© The Authors 2023

1 Introduction

Methanol (CH₃OH) is the simplest O-bearing complex organic molecule (COM) in the interstellar medium and an important precursor of saturated, more massive COMs that are formed in the gas phase. Previous work, however, showed theoretically (Garrod et al. 2006) and experimentally (Geppert et al. 2005) that the presumed gas-phase formation route is very inefficient and unable to account for the observed gas-phase abundances. Simultaneously, the formation of methanol on the surface of dust grains was proposed and experimentally investigated in an extensive manner in several independent projects (e.g. Watanabe et al. 2002; Fuchs et al. 2009). It was concluded that formaldehyde and methanol can be produced by multiple successive addition reactions of CO with diffusive hydrogen atoms. The measured significant isotope effect for hydrogen and deuterium suggests that hydrogenation and deuteration on the surface proceeds via quantum tunnelling reactions (Hidaka et al. 2007) and therefore progresses at a much higher rate than in the gas phase. Additionally, the existence of so-called abstraction reactions, which remove H from the molecule and thereby reverse the addition reaction, was postulated and found in various laboratory experiments (e.g. Hidaka et al. 2009; Minissale et al. 2016c). However, there seems to be some disagreement about the exact reaction scheme and the magnitude of the reaction rates of the different experimental approaches.

For the production of more advanced COMs in the gas phase, methanol needs to desorb from the surface of the dust grains. In hot cores and corinos, the desorption of the molecular contents of the surface phase does proceed very efficiently via thermal evaporation and photoevaporation. In the cold and dark environment of pre-stellar cores, however, these mechanisms are negligible. Therefore, only very low abundances of methanol were expected to exist in the gas phase of pre-stellar cores. Surprisingly, several surveys, for example Bacmann et al. (2012), Cernicharo et al. (2012), and Jiménez-Serra et al. (2016), conducted towards dark molecular clouds found comparatively high abundances of methanol and other COMs.

In order to explain these unexpected findings, astrochemical models need to employ other mechanisms for the evaporation of surface molecules. One promising candidate, reactive desorption, is based on the energy released in an exothermic chemical reaction, which can lead to the desorption of the reaction product(s). Nowadays, most chemical codes include a simple treatment of reactive desorption following the recipe described by Garrod et al. (2007). According to those authors, a (nearly) constant reactive desorption efficiency, typically 1%, is employed independent of the desorbing molecule. However, recent laboratory experiments have found that the reactive desorption efficiency is strongly dependent on the chemical reaction and the type of underlying surface (Minissale et al. 2016b; Chuang et al. 2018).

In this paper we developed an updated version of the description of reactive desorption presented by Vasyunin et al. (2017), itself based on the experiments by Minissale et al. (2016b). Our work is motivated by the column density maps of CH₃OH and CH₂DOH and the theoretical predictions presented by Chacón-Tanarro et al. (2019). They carried out single-dish observations with the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope (also Bizzocchi et al. 2014; Vastel et al. 2014; Spezzano et al. 2016), showing that CH₃OH peaks in an asymmetric ring around the dust peak, with the strongest emission in the northern part of the pre-stellar core. Chacón-Tanarro et al. (2019) analysed the formaldehyde and methanol column densities along a cut set by the position of the dust and the offset methanol peak.

An important aim of this work is to improve upon the theoretical column density profiles made by the chemical code pyRate (hereafter model S16; Sipilä et al. 2015a, Sipilä et al. 2019b), particularly to get more accurate predictions about the deuteration of methanol. Additionally, we compared our predictions with the results of the chemical code presented in Vasyunin et al. (2017, hereafter the V17 model). These two models were used in Chacón-Tanarro et al. (2019), who found that the V17 model produced acceptable results for non-deuterated methanol (CH₃OH). However, the column density profile for singly deuterated methanol (CH₂DOH) had to be obtained using the N(CH₂DOH)/N(CH₃OH) ratio predicted by pyRate since the MONACO code (Vasyunin et al. 2017) does not treat deuteration, whereas pyRate includes an extensive description of deuterium chemistry. To carry out a similar analysis in a consistent manner, we implemented an updated version of the reactive desorption mechanism used in Vasyunin et al. (2017) into pyRate.

Vasyunin et al. (2017) and Chacón-Tanarro et al. (2019) both find that diffusion via the quantum tunnelling of atomic and molecular hydrogen had to be employed in order to explain the observed methanol column density profile. In the present paper, we employed tunnelling diffusion for hydrogen and deuterium atoms for our fiducial model, neglecting the tunnelling diffusion of molecular hydrogen. However, we also explore the option to promote thermal diffusion by decreasing the diffusion-to-binding energy ratio, E_d/E_b, which determines the threshold for the diffusion of molecules over the surface of an interstellar dust grain via thermal hopping, from 0.55 to the lowest value reported in the experimental literature (0.2; Furuya et al. 2022). Moreover, we test several alternative sets of input parameter to see how well they are able to explain the observed deuterium fraction profile.

The paper is structured as follows: Sect. 2 describes the new reactive desorption mechanism and the chemical and physical model in detail. Section 3 describes the results of the fiducial model and compares them with the observationally obtained column density and deuterium fraction profiles from Chacón-Tanarro et al. (2019) and with the results of the V17 model Vasyunin et al. (2017). Section 4 discusses multiple modifications to the chemical and physical parameters and their effects on the results. Section 5 presents our conclusions. Appendices A and B provide additional information on the chemical parameters and reaction schemes used.

2 Model

2.1 Treatment of reactive desorption

For the extension of the reactive desorption mechanism, we adopted the physical scenario proposed and experimentally justified by Minissale et al. (2016b). The mechanism is implemented following Vasyunin et al. (2017), but it is modified for this paper, as detailed below.

Minissale et al. (2016b) developed a formula that expresses the reactive desorption efficiency, RD, as a function of reaction and surface-dependent properties: $$RD=\exp \left(-\frac{E_{b}N}{ϵ\Delta H}\right),$$(1)

where ∆H is the reaction enthalpy, E_b the binding energy, N the number of the degrees of freedom (translational, rotational and vibrational) of the reaction product and ϵ the fraction of kinetic energy retained by the reaction product. We assumed that the available reaction enthalpy is distributed equally into all the degrees of freedom, N, and that only the energy going into the vertical translational degree of freedom is used for the desorption of an atom or a molecule from the surface of a dust grain. Therefore, for one-product reactions, 1/N of the reaction enthalpy is distributed into motion perpendicular to the surface of the dust grain. The number of degrees of freedom, N, of the product can be simply derived as N = 3n_atoms, with n_atoms being the number of atoms of the reaction product. This approach has the effect that the less complex molecules, consisting of fewer atoms, can use a larger share of the reaction enthalpy for vertical motion off the surface. For the more complex molecules, the total number of degrees of freedom increases quickly with increasing number of atoms, due to the fact that larger molecules have more possibilities for vibrational excitation. However, the number of translational degrees of freedom stays the same. As a consequence, a smaller share of the reaction enthalpy is converted into vertical motion. For example, a forming CO molecule, with only 2 atoms, receives 1/6 of the reaction enthalpy for vertical motion, while a forming CH₃OH molecule, with 6 atoms, receives only 1/18.

The dependence of the reactive desorption efficiency on the type of surface is expressed by the fraction of kinetic energy, ϵ, received by the reaction product with mass m when colliding with a surface element with an effective mass M: $$ϵ={\left(\frac{M-m}{M+m}\right)}^2.$$(2)

In the experiments by Minissale et al. (2016b), the effective mass, M, is a parameter that needs to be fitted. It was found to be typically much larger than the mass of a single atom or molecule of the surface species, but is more consistent with a collective behaviour of also the neighbouring molecules that is induced by the rigidity of the surface. The more rigid the surface, the higher the effective mass of the surface element M and the easier it is for the reaction products to bounce off from the surface into the gas phase. For the effective masses we mostly kept the values adopted by Vasyunin et al. (2017), namely M = 48 amu for a H₂O surface and M = 100 amu for a CO-surface, with the exception of the effective mass for bare grain, where we took the somewhat lower value of M = 120 amu that has been suggested as a common value for carbonaceous and silicate grains following the experiments conducted by Minissale et al. (2016b). Every surface that does not consist of H₂O or bare grain is treated as if it were CO. Considering that grain surfaces in the inner, very cold regions of pre-stellar cores are covered, aside from CO itself, by species with a similar molecular weight (e.g. N₂, H₂CO and CH₃OH), this seems to be a reasonable approximation given the lack of more detailed measurements. To consider that a grain surface might be covered by a mixture of bare surface, H₂O or CO, we followed the average surface composition of the dust grains over time, starting with a bare grain surface, which is then quickly covered by water ice and later on CO. Then, we scaled the individual reactive desorption efficiencies for a reaction i by the fraction of the surface sites that are covered by the particular surface type, j: $$R{D}_{\text{tot}}\left(i,t\right)={\displaystyle \sum _{j}R{D}_{\text{j}}\left(i\right)}\cdot \frac{n_j^{*}\left(t\right)}{n_{\text{tot}}^{*}\left(t\right)},$$(3)

