# Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544

W. Riedel, O. Sipilä, E. Redaelli, P. Caselli, A. Vasyunin, Francois Dulieu, N.

Watanabe

► To cite this version:

W. Riedel, O. Sipilä, E. Redaelli, P. Caselli, A. Vasyunin, et al., Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544. Astronomy and Astrophysics - A&A, 2023, 680, pp.A87. 10.1051/0004-6361/202245367. [hal-04626717](https://hal.science/hal-04626717)

HAL Id: hal-04626717

https://hal.science/hal-04626717

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544

W. Riedel, O. Sipilä, E. Redaelli, P. Caselli, A. Vasyunin, Francois Dulieu, N.

Watanabe

► To cite this version:

W. Riedel, O. Sipilä, E. Redaelli, P. Caselli, A. Vasyunin, et al. Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544. Astronomy & Astrophysics - A&A, 2023, 680, pp.A87. 10.1051/0004-6361/202245367. hal-04626717

HAL Id: hal-04626717

https://hal.science/hal-04626717v1

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544

W. Riedel<fnref n="1" />, O. Sipilä<fnref n="1" />, E. Redaelli<fnref n="1" />, P. Caselli<fnref n="1" />, A. I. Vasyunin<fnref n="2" />, F. Dulieu<fnref n="3" />, and N. Watanabe<fnref n="4" />

Received 3 November 2022 / Accepted 11 October 2023

## ABSTRACT

**Context.** In the extremely cold and dark environments of pre-stellar cores, methanol is formed on the surface of interstellar dust grains and released into the gas phase via non-thermal desorption mechanisms. Gaseous methanol constitutes the starting point for the formation of many massive complex organic molecules and is therefore of utmost importance for the build-up of chemical complexity.

**Aims.** We aim to improve upon a previous model for the prediction of column densities and deuterium fractions of non-deuterated and singly deuterated methanol. Thereby, we try to identify crucial chemical and physical parameters for which the study of deuteration could provide valuable additional constraints.

**Methods.** We employed a gas-grain chemical code to devise a model that is in agreement with the observed column density and deuterium fraction profiles of the innermost region of the pre-stellar core L1544. For that purpose, we developed a new treatment of reactive desorption, deriving an individual reactive desorption efficiency for every product species in a chemical reaction that depends on the reaction enthalpy and type of the underlying surface. Furthermore, we explored several options to promote the diffusion of hydrogen and deuterium atoms over the surface of interstellar dust grains in order to increase methanol formation.

**Results.** Our fiducial model employs diffusion via the quantum tunnelling of hydrogen and deuterium atoms, resulting in CH₃OH and CH₂DOH column densities that are approximately an order of magnitude lower than the observed values, which is an improvement over the results of the previous model by a factor of 10. The N(CH₂DOH)/N(CH₃OH) ratio is reproduced within a factor of 1.2 for the centre and 1.8 for the position of the methanol peak. Given the large uncertainties that chemical models typically have, we consider our predictions to be in agreement with the observations. In general, we conclude that a diffusion process with a high diffusion rate needs to be employed to obtain methanol column densities that are in accordance with the observed values. Also, we find that the introduction of abstraction reactions into the methanol formation scheme suppresses deuteration when used in combination with a high diffusion rate.

**Key words.** ISM: abundances – ISM: clouds – ISM: molecules – astrochemistry – stars: formation – ISM: individual objects: L1544

## 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

<footnote n="1">¹ Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching bei München, Germany
e-mail: riedel@mpe.mpg.de</footnote>

<footnote n="2">² Ural Federal University, 620002, 19 Mira street, Yekaterinburg, Russia</footnote>

<footnote n="3">³ CY Cergy Paris Université, Sorbonne Université, Observatoire de Paris, PSL University, CNRS, LERMA, 95000 Cergy, France</footnote>

<footnote n="4">⁴ Institute of Low Temperature Science, Hokkaido University, Sapporo, Hokkaido 060-0819, Japan</footnote>
# Modelling deuterated isotopologues of methanol towards the pre-stellar core L1544

W. Riedel<fnref n="1" />, O. Sipilä<fnref n="1" />, E. Redaelli<fnref n="1" />, P. Caselli<fnref n="1" />, A. I. Vasyunin<fnref n="2" />, F. Dulieu<fnref n="3" />, and N. Watanabe<fnref n="4" />

Received 3 November 2022 / Accepted 11 October 2023

## ABSTRACT

**Context.** In the extremely cold and dark environments of pre-stellar cores, methanol is formed on the surface of interstellar dust grains and released into the gas phase via non-thermal desorption mechanisms. Gaseous methanol constitutes the starting point for the formation of many massive complex organic molecules and is therefore of utmost importance for the build-up of chemical complexity.

**Aims.** We aim to improve upon a previous model for the prediction of column densities and deuterium fractions of non-deuterated and singly deuterated methanol. Thereby, we try to identify crucial chemical and physical parameters for which the study of deuteration could provide valuable additional constraints.

**Methods.** We employed a gas-grain chemical code to devise a model that is in agreement with the observed column density and deuterium fraction profiles of the innermost region of the pre-stellar core L1544. For that purpose, we developed a new treatment of reactive desorption, deriving an individual reactive desorption efficiency for every product species in a chemical reaction that depends on the reaction enthalpy and type of the underlying surface. Furthermore, we explored several options to promote the diffusion of hydrogen and deuterium atoms over the surface of interstellar dust grains in order to increase methanol formation.

**Results.** Our fiducial model employs diffusion via the quantum tunnelling of hydrogen and deuterium atoms, resulting in CH₃OH and CH₂DOH column densities that are approximately an order of magnitude lower than the observed values, which is an improvement over the results of the previous model by a factor of 10. The N(CH₂DOH)/N(CH₃OH) ratio is reproduced within a factor of 1.2 for the centre and 1.8 for the position of the methanol peak. Given the large uncertainties that chemical models typically have, we consider our predictions to be in agreement with the observations. In general, we conclude that a diffusion process with a high diffusion rate needs to be employed to obtain methanol column densities that are in accordance with the observed values. Also, we find that the introduction of abstraction reactions into the methanol formation scheme suppresses deuteration when used in combination with a high diffusion rate.

**Key words.** ISM: abundances – ISM: clouds – ISM: molecules – astrochemistry – stars: formation – ISM: individual objects: L1544

## 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

<footnote n="1">¹ Max-Planck-Institut für extraterrestrische Physik, Gießenbachstraße 1, 85748 Garching bei München, Germany
e-mail: riedel@mpe.mpg.de</footnote>

<footnote n="2">² Ural Federal University, 620002, 19 Mira street, Yekaterinburg, Russia</footnote>

<footnote n="3">³ CY Cergy Paris Université, Sorbonne Université, Observatoire de Paris, PSL University, CNRS, LERMA, 95000 Cergy, France</footnote>

<footnote n="4">⁴ Institute of Low Temperature Science, Hokkaido University, Sapporo, Hokkaido 060-0819, Japan</footnote>
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$_3$OH and CH$_2$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$_3$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, 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$_3$OH). However, the column density profile for singly deuterated methanol (CH$_2$DOH) had to be obtained using the N(CH$_2$DOH)/N(CH$_3$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}{\epsilon \Delta H}\right), \qquad (1)$$

