© The Authors 2024

1 Introduction

At the low temperatures of the dense regions of the interstellar medium (10–20 K), and in the absence of any external energy source, all interstellar gas containing heavy elements should be depleted in interstellar dust grains (Willacy & Millar 1998). While such a picture is ideal, it highlights the immense importance of interstellar dust grains in the chemical evolution of our cosmos. At the temperatures mentioned above, observations show that dust grains are covered by a thick ice mantle made of different molecules (Boogert et al. 2015) such as H₂O, CO, CO₂, and CH₃OH, along with minor fractions of the so-called complex organic molecules (COMs), broadly defined as organic molecules with six atoms or more (Herbst & van Dishoeck 2009). These COMs are fundamental to understanding the evolution of the chemical complexity taking place in the very early stages of star and planet formation.

The formation of the ice mantle atop dust grains is sequential. Although H₂O is formed via surface reactions even in an early stage of molecular clouds (Tielens 1982; Dulieu et al. 2010; Ioppolo et al. 2008; Lamberts et al. 2014, 2013; Meisner et al. 2017; Molpeceres et al. 2019), CO is formed in the gas phase and accretes to the grains when the cloud reaches a critical density (Boogert et al. 2015). This phenomenon, known as catastrophic CO freeze-out, is demonstrated by astronomical observations of prestellar cores by an abrupt decay of gas CO abundance (Caselli et al. 1999; Cazaux et al. 2011; Taquet et al. 2012). In fact, in the traditional picture of the formation of COMs, the onset of organic chemistry is the hydrogenation of accreted CO to CH₃OH (Watanabe & Kouchi 2002; Fuchs et al. 2009; Simons et al. 2020; Ferrero et al. 2023). This picture is still prevalent today, although it is complemented by other theories of COMs’ formation, such as carbon atom accretion (Molpeceres et al. 2021a; Potapov et al. 2021; Ferrero et al. 2023, 2024) or gas-phase formation (Redondo et al. 2017; Cooke & Sims 2019; Ceccarelli et al. 2022).

One conundrum of interstellar surface chemistry is the impossibility of heavy species roaming the surface of interstellar ice. Although various COMs are detected in cold prestellar cores (e.g., Bacmann & Faure 2014; Jiménez-Serra et al. 2016), diffusive reaction mechanisms such as Langmuir-Hinshelwood are limited by thermal diffusion of the radicals, which can be negligibly low at ~10 K. Several alternative mechanisms may explain the presence of COMs in ice mantles and their subsequent return to the gas phase. Some of them are bimolecular reactions through the Eley-Rideal mechanism (Le Bourlot et al. 2012; Ruaud et al. 2015; Ferrero et al. 2023) and nondiffusive chemistry (Ioppolo et al. 2021; Jin & Garrod 2020; Simons et al. 2020; Garrod et al. 2022) or energetic mechanisms such as UV or e⁻ irradiation (Munñoz Caro et al. 2002; Bennett et al. 2006; Muñoz Caro et al. 2014; Turner et al. 2018; Shingledecker et al. 2017; Ishibashi et al. 2021).

The explanation for the emergence of COMs we investigate in this paper is the diffusive chemistry on CO ices. We have kept the Langmuir-Hinshelwood mechanism but modified the substrate. On CO ices, the intermolecular forces that bind the adsorbates to the surface are of a lower magnitude due to the very low value of the dipole moment of the CO molecule |μ| = 0.122 D (Muenter 1975), in stark contrast to the one of H₂O, 1.85 D. Lower forces in CO translate into weaker binding and enhanced diffusion of adsorbates on its surface, compared to H₂O ice. While Bergin et al. (1995) investigated the effect of weaker binding on CO-covered ices on molecular depletion in molecular clouds, in the present work we investigated the effect of enhanced diffusion on the formation of COMs. However, due to its weak binding energy (BE), pure CO ice desorbs thermally at ~18 K (Bisschop et al. 2006; Fuchs et al. 2006). Thus, the CO-dominant ice exists only in a narrow temperature range, in which only specific radicals and molecules can diffuse. Nevertheless, laboratory experiments do show rich COM chemistry on CO ice, mainly related to its hydrogenation sequence (Fedoseev et al. 2015, 2016, 2017; Simons et al. 2020; Chuang et al. 2020).

For this work, we ran astrochemical models to investigate the effect of the enhanced diffusion of certain species on CO ice in the aforementioned narrow temperature window. We performed this task using an average proportionality factor ϵ_CO between BE in CO and H₂O that is dynamically updated as a function of the dust surface composition. Although most molecules are unaffected by the update of binding energies as a result of a low temperature and competition with CO desorption, the diffusion of some key radicals is enabled as the CO coverage increases. In Furuya et al. (2022), we determined that an increase in the diffusion of six radicals, namely CH₃, CH₂, HCO, O, N, and NH, affected the chemistry of cold cores significantly. For these radicals, we derived an individualized scaling factor. We investigated not only the specific influence of CO ice on diffusive chemistry but also how the chemistry couples with other physical magnitudes of the prestellar cores, such as the dust temperature and the cosmic ray ionization rate ζ. Our main aim is to qualitatively show how the surface composition of interstellar ice can modulate molecular abundances in the interstellar medium (ISM). This work serves as a stepping stone for the development of models that recognize the heterogeneous composition of interstellar ice.

Our paper is structured as follows. Section 2 discusses our methodology for astrochemical modeling. We present our results of the grid of astrochemical models of prestellar cores in Sec. 3, and discuss how the surface composition affects the abundances of certain species. We summarize our main findings in Sec. 4 and contextualize our results in the tapestry of modern astrochemical observations.

2 Astrochemical models

We used our in-house astrochemical model ROKKO described in Furuya et al. (2015) and Furuya & Persson (2018) using the reaction network and grain surface parameters presented in Garrod (2013) and Belloche et al. (2014). Briefly, the Garrod (2013) reaction network comprises almost 1300 species considering gas, surface, and mantle species, coupled through 11,800 reactions. We made several modifications to the network model. The most relevant are setting the BE of C and H atoms to 10000 K (Wakelam et al. 2017; Shimonishi et al. 2018; Duflot et al. 2021; Molpeceres et al. 2021a) and 400 K (Wakelam et al. 2017, M06-2X values), setting the formation of CO₂ through the HOCO radical (Molpeceres et al. 2023), and expanding the recombination of cations and grains upon collision with negatively charged grains to all cations. For convenience, the BE of H₂ is set equal to the one of H, that is, 400 K. We note that $E_{\mathrm{b}\mathrm{i}\mathrm{n},{\mathrm{H}}_2}^{\mathrm{C}\mathrm{O}}$ is lower than ${\mathrm{E}}_{\text{bin, }{\mathrm{H}}_2}^{{\mathrm{H}}_2\mathrm{O}}$, as described below, while we set the same BE for H and H₂ (i.e., 400 K) on H₂O ice. Finally, we included the activation energies for radical recombinations in Enrique-Romero et al. (2022). In Enrique-Romero et al. (2022) the authors report activation energies for some radical-radical recombinations on H₂O ice at different binding sites and in the gas phase. Among other conclusions, the activation energies reported in Enrique-Romero et al. (2022) are binding site dependent. Due to the absence of an H-bond network, in CO ices these activation energies should vanish, as reported by Lamberts et al. (2019) for CH₃CHO. The applicable activation barriers, when present, will be closer to the gas phase values. Therefore, we have adopted the gas phase values of Table 3 of Enrique-Romero et al. (2022) in our rate equation model.

In order to investigate how chemistry varies with the ice surface composition, the BE of a species is calculated at each time step as follows. $$E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{i}}=(1-\theta )E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{i}}^{{\mathrm{H}}_2\mathrm{O}}+\theta E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{i}}^{\mathrm{C}\mathrm{O}},$$(1)

with $E_{\mathrm{b}\mathrm{i}\mathrm{n}}^{{\mathrm{H}}_2\mathrm{O}}$ the BE of the adsorbate on amorphous solid water, θ the CO coverage at a given time and $E_{\mathrm{b}\mathrm{i}\mathrm{n}}^{\mathrm{C}\mathrm{O}}$ the CO ice counterpart. Surface coverage, θ, is defined as the fractional abundance of CO on the surface of the ice, which in our models corresponds to 4 monolayers; see below. This approach was used in previous works (Garrod & Pauly 2011; Taquet et al. 2014; Furuya & Persson 2018).

In general, $E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{i}}^{\mathrm{C}\mathrm{O}}$ is unknown because the available data are much more scarce, although several studies exist (Bisschop et al. 2006; Fuchs et al. 2006, 2009; Lamberts et al. 2019; Mondal et al. 2021; Sil et al. 2021; Molpeceres et al. 2021b; Ferrari et al. 2023; Ferrero et al. 2023). For pure diffusion energies, the only available comparison is the one between Hama et al. (2012) (on H₂O) and Kimura et al. (2018) (on CO) for the H atom, finding similar values in the two substrates (See below).

To estimate the influence of CO ices in the chemistry of the ISM, we define the binding energies on CO in our model as: $$E_{\text{bin},\mathrm{i}}^{\mathrm{C}\mathrm{O}}={ϵ}_{\mathrm{C}\mathrm{O}}{E}_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{i}}^{{\mathrm{H}}_2\mathrm{O}}.$$(2)