where RD_j is the individual reactive desorption efficiency for either bare grain H₂O or CO, $n_j^{*}$ the number of surface sites inhabited by surface type j and $n_{\text{tot}}^{*}$ the total number of surface sites. The largest difference between the procedure applied in this paper and the one presented in Vasyunin et al. (2017) is the treatment of reactions with more than one reaction product. Vasyunin et al. (2017) treated multi-product reactions identically to one-product reactions, so that every reaction product receives the entire available reaction enthalpy. This approach clearly violates energy conservation. However, this simplification was considered to be negligible for the chemical network that was used in Vasyunin et al. (2017) as it includes mostly two-product reactions involving heavier reactants that proceed inefficiently at the low grain temperatures typical of pre-stellar cores.

In contrast, the method presented in this paper includes a recipe for the mass-dependent partitioning of the total reaction enthalpy in the case of a two-product reaction. We assumed that the desorbing molecules are first isolated from the surface and then undergo an elastic collision with the surface of the grain to gain velocity in the vertical direction. Considering the conservation of linear momentum, one can derive an expression for the kinetic energy, E_kin, received by a product species i: $$E_{\text{kin}}\left(\text{i}\right)=E_{\text{kin}}\cdot \frac{m_{\text{j}}}{m_{\text{i}}+m_{\text{j}}}=\frac{N_{\text{trans}}}{N_{\text{tot}}}\cdot \Delta H\cdot \frac{m_{\text{j}}}{m_{\text{i}}+m_{\text{j}}}.$$(4)

The kinetic energy received by species i is scaled by the fraction of the mass of the other reaction partner (species j) and the combined mass of both reaction products. Additionally, we derived the share of the reaction enthalpy that goes into kinetic energy as the fraction of the number of translational degrees of freedom to the total number of degrees of freedom. We note that the number of translational degrees of freedom is the combined number for both reactants: N_trans = 3n_prod, where n_prod is the number of products. The total number of degrees of freedom is similarly defined as the sum of the individual degrees of freedom for both reaction partners: N_tot = ∑_i 3n_atoms,i, with n_atoms,i being the number of atoms of reaction product i. Here, we considered that the vertical motion that is responsible for the desorption is only one of three translational degrees of freedom. This introduces an additional factor of 1/3. Finally, we could derive an equation for the reactive desorption efficiency of a two-product reaction (5), which simplifies to Eq. (1) if the occurring reaction has only one reaction product: $$R{D}_{\text{i}}=\exp \left(-\frac{1}{ϵ}\left(\frac{E_b}{\Delta H}\right){\left(\frac{1}{3}\frac{N_{\text{trans}}}{N_{\text{tot}}}\frac{m_{\text{j}}}{m_{\text{i}}+m_{\text{j}}}\right)}^{-1}\right).$$(5)

The new approach has the effect that the lighter of the two products gets the larger share of the kinetic energy and consequently has an increased reactive desorption efficiency, while the one for the heavier product is decreased. For example, considering the reaction $${\text{H}}_2\text{CO}+\text{H}\to \text{HCO}+{\text{H}}_2,$$(6)

with HCO being approximately 14 times more massive than H₂, the HCO molecule receives only 1/15 of the kinetic energy, whereas the molecular hydrogen gets 14/15. This results in a reactive desorption efficiency for a CO surface of 54.0% for the H₂ molecule and only 1.5 × 10⁻³⁶% for the HCO molecule. The procedure presented in Vasyunin et al. (2017) would provide an efficiency of 63.1% for H₂ and 0.1% for the HCO molecule. However, there are no similar reactions efficient at low temperatures in the Vasyunin et al. (2017) model.

Fig. 1 Reaction scheme for the formation of non-deuterated CH3OH by successive hydrogenation. Hydrogen molecules can also be segregated from CH3OH or its precursors via abstraction reactions.

2.2 Chemical model

We incorporated the new description of reactive desorption into the rate-equation-based chemical code pyRate, described in more detail in Sipilä et al. (2015a, Sipilä et al 2019b). It tracks the chemical evolution both in the gas phase and on the grain surface. The basis of the chemical network is the 2014 public release of the Kinetic Database for Astrochemistry (KIDA) gas-phase network (kida.uva.2014, Wakelam et al. 2015), which was extended by deuterium chemistry for molecules with up to seven atoms. The code also tracks the various spin states of the light hydrogen-bearing species H₂, ${\text{H}}_2^{+}$ and ${\text{H}}_3^{+}$ and their deuterated isotopologues, as well as multiply protonated or deuterated species involved in the water and ammonia formation networks. All together, the network includes ≈74 000 gas-phase reactions and ≈2100 grain surface reactions. For the inclusion of abstraction reactions into the methanol formation pathway in the models presented in this work, we refer to the chemical network proposed by Hidaka et al. (2009) and depicted in Fig. 1 (representation with H) and Fig. B.1 (more extensive representation with H and D). We adopted atomic initial abundances (see Table 1) taken from Semenov et al. (2010), as they were also used for the S16 model, whose improvement is the main aim of this work. Also, we employed a three-phase model, consisting of a gas phase, a chemically active surface phase, and a chemically inert mantle phase. The dust grains were assumed to be spherically symmetric with a radius of 0.1 μm.

For exploring our new treatment of reactive desorption, the choice of binding energies E_b and formation enthalpies H_form is crucial. The values of E_b and H_form are displayed in Table A.1. The binding energies are taken from Semenov et al. (2010). Most of the formation enthalpies are adopted from Du et al. (2012). For the other remaining species, the data sources are marked in Table A.1. Unfortunately, experimental values for deuterated molecules are quite scarce. For this reason, we applied the same values as for the non-deuterated isotopologues, with the exception of the species marked with a star in Table A.1, for which we found individual values in the NIST Chemistry WebBook¹.

Fig. 2 Physical model developed in Keto & Caselli (2010) that yields static radial profiles of the H2 density, n(H2) (blue, in logarithmic scale), the visual extinction, AV (red), the gas temperature, Tgas (orange), and the dust temperature, Tdust (green).

2.3 Physical model

