© The Authors 2025

1 Introduction

Currently, over 330 molecules¹ have been discovered in the interstellar medium, and their numbers are rapidly growing as observers make use of today’s outstanding telescope facilities in the radio, millimetre, and far-infrared regime. Around 47% of the detected molecules are considered to be large, consisting of six atoms or more. Interestingly, all of these larger molecules contain at least one carbon atom, which is why they are commonly called complex organic molecules (COMs; Herbst & van Dishoeck 2009) or interstellar COMs (iCOMs; Ceccarelli et al. 2017). The apparent dominance of carbon chemistry, both in the interstellar medium and on Earth, together with other chemical markers (e.g. D/H ratios), strongly suggest that the precursors of the building blocks of organic life on our planet were formed – and possibly conserved to some extent - in the earliest stages of star formation (Caselli & Ceccarelli 2012).

The consensus is that methanol is formed almost exclusively on the surface of dust grains by successive H-addition reactions to CO (e.g. Watanabe & Kouchi 2002; Fuchs et al. 2009; Santos et al. 2022). The formation scenarios of most other COMs, however, are still under debate. The present and ongoing discussion concerns the degree to which either gas-phase or grain-surface and ice mechanisms dominate the production of COMs in general. In the former case, this would occur as common grain mantle components (e.g. H₂O, CO₂, NH₃, CH₄, and CH₃OH; Boogert et al. 2015), which are formed mostly by the hydrogenation f atoms or simple molecules, are released into the gas phase by thermal and nonthermal desorption processes and then further processed by gas-phase reactions to form more complex molecules. In the case in which grain-surface chemistry is dominant, production could occur through diffusive or non-diffusive reactions, or through energetic processing (e.g. UV photolysis). The products would then be desorbed into the gas via thermal or nonthermal means, depending on the environment. Although these two possibilities are sometimes presented as either-or options for all COMs (e.g. Ceccarelli et al. 2023), it is apparent from astrochemical models that certain COMs have highly effective gas-phase production mechanisms, while others do not, with many still unexplored. Understanding the early formation mechanisms of methanol is extremely relevant in both formation scenarios, as it plays an important role in both of them. In a scenario where grain surface mechanisms dominate, COMs are usually assumed to form by radical-radical or radical-neutral reactions in which HCO, CH₃O, and CH₂OH, as intermediate products of the methanol formation route, play a fundamental role. In contrast, in the gas-dominated scenario, methanol is desorbed into the gas phase subsequent to its formation on the grain surface and acts as an important precursor for larger COMs.

In the case of grain surface formation, COMs were formerly believed to form only in a lukewarm temperature regime (30–100 K) during a gradual warm-up phase of the molecular clouds. In this warm-up scenario, grain surface radicals are first produced by the irradiation of UV photons and cosmic rays. Once the grain temperature increases to about 30 K, the radicals become more mobile, diffuse over the grain surface, meet each other, and react (Garrod & Herbst 2006; Garrod et al. 2008). However, detections of multiple COMs in even earlier, colder (< 15 K) stages of star formation (e.g. Öberg et al. 2010; Bacmann et al. 2012; Cernicharo et al. 2012; Jiménez-Serra et al. 2016 ) bring doubt regarding the warm-up scenario as the sole formation process. Alternatively, or in addition to the warm-up scenario, the introduction of non-diffusive reaction mechanisms has been proposed in various works. Hereby, the term non-diffusive chemistry consists of various reaction mechanisms, which have in common that they do not require the diffusion of species heavier than commonly present atoms (e.g. H, O, C, and N). Instead, they rely on the (small) probability that both reaction partners are formed in close proximity to each other. At the low temperatures of pre-stellar cores (6 K ≤ T_grain ≤ 15 K), non-diffusive reaction mechanisms could be the most relevant process for the formation of COMs, as only the lightest species are expected to diffuse over the grain surface, likely excluding a diffusion-driven COM formation scenario for these conditions. Garrod & Pauly (2011) introduced a single three-body reaction, very similar to the one used in the present work, between the reactants H, O, and CO to explain the formation of CO₂ on the surface of dust grains under dark-cloud conditions. Ruaud et al. (2015) showed that Eley-Rideal reactions in combination with a complex induced reaction mechanism of gas phase carbon atoms with the grain surface are able to increase the formation rate of multiple O-bearing COMs. Chang & Herbst (2016) added a chain-reaction mechanism to their unified microscopic-macroscopic Monte Carlo code. Conceptually, their mechanism is very similar to the three-body reactions used in this work; however, it is not straightforwardly applicable to a rate-equation based code. Similar to Garrod & Pauly (2011), Dulieu et al. (2019) introduced a single three-body reaction, between H₂NO and H₂CO to reproduce the formation of formamide (NH₂CHO) in co-deposition experiments of H, NO, and H₂CO at 10 K. Their approach tracks two types of H₂NO molecules: the ones that are formed next to an active H₂CO molecule and the ones that are formed next to other molecules. Although different from the one used in this work, their treatment should produce the same results. Shingledecker et al. (2019) showed that radiolysis can efficiently form radicals and other reactive species in the grain’s bulk, which are then able to react non-diffusively with neighbouring molecules and drive a rich chemistry there. The most extensive introduction and the one adopted in the present work, however, is presented in Jin & Garrod (2020) and Garrod et al. (2022).