In the above definition, ϵ_CO is the proportionality factor between the binding energies on pure amorphous solid water (ASW) and pure CO ices. In Appendix A we show a derivation of ϵ_CO for the six radicals deemed fundamental for molecular cloud chemistry, that is, CH₃, HCO, N, O, CH₂, NH (Furuya et al. 2022), using quantum chemical calculations. A particular case is the BE of the H atom, taken from the study of adsorption on a dimer by Wakelam et al. (2017). In all our models, the BE of H is independent of θ, e.g., it is not updated according to Eq. (1) based on the similarity between the values found in Wakelam et al. (2017) for H₂O and in Fuchs et al. (2009) for CO ice, also found for diffusion barriers comparing the works of Hama et al. (2012) and Kimura et al. (2018). Therefore, we set $E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{H}}^{\mathrm{C}\mathrm{O}}=E_{\mathrm{b}\mathrm{i}\mathrm{n},\mathrm{H}}^{{\mathrm{H}}_2\mathrm{O}}=400\mathrm{K}$ in our models. For all the remaining adsorbates not included in the above list (and H₂, due to its importance), we set an ϵ_CO=0.66 (Furuya & Persson 2018, and references therein), i.e. the same as the ratio between the BE of CO on H₂O ice and CO ice. In Appendix B, we show that the selection of a universal ϵ_CO for these species does not affect the subsequent conclusions of our calculations. However, we note that 0.66 is already a conservative value and that our quantum chemical estimations suggest that ϵ_CO should be lower for a wide amount of molecules (Appendix A). The values of ϵ_CO are summarized in Table 1. In our chemical models, we compared models that enable Eq. (1) (i.e. considering ϵ_CO) with models where the BE is constant and equivalent to the binding energies on ASW. We track how ϵ_CO and θ affect the evolution of molecules in prestellar cores.

Several nonthermal mechanisms are included in the models. The reactive desorption fraction is given by the scheme of Garrod et al. (2007): $$f=\frac{ap}{1+ap}$$(3) $$p={[1-\frac{E_{\mathrm{b}\mathrm{i}\mathrm{n}}}{E_r}]}^{s-1},$$(4)

with an a parameter of 0.01, and a reaction energy, E_r. The parameter a corresponds to the fraction of the reaction energy that is inoculated into the desorption vibration mode. Contrary to thermal diffusion or desorption, which we made dependent on the surface coverage of CO through Eq. (1), we did not scale f accordingly. The reason for this is that the literature does not contain data on a on CO ice. Other nonthermal mechanisms in our model are cosmic-ray induced desorption that follows the Hasegawa-Herbst formulation (Hasegawa & Herbst 1993) sampling cosmic ray ionization rates (ζ) in the range between 1×10⁻¹⁸−1×10⁻¹⁵ s⁻¹ (50 logarithmic equispaced values). Photodesorption is included using a constant desorption yield of 1×10⁻³ per photon with a constant extinction coefficient A_v=100 mag, so CR-induced UV photons dominate the UV radiation field. The gas particle density is set to 2×10⁵ cm⁻³. The temperatures of gas and dust (T_g and T_d) are assumed to be the same and are sampled in the range between 8–18 K (50 values). Above that temperature, the thermal desorption of CO becomes important, and the impact of CO coverage on COMs formation is diminished. Combining the sampling of ζ and T_d, our grid of models consists of 2500 models, that are carried out varying and without varying ϵ_CO for comparison, for a total of 5000 models.

Following the previous models of the gas–grain reaction network, the relationship between E_bin and E_diff (the diffusion energy) is assumed to be E_diff/E_bin = 0.4, e.g. a constant value. The drawbacks of this approach are discussed in Furuya et al. (2022), but this choice allows us to keep the parameter space search (e.g. ζ and T_d) tractable. While quantitative conclusions are indeed affected by a different choice of E_diff/E_bin, our qualitative conclusions about the enhancement of the chemistry in CO are not affected by this model choice as we show in Appendix C. We considered the reaction-diffusion competition scheme in all our models (Chang et al. 2007). Finally, in our models, we set the four topmost (N_act = 4) layers as chemically active and regarded as surface species. On the contrary, the mantle chemistry below these four monolayers was considered inactive. A summary of the physical parameters of our simulations is presented in Table 2.

The initial conditions of the models are selected to maximize the impact of CO chemistry; we started with a significant fraction of C and O nuclei locked in gas-phase CO (before freeze-out) and H₂O ice. The elemental abundances and the initial molecular/atomic abundances are presented in Table 3. The elemental abundances are taken from Aikawa & Herbst (1999), low-metallicity values. We considered a single population of interstellar grains of size 0.1 μm with approximately 10⁶ binding sites per grain.

We defined our reference model with parameters of ϵ_CO value of 1.0, T_d of 10.0 K and ζ of 1.3×10⁻¹⁷ s⁻¹. The temporal variations of molecular abundances are shown in Fig. 1 for the main ice species. In the bottom panel of Fig. 1, we show the fractional surface coverage of the main ice species. In the model, we observe a clear and abrupt depletion of CO gas around 5×10⁴ yr and the formation of an apolar-rich ice surface in a time window between 4×10⁴ and 3×10⁵ yr. Thus, in Sec. 3, we present molecular abundances at 2×10⁵ yr for a grid of models. This time corresponds to an average lifetime of a pre-stellar core (André et al. 2014) and also a time after most of the catastrophic CO freeze-out has taken place when the effects of the CO chemistry are more apparent. The surface coverage of CO (θ) peaks at 6×10⁴ yr, θ ~0.65, with CH₄, N₂ and H₂CO as other main components. Because the freeze-out time scale depends only slightly on the temperature at 10–20 K, the choice of 2×10⁵ yr is also justified for the rest of the surface temperatures.

Fig. 1 Reference model. Top: evolution of the abundances of main ice species, surface H2O and CO, and gas CO in our reference model. Bottom: the fractional surface coverage of the main ice species corresponds to the fraction of molecules in the four topmost layers of a dust grain.

3 Results

3.1 Initial analysis

3.1.1 Diffusion times as a function of time and T_d

Qualitative differences in the chemistry on H₂O-dominant ices and CO-dominant ices would be caused by species whose diffusion rate is negligible on H₂O ice but significant on CO ice. Hence it is beneficial to check the diffusion timescale (τ, i.e., inverse of the rate constant for diffusion): $$\tau =1/(\nu \exp \left[\frac{-0.4((1-\theta )E_{\text{bin,i}}^{{\mathrm{H}}_2\mathrm{O}}+\theta E_{\text{bin},\mathrm{i}}^{\mathrm{C}\mathrm{O}}}{T_{\mathrm{d}}}\right]),$$(5)

of the key radicals considered in this work on pure H₂O ice and pure CO ice (Fig. 2), as a function of T_d and θ_CO. In Eq. (5) we assume an E_diff/E_bin ratio of 0.4 and a preexponential factor of diffusion (ν) of 1×10¹² s⁻¹. Then τ must be confronted against the competitive process that is the hydrogenation of the radicals, whose timescale depends on the H accretion (and diffusion rate). The H accretion frequency can be determined as $$f_{\mathrm{H}}=\frac{1}{\sigma v_{\mathrm{H}}{n}_{{\mathrm{H}}_{\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}}}},$$(6)

where σ is the cross-section of a dust grain, υ_H is the average thermal velocity of the H atom at 10 K and $n_{{\mathrm{H}}_{\text{atom }}}$ is the number density of H atoms (~1 cm⁻³ for ζ=1.3×10⁻¹⁷ s⁻¹, Goldsmith & Li 2005). This timescale is depicted in Fig. 2 as a dotted black line. Because f_H depends strongly on ζ, we consider $n_{H_{\text{atom }}}$ in a region that spans a factor 5 from the above value (represented as a shaded gray box in Fig. 2).

Figure 2 (top panel) indicates that below 12 K and in a θ=0.65 CO mixed ice, only species with effective E_diff below 800 K may experience rapid diffusion, for example N, O or CH₂. We note that mobility of radicals such as CH₂ or O is not negligible even below 10 K. Slow but plausible diffusion of species with binding energies close to 1500 K is also enabled (i.e. CH₃) close to 12 K. Above 14 K, the diffusion of molecules with binding energies around 1500 K is fast. Between 12–14 K, radicals with binding energies higher than 2000 K start to diffuse slowly (HCO, NH). The bottom panel of Fig. 2 is informative to rationalize a “critical” CO coverage at which the chemistry on CO is truly enabled. This surface coverage is roughly 0.60, which is close to the maximum ~0.65 found in our reference model (Fig. 1, bottom panel). We note that owing to the differences in physical conditions within the grid of models to be presented in this section, the maximum θ can vary between models; see below. However, a valuable piece of information derived from this analysis is the threshold concentration of CO molecules on the surface to activate new chemistry.

In terms of chemistry, including a coverage-dependent term in the diffusion allows for rapid diffusion of molecules below 12 K, as mentioned above. Therefore, the abundance of molecules whose parent species are CH₃ (and CH₂ as the parent of CH₃), or O and N atoms, are susceptible to their formation of the apolar ice layer. Above 12 K, HCO or NH also begin to diffuse, competing with the rapid consumption of CO and the consequential increase in E_bin. Therefore, mainly due to the delicate balance established by CO consumption through several mechanisms (hydrogenation, desorption), the species with the most likely impact on the chemistry of apolar ices at very low temperatures are CH₃, CH₂, N and O because their diffusion rate is significant considering significant CO ice coverage while they are negligible with θ=0.