As mentioned in Sect. 1, one of the main aims of this work is to improve the theoretical predictions made by the chemical code pyRate for the column density profiles of CH₃OH and CH₂DOH and to compare these again to the observational (Chacón-Tanarro et al. 2019 and theoretical profiles from the V17 model (Vasyunin et al. 2017). For that purpose, we extensively explored the chemical evolution in the pre-stellar core L1544 with a one-dimensional physical model. It was derived from the one presented in Keto & Caselli (2010) and described in more detail in Sipilä et al. (2019a). The model provides radius-dependent, but time-independent, values for the H₂-density n(H₂), the gas temperature T_gas, the dust temperature T_dust and the visual extinction A_V as shown in Fig. 2. The core model consists of 35 concentric shells spanning the core radius of 0.32 pc. The chemistry is solved separately for each shell, yielding a spherically symmetric spatio-temporal evolution of molecular abundances. We calculated column densities by integrating along the line of sight for different impact factors from the core centre. Afterwards, the column density distribution was convolved with a 30” Gaussian beam, corresponding to the angular resolution of the observations by Chacón-Tanarro et al. (2019).

Taking the core model as a basis, we ran several simulations varying multiple chemical and physical parameters. All presented models include the new treatment for reactive desorption with individual efficiencies for every exothermic surface reaction and various surface types, as described in Sect. 2.1. An overview of the various models is presented in Table 2.

3 Results

3.1 Fiducial model

We selected the 1D-4 model as our fiducial model because it is the closest to V17 in terms of parameter space. In the fiducial model, the diffusion of H and D atoms by quantum tunnelling is enabled, while abstraction reactions as shown in Fig. 1 (or their deuterated analogues) are not included. This choice is discussed in more detail in Sect. 4.1.1.

We calculated column density profiles for all species observed by Chacón-Tanarro et al. (2019), for several time steps in the range of 10⁵ yr to 10⁶ yr in the fiducial model. Figure 3 shows the molecular abundances with respect to H₂ and Fig. 4 shows the column density profiles for H₂CO, CH₃OH and CH₂DOH for four different time steps (1 × 10⁵ yr, 3 × 10⁵ yr, 5 × 10⁵ yr, and 1 × 10⁶ yr). We note that these species freeze out onto dust grains in the very cold centre of the pre-stellar core and then peak at a density of n(H₂) ≈ 10⁴ cm⁻³, which corresponds to a radius of ≈5200 AU in the theoretical profiles, offset from the position of the dust peak. For a time of t = 3.5 × 10⁵ yr, the CH₃OH abundance reaches its maximum value of approximately 4 × 10⁻¹⁰, which is in good agreement with values observed in various pre-stellar cores (e.g. Scibelli & Shirley 2020; Harju et al. 2020; Spezzano et al. 2020; Punanova et al. 2022).

In Chacón-Tanarro et al. (2019), the time step in the S16 model was chosen such that the simulated CO column density approximately matched the observed one (Caselli et al. 1999), which occurred at a very early time of 3 × 10⁴ yr. However, S16 did not include CO self-shielding, and is therefore probably underestimating the total amount of CO. For this very early time step, the S16 model produces a centrally flat CH₃OH column density profile with an amplitude of N ≈ 1 × 10¹¹ cm⁻², underproducing the observed values for methanol of 3.9 × 10¹³ cm⁻² at the dust peak and 5.9 × 10¹³ cm⁻² at the methanol peak by roughly two orders of magnitude. On the other hand, the N(CH₂DOH)/N(CH₃OH) ratio in the S16 model approaches unity, thereby overestimating the observed deuterium fraction of 0.07 by a factor of 10. We note that the chemical network that was used for the S16 model includes reactions like CH₂DOH + H₃O⁺ → CH₃OHD⁺ + H₂O, or analogues for other deuterated forms of methanol, with the effect of increasing the deuterium fraction. We do not permit this sort of exchange to happen in the revised chemical network, as this would require the addition of a (hydrogen) atom as well as the exchange of a hydrogen and a deuterium atom between the different functional groups of methanol, which we consider unlikely.

To estimate a new best-fit time for the fiducial model, we calculated the χ² values of the observed central column density versus the corresponding values in the fiducial model for different time steps and species. Subsequently the time step where the sum of the χ² values of all species is minimised was taken as the best-fit time. For the 1D-4 model, this corresponds to a best-fit time of t = 3 × 10⁵ yr, which is roughly consistent with estimations from other molecules (e.g. Redaelli et al. 2019a, Redaelli et al. 2021b). Additionally, it coincides with the occurrence of the highest methanol column density in the probed time frame. For t = 3.0 × 10⁵ yr, the fiducial model reaches a column density of methanol of ≈1.5 × 10¹²cm⁻², which is roughly an order of magnitude lower than the column densities determined by the observations. Given the large uncertainties that chemical modelling typically experiences, this is an acceptable agreement (Vasyunin et al. 2004, Vasyunin et al. 2008, Wakelam et al. 2010, and references therein).

Figure 5 compares the modelled deuterium fraction N(CH₂DOH)/N(CH₃OH) with the observed one. At the early time step (t = 1 × 10⁵ yr), the deuterium fraction profile is slightly higher than the observed one, but the shape of the profiles are nearly identical. The modelled deuterium fraction flattens with time. This behaviour can probably be explained by the use of the static physical model that is employed here. At the best-fit time of t = 3 × 10⁵ yr, the deuterium fraction assumes an almost completely flat profile with a value of ≈0.08. At the intermediate (t = 5 × 10⁵ yr) and late (t = 1 × 10⁶ yr) time steps, the N(CH₂DOH)/N(CH₃OH) ratio starts to decrease slightly, but stays at a level that is in agreement with the observed one.

Fig. 3 Gas-phase abundance profiles of H2CO, CH3OH, and CH2DOH for the 1D-4 model for four different time steps ranging from 105 yr (left) to 106 yr (right). The best-fit time is t = 3 × 105 yr.

3.2 Comparison with V17

The fiducial model (1D-4) shows a much better agreement for the column densities of non-deuterated methanol than the S16 model, which was produced by an earlier version of pyRate, when compared to the V17 model produced by the chemical code MONACO and presented in Vasyunin et al. (2017) and Chacón-Tanarro et al. (2019, see our Fig. 6). We set the values of the various physical parameters to correspond to those in V17 as closely as possible. Still, several differences remain between MONACO and pyRate, and these cannot be fully understood without a detailed direct comparison of the two models, which is out of the scope of this paper.

Here, we point out some of the most noticeable differences. The chemical models are not identical, as we used different chemical networks. MONACO is specialised to describe the formation of COMs, while pyRate concentrates on the description of deuteration. V17 does not consider deuteration at all. Also, pyRate’s chemical network contains ‘backward’ abstraction reactions adopted from Hidaka et al. (2009). We removed them for several models presented in this paper, including the fiducial model 1D-4, in order to compare our model to the V17 model. In V17, these reactions are deliberately left out due to their badly constrained reaction rates.

Both codes use a so-called three-phase grain model with multiple layers as opposed to a simpler two-phase model. Three-phase models usually consist of three distinct phases: the gas, the ice surface and the bulk phase, which again can be subdivided into individual layers. Two-phase models only distinguish between gas and surface phase. In most astrochemical codes, the bulk is simply a chemically inert storage of accreted molecules, while the chemical reactions occur solely in the gas and surface phase. A two-phase model has no such storage. Every molecule on the grain can react or desorb at any time. The introduction of layers into a grain model enables the storage of molecules in the bulk in the order that they are accreted. Gas-phase molecules can only accrete to the surface phase, where they can react with another surface molecule. The surface phase usually consists of only one layer. If all binding sites in the surface phase are filled, molecules are transferred continuously into the bulk phase, which keeps growing.

The MONACO code has a more advanced description of the physical processes taking place on dust grains than most other astrochemical codes. In MONACO, the bulk experiences a slow type of diffusion and species are therefore able to meet one another and react, instead of just being stored away. Moreover, the surface phase consists not only of the uppermost layer, but of the first four layers in order to also allow atoms to be diffused in the vertical direction. The fiducial model presented in this paper also uses a multilayer dust model, but with a chemically inactive bulk phase and only one layer in the surface phase, which reduces the surface area on which methanol and its precursors can be hydrogenated.

To evaluate the consequences of our choice to use only one layer in the chemically active surface phase in our models, we also tested a modification of our fiducial model, 1D-4-1, where we set the number of layers in the surface phase to four and, as a reference, a two-phase model of the fiducial model setup, 1D-4-2, where the entirety of the ice is available for reactions. For 1D-4-1, we find that both the CH₃OH and the CH₂DOH column densities are increased by a factor of a few in comparison to the fiducial model in the entire considered time frame. This makes sense as with the increase in layers in the surface phase the model approaches the two-phase model, showing consistently higher column densities for both isotopologues of methanol. This behaviour is credited both to an increased overall production of methanol and an increased reactive desorption efficiency in the very centre of the pre-stellar core, due to a higher coverage of the surface with CO. However, the deuterium fraction of the two-phase model is quite low (see also Sect. 4.1.2), as the increase in column densities for CH₂DOH is lower than for CH₃OH. Surprisingly, this is not the case for the model with four layers in the surface phase, as is depicted in Fig. 5. Instead, we find a considerable increase in the deuterium fraction compared to the fiducial model for early time steps, reaching a value of up to ≈ 0.16 at t = 1 × 10⁵ yr. The deuterium fraction of the model with four layers in the surface phase drops slightly below that of the fiducial model for time steps around when the methanol column density peaks, but increases again above it for later time steps.

Both models, the fiducial and V17, employ a similar form of diffusion. The diffusion-to-binding energy E_d/E_b is set to a constant value of 0.55, as suggested by Minissale et al. (2016a), for all surface species. Additionally, both codes consider the inclusion of diffusion by quantum tunnelling. Although, while pyRate only follows the tunnelling diffusion of H (and D), MONACO also traces that of H₂.

The mechanism for reactive desorption in the two codes is similar: it is altered here for chemical reactions with multiple products to ensure that it adheres to the conservation of energy, and it is identical for reactions with a single product. If one does not introduce abstraction reactions in the formation scheme of methanol, reactions with two reaction products are considered to be negligible, as the formation just proceeds by successive addition reactions of H. We note that even though the procedure for single product reactions is the same in both codes, they can produce unequal reactive desorption efficiencies owing to the fact that we used different sets of enthalpies and binding energies. The binding energies and formation enthalpies for the models presented in this work and in the S16 model are adopted from Semenov et al. (2010), Du et al. (2012), the NIST Chemistry WebBook² and KIDA³. Detailed information can be found in Table A.1. As little modification as possible has been made to the list of these values, as the main aim of this work is to improve the result of the S16 model. Only if necessary for the calculation of the reactive desorption efficiencies, values have been added. This concerns all values adopted from the NIST Chemistry WebBook and KIDA. For the V17 model, the enthalpies are also mostly taken from Du et al. (2012) and the binding energies from Minissale et al. (2016b). Apparently, discrepancies between the two codes also exist for some of the reactions involved in the hydrogenation chain towards methanol: $$\text{HCO}+\text{H}\underset{\text{py\hspace{0.17em}Rate:15.4\%}}{\overset{\text{MONACO}:5.4\%}{\to }}{\text{H}}_2\text{CO}$$(7) $${\text{CH}}_2{\text{OH/CH}}_3\text{O}+\text{H}\underset{\text{py\hspace{0.17em}Rate}:0.05\%/0.08\%}{\overset{\text{MONACO}:0.64\%}{\to }}{\text{CH}}_3\text{OH}.$$(8)

For example for reaction (7) and (8) the two codes yield very different reactive desorption efficiencies. In pyRate, H₂CO is desorbed more efficiently before it has the opportunity to react further and eventually form CH₃OH, which is also less likely to desorb from the surface of the dust grain compared to the formation scenario in MONACO.

Figure 6 presents the column density profiles of CO, H₂CO, CH₃OH and CH₂DOH of models 1D-4, S16 and V17. In the case of V17, the CH₂DOH column density was derived by scaling the CH₃OH column density with the respective deuteration ratio from S16, as MONACO does not include a description of deuterium chemistry. We note that the column density profiles are shown at different time steps. The V17 model and the S16 model presented in Chacón-Tanarro et al. (2019) showed the time steps when the CO column density in the respective model is comparable to the observed value. For V17 this corresponds to t = 1.6 × 10⁵ yr, which coincides with the peak of the COM abundances. The same estimation for S16 yields a very early time step of t = 3.0 × 10⁴ yr. The resulting column densities are shown in Fig. 6 as the dark red line. We have derived a new best-fit time corresponding to the lowest χ² value of the observed column densities versus the corresponding values in the 1D-4 model: t = 3.0 × 10⁵ yr. This time step corresponds roughly with the peak of the methanol abundance in the 1D-4 model and is also in good agreement with the one estimated by the V17 model, varying only by a factor of 2. In Fig. 6 we show the results for the S16 model and the 1D-4 model at the new best-fit time step as bright red or green line, respectively. The S16 model produces less non-deuterated and singly deuterated methanol at this later time step, which is caused by gas-phase chemical reactions where atom exchanges between the two functional groups of methanol were allowed. Such reactions are not allowed in the present work. The 1D-4 model, on the other hand, is more consistent with the results of the V17 model, though it still underestimates the column densities of methanol by roughly an order of magnitude. The V17 model overestimates the CH₃OH column density and as consequence also CH₂DOH column density, which was credited to a likely overestimation of the reactive desorption efficiency. All models overestimate the amount of gas-phase H₂CO. The V17 and S16 models, which are evaluated at a time step to match the observed CO column densities, are off by more than an order of magnitude, whereas S16 at t = 3.0 × 10⁵ yr and 1D-4 produce only twice the observed column density. Although the V17 model and the newly developed fiducial model, 1D-4, still differ for the CH₃OH column density by almost two orders of magnitude (one order of magnitude of overestimation by V17 and one order of magnitude underestimation by 1D-4), we are able to confirm qualitatively some results of Vasyunin et al. (2017). The most important conclusion is that we are only able to reconcile our models with the observed values if we consider some form of enhanced diffusion on the surface of dust grains. One possibility to get a higher diffusion rate is to enable diffusion by quantum tunnelling of H and D atoms in the models. Other options are discussed below. Also, similar to Vasyunin et al. (2017), we conclude that other forms of non-thermal desorption, for example cosmic-ray-induced desorption (CRD) and photo-desorption, seem to have a negligible impact on the release of methanol and its precursors into the gas phase.

Fig. 4 Column density profiles of H2CO, CH3OH, and CH2DOH for the 1D-4 model for four different time steps ranging from 105 yr to 106 yr. The lines with markers show the modelled results, integrated along the line of sight and convolved with a 30″ beam. The solid lines show the observationally obtained column density profiles obtained by Chacón-Tanarro et al. (2019) by taking a cut through the dust and methanol peaks. The grey-shaded areas indicate the error bars of the column densities. The position of the dust peak is at r = 0 AU, while r > 0 AU is the direction towards the methanol peak.

Fig. 5 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for the 1D-4 model for four different time steps ranging from 105 yr (top left) to 106 yr (bottom right). Additionally, we show two variations of this model: one with four layers in the chemically active surface phase instead of one (1D-4-1), and one with a two-phase model (1D-4-2; see also Table 2 and Sect. 3.2 for a more detailed explanation). The coloured lines show the column density ratio of the models, while the black line indicates the observed ratio (errors as grey-shaded areas).

4 Discussion

In addition to the comparison with the V17 model, we explored various alterations of our model, mainly to investigate the effect on the deuteration of methanol and possibly improve the agreement between the predicted and observed N (CH₂DOH)/N (CH₃OH) ratio. An overview of the various models and their modified chemical and physical parameters is presented in Table 2. The 1D-4 model was picked as the fiducial model because of its closeness to V17 in terms of parameter space. However, for most of the variations of the other physical parameters, we used the 1D-3 model as the point of reference, comprising the abstraction reactions, as most other modelling works include them. In Table 2, models that are derived by modifying the fiducial model 1D-4 are indicated with a star (^⋆) and those derived from 1D-3 with a dagger (^†). We only varied one physical parameter at a time, in order to undoubtedly ascribe the altered results to the made modifications. For the discussion of the effects of the variations, we roughly divided them into chemical and physical variations. We collected an overview of the resulting N (CH₂DOH)/N (CH₃OH) ratios for all presented models at four different time steps for the positions of the dust peak and of the methanol peak in Table 3.

Fig. 6 Comparison of the models presented in Chacón-Tanarro et al. (2019) – S16 (obtained with pyRate) and V17 (obtained with MONACO) – to the 1D-4 model (fiducial model) computed with an updated version of pyRate. The column density profiles are shown at different time steps: t = 1.6 × 105 yr for V17, t = 3.0 × 104 yr and t = 3.0 × 105 yr for S16, and t = 3.0 × 105 yr for 1D-4. The observed profiles are depicted in black (errors as grey-shaded areas). The CO column densities were obtained using observations of C17O and the average isotopic ratios of 16O/18O = 557 and 18O/17O = 3.6 in the local interstellar medium (Wilson et al. 1999).

4.1 Chemical variation

4.1.1 1D-1 to 1D-6: Combining enhanced diffusion and abstraction reactions

As already pointed out in Sect. 3.2, to explain the observed magnitudes of molecular abundances and column densities, we had to introduce some type of enhanced diffusion of atoms on the grain surface. In V17 and the fiducial model, this is realised by enabling the diffusion of H (and D) via quantum tunnelling. Laboratory experiments on an amorphous solid water surface (Hama et al. 2012) as well as a CO surface (Kimura et al. 2018) actually suggest that the diffusion of H and D is dominated by thermal hopping, even at temperatures around 10 K. This was concluded since a significant isotope effect, which is expected if the diffusion proceeds mainly via quantum tunnelling, could not be observed. Therefore, we explored both the option for an increased, fast thermal diffusion and the diffusion of H and D atoms by quantum tunnelling.

For the former, a value of the diffusion-to-binding energy, E_d/E_b, has to be chosen. Due to the lack of detailed experimental data, most chemical codes apply only one value for the diffusion-to-binding energy to every surface species. In reality, it is likely that different species have an individual diffusion-to-binding energy that extends over a wider range of values, depending on the various types of potential wells present. Reasonable values of the diffusion-to-binding energy cover the range of 0.2 to 0.7 (Furuya et al. 2022). As a first approximation, we assumed that the comparatively large number of hydrogen and deuterium atoms might fill up the deeper potential wells quite quickly, making the shallower potential wells more relevant for the diffusion process. For that reason, we adopted the lowest debated value of 0.2 for the diffusion-to-binding energy. As a reference, we also tested the option of introducing no enhanced diffusion – neither fast thermal diffusion nor tunnelling diffusion – only employing a typically assumed value for the diffusion-to-binding energy of 0.55, leading to slow thermal hopping. Additionally, we decided to run two models for every explored type of diffusion process: one with only addition reactions and one with addition and abstraction reactions.

Models 1D-1 and 1D-2, employing slow thermal diffusion with a diffusion-to-binding energy E_d/E_b of 0.55, serve as references to the models with enhanced diffusion. Both produce CH₃OH column densities around 10⁷ cm⁻², which is several orders of magnitude lower than the observed value of 10¹³ cm⁻². Models 1D-3 and 1D-4, additionally employing the diffusion via tunnelling for H and D atoms, as well as models 1D-5 and 1D-6, relying on fast thermal hopping, show significantly higher column densities of the order of 10¹² cm⁻². These results deviate only within a factor of 10 from the observed values and are thereby matching the observation. Therefore, we conclude that some form of enhanced diffusion of H and D atoms over the grain surface has to take place in order to reach similar column densities as the ones measured in the pre-stellar core L1544.

The column densities in the models that employ fast thermal hopping are higher for early time steps (t = 1 × 10⁵ yr) and intermediate time steps (t = 5 × 10⁵ yr) – up to a factor of 6 for CH₃OH and up to a factor of 4 for CH₂DOH – but decline earlier than in the models with tunnelling diffusion. Moreover, we find that the CH₃OH column density profiles in the 1D-4 model including only addition reactions are lower (by up to a factor of 3) than for the 1D-3 model comprising both addition and abstraction reactions. The CH₂DOH column density profiles, on the other hand, are higher (by up to a factor of 2.5) in the models with addition reactions only. For the 1D-5 and 1D-6 model, employing fast thermal hopping, this effect is even more pronounced with a factor of up to 10 for CH₃OH and up to eight for CH₂DOH.

In addition to the order of magnitude of column densities, we attempted to reproduce the deuterium fraction N (CH₂DOH)/N (CH₃OH) as closely as possible. This ratio is likely less affected by the individual modelling choices, as the effects on CH₃OH and CH₂DOH are probably similar for most parameter selections. Figure 7 depicts the ratio of the column densities of singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH). Models 1D-1 and 1D-2 show a similar level of deuteration reaching a central value of ≈0.05 for the best-fit time of t = 3.0 × 10⁵yr. 1D-1, the model that includes addition and abstraction reactions, has a slightly higher deuterium fraction than the one with only addition reactions. The slightly higher deuterium fraction is expected theoretically and can be explained by the following: reaction (9) has a much higher activation energy than the competing reaction (10), leading to its isomer. Therefore, the hydrogenation of H₂CO preferably proceeds via the latter reaction: $${\text{H}}_2\text{CO}+\text{H}\stackrel{E_{\text{A}}=5.16\times {10}^3\text{K}}{\to }{\text{CH}}_2\text{OH}$$(9) $${\text{H}}_2\text{CO}+\text{H}\stackrel{E_{\text{A}}=2.00\times {10}^3\text{K}}{\to }{\text{CH}}_3\text{O.}$$(10)