The importance of reactions between heavy reaction partners under dark cloud conditions is supported by several co-deposition experiments of simple and abundant molecules. In experiments with CO molecules and H atoms at 13 K, Fedoseev et al. (2015) found the formation of glycoalde-hyde (HC(O)CH₂OH) and ethylene glycol (H₂C(OH)CH₂OH) in addition to formaldehyde (H₂CO) and methanol (CH₃OH) formation. The result was confirmed and extended to methyl formate (HC(O)OCH₃) formation in CO+H+H₂CO and CO+H+CH₃OH co-deposition experiments at 15 K by Chuang et al. (2016). There, H₂CO and CH₃OH were only deposited simultaneously to artificially enhance their abundances to provide sufficiently high methyl formate yields to meet the detection sensitivity limits of the experimental setup. In a similar manner, Fedoseev et al. (2017) detected glycerol (HOCH₂CH(OH)CH₂OH and tentatively detected glyceraldehyde (HOCH₂CH(OH)CHO) upon co-deposition of CO+H+HOCH₂ CHO at 13 K. Qasim et al. (2019) showed that the formation of propanal (H₃CCH₂CHO) in co-deposition of C₂H₂+CO+H at 10 K and its subsequent hydrogenation to 1-propanol (CH₃CH₂CH₂OH) is possible. Moreover, the formation of the simplest amino acid, glycine (NH₂CH₂COOH) has been shown in co-deposition of CO+NH₂CH₃+O₂+H at 13 K by Ioppolo et al. (2021). Their proposed formation routes are supported by astrochemical models, performed with MAGICKAL (Garrod 2013), including non-diffusive reaction mechanisms. Fedoseev et al. (2022) find that ketene (CH₂CO) is formed in co-deposition of H₂O + C + CO + H at 10 K with possible further hydrogenation steps leading to acetaldehyde (CH₃CHO) and ethanol (CH₃CH₂OH) formation. However, the most relevant finding for the present work are the co-deposition experiments of H₂CO+H performed in a temperature range from 10 K to 16 K by Santos et al. (2022). The authors confirmed the predominant formation of methanol by the radical-molecule reaction H₂CO + CH₃O → CH₃OH + HCO over the simple H-addition reaction, which was formerly suggested by microscopic kinetic Monte Carlo simulations (Simons et al. 2020. What all the proposed formation pathways have in common is that they involve at least one reaction step where larger surface radicals or molecules react with one another, possibly preceded and-or followed by multiple hydrogenation steps. Before, such reaction steps were deemed inefficient, as both reaction partners are considered to be immobile on the surface of dust grains in the extremely cold conditions of pre-stellar cores. Their inclusion makes invoking non-diffusive reaction mechanisms into astro-chemical models necessary. Otherwise, the efficiency of the proposed reactions could be severely underestimated. Interestingly, the found pathways do not require that the grain mantle products are energetically processed by UV photons or cosmic rays; most of the radicals are instead formed by simple surface hydrogenation addition and abstraction reactions.

The inclusion of non-diffusive reaction mechanisms into rate-equation based astrochemical codes is still under scrutiny. While microscopic Monte Carlo codes account for the exact position of each molecule at all times and can therefore accurately determine the amount of respective reactants sitting in close proximity to each other, rate-equation based codes compute only the average formation and destruction rates of the molecules included in their chemical network. Unfortunately, this large amount of information comes with a high computational cost, and hence, studies conducted with a microscopic Monte Carlo code usually use chemical networks including only a small number of species and reactions. To obtain a broad overview over the chemistry occurring in any interstellar region, rate-equation based codes are typically better suited accounting for the multitude of different species and reactions present. In the formulation of non-diffusive chemistry proposed by Jin & Garrod (2020) and Garrod et al. (2022), the probability of the second reactant being close by when the first one forms is simply estimated to be N_i/N_S. Here, N_i is the average number of atoms or molecules of species A or B on the grains, and N_S is the total number of binding sites on a grain. While this recipe is simple to integrate into a rate-equation based code, it can lead to an overestimation of reactants being formed close to each other in the case only very small numbers of the respective reactants are present on the grain. Additionally, rate-equation based codes are unable to account for clustering of surface species, which could potentially make certain non-diffusive reactions more efficient than would otherwise be the case for a random distribution of reactants on the grain surface.

Here, we present the inclusion of non-diffusive reaction mechanisms following the formulation by Jin & Garrod (2020) and Garrod et al. (2022) into an astrochemical code capable of simulating deuterium chemistry (pyRate; Sipilä et al. 2015a, 2019b). The addition of non-diffusive chemistry to pyRate provides valuable additional information about the deuteration of COMs that could formerly not be obtained. We applied the added reaction mechanisms to methanol, the simplest of the O-bearing COMs, as it is a key species in the two major formation scenarios for COMs discussed in the current literature. To connect to our previous works (Vasyunin et al. 2017; Chacón-Tanarro et al. 2019; Riedel et al. 2023), all the necessary modifications were performed in a step-by-step manner to enable the reader to separately evaluate the effects of each step. In a subsequent work, we will extend our chemical network to larger species and investigate the formation and deuteration of other COMs.

This paper is structured as follows: in Section 2 we describe the physical model and the updated chemical model in detail. In Section 3 we add various non-diffusive reaction mechanisms to the base model individually and in combination with each other. We discuss their importance for the formation and deuteration of methanol and compare them to more conventional modelling approaches and the observationally obtained column density and deuterium fraction profiles from Chacón-Tanarro et al. (2019). Section 4 discusses additional modifications to the chosen chemical and physical parameters. Section 5 presents our conclusions. Appendices A and B provide more detailed information for some minor intermediate modification steps.

2 Methods

2.1 Introducing non-diffusive reaction mechanisms

PyRate was updated to include the non-diffusive mechanisms presented in Jin & Garrod (2020) and Garrod et al. (2022). In general, the present work follows their formulation closely, and we refer to these works for an extensive explanation of the adopted mechanisms and their respective formulation. However, for some of our purposes, namely to ensure comparability to our previous work (Riedel et al. 2023), it proved necessary to modify some equations slightly. Therefore, we lay out very briefly relevant equations below to explain these minor modifications. In Section 3.2, following the comparison to previous results, we switch to the original formulation of Jin & Garrod (2020) and Garrod et al. (2022) and use it for the remainder of the paper to ensure comparability to the literature.

According to the standard formulation, as detailed in Hasegawa et al. (1992), the reaction rate R_AB of reactants A and B in chemical reactions proceeding diffusively, is given by the following expression: $$\begin{array}{c}{R}_{AB}=f_{\text{act}}\left(AB\right)\left[k_{\text{hop}}\left(A\right)N\left(A\right)\right]\frac{N\left(B\right)}{N_{\text{S}}}\\ +f_{\text{act}}\left(AB\right)\left[k_{\text{hop}}\left(B\right)N\left(B\right)\right]\frac{N\left(A\right)}{N_{\text{S}}},\end{array}$$(1)

where k_hop(A) or k_hop(B) are the rates of thermal hopping (occasionally modified to account for diffusion by quantum tunneling) of molecule A or B, respectively, N(A) or N(B) the abundance of species A or B on the grain (here: average number of molecules on individual grain) and N_S total number of binding sites on the grain. The quantity ƒ_act denotes the efficiency of a chemical reaction and assumes a value between 0 and 1.

Reaction rates R_AB for non-diffusive reactions, as detailed in Jin & Garrod (2020) and Garrod et al. (2022), follow a similar structure as Eq. (1): $$\begin{array}{c}{R}_{AB}=f_{\text{act}}\left(AB\right)R_{\text{comp}}\left(A\right)\frac{N\left(B\right)}{N_{\text{S}}}\\ +f_{\text{act}}\left(AB\right)R_{\text{comp}}\left(B\right)\frac{N\left(A\right)}{N_{\text{S}}},\end{array}$$(2)

where the diffusion rate [k_hop(A)N(A)] of reactant A or the diffusion rate [k_hop(B)N(B)] of reactant B is replaced by the completion rate R_comp(A) or R_comp(B) that describe the rate of ‘appearance’ of reactant A and subsequent reaction with reactant B (or vice versa). The completion rates are derived by Eqs. (3) and (4): $$R_{\text{comp}}\left(A\right)=\frac{1}{1/R_{\text{app}}\left(A\right)+t_{\text{AB}}}$$(3) $$R_{\text{comp}}\left(B\right)=\frac{1}{1/R_{\text{app}}\left(B\right)+t_{\text{AB}}},$$(4)

where R_app(A) and R_app(B) are the appearance rates, that depend on the specific non-diffusive mechanism used, and t_AB is the average lifetime of both reaction partners staying in a state in which they are ready to react.

Non-diffusive mechanisms are always allowed to act on the entire chemical network and in addition to the well-established diffusive reactions. The added mechanisms include the Eley–Rideal mechanism, photodissociation-induced mechanisms and three-body reactions. Eley–Rideal reactions are reactions in which a molecule from the gas phase adsorbs onto the grain surface and immediately finds a reaction partner in the vicinity of its adsorption spot. The appearance rate for a reactant is simply equal to the adsorption rate of this species. Photodissociation-induced reactions are initiated by a photon dissociating a larger molecule on the dust grain surface into two fragments, followed by a reaction of one of those fragments with another surface molecule in close proximity. The species appearance rate is calculated as the sum of photodissociation rates producing this species. Here, we adopt the description in Jin & Garrod (2020) (or the PDI version in Garrod et al. 2022), as the proposed refinement that considers that unsuccessful immobile photo-products recombine immediately is not of importance for surface species, since they can diffuse away from their production site before they would recombine again. In three-body reactions, a molecule is initially formed either diffusively, in a reaction with at least one mobile reaction partner (e.g. H), or via the above described non-diffusive reaction mechanisms. Subsequently, the newly formed reaction product can then spontaneously react with another adjacent molecule. Here, the appearance rate is the sum of the reaction rates of all diffusive reactions producing this species. In case they are switched on, the summed reaction rates of Eley–Rideal and photodissociation-induced reactions are also added to the appearance rate. Additionally, the three-body reaction mechanism has two variations. It is possible to perform several successive rounds of three-body reactions. While the first round can be initiated by diffusive reactions, and possibly Eley-Rideal and photodissociation-induced reactions, the following rounds are initiated by the respective preceding round of three-body reactions. The cost of running multiple rounds of three-body reactions with a chemical network comprising 75 000 reactions is extremely large and the effect expected to diminish with each round, which is why we limit ourselves to run models with two rounds of three-body reactions. Another option, three-body reactions with excited formation, is to assume that the reaction enthalpy released by the initiating reaction can (partially) help to overcome the activation-energy barrier of the follow-up reaction. This mechanism is only sensible in combination with the basic three-body reactions.

The above described non-diffusive reaction mechanisms can, in principle, be applied in arbitrary combination with each other. However, as running every possible combination would result in an unreasonable amount of models that also not necessarily provide more insights, we restrict ourselves to six models, involving non-diffusive chemistry. We perform a model run, where Eley–Rideal reactions (ND1/ND7), photodissociation-induced reactions (ND2/ND8) and three-body reactions (ND3/ND9) each are enabled separately. Moreover, we test a model with two rounds of three-body reactions (ND4/ND10) and a model where three-body reactions with excited formation are enabled alongside the basic version of three-body reactions (ND5/ND11). At last, we also run a model, where we combine Eley–Rideal, photodissociation-induced and three-body reactions (ND6/ND12), but none of the variations.

2.2 Chemical model

The chemical evolution of molecular abundances is calculated with the gas-grain astrochemical code pyRate, capable of simulating deuterium chemistry and described in more detail in Sipilä et al. (2015a, 2019b). The overall chemical network is based on the 2014 public release of the Kinetic Database for Astro-chemistry (KIDA) gas-phase network (kida.uva.2014, Wakelam et al. 2015). Respective reactions were cloned to include deuterium chemistry for molecules with up to seven atoms and spin-state chemistry for H₂, ${\text{H}}_2^{+}$, and ${\text{H}}_3^{+}$, as well as multiply pro-tonated species involved in the water and ammonia formation networks and their respective deuterated isotopologues. More details can be found in Sipilä et al. (2015a,b) and Sipilä et al. (2019b). In total the network includes a total of ≈75000 gas-phase and grain surface reactions, making it substantially larger than commonly used chemical networks that do not treat the chemistry of isotopologues explicitly. The methanol formation network is the same as presented in detail in Riedel et al. (2023). If not otherwise indicated, a network without abstraction reactions is adopted. However, in Section 3.3, the experimentally verified (Santos et al. 2022) methanol formation pathway H₂CO + CH₃O → CH₃OH + HCO is added to the network. Moreover, in Section 4.2, we also present and discuss selected models including the abstraction scheme proposed by Hidaka et al. (2009).

Surface reactions are allowed to proceed through the Langmuir-Hinshelwood mechanism relying purely on thermal diffusion. Diffusion by quantum tunneling of hydrogen and deuterium atoms is not enabled in most models. We show, however, some models including both types of diffusion for the purpose of comparison to previous works (Vasyunin et al. 2017; Chacón-Tanarro et al. 2019; Riedel et al. 2023) in Section 3.1. The diffusion-to-binding energy E_d/E_b is set to 0.55, a value suggested for atoms by Minissale et al. (2016a) based on their measurements of diffusion and desorption energies of O and N atoms.

Dust grains are assumed to be spherical, with a radius of 0.1 µm and a surface site density of 1.5 × 10¹⁵ cm⁻². We employed a three-phase grain model, including a gas phase, a chemically active surface phase and an inert mantle phase. In contrast to Jin & Garrod (2020), diffusive and non-diffusive mechanisms are carried out solely on the surface of the dust grains. Chemical reactions or any form of bulk diffusion is not included. The impact of bulk diffusion and reactions is probably small in the physical conditions of pre-stellar cores, where the grains undergo net adsorption rather than net desorption and as long as no invasive desorption mechanisms, that partially (or entirely) destroy the grain mantle, are applied. It can, however, become important if one considers the effect of cosmic-ray-driven radiation chemistry in combination with non-diffusive bulk reactions as was done in the works by Shingledecker et al. (2019), Shingledecker et al. (2020). This is beyond the scope of the present paper.

Initially, we adopted atomic initial abundances (see Table 1) taken from Semenov et al. (2010), as they were used for previous works (Chacón-Tanarro et al. 2019; Riedel et al. 2023). This choice is revisited in Section 4.3.

Besides thermal desorption, pyRate contains several non-thermal desorption mechanisms. In this work, we only apply cosmic-ray induced desorption and reactive desorption as non-thermal desorption mechanisms, since those can efficiently desorb methanol from the grain surface with reactive desorption presumed to be the dominant one. Photo-desorption of methanol has been shown experimentally to have negligible impact on the release of methanol from the grain surface (Bertin et al. 2016; Cruz-Diaz et al. 2016). PyRate contains two different options for the treatment of reactive desorption. The first option is a basic mechanism where the reaction products desorb with an efficiency of a specified, but constant value. The second option derives an individual reactive desorption efficiency for every product species that depends on the reaction enthalpy and type of underlying surface Minissale et al. (2016b) and was previously tested in Riedel et al. (2023). Most presented models apply the first option with an efficiency of 1%, as this provides a better comparability to the literature. However in Section 4.2, we also present the effects of the experiment-based reactive desorption mechanism for a small selection of models.

During the test phase, we found that when using the Eley–Rideal mechanism, the code will not converge to a solution, unless there is some method in place to get rid of an unphysical buildup of H₂ on the surface of the dust grains. In Section 4.1, we test various options to lower the H₂ coverage, including scaling the binding energies of H₂ and its isotopologues by a factor of 0.1 and the modified binding energy treatment presented in Garrod & Pauly (2011) applied to all species or solely to H, H₂ and its isotopologues. The overall effects on the column densities of CH₃OH are small, which is why we have used the simple scaling option on all models presented in the main body of the paper.

Many chemical reactions both in the gas phase and on the surface of dust grains are exothermic but still require the reactants to overcome activation-energy barriers to react with each other. The formation scheme of methanol comprises at least two hydrogenation steps possessing considerable barriers (CO → HCO; E_A = 1.76 × 10³ K, H₂CO → CH₃O/CH₂OH; 2.00 × 10³ K/5.16 × 10³ K). However, the extremely cold environments of pre-stellar cores pose a significant obstacle for those reactions to proceed, as the reactants typically have only very low energies. Hence, the details of calculating the efficiency of reactions with activation-energy barriers are crucial for the amount of methanol (and other COMs) produced.

In previous models (see Chacón-Tanarro et al. 2019; Riedel et al. 2023), a single collision reaction probability (hereafter: SC), which was proposed in Hasegawa et al. (1992), was applied. It derives the probability κ for the reaction to happen upon encounter of the reactants as either a Boltzmann factor, κ = exp(−E_A/T_d), where E_A is the activation energy of the reaction and T_d the dust temperature, or a tunneling probability, $\kappa =\exp \left(-2\left(a/ħ\right)\sqrt{2\mu k_{\text{B}}{E}_{\text{A}}}\right)$, where a is the thickness of a rect-angular barrier and µ the reduced mass. In Riedel et al. 2023 and the present work, reactive tunneling with a barrier thickness of 1 Å is used.

However, the formulation of non-diffusive mechanisms proposed in Jin & Garrod (2020) assumes that rates for reactions containing an activation-energy barrier are calculated using reaction probabilities derived from the reaction-diffusion competition (RDC) model proposed by Chang et al. (2007) (see also Garrod & Pauly 2011). It considers that the reaction partners are confined in the same binding site until one of them diffuses away again. Therefore, they have multiple opportunities to react with each other as opposed to just one in the SC model. In the RDC model, a reaction with an activation-energy barrier is typically significantly more likely to happen as compared to the SC model, as the reactants can undergo multiple attempts to react. Therefore, reactants have to meet less frequently to produce the same amount of successful reactions. In recent years, many modellers adopted the RDC model into their chemical models, as it can provide higher gas phase abundances without assuming a highly efficient diffusion process on the surface or an unlikely high desorption rate.

First, in Section 3.1, to present a consistent approach for both diffusive and non-diffusive reactions and to be fully able to compare the effects of non-diffusive chemistry to the previous models, we modified the equations presented in Jin & Garrod (2020) and Garrod et al. (2022) accordingly, which can be simply done by setting ƒ_act to κ. We, also, removed the t_AB-term from Eqs. (3) and (4), as the SC model assumes that the reaction either happens on first encounter of the reactants or not at all. In Section 3.2, we switch from using the SC model to using the RDC model for the remainder of the paper. For that, the efficiency factor ƒ_act is adjusted to be $$f_{\text{act}}\left(AB\right)=\frac{v_{AB}{\kappa }_{AB}}{v_{AB}{\kappa }_{AB}+k_{\text{hop}}\left(A\right)+k_{\text{hop}}\left(B\right)},$$(5)

where κ_AB is, as described above, either a Boltzmann factor or a tunneling probability. k_hop(A) and k_hop(B) are the rates of thermal hopping of species A or B, respectively. ν_AB denotes the frequency of collision of both reactants, which is here taken to be the larger of the characteristic frequencies of the reactants derived by the respective equation that is presented in Hasegawa et al. (1992).

2.3 Physical model

We present a set of static models using a one dimensional physical model (see Figure 1) for H₂ density, gas temperature T_gas, dust temperature T_d and visual extinction A_V, which was derived from the one presented in Keto & Caselli (2010), with a more detailed description in Sipilä et al. (2019a). The chemical evolution is solved separately for 35 concentric shells spanning the entire core radius of 0.32 pc, resulting in spherical symmetric spatio-temporal evolution of molecular abundances.

The column density distribution is calculated by an integration along the line of sight for different impact factors and is subsequently convolved with a 30″ Gaussian beam to be comparable to the observations by Chacón-Tanarro et al. (2019).

3 Results

In the following, various modifications are introduced in a step-by-step manner to make it possible to follow their effects separately. All presented models are listed in the order of their appearance in Table 2. Modified physical and chemical parameters are indicated as well.

3.1 Applying a single collision model

The upper panel of Figure 2 presents the column density profiles of CH₃OH (left) and N(CH₂DOH)/N(CH₃OH) profiles (right) of various models including non-diffusive reaction mechanisms at the best-fit time of t = 3.0 × 10⁵ yr. The lower panel sets these results in comparison to previous results (Riedel et al. 2023) and the available gas phase observations by Chacón-Tanarro et al. (2019). The single-dish observations, carried out with the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope, show that CH₃OH peaks in an asymmetric ring (see also Bizzocchi et al. 2014; Vastel et al. 2014; Spezzano et al. 2016). The position of strongest emission (methanol peak) is located in the northeastern part of the ring offset from the position of the dust peak. The best-fit time is taken to be the time-step, where the χ² value of the observed central column density versus the corresponding value in the base model for H₂CO and CH₃OH (and their observed deuterated isotopologues) are minimised.

A purely diffusive model with a diffusion-to-binding energy ratio of E_d/E_b of 0.55 (slow thermal diffusion; model D1) is always shown as a reference. We find that applying the non-diffusive mechanisms as proposed by Jin & Garrod (2020) and Garrod et al. (2022) and modified to be consistent with the SC model as detailed above, is only adding minor contributions to the gas phase methanol reservoir, independent of the respective mechanism. While the largest impact on the CH₃OH and N(CH₂DOH)/N(CH₃OH) profiles is seen when adding three-body reactions (model ND3), a further addition of a second round of three-body reactions (ND4) or three-body reactions with excited formation (model ND5) is not increasing the methanol yield much more. Also, Eley–Rideal reactions (model ND1) only cause a slight increase in methanol production.

Photodissociation-induced reactions (model ND2) have a negligible impact on the methanol reservoir. The increase for CH₃OH is typically larger than for its deuterated counterpart CH₂DOH, consequently reducing the N(CH₂DOH)/N(CH₃OH) slightly. Most of the increase in methanol can be accredited to opening up an additional pathway for simple H-addition reactions. This is particularly the case for the ND1 model, where the Eley–Rideal mechanism contributes up to 44% at the dust peak (densest part) and 13% at the methanol peak. Reactions between heavier reactants (e.g. OH + CH₃ → CH₃OH or OH + CH₂D→CH₂DOH), which occur predominantly via non-diffusive reaction mechanisms, start to play a role as well and are even the dominant formation route on occasion. Especially for the formation of the methanol precursors – CH₃O and CH₂OH – the reactions O + CH₃ → CH₃O and OH+CH₂ → CH₂OH dominate the formation (for CH₃O) or are comparable to the hydrogenation reaction (for CH₂OH). These formation routes are likely favoured as they do not possess an activation-energy barrier, while the hydrogenation reactions H+H₂CO→CH₃O (E_A = 2 × 10³ K) and H+H₂CO→CH₂OH (E_A = 5.16 × 10³ K) have a substantial one. For Eley–Rideal reactions, we typically find an increased importance to the innermost part of the core, where the higher gas phase abundances result in larger adsorption rates. The photodissociation-induced and three-body reactions, however, show higher contributions at the position of the methanol peak.

Riedel et al. (2023), in agreement with Vasyunin et al. (2017) and Chacón-Tanarro et al. (2019), conclude that an increase of the diffusion rates of hydrogen and deuterium atoms over the grain is necessary to explain the observations. They tested two available options: either allowing the diffusion to take place both by a slow thermal hopping process (E_d/E_b = 0.55) and additional quantum tunneling of hydrogen and deuterium atoms through a rectangular barrier with a width of 1 Å or by a fast thermal hopping process (E_d/E_b = 0.2). Additionally, a reference model with solely a slow thermal hopping process (E_d/E_b = 0.55) was performed. Both options with increased diffusion rates are able to increase the methanol column densities by several orders of magnitude compared to the reference model, but they still underproduce the observed methanol column densities by roughly an order of magnitude. However, most models presented in Vasyunin et al. (2017), Chacón-Tanarro et al. (2019), and Riedel et al. (2023) apply an experiment-based formulation of reactive desorption proposed by Minissale et al. (2016b), which derives an individual reactive desorption efficiency for every product species in a chemical reaction, depending on the reaction enthalpy and type of underlying surface. Since this version of reactive desorption typically produces reactive desorption efficiencies much lower than 1% for heavier reaction products, this underproduction could potentially be caused by an underestimation of the desorption process. For a full comparison of the models using a form of increased diffusion rate with the models comprising non-diffusive reaction mechanisms presented here, we performed the former ones again, but with scaled H₂ binding energies and a constant reactive desorption efficiency of 1% (see Appendices A and B for more details).

In the framework of the SC model, we conclude, when using solely a slow diffusion process (E_d/E_b = 0.55), the diffusion rate of hydrogen (and deuterium) atoms is not efficient enough. Applying the non-diffusive mechanisms is only adding a minor contribution to the gas phase methanol reservoir. Models increasing the diffusion rate, especially of hydrogen (and deuterium) atoms, tested in the previous works, are more successful in that regard and are able, when using a constant reactive desorption efficiency of 1%, to increase column densities to almost the observed order of magnitude and simultaneously produce a N(CH₂DOH)/N(CH₃OH) profile, which is well within the area of uncertainty of the observations.

However, available laboratory experiments (Hama et al. 2012 on an ASW surface, Kimura et al. 2018 on a CO surface) found that the diffusion of hydrogen (and deuterium) atoms at 10 K is not dominated by quantum tunneling but rather thermal hopping, as a significant isotope effect between hydrogen and deuterium atoms has not been observed. These results are also backed up by theoretical work (Senevirathne et al. 2017).

Also, measurements for diffusion energies E_d or more commonly diffusion-to-binding energy ratios E_d/E_b are quite sparse. They are found to be species-dependent with E_d/E_b = 0.2 (as used for the model with a fast thermal hopping process) being the lowest value debated in the literature (Furuya et al. 2022). However, these values were measured for CO and CH₄ and not for H. Moreover, experiments (Watanabe et al. 2010; Hama et al. 2012) conclude that potential wells on the surface of the dust grain can be divided into groups according to their depth, with Hama et al. (2012) identifying at least three different ones. Most astrochemical codes are not able to model such complicated grain surface properties, with only a few exceptions (e.g. Grassi et al. 2020), but only use a uniform diffusion-to-binding energy ratio. Additionally, Jin & Garrod (2020) find that a low value for the diffusion-to-binding energy ratio (E_d/E_b = 0.35 in their case) causes a quicker hydrogenation of radicals, which shortens their lifetimes and consequently could make proposed surface formation pathways for COMs less efficient.

Fig. 1 Physical model first presented in Keto & Caselli (2010) and described in more detail in Sipilä et al. (2019a) providing static radial profiles of the H2 volume density, n(H2) (blue, logarithmic scale), the visual extinction, AV (red), Tgas (orange), and the dust temperature, Tdust (green).

3.2 Applying a reaction-diffusion competition model

Figure 3 presents the column density profiles of CH₃OH (left) and N(CH₂DOH)/N(CH₃OH) profiles (right) of the diffusive model with a slow thermal hopping process (E_d/E_b = 0.55; D4) and various non-diffusive models (ND7-ND12), derived by assuming the RDC model. They are presented at the best-fit time of t = 3.0 × 10⁵ yr. We find that the column densities of CH₃OH and CH₂DOH in model D4 are increased by several orders of magnitude compared to the respective model derived with the SC model (D1) due to the higher reaction probabilities of the reactions with activation-energy barriers (CO + H → HCO and CH₂OH/CH₃O + H → CH₃OH). In fact, the CH₃OH column density profile matches the observationally constrained profiles within a factor of around two, which is considered to be a good agreement for astrochemical models.

Interestingly, the N(CH₂DOH)/N(CH₃OH) in model D4 is decreased by a factor of two to three compared to model D1. This decrease is caused by the change in the reaction probabilities for important formation pathways of CH₂DOH when applying the RDC model. CH₂DOH can be formed via three different pathways, starting from H₂CO: (I) the hydrogenation of H₂CO to CH₂OH and its subsequent deuteration (see reaction (6)), (II) the deuteration of H₂CO to CHDOH and its subsequent hydrogenation (see reaction (7)) and (III) the deuteration of H₂CO to CH₂DO and its subsequent hydrogenation (see reaction (8)). The CH₂DOH formation from CH₃O or CH₂OD is excluded in the selected formation scheme. $${\text{H}}_2\text{CO}\stackrel{+\text{H}}{\to }{\text{CH}}_2\text{OH}\stackrel{+\text{D}}{\to }{\text{CH}}_2\text{DOH}$$(6) $${\text{H}}_2\text{CO}\stackrel{+\text{D}}{\to }\text{CHDOH}\stackrel{+\text{H}}{\to }{\text{CH}}_2\text{DOH}$$(7) $${\text{H}}_2\text{CO}\stackrel{+\text{D}}{\to }{\text{CH}}_2\text{DO}\stackrel{+\text{H}}{\to }{\text{CH}}_2\text{DOH}$$(8)

When applying the SC model for the calculation of reaction probabilities for reactions with an activation-energy barrier, the hydrogenation reaction of H₂CO will much more often result in the formation of the isomer CH₃O as opposed to CH₂OH due to its much higher activation-energy barrier (E_A = 5.16 × 10³ K versus E_A = 2.00 × 10³ K). Consequently CH₃O is around 10⁴x more abundant than CH₂OH at the dust peak for the best-fit time. Without the inclusion of abstraction reactions (see also Section 4.2), a deuteration to CH₂DOH is not possible any more once CH₃O is formed. Since abstraction reactions are not considered in the first part of this work, the formation of CH₂DOH mainly proceeds via reaction (8) as it has the lowest activation-energy barrier: E_A = 1.28 × 10³ K compared to E_A = 5.16 × 10³ K for the other two options. When applying the RDC model, however, the formation of CH₂OH is not hindered significantly by the reaction’s high activation-energy barrier, as both reactants can undergo numerous attempts to react until one of them diffuses away again. In the extremely cold environment of pre-stellar cores the diffusion rates, even of hydrogen atoms, are much lower than the typical collision frequencies ν_AB. Consequently, the efficiency factor ƒ_act is much larger than in the SC model, resulting in an almost equal abundance of CH₂OH and CH₃O, while the abundance of CH₂DO is decreased. Additionally, the D/H ratio on the grain surface in the RDC model is lower – around two order of magnitudes at the dust peak – than in the SC model. Combined with the lower D/H ratio and the slower hopping rate of D atoms, the higher amount of CH₂OH decreases the relative CH₂DOH formation as a hydrogenation reaction of CH₂OH to form CH₃OH is more likely than the respective deuteration reaction to form CH₂DOH.

The impact of applying the various non-diffusive reaction mechanisms is again quite minor, which is in agreement with the results presented in Jin & Garrod (2020) (see Figure 11). In the models described above, methanol forms mainly by successive H/D-addition reactions to CO. Assuming a more efficient diffusion process is not necessary, when applying the RDC model, as the activation-energy barriers can be overcome more often due to the possibility of multiple reaction attempts, thereby reducing the need for a high amount of meeting events between the reactants. Reactions between two heavier reactants, such as CH₃ + OH → CH₃OH or similar reactions for precursors of methanol, do not contribute significantly.

Fig. 2 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various combinations of the introduced non-diffusive mechanisms (model ND1-ND6; upper panel) and of diffusive models, including the diffusion of hydrogen (and deuterium) atoms by quantum tunneling (model D2; red) or a diffusion-to-binding energy ratio of Ed/Eb of 0.2 (fast thermal diffusion; model D3; orange; lower panel). A diffusive model with a diffusion-to-binding energy ratio of Ed/Eb of 0.55 (slow thermal diffusion; model D1; grey) is shown as a reference. All models are performed within the framework of the single collision model. The observed profiles, ranging from the dust peak into the direction of the methanol peak, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas). The results are presented for the best-fit time of t = 3.0 × 105 yr.

Fig. 3 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various combinations of the introduced non-diffusive mechanisms (model ND7-ND12). A diffusive model with a diffusion-to-binding energy ratio of Ed/Eb of 0.55 (slow thermal diffusion; model D4; grey) is shown as a reference. All models are performed within the framework of the RDC model. The observed profiles are depicted in black (errors as grey-shaded areas). The results are presented for the best-fit time of t = 3 × 105 yr.

3.3 Introducing H₂CO + CH₃O → CH₃OH + HCO

Typically, the methanol formation pathway solely comprises successive H-addition reactions (or H-abstraction reactions) of atomic H to CO (e.g. Watanabe & Kouchi 2002; Fuchs et al. 2009). The hydrogen atom and its deuterated counterpart are both light enough to efficiently diffuse over the surface of the dust grain even at temperatures ranging around 10 K. If at least one of the reactants can diffuse efficiently, the contribution of non-diffusive reaction mechanisms to the formation yields of a reaction is generally small.

However, as Jin & Garrod (2020) showed, non-diffusive reaction mechanisms can have major contributions for reactions where both reactants are too heavy to efficiently diffuse over the surface. Several such reactions for the formation of methanol were proposed and tested in recent years. In their microscopic kinetic Monte Carlo (KMC) code, Simons et al. (2020) included several H-abstraction reactions between heavy intermediate products of the methanol formation pathway – HCO, H₂CO and CH₃O – with each other. While most of those reactions seem to produce only small methanol yields, Simons et al. (2020) surprisingly derived a contribution of 60–90%, depending on the physical conditions, for the reaction between H₂CO and CH₃O: $${\text{H}}_2\text{CO}+{\text{CH}}_3\text{O}\to {\text{CH}}_3\text{OH}+\text{HCO},$$(9)

making it the dominant pathway for methanol formation over the simple H-addition reaction. This result was later confirmed experimentally by Santos et al. (2022). These authors performed several H₂CO and H (or other combinations with D₂CO and D) co-deposition experiments with H/H₂CO ratios of either 10 or 30 in a temperature range between 10 K and 16 K. Similar to Simons et al. (2020), they find reaction (9) to be the dominant pathway with a contribution of at least 83% independent of the H/H₂CO ratios and temperatures. Additionally, they conclude that forming methanol predominantly by reaction (9) would result in a decrease of the expected deuterium fractionation due to the kinetic isotope effect.

Although these findings would signify a major paradigm shift in the methanol formation scheme, most rate-equation codes, unlike microscopic Monte Carlo codes, are not equipped to properly test the efficiency of chemical reactions between two heavy reactants in low temperature regimes, as they do not account for the possibility that the reactants could be already in close proximity by chance. Finally, with the inclusion of non-diffusive reaction mechanisms, we are in a position to test the efficiency of reaction (9).

We adopted the branching ratio of 0.4 from Simons et al. (2020) for reaction (9). They also include the reaction channel H₂CO+CH₃O -> CH₃OCH₂O and a nonreactive channel to account for cases where the relative geometry of the reactants is not suitable for a reaction and the reactants are unable to reorient. In the present work, these other reaction channels are only implicitly accounted for by using a branching ratio unequal one. As reaction (9) is a radical-molecule reaction, it has an activation-energy barrier, which was estimated to be 2670 K by Álvarez-Barcia et al. (2018). Subsequently, reaction (9) is cloned to generate the reactions involving the various deuterated counterparts of the reactant and product species, resulting in 16 distinct versions. After adding the resulting reaction list to the chemical network used above, we run the diffusive model D4, resulting in model NM1, and the non-diffusive model including three-body reactions ND9, resulting in model NM2. The results are depicted in Figure 4. The former run is performed to verify that reaction (9) is inefficient when assuming purely diffusive reactions to occur, which can be indeed confirmed, as the column density profiles of both CH₃OH and CH₂DOH of model D4 and NM1 are almost identical. The latter run is performed, because non-diffusive reaction mechanisms had the largest contribution for models including three-body reactions in the previous models presented above. In contrast to the results of Simons et al. (2020) and Santos et al. (2022), we only find negligible contributions to the methanol formation yields from reaction (9), the majority of methanol is still formed by the addition of atomic hydrogen to CH₃O or to CH₂OH. Our explanation for these contradicting results is that the occurrence rate of either H₂CO or CH₃O being produced in close proximity to the other reactant is much smaller than a hydrogen atom diffusing to a binding site, where a CH₃O molecule is sitting.

Apparently, this is not the case for the other two works by Simons et al. (2020) and Santos et al. (2022). Both groups of authors argue that reaction (9) might be the dominant route in their works due to a higher availability of H₂CO and CH₃O close to each other as opposed to hydrogen atoms on the grain’s surface. The chemical network in Simons et al. (2020) is limited to a hydrogen addition chain from CO to form methanol that was extended to also account for some recent experimental results on the formation of larger COMs (see Fedoseev et al. (2015) and Chuang et al. (2016)). Considering significantly more chemical reactions would be prohibitively computationally costly for a microscopic Monte Carlo code to perform. Naturally, the surface of their simulated dust grains only consists of the limited amount of species included in their chemical network, which could cause an overestimation of the occurrence rate of H₂CO and CH₃O molecules in close proximity to each other. Chemical codes that are based on the rate-equation approach are not that severely limited in the same sense. The simulations presented in this work include over 1000 species that are coupled by ≈75 000 reactions, resulting in a more realistic surface composition for interstellar dust grains. Similarly is a H/H₂CO ratio of 30 or 10 used in the H₂CO and H codeposition experiments by Santos et al. (2022) quite low as compared to the interstellar conditions, which could again result in an overestimation of the availability of H₂CO close to CH₃O molecules. In a more recent work by Jiménez-Serra et al. (2025), an adaption of the KMC code used in Simons et al. (2020) was employed to model the ice chemistry in the Chamaeleon I cloud. The chemical network used in this work also includes relevant species and reactions for the formation of CO₂, CH₄, NH₃, and H₂O, in addition to the network used in Simons et al. (2020). Interestingly, with this larger network Jiménez-Serra et al. (2025) derive only a contribution of ≈ 10% for reaction (9) and ≈90% for the hydrogenation of CH₂OH and CH₃O. However, since the physical conditions in these two works are not identical, the comparison of both results has to be regarded with some caution.

On the other hand, the specific properties of hydrogen diffusion, its dependency on ice structure, binding sites and - energies and how to include those into chemical codes is still under debate. An overestimation of the hydrogen diffusion rate on the surface could also result in an overestimation of the hydrogenation reaction’s efficiency and consequently lower the contribution of reaction (9).

Fig. 4 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of the diffusive model (D4) and a model including three-body reactions in addition to diffusive reactions (ND9). Both models use the methanol formation scheme presented in Riedel et al. (2023). Subsequently, they were redone with a modified chemical network, including the reaction H2CO + CH3O → CH3OH + HCO (NM1 and NM2). The observed profiles are depicted in black (errors as grey-shaded areas).

4 Discussion

In addition to the introduction of non-diffusive reaction mechanisms, we explored various parameter modifications. They proved to be either necessary for the inclusion of non-diffusive chemistry or showed large effects on the CH₃OH column density or N(CH₂DOH)/N(CH₃OH) ratio.

4.1 Removing surface H₂

The buildup of an unreasonable amount of H₂ on the surface of dust grains in cold dense cores is a common issue for astrochemical codes. The H₂ adsorption rate is large due to its high abundance in the gas phase, while desorption rates in the extremely cold and highly shielded environments are quite low. The H₂ problem becomes more severe in three-phase models, as H₂ molecules are incorporated into the mantle and locked in, resulting in a very high amount of mantle layers. Additionally, we found in this work that when applying the Eley–Rideal mechanism, the code will not converge to a solution, unless there is also a mechanism in place that removes the excess H₂ from the surface of the dust grains.

The H₂ excess is likely caused by the usage of binding energies appropriate for a pure water ice surface, while a realistic grain surface is composed of mostly H₂O, CO, CO₂, and ill constrained amounts of H₂. Cuppen et al. (2009) proposed that the binding energy to an H₂ surface could be around 10 times weaker than to other surface types. This motivated Garrod & Pauly (2011) to introduce effective binding E_b,eff and diffusion energies E_d,eff considering the fractional H₂ surface coverage θ(H₂): $$E_{\text{b,eff}}=E_{\text{b}}\left[1-\theta \left({\text{H}}_2\right)\right]+0.1E_{\text{b}}\theta \left({\text{H}}_2\right).$$(10)

In Equation (10), the surface fraction θ(H₂) that is covered by H₂ contributes only 0.1 times as much to the effective desorption energy as the remaining surface fraction. The main effect of the method is a significant reduction of the amount of H₂ on the surface and in the mantle of the dust grain. However, in Garrod & Pauly (2011), the method is applied to every grain species and consequently also causes a decrease of the binding energies of other species. The same is true for the diffusion energies E_d, as they are usually coupled to the binding energies by the diffusion-to-binding energy ratio E_d/E_b, which is typically due to the lack of sufficient experimental constraints, simply a single, constant value. Therefore, applying the above-described method causes also a decrease in the diffusion energies, increasing the grain species’ mobility.

Slightly more complex variations of this method have also been adopted in several other works. Taquet et al. (2014) considers not only H₂ as a surface constituent with distinctly different binding energies from H₂O, but also bare grain and pure ice; Furuya & Persson (2018) also considered CO, CO₂ and CH₃OH. Also, Garrod et al. (2022) propose a refinement of their original formulation, where they assume that every surface species is bound to four adjacent binding partners that collectively contribute to the binding energy. Each of the binding partners can either be an H₂ molecule or not, introducing five types of binding sites of varying strength (1.0 for no H₂ molecules; 0.1 for 4 H₂ molecules). Considering this more intricate version mainly results in a further reduction of the H₂ surface coverage and a lessened influence of the H₂ content on the diffusion and binding energies of other species.

We tested different methods to remove the excess surface H₂: scaling the binding energies of H₂ and its isotopologues by a factor of 0.1 (model D4) and several versions of the original formulation by Garrod & Pauly (2011) applied to all species (model BE2); only to H₂ and its deuterated isotopologues (model BE3); or only to H₂, H, and their isotopologues (model BE4). Model BE3 and BE4 were tested, as heavier species should be able to penetrate through the single monolayer of H₂ and bind as usual to the underlying ice, which is backed by laboratory experiments that do not find evidence for binding energy changes in experiments with H₂ presence in the apparatus (priv. comm. with Gleb Fedoseev). All models presented in Section 4.1 apply only diffusive reaction mechanisms. Non-diffusive reaction mechanisms are not included, as we wanted to be able to compare the effects of the various H₂ removal methods to a reference model (model BE1) without one. Figure 5 presents the column density profiles of CH₃OH (left) and N(CH₂DOH)/N(CH₃OH) profiles (right), while Figure 6 shows the surface coverage of the grain as a function of time (x-axis) and radius (y-axis). In general, the CH₃OH column densities and N(CH₂DOH)/N(CH₃OH) of the modified models (D4 and BE2-BE4) differ only by a factor of around two or less when compared to the reference model (BE1). The largest difference to the reference model is seen when using scaled H₂-binding energies (model D4), which is likely due to the fact that it represents the most extreme H₂ removal method. Its application removes the surface H₂ almost completely, decreases the CH₃OH column density by approximately a factor of two and increases the N(CH₂DOH)/N(CH₃OH) by the same amount. It is noteworthy, that we partially see the opposite trend when this method is applied in the context of the single collision model (see Appendix A). The other H₂ removal methods (models BE2-BE4) produce CH₃OH column densities and N(CH₂DOH)/N(CH₃OH) ratios that are closer to the reference model and also each other. For model BE2, where the Garrod & Pauly (2011) method is applied to all grain species, are the CH₃OH column densities slightly higher than in the reference model, while for model BE3 and BE4 the difference to the reference model is negligible. Models BE2-BE4 show a similar level of H₂ removal, with the H₂ coverage being altered from 99.9% to 56% for the best-fit time of t = 3.0 × 10⁵ yr at the position of the dust peak. Interestingly, the N(CH₂DOH)/N(CH₃OH) ratio of model BE1 and model BE2-BE4 is lower than that of model D4 with only slight differences between the former models. Overall, Model BE3, where the Garrod & Pauly (2011) method is applied to only H₂ and its isotopologues, shows the smallest difference to the reference model. It is, however, difficult to decide which method is the closest to a realistic description of the grain surface, as we currently still lack good constraints on the amount of surface H₂ in pre-stellar cores. A removal of the excess surface H₂ is mainly necessary to be able to run the Eley-Rideal mechanism, as otherwise the code will not converge to a solution. In principle, all tested mechanisms are able to remove this convergence problem, but the method of scaling the H₂-binding energies (model D4) shows the lowest computative cost and the least issues with convergence. Therefore, we adopted it for the majority of this work.

Fig. 5 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various diffusive models, that we used to test different approaches to prohibit an unphysical buildup of H2 on the surface of dust grains. A reference model (BE1) without a H2 removal method is shown. Additionally, we present a model (D4), where the binding energy of H2 and its isotopologues were scaled by 0.1 and several models adopting the H2 removal method proposed by Garrod & Pauly (2011) applied to all species (BE2), only to H2 (BE3) and only to H and H2 (BE4). The observed profiles are depicted in black (errors as grey-shaded areas).

Fig. 6 Coverage of dust grain as a function of time (x-axis) and radius (y-axis). The left column shows the percentage of empty binding sites, which are typically filled quickly by absorbed and newly formed molecules. The right column shows the percentage of H2.

4.2 Combining H-abstraction reactions and chemical desorption

Various laboratory experiments (Hidaka et al. 2009, Chuang et al. 2016; Minissale et al. 2016c) agree that H-abstraction reactions play an important role in the methanol formation scheme. Theoretically, the inclusion of H-abstraction should have two important effects on the gas phase reservoirs of CH₃OH and CH₂DOH. (I) The reactive desorption rate of methanol increases due to cyclic H-abstraction events (e.g. Jin & Garrod 2020): when a CH₃OH molecule (or a CH₂DOH molecule) is formed, it can be desorbed into the gas phase with a certain probability. However, in the majority of cases, it stays on the grain surface, where it can be subject of H-abstraction reactions that reduce it back to CH₃O or CH₂OH, which can be hydrogenated to methanol again and thereby desorbing with some probability into the gas phase. In principle, the described addition and abstraction steps can be repeated multiple times, thereby amplifying the reactive desorption process beyond what is suggested by the reactive desorption efficiency. (II) The N(CH₂DOH)/N(CH₃OH) is increased due to preferential abstraction of CH₂DOH to CH₂DO and CHDOH (Hidaka et al. 2009): the reaction CH₂DOH + H has three possible out-comes in our chemical network. It can either reduce to CH₂OH, from where a hydrogenation to CH₃OH is more likely to happen than a deuteration to CH₂DOH again, or it can reduce to CHDOH and CH₂DO, from where a hydrogenation to CH₂DOH is more likely to happen than a deuteration to CHD₂OH. The branching ratio (0.17) for the former reaction, likely reversing the deuteration, is lower than for the latter two possibilities (0.33 and 0.5, respectively), preserving the deuteration. It is noteworthy that the likelihood of hydrogenation versus deuteration reactions depends on the ratio of deuterium to hydrogen atoms and the assumed mode of diffusion on the surface of the dust grains. In case the D/H ratio is sufficient to compensate the slower diffusion of deuterium atoms, deuteration reactions can become more likely than hydrogenation reactions. However, this is not the case for the models presented in this Section.

Figure 7 shows the impact of the inclusion of the H-abstraction reactions following the Hidaka et al. 2009-scheme, and their combination with either a simple constant reactive desorption efficiency of 1% (model D5) or the experiment-based reactive desorption efficiencies following Minissale et al. (2016b) (model D7). For the purpose of comparison, we also show the respective models solely including H-addition reactions (model D4 and D6).

Models D4 and D5, using a constant reactive desorption efficiency of 1% show both theoretically expected effects. Model D5, which includes H-abstraction reactions, has higher column densities of both CH₃OH and CH₂DOH. For CH₃OH the increase can be attributed to cyclic H-abstraction (effect I), whereas for CH₂DOH it is a combination of both described effects (effects I and II). In alignment with our expectations, we find that the increase in CH₂DOH is higher than the increase in CH₃OH, resulting in an approximately two percentage points higher N(CH₂DOH)/N(CH₃OH) ratio.

Models D6 and D7 apply an experiment-based desorption mechanism following the proposed equation by Minissale et al. (2016b), which produces an individual desorption efficiency for every product species and surface coverage. In general, they produce lower column densities than model D4 and D5 for both CH₃OH and CH₂DOH, which is expected as the experiment-based reactive desorption mechanism yields much lower desorption efficiencies than 1% for heavier molecules. Given the used values for reaction enthalpies H_form and binding energies E_b (presented in Table A1 of Riedel et al. 2023), the desorption efficiency is 0.05 for CH₂OH→CH₃OH and 0.08% for CH₃O → CH₃ OH on a pure CO surface.

Models D6 and D7 do not follow the same trend as described above for model D4 and D5. The column densities of CH₃OH and CH₂DOH are both increased due to the cyclic H-abstraction (effect I) in model D7. However, the increase for CH₂DOH is lower than for CH₃OH for radii above 4000 AU, consequently resulting in a lower N(CH₂DOH)/N(CH₃OH) ratio for model D7 compared to model D6. The smaller increase in CH₂DOH is caused by different reactive desorption efficiencies for the different formation reactions of CH₂DOH. Reaction (6) (CH₂OH → CH₂DOH) and reaction (7) (CHDOH → CH₂DOH) both have a reactive desorption efficiency of around 0.05%, while reaction (8) (CH₂DO → CH₂DOH) has a reactive desorption efficiency of around 0.08%. For the inner radii, CH₂DOH formation is dominated by reaction (8). For the outer radii, however, the efficiency of reaction (6) increases, as the D atoms mobility is increased due to the higher dust temperatures. Therefore, the effective CH₂DOH reactive desorption efficiency is decreased, since this formation route contributes less of its products to the gas phase.

In our previous works, abstraction reactions were left out entirely or were only used for selected parameter sets. Vasyunin et al. (2017) argued that branching ratios between H-addition and H-abstraction reactions were not measured accurately enough to lower the overall uncertainties of their presented model. Therefore, the authors opted to use a simpler network solely including H-addition reactions. Riedel et al. (2023) included H-abstraction reactions following the scheme proposed by Hidaka et al. (2009) but also tested selected parameter sets with a simpler network comprising only the H-addition reactions. They found the inclusion of H-abstraction reactions to cause a substantial increase in the CH₃OH column densities, while the CH₂DOH column densities were found to be decreased, consequently lowering the N(CH₂DOH)/N(CH₃OH) to values that are in worse agreement with observations. Admittedly, this result was impacted by the fact that the presented models also used the experiment-based mechanism for reactive desorption proposed by Minissale et al. (2016b) and therefore experienced the same effect as described above for models D6 and D7. The effect could be seen even more pronounced in the models presented in Riedel et al. (2023), as the use of the SC model causes CH₃OH to be formed predominantly by hydrogenation of CH₃O, which has a higher reactive desorption efficiency as the competing reaction involving its isome CH₂OH. Hence, the effective CH₃OH reactive desorption efficiency is higher than in the RDC model.

Recently, Punanova et al. (2025) pointed out that H₂CO, an important intermediary in the methanol formation chain, is typically strongly overproduced, resulting in an inverse H₂CO:CH₃OH abundance of the chemical models compared to observations. Although formaldehyde is formed efficiently on the grain, its main formation reaction is the gas phase reaction CH₃ + O → H₂CO + H. Punanova et al. (2025) propose to remedy the overproduction of H₂CO by introducing the additional reaction channel CH₃ + O → CO + H₂ + H and constraining its branching ratio (CO+H₂+H:H₂CO+H = 8:1) observationally. Since the publication of this work is very recent, its proposed constraint has not been tested by us. However, we find that it is unlikely to have large consequences for the formation of methanol, as it concerns the gas phase formation of formaldehyde. We also would like to point out that for the models presented in the present work, the H₂CO:CH₃OH ratio strongly depends on the model parameters. Indeed, the H₂CO:CH₃OH ratio for model D5 is not inverse to the observed one.

Fig. 7 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various diffusive models, either including solely H-addition reactions (round markers) or H-addition and H-abstraction reactions (diamond markers) combined with either a constant reactive desorption efficiency of 1% (grey) or an experiment-based reactive desorption efficiency based on Minissale et al. (2016b) (navy). The observed profiles are depicted in black. Errors are shown as grey-shaded areas.

Fig. 8 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of the diffusive model D4 (grey), with an initial atomic hydrogen abundance of 10−8, and diffusive model D8 (green), with an initial atomic hydrogen abundance of 0.5. The observed profiles are depicted in black (errors as grey-shaded areas).

4.3 Varying the initial H abundances

Unfortunately, a direct measurement of the H₂ abundance in cold dense cores is not possible, as the highly symmetric H₂ molecule does not possess a permanent dipole moment and is therefore only directly accessible in the UV and near- and far-infrared regime. Therefore, a commonly used tool is the use of ¹³CO or ¹⁸CO abundances combined with a conversion relation of the respective CO isotopologue’s abundance to the one of molecular hydrogen. Simultaneously, the atomic hydrogen abundance can be obtained by accompanying HI narrow self-absorption observations, as was done by Li & Goldsmith (2003). They determined [HI]/[H₂] number density ratios between 10⁻⁴ and 8 × 10⁻³ for over 30 dark clouds in the Taurus and Perseus region. Goldsmith & Li (2005) later studied four dark clouds in more detail. There, they obtained an [HI]/[H₂] ratio of 10⁻³ specifically for L1544, resulting in an H abundance of 5 × 10⁻⁴ with respect to the total proton density n_H= n(H) + 2n(H₂). Although it is a common modelling approach to start the simulation with no (or very small) atomic hydrogen abundances, in recent years several modellers adopted initial abundances for atomic hydrogen following the values suggested by the above mentioned works (e.g. Jin & Garrod (2020): 5 × 10⁻⁴ for L1544; Wakelam et al. (2021): 5 × 10⁻⁴, 1 × 10⁻³ and 1 × 10⁻² for an averaged physical structure representative of TMC1).

In this section, we investigate the effect of a variation in the initial abundances of atomic and molecular hydrogen on the CH₃OH column density and N(CH₂DOH)/N(CH₃ OH) ratio profile (see Figure 8). Here, we adopt an even higher value of 0.5 for atomic hydrogen and 0.25 for molecular hydrogen, which is a ratio more typical of translucent clouds. The initial abundances for pH₂ and oH₂ are adjusted accordingly to preserve an ortho-to-para ratio of 10⁻³. The HD abundance is divided in half, where one half remains HD and the other is thought to be in the form of H and D atoms, respectively, thereby increasing the D abundance to the same value as the HD abundance. Consequently, the initial D/H ratio is lowered from a value of 1.0 to 1.6 × 10⁻⁵. The modified initial abundances are summarised in Table 3. The diffusive model D4, not including H-abstraction reactions, is run again with the modified initial abundances, resulting in model D8.

An increased atomic hydrogen abundance is a priori expected to increase the production of molecules that are formed by hydrogenation reactions of common ice species, such as CH₃OH. Indeed, we find that the CH₃OH column density is increased by roughly a factor of 2.5, while the CH₂DOH column density is reduced by that factor, resulting in a slightly lower N(CH₂DOH)/N(CH₃OH) ratio. The decrease of the CH₂DOH column density is caused by the lower initial D/H ratio. The N(CH₂DOH)/N(CH₃OH) in model D8 is eventually able to reach values similar to those in model D4; however, the two models only converge for evolutionary times past 10⁶ yr.

A similar result is obtained in Sipilä et al. (2024) (see Appendix B). There, the authors found a highly increased H₂CO and CH₃OH formation for larger radii, while the variation for inner radii is small. Since the volume density for the outer radii is low, their contribution to the CH₃OH column density is low as well. Therefore, the CH₃OH column density of the 50:50 H/H₂ model is only increased by a factor 1.8 compared to reference model with initially no atomic hydrogen.

5 Conclusion

Surprisingly high abundances of methanol and other COMs were found in several surveys targeting dark molecular clouds (e.g. Bacmann et al. 2012; Cernicharo et al. 2012; Jiménez-Serra et al. 2016). Additionally, multiple experiments (Fedoseev et al. 2015; Chuang et al. 2016; Fedoseev et al. 2017; Qasim et al. 2019; Ioppolo et al. 2021; Fedoseev et al. 2022; Santos et al. 2022) conducted at low temperatures (T = 10 − 16 K) have proposed grain surface COM formation routes involving (at least) one reaction step with two heavy immobile reactants. These kinds of reactions pose a severe bottleneck in a diffusive-only grain chemistry and likely render the formation route inefficient. Including non-diffusive grain chemistry evades the obstacle of immobility of reactants by taking into consideration that in a (small) number of cases, the reactants are already formed in close proximity to each other.

We updated the astrochemical code pyRate by including various non-diffusive reaction mechanisms, including Eley-Rideal reactions, photodissociation-induced reactions, and three-body reactions, thereby following closely the proposed formulation by Jin & Garrod (2020) and Garrod et al. (2022). We investigated methanol and its deuterated isotopologues, as methanol and its precursors are considered to be fundamental building blocks for both gas phase and grain surface formation routes of COMs. We presented several models for the prediction of abundance, column density, and deuterium fraction profiles of CH₃OH and CH₂DOH in the prototypical pre-stellar core L1544. Subsequently, we compared our model results to our previous work (Riedel et al. 2023), exploring various methods to increase the efficiency of hydrogen and deuterium atom diffusion and with available single-dish observations presented in Chacón-Tanarro et al. (2019).

Our main conclusions are as follows. When applying the single collision model (SC model) for the derivation of reaction probabilities of reactions with an activation-energy barrier:

1. Models with a slow diffusion process (E_d/E_b = 0.55) cannot reproduce the observed column density profiles;

2. The addition of non-diffusive reaction mechanisms provides only minor contributions to methanol formation, independent of the selected mechanisms. Although the majority of CH₃OH and CH₂DOH is still produced by successive addition of H (or D) atoms, reactions between heavier products (e.g. OH+CH₃ → CH₃OH, O+CH₃ → CH₃O and OH+CH₂ → CH₂OH) can contribute significantly in certain regions and/or at certain time steps;

3. Increasing the diffusion rate of hydrogen and deuterium atoms by either a slow thermal hopping process (E_d/E_b = 0.55) that is accompanied by diffusion by quantum tunneling of H and D atoms or a fast thermal hopping process (E_d/E_b = 0.2) results in significantly better agreement with observed profiles.

   When applying the reaction-diffusion competition model (RDC model) for the derivation of reaction probabilities:

4. In the RDC model, reactions with an activation-energy barrier are significantly more likely to happen when compared to the SC model, reducing the importance of an efficient diffusion process. Models with a slow diffusion process (E_d/E_b = 0.55) match the observed column density profile for CH₃OH within a factor of around two, which is a good agreement;

5. The N(CH₂DOH)/N(CH₃OH) ratio derived with the RDC model is lower by a factor of two to three in comparison to the one derived with the SC model and the observed ratio. This is due to a combination of an increased importance of H₂CO→CH₂OH→CH₂DOH relative to H₂CO→ CH₂DO→ CH₂DOH and a lowered D/H ratio;

6. Similar to the results derived with the SC model, the addition of non-diffusive reaction mechanisms provides only minor contributions to methanol formation, independent of the selected mechanisms;

7. The reaction H₂CO+CH₃O → CH₃OH + HCO (reaction (9)) is theoretically (Simons et al. 2020) and experimentally (Santos et al. 2022) predicted to provide a contribution of 60–90% to the formation of CH₃OH. Following the introduction of reaction (9) to our chemical network, we derived only a negligible contribution for reaction (9) due to the fact that the occurrence rate of either H₂CO or CH₃O produced in close proximity to the other reactant is quite low, resulting in a very small reaction rate. Therefore, we cannot confirm the result;

8. Unlike other COMs, methanol can be considered a unique case, as its formation does not include a bottleneck reaction between two heavier reactants, which is the case for other proposed COM formation routes. If at least one of the reactants, for example, hydrogen (or deuterium) atoms, can diffusive efficiently, the contribution of non-diffusive reactions mechanisms is generally small;

9. The inclusion of H-abstraction reactions has two important effects on the gas phase reservoir: (I) The desorption rate of CH₃OH and CH₂DOH increases due to cyclic H-abstraction events, and (II) the N(CH₂DOH)/N(CH₃OH) ratio is increased due to preferential abstraction of CH₂DOH to CH₂DO and CHDOH, instead of CH₂OH. Both effects result in higher CH₃OH column densities and higher N(CH₂DOH)/N(CH₃OH) ratios;

10. A combination of H-abstraction reactions with an experiment-based reactive desorption mechanism (Minissale et al. 2016b), as opposed to a constant 1% reactive desorption efficiency, also increases the formation of CH₃OH. The N(CH₂DOH)/N(CH₃OH) ratio, however, is lower than in the model without abstraction reactions due to lower effective reactive desorption efficiencies;

11. Increasing the initial amount of atomic H from a very low value (10⁻⁸) to 50% of the total H content results in an increase of the CH₃OH column density by roughly a factor of 2.5, while the CH₂DOH column density is reduced by that factor, lowering the N(CH₂DOH)/N(CH₃OH) ratio.

The present work represents the first step in a set of projects exploring the deuteration of COMs. In a future paper, we will expand the chemical network to include COMs larger than methanol (e.g. acetaldehyde, methyl formate, and glycoalde-hyde) and constrain their expected level of deuteration. In addition, a closer look into the chemistry of doubly or even triply deuterated forms of methanol could provide valuable additional constraints for the formation of methanol and its deuterated isotopologues.

Acknowledgements

W.R., O.S. and P.C. gratefully acknowledge the funding by the Max Planck Society. R.T.G. thanks the National Science Foundation for funding through the Astronomy & Astrophysics program (grant number 2206516). The work by A.I.V. is supported via the State Assignment contract FEUZ-2025-0003.

Appendix A Scaling the H₂ binding energies

Fig. A.1 Modelled column density profiles. Comparison of the main models presented in Riedel et al. (2023) (solid lines) with a respective version of those models with a scaling factor of 0.1 applied to the binding energy of H2 and its deuterated isotopologues (dashed lines). The results are presented for the best-fit time of t = 3.0 × 105 yr. The observed profiles, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas).