In the following sections, we present the results of our chemical model for some selected species and discuss how the CO surface coverage affects these specific molecules. We discuss the chemical reactions that are enabled on CO ice and that lead to the different abundances of the target species. The selection of specific molecules affected by the different chemistry on CO ice does not imply that these are the only molecules susceptible to enhanced chemistry on apolar ice. However, selecting key molecules presented in the extensive Garrod (2013) chemical network is necessary to keep our discussion as general as possible.

Fig. 2 Initial analysis of diffusion timescales. Top: diffusion timescale as a function of the dust temperature (Td) for θ = 0.65 Bottom: as a function of the surface coverage of CO (θ) at 12 K for the main radicals considered (Furuya et al. 2022). The dotted black line and the shaded region around it show the H accretion timescale for different $n_{{\mathrm{H}}_{\text{atom }}}$ (see Sec. 3.1).

3.1.2 The importance of ζ

As described in Sec. 2 and Table 2, we vary two physical parameters in our grid of models: surface temperature, T_d, and cosmic ray ionization rate, ζ. Although the effect of T_d is evident from Eq. (5), the importance of ζ may not be obvious. In fact, ζ is the most important factor affecting CO coverage and CO-related chemistry on the ice surface.

At low temperatures, CO is one of the few abundant closed shell species that can be chemically converted through the chain of (forward and backward) reactions (Watanabe & Kouchi 2002; Fuchs et al. 2009): $$\mathrm{C}\mathrm{O}\stackrel{+\mathrm{H}}{⟶}\mathrm{H}\mathrm{C}\mathrm{O}\stackrel{+\mathrm{H}}{⟶}{\mathrm{H}}_2\mathrm{C}\mathrm{O}\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}\mathrm{H}}_3\mathrm{O}\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}$$(R1) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\stackrel{+\mathrm{H},-{\mathrm{H}}_2}{⟶}{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}\stackrel{+\mathrm{H},-{\mathrm{H}}_2}{⟶}{\mathrm{H}}_2\mathrm{C}\mathrm{O}\stackrel{+\mathrm{H},-{\mathrm{H}}_2}{⟶}\mathrm{H}\mathrm{C}\mathrm{O}\stackrel{+\mathrm{H},-{\mathrm{H}}_2}{⟶}\mathrm{C}\mathrm{O}.$$(R2)

Therefore, the abundance of CO on the surface is inversely correlated with the presence of H atoms. The steady-state H abundance is determined by ζ (Goldsmith & Li 2005), because cosmic rays dissociate H₂ molecules very effectively into two H atoms that later accrete on the ice and reduce the abundance of CO. In summary, in models with higher ζ, CO-coverage is lower, reducing COM formation through efficient diffusion on apolar ices. This fact is slightly counterintuitive because higher ζ is known to produce more radicals on H₂O due to the fragmentation of molecules. However, in CO ices, removing the CO layer by H atoms overwhelms the latter effect. As we show in the next section CH₃OH is an exception to this rule, precisely because it is the endpoint of the CO hydrogenation sequence.

The other magnitude that depends on ζ is the cosmic-ray-induced desorption that returns molecules to the gas phase via stochastic heating after the collision (Hasegawa & Herbst 1993). This effect indirectly depends on the binding energy and synergizes with chemical desorption. It will be briefly discussed in Sec. 3.4. However, we found it to be less critical in the context of the chemistry of CO ices than in the production of H atoms discussed above.

3.2 Chemistry on CO activated at ≤ 12K

We start our analysis with an examination of the abundances of molecules originating from diffusion below 12 K (Fig. 2). As we introduced at the end of Sec. 2, our analysis of molecular abundances is performed at 2×10⁵ yr. Disentangling the influence of surface composition requires careful examination of many factors that influence surface chemistry in interstellar environments. These include not only the binding energies updated by the chemical composition of the surface but also how this effect couples with the cosmic ray-induced processes. It is customary to mention that E_diff/E_bin is an important factor in the model (Furuya et al. 2022). We assume a constant value for all of our simulations. This is a caveat of our models. However, we note that we can estimate at least qualitatively which molecules and radicals would play a critical role in the CO-dominated ice as described in Sec. 3.1 when a different ratio of E_diff/E_bin is adopted (Appendix C).

Figure 3 shows the abundances of assorted molecules at t=2×10⁵ yr as a function of ζ and T_d for models of ϵ_CO ≠ 1.0 and ϵ_CO=1.0. The latter corresponds to the classical model, in which the adsorption and diffusion energies do not depend on the surface composition of the ice. We selected molecules directly related to the radicals analyzed in Sec. 3.1 whose diffusion timescale is significantly affected by the composition of the surface ice. In the following, we highlight CH₃CHO, CH₃OCH₃, C₂H₅OH, CH₃COOH, C₂H₆, and CH₃NH₂. Although the first four molecules are sensitively dependent on ϵ_CO, the last two species are less affected by a variety of factors that we detail in the following paragraphs. Other molecules, such as CH₃OH, CH₃CN, or C₂H₄, are also described as examples of unaffected molecules.

CH₃CHO (Acetaldehyde). Acetaldehyde is the molecule most affected by the update of the binding energies. Introducing the θ dependence in our models can increase the abundance of CH₃CHO up to four orders of magnitude in the temperature range investigated in the present work. There are three related factors.

First, diffusion of CH₃ and CH₂ on a CO-dominated surface enables the surface reactions: $${\mathrm{C}\mathrm{H}}_3+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{H}\mathrm{O}$$(R3)