However, the deuteration of CH₃O is not able to produce CH₂DOH in our chemical network, as this would require an exchange of atoms between the two functional groups. It is only able to react to CH₃OD in the presence of D. CH₂DOH can only be produced by deuteration of CH₂OH. Including abstraction reactions into the chemical network opens up another channel for the formation of CH₂DOH via the abstraction of non-deuterated methanol: $${\text{CH}}_3\text{OH}+\text{H/D}\stackrel{E_{\text{A}}=3.62\times {10}^3/3.24\times {10}^3\text{K}}{\to }{\text{CH}}_2\text{OH}+{\text{H}}_2\text{/HD}$$(11) $${\text{CH}}_3\text{OH}+\text{H/D}\stackrel{E_{\text{A}}=5.56\times {10}^3/5.16\times {10}^3\text{K}}{\to }{\text{CH}}_3\text{O}+{\text{H}}_2\text{/HD.}$$(12)

Reaction (11) is favoured in comparison to the analogue including its isomer (12). A similar effect on the deuterium fraction of the molecular abundances was discussed for the static 0D-models in Taquet et al. (2012). They found large enhancements for the N (CH₂DOH)/N (CH₃OH) ratio in the high density case (n_H ≥ 1 × 10⁶ cm⁻³) eventually reaching values above unity, while there was no significant increase for the low density case (n_H = 1 × 10⁴ cm⁻³ − 1 × 10⁵ cm⁻³). Aikawa et al. (2012), on the other hand, could not fully confirm these results with their 1D radiative hydrodynamics model.