An important aim of this work is to ensure comparability of the models presented in this work to our previous work (Riedel et al. 2023). Therefore, we decided to introduce all modifications in a step-by-step manner and also present the results of all intermediate steps to enable the reader to follow each modification and its effects. The models, presented in Riedel et al. (2023) were performed by applying a single collision model (Hasegawa et al. 1992) for the derivation of reaction probabilities of reactions with an activation-energy barrier. It was found that a model with only a slow thermal hopping process (E_d/E_b = 0.55, model I1) is not able to reproduce the observed methanol column densities, as it underproduces them by several orders of magnitude. Instead, it proved necessary to introduce a highly efficient diffusion process for H (and D) atoms to ensure that the reactants ‘meet’ often enough to lead to a sufficient amount of successful reactions. We tested two different options: (I) a model, using a slow thermal hopping process (E_d/E_b = 0.55) accompanied by diffusion by quantum tunneling of H and D atoms (model I2) and (II) a model, using a fast thermal hopping process (E_d/E_b = 0.2; model I3). Since we test each of those models - with slightly different parameters – in this work again, we used them as a starting point for our modifications.

To get models that include the Eley-Rideal mechanism to work, we needed to reduce the amount of H₂ built up on the surface of dust grains. For that purpose, we adopted the scaling of the binding energies of H₂ and its deuterated isotopologues by a factor of 0.1. In Figure A.1, we present the column density profiles of CH₃OH (left) and N(CH₂DOH)/N(CH₃OH) profiles (right) of the most important models (I1, I2, I3) as presented in Riedel et al. (2023) and with scaled H₂ binding energies (I4, I5, I6). The CH₃OH column density of models I4, I5 and I6 is only varied within a factor of less than two compared to the models without the scaling of H₂ binding energies. The N(CH₂DOH)/N(CH₃OH) ratio, however, is decreased by a factor of two to four, depending on the type of diffusion applied to H and D atoms. We find that models I1, I2 and I3 present a larger D/H ratio on the grain surface for the time interval most relevant for the production of CH₃OH and CH₂DOH. Atomic hydrogen (and deuterium) on the surface of the dust grain is formed either by absorption from the gas phase or by surface reactions involving H₂ or its deuterated isotopologues (e.g. H₂ + OH → H₂O + H, H₂ +CH₂ → CH₃ + H). Naturally, the removal of large quantities of surface H₂ in models I4, I5 and I6 is reducing the importance of the second pathway. Instead, proportionally more of the available surface H and D atoms are absorbed from the gas phase with a typically significantly lower D/H ratio than the one on the grain surface.