$${\mathrm{C}\mathrm{H}}_2+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{H}\mathrm{O}\stackrel{\mathrm{H}}{⟶}{\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{H}\mathrm{O},$$(R4)

Reaction R3 was investigated on pure H₂O ice by Enrique-Romero et al. (2021), finding a critical dependence on the E_diff/E_bin, with respect to the competitive channel: $${\mathrm{C}\mathrm{H}}_3+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_4+\mathrm{C}\mathrm{O}.$$(R5)

Enrique-Romero et al. (2021) find a dominance of abstraction (Reaction R5) over addition (Reaction R3) for the reaction at a cavity on H₂O ice with E_diff/E_bin=0.3, due to the interactions of the adsorbates with the H₂O substrate in a particular binding mode. This dominance vanishes for E_diff/E_bin=0.4. In general, the competition between addition and abstraction on H₂O ice, not only for CH₃CHO but for many COMs, depends on the binding site (Enrique-Romero et al. 2016; Butscher et al. 2017; Enrique-Romero et al. 2022). Mixed interstellar ice or CO ices have weaker binding sites than pure H₂O ice so the contribution of the particular binding mode becomes less significant. Since both reactions R3 and R5 do not have an intrinsic activation barrier (i.e., they are barrierless in the gas phase), we expect an equal branching ratio, that is, 0.5 each. The lack of intrinsic barriers was also found by Lamberts et al. (2019) for the formation of CH₃CHO on the CO ice. Although Enrique-Romero et al. (2021) relies on the uncertainty of the diffusion energies of CH₃ and HCO to postulate the possibility of formation of CH₃CHO, we consider that the surface coverage of CO plays an active role in promoting the reaction. Our conclusions are complementary and have different applicability, depending on the evolutionary stage of the molecular cloud, for example, depending on the CO/H₂O ratio on the ice surface.

The second reason for the enhanced formation of CH₃CHO is the abundance of HCO on the surface of the grains. Compared to other radicals, HCO is relatively abundant thanks to CO hydrogenation (Eq. (R1)).

The third reason is that, at slightly higher temperatures, HCO also starts to be relatively mobile, enhancing the reaction (R3). All these three factors allow the formation of acetaldehyde at low and intermediate ζ, in the entire temperature range considered in this study.

CH₃OCH₃ (dimethyl ether). The arguments presented for CH₃CHO can be directly transferred to CH₃OCH₃. Therefore, we observe a remarkable similarity between the graphs corresponding to these molecules in Fig. 3. The difference between these two species is the second radical that participates in the reaction, that for CH₃OCH₃ is CH₃O: $${\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_3\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_3{\mathrm{O}\mathrm{C}\mathrm{H}}_3,$$(R6)

with the competing channel with an activation barrier of 204 K (Enrique-Romero et al. 2022): $${\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_3\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_4+{\mathrm{H}}_2\mathrm{C}\mathrm{O}.$$(R7)

Therefore, the lower abundance of CH₃OCH₃ to CH₃CHO is due to the lower abundance of its parent radical CH₃O than HCO, about one order of magnitude. The higher BE of CH₃O also makes it effectively immobile, unlike what happens with HCO, reducing the relative abundance with respect to CH₃CHO at higher T_d. Both HCO and CH₃O are formed through the CO hydrogenation sequence (Watanabe & Kouchi 2002; Fuchs et al. 2009; Simons et al. 2020).

C₂H₅OH (ethanol). The set of the three molecules most affected by CH₃ interacting with hydrogenated CO products concludes with ethanol, C₂H₅OH. Like in the two previous molecules, the enhanced reaction on CO ice is related to the hydrogenation sequence on CO ice (Eq. (R1)). The reaction that enhances the formation of C₂H₅OH is: $${\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}⟶{\mathrm{C}}_2{\mathrm{H}}_5\mathrm{O}\mathrm{H},$$(R8)

that again competes with the following reaction with an 1170 K barrier: $${\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_4+{\mathrm{H}}_2\mathrm{C}\mathrm{O}.$$(R9)

The very similar profiles found for CH₃CHO, CH₃OCH₃, and C₂H₅OH in Fig. 3 confirm this explanation, and, as discussed for CH₃OCH₃, variations between them can be attributed to variations in the abundance of the parent species.

CH₃COOH (acetic acid). Acetic acid abundance is greatly enhanced by CO ices. In contrast with the ternary of previous molecules, the reason for its increase is due to the enabling of the reaction $$\mathrm{C}\mathrm{O}+\mathrm{O}\mathrm{H}⟶\mathrm{H}\mathrm{O}\mathrm{C}\mathrm{O},$$(R10)

following the recent study of Molpeceres et al. (2023). Reaction (R10) was previously considered to form directly CO₂ (Ioppolo et al. 2011; Garrod & Pauly 2011; Oba et al. 2010; Terwisscha van Scheltinga et al. 2022), but in our recent study we showed otherwise. Hence, the chemistry of HOCO can be initiated on the surface of the CO-rich ice. Associated with the selective production of HOCO, we also find an increase in the diffusion of CO (ϵ=0.66) in mixed ices (θ ≤ 0.5) and self-diffusion in CO-rich ices (θ>0.5), further increasing the production of HOCO. With a high concentration of HOCO radicals on the surface, the formation of CH₃COOH can proceed through: $$\mathrm{H}\mathrm{O}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_3⟶{\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}.$$(R11)

Both reactions (R10) and (R11) are facile on CO ice in our models, in contrast to H₂O ice, which requires a threshold temperature of (~ 12 K) for the diffusion of CO (Garrod & Pauly 2011). On CO ices, the lack of this threshold summed to the easy diffusion of CH₃ is what causes the dramatic increase in the production of this molecule. However, we note the subtle differences in the HOCO formation barriers on H₂O and CO ices through reaction (R10) (Molpeceres et al. 2023). Regretfully, we cannot capture these differences with our models yet.

C₂H₆ (ethane). Ethane is an interesting molecule in this context because it presents three different regimes in its abundance heat maps in Fig. 3. The formation of ethane on CO ice is enhanced by the surface reaction: $${\mathrm{C}\mathrm{H}}_3+{\mathrm{C}\mathrm{H}}_3⟶{\mathrm{C}}_2{\mathrm{H}}_6\text{, }$$(R12)

primarily around 12 K with deviations between ϵ_CO=1.0 and ϵ_CO ≠ 1.0 of roughly a factor 3. However, above 12 K, we observe an increase in the production of C₂H₆ that cannot be explained by the recombination of CH₃ alone. A careful inspection of the reaction rates shows that the reaction: $${\mathrm{C}}_2{\mathrm{H}}_5+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{C}}_2{\mathrm{H}}_6+\mathrm{C}\mathrm{O}.$$(R13)

is responsible for the increase in C₂H₆ abundance at medium temperatures (and medium ζ), promoted by the diffusion of HCO above 12 K. The less marked influence of the CO ice for this molecule is due to the importance of the reaction $${\mathrm{C}}_2{\mathrm{H}}_5+\mathrm{H}⟶{\mathrm{C}}_2{\mathrm{H}}_6$$(R14)

across our parameter space, which is the last step in the hydrogenation of C₂H₂ (Kobayashi et al. 2017; Molpeceres & Rivilla 2022). Although reactions R12 and R13 do enhance the formation of C₂H₆, they are secondary channels. In this study, we found that, when a molecule can be formed by hydrogenation of abundant species, the influence of CO ice chemistry is minor or null and only molecules without a clear hydrogenation route have a significant increase in abundance (at low ζ).

CH₃NH₂ (methylamine). In Fig. 3, we observe that the models for ϵ_CO ≠1.0 and ϵ_CO=1.0 present a similar profile but are slightly vertically displaced, meaning that the same chemistry begins to operate, but at a lower temperature for ϵ_CO ≠ 1.0. Reactions that are common to the ϵ_CO ≠ 1.0 and ϵ_CO=1.0 models are $${\mathrm{C}\mathrm{H}}_3\mathrm{N}\mathrm{H}+\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_3{\mathrm{N}\mathrm{H}}_2$$(R15) $${\mathrm{C}\mathrm{H}}_2{\mathrm{N}\mathrm{H}}_2+\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_3{\mathrm{N}\mathrm{H}}_2\text{. }$$(R16)

The enhanced formation of CH₃NH₂ at low temperatures on CO ice is a consequence of the following reactions: $${\mathrm{C}\mathrm{H}}_3+\mathrm{N}⟶{\mathrm{C}\mathrm{H}}_2\mathrm{N}\mathrm{H}$$(R17) $${\mathrm{C}\mathrm{H}}_3+\mathrm{N}\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_3\mathrm{N}\mathrm{H}$$(R18) $${\mathrm{C}\mathrm{H}}_2+{\mathrm{N}\mathrm{H}}_2⟶{\mathrm{C}\mathrm{H}}_2{\mathrm{N}\mathrm{H}}_2$$(R19) $${\mathrm{C}\mathrm{H}}_2+\mathrm{N}\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_2\mathrm{N}\mathrm{H}$$(R20) $${\mathrm{C}\mathrm{H}}_3+{\mathrm{N}\mathrm{H}}_2⟶{\mathrm{C}\mathrm{H}}_3{\mathrm{N}\mathrm{H}}_2.$$(R21)

The first four products are further hydrogenated to CH₃NH₂. The heatmap in Fig. 3 is vertically displaced because the diffusion of N atom is activated early on water ice and CO ice only plays a role at a very low temperature. At higher ζ, NH₃ formation is favored, and NH_x (x≤2) abundance is lower on the surface. At high T_d, the population of H atoms on surfaces is lower due to thermal desorption, which decays the abundance of CH₃NH₂. Thus, CH₃NH₂ formation is slightly enhanced on CO ices at low temperatures and ζ. However, the effect of the CO ice in the predicted abundance of CH₃NH₂, is smaller than in the O-bearing COMs described above.

Unaffected molecules: CH₃OH (methanol), C₂H₄ (ethylene), CH₄ (methane) and CH₃CN (methyl cyanide). The abundances of CH₃OH, CH₃CN, and CH₄ at time 2×10⁵ yr for ϵ_CO=1.0 and ϵ_CO ≠ 1.0 are presented in Fig. 4. We observe an almost null effect of the CO ice layer on the abundances, CH₃OH, C₂H₄, and CH₃CN, or a decrease in the production of the molecule on CO ice (CH₄). We found common reasons for all the cases: the competition between hydrogenations and reactions with heavier radicals. For CH₃OH, the CO hydrogenation (sequence reactions (R1)) dominates the whole chemistry of the molecule and alternative reactions, like $${\mathrm{C}\mathrm{H}}_3+\mathrm{O}\mathrm{H}⟶{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}$$(R22)

are minor. Similar arguments can be found for C₂H₄ and CH₃CN, whose chemistry is dominated by the chain hydrogenation reactions: $${\mathrm{C}}_2{\mathrm{H}}_2\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}}_2{\mathrm{H}}_3\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}}_2{\mathrm{H}}_4\left(\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}}_2{\mathrm{H}}_5\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}}_2{\mathrm{H}}_6\right)$$(R23) $${\mathrm{C}}_2\mathrm{N}\stackrel{+\mathrm{H}}{⟶}\mathrm{H}\mathrm{C}\mathrm{C}\mathrm{N}\stackrel{+\mathrm{H}}{⟶}{\mathrm{H}}_2\mathrm{C}\mathrm{C}\mathrm{N}\stackrel{+\mathrm{H}}{⟶}{\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{N}.$$(R24)

Alternative reactions like $${\mathrm{C}\mathrm{H}}_2+{\mathrm{C}\mathrm{H}}_2⟶{\mathrm{C}}_2{\mathrm{H}}_4$$(R25) $${\mathrm{C}\mathrm{H}}_3+\mathrm{C}\mathrm{N}⟶{\mathrm{C}\mathrm{H}}_3\mathrm{C}\mathrm{N},$$(R26)

are minor. Interestingly, and as mentioned in the previous section, C₂H₆ has a small increase in its abundance on CO ices, in comparison with the related C₂H₄ molecule. The reason is that reaction R12 contributes to a small extent to the chemistry of C₂H₆, possibly due to the larger abundance of CH₃ radicals on the surface. Finally, the formation of CH₄ is enhanced on H₂O ice over CO ice owing to the lower abundance of CH₃ radicals on the surface that are used in all CO enhanced reactions presented above.

The absence of important alternative routes of formation on CO ice for the molecules in this section makes us conclude that enhanced chemistry is only applicable to molecules that cannot be formed by hydrogenation. When hydrogenation is possible, and because the diffusion rate of the H atom on ice does not depend on the ice’s surface composition in our models, both the updated model and the classical one yield the same results.

Fig. 3 Ice abundances at t=2×105 yr in models of various ζ and Td. The selected molecules are CH3CHO, CH3OCH3, C2H5OH, CH3COOH, C2H6, and CH3NH2.

Fig. 4 Same as Fig. 3 but for CH3OH, C2H4, CH4, and CH3CN.

3.3 Chemistry on CO activated at ≥ 14 K

The second grand set of molecules affected by the surface composition of the ice substrate refers to molecules related to HCO and NH that slowly diffuse below 12 K, but start to be rather mobile at ~14 K. The impact of CO ice on the chemistry of these molecules is smaller than those presented in Sec. 3.2 because HCO and NH need to compete with the chemistry initiated by CH₃, CH₂, O and N that have lower diffusion energies and especially by H because both HCO and NH are easily hydrogenated. The effect of the CO ice is, however, noticeable. We explicitly list some molecules here, namely NH₂CHO (formamide), HCOOCH₃ (methyl formate), OHCCHO (glyoxal), HOCH₂CHO (glycolaldehyde), and NH₂CH₂COOH (glycine), whose abundances are pictured in Fig. 5.

The explanation for the enhanced formation on the CO ice is, in this case, common to all molecules except glycine and begins to apply mostly above 12–13 K, coinciding with the start of a slow but noticeable diffusion of HCO and NH (Fig. 2). This makes an itemized description of the reactivity, as in the previous section, reiterative. The enhanced abundance of these molecules arises from radical couplings between the parent radicals and HCO/NH. For the species portrayed in Fig. 5 these reactions are, where the competing reactions are in parentheses: $${\mathrm{N}\mathrm{H}}_2+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{N}\mathrm{H}}_2\mathrm{C}\mathrm{H}\mathrm{O}$$(R27) $$({\mathrm{N}\mathrm{H}}_2+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{N}\mathrm{H}}_3+\mathrm{C}\mathrm{O})$$(R28) $$\mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_3\mathrm{O}⟶{\mathrm{H}\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{C}\mathrm{H}}_3$$(R29) $$(\mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_3\mathrm{O}⟶{\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{C}\mathrm{O})$$(R30) $$\mathrm{H}\mathrm{C}\mathrm{O}+\mathrm{H}\mathrm{C}\mathrm{O}⟶\mathrm{O}\mathrm{H}\mathrm{C}\mathrm{C}\mathrm{H}\mathrm{O}$$(R31) $$(\mathrm{H}\mathrm{C}\mathrm{O}+\mathrm{H}\mathrm{C}\mathrm{O}⟶{\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{C}\mathrm{O})$$(R32) $$\mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}⟶{\mathrm{H}\mathrm{O}\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{H}\mathrm{O}$$(R33) $$(\mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}⟶2{\mathrm{H}}_2\mathrm{C}\mathrm{O}).$$(R34)

The dependence on T_d is now clearer to explain; all chemistry is activated by the diffusion of HCO/NH₂. On the other hand, the ζ dependence is precisely the same as in Sec. 3.2. With more available H atoms, diffusive chemistry on CO ices converges to diffusive chemistry on H₂O¹, because heavy radical diffusion is hampered by the hydrogenation of the radicals. The small differences between the molecules in Fig. 5 are attributed to subtle differences in the abundances of parent molecules and precursors at the given T_d. It is more difficult to pinpoint the exact reason for each case.

Even with the increased formation of COMs on CO ice, the species whose parent species is HCO have lower abundances than the ones related with CH₃/CH₂. This is especially true for HCOOCH₃ and HOCH₂CHO. For NH₂CHO, an important prebiotic precursor (López-Sepulcre et al. 2019), the abundances predicted by the CO-rich model bring the ice abundance of NH₂CHO to values up to ~ 5×10⁻⁹ at 14 K and canonical ζ (1.3×10⁻¹⁷ s⁻¹). We found a narrow maximum at 18 K and ζ=1×10⁻¹⁶ s⁻¹ with fractional abundances of 5×10⁻⁸, which indicate the ideal conditions for the production of this molecule. However, at the maximum, the T_d is rather close to the CO desorption temperature (Bisschop et al. 2006; Fuchs et al. 2006). The high reference binding energy of NH₂CHO, 5400 K in our model and in reality likely even higher, challenges the return of the molecule to the gas phase in prestellar cores. Therefore, we encourage non thermal desorption studies of NH₂CHO on CO ice to confirm the viability of releasing this molecule to the gas phase or search it JWST observations.

In the case of glycine, one of the interstellar molecules without positive detection that attracts more attention as the simplest amino acid, it has a low-temperature synthetic route on CO: $$\mathrm{H}\mathrm{O}\mathrm{C}\mathrm{O}+{\mathrm{C}\mathrm{H}}_2⟶{\mathrm{H}}_2\mathrm{C}\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}$$(R35) $${\mathrm{H}}_2\mathrm{C}\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}+\mathrm{N}\mathrm{H}⟶{\mathrm{N}\mathrm{H}\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}$$(R36) $${\mathrm{N}\mathrm{H}\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}+\mathrm{H}⟶{\mathrm{N}\mathrm{H}}_2{\mathrm{C}\mathrm{H}}_2\mathrm{C}\mathrm{O}\mathrm{O}\mathrm{H}$$(R37)