In contrast to the models with slow thermal diffusion, the models that include a form of enhanced diffusion – either fast thermal hopping or tunnelling diffusion – show the opposite behaviour when it comes to the inclusion of abstraction reactions (see Fig. 7). The models with only addition, 1D-4 and 1D-6, have an up to 8 times higher N(CH₂DOH)/N (CH₃OH) ratio at certain time steps than the one where addition and abstraction reactions are included. The downside of this increased deuterium fraction is that the increased amount of CH₂DOH apparently causes a decrease in the amount of CH₃OH that is produced.

Models 1D-3 and 1D-5, comprising addition and abstraction reactions, are not able to reproduce the observed N (CH₂ DOH)/N (CH₃OH) ratio as well. They only match the level of the observations for the innermost part of the core at early time steps (t = 1 × 10⁵yr), but decline too quickly with increasing radius. For late time steps (t = 1 × 10⁶ yr), the deuterium fraction is increasing again to a level, where the intermediate radii (r = 2000 AU–r = 6000) match the observational profile. The innermost part, however, shows less deuteration than both the outer ring and the observations.

By looking at the reaction rates of for example model 1D-3, one can see that at a time step close to the peak in methanol formation (≈t = 5 × 10⁵ yr), the magnitude of the addition reaction rates and the one from the abstraction reactions approach each other and basically become almost equal in value, meaning that the net formation of methanol does proceed much more slowly. We therefore suspect that the reaction rates for the abstraction reactions are too high/much higher than in reality.

Having assessed that some form of enhanced diffusion is needed, it is not possible to decide based on our modelling results which of the processes of enhanced diffusion matches the observations better. For times before the time step when the methanol column density peaks, the fiducial model, employing tunnelling diffusion, has a higher N (CH₂ DOH)/N (CH₃OH) ratio profile than the 1D-6 model, relying on thermal hopping. The former is slightly above the observed values, but reproduces the observed shape of the profile very well. The latter is well within the area of uncertainty of the observations. It has, however, a much flatter profile than is observed. The profile of the fiducial model flattens with time, until both models become nearly identical at the best-fit time of 3.0 × 10⁵ yr. After the time step when the methanol column density peaks, the N(CH₂DOH)/N(CH₃OH) ratio of the fiducial model decreases and begins to match the observed values quite closely for time steps beyond t = 5.0 × 10⁵ yr. The N(CH₂DOH)/N(CH₃OH) ratio of the 1D-6 model starts increasing far above the observed values, until it reaches values of almost 0.14 for late times (t = 1 × 10⁶ yr).

Fig. 7 Modelled ratio between singly deuterated methanol (CH2 DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105 yr to 106 yr. The best-fit time is t = 3 × 105 yr. We show two models with slow thermal hopping (Ed/Eb = 0.55; blue lines), two models with slow thermal hopping and tunnelling diffusion of H and D (red lines), and two models with fast thermal hopping (Ed/Eb = 0.2; orange lines). The solid lines indicate models with addition and abstraction reactions, while the dashed lines indicate models with only addition reactions. The black line shows the observed ratio (errors as grey-shaded areas).

Fig. 8 Modelled column density profiles of non-deuterated methanol (CH3OH) and singly deuterated methanol (CH2DOH) for several models at four different time steps ranging from 105yr to 106 yr. We present the varied models that show the largest effects in terms of the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product), and 1D-12 (location-dependent cosmic-ray ionisation rate). The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH3OH column densities, and the dashed lines indicate the CH2DOH column densities.

4.1.2 1D-7: Variation of the grain model

As described in Sect. 2.2, the grain model is a three-phase model, consisting of a gas phase, a chemically active surface phase and a chemically inert mantle phase, for most of the models presented here (see also Table 2). Other modelling tasks performed with pyRate (Sipilä et al. 2016a) show that the choice of the grain model can have a significant effect on the magnitude of the deuterium fraction. Therefore, we decided to run a simulation with a simpler two-phase model, consisting only of a gas phase and a chemically active surface phase. A mantle phase, where frozen-out molecules are stored without chemical alteration is not considered in this model. The 1D-7 model presented in this section is identical to the 1D-4-2 model, but for the inclusion of abstraction reactions for methanol in 1D-7 in contrast to 1D-4-2.

The resulting column density profiles and the N(CH₂DOH)/N(CH₃OH) ratio are depicted in Figs. 8 and 9, respectively. The two phase model, 1D-7, does produce a much lower deuterium fraction than the three-phase model. The deuterium fraction in the two phase model is of the order 10⁻³ to 10⁻⁴ compared to 10⁻² in the fiducial model, in model 1D-3 and in the observed N(CH₂ DOH)/N (CH₃OH) ratio. This finding is consistent with previous results that deuteration on the surface is hindered by using only two phases (Sipilä et al. 2016a). The two-phase model produces higher column densities than model 1D-3 or model 1D-4 for CH₃OH. They are of the order of 10¹³ cm⁻², differing only a factor of a few from the observed values. The column densities of CH₂DOH on the other hand are one or two orders of magnitude lower than in the three phase model and CH₂DOH is thereby severely under-produced.

An explanation for this behaviour is that in two phase models, deuterium atoms on the surface of dust grains tend to get locked into deuterated forms of water, ammonia and methane, forming quite stable bonds that are not easily broken up again. Additionally, the aforementioned molecules are not readily desorbed into the gas phase, which hinders the release of deuterium atoms by the dissociation of gas-phase species. As a consequence, the majority of deuterium is trapped in the molecular ice contents and deuteration of other surface species, including methanol, is suppressed in two phase models compared to the more advanced multilayer models.

4.1.3 1D-8: Variation of the reactive desorption mechanism

The reactive desorption mechanism that was in place in pyRate before the extension to the more advanced treatment, described in Sect. 2.1, allowed only the reactive desorption of exothermic surface reactions with a single reaction product. Additionally, the mechanism applied the same reactive desorption efficiency, typically 1%, to every eligible reaction and did not distinguish between surface types. In fact that is the case for most other reactive desorption mechanisms used in the literature, except for the one presented here and the one in Vasyunin et al. (2017).