Appendix B Switching to a simpler reactive desorption mechanism

Fig. B.1 Modelled column density profiles. Comparison of the models presented in Appendix A (solid lines) with a respective version of those models with a reactive desorption efficiency of 1% independent of reaction and surface coverage (dashed lines). The results are presented for the best-fit time of t = 3.0 × 105 yr. The observed profiles, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas).

In Riedel et al. (2023), we adopted an experiment-based reactive desorption mechanism following the semi-empirical equation suggested in Minissale et al. (2016b). It derives an individual reactive desorption efficiency for every product species in a reaction, that depends on the reaction enthalpy and underlying surface. Similar versions of this mechanism were also added in several other works (e.g. Vasyunin et al. 2017; Wakelam et al. 2021). The largest difference to those works is the treatment of multi-product reactions. The method used in Riedel et al. (2023), and partially in this work, also includes a recipe for the mass-dependent partitioning of the total reaction enthalpy in the case of a two-product reaction. However, since most other works use a constant reactive desorption efficiency for every reaction and independent of the surface coverage, we chose, for the sake of comparability to the literature, to switch to this simpler version of reactive desorption. In most models of this work, we adopted a reactive desorption efficiency of 1%.

Figure B.1 presents the column density profiles of CH₃OH (left) and N(CH₂DOH)/N(CH₃OH) profiles (right) of the most important models (I4, I5, I6) presented in Riedel et al. (2023), with scaled H₂ binding energies (see Appendix A) and the experiment-based reactive desorption mechanism, and the same models with a constant 1% reactive desorption efficiency (models D1,D2 and D3). The CH₃OH column density of models D1, D2, and D3, using a constant 1% reactive desorption efficiency, is increased by roughly an order of magnitude in comparison to models I4, I5 and I6, using the experiment-based reactive desorption mechanism. The increase is caused by the comparably low reactive desorption efficiencies of the experiment-based reactive desorption mechanism for heavier products. When applying the SC model, CH₃OH is predominantly produced by CH₃O+H→CH₃OH, which has a reactive desorption efficiency of only 0.08% on a pure CO surface, given the reaction enthalpies and binding energies adopted in this work. Additionally, the N(CH₂DOH)/N(CH₃OH) ratio is increased slightly when switching from models I4, I5, and I6 to models D1, D2, and D3, as the formation reactions of CH₂DOH partially have lower reactive desorption efficiencies than the one leading to CH₃OH (CH₂DO+H→CH₂DOH: 0.08%; CHDOH+H→CH₂DOH: 0.05%, and CH₂OH+D→CH₂DOH: 0.05% on a pure CO surface). Consequently, the higher the relative reaction rates of reactions CHDOH+H→CH₂DOH and CH₂OH+D→CH₂DOH as compared to CH₂DO+H→CH₂DOH, the lower the N(CH₂DOH)/N(CH₃OH) ratio.