e.g. a three-step synthesis. In this case, the formation of the molecule is allowed through the concomitant diffusion of the reactants and the formation of HOCO (R10). Despite an increase in the formation of glycine by many orders of magnitude, the absolute abundance of glycine in ice remains too small for potential detection. However, we refer to Ioppolo et al. (2021) for alternative routes relying on nonthermal schemes that could potentially bring the abundances closer to detectable amounts.

Finally, we summarize and merge all the results discussed for ice abundances in Table 4 for a semiquantitative estimation of the effect of CO ice in the chemistry of the molecules considered in these sections. The table reinforces our general conclusion that COM abundances are enhanced in the models of ϵ_CO ≠1.0 compared with those of ϵ_CO=1.0 in most of the combinations of ζ and T_d. The enhancement, however, diminishes in the case of high ζ and low T_d.

Fig. 5 Same as Fig. 3 but for NH2CHO, HCOOCH3, OHCCHO, HOCH2CHO, and glycine, NH2CH2COOH.

3.4 Impact on gas phase abundances

We briefly discuss how CO ice affects the desorption of some molecules into the gas phase. Lower binding energies on CO-dominated ice surfaces can potentially enhance a series of thermal and nonthermal desorption mechanisms. From the three nonthermal mechanisms considered in our models, photodesorption is the only formally independent of the adsorbates’ binding energy, at least in our models. The remaining two, cosmic-ray-induced and chemical desorption, explicitly depend on the binding energy. The former is simulated in the Hasegawa-Herbst formulation (Hasegawa & Herbst 1993), which means that every grain is briefly heated stochastically to 70 K (10 μs) and the update in BE affects the desorption rates of molecules. With respect to chemical desorption, in Sec. 2 we discuss that our models use the same chemical desorption probability on H₂O ice and CO ice because the literature does not have experimental or theoretical data on the partition of the reaction energy into CO ices. Furthermore, the chemical desorption fraction (f, Eq. (3)) will also depend on the molecule under consideration (Fredon et al. 2021). Therefore, there are severe uncertainties in the probability of chemical desorption, and we prefer not to enlarge them in our models. However, we note that if f is calculated as a function of CO surface coverage, the differences in gaseous molecular abundances between models with ϵ_CO=1.0 and ϵ_CO ≠1.0 should be larger.

We study the effect on gas phase abundances of two small molecules (CH₄ and CO₂) and three COMs (C₂H₆, CH₃CHO and CH₃OH). Figure 6 shows the gas-phase abundances of these molecules at t=2× 10⁵ yr as a function of ζ and T_d. In contrast to ices at low temperatures, where molecules are stored and protected in the mantle, the gas is in a harsher environment where molecules do not reach a steady state, and are continuously formed and destroyed. Therefore, there is no unique choice of time to collect our data, and t=2×10⁵ yr does not reflect the ideal time for detecting the molecule, and the abundances presented here should not be considered as predictions, but rather as a comparison to show the differences in the models when desorption and diffusion rates vary with the composition of the ice surface. This is the reason, combined with the complexities in the gas-phase reaction network, why we refrain from providing a Table similar to Table 4 for gas-phase species.

The increased abundance of gas-phase molecules is hardly related to a single physical or chemical process; for example, gas-phase abundances are much more sensitive to destruction by gas-phase reactions as mentioned above. We first observe a smaller dust temperature dependence of gaseous molecular abundances to ice abundances, except for CO at high T_d. This is indirect evidence for the dominance of nonthermal effects in desorption at 8–18 K. For example, with low ζ, the gaseous CH₄ abundance at T_d < 18 K is higher in the model with ϵ_CO ≠1.0 due to the lower BE, which improves CR-induced desorption. However, grain contributions compete with the following gas-phase destruction reaction: $${\mathrm{C}\mathrm{H}}_4+{\mathrm{H}}_3^{+}⟶{\mathrm{C}\mathrm{H}}_5^{+}+{\mathrm{H}}_2\text{, }$$(R38)