The new reactive desorption mechanism extends its application to chemical reactions with two reaction products. It is now interesting to quantify the consequences of this change. We note that the newly developed mechanism yields a larger reactive desorption efficiency for the lighter reaction partner and a lower efficiency for the heavier one, depending on their mass ratio. In model 1D-8, we apply the modified reactive desorption mechanism only to reactions with one product. We anticipate that fewer of the light species are expelled from the surface of the dust grains. Although the reactive desorption mechanism is in place for all exothermic surface reactions, we also expect to hinder a specific effect caused by the existence of the abstraction reactions. For example for reaction (13): $${\text{H}}_2\text{CO}+\text{H}\to \text{HCO}+{\text{H}}_2,$$(13)

the new reactive desorption mechanism yields on a CO-surface an efficiency of 54% for the H₂ and of 1.5 × 10⁻³⁶ % for HCO. The abstraction reactions cause the expulsion of significant amounts of H₂, HD, and D₂, while the desorption of the larger reaction partner is negligible.

The column densities of CH₃OH and CH₂DOH in the 1D-8 model are increased compared to the 1D-3 model (see Fig. 8). The difference between the two models grows with time, from a factor of ≈ 2 at early time steps (t = 1 × 10⁵ yr) to a factor of ≈ 16 at late time steps (t = 1 × 10⁶ yr) for CH₃OH or from ≈ 2.5 to ≈ 12 for CH₂DOH. This results in a N(CH₂DOH)/N(CH₃OH) profile (see Fig. 9) that is very similar in shape to the 1D-3 model, but slightly higher than that of the 1D-3 model around the time steps before the methanol column density peaks in the 1D-3 model, as the CH₂DOH column density increases more quickly than the one of CH₃OH in the 1D-8 model. This could indicate that the determined reactive desorption efficiencies do not describe the physical reality very well, and in particular that our mechanism is overestimating the efficiency with which lighter particles are expelled from the surface. However, light particles, as for example H, H₂ and their deuterated isotopologues, are especially important for the formation of methanol on dust grains, as it proceeds by successive addition of hydrogen and deuterium atoms. Indeed, there are hints that some of the assumptions made to set up the mechanism might not be fulfilled. For example, Fredon et al. (2021) pointed out that the equal distribution of energy into all the degrees of freedom of the reaction product is unlikely to occur, and it is more likely that one or multiple degrees of freedom are favoured against the others. However, since the present work represents the first step to a more advanced treatment, we kept our assumptions simple and as general as possible. Further work could investigate different options for the partitioning of the available reaction enthalpy.

Fig. 9 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105 yr to 106 yr. We present the varied models that show the largest effects in terms of the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product), and 1D-12 (location-dependent cosmic-ray ionisation rate). The black line shows the observed ratio (errors as grey-shaded areas).

4.1.4 1D-9: Variation of the activation energy

The formation of methanol proceeds by the successive hydrogenation of CO along HCO, H₂CO and CH₂OH/CH₃O to CH₃OH. Reactions (14)–(16), $$\text{CO}+\text{H}\stackrel{E_{\text{A}}=1.76\times {10}^3\text{K}}{\to }\text{HCO}$$(14) $${\text{H}}_2\text{CO}+\text{H}\stackrel{E_{\text{A}}=5.16\times {10}^3\text{K}}{\to }{\text{CH}}_2\text{OH}$$(15) $${\text{H}}_2\text{CO\hspace{0.17em}+\hspace{0.17em}H}\stackrel{E_{\text{A}}=2.00\times {10}^3\text{\hspace{0.17em}K}}{\to }{\text{CH}}_3\text{O},$$(16)

have an activation barrier. Their corresponding activation energies E_A are indicated on top of the arrows. The remaining reactions are barrier-less. The reactions leading to deuterated analogues of these species have similar activation energies, at times with somewhat lower values. A complete overview is shown in Appendix B for both addition and abstraction reactions.

The exact values of the activation energies, especially the difference for reactions leading to non-deuterated and deuterated isotopologues, could potentially have a large impact on the N(CH₂DOH)/N(CH₃OH) ratio. Therefore, we explored how the results are affected by a small variation in the activation energy. Specifically, we aimed to test if decreasing the activation energies for reactions producing deuterated isotopologues could lead to a significant increase in the N(CH₂DOH)/N(CH₃OH) ratio. Hence, we decided to decrease the activation energy for those reactions by 200 K. We only varied the activation energy of specific addition reactions (marked with a star in the overview in Appendix B). The abstraction reactions were left untouched.

Undertaking this variation produces a slightly higher CH₂DOH column density (see Fig. C.1). For early times (t = 1 × 10⁵ yr), the 1D-9 model shows a twice as high CH₂DOH column density as compared to the 1D-3 model. The difference between the two models decreases towards the time step when the methanol column density peaks at t = 3.0 × 10⁵ yr. At this time step, the 1D-8 model has a slightly lower CH₂DOH column density profile as compared to the 1D-3 model. After the temporal methanol peak, the difference increases for intermediate (t = 5.0 × 10⁵ yr) and late (t = 1 × 10⁶ yr) time steps to a factor between 1 and 2. The CH₃OH column densities differ by a negligible amount between the two models (see Fig. C.1). Consequently, the N(CH₂DOH)/N(CH₃OH) profile (see Fig. C.2) experiences an increase as well. For early time steps the N(CH₂DOH)/N(CH₃OH) ratio reaches high values of up to 16.6 at t = 1.15 × 10⁵ yr and is for several other time steps well within the area of uncertainty of the observed profile for the very inner centre of the core (0.055–0.091) but is declining significantly more steeply at larger radii. The ratio decreases rapidly at the time when the methanol column density peaks, and increases again at late times (t = 1 × 10⁶ yr), hitting the area of uncertainty, but presenting only a small upward shift compared to the 1D-3 model.

4.1.5 1D-10: Variation of the formation enthalpies

The incorporation of the more sophisticated reactive desorption mechanism required to expand the list of tabulated formation enthalpies H_form, which are necessary to compute the reaction enthalpies ∆H. The complete list of formation enthalpies (and binding energies) of species involved in exothermic surface reactions is shown in Appendix A. Since experimental data for deuterated isotopologues are scarce, we have for the most part adopted the same formation enthalpy values for the deuterated isotopologues as for their non-deuterated counterparts. In a few cases, however, we were able to find experimentally measured values for the deuterated analogues in the NIST Chemistry WebBook⁴. The adopted values are marked with a star in Table A.1. The formation enthalpies of the non-deuterated molecules do not differ strongly from the values of their deuterated analogues. For most species, there is a difference of approximately 4 kJ mol⁻¹ or less for singly deuterated, up to 7 kJ mol⁻¹ for doubly deuterated and up to 13 kJ mol⁻¹ for triply deuterated analogues. This list of formation enthalpies and binding energies was used for all the presented models.

In order to secure that changing the formation enthalpies for only some of the deuterated isotopologues has no significant effect, we ran a model in which we adopted the same values for non-deuterated and deuterated isotopologues. The column density profiles for CH₃OH and CH₂DOH of this model are shown in Fig. C.1 and the deuterium fraction profiles are shown in Fig. C.2. The effect on both the CH₃OH and CH₂DOH column densities and on the N (CH₂ DOH)/N (CH₃OH) ratio is vanishingly small. The difference with respect to the reference model, 1D-3, ranges around 4% for the best-fit time (t = 3.0 × 10⁵ yr.)

4.2 Physical variation

4.2.1 1D-11: Variation of the gas temperature

Based on observations of the NH₃ (1,1) and (2,2) lines and especially their relative strengths, there is reason to believe that the determined gas temperatures in the used physical model by Keto & Caselli (2010) are too high at the intermediate densities, where the maximum of the CO desorption and methanol formation occurs. In principle, lower temperatures should help promote the deuteration process. Therefore, we decided to test a model, in which we decreased the gas temperature throughout the entire core by 1 K, as a first approximation for a revised temperature profile.

The obtained column density profiles for both isotopologues of methanol, non-deuterated and singly deuterated (see Fig. C.1), are quite close to the reference model 1D-3. The CH₃OH column densities of the 1D-3 model are a bit higher until the time step when the methanol column density peaks, which is reversed after the methanol column densities start decreasing again. The CH₂DOH column densities in the 1D-11 model are a little higher than in the 1D-3 model around the methanol peak, but are lower before and almost identical after the peak. The deuterium fraction in the 1D-11 model (see Fig. C.2) increases less quickly, but the maximum value does reach a slightly higher maximum value than the 1D-3 model at a later time step. After the time step of the methanol peak the N(CH₂DOH)/N(CH₃OH) ratios become very similar.

4.2.2 1D-12: Variation of the cosmic-ray ionisation rate