References

All Tables

All Figures

	Fig. 1 Physical model first presented in Keto & Caselli (2010) and described in more detail in Sipilä et al. (2019a) providing static radial profiles of the H2 volume density, n(H2) (blue, logarithmic scale), the visual extinction, AV (red), Tgas (orange), and the dust temperature, Tdust (green).
In the text

Fig. 2 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various combinations of the introduced non-diffusive mechanisms (model ND1-ND6; upper panel) and of diffusive models, including the diffusion of hydrogen (and deuterium) atoms by quantum tunneling (model D2; red) or a diffusion-to-binding energy ratio of Ed/Eb of 0.2 (fast thermal diffusion; model D3; orange; lower panel). A diffusive model with a diffusion-to-binding energy ratio of Ed/Eb of 0.55 (slow thermal diffusion; model D1; grey) is shown as a reference. All models are performed within the framework of the single collision model. The observed profiles, ranging from the dust peak into the direction of the methanol peak, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas). The results are presented for the best-fit time of t = 3.0 × 105 yr.

In the text

Fig. 3 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various combinations of the introduced non-diffusive mechanisms (model ND7-ND12). A diffusive model with a diffusion-to-binding energy ratio of Ed/Eb of 0.55 (slow thermal diffusion; model D4; grey) is shown as a reference. All models are performed within the framework of the RDC model. The observed profiles are depicted in black (errors as grey-shaded areas). The results are presented for the best-fit time of t = 3 × 105 yr.