at intermediate ζ. At high ζ, the higher abundance of CH₄ desorbed from ice (see Sec. 3.2) dominates over the above reaction. This is an example of the complex reasons behind the gas-phase molecular abundances. Similar arguments can be given for CO₂, where competing reactions are at play in the predicted abundances. The increased abundance of the gas phase CO₂ in models at low ζ and ϵ_CO ≠1.0 is not due to the desorption of CO₂ from the grain, but rather to the increased formation in the gas phase through: $${\mathrm{H}}_3^{+}+\mathrm{C}\mathrm{O}⟶{\mathrm{H}\mathrm{C}\mathrm{O}}^{+}+{\mathrm{H}}_2$$(R39) $${\mathrm{H}\mathrm{C}\mathrm{O}}^{+}+\mathrm{O}\mathrm{H}⟶{\mathrm{H}\mathrm{C}\mathrm{O}}_2^{+}+\mathrm{H}$$(R40) $${\mathrm{H}\mathrm{C}\mathrm{O}}_2^{+}+\mathrm{C}\mathrm{O}⟶{\mathrm{H}\mathrm{C}\mathrm{O}}^{+}+{\mathrm{C}\mathrm{O}}_2$$(R41) $${\mathrm{H}\mathrm{C}\mathrm{O}}_2^{+}+{\mathrm{e}}^{-}⟶\mathrm{H}+{\mathrm{C}\mathrm{O}}_2$$(R42) $${\mathrm{H}\mathrm{C}\mathrm{O}}_2^{+}+{\mathrm{N}\mathrm{H}}_3⟶{\mathrm{N}\mathrm{H}}_4^{+}+{\mathrm{C}\mathrm{O}}_2$$(R43)

so it depends indirectly on the amount of CO in the gas, which increases thanks to cosmic ray induced desorption at low and intermediate ζ (e.g., when hydrogenation of CO ice is slower), as shown above for CH₄. This highlights once again the complex and hardly predictable grain-gas interaction as well as the importance of chemical models in addressing this complex chemistry.

For most of the COMs, the main non-thermal desorption mechanism is chemical desorption, as explained above. For C₂H₆, cosmic ray-induced desorption and chemical desorption are on par, owing to the low binding energy of C₂H₆. In the case of CH₃CHO, whose ice abundance is most affected by the binding energy update on the CO ice (see Sec. 3.2), the explanation for the enhanced gas-phase abundance with ϵ_CO ≠1.0 models is straightforward. Chemical desorption is efficient because CH₃CHO has relatively low binding energy (Molpeceres et al. 2022; Ferrero et al. 2022), well represented by the Garrod (2013) reaction network (2500 K). The increased abundance on the ices presented in Sec. 3.2 directly correlates with the observed abundance in the gas, considering that our chemical desorption scheme does not change with ϵ_CO. Gas abundances of molecules with intermediate binding energies, for example, CH₃OCH₃, are also enhanced. However, we note that the abundances of CH₃CHO are low for its detection, but that is a consequence of the choice of time for the extraction of the data, as mentioned above. As an example, the global maximum in gas CH₃CHO abundances in our models corresponds to the model with ϵ_CO ≠1.00, T_d=10.7 K and ζ=1×10⁻¹⁸ s⁻¹, at t=3.2×10⁴ yr, with a fractional abundance of 1.2×10⁻¹⁰ with respect to H, e.g. a detectable amount consistent with prestellar core observations for this molecule (Jiménez-Serra et al. 2016; Megías et al. 2022).

Finally, it is worth indicating that our chemical models for ϵ_CO ≠1.00 do not show a notable increase in the abundance of CH₃OH in the gas phase (Fig. 6, bottom panel), as expected from the similar abundances of ice on the H₂O and CO ices reported in Sec. 3.2. An in-depth study of chemical desorption of CH₃OH from CO ices is desirable, to improve our description of chemical desorption and complement other desorption mechanisms available in the literature, e.g., photodesorption, co-desorption, and recently resonant infrared photodesorption; Öberg et al. (2009); Bertin et al. (2016); Martín-Doménech et al. (2016); Ligterink et al. (2018); Santos et al. (2023).

Fig. 6 Gas abundances at t=2×105 yr (see text for an explanation of the chemical timescales for gaseous molecules) in models of various ζ and Td. The selected molecules are CH4, CO2, C2H6, CH3CHO, and CH3OH.

4 Discussion and conclusions

Using state-of-the-art astrochemical models, we have investigated the role of reduced binding energies of adsorbates on CO-dominated ices in diffusion-mediated astrochemical reaction networks. The enhanced chemistry we found is attributed solely to the change in the intermolecular forces that drive diffusion, and new reactivity is not explicitly addressed (see, for example, the exclusive CO ice Eley-Rideal reaction, $\mathrm{C}+\mathrm{C}\mathrm{O}⟶\mathrm{C}\mathrm{C}\mathrm{O}$; Fedoseev et al. 2022; Ferrero et al. 2023). Diffusive chemistry is affected by CO ices, although non-diffusive mechanisms can also play a significant role in the promotion of chemistry in interstellar grains (Jin & Garrod 2020; Garrod et al. 2022). In fact,our models share some trends with the non-diffusive scenario, where Jin & Garrod (2020) report efficient production of CH₃CHO and CH₃OCH₃. Although our proposed routes differ, our findings are complementary for CH₃CHO and CH₃OCH₃. The non-diffusive mechanism relies upon the reactions CH₃* + CO and CH₃* + H₂CO (where the asterisk denotes a translational excited radical) for reactions to form CH₃CO and CH₃OCH₂ that can be easily hydrogenated to CH₃CHO and CH₃OCH₃.

Our chemical models reveal that the importance of CO-dominated ice is highly dependent on the physical conditions of the core, mainly through the dust temperature T_d and the cosmic ray ionization rate ζ. From these two, the critical parameter for enhanced chemistry on CO ices is ζ, which can vary several orders of magnitude (Dalgarno 2006; Padovani et al. 2009; Padovani et al. 2018b; Padovani et al. 2018a). For example, in the well studied prestellar core L1527, the cosmic ray ionization rate is likely to be closer to the canonical value of ζ=1.3×10⁻¹⁷ s⁻¹ (Redaelli et al. 2021) while in the giant molecular cloud G+0.693–0.027 in the center of our galaxy, ζ may present values between 100 and 1000 higher than the canonical value (Goto et al. 2013; Goto 2014). In general, COM formation is enhanced on CO-dominated ice with low ζ, while chemistry on CO-dominated ices and H₂O-dominated ices converges at very high ζ and low T_d. However, the effect at intermediate values of ζ is not trivial to predict, and astrochemical models remain instrumental in unraveling each specific case.

Our work comes at a good time for the astrochemistry canvas, with JWST observations indicating an important component of COM before forming a protostar (McClure et al. 2023). However, the definitive identification of COMs other than CH₃OH in ice observations has proven challenging Nazari et al. (2024). Our results align with the observed segregation of COMs in N-and O-bearing COMs presented in Jiménez-Serra et al. (2016); Jiménez-Serra et al. (2021); Megías et al. (2022) in the starless cores L1544, L1498, and L1517B². These observations show that O-bearing molecules and specifically CH₃CHO abundances are located toward the center of the core. On the contrary, N molecules are located at the edges of the cloud, making them more efficiently formed in earlier stages, possibly because of gas-phase nitrogen chemistry. Based on the evidence of efficient cyanation reactions (Balucani et al. 1999; Vazart et al. 2015; Puzzarini et al. 2020; Heitkämper et al. 2022) in the gas phase and the evidence collected in this work, the enhanced chemistry in CO appears to be more than plausible for the increased abundance of O-bearing COMs. This is in line with our models, where the impact of CO chemistry in the abundance of N-bearing COMs is less pronounced than that of O-bearing COMs. In the models presented in Sects. 3.2 and 3.3, we show abundances of CH₃NH₂, CH₃CN, and NH₂CHO, but this observation also applies to other molecules such as C₂H₃CN or C₂H₅CN. We note that the abundance of CH₃NH₂ is high in the grains (Fig. 3), but its abundance is also high in models that do not consider CO ice. Therefore, its chemistry on the CO ice is not as enhanced as in the case of O-bearing COMS. The reason for this chemical segregation is the importance of CO hydrogenation products in the promotion of COM formation. These products (HCO, CH₂OH, CH₃O) are precursors to O-bearing COM.

4.1 Required inputs to improve the models

It is evident from our simulations that CO-dominated ice promotes richer chemistry than the one present on H₂O ice, which is wholly dominated by hydrogenations at low temperatures, ≪20 K. Considering the relatively fast diffusion of radicals on CO ice, the diffusive constraints for several debated grain routes for the formation of COMs on ice are lifted, such as NH₂CHO (Noble et al. 2015; Song & Kästner 2016; Rimola et al. 2018; López-Sepulcre et al. 2019) or CH₃CHO (Enrique-Romero et al. 2016; Lamberts et al. 2019; Enrique-Romero et al. 2021). Despite models represent a theoretical effort to characterize the phenomenon on a qualitative level, physicochemical input from different sources can increase the predictive potential of the models, especially with regard to gas-phase species and comparison with observations.

The uncertainties in the rate equation model are mostly related to the E_diff/E_bin ratio (See Appendix C). Because CO chemistry occurs in a narrow range between 8–18 K, we can easily identify which radicals are vital in the production on COMs on CO ice, most notably CH₃ or HCO, depending on T_d (see Fig. 2). Should experiments or dedicated simulations provide better constraints on E_diff/E_bin for these species, a parameter that we encourage to search for, it will be easy to update the values derived here, seeking a better agreement between observations and models of prestellar cores.