UV photons are not able to penetrate the inner, denser parts of molecular clouds with visual extinctions A_V ≥ 1, as they are already efficiently absorbed by the outer layers of the cloud. Therefore, cosmic rays take their place as the main ionising agents in the central parts, constituting the start of the ion-molecule chemistry. The penetrating cosmic rays will ionise molecular hydrogen to form ${\text{H}}_2^{+}$, which then in turn reacts again with the large reservoir of hydrogen molecules, thereby forming the ${\text{H}}_3^{+}$ ion. This particular ion can react with deuterated molecular hydrogen HD in the following reversible reaction: $${\text{H}}_3^{+}+\text{HD}\rightleftharpoons {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2.$$(17)

The direction of the reaction from left to right is strongly favoured due to the lower zero-point energy of H₂D⁺ for temperatures below 30 K, which pre-stellar cores usually exhibit (strictly true only if all reactants and products are in para form Pagani et al. 1992). Additionally, CO, the main destroyer of ${\text{H}}_3^{+}$ is mostly frozen out on dust grains in the very inner part of the pre-stellar core. These convenient conditions promote a very efficient deuteration process in the inner parts of the core, which also quickly translates the high D/H ratios to more complex molecules.

While we used the canonical value of the cosmic-ray ionisation rate per hydrogen molecule ζ(H) = 1.3 × 10⁻¹⁷ s⁻¹ for most of the models presented here, for the 1D-12 model we used a physical model that is attenuating the cosmic-ray ionisation rate depending on its distance from the centre. The ℒ-model, presented in Padovani et al. (2018) and already tested in the context of the pre-stellar core L1544 by Redaelli et al. (2021b), increases the cosmic-ray ionisation rate from ζ(H) = 2.02 × 10⁻¹⁷ s⁻¹ in the innermost cell to ζ(H) = 4.89 × 10⁻¹⁷ s⁻¹ at the outer boundary of the core.

The 1D-12 model exhibits higher CH₃OH and CH₂DOH column densities as compared to the 1D-3 model (see Fig. 8). Especially for early times (t = 1 × 10⁵ yr), the CH₃OH column densities are a factor of 13 higher than in the 1D-3 model. However, the difference to the 1D-3 model is decreasing over time: at intermediate times (t = 5 × 10⁵ yr) it is approximately a factor of 2 and even lower at late times (t = 1 × 10⁶ yr). The shape of the column density profile (see Fig. 9) is almost identical between the two models, with small deviations at early times. For the increase in the singly deuterated methanol, we see a time-delayed behaviour compared to the non-deuterated isotopo-logue. CH₂DOH is amplified by a factor of 3 at early times (t = 1 × 10⁵ yr), but it never gets to the values observed for CH₃OH. Nevertheless, the amplification for CH₂DOH overtakes the one for CH₃OH after the time step when the methanol column density peaks, which results in a larger deuterium fraction magnitude than in the 1D-3 model for the later time steps. The largest N(CH₂DOH)/N(CH₃OH) ratio of ≈0.74 is reached at late times (t = 1 × 10⁶ yr).

4.2.3 1D-13: Variation of the cosmic-ray desorption mechanism

A frequently adopted model for the CRD is the one laid out in Hasegawa & Herbst (1993). It is also used for all the models presented in this work so far. The model assumes that cosmic rays in the 20−70 MeV nucleon⁻¹ energy range deposit 0.4 MeV of energy into the dust grain, heating it up to a temperature T_max of 70 K. The CRD rate coefficient for molecule i is calculated as the product of the thermal desorption rate, k_max(i, T_max), of i at the temperature T_max and an efficiency term f (a, T_max) for a grain of radius a: $$k_{\text{CR}}\left(i\right)=f\left(a,T_{\max }\right)k_{\text{therm}}\left(i,T_{\max }\right)$$(18)

The efficiency factor f (a, T_max) is determined as the ratio between the cooling time of the grains τ_cool to the heating interval τ_heat. The Hasegawa & Herbst (1993) model adopts constant values for these quantities, for example f (a, T_max) = 1 × 10⁻⁵ s/3.16 × 10¹³ s = 3.16 × 10⁻¹⁹ for a grain of 0.1 µm. A revised version of CRD presented by Sipilä et al. (2021) refines the description of the process by introducing two major modifications to the established scheme. On one hand, the grain cooling time τ_cool is determined now by a dynamic mechanism taking into account the individual sublimation rates of the surface molecules as a function of their time-dependent ice abundances. On the other hand, several different cosmic-ray fluxes can be considered for the calculation of the heating intervals τ_heat.

In order to test how this new mechanism affects the formation of methanol and its deuterated isotopologues, we chose the cosmic-ray flux presented in Léger et al. (1985), as this is the one most consistent with the canonical value of the cosmic-ray ion-isation rate ζ(H) = 1.3 × 10⁻¹⁷ s⁻¹, which we used for the other models. The impact on singly and non-deuterated methanol formation seems to be minor compared to the 1D-3 model. This results fits well with the finding that the desorption of species involved in the methanol formation scheme is dominated by reactive desorption, rather than CRD. The CH₃OH column density profile (see Fig. C.1) is slightly decreased, especially for the earlier time steps, while the CH₂DOH profile is not significantly changed, resulting in a little higher N(CH₂DOH)/N(CH₃OH) ratio (see Fig. C.2). It reaches its highest value of 0.088 at t = 1 × 10⁵ yr in the very inner centre. However, the decline is again much steeper than is observed.

5 Conclusion

We have presented several models for predicting column densities and deuterium fractions of methanol and its deuterated isotopologues in pre-stellar cores. As a comparison to observed quantities, we used single-dish observations of H₂CO and CH₃OH and some of their deuterated isotopologues towards the pre-stellar core L1544 conducted and analysed by Chacón-Tanarro et al. (2019).

All our models use a novel treatment of reactive desorp-tion of molecules from the surface of interstellar dust grains. The treatment is experimentally justified (Minissale et al. 2016b) and derives an individual reactive desorption efficiency for every species, depending on the forming chemical reaction(s) and the type of underlying surface.

The results of our fiducial model were compared to the results of the models V17 (MONACO) and S16 (pyRate) presented in Chacón-Tanarro et al. (2019). The fiducial model includes thermal diffusion (diffusion-to-binding energy ratio E_d/E_b = 0.55) as well as the diffusion of hydrogen and deuterium atoms via quantum tunnelling. The chemical network does not comprise abstraction reactions for the methanol reaction scheme. We estimated a best-fit time that coincides with the occurrence of the methanol peak at 3.0 × 10⁵ yr. At this time step, we find a better agreement with the observations than for the S16 model. The column densities of CH₃OH and CH₂DOH are still underestimated. However, instead of a more than two orders of magnitude deviation, we are able to reduce the difference to approximately an order of magnitude. This improvement is not possible without increasing the diffusion rate on the surface of the dust grain. The observed N(CH₂DOH)/N(CH₃OH) can be reproduced quite closely. Additionally, we find that increasing the number of layers in the chemically active surface phase from one to four, which allows the atoms to also diffuse in the vertical direction as in the V17 model, increases the column densities of non-deuterated and singly deuterated methanol as well as their D/H ratio in the time frame in question.

Previous work by Vasyunin et al. (2017) and Chacón-Tanarro et al. (2019) shows that to reproduce the observed order of magnitude for the methanol abundances, it is necessary to both employ an increased rate of surface diffusion and to disregard abstraction reactions from the reaction scheme. Therefore, we also explored various types of diffusion processes: slow thermal hopping (E_d/E_b = 0.55), fast thermal hopping (E_d/E_b = 0.2), and slow thermal hopping (E_d/E_b = 0.55) combined with the diffusion of H and D atoms via quantum tunnelling. From these tests, we conclude that a form of enhanced diffusion over the surface needs to take place to explain the observational results. Only employing slow thermal hopping produces CH₃OH and CH₂DOH column densities that are several orders of magnitudes below the observed ones. However, we cannot decide based solely on our models which enhanced diffusion process – fast thermal hopping or tunnelling diffusion – matches the observed N(CH₂DOH)/N(CH₃OH) ratio profiles better. The two produce D/H ratios of a similar level that are in good agreement with the observed values. For the best-fit time, the two options have nearly identical profiles.

We also tested both options, employing either addition and abstraction reactions or only addition reactions, for every explored diffusion process. In general, we conclude that including abstraction reactions following Hidaka et al. (2009), in combination with an increased diffusion rate, leads to N(CH₂DOH)/N(CH₃OH) ratios that are a factor of a few lower than in the models that only include addition reactions. We ascribe this behaviour to the fact that the reaction rates of the abstraction reactions become comparable to the addition reactions when combined with enhanced diffusion processes.

Furthermore, we explored other modifications to our model that we suspected would have an effect on the N(CH₂DOH)/N(CH₃OH) ratio:

• A two-phase model resulted in higher CH₃OH and CH₂DOH column densities; the deuterium fraction is severely underestimated, by more than one order of magnitude;

• Reactive desorption applied only to reactions with one reaction product resulted in higher CH₃OH and CH₂DOH column densities and thus a slightly better agreement with the observations;

• A location-dependent cosmic-ray ionisation rate resulted in higher CH₃OH and CH₂DOH column densities that only differ from the observations by a factor of 2 to 3.

Only small effects on the CH₃OH and CH₂DOH column density and deuterium fraction profiles were found for the following model changes:

• A decrease in the activation energy by 200 K leading to deuterated isotopologues (as opposed to non-deuterated species);

• Inclusion of individual formation enthalpies for some deuter-ated isotopologues as opposed to using the same formation enthalpies for hydrogenated and deuterated isotopologues;

• Decrease in the gas temperature by 1 K throughout the entire core;

• Refinement of the used cosmic-ray desorption mechanism following Sipilä et al. (2021).

Further work needs to be carried out on quantifying the reactive desorption mechanism. On the one hand, there is reason to question the assumption of equal partitioning of energy into all the degrees of freedom. Also, some of our results could hint at the fact that light particles are desorbed too easily into the gas phase with the employed mechanism. On the other hand, it would be interesting to investigate how the reactive desorption mechanism influences other species that are formed on the surface of dust grains.

Additionally, a closer look into the intricacies of the surface diffusion processes, especially of H and D, is needed. Particularly interesting would be to explore the effects of introducing different types of potential wells in which species can be trapped. Other reaction mechanisms, such as the Eley-Rideal mechanism, which is not treated by many chemical codes, or non-diffusive chemistry could play an important role in the formation and deuteration of methanol. These mechanisms will be the subject of a future paper.

We conclude that to obtain a reasonable match with the observational column density and deuterium fraction profiles, it is necessary to employ a form of enhanced diffusion process – either fast thermal hopping or diffusion via quantum tunnelling. Furthermore, the inclusion of abstraction reactions in the methanol formation scheme while also using a fast diffusion process leads to N(CH₂DOH)/N(CH₃OH) ratios that are a factor of a few lower than without the abstraction reactions.

Acknowledgements

The work by A.V. is supported by the Russian Ministry of Science and Higher Education via the Project FEUZ-2020-0038.

Appendix A Formation enthalpies and binding energies

In Table A.1 we present the formation enthalpies H_form and binding energies E_b for all species involved in exothermic surface reactions, to which we applied the newly developed reactive desorption mechanism laid out in Sect. 2.1.

Appendix B Applied reaction scheme for the formation of methanol and deuterated isotopologues

Figure B.1 depicts the reaction scheme for the formation of methanol that we employed for all the presented models.

Fig. B.1 Reaction scheme for the formation of CH3OH and its deuterated isotopologues. The upper half shows the applied ‘forward reactions’ (addition reactions), either adding H (in the horizontal direction) or D (in the vertical direction). The lower half shows the employed ‘backward reactions’ (abstraction reactions), reacting with H or D and thereby removing a H2, HD, or D2 molecule. The values on top of the arrows indicate, if existing, the activation energies, EA, in 103 K. There are always two values for the abstraction reactions. The first one gives the value for the reaction with H, the second for the reaction with D. Note that our reaction scheme also includes two ‘substitution reactions’, exchanging one D atom for a H atom in H2CO and HDCO, which are not shown here for the sake of clarity.

Appendix C Column density and deuterium fraction profiles

Figures C.1 and C.2 respectively show the modelled column density and deuterium fraction profiles of the models with only small effects. We present the models 1D-9 (decreased activation energy for deuterated isotopologues), 1D-10 (individual formation enthalpies for deuterated isotopologues), 1D-11 (decrease in the gas temperature) and 1D-13 (refinement of cosmic-ray desorption mechanism).

Fig. C.1 Modelled column density profiles of non-deuterated methanol (CH3OH) and singly deuterated methanol (CH2DOH) for several models at four different time steps ranging from 105yr to 106yr. The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH3OH column densities, and the dashed lines indicate the CH2DOH column densities.

Fig. C.2 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105yr to 106yr. The black line shows the observed ratio (errors as grey-shaded areas).

References

All Tables

All Figures

	Fig. 1 Reaction scheme for the formation of non-deuterated CH3OH by successive hydrogenation. Hydrogen molecules can also be segregated from CH3OH or its precursors via abstraction reactions.
In the text

	Fig. 2 Physical model developed in Keto & Caselli (2010) that yields static radial profiles of the H2 density, n(H2) (blue, in logarithmic scale), the visual extinction, AV (red), the gas temperature, Tgas (orange), and the dust temperature, Tdust (green).
In the text

	Fig. 3 Gas-phase abundance profiles of H2CO, CH3OH, and CH2DOH for the 1D-4 model for four different time steps ranging from 105 yr (left) to 106 yr (right). The best-fit time is t = 3 × 105 yr.
In the text

Fig. 4 Column density profiles of H2CO, CH3OH, and CH2DOH for the 1D-4 model for four different time steps ranging from 105 yr to 106 yr. The lines with markers show the modelled results, integrated along the line of sight and convolved with a 30″ beam. The solid lines show the observationally obtained column density profiles obtained by Chacón-Tanarro et al. (2019) by taking a cut through the dust and methanol peaks. The grey-shaded areas indicate the error bars of the column densities. The position of the dust peak is at r = 0 AU, while r > 0 AU is the direction towards the methanol peak.

In the text

Fig. 5 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for the 1D-4 model for four different time steps ranging from 105 yr (top left) to 106 yr (bottom right). Additionally, we show two variations of this model: one with four layers in the chemically active surface phase instead of one (1D-4-1), and one with a two-phase model (1D-4-2; see also Table 2 and Sect. 3.2 for a more detailed explanation). The coloured lines show the column density ratio of the models, while the black line indicates the observed ratio (errors as grey-shaded areas).

In the text

Fig. 6 Comparison of the models presented in Chacón-Tanarro et al. (2019) – S16 (obtained with pyRate) and V17 (obtained with MONACO) – to the 1D-4 model (fiducial model) computed with an updated version of pyRate. The column density profiles are shown at different time steps: t = 1.6 × 105 yr for V17, t = 3.0 × 104 yr and t = 3.0 × 105 yr for S16, and t = 3.0 × 105 yr for 1D-4. The observed profiles are depicted in black (errors as grey-shaded areas). The CO column densities were obtained using observations of C17O and the average isotopic ratios of 16O/18O = 557 and 18O/17O = 3.6 in the local interstellar medium (Wilson et al. 1999).

In the text

Fig. 7 Modelled ratio between singly deuterated methanol (CH2 DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105 yr to 106 yr. The best-fit time is t = 3 × 105 yr. We show two models with slow thermal hopping (Ed/Eb = 0.55; blue lines), two models with slow thermal hopping and tunnelling diffusion of H and D (red lines), and two models with fast thermal hopping (Ed/Eb = 0.2; orange lines). The solid lines indicate models with addition and abstraction reactions, while the dashed lines indicate models with only addition reactions. The black line shows the observed ratio (errors as grey-shaded areas).

In the text

Fig. 8 Modelled column density profiles of non-deuterated methanol (CH3OH) and singly deuterated methanol (CH2DOH) for several models at four different time steps ranging from 105yr to 106 yr. We present the varied models that show the largest effects in terms of the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product), and 1D-12 (location-dependent cosmic-ray ionisation rate). The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH3OH column densities, and the dashed lines indicate the CH2DOH column densities.

In the text

Fig. 9 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105 yr to 106 yr. We present the varied models that show the largest effects in terms of the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product), and 1D-12 (location-dependent cosmic-ray ionisation rate). The black line shows the observed ratio (errors as grey-shaded areas).

In the text

Fig. B.1 Reaction scheme for the formation of CH3OH and its deuterated isotopologues. The upper half shows the applied ‘forward reactions’ (addition reactions), either adding H (in the horizontal direction) or D (in the vertical direction). The lower half shows the employed ‘backward reactions’ (abstraction reactions), reacting with H or D and thereby removing a H2, HD, or D2 molecule. The values on top of the arrows indicate, if existing, the activation energies, EA, in 103 K. There are always two values for the abstraction reactions. The first one gives the value for the reaction with H, the second for the reaction with D. Note that our reaction scheme also includes two ‘substitution reactions’, exchanging one D atom for a H atom in H2CO and HDCO, which are not shown here for the sake of clarity.

In the text

	Fig. C.1 Modelled column density profiles of non-deuterated methanol (CH3OH) and singly deuterated methanol (CH2DOH) for several models at four different time steps ranging from 105yr to 106yr. The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH3OH column densities, and the dashed lines indicate the CH2DOH column densities.
In the text

	Fig. C.2 Modelled ratio between singly deuterated methanol (CH2DOH) and non-deuterated methanol (CH3OH) for several models at four different time steps ranging from 105yr to 106yr. The black line shows the observed ratio (errors as grey-shaded areas).
In the text