where $\Delta 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 $\epsilon$ 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_{\text{atoms}}$, with $n_{\text{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$_3$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, $\epsilon$, received by the reaction product with mass $m$ when colliding with a surface element with an effective mass $M$:

$$\epsilon = \left( \frac{M - m}{M + m} \right)^2. \qquad (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$_2$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
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$_3$OH and CH$_2$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$_3$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, 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$_3$OH). However, the column density profile for singly deuterated methanol (CH$_2$DOH) had to be obtained using the N(CH$_2$DOH)/N(CH$_3$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}{\epsilon \Delta H}\right), \qquad (1)$$

where $\Delta 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 $\epsilon$ 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_{\text{atoms}}$, with $n_{\text{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$_3$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, $\epsilon$, received by the reaction product with mass $m$ when colliding with a surface element with an effective mass $M$:

$$\epsilon = \left( \frac{M - m}{M + m} \right)^2. \qquad (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$_2$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
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

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 H2O 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. N2, H2CO and CH3OH), 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, H2O 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:

$$RD_{\text{tot}}(i, t) = \sum_{j} RD_{j}(i) \cdot \frac{n_{j}^{*}(t)}{n_{\text{tot}}^{*}(t)}, \qquad (3)$$

where RDj is the individual reactive desorption efficiency for either bare grain H2O or CO, n*j the number of surface sites inhabited by surface type j and n*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, Ekin, received by a product species i:

$$E_{\text{kin}}(i) = E_{\text{kin}} \cdot \frac{m_j}{m_i + m_j} = \frac{N_{\text{trans}}}{N_{\text{tot}}} \cdot \Delta H \cdot \frac{m_j}{m_i + m_j}. \qquad (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: Ntrans = 3nprod, where nprod 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: Ntot = Σi 3natoms,i, with natoms,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:

$$RD_i = \exp\left(-\frac{1}{\epsilon}\left(\frac{E_b}{\Delta H}\right)\left(\frac{1}{3}\frac{N_{\text{trans}}}{N_{\text{tot}}}\frac{m_j}{m_i+m_j}\right)^{-1}\right). \qquad (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

$$H_2CO + H \rightarrow HCO + H_2, \qquad (6)$$

with HCO being approximately 14 times more massive than H2, 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 H2 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 H2 and 0.1% for the HCO molecule. However, there are no similar reactions efficient at low temperatures in the Vasyunin et al. (2017) model.

## 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, 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 H2, H2+ and H3+ 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
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

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 H2O 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. N2, H2CO and CH3OH), 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, H2O 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:

$$RD_{\text{tot}}(i, t) = \sum_{j} RD_{j}(i) \cdot \frac{n_{j}^{*}(t)}{n_{\text{tot}}^{*}(t)}, \qquad (3)$$

where RDj is the individual reactive desorption efficiency for either bare grain H2O or CO, n*j the number of surface sites inhabited by surface type j and n*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, Ekin, received by a product species i:

$$E_{\text{kin}}(i) = E_{\text{kin}} \cdot \frac{m_j}{m_i + m_j} = \frac{N_{\text{trans}}}{N_{\text{tot}}} \cdot \Delta H \cdot \frac{m_j}{m_i + m_j}. \qquad (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: Ntrans = 3nprod, where nprod 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: Ntot = Σi 3natoms,i, with natoms,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:

$$RD_i = \exp\left(-\frac{1}{\epsilon}\left(\frac{E_b}{\Delta H}\right)\left(\frac{1}{3}\frac{N_{\text{trans}}}{N_{\text{tot}}}\frac{m_j}{m_i+m_j}\right)^{-1}\right). \qquad (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

$$H_2CO + H \rightarrow HCO + H_2, \qquad (6)$$

with HCO being approximately 14 times more massive than H2, 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 H2 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 H2 and 0.1% for the HCO molecule. However, there are no similar reactions efficient at low temperatures in the Vasyunin et al. (2017) model.

## 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, 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 H2, H2+ and H3+ 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
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Physical model developed in Keto & Caselli (2010) that yields static radial profiles of the H₂ density, *n*(*H*₂) (blue, in logarithmic scale), the visual extinction, *A*<sub>V</sub> (red), the gas temperature, *T*<sub>gas</sub> (orange), and the dust temperature, *T*<sub>dust</sub> (green).</figcaption></figure>

<table><caption>**Table 1.** Initial chemical abundances with respect to *n*<sub>H</sub>.</caption><thead><tr><th>Species</th><th>Initial abundance</th><th>Species</th><th>Initial abundance</th></tr></thead><tbody><tr><td>He</td><td>9.00 × 10<sup>-2</sup></td><td>S<sup>+</sup></td><td>8.00 × 10<sup>-8</sup></td></tr><tr><td>pH<sub>2</sub></td><td>5.00 × 10<sup>-1</sup></td><td>Si<sup>+</sup></td><td>8.00 × 10<sup>-9</sup></td></tr><tr><td>oH<sub>2</sub></td><td>5.00 × 10<sup>-4</sup></td><td>Na<sup>+</sup></td><td>2.00 × 10<sup>-9</sup></td></tr><tr><td>HD</td><td>1.60 × 10<sup>-5</sup></td><td>Mg<sup>+</sup></td><td>7.00 × 10<sup>-9</sup></td></tr><tr><td>H</td><td>1.00 × 10<sup>-8</sup></td><td>Fe<sup>+</sup></td><td>3.00 × 10<sup>-9</sup></td></tr><tr><td>D</td><td>1.00 × 10<sup>-8</sup></td><td>P<sup>+</sup></td><td>2.00 × 10<sup>-10</sup></td></tr><tr><td>C<sup>+</sup></td><td>1.20 × 10<sup>-4</sup></td><td>Cl<sup>+</sup></td><td>1.00 × 10<sup>-9</sup></td></tr><tr><td>N</td><td>7.60 × 10<sup>-5</sup></td><td>F</td><td>2.00 × 10<sup>-9</sup></td></tr><tr><td>O</td><td>2.56 × 10<sup>-4</sup></td><td></td><td></td></tr></tbody></table>

**Notes.** (a) From Semenov et al. (2010).

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*<sub>b</sub> and formation enthalpies *H*<sub>form</sub> is crucial. The values of *E*<sub>b</sub> and *H*<sub>form</sub> 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<fnref n="1" />.

## 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*<sub>gas</sub>, the dust temperature *T*<sub>dust</sub> and the visual extinction *A*<sub>V</sub> 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

<footnote n="1">¹ https://webbook.nist.gov/chemistry</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Physical model developed in Keto & Caselli (2010) that yields static radial profiles of the H₂ density, *n*(*H*₂) (blue, in logarithmic scale), the visual extinction, *A*<sub>V</sub> (red), the gas temperature, *T*<sub>gas</sub> (orange), and the dust temperature, *T*<sub>dust</sub> (green).</figcaption></figure>

<table><caption>**Table 1.** Initial chemical abundances with respect to *n*<sub>H</sub>.</caption><thead><tr><th>Species</th><th>Initial abundance</th><th>Species</th><th>Initial abundance</th></tr></thead><tbody><tr><td>He</td><td>9.00 × 10<sup>-2</sup></td><td>S<sup>+</sup></td><td>8.00 × 10<sup>-8</sup></td></tr><tr><td>pH<sub>2</sub></td><td>5.00 × 10<sup>-1</sup></td><td>Si<sup>+</sup></td><td>8.00 × 10<sup>-9</sup></td></tr><tr><td>oH<sub>2</sub></td><td>5.00 × 10<sup>-4</sup></td><td>Na<sup>+</sup></td><td>2.00 × 10<sup>-9</sup></td></tr><tr><td>HD</td><td>1.60 × 10<sup>-5</sup></td><td>Mg<sup>+</sup></td><td>7.00 × 10<sup>-9</sup></td></tr><tr><td>H</td><td>1.00 × 10<sup>-8</sup></td><td>Fe<sup>+</sup></td><td>3.00 × 10<sup>-9</sup></td></tr><tr><td>D</td><td>1.00 × 10<sup>-8</sup></td><td>P<sup>+</sup></td><td>2.00 × 10<sup>-10</sup></td></tr><tr><td>C<sup>+</sup></td><td>1.20 × 10<sup>-4</sup></td><td>Cl<sup>+</sup></td><td>1.00 × 10<sup>-9</sup></td></tr><tr><td>N</td><td>7.60 × 10<sup>-5</sup></td><td>F</td><td>2.00 × 10<sup>-9</sup></td></tr><tr><td>O</td><td>2.56 × 10<sup>-4</sup></td><td></td><td></td></tr></tbody></table>

**Notes.** (a) From Semenov et al. (2010).

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*<sub>b</sub> and formation enthalpies *H*<sub>form</sub> is crucial. The values of *E*<sub>b</sub> and *H*<sub>form</sub> 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<fnref n="1" />.

## 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*<sub>gas</sub>, the dust temperature *T*<sub>dust</sub> and the visual extinction *A*<sub>V</sub> 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

<footnote n="1">¹ https://webbook.nist.gov/chemistry</footnote>
<table><caption>**Table 2.** Overview of the various models investigated in this work.</caption><thead><tr><th>Model</th><th>Tunnelling diffusion</th><th>Thermal diffusion</th><th>E<sub>d</sub>/E<sub>b</sub><sup>(a)</sup></th><th>Abstraction reactions</th><th>Gas-grain model</th><th>Cosmic-ray ionisation rate</th><th>Other modifications</th></tr></thead><tbody><tr><td colspan="8"><strong>Chemical variation</strong></td></tr><tr><td>1D-1<sup>(†)</sup></td><td>✗</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-2<sup>(*)</sup></td><td>✗</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-3</td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-4</td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-4-1<sup>(*)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>4 active surface layers</td></tr><tr><td>1D-4-2<sup>(*)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>2 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-5<sup>(†)</sup></td><td>✗</td><td>✓</td><td>0.2</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-6<sup>(*)</sup></td><td>✗</td><td>✓</td><td>0.2</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-7<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>2 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-8<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>No reactive desorption for reactions with 2 products</td></tr><tr><td>1D-9<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>Reduced E<sub>A</sub><sup>(b)</sup> for deuterated species by 200 K</td></tr><tr><td>1D-10<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>H<sub>form</sub>(XH) = H<sub>form</sub>(XD)<sup>(c)</sup> (see Table A.1)</td></tr><tr><td colspan="8"><strong>Physical variation</strong></td></tr><tr><td>1D-11<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>Decreased T<sub>gas</sub> by 1 K throughout the core</td></tr><tr><td>1D-12<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>L-model<sup>(d)</sup></td><td></td></tr><tr><td>1D-13<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td></td><td>Dynamic mechanism for CR- desorption<sup>(e)</sup></td></tr></tbody></table>

**Notes.** Only the chemical and physical properties that vary between models are listed. The 1D-4 model is the fiducial model presented in Sect. 3.1. <sup>(*)</sup> Models that take the 1D-4 model as starting point. <sup>(†)</sup> Models that take the 1D-3 model, including the abstraction reactions, as a starting point. <sup>(a)</sup> *E*<sub>d</sub> is the diffusion energy and *E*<sub>b</sub> the binding energy. <sup>(b)</sup> *E*<sub>A</sub> is the activation energy of the reaction. <sup>(c)</sup> *H*<sub>form</sub>(XH) is the formation enthalpy of a non-deuterated species and *H*<sub>form</sub>(XD) is the formation enthalpy of a deuterated species. <sup>(d)</sup> From Padovani et al. (2018). <sup>(e)</sup> From Sipilä et al. (2021).

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 \times 10^5$ yr, the CH<sub>3</sub>OH abundance reaches its maximum value of approximately $4 \times 10^{-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 \times 10^4$ 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<sub>3</sub>OH column density profile with an amplitude of $N \approx 1 \times 10^{11}$ cm<sup>-2</sup>, under-producing the observed values for methanol of $3.9 \times 10^{13}$ cm<sup>-2</sup> at the dust peak and $5.9 \times 10^{13}$ cm<sup>-2</sup> at the methanol peak by roughly two orders of magnitude. On the other hand, the N(CH<sub>2</sub>DOH)/N(CH<sub>3</sub>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<sub>2</sub>DOH + H<sub>3</sub>O<sup>+</sup> → CH<sub>3</sub>OHD<sup>+</sup> + H<sub>2</sub>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 $\chi^2$ 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 $\chi^2$ 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 \times 10^5$ yr, which is roughly consistent with estimations from other molecules (e.g. Redaelli et al. 2019a, 2021b). Additionally, it coincides with the occurrence of the highest methanol column density in the probed time frame. For $t = 3.0 \times 10^5$ yr, the fiducial model reaches a column density of methanol of ≈1.5 × 10<sup>12</sup> cm<sup>-2</sup>, 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, 2008, Wakelam et al. 2010, and references therein).

Figure 5 compares the modelled deuterium fraction N(CH<sub>2</sub>DOH)/N(CH<sub>3</sub>OH) with the observed one. At the early time step ($t = 1 \times 10^5$ 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
<table><caption>**Table 2.** Overview of the various models investigated in this work.</caption><thead><tr><th>Model</th><th>Tunnelling diffusion</th><th>Thermal diffusion</th><th>E<sub>d</sub>/E<sub>b</sub><sup>(a)</sup></th><th>Abstraction reactions</th><th>Gas-grain model</th><th>Cosmic-ray ionisation rate</th><th>Other modifications</th></tr></thead><tbody><tr><td colspan="8"><strong>Chemical variation</strong></td></tr><tr><td>1D-1<sup>(†)</sup></td><td>✗</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-2<sup>(*)</sup></td><td>✗</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-3</td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-4</td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-4-1<sup>(*)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>4 active surface layers</td></tr><tr><td>1D-4-2<sup>(*)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✗</td><td>2 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-5<sup>(†)</sup></td><td>✗</td><td>✓</td><td>0.2</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-6<sup>(*)</sup></td><td>✗</td><td>✓</td><td>0.2</td><td>✗</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-7<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>2 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td></td></tr><tr><td>1D-8<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>No reactive desorption for reactions with 2 products</td></tr><tr><td>1D-9<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>Reduced E<sub>A</sub><sup>(b)</sup> for deuterated species by 200 K</td></tr><tr><td>1D-10<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>H<sub>form</sub>(XH) = H<sub>form</sub>(XD)<sup>(c)</sup> (see Table A.1)</td></tr><tr><td colspan="8"><strong>Physical variation</strong></td></tr><tr><td>1D-11<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>1.3 × 10<sup>-17</sup> s<sup>-1</sup></td><td>Decreased T<sub>gas</sub> by 1 K throughout the core</td></tr><tr><td>1D-12<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td>L-model<sup>(d)</sup></td><td></td></tr><tr><td>1D-13<sup>(†)</sup></td><td>✓</td><td>✓</td><td>0.55</td><td>✓</td><td>3 phases</td><td></td><td>Dynamic mechanism for CR- desorption<sup>(e)</sup></td></tr></tbody></table>

**Notes.** Only the chemical and physical properties that vary between models are listed. The 1D-4 model is the fiducial model presented in Sect. 3.1. (*) Models that take the 1D-4 model as starting point. (†) Models that take the 1D-3 model, including the abstraction reactions, as a starting point. (a) $E_d$ is the diffusion energy and $E_b$ the binding energy. (b) $E_A$ is the activation energy of the reaction. (c) $H_{form}(XH)$ is the formation enthalpy of a non-deuterated species and $H_{form}(XD)$ is the formation enthalpy of a deuterated species. (d) From Padovani et al. (2018). (e) From Sipilä et al. (2021).

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 \times 10^5$ yr, the CH$_3$OH abundance reaches its maximum value of approximately $4 \times 10^{-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 \times 10^4$ 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$_3$OH column density profile with an amplitude of $N \approx 1 \times 10^{11}$ cm$^{-2}$, under-producing the observed values for methanol of $3.9 \times 10^{13}$ cm$^{-2}$ at the dust peak and $5.9 \times 10^{13}$ cm$^{-2}$ at the methanol peak by roughly two orders of magnitude. On the other hand, the $N(CH_2DOH)/N(CH_3OH)$ 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$_2$DOH + H$_3$O$^+$ → CH$_3$OHD$^+$ + H$_2$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 $\chi^2$ 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 $\chi^2$ 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 \times 10^5$ yr, which is roughly consistent with estimations from other molecules (e.g. Redaelli et al. 2019a, 2021b). Additionally, it coincides with the occurrence of the highest methanol column density in the probed time frame. For $t = 3.0 \times 10^5$ yr, the fiducial model reaches a column density of methanol of ≈1.5 × 10$^{12}$ cm$^{-2}$, 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, 2008, Wakelam et al. 2010, and references therein).

Figure 5 compares the modelled deuterium fraction $N(CH_2DOH)/N(CH_3OH)$ with the observed one. At the early time step ($t = 1 \times 10^5$ 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
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Gas-phase abundance profiles of H₂CO, CH₃OH, and CH₂DOH for the 1D-4 model for four different time steps ranging from 10⁵ yr (left) to 10⁶ yr (right). The best-fit time is t = 3 × 10⁵ yr.</figcaption></figure>

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.

## 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
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Gas-phase abundance profiles of H₂CO, CH₃OH, and CH₂DOH for the 1D-4 model for four different time steps ranging from 10⁵ yr (left) to 10⁶ yr (right). The best-fit time is t = 3 × 10⁵ yr.</figcaption></figure>

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.

## 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
<figure><img src="image_5.png" /><figcaption>Fig. 4. Column density profiles of H₂CO, CH₃OH, and CH₂DOH for the 1D-4 model for four different time steps ranging from 10⁵ yr to 10⁶ 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.</figcaption></figure>

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<sub>d</sub>/E<sub>b</sub> 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<fnref n="2" /> and KIDA<fnref n="3" />. 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:

$$HCO + H \xrightarrow[\text{pyRate:15.4\%}]{\text{MONACO:5.4\%}} H_2CO \quad (7)$$

$$CH_2OH/CH_3O + H \xrightarrow[\text{pyRate:0.05\%/0.08\%}]{\text{MONACO:0.64\%}} CH_3OH. \quad (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

<footnote n="2">² https://webbook.nist.gov/chemistry</footnote>

<footnote n="3">³ https://kida.astrochem-tools.org</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Column density profiles of H₂CO, CH₃OH, and CH₂DOH for the 1D-4 model for four different time steps ranging from 10⁵ yr to 10⁶ 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.</figcaption></figure>

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<sub>d</sub>/E<sub>b</sub> 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<fnref n="2" /> and KIDA<fnref n="3" />. 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:

$$HCO + H \xrightarrow[\text{pyRate:15.4\%}]{\text{MONACO:5.4\%}} H_2CO \quad (7)$$

$$CH_2OH/CH_3O + H \xrightarrow[\text{pyRate:0.05\%/0.08\%}]{\text{MONACO:0.64\%}} CH_3OH. \quad (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

<footnote n="2">² https://webbook.nist.gov/chemistry</footnote>

<footnote n="3">³ https://kida.astrochem-tools.org</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for the 1D-4 model for four different time steps ranging from 10⁵ yr (top left) to 10⁶ 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).</figcaption></figure>

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 \times 10^5$ 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 \times 10^4$ 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 $\chi^2$ value of the observed column densities versus the corresponding values in the 1D-4 model: $t = 3.0 \times 10^5$ 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 \times 10^5$ 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.

## 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
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for the 1D-4 model for four different time steps ranging from 10⁵ yr (top left) to 10⁶ 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).</figcaption></figure>

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 \times 10^5$ 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 \times 10^4$ 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 $\chi^2$ value of the observed column densities versus the corresponding values in the 1D-4 model: $t = 3.0 \times 10^5$ 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 \times 10^5$ 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.

## 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
<figure><img src="image_7.png" /><figcaption>**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 × 10⁵ yr for V17, *t* = 3.0 × 10⁴ yr and *t* = 3.0 × 10⁵ yr for S16, and *t* = 3.0 × 10⁵ 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 C¹⁷O and the average isotopic ratios of ¹⁶O/¹⁸O = 557 and ¹⁸O/¹⁷O = 3.6 in the local interstellar medium (Wilson et al. 1999).</figcaption></figure>

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.

## 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*<sub>d</sub>/*E*<sub>b</sub>, 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
<figure><img src="image_7.png" /><figcaption>**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 × 10⁵ yr for V17, *t* = 3.0 × 10⁴ yr and *t* = 3.0 × 10⁵ yr for S16, and *t* = 3.0 × 10⁵ 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 C¹⁷O and the average isotopic ratios of ¹⁶O/¹⁸O = 557 and ¹⁸O/¹⁷O = 3.6 in the local interstellar medium (Wilson et al. 1999).</figcaption></figure>

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.

## 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*<sub>d</sub>/*E*<sub>b</sub>, 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
<table><caption>**Table 3.** N(CH₂DOH)/N(CH₃OH) ratio at the centre of the pre-stellar core, r_cen, and the radius of the methanol peak, r_max = 5200 AU.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">t = 1 × 10<sup>5</sup> yr</th><th colspan="2">t = 3 × 10<sup>5</sup> yr</th><th colspan="2">t = 5 × 10<sup>5</sup> yr</th><th colspan="2">t = 1 × 10<sup>6</sup> yr</th></tr><tr><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th></tr></thead><tbody><tr><td>1D-1</td><td>3.8</td><td>3.3</td><td>5.4</td><td>4.6</td><td>6.2</td><td>5.0</td><td>6.6</td><td>5.1</td></tr><tr><td>1D-2</td><td>3.8</td><td>3.3</td><td>5.2</td><td>4.5</td><td>5.9</td><td>4.8</td><td>6.4</td><td>4.9</td></tr><tr><td>1D-3</td><td>3.0</td><td>0.5</td><td>3.3</td><td>1.1</td><td>2.4</td><td>1.6</td><td>4.8</td><td>4.2</td></tr><tr><td>1D-4*</td><td>8.8</td><td>6.2</td><td>8.2</td><td>8.0</td><td>7.1</td><td>6.9</td><td>6.3</td><td>6.1</td></tr><tr><td>1D-4-1</td><td>16.1</td><td>4.5</td><td>8.0</td><td>7.3</td><td>8.1</td><td>7.5</td><td>7.7</td><td>7.2</td></tr><tr><td>1D-4-2</td><td>2.5</td><td>2.3</td><td>3.0</td><td>2.5</td><td>3.0</td><td>2.2</td><td>2.7</td><td>2.0</td></tr><tr><td>1D-5</td><td>1.1</td><td>0.2</td><td>1.3</td><td>0.5</td><td>1.7</td><td>1.3</td><td>2.0</td><td>2.3</td></tr><tr><td>1D-6</td><td>5.8</td><td>4.9</td><td>8.0</td><td>7.7</td><td>11.4</td><td>11.3</td><td>13.6</td><td>13.5</td></tr><tr><td>1D-7</td><td>0.1</td><td>0.04</td><td>0.2</td><td>0.04</td><td>0.2</td><td>0.04</td><td>0.1</td><td>0.03</td></tr><tr><td>1D-8</td><td>4.3</td><td>0.8</td><td>6.8</td><td>2.3</td><td>5.7</td><td>3.0</td><td>3.5</td><td>1.4</td></tr><tr><td>1D-9</td><td>7.0</td><td>1.4</td><td>2.7</td><td>1.1</td><td>3.4</td><td>2.2</td><td>6.3</td><td>5.5</td></tr><tr><td>1D-10</td><td>3.0</td><td>0.5</td><td>3.2</td><td>1.1</td><td>2.3</td><td>1.5</td><td>4.6</td><td>4.0</td></tr><tr><td>1D-11</td><td>3.0</td><td>0.5</td><td>4.3</td><td>1.4</td><td>2.2</td><td>1.4</td><td>4.9</td><td>4.2</td></tr><tr><td>1D-12</td><td>0.8</td><td>0.1</td><td>2.3</td><td>1.1</td><td>3.1</td><td>2.2</td><td>7.4</td><td>6.9</td></tr><tr><td>1D-13</td><td>3.9</td><td>0.7</td><td>2.7</td><td>1.0</td><td>2.3</td><td>1.5</td><td>4.3</td><td>3.7</td></tr></tbody></table>

**Notes.** The position of the methanol peak is determined by the profile of the fiducial model at the best-fit time. An overview of the chemical and physical properties of the various models is given in Table 2. The model marked with a star is the fiducial model and $t = 3 \times 10^5$ yr is the best-fit time obtained as described in Sect. 3.1. The observed values are N(CH₂DOH)/N(CH₃OH)(r_cen) = 7.1% ± 1.9% for the centre and N(CH₂DOH)/N(CH₃OH)(r_max) = 4.5% ± 1.2% for a radius of 5200 AU.

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^7$ cm⁻², which is several orders of magnitude lower than the observed value of $10^{13}$ 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^{12}$ 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 \times 10^5$ yr) and intermediate time steps ($t = 5 \times 10^5$ 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_2DOH)/N(CH_3OH)$ 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 \times 10^5$ 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:

$$H_2CO + H \xrightarrow{E_A=5.16 \times 10^3 \text{ K}} CH_2OH \quad (9)$$

$$H_2CO + H \xrightarrow{E_A=2.00 \times 10^3 \text{ K}} CH_3O. \quad (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:

$$CH_3OH + H/D \xrightarrow{E_A=3.62 \times 10^3/3.24 \times 10^3 \text{ K}} CH_2OH + H_2/HD \quad (11)$$

$$CH_3OH + H/D \xrightarrow{E_A=5.56 \times 10^3/5.16 \times 10^3 \text{ K}} CH_3O + H_2/HD. \quad (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
<table><caption>**Table 3.** N(CH₂DOH)/N(CH₃OH) ratio at the centre of the pre-stellar core, r_cen, and the radius of the methanol peak, r_max = 5200 AU.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">t = 1 × 10<sup>5</sup> yr</th><th colspan="2">t = 3 × 10<sup>5</sup> yr</th><th colspan="2">t = 5 × 10<sup>5</sup> yr</th><th colspan="2">t = 1 × 10<sup>6</sup> yr</th></tr><tr><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th><th>r<sub>cen</sub><br>(%)</th><th>r<sub>max</sub><br>(%)</th></tr></thead><tbody><tr><td>1D-1</td><td>3.8</td><td>3.3</td><td>5.4</td><td>4.6</td><td>6.2</td><td>5.0</td><td>6.6</td><td>5.1</td></tr><tr><td>1D-2</td><td>3.8</td><td>3.3</td><td>5.2</td><td>4.5</td><td>5.9</td><td>4.8</td><td>6.4</td><td>4.9</td></tr><tr><td>1D-3</td><td>3.0</td><td>0.5</td><td>3.3</td><td>1.1</td><td>2.4</td><td>1.6</td><td>4.8</td><td>4.2</td></tr><tr><td>1D-4*</td><td>8.8</td><td>6.2</td><td>8.2</td><td>8.0</td><td>7.1</td><td>6.9</td><td>6.3</td><td>6.1</td></tr><tr><td>1D-4-1</td><td>16.1</td><td>4.5</td><td>8.0</td><td>7.3</td><td>8.1</td><td>7.5</td><td>7.7</td><td>7.2</td></tr><tr><td>1D-4-2</td><td>2.5</td><td>2.3</td><td>3.0</td><td>2.5</td><td>3.0</td><td>2.2</td><td>2.7</td><td>2.0</td></tr><tr><td>1D-5</td><td>1.1</td><td>0.2</td><td>1.3</td><td>0.5</td><td>1.7</td><td>1.3</td><td>2.0</td><td>2.3</td></tr><tr><td>1D-6</td><td>5.8</td><td>4.9</td><td>8.0</td><td>7.7</td><td>11.4</td><td>11.3</td><td>13.6</td><td>13.5</td></tr><tr><td>1D-7</td><td>0.1</td><td>0.04</td><td>0.2</td><td>0.04</td><td>0.2</td><td>0.04</td><td>0.1</td><td>0.03</td></tr><tr><td>1D-8</td><td>4.3</td><td>0.8</td><td>6.8</td><td>2.3</td><td>5.7</td><td>3.0</td><td>3.5</td><td>1.4</td></tr><tr><td>1D-9</td><td>7.0</td><td>1.4</td><td>2.7</td><td>1.1</td><td>3.4</td><td>2.2</td><td>6.3</td><td>5.5</td></tr><tr><td>1D-10</td><td>3.0</td><td>0.5</td><td>3.2</td><td>1.1</td><td>2.3</td><td>1.5</td><td>4.6</td><td>4.0</td></tr><tr><td>1D-11</td><td>3.0</td><td>0.5</td><td>4.3</td><td>1.4</td><td>2.2</td><td>1.4</td><td>4.9</td><td>4.2</td></tr><tr><td>1D-12</td><td>0.8</td><td>0.1</td><td>2.3</td><td>1.1</td><td>3.1</td><td>2.2</td><td>7.4</td><td>6.9</td></tr><tr><td>1D-13</td><td>3.9</td><td>0.7</td><td>2.7</td><td>1.0</td><td>2.3</td><td>1.5</td><td>4.3</td><td>3.7</td></tr></tbody></table>

**Notes.** The position of the methanol peak is determined by the profile of the fiducial model at the best-fit time. An overview of the chemical and physical properties of the various models is given in Table 2. The model marked with a star is the fiducial model and $t = 3 \times 10^5$ yr is the best-fit time obtained as described in Sect. 3.1. The observed values are N(CH₂DOH)/N(CH₃OH)(r_cen) = 7.1% ± 1.9% for the centre and N(CH₂DOH)/N(CH₃OH)(r_max) = 4.5% ± 1.2% for a radius of 5200 AU.

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^7$ cm⁻², which is several orders of magnitude lower than the observed value of $10^{13}$ 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^{12}$ 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 \times 10^5$ yr) and intermediate time steps ($t = 5 \times 10^5$ 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 \times 10^5$ 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:

$$H_2CO + H \xrightarrow{E_A=5.16 \times 10^3 \text{ K}} CH_2OH \quad (9)$$

$$H_2CO + H \xrightarrow{E_A=2.00 \times 10^3 \text{ K}} CH_3O. \quad (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:

$$CH_3OH + H/D \xrightarrow{E_A=3.62 \times 10^3/3.24 \times 10^3 \text{ K}} CH_2OH + H_2/HD \quad (11)$$

$$CH_3OH + H/D \xrightarrow{E_A=5.56 \times 10^3/5.16 \times 10^3 \text{ K}} CH_3O + H_2/HD. \quad (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
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵ yr to 10⁶ yr. The best-fit time is t = 3 × 10⁵ yr. We show two models with slow thermal hopping (*E*<sub>d</sub>/*E*<sub>b</sub> = 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 (*E*<sub>d</sub>/*E*<sub>b</sub> = 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).</figcaption></figure>

0D-models in Taquet et al. (2012). They found large enhance-
ments for the N(CH₂DOH)/N(CH₃OH) ratio in the high density
case (nH ≥ 1 × 10⁶ cm⁻³) eventually reaching values above unity,
while there was no significant increase for the low density case
(nH = 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 abstrac-
tion 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 inter-
mediate 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 for-
mation (≈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 obser-
vations 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 val-
ues 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 increas-
ing far above the observed values, until it reaches values of
almost 0.14 for late times (t = 1 × 10⁶ yr).
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵ yr to 10⁶ yr. The best-fit time is t = 3 × 10⁵ yr. We show two models with slow thermal hopping (*E*<sub>d</sub>/*E*<sub>b</sub> = 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 (*E*<sub>d</sub>/*E*<sub>b</sub> = 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).</figcaption></figure>

0D-models in Taquet et al. (2012). They found large enhance-
ments for the N(CH₂DOH)/N(CH₃OH) ratio in the high density
case (nH ≥ 1 × 10⁶ cm⁻³) eventually reaching values above unity,
while there was no significant increase for the low density case
(nH = 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 abstrac-
tion 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 inter-
mediate 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 for-
mation (≈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 obser-
vations 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 val-
ues 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 increas-
ing far above the observed values, until it reaches values of
almost 0.14 for late times (t = 1 × 10⁶ yr).
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Modelled column density profiles of non-deuterated methanol (CH₃OH) and singly deuterated methanol (CH₂DOH) for several models at four different time steps ranging from 10⁵yr to 10⁶ 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 CH₃OH column densities, and the dashed lines indicate the CH₂DOH column densities.</figcaption></figure>

### 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
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Modelled column density profiles of non-deuterated methanol (CH₃OH) and singly deuterated methanol (CH₂DOH) for several models at four different time steps ranging from 10⁵yr to 10⁶ 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 CH₃OH column densities, and the dashed lines indicate the CH₂DOH column densities.</figcaption></figure>

### 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
<figure><img src="image_10.png" /><figcaption>Fig. 9. Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵ yr to 10⁶ 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).</figcaption></figure>

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

$$H_2CO + H \rightarrow HCO + H_2, \qquad (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.

#### 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
<figure><img src="image_10.png" /><figcaption>Fig. 9. Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵ yr to 10⁶ 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).</figcaption></figure>

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

$$H_2CO + H \rightarrow HCO + H_2, \qquad (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.

#### 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$_3$OH. Reactions (14)–(16),

$$CO + H \xrightarrow{E_A=1.76 \times 10^3 \text{ K}} HCO \quad (14)$$

$$H_2CO + H \xrightarrow{E_A=5.16 \times 10^3 \text{ K}} CH_2OH \quad (15)$$

$$H_2CO + H \xrightarrow{E_A=2.00 \times 10^3 \text{ K}} CH_3O, \quad (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_2DOH)/N(CH_3OH)$ 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_2DOH)/N(CH_3OH)$ 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$_2$DOH column density (see Fig. C.1). For early times ($t = 1 \times 10^5$ yr), the 1D-9 model shows a twice as high CH$_2$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 \times 10^5$ yr. At this time step, the 1D-8 model has a slightly lower CH$_2$DOH column density profile as compared to the 1D-3 model. After the temporal methanol peak, the difference increases for intermediate ($t = 5.0 \times 10^5$ yr) and late ($t = 1 \times 10^6$ yr) time steps to a factor between 1 and 2. The CH$_3$OH column densities differ by a negligible amount between the two models (see Fig. C.1). Consequently, the $N(CH_2DOH)/N(CH_3OH)$ profile (see Fig. C.2) experiences an increase as well. For early time steps the $N(CH_2DOH)/N(CH_3OH)$ ratio reaches high values of up to 16.6 at $t = 1.15 \times 10^5$ 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 \times 10^6$ 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 $\Delta 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<fnref n="4" />. 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$_3$OH and CH$_2$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$_3$OH and CH$_2$DOH column densities and on the $N(CH_2DOH)/N(CH_3OH)$ 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 \times 10^5$ yr.)

## 4.2. Physical variation

### 4.2.1. 1D-11: Variation of the gas temperature

Based on observations of the NH$_3$ (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$_3$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$_2$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_2DOH)/N(CH_3OH)$ 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 \ge 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 $H_2^+$, which then in turn reacts again with the large reservoir of hydrogen molecules, thereby forming the $H_3^+$ ion. This particular ion can react with deuterated molecular hydrogen HD in the following reversible reaction:

$$H_3^+ + HD \rightleftharpoons H_2D^+ + H_2. \quad (17)$$

<footnote n="4">⁴ https://webbook.nist.gov/chemistry/</footnote>
CH$_3$OH. Reactions (14)–(16),

$$CO + H \xrightarrow{E_A=1.76 \times 10^3 \text{ K}} HCO \quad (14)$$

$$H_2CO + H \xrightarrow{E_A=5.16 \times 10^3 \text{ K}} CH_2OH \quad (15)$$

$$H_2CO + H \xrightarrow{E_A=2.00 \times 10^3 \text{ K}} CH_3O, \quad (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_2DOH)/N(CH_3OH)$ 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_2DOH)/N(CH_3OH)$ 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$_2$DOH column density (see Fig. C.1). For early times ($t = 1 \times 10^5$ yr), the 1D-9 model shows a twice as high CH$_2$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 \times 10^5$ yr. At this time step, the 1D-8 model has a slightly lower CH$_2$DOH column density profile as compared to the 1D-3 model. After the temporal methanol peak, the difference increases for intermediate ($t = 5.0 \times 10^5$ yr) and late ($t = 1 \times 10^6$ yr) time steps to a factor between 1 and 2. The CH$_3$OH column densities differ by a negligible amount between the two models (see Fig. C.1). Consequently, the $N(CH_2DOH)/N(CH_3OH)$ profile (see Fig. C.2) experiences an increase as well. For early time steps the $N(CH_2DOH)/N(CH_3OH)$ ratio reaches high values of up to 16.6 at $t = 1.15 \times 10^5$ 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 \times 10^6$ 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 $\Delta 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<fnref n="4" />. 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$_3$OH and CH$_2$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$_3$OH and CH$_2$DOH column densities and on the $N(CH_2DOH)/N(CH_3OH)$ 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 \times 10^5$ yr.)

## 4.2. Physical variation

### 4.2.1. 1D-11: Variation of the gas temperature

Based on observations of the NH$_3$ (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$_3$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$_2$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_2DOH)/N(CH_3OH)$ 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 \ge 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 $H_2^+$, which then in turn reacts again with the large reservoir of hydrogen molecules, thereby forming the $H_3^+$ ion. This particular ion can react with deuterated molecular hydrogen HD in the following reversible reaction:

$$H_3^+ + HD \rightleftharpoons H_2D^+ + H_2. \quad (17)$$

<footnote n="4">⁴ https://webbook.nist.gov/chemistry/</footnote>
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 H₃⁺ 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 \times 10^{-17} s^{-1}$ 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 L-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 \times 10^{-17} s^{-1}$ in the innermost cell to $ζ(H) = 4.89 \times 10^{-17} s^{-1}$ 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 \times 10^5$ 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 \times 10^5$ yr) it is approximately a factor of 2 and even lower at late times ($t = 1 \times 10^6$ 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 isotopologue. CH₂DOH is amplified by a factor of 3 at early times ($t = 1 \times 10^5$ 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 \times 10^6$ 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_{CR}(i) = f(a, T_{max})k_{therm}(i, T_{max}). \quad (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 \times 10^{-5} s/3.16 \times 10^{13} s = 3.16 \times 10^{-19}$ 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 ionisation rate $ζ(H) = 1.3 \times 10^{-17} s^{-1}$, 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 \times 10^5$ 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 desorption 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 \times 10^5$ 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$),
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 H₃⁺ 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 \times 10^{-17} s^{-1}$ 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 L-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 \times 10^{-17} s^{-1}$ in the innermost cell to $ζ(H) = 4.89 \times 10^{-17} s^{-1}$ 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 \times 10^5$ 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 \times 10^5$ yr) it is approximately a factor of 2 and even lower at late times ($t = 1 \times 10^6$ 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 isotopologue. CH₂DOH is amplified by a factor of 3 at early times ($t = 1 \times 10^5$ 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 \times 10^6$ 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_{CR}(i) = f(a, T_{max})k_{therm}(i, T_{max}). \quad (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 \times 10^{-5} s/3.16 \times 10^{13} s = 3.16 \times 10^{-19}$ 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 ionisation rate $ζ(H) = 1.3 \times 10^{-17} s^{-1}$, 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 \times 10^5$ 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 desorption 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 \times 10^5$ 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$_3$OH and CH$_2$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$_2$DOH)/N(CH$_3$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$_2$DOH)/N(CH$_3$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$_2$DOH)/N(CH$_3$OH) ratio:

- A two-phase model resulted in higher CH$_3$OH and CH$_2$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$_3$OH and CH$_2$DOH column densities and thus a slightly better agreement with the observations;

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

Only small effects on the CH$_3$OH and CH$_2$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 deuterated 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$_2$DOH)/N(CH$_3$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.

## References

Aikawa, Y., Wakelam, V., & Hersant, F. 2012, ApJ, 760, 40

Bacmann, A., Tacquet, V., Faure, A., et al. 2012, A&A, 541, A12

Bizzocchi, L., Caselli, P., Spezzano, S., et al. 2014, A&A 569, A27

Caselli, P., Walmsley, C. M., Tafalla, M., et al. 1999, ApJ, 523, L165

Cernicharo, J., Marcelino, N., Roueff, E., et al. 2012, ApJ, 759, L43

Chacón-Tanarro, A., Caselli, P., Bizzocchi, L., et al. 2019, A&A, 622, A141

Chuang, K.-J., Fedoseev, G., Qasium, D., et al. 2018, ApJ, 853, 102

Du, F., Parise, B., & Bergman, P. 2012, A&A, 538, A91

Fredon, A., Radchenko, A. K., & Cuppen, H. M. 2021, Acc. Chem. Res., 54, 745

Fuchs, G. W., Cuppen, H. W., & Ioppolo, S. 2009, A&A, 505, 629

Furuya, K., Hama, T., Oba, Y., et al. 2022, ApJ, 933, L16

Garrod, R. T., Park, I. H., Caselli, P., et al. 2006, Faraday Discuss., 133, 51

Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103

Geppert, W. D., Hellberg, F., Österdahl, F., et al. 2005, Proc. IAU Symp., 231, 117

Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, ApJ, 757, 185

Harju, J., Pineda, J. E., Vasyunin, A. I., et al. 2020, ApJ, 895, 101

Hasegawa, T., & Herbst, E. 1993, MNRAS, 261, 83

Hidaka, H., Kouchi, A., & Watanabe, N. J. 2007, Chem. Phys., 126, 204707

Hidaka, H., Watanabe, M., Kouchi, A., et al. 2009, ApJ, 702, 291

Jiménez-Serra, I., Vasyunin, A. I., Caselli, P., et al. 2016, ApJ, 830, L6

Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625

Kimura, Y., Tsuge, M., Pirronello, V., et al. 2018, ApJ, 858, L23

Léger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147

Minissale, M., Congiu, E., & Dulieu, F. 2016a, A&A, 585, A146

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016b, A&A, 585, A24

Minissale, M., Moudens, A., Baouche, S., et al. 2016c, MNRAS, 458, 2953

Padovani, M., Ilev, A. V., & Galli, D. 2018, A&A, 614, A111

Pagani, L., Salez, M., & Wannier, P. G. 1992, A&A, 258, 479

Punanova, A., Vasyunin, A. I., Caselli, P., et al. 2022, ApJ, 927, 213

Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15

Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109

Scibelli, S., & Shirley, Y. 2020, ApJ, 891, 73

Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42

Sipilä, O., Caselli, P., & Harju, J. 2015, A&A, 578, A55

Sipilä, O., Caselli, P., & Taquet, V. 2016, A&A, 591, A9

Sipilä, O., Caselli, P., & Redaelli, E., et al. 2019a, MNRAS, 487, 1269

Sipilä, O., Caselli, P., & Harju, J. 2019b, A&A, 631, A63

Sipilä, O., Silsbee, K., & Caselli, P. 2021, ApJ, 922, 126

Spezzano, S., Bizzocchi, L., Caselli, P., et al. 2016, A&A, 592, L11

Spezzano, S., Caselli, P., Pineda, J. E., et al. 2020, A&A, 643, A60

Taquet, V., Ceccarelli, C., & Kahane, C. 2012, ApJ, 748, L3

Vastel, C., Ceccarelli, C., Lefloch, B., et al. 2014, ApJ, 795, L2

Vasyunin, A. I., Sobolev, A. M., Wiebe, D. S., et al. 2004, Astron. Lett., 30, 566

Vasyunin, A. I., Semenov, D., Henning, Th., et al. 2008, ApJ, 672, 629

Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33

Wakelam, V., Herbst, E., Le Bourlot, J., et al. 2010, A&A, 517, A21

Wakelam, V., Loison, J. C., Herbst, E., et al. 2015, ApJ, 217, 20

Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173

Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143
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$_3$OH and CH$_2$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$_2$DOH)/N(CH$_3$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$_2$DOH)/N(CH$_3$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$_2$DOH)/N(CH$_3$OH) ratio:

- A two-phase model resulted in higher CH$_3$OH and CH$_2$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$_3$OH and CH$_2$DOH column densities and thus a slightly better agreement with the observations;

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

Only small effects on the CH$_3$OH and CH$_2$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 deuterated 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$_2$DOH)/N(CH$_3$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.

## References

Aikawa, Y., Wakelam, V., & Hersant, F. 2012, ApJ, 760, 40

Bacmann, A., Tacquet, V., Faure, A., et al. 2012, A&A, 541, A12

Bizzocchi, L., Caselli, P., Spezzano, S., et al. 2014, A&A 569, A27

Caselli, P., Walmsley, C. M., Tafalla, M., et al. 1999, ApJ, 523, L165

Cernicharo, J., Marcelino, N., Roueff, E., et al. 2012, ApJ, 759, L43

Chacón-Tanarro, A., Caselli, P., Bizzocchi, L., et al. 2019, A&A, 622, A141

Chuang, K.-J., Fedoseev, G., Qasium, D., et al. 2018, ApJ, 853, 102

Du, F., Parise, B., & Bergman, P. 2012, A&A, 538, A91

Fredon, A., Radchenko, A. K., & Cuppen, H. M. 2021, Acc. Chem. Res., 54, 745

Fuchs, G. W., Cuppen, H. W., & Ioppolo, S. 2009, A&A, 505, 629

Furuya, K., Hama, T., Oba, Y., et al. 2022, ApJ, 933, L16

Garrod, R. T., Park, I. H., Caselli, P., et al. 2006, Faraday Discuss., 133, 51

Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103

Geppert, W. D., Hellberg, F., Österdahl, F., et al. 2005, Proc. IAU Symp., 231, 117

Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, ApJ, 757, 185

Harju, J., Pineda, J. E., Vasyunin, A. I., et al. 2020, ApJ, 895, 101

Hasegawa, T., & Herbst, E. 1993, MNRAS, 261, 83

Hidaka, H., Kouchi, A., & Watanabe, N. J. 2007, Chem. Phys., 126, 204707

Hidaka, H., Watanabe, M., Kouchi, A., et al. 2009, ApJ, 702, 291

Jiménez-Serra, I., Vasyunin, A. I., Caselli, P., et al. 2016, ApJ, 830, L6

Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625

Kimura, Y., Tsuge, M., Pirronello, V., et al. 2018, ApJ, 858, L23

Léger, A., Jura, M., & Omont, A. 1985, A&A, 144, 147

Minissale, M., Congiu, E., & Dulieu, F. 2016a, A&A, 585, A146

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016b, A&A, 585, A24

Minissale, M., Moudens, A., Baouche, S., et al. 2016c, MNRAS, 458, 2953

Padovani, M., Ivlev, A. V., & Galli, D. 2018, A&A, 614, A111

Pagani, L., Salez, M., & Wannier, P. G. 1992, A&A, 258, 479

Punanova, A., Vasyunin, A. I., Caselli, P., et al. 2022, ApJ, 927, 213

Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15

Redaelli, E., Sipilä, O., Padovani, M., et al. 2021, A&A, 656, A109

Scibelli, S., & Shirley, Y. 2020, ApJ, 891, 73

Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42

Sipilä, O., Caselli, P., & Harju, J. 2015, A&A, 578, A55

Sipilä, O., Caselli, P., & Taquet, V. 2016, A&A, 591, A9

Sipilä, O., Caselli, P., & Redaelli, E., et al. 2019a, MNRAS, 487, 1269

Sipilä, O., Caselli, P., & Harju, J. 2019b, A&A, 631, A63

Sipilä, O., Silsbee, K., & Caselli, P. 2021, ApJ, 922, 126

Spezzano, S., Bizzocchi, L., Caselli, P., et al. 2016, A&A, 592, L11

Spezzano, S., Caselli, P., Pineda, J. E., et al. 2020, A&A, 643, A60

Taquet, V., Ceccarelli, C., & Kahane, C. 2012, ApJ, 748, L3

Vastel, C., Ceccarelli, C., Lefloch, B., et al. 2014, ApJ, 795, L2

Vasyunin, A. I., Sobolev, A. M., Wiebe, D. S., et al. 2004, Astron. Lett., 30, 566

Vasyunin, A. I., Semenov, D., Henning, Th., et al. 2008, ApJ, 672, 629

Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33

Wakelam, V., Herbst, E., Le Bourlot, J., et al. 2010, A&A, 517, A21

Wakelam, V., Loison, J. C., Herbst, E., et al. 2015, ApJ, 217, 20

Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173

Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143
# Appendix A: Formation enthalpies and binding energies

In Table A.1 we present the formation enthalpies H<sub>form</sub> and binding energies E<sub>b</sub> for all species involved in exothermic surface reactions, to which we applied the newly developed reactive desorption mechanism laid out in Sect. 2.1.

<table><caption>**Table A.1.** Formation enthalpies, H<sub>form</sub>, and binding energies, E<sub>b</sub>.</caption><thead><tr><th>Species</th><th>H<sub>form</sub><br>[kJ mol<sup>-1</sup>]</th><th>E<sub>b</sub><br>[K]</th><th>Species</th><th>H<sub>form</sub><br>[kJ mol<sup>-1</sup>]</th><th>E<sub>b</sub><br>[K]</th></tr></thead><tbody><tr><td>C</td><td>716.70</td><td>800.0</td><td>O<sub>2</sub></td><td>0.00</td><td>1000.0</td></tr><tr><td>CH</td><td>594.10</td><td>925.0</td><td>O<sub>2</sub>H</td><td>2.10</td><td>3650.0</td></tr><tr><td>CH<sub>2</sub></td><td>386.40</td><td>1050.0</td><td>O<sub>3</sub></td><td>142.70</td><td>1800.0</td></tr><tr><td>CH<sub>3</sub></td><td>145.70</td><td>1175.0</td><td>OCN</td><td>159.40</td><td>2400.0</td></tr><tr><td>CH<sub>3</sub>O</td><td>17.00</td><td>3800.0</td><td>OCS</td><td>138.40</td><td>2888.0</td></tr><tr><td>CH<sub>2</sub>OH</td><td>-9.00</td><td>5084.0</td><td>OH</td><td>39.00</td><td>2850.0</td></tr><tr><td>CH<sub>3</sub>OH</td><td>-201.20</td><td>5534.0</td><td>OD</td><td>36.60<sup>a</sup></td><td>2850.0</td></tr><tr><td>CH<sub>4</sub></td><td>-74.90</td><td>1300.0</td><td>S</td><td>277.00</td><td>1100.0</td></tr><tr><td>CN</td><td>435.10</td><td>1600.0</td><td>SO</td><td>5.00</td><td>2600.0</td></tr><tr><td>CO</td><td>-110.50</td><td>1150.0</td><td>SO<sub>2</sub></td><td>-296.80</td><td>3405.0</td></tr><tr><td>CO<sub>2</sub></td><td>-393.50</td><td>2575.0</td><td>C<sub>2</sub></td><td>837.74</td><td>1600.0</td></tr><tr><td>CS</td><td>280.30</td><td>1900.0</td><td>CCH</td><td>476.98</td><td>2137.0</td></tr><tr><td>H</td><td>218.00</td><td>450.0</td><td>C<sub>2</sub>H<sub>2</sub></td><td>226.73</td><td>2587.0</td></tr><tr><td>D*</td><td>221.72<sup>a</sup></td><td>450.0</td><td>C<sub>3</sub></td><td>820.06<sup>a</sup></td><td>2400.0</td></tr><tr><td>H<sub>2</sub></td><td>0.0</td><td>500.0</td><td>C<sub>4</sub></td><td>970.69<sup>a</sup></td><td>3200.0</td></tr><tr><td>HD*</td><td>0.32<sup>a</sup></td><td>500.0</td><td>C<sub>5</sub></td><td>979.06<sup>a</sup></td><td>4000.0</td></tr><tr><td>H<sub>2</sub>CO</td><td>-115.90</td><td>2050.0</td><td>C<sub>2</sub>H<sub>3</sub></td><td>299.00<sup>a</sup></td><td>3037.0</td></tr><tr><td>H<sub>2</sub>O</td><td>-241.80</td><td>5700.0</td><td>C<sub>4</sub>H<sub>2</sub></td><td>464.00<sup>a</sup></td><td>4187.0</td></tr><tr><td>HDO*</td><td>-245.37<sup>a</sup></td><td>5700.0</td><td>C<sub>2</sub>H<sub>4</sub></td><td>52.40<sup>a</sup></td><td>3487.0</td></tr><tr><td>D<sub>2</sub>O*</td><td>-249.20<sup>a</sup></td><td>5700.0</td><td>C<sub>2</sub>H<sub>5</sub></td><td>119.00<sup>a</sup></td><td>3937.0</td></tr><tr><td>H<sub>2</sub>O<sub>2</sub></td><td>-135.80</td><td>5700.0</td><td>HC<sub>3</sub>N</td><td>354.00<sup>a</sup></td><td>4580.0</td></tr><tr><td>H<sub>2</sub>S</td><td>-20.5</td><td>2743.0</td><td>H<sub>2</sub>CS</td><td>118.00<sup>a</sup></td><td>2700.0</td></tr><tr><td>D<sub>2</sub>S*</td><td>-23.89<sup>a</sup></td><td>2743.0</td><td>MgH</td><td>169.03<sup>a</sup></td><td>5750.0</td></tr><tr><td>HCN</td><td>135.10</td><td>2050.0</td><td>NaH</td><td>124.27<sup>a</sup></td><td>12250.0</td></tr><tr><td>HNC</td><td>135.10</td><td>2050.0</td><td>PH</td><td>253.55<sup>a</sup></td><td>5000.0</td></tr><tr><td>HCO</td><td>43.50</td><td>1600.0</td><td>PH<sub>2</sub></td><td>125.94<sup>a</sup></td><td>5000.0</td></tr><tr><td>HCS</td><td>296.20</td><td>2350.0</td><td>PH<sub>2</sub></td><td>125.94<sup>a</sup></td><td>5000.0</td></tr><tr><td>HNO</td><td>99.60</td><td>2050.0</td><td>SiH</td><td>376.66<sup>a</sup></td><td>3150.0</td></tr><tr><td>HS</td><td>139.30</td><td>1450.0</td><td>SiH<sub>4</sub></td><td>34.31<sup>a</sup></td><td>4500.0</td></tr><tr><td>DS*</td><td>138.49<sup>a</sup></td><td>1450.0</td><td>NS</td><td>263.59<sup>a</sup></td><td>1900.0</td></tr><tr><td>N</td><td>472.70</td><td>800.0</td><td>CCO</td><td>286.60<sup>a</sup></td><td>1950.0</td></tr><tr><td>N<sub>2</sub></td><td>0.00</td><td>1000.0</td><td>C<sub>4</sub>H</td><td>775.02<sup>b</sup></td><td>3737.0</td></tr><tr><td>N<sub>2</sub>H</td><td>245.20</td><td>1450.0</td><td>C<sub>6</sub></td><td>1261.02<sup>b</sup></td><td>4800.0</td></tr><tr><td>N<sub>2</sub>H<sub>2</sub></td><td>213.00</td><td>4756.0</td><td>C<sub>6</sub>H</td><td>991.80<sup>b</sup></td><td>5337.0</td></tr><tr><td>N<sub>2</sub>D<sub>2</sub></td><td>207.11<sup>a</sup></td><td>4756.0</td><td>C<sub>7</sub></td><td>1309.34<sup>b</sup></td><td>5600.0</td></tr><tr><td>NH</td><td>376.60</td><td>2378.0</td><td>I-C<sub>3</sub>H</td><td>714.09<sup>b</sup></td><td>2937.0</td></tr><tr><td>ND*</td><td>375.31<sup>a</sup></td><td>2378.0</td><td>c-C<sub>3</sub>H</td><td>714.09<sup>b</sup></td><td>2937.0</td></tr><tr><td>NH<sub>2</sub></td><td>190.40</td><td>3956.0</td><td>C<sub>5</sub>H<sub>2</sub></td><td>690.36<sup>b</sup></td><td>4987.0</td></tr><tr><td>ND*</td><td>185.35<sup>a</sup></td><td>3956.0</td><td>H<sub>2</sub>CN</td><td>242.23<sup>b</sup></td><td>2400.0</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>-186.00</td><td>5556.0</td><td>I-C<sub>3</sub>H<sub>2</sub></td><td>650.36<sup>b</sup></td><td>3387.0</td></tr><tr><td>NH<sub>3</sub></td><td>-45.90</td><td>5534.0</td><td>c-C<sub>3</sub>H<sub>2</sub></td><td>477.96<sup>b</sup></td><td>3387.0</td></tr><tr><td>ND*</td><td>-58.58<sup>a</sup></td><td>5534.0</td><td>C<sub>4</sub>H<sub>3</sub></td><td>545.65<sup>b</sup></td><td>4637.0</td></tr><tr><td>NO</td><td>90.30</td><td>1600.0</td><td>SiH<sub>2</sub></td><td>275.00<sup>b</sup></td><td>3600.0</td></tr><tr><td>NO<sub>2</sub></td><td>33.10</td><td>2400.0</td><td>SiH<sub>3</sub></td><td>204.09<sup>b</sup></td><td>4050.0</td></tr><tr><td>O</td><td>249.20</td><td>1390.0</td><td>HOOH</td><td>-129.89<sup>b</sup></td><td>5700.0</td></tr></tbody></table>

**Notes.** Species marked with a star (*) are newly added formation enthalpies for deuterated isotopologues. If not stated otherwise the formation enthalpies are adopted from Du et al. (2012). The values marked with a<sup>a</sup> are adopted from the NIST Chemistry WebBook<sup>(1)</sup> and the ones marked with b<sup>b</sup> are from the Kinetic Database for Astrochemistry<sup>(2)</sup>. The binding energies are taken from Semenov et al. (2010).<sup>(1)</sup>https://webbook.nist.gov/chemistry<sup>(2)</sup>https://kida.astrochem-tools.org
## 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.

<table><caption>**Table A.1.** Formation enthalpies, $H_{form}$, and binding energies, $E_b$.</caption><thead><tr><th>Species</th><th>H<sub>form</sub><br>[kJ mol<sup>-1</sup>]</th><th>E<sub>b</sub><br>[K]</th><th>Species</th><th>H<sub>form</sub><br>[kJ mol<sup>-1</sup>]</th><th>E<sub>b</sub><br>[K]</th></tr></thead><tbody><tr><td>C</td><td>716.70</td><td>800.0</td><td>O<sub>2</sub></td><td>0.00</td><td>1000.0</td></tr><tr><td>CH</td><td>594.10</td><td>925.0</td><td>O<sub>2</sub>H</td><td>2.10</td><td>3650.0</td></tr><tr><td>CH<sub>2</sub></td><td>386.40</td><td>1050.0</td><td>O<sub>3</sub></td><td>142.70</td><td>1800.0</td></tr><tr><td>CH<sub>3</sub></td><td>145.70</td><td>1175.0</td><td>OCN</td><td>159.40</td><td>2400.0</td></tr><tr><td>CH<sub>3</sub>O</td><td>17.00</td><td>3800.0</td><td>OCS</td><td>138.40</td><td>2888.0</td></tr><tr><td>CH<sub>2</sub>OH</td><td>-9.00</td><td>5084.0</td><td>OH</td><td>39.00</td><td>2850.0</td></tr><tr><td>CH<sub>3</sub>OH</td><td>-201.20</td><td>5534.0</td><td>OD</td><td>36.60<sup>a</sup></td><td>2850.0</td></tr><tr><td>CH<sub>4</sub></td><td>-74.90</td><td>1300.0</td><td>S</td><td>277.00</td><td>1100.0</td></tr><tr><td>CN</td><td>435.10</td><td>1600.0</td><td>SO</td><td>5.00</td><td>2600.0</td></tr><tr><td>CO</td><td>-110.50</td><td>1150.0</td><td>SO<sub>2</sub></td><td>-296.80</td><td>3405.0</td></tr><tr><td>CO<sub>2</sub></td><td>-393.50</td><td>2575.0</td><td>C<sub>2</sub></td><td>837.74</td><td>1600.0</td></tr><tr><td>CS</td><td>280.30</td><td>1900.0</td><td>CCH</td><td>476.98</td><td>2137.0</td></tr><tr><td>H</td><td>218.00</td><td>450.0</td><td>C<sub>2</sub>H<sub>2</sub></td><td>226.73</td><td>2587.0</td></tr><tr><td>D*</td><td>221.72<sup>a</sup></td><td>450.0</td><td>C<sub>3</sub></td><td>820.06<sup>a</sup></td><td>2400.0</td></tr><tr><td>H<sub>2</sub></td><td>0.0</td><td>500.0</td><td>C<sub>4</sub></td><td>970.69<sup>a</sup></td><td>3200.0</td></tr><tr><td>HD*</td><td>0.32<sup>a</sup></td><td>500.0</td><td>C<sub>5</sub></td><td>979.06<sup>a</sup></td><td>4000.0</td></tr><tr><td>H<sub>2</sub>CO</td><td>-115.90</td><td>2050.0</td><td>C<sub>2</sub>H<sub>3</sub></td><td>299.00<sup>a</sup></td><td>3037.0</td></tr><tr><td>H<sub>2</sub>O</td><td>-241.80</td><td>5700.0</td><td>C<sub>4</sub>H<sub>2</sub></td><td>464.00<sup>a</sup></td><td>4187.0</td></tr><tr><td>HDO*</td><td>-245.37<sup>a</sup></td><td>5700.0</td><td>C<sub>2</sub>H<sub>4</sub></td><td>52.40<sup>a</sup></td><td>3487.0</td></tr><tr><td>D<sub>2</sub>O*</td><td>-249.20<sup>a</sup></td><td>5700.0</td><td>C<sub>2</sub>H<sub>5</sub></td><td>119.00<sup>a</sup></td><td>3937.0</td></tr><tr><td>H<sub>2</sub>O<sub>2</sub></td><td>-135.80</td><td>5700.0</td><td>HC<sub>3</sub>N</td><td>354.00<sup>a</sup></td><td>4580.0</td></tr><tr><td>H<sub>2</sub>S</td><td>-20.5</td><td>2743.0</td><td>H<sub>2</sub>CS</td><td>118.00<sup>a</sup></td><td>2700.0</td></tr><tr><td>D<sub>2</sub>S*</td><td>-23.89<sup>a</sup></td><td>2743.0</td><td>MgH</td><td>169.03<sup>a</sup></td><td>5750.0</td></tr><tr><td>HCN</td><td>135.10</td><td>2050.0</td><td>NaH</td><td>124.27<sup>a</sup></td><td>12250.0</td></tr><tr><td>HNC</td><td>135.10</td><td>2050.0</td><td>PH</td><td>253.55<sup>a</sup></td><td>5000.0</td></tr><tr><td>HCO</td><td>43.50</td><td>1600.0</td><td>PH<sub>2</sub></td><td>125.94<sup>a</sup></td><td>5000.0</td></tr><tr><td>HCS</td><td>296.20</td><td>2350.0</td><td>PH<sub>2</sub></td><td>125.94<sup>a</sup></td><td>5000.0</td></tr><tr><td>HNO</td><td>99.60</td><td>2050.0</td><td>SiH</td><td>376.66<sup>a</sup></td><td>3150.0</td></tr><tr><td>HS</td><td>139.30</td><td>1450.0</td><td>SiH<sub>4</sub></td><td>34.31<sup>a</sup></td><td>4500.0</td></tr><tr><td>DS*</td><td>138.49<sup>a</sup></td><td>1450.0</td><td>NS</td><td>263.59<sup>a</sup></td><td>1900.0</td></tr><tr><td>N</td><td>472.70</td><td>800.0</td><td>CCO</td><td>286.60<sup>a</sup></td><td>1950.0</td></tr><tr><td>N<sub>2</sub></td><td>0.00</td><td>1000.0</td><td>C<sub>4</sub>H</td><td>775.02<sup>b</sup></td><td>3737.0</td></tr><tr><td>N<sub>2</sub>H</td><td>245.20</td><td>1450.0</td><td>C<sub>6</sub></td><td>1261.02<sup>b</sup></td><td>4800.0</td></tr><tr><td>N<sub>2</sub>H<sub>2</sub></td><td>213.00</td><td>4756.0</td><td>C<sub>6</sub>H</td><td>991.80<sup>b</sup></td><td>5337.0</td></tr><tr><td>N<sub>2</sub>D<sub>2</sub></td><td>207.11<sup>a</sup></td><td>4756.0</td><td>C<sub>7</sub></td><td>1309.34<sup>b</sup></td><td>5600.0</td></tr><tr><td>NH</td><td>376.60</td><td>2378.0</td><td>I-C<sub>3</sub>H</td><td>714.09<sup>b</sup></td><td>2937.0</td></tr><tr><td>ND*</td><td>375.31<sup>a</sup></td><td>2378.0</td><td>c-C<sub>3</sub>H</td><td>714.09<sup>b</sup></td><td>2937.0</td></tr><tr><td>NH<sub>2</sub></td><td>190.40</td><td>3956.0</td><td>C<sub>5</sub>H<sub>2</sub></td><td>690.36<sup>b</sup></td><td>4987.0</td></tr><tr><td>ND*</td><td>185.35<sup>a</sup></td><td>3956.0</td><td>H<sub>2</sub>CN</td><td>242.23<sup>b</sup></td><td>2400.0</td></tr><tr><td>NH<sub>2</sub>CHO</td><td>-186.00</td><td>5556.0</td><td>I-C<sub>3</sub>H<sub>2</sub></td><td>650.36<sup>b</sup></td><td>3387.0</td></tr><tr><td>NH<sub>3</sub></td><td>-45.90</td><td>5534.0</td><td>c-C<sub>3</sub>H<sub>2</sub></td><td>477.96<sup>b</sup></td><td>3387.0</td></tr><tr><td>ND*</td><td>-58.58<sup>a</sup></td><td>5534.0</td><td>C<sub>4</sub>H<sub>3</sub></td><td>545.65<sup>b</sup></td><td>4637.0</td></tr><tr><td>NO</td><td>90.30</td><td>1600.0</td><td>SiH<sub>2</sub></td><td>275.00<sup>b</sup></td><td>3600.0</td></tr><tr><td>NO<sub>2</sub></td><td>33.10</td><td>2400.0</td><td>SiH<sub>3</sub></td><td>204.09<sup>b</sup></td><td>4050.0</td></tr><tr><td>O</td><td>249.20</td><td>1390.0</td><td>HOOH</td><td>-129.89<sup>b</sup></td><td>5700.0</td></tr></tbody></table>

**Notes.** Species marked with a star (*) are newly added formation enthalpies for deuterated isotopologues. If not stated otherwise the formation enthalpies are adopted from Du et al. (2012). The values marked with a are adopted from the NIST Chemistry WebBook (1) and the ones marked with b are from the Kinetic Database for Astrochemistry (2). The binding energies are taken from Semenov et al. (2010). (1) https://webbook.nist.gov/chemistry (2) https://kida.astrochem-tools.org
**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.

<figure><img src="image_11.png" /><figcaption>**Fig. B.1.** Reaction scheme for the formation of CH<sub>3</sub>OH 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 H<sub>2</sub>, HD, or D<sub>2</sub> molecule. The values on top of the arrows indicate, if existing, the activation energies, E<sub>A</sub>, in 10<sup>3</sup> 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 H<sub>2</sub>CO and HDCO, which are not shown here for the sake of clarity.</figcaption></figure>
# 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.

<figure><img src="image_11.png" /><figcaption>**Fig. B.1.** Reaction scheme for the formation of CH₃OH 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 H₂, HD, or D₂ molecule. The values on top of the arrows indicate, if existing, the activation energies, Eₐ, in 10³ 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 H₂CO and HDCO, which are not shown here for the sake of clarity.</figcaption></figure>
**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).

<figure><img src="image_12.png" /><figcaption>**Fig. C.1.** Modelled column density profiles of non-deuterated methanol (CH₃OH) and singly deuterated methanol (CH₂DOH) for several models at four different time steps ranging from 10⁵yr to 10⁶yr. The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH₃OH column densities, and the dashed lines indicate the CH₂DOH column densities.</figcaption></figure>
**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).

<figure><img src="image_12.png" /><figcaption>**Fig. C.1.** Modelled column density profiles of non-deuterated methanol (CH₃OH) and singly deuterated methanol (CH₂DOH) for several models at four different time steps ranging from 10⁵yr to 10⁶yr. The black lines show the observed profiles (errors as grey-shaded areas). The solid lines indicate the CH₃OH column densities, and the dashed lines indicate the CH₂DOH column densities.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. C.2.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵yr to 10⁶yr. The black line shows the observed ratio (errors as grey-shaded areas).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. C.2.** Modelled ratio between singly deuterated methanol (CH₂DOH) and non-deuterated methanol (CH₃OH) for several models at four different time steps ranging from 10⁵yr to 10⁶yr. The black line shows the observed ratio (errors as grey-shaded areas).</figcaption></figure>