Although less critical than the E_diff/E_bin ratio, two other factors could improve our models. First, to improve the description of chemical desorption and our predictions concerning the abundances of the gas phase, information on the energy transfer of chemical energy to a CO ice (compared to a H₂O ice) is required. Although intuition suggests that chemical desorption must be enhanced on CO ice, due to the lower BE of the adsorbates on the surface, it is not clear how chemical energy transfer occurs in these ices. In our models, energy dissipation is governed by a (Eq. (3)), and is assumed to be equal on CO ices and H₂O ices. The lack of knowledge about the differences in chemical desorption probabilities on different substrates makes us consider a similar desorption probability, something that can be easily improved when new investigations address this particular issue. The second input that could also slightly improve our models is the inclusion of exclusive CO chemistry, that is, reactions that can only occur on CO ice. To our knowledge, the only example to date involves the condensation of atomic carbon on CO ice (Fedoseev et al. 2022; Ferrero et al. 2023), a reaction whose modeling is the subject of a subsequent work. When new reactions are postulated, our models can easily incorporate them.

Finally, it is vital to indicate that our models’ initial conditions are tailored to produce CO-rich ice, allowing us to investigate the maximum impact of CO-dominated ice. Today, there is no clear consensus on the morphology of the CO/H₂O ice in interstellar environments, for example, whether it is layered (Boogert et al. 2015) or, on the contrary, the different components of the ice are mixed, forming islands on the core of the grain (Potapov et al. 2020). Recent JWST results support the presence of pure CO ice without significant mixing with H₂O ice (McClure et al. 2023). Even assuming the creation of layered ice, it is not straightforward to confirm whether interstellar CO completely wets the surface or not (Noble et al. 2012; Kouchi et al. 2021) or whether CO ice is crystalline or amorphous (Kouchi et al. 2021; He et al. 2021). These questions need to be answered by astronomical observations of interstellar ice. Our models coarse-grain all these unknowns via a continuous function of BE as a function of the CO surface coverage (Eq. (1)). We will continue monitoring the recent discoveries in interstellar ice structure through JWST observations and adapt our models accordingly.

4.2 Conclusions

We conclude this work with a series of bullet points that can be extracted from this work:

• 1.

  The diffusion of relatively heavy radicals on CO-ice (CH₃, HCO, N, O, CH₂ and NH) enhances alternative reactions to the hydrogenations on said ice. In particular, the formation of CH₃CHO, CH₃OCH₃, and C₂H₅OH appears to be significantly enhanced on CO ices;

• 2.

  The dust temperature contribution to the chemistry on CO-dominated ice is hardly predictable at first sight and needs to be considered on a per-reaction basis, although the competition between H atom diffusion and desorption versus other radicals diffusion is the main driving mechanisms for the observed changes;

• 3.

  The enhanced COM formation on CO-dominated ice is inhibited when ζ is high (normally >5×10⁻¹⁷ s⁻¹) and T_d is low. High ionization rates enhance H atom abundances and thus hydrogenation of radicals. Molecules that can be formed by hydrogenation of simpler precursors like CO, C₂N or C₂H₂ (yielding CH₃OH, CH₃CN, C₂H₄) do not see an appreciable increase in their abundances due to competition between heavy radical diffusion and H diffusion on CO ices;

• 4.

  There are two regimes for the enhancement of chemistry on CO ice. The first occurs for T_d ≤12 K, where we observe an enhancement of the chemistry initiated by CH₃, CH₂, O or N. A second regime for enhanced COMs formation is revealed above ≥14 K, initiated by the chemistry of HCO and NH. This second regime is less impactful due to the strong competition with hydrogenation. Still, the increase in abundance of COMs like NH₂CHO, HCOOCH₃, OHCCHO, HOCH₂CHO and NH₂CH₂COOH on CO dominated ices is notable, at low ζ;

• 5.

  The desorption of species from CO-dominated ices is also promoted under several conditions, such as an overall small binding energy of the desorbed molecule or a greater COM formation compared to the H₂O-dominated ices. However, extracting quantitative gas abundances becomes even more complicated when considering gas-phase chemistry in addition to the uncertainties in the model;

• 6.

  Our models can be enhanced once several unknowns are resolved from experiments or quantum chemical calculations. For example, we highlight the determination of E_diff/E_bin on CO ices in comparison with H₂O ices, the rate of chemical energy transfer for chemical desorption on CO ices, and the study of reactions taking place exclusively on CO ice.

The results presented in this paper serve as a step toward improving the description of the chemistry on interstellar dust grains and complement other contemporary advances, such as nondiffusive chemistry. We will continue improving our models, based on new observational evidence and theoretical advances. An example of these new techniques includes the recent approach presented in Garrod et al. (2022) of considering different binding sites on a substrate based on a probabilistic distribution of binding partners.

Our results are timely in the context of JWST observations. We will also continue expanding the quantum chemical dataset in search of binding energy outliers and adapting our methodology to simulate realistic pre-stellar cores, including the physical evolution of the objects in the astrochemical models.

Acknowledgements

We thank Dr. Joan Enrique-Romero for his insightful comments. G.M. thanks the Japan Society for the Promotion of Science (JSPS International Fellow P22013, and Grant-in-aid JP22F22013) for its support. G.M. also acknowledges support from the grant RYC2022-035442-I funded by MCIN/AEI/10.13039/501100011033 and ESF+. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no. INST 40/575-1 FUGG (JUSTUS 2 cluster) and the Research Center for Computational Science, Okazaki, Japan (Projects: 22-IMS-C301, 23-IMS-C128). K.F. acknowledges support from JSPS KAKENHI grant Numbers 20H05847 and 21K13967. Y.A. acknowledges support by Grant-in-Aid for Scientific Research (S) 18H05222 and Grant-in-Aid for Transformative Research Areas (A) grant Nos. 20H05847.

Appendix A Quantum chemical determination of ϵ_CO for CH₃, HCO, N, O, CH₂, NH and H₂

As discussed in the main text, our chemical models approximate BE on pure and mixed CO ice through the scaling of BE on water (Equation 2, ϵ). The reason for this is the absence of data on the former two with respect to the latter, due to a very high computational cost of computing binding energy distributions on realistic CO ice analogs. Furthermore, because BE is varied depending on θ_CO at each time step, we need a formulation applicable to intermediate situations that do not contain pure CO ices or H₂O ices (Equation 1). We use quantum-chemical calculations to derive ϵ_CO. For CO ice clusters, we follow a similar scheme as in (Ferrari et al. 2023), creating 5 clusters with 10 CO molecules each and sampling 8 adsorption positions for each cluster, for a total of 40 adsorption geometries per adsorbate+cluster combination. We chose to use several clusters and different adsorption geometries per cluster to account for a wide variety of binding sites for our reduced-size cluster models. For H₂O we repeat the same sampling for the same number of molecules. All cluster geometries are shown in Figure A.1. We then place the adsorbates (CH₃, HCO, N, O, CH₂, NH, and H₂) randomly in equispaced points on the sphere surrounding the clusters, projecting the adsorbates to ensure a minimum distance of ~3 Å between the center of mass of the adsorbate and the nearest neighbor on the cluster. Binding energies are then computed as the difference between the energy of the complex adsorbate+cluster and the sum of the separated components (adsorbate + cluster). The ωB97 family of functionals showed excellent performance in the study of binding energies, both on CO and H₂O ice (Ferrero et al. 2023; Hendrix et al. 2024; Ferrari et al. 2023). Therefore, for this work we used ωB97M-D4/def2-TZVP (Mardirossian & Head-Gordon 2016; Weigend & Ahlrichs 2005; Caldeweyher et al. 2019; Najibi & Goerigk 2020). All calculations shown in this Appendix are performed using ORCA 5.0.4 (Neese et al. 2020).

Zero-point energy corrections to the binding energy are introduced via a scaling factor, because of the high computational cost associated with deriving the molecular Hessian of every sampling point. We obtain the proportionality factor, α, using the same theoretical method (ωB97M-D4/def2-TZVP) and connecting non-ZPVE inclusive and ZPVE inclusive BE using a model system containing 4 molecules (either H₂O or CO) and an adsorbate test set consisting of CH₃, HCO, N, O, CH₂, NH, CH₂OH, CH₃OH, CH₄, CO, CO₂, H₂CO, H₂O, H₃CO, N₂, NH₂, NH₃, and OH. This approach has been successfully applied in computational works dedicated to the determination of binding energies of collections of molecules (Ferrero et al. 2020; Bovolenta et al. 2022; Perrero et al. 2022). The results of this study are shown in Figure A.2 and the models used for the ice substrate are shown in Figure A.1. From the graph, we observe a correlation between the non-ZPVE corrected BE (BE) and the ZPVE corrected, BE(0), with a correlation coefficient in line with other works in the literature. Then the binding energies used to derive ϵ_CO are obtained as: $$\mathrm{B}\mathrm{E}(0,{\mathrm{H}}_2\mathrm{O})=0.7921\times \mathrm{B}\mathrm{E}\left({\mathrm{H}}_2\mathrm{O}\right)$$(A.1) $$\mathrm{B}\mathrm{E}(0,\mathrm{C}\mathrm{O})=0.7183\times \mathrm{B}\mathrm{E}(\mathrm{C}\mathrm{O}).$$(A.2)