In the text

Fig. 4 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of the diffusive model (D4) and a model including three-body reactions in addition to diffusive reactions (ND9). Both models use the methanol formation scheme presented in Riedel et al. (2023). Subsequently, they were redone with a modified chemical network, including the reaction H2CO + CH3O → CH3OH + HCO (NM1 and NM2). The observed profiles are depicted in black (errors as grey-shaded areas).

In the text

Fig. 5 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various diffusive models, that we used to test different approaches to prohibit an unphysical buildup of H2 on the surface of dust grains. A reference model (BE1) without a H2 removal method is shown. Additionally, we present a model (D4), where the binding energy of H2 and its isotopologues were scaled by 0.1 and several models adopting the H2 removal method proposed by Garrod & Pauly (2011) applied to all species (BE2), only to H2 (BE3) and only to H and H2 (BE4). The observed profiles are depicted in black (errors as grey-shaded areas).

In the text

	Fig. 6 Coverage of dust grain as a function of time (x-axis) and radius (y-axis). The left column shows the percentage of empty binding sites, which are typically filled quickly by absorbed and newly formed molecules. The right column shows the percentage of H2.
In the text

Fig. 7 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of various diffusive models, either including solely H-addition reactions (round markers) or H-addition and H-abstraction reactions (diamond markers) combined with either a constant reactive desorption efficiency of 1% (grey) or an experiment-based reactive desorption efficiency based on Minissale et al. (2016b) (navy). The observed profiles are depicted in black. Errors are shown as grey-shaded areas.

In the text

	Fig. 8 Column density profiles of CH3OH (left) and N(CH2DOH)/N(CH3OH) profiles (right) of the diffusive model D4 (grey), with an initial atomic hydrogen abundance of 10−8, and diffusive model D8 (green), with an initial atomic hydrogen abundance of 0.5. The observed profiles are depicted in black (errors as grey-shaded areas).
In the text

Fig. A.1 Modelled column density profiles. Comparison of the main models presented in Riedel et al. (2023) (solid lines) with a respective version of those models with a scaling factor of 0.1 applied to the binding energy of H2 and its deuterated isotopologues (dashed lines). The results are presented for the best-fit time of t = 3.0 × 105 yr. The observed profiles, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas).

In the text

Fig. B.1 Modelled column density profiles. Comparison of the models presented in Appendix A (solid lines) with a respective version of those models with a reactive desorption efficiency of 1% independent of reaction and surface coverage (dashed lines). The results are presented for the best-fit time of t = 3.0 × 105 yr. The observed profiles, presented first in Chacón-Tanarro et al. (2019), are depicted in black (errors as grey-shaded areas).

In the text