We observe a higher influence of ZPVE on CO ices than in H₂O ices, which results in a lower value of ϵ_CO, further facilitating diffusive chemistry. We note that even though in this appendix we differentiate between BE and BE(0), in the main text BE refers to BE(0) throughout.

Fig. A.1 Structures of the clusters used for the sampling of BE and BE(0) ZPVE corrections in this appendix. Further molecular geometries can be provided on request.

Finally, after calculating BE(0) using the 40 adsorption site sampling in the 10 CO/H₂O cluster and α, we can extract a representative binding energy as the average BE(0) considering the different adsorption sites. The results of this analysis, along with standard deviations, are collated in Figure A.3. It is worth mentioning that higher deviation of BE is found for H₂O clusters, due to the variability in H-bond topologies in the different binding sites. The BE(0) ratios shown in Figure A.3 are the ones presented in Table 1 of the main text.

The values in Table 1 of the main manuscript and Figure A.3 show that, in general, ϵ_CO lower than 0.40 describe well the binding of the adsorbates on CO to H₂O. This is in good agreement with Mondal et al. (2021), that reports 0.35–0.40 for molecules with permanent dipole moment or significant polarizability. Both our calculations and Mondal et al. (2021) show that molecules that form H-bonds, dipole-dipole, or effective dipole-induced dipole (as in CH₃) interactions will have ϵ_CO below 0.4. Among the molecules that we considered in Table 1, exceptions are N, CO, and H and H₂; they are apolar, small, and thus hardly polarizable. In the absence of a permanent dipole to form dipoledipole interactions, adsorbate surface interactions are driven by weaker forces, depending on the polarizability of the adsorbate. This allows us to rationalize the unity value found for H atoms (Fuchs et al. 2009; Hama et al. 2012; Kimura et al. 2018), because the interactions affecting the H atom, the smallest possible chemical adsorbate, are essentially the same on H₂O and CO ices. Likewise, a similar case can be found for the N atom, but here the species is larger, and a certain degree of polarizability is likely admitted. In fact, the phosphorous atom (P) Sil et al. (2021) finds an ϵ_CO of 0.53, using a monomer model, which follows the polarizability trend that we discussed in this paragraph. Not surprisingly, different P-bearing molecules have different ϵ_CO depending on the interactions at hand. Because P-bearing molecules are out of the scope of this modeling effort we refer the reader to the work of Sil et al. (2021). Continuing with the adsorbates with large ϵ_CO we have H₂, that according to Molpeceres & Kästner (2020) and Molpeceres et al. (2021b), should have ϵ_CO of around 0.5, and CO, that according to (Bisschop et al. 2006; Fuchs et al. 2009; Furuya & Persson 2018; Ferrari et al. 2023) should have values above 0.60 (0.66). The difference between these apolar, hardly polarizable, small molecules stems most likely from their quadrupole moments. Carbon monoxide having a larger quadrupole moment is reflected by a stronger binding with CO ice, in contrast with H₂ whose quadrupole moment is significantly lower.

Fig. A.2 Structures of the clusters used for the sampling of BE and BE(0) ZPVE corrections in this appendix. Further molecular geometries can be provided on request.

Fig. A.3 Average of the binding energies obtained using a ten molecule cluster for the key adsorbates considered in this work. In blue, binding energies on CO clusters, in orange, binding energies on H2O. Error bars represent standard deviation.

Appendix B Chemical models varying ϵ_CO for the majority of radicals

Our production models update all surface species not included in Table 1 with a constant scaling using ϵ_CO=0.66. While this selection is made to ensure that all species in the reaction network can have their binding energies updated, the selection of ϵ_CO=0.66 is based on the ratio found for CO (Furuya & Persson 2018, and references therein) and may be different for other species. To confirm that such a choice does not significantly alter our results, we ran additional models using a reduced grid for T_d and ζ (e.g. 25x25), considering ϵ_CO = 0.66 for all adsorbates not explicitly mentioned in Table 1 and ϵ_CO=1.00, e.g. not updating BE for those molecules. The results are visually shown for CH₃CHO, CH₃OH and HCOOCH₃ in Figure B.1 where the color scheme is selected differently from all the other heatmaps to highlight the comparison. A visual inspection of Figure B.1 shows that although there are differences between the models, these are only minor, as shown by the similar profiles in the heatmaps, especially compared with a model without a BE update (See Figures 3 and 5). These small changes demonstrate that the diffusion of the selected species in Table 1 is the main reason behind the chemistry enhancement on CO ice.

Fig. B.1 Comparison of abundances for a reduced grid of models considering an update of BE through ϵCO of 0.66 (labeled “UPDATE”) or 1.00 (labeled “NO UPDATE”) for all species labeled as “Rest of surface species” in 1. The color scheme is chosen deliberately different to the main text.

Appendix C The effect of E_diff/E_bin

As we discussed in the main text, the biggest unknown of our models is the diffusion vs. desorption ratio (E_diff/E_bin), which is the main parameter determining the mobility of chemical species on ices. To show how E_diff/E_bin affects our results, we run models for fixed T_d and ζ for different values E_diff/E_bin in comparison to a model without updating ϵ_CO and E_diff/E_bin=0.5. In the main text, the fiducial model considers E_diff/E_bin=0.4, but since in this appendix the maximum E_diff/E_bin is 0.5 the comparison needs to be held taking this value as a reference. The steady-state abundances for each of these conditions at 10, 12, and 14 K are shown in Figure C.1 for CH₃CHO, CH₃OCH₃ and C₂H₅OH, which are much affected by the update of BE in section 3.2. The effect of E_diff/E_bin shown in Figure C.1 is evident, especially at 10 K where the steady-state abundances can vary by several orders of magnitude. This effect is diminished at higher T_d. While better constraints on E_diff/E_bin on CO ices would surely allow for a more quantitative prediction of ice and gas phase abundances, it is clear from Figure C.1 that models that update the BE on CO ices yield systematically higher abundances than models without such an update (Shown by black lines in Figure C.1). Such a conclusion demonstrates that our hypothesis in this work, that is CO ices promote COMs chemistry, holds for different E_diff/E_bin values.

Fig. C.1 Steady state abundances (i.e.,1×106) yr for CH3CHO, CH3OCH3, and C2H5OH for fixed Td, and ζ=1.3×10−17 s−1. At each Td we vary the Ediff/Ebin between 0.3–0.5 for the models with updated BE on CO ice. For comparison, we also show the models for Ediff/Ebin=0.5 without the update of BE (labeled 0.5, H2O in the graphic’s legend).

References

All Tables

All Figures

	Fig. 1 Reference model. Top: evolution of the abundances of main ice species, surface H2O and CO, and gas CO in our reference model. Bottom: the fractional surface coverage of the main ice species corresponds to the fraction of molecules in the four topmost layers of a dust grain.
In the text

	Fig. 2 Initial analysis of diffusion timescales. Top: diffusion timescale as a function of the dust temperature (Td) for θ = 0.65 Bottom: as a function of the surface coverage of CO (θ) at 12 K for the main radicals considered (Furuya et al. 2022). The dotted black line and the shaded region around it show the H accretion timescale for different $n_{{\mathrm{H}}_{\text{atom }}}$ (see Sec. 3.1).
In the text

	Fig. 3 Ice abundances at t=2×105 yr in models of various ζ and Td. The selected molecules are CH3CHO, CH3OCH3, C2H5OH, CH3COOH, C2H6, and CH3NH2.
In the text

	Fig. 4 Same as Fig. 3 but for CH3OH, C2H4, CH4, and CH3CN.
In the text

	Fig. 5 Same as Fig. 3 but for NH2CHO, HCOOCH3, OHCCHO, HOCH2CHO, and glycine, NH2CH2COOH.
In the text

	Fig. 6 Gas abundances at t=2×105 yr (see text for an explanation of the chemical timescales for gaseous molecules) in models of various ζ and Td. The selected molecules are CH4, CO2, C2H6, CH3CHO, and CH3OH.
In the text

	Fig. A.1 Structures of the clusters used for the sampling of BE and BE(0) ZPVE corrections in this appendix. Further molecular geometries can be provided on request.
In the text

	Fig. A.2 Structures of the clusters used for the sampling of BE and BE(0) ZPVE corrections in this appendix. Further molecular geometries can be provided on request.
In the text

	Fig. A.3 Average of the binding energies obtained using a ten molecule cluster for the key adsorbates considered in this work. In blue, binding energies on CO clusters, in orange, binding energies on H2O. Error bars represent standard deviation.
In the text

	Fig. B.1 Comparison of abundances for a reduced grid of models considering an update of BE through ϵCO of 0.66 (labeled “UPDATE”) or 1.00 (labeled “NO UPDATE”) for all species labeled as “Rest of surface species” in 1. The color scheme is chosen deliberately different to the main text.
In the text

	Fig. C.1 Steady state abundances (i.e.,1×106) yr for CH3CHO, CH3OCH3, and C2H5OH for fixed Td, and ζ=1.3×10−17 s−1. At each Td we vary the Ediff/Ebin between 0.3–0.5 for the models with updated BE on CO ice. For comparison, we also show the models for Ediff/Ebin=0.5 without the update of BE (labeled 0.5, H2O in the graphic’s legend).
In the text