Modelling Deuterated Isotopologues of Methanol toward the Pre-Stellar Core L1544

Context. In the extremely cold and dark environments of pre-stellar cores, methanol is formed on the surface of interstellar dust grains and released to the gas phase by 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 a previous model for the prediction of column densities and deuterium fractions of non-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 e ffi ciency for every product species in a chemical reaction, that depends on the reaction enthalpy and type of underlying surface. Furthermore, we explored several options to promote the di ff usion of hydrogen and deuterium atoms over the surface of interstellar dust grains, in order to increase methanol formation. Results. Our fiducial model employs di ff usion by quantum tunneling of hydrogen and deuterium atoms, resulting in CH 3 OH and CH 2 DOH column densities that are approximately an order of magnitude lower than the observed values, which improves the results compared to the previous model by a factor 10. The N (CH 2 DOH) / N (CH 3 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 di ff usion process with a high di ff usion 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 di ff usion rate.


Introduction
Methanol (CH 3 OH) is the simplest O-bearing complex organic molecule (COM) in the interstellar medium and an important precursor for saturated, more massive, COMs that are formed in the gas phase.Previous works, 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 investigated in an extensive manner experimentally 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 tunneling reactions (Hidaka et al. 2007) and is therefore progressing with a much higher rate than in the gas phase.Additionally, the existence of so-called abstraction reactions, that remove H 2 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 between 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 by 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, e.g.: Bacmann et al. (2012), Cernicharo et al. (2012) and Jiménez-Serra et al. (2016), conducted toward dark molecular clouds, found comparatively high abundances of methanol and other complex organic molecules.
In order to explain these unexpected findings, astrochemical models need to employ other mechanisms for the evaporation of surface molecules.One promising candidate, reactive desorption, is based on the energy released in an exothermic chemical reaction, which can lead to the desorption of the reaction product(s).Nowadays, most chemical codes include a simple treat-ment of reactive desorption following the recipe described by Garrod et al. (2007).According to those authors, a (nearly) constant reactive desorption efficiency, of typically 1%, is employed independent of the desorbing molecule.However, recent laboratory experiments found a strong dependence of the reactive desorption efficiency on the chemical reaction and the type of underlying surface (Minissale et al. 2016b;Chuang et al. 2018).
In this paper, we develop 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).The motivation for our work are the column density maps of CH 3 OH and CH 2 DOH and theoretical predictions presented by Chacón-Tanarro et al. (2019).They carried out single-dish observations with the 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 the theoretical column density profiles (hereafter model S16) made by the chemical code pyRate (Sipilä et al. 2015a(Sipilä et al. , 2019b)), particularly to get more accurate predictions about the deuteration of methanol.Additionally, we compare our predictions with the results of the chemical code presented in Vasyunin et al. (2017) (hereafter V17 model).These two models were used in the work of Chacón-Tanarro et al. (2019), where it was 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 by using the N(CH 2 DOH)/N(CH 3 OH) ratio predicted by pyRate, as 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 implement here 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 found that diffusion by quantum tunneling of atomic and molecular hydrogen had to be employed in order to explain the observed methanol column density profile.In the present paper, we employ tunneling diffusion for hydrogen and deuterium atoms for our fiducial model, neglecting the tunneling diffusion of molecular hydrogen.However, we also explore the option to promote thermal diffusion by decreasing the diffusionto-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 and how well they are able to explain the observed deuterium fraction profile.
The paper is structured as follows: Section 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 it with the observationally obtained column density and deuterium fraction profiles by Chacón-Tanarro et al. (2019) and with the model presented in 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 used chemical parameters and reaction schemes.

Treatment of Reactive Desorption
For the extension of the reactive desorption mechanism, we adopt 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: where ∆H is the reaction enthalpy, E b the binding energy, N the number of the degrees of freedom (translational, rotational and vibrational) of the reaction product and ϵ the fraction of kinetic energy retained by the reaction product.We assume 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 by N = 3n atoms with n atoms being the number of atoms of the reaction product.This approach has the effect that the less complex molecules, consisting of fewer atoms, can use a larger share of the reaction enthalpy for vertical motion off the surface.For the more complex molecules, the total number of degrees of freedom increases quickly with increasing number of atoms, due to the fact that larger molecules have more possibilities for vibrational excitation.However, the number of translational degrees of freedom stays the same.As a consequence, a smaller share of the reaction enthalpy is converted into vertical motion.For example, a forming CO molecule, with only 2 atoms, receives 1/6 of the reaction enthalpy for vertical motion, while a forming CH 3 OH molecule, with 6 atoms, receives only 1/18.The dependency of the reactive desorption efficiency on the type of surface is expressed by the fraction of kinetic energy ϵ received by the reaction product with mass m when colliding with a surface element with an effective mass M: 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 have 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 have taken a somewhat lower value of M = 120 amu that has been suggested as a common value for carbonaceous and silicate grains following the experiments conducted by Minissale et al. (2016b).Every surface which does not consist of H 2 O or bare grain is treated as if it were CO.Considering that grain surfaces in the inner, very cold regions of pre-stellar cores are covered, aside from CO itself, by species with a similar molecular weight (e.g.: N 2 , H 2 CO and CH 3 OH), this seems to be a reasonable approximation given the lack of more detailed measurements.To consider that a grain surface might be covered by a mixture of bare surface, H 2 O or CO, we follow 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 scale 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: where RD j is the individual reactive desorption efficiency for either bare grain, H 2 O 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 oneproduct 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 assume that the desorbing molecules are first isolated from the surface and then undergo an elastic collision with the surface of the grain to gain velocity in the vertical direction.Considering the conservation of linear momentum, one can derive an expression for the kinetic energy E kin received by a product species i: 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 derive 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.Note that the number of translational degrees of freedom is the combined number for both reactants: N trans = 3n prod , where n prod is the number of products.The total number of degrees of freedom is defined similarly as the sum of the individual degrees of freedom for both reaction partners: N tot = i 3n atoms,i , with n atoms,i being the number of atoms of reaction product i.Here, we consider 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 can derive an equation for the reactive desorption efficiency of a two-product reaction (5), that simplifies to Eq. ( 1) in the case that the occurring reaction has only one reaction product: 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:

Chemical Model
We incorporate the new description of reactive desorption into the rate-equation based chemical code pyRate, described in more detail in Sipilä et al. (2015aSipilä et al. ( , 2019b)).It tracks the chemical evolution both in the gas phase and on the grain surface.The basis of the chemical network is the 2014 public release of the KIDA gas phase network (kida.uva.2014(kida.uva. , Wakelam et al. 2015)), which was extended by deuterium chemistry for molecules with up to seven atoms.The code also tracks the various spin states of the light hydrogen-bearing species H 2 , H + 2 and H + 3 and their deuterated isotopologues, as well as multiply protonated or deuterated species involved in the water and ammonia formation networks.All together, the network includes ≈ 74000 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 Figure 1 (representation with H) and Figure B.1 (more extensive representation with H and D).We adopt 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 employ a three-phase model, consisting of a gas phase, a chemically-active surface phase and a chemically-inert mantle phase.The dust grains are assumed to be spherically symmetric with a radius of 0.1 µm.
For exploring our new treatment of reactive desorption, the choice of binding energies E b and formation enthalpies H form is crucial.The values of E b and H form are displayed in Table (A.1).The binding energies are taken from Semenov et al. (2010).Notes. (a) 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 in the appendix.Unfortunately, experimental values for deuterated molecules are quite scarce.For this reason, we apply 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 1 .

Physical Model
As mentioned in Section 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 3 OH and CH 2 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 2density n(H 2 ), the gas temperature T gas , the dust temperature T dust and the visual extinction A V as shown in Figure 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 calculate column densities by integrating along the line-of-sight for different impact factors from the core centre.Afterwards, the column density distribution is 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 have run 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 Section 2.1.An overview of the various models is presented in Table 2.

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 tunneling is enabled, while, abstraction reactions as shown in Figure 1 (or their deuterated analogues) are not included.This choice is discussed in more detail in Section 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 5 yr to 10 6 yr in the fiducial model.Figure 3 shows the molecular abundances with respect to H 2 and Figure 4 shows the column density profiles for H 2 CO, CH 3 OH and CH 2 DOH for four different time steps (1 × 10 5 yr, 3 × 10 5 yr, 5 × 10 5 yr and 1 × 10 6 yr).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 2 ) ≈ 10 4 cm −3 , which corresponds to a radius of ≈ 5200 AU in the theoretical profiles, offset from the position of the dust peak.For a time of t = 3.5 × 10 5 yr, the CH 3 OH abundance reaches its maximum value of approximately 4 × 10 −10 , which is in good agreement with values observed in various prestellar 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 matched approximately the observed one (Caselli et al. 1999), which occurred at a very early time of 3 × 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 ≈ 1 × 10 11 cm −2 , underproducing the observed values for methanol of 3.9 × 10 13 cm −2 at the dust peak and 5.9 × 10 13 cm −2 at the methanol peak by roughly two orders of magnitude.On the other hand, the N(CH 2 DOH)/N(CH 3 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 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 to have a low probability.
To estimate a new best-fit time for the fiducial model, we calculated the χ 2 -values of the observed central column density vs. the corresponding values in the fiducial model for different time steps and species.Subsequently the time step where the sum of the χ 2 -values of all species is minimised was taken as the bestfit time.For the 1D-4 model, this corresponds to a best-fit time of t = 3 × 10 5 yr, which is roughly consistent with estimations from other molecules (e.g.Redaelli et al. 2019aRedaelli et al. , 2021b)).Additionally, it coincides with the occurrence of the highest methanol column density in the probed time frame.For t = 3.0 × 10 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).2010) yielding static radial profiles of the H 2 density n(H 2 ) (blue, in logarithmic scale), the visual extinction A V (red), the gas temperature T gas (orange) and the dust temperature T dust (green).
time step (t = 1×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 time of t = 3 × 10 5 yr, the deuterium fraction assumes an almost completely flat profile with a value of ≈ 0.08.At the intermediate (t = 5 × 10 5 yr) and late (t = 1 × 10 6 yr) time steps, the N(CH 2 DOH)/N(CH 3 OH) ratio starts to decrease slightly, but stays at a level that is in agreement with the observed one.

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 Figure 6).We have 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 are using different chemical networks.MONACO is specialised to describe the formation of complex organic molecules, 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.Threephase 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 into 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.In the case that 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 being just 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 vertical direction.The fiducial model presented in this paper also uses a multilayer dust model, but with a chemically inactive bulk phase and only one layer in the surface phase, which reduces the surface area on which methanol and its precursors can be hydrogenated.
To evaluate the consequences of our choice to use only one layer in the chemically active surface phase in our models, we have 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 set up, 1D-4-2, where the entire ice is available for reactions.For 1D-4-1, we find that both the CH 3 OH and the CH 2 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 of 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, but also 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 Section 4.1.2),as the increase in column densities for CH 2 DOH is lower than for CH 3 OH.Surprisingly, this is not the case for the model with four layers in the surface phase, as is depicted in Figure 5. In- stead, we find a considerable increase of the deuterium fraction compared to the fiducial model for early time steps, reaching a value of up to ≈ 0.16 at t = 1 × 10 5 yr.The deuterium fraction of the model with four layers in the surface phase drops slightly below that of the fiducial model for time steps around when the methanol column density peaks, but increases again above it for later time steps.
Both models, the fiducial and V17, employ a similar form of diffusion.The diffusion-to-binding energy E d /E b is set to a constant value of 0.55, as suggested by Minissale et al. (2016a), for all surface species.Additionally, both codes consider the inclusion of diffusion by quantum tunneling.Although, while pyRate only follows the tunneling diffusion of H (and D), MONACO also traces that of H 2 .
The mechanism for reactive desorption in the two codes is similar: altered here for chemical reactions with multiple products, to ensure that it adheres to conservation of energy, and 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.Note that even though the procedure for single product reactions is the same in both codes, they can produce unequal re-active desorption efficiencies owing to the fact that we use 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 WebBook2 and KIDA3 .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 both codes also exist for some of the reactions involved in the hydrogenation chain towards methanol: For example for reaction 7 and 8 the two codes yield very different reactive desorption efficiencies.In pyRate, H 2 CO is desorbed more efficiently before it has the opportunity to react further and eventually form CH 3 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 2 CO, CH 3 OH and CH 2 DOH of models 1D-4, S16 and V17.In the case of V17, the CH 2 DOH column density was derived by scaling the CH 3 OH column density with the respective deuteration ratio from S16, as MONACO does not include a description of deuterium chemistry.Note that the column density profiles are shown at different time steps.The V17 model and the S16 model presented in Chacón-Tanarro et al. (2019) showed the time steps when the CO column density in the respective model is comparable to the observed value.For V17 this corresponds to t = 1.6 × 10 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 × 10 4 yr.The resulting column densities are shown in Figure 6 as the dark red line.We have derived a new best fit time corresponding to the lowest χ 2 value of the observed column densities vs. the corresponding values in the 1D-4 model, yielding t = 3.0 × 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 2. In Figure 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 3 OH column density and as consequence also CH 2 DOH column density, which was credited to a likely overestimation of the reactive desorption efficiency.All models overestimate the amount of gas phase H 2 CO.The V17 and S16 model, that are evaluated at a time step to match the observed CO column densities, are off by more than an order of magnitude, whereas S16 at t = 3.0 × 10 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 3 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 tunneling 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 e.g.: cosmic-ray induced desorption and photodesorption seem to have a negligible impact on the release of methanol and its precursors into the gas phase.

Discussion
In addition to the comparison with the V17 model, we have explored various alterations of our model, mainly to investigate the effect on the deuteration of methanol and possibly t = 10 6 yr Fig. 5. Modelled ratio between singly-deuterated methanol (CH 2 DOH) and non-deuterated methanol (CH 3 OH) for the 1D-4 Model for four different time steps ranging from 10 5 yr (top left) to 10 6 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 Section 3.2) for a more detailed explanation.)Thecoloured lines shows the column density ratio of the models, while the black line indicates the observed ratio (errors as grey-shaded areas).
improve the agreement between the predicted and observed N(CH 2 DOH)/N(CH 3 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 have 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 modifying the fiducial model 1D-4 are indicated by an asterisk( ⋆ ), whilst those derived from 1D-3 are highlighted by a dagger( † ).We have 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 have roughly divided them into chemical and physical variations.We collected an overview of the resulting N(CH 2 DOH)/N(CH 3 OH) ratios for for all presented models at four different time steps for the positions of the dust peak and of the methanol peak in Table 3.

1D-1 to 1D-6: Combining Enhanced Diffusion and Abstraction Reactions
As already pointed out in Section 3.2, to explain the observed magnitudes of molecular abundances and column densities, we have 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) by quantum tunneling.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 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 in the case that the diffusion proceeds mainly via quantum tunneling, could not be observed.Therefore, we have explored both the option for an increased, fast thermal diffusion as well as the diffusion of H and D atoms by quantum tunneling.
For the former, a value of the diffusion-to-binding energy E d /E b has to be chosen.Due to the lack of detailed experimental data, most chemical codes apply only one value for the diffusionto-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 assume 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 have adopted the lowest debated value of 0.2 for the diffusion-to-binding energy.As a reference, we have also tested the option to introduce no enhanced diffusion -neither fast thermal diffusion nor tunneling diffusion, but only to employ 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 ref-erences to the models with enhanced diffusion.Both produce CH 3 OH column densities around 10 7 cm −2 , which is several orders of magnitude lower than the observed value of 10 13 cm −2 .Models 1D-3 and 1D-4, additionally employing the diffusion via tunneling 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 −2 .These results deviate only within a factor of 10 from the observed values and are thereby matching the observation.Therefore, we conclude that some form of enhanced diffusion of H and D atoms over the grain surface has to take place in order to reach similar column densities as the ones measured in the pre-stellar core L1544.
The column densities in the models that employ fast thermal hopping are higher for early time steps (t = 1 × 10 5 yr) to intermediate time steps (t = 5 × 10 5 yr) -up to a factor of six for CH 3 OH and up to a factor of four for CH 2 DOH -, but decline earlier than in the models with tunneling diffusion.Moreover, we find that the CH 3 OH column density profiles in the 1D-4 model including only addition reactions are lower (by up to a factor of three) than for the 1D-3 model comprising both addition and abstraction reactions.The CH 2 DOH column density profiles, on  Notes.Listed are only chemical and physical properties that were varied between the models.The 1D-4 model is the fiducial model presented in Section 3.1.Models taking the 1D-4 model as starting point are marked with an asterisk( ⋆ ).Models taking the 1D-3 model, including the abstraction reactions, as a starting point are marked with a dagger( † ).  3. N(CH 2 DOH)/N(CH 3 OH) ratio at the centre of the pre-stellar core r cen and the radius of the methanol peak r max = 5200 AU.

Model
t = 1 × 10 5 yr t = 3 × 10 5 yr t = 5 × 10 5 yr t = 1 × 10 6 yr r cen r max r cen r max r cen r max r cen r max 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 the asterisk is the fiducial model and t = 3 × 10 5 yr is the best-fit time obtained as described in Section 3.1.The observed values are N(CH 2 DOH)/N(CH 3 OH)(r cen ) = 7.1% ± 1.9% for the centre and N(CH 2 DOH)/N(CH 3 OH)(r max ) = 4.5% ± 1.2% for a radius of 5200 AU.
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 ten for CH 3 OH and up to eight for CH 2 DOH.
In addition to the order of magnitude of column densities, we attempt to reproduce the deuterium fraction N(CH 2 DOH)/N(CH 3 OH) as closely as possible.This ratio is likely less affected by the individual modelling choices, as the effects on CH 3 OH and CH 2 DOH are probably similar for most parameter selections.Figure 7 depicts the ratio of the column densities of singly deuterated methanol (CH 2 DOH) and nondeuterated methanol (CH 3 OH).Models 1D-1 and 1D-2 show a similar level of deuteration reaching a central value of ≈ 0.05 for the best fit time of t = 3.0 × 10 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 2 CO proceeds preferably by the latter reaction: However, the deuteration of CH 3 O is not able to produce CH 2 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 3 OD in the presence of D. CH 2 DOH can only be produced by deuteration of CH 2 OH.Including abstraction reactions into the chemical network opens up another channel for the formation of CH 2 DOH via the abstraction of non-deuterated methanol: Reaction 11 is favoured in comparison to the analogue including its isomer 12.A similar effect on the deuterium fraction of the molecular abundances was discussed for the static 0D-models in Taquet et al. (2012).They found large enhancements for the N(CH 2 DOH)/N(CH 3 OH) ratio in the high density case (n H ≥ 1 × 10 6 cm −3 ) eventually reaching values above unity, while there was no significant increase for the low density case (n H = 1 × 10 4 cm −3 -1 × 10 5 cm −3 ).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 tunneling diffusion -show the opposite behaviour when it comes to the inclusion of abstraction reactions (see Figure 7).The models with only addition, 1D-4 and 1D-6, have an up to eight times higher N(CH 2 DOH)/N(CH 3 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 2 DOH apparently causes a decrease of the amount of CH 3 OH that is produced.Models 1D-3 and 1D-5, comprising addition and abstraction reactions, are not able to reproduce the observed N(CH 2 DOH)/N(CH 3 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 5 yr), but decline too quickly with increasing radius.For late time steps (t = 1×10 6 yr), the deuterium fraction is increasing again to a level, where the intermediate radii (r=2000AU -r=6000) match the observational profile.The innermost part, however, shows less deuteration than the outer ring and also the observations.By looking at the reaction rates of for example model 1D-3, one can see that at a time step close to the peak in methanol formation (≈ t = 5 × 10 5 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 slower.We therefore suspect that the reaction rates for the abstraction reactions are too high/much higher than in reality.
Having assessed that some form of enhanced diffusion is needed, it is not possible to decide based on our modelling results, which of the processes of enhanced diffusion matches the observations better.For times before the time step when the methanol column density peaks, the fiducial model, employing tunneling diffusion, has a higher N(CH 2 DOH)/N(CH 3 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 5 yr.After the time step when the methanol column density peaks, the N(CH 2 DOH)/N(CH 3 OH) ratio of the fiducial model decreases and begins to match the observed values quite closely for time steps beyond t = 5.0 × 10 5 yr.The N(CH 2 DOH)/N(CH 3 OH) ratio of the 1D-6 model starts increasing far above the observed values, until it reaches values of almost 0.14 for late times (t = 1 × 10 6 yr).

1D-7: Variation of the Grain Model
As described in Section 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) have shown 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 2 DOH)/N(CH 3 OH) ratio are depicted in Figure 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 −3 to 10 −4 compared to 10 −2 in the fiducial model, in model 1D-3 and in the observed N(CH 2 DOH)/N(CH 3 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 3 OH.They are of the order of 10 13 cm −2 , differing only a factor of a few from the observed values.The column densities of CH 2 DOH on the other hand are one or two orders of magnitude lower than in the three phase model and CH 2 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.

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 Section 2.1, allowed only the reactive desorption of exother-mic surface reactions with a single reaction product.Additionally, the mechanism applied the same reactive desorption efficiency, typically 1%, to every eligible reaction and did not distinguish between surface types.In fact that is the case for most other reactive desorption mechanisms used in the literature, except for the one presented here and the one in Vasyunin et al. (2017).
The new reactive desorption mechanism extends its application to chemical reactions with two reaction products.It is now interesting to quantify the consequences of this change.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 less 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: t = 10 6 yr Fig. 8. Modelled column density profiles of non-deuterated methanol (CH 3 OH) and singly-deuterated methanol (CH 2 DOH) for several models at four different time steps ranging from 10 5 yr to 10 6 yr.We present the varied models showing the largest effects on the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product) and 1D-12 (radius-dependent cosmic-ray ionisation rate).The black lines show the observed profiles (errors as grey-shaded areas).The solid lines indicate the CH 3 OH column densities, the dashed lines indicate the CH 2 DOH column densities respectively.
the new reactive desorption mechanism yields on a CO-surface an efficiency of 54% for the H 2 and of 1.5×10 −36 % for HCO.The abstraction reactions cause the expulsion of significant amounts of H 2 , HD and D 2 , while the desorption of the larger reaction partner is negligible.
The column densities of CH 3 OH and CH 2 DOH in the 1D-8 model are increased compared to the 1D-3 model (see Figure 8).The difference between the two models grows with time, beginning with a factor of ≈ 2 at early time steps (t = 1 × 10 5 yr) to a factor of ≈ 16 at late time steps (t = 1 × 10 6 yr) for CH 3 OH or from ≈ 2.5 to ≈ 12 for CH 2 DOH respectively.This results in a N(CH 2 DOH)/N(CH 3 OH) profile (see Figure 9) which 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 2 DOH column density increases more quickly than the one of CH 3 OH in the 1D-8 model.This could indicate that the determined reactive desorption efficiencies do not describe the physical reality very well.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 2 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 keep our assumptions simple and as general as possible.Further work could investigate different options for the partitioning of the available reaction enthalpy.t = 10 6 yr Fig. 9. Modelled ratio between singly-deuterated methanol (CH 2 DOH) and non-deuterated methanol (CH 3 OH) for several models at four different time steps ranging from 10 5 yr to 10 6 yr.We present the varied models showing the largest effects on the amount of methanol in the gas phase: 1D-7 (two-phase model), 1D-8 (only reactive desorption with one product) and 1D-12 (radius-dependent cosmic-ray ionisation rate).The black line shows the observed ratio (errors as grey-shaded areas).

1D-9: Variation of the Activation Energy
The formation of methanol proceeds by the successive hydrogenation of CO along HCO, H 2 CO and CH 2 OH/CH 3 O to CH 3 OH.Reactions 14 to 16: have an activation barrier.Their corresponding activation energies E A are indicated on top of the arrows.The remaining reactions are barrierless.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 2 DOH)/N(CH 3 OH) ratio.Therefore we explore how the results are affected by a small variation of the activation energy.Specifically, we aim to test if decreasing the activation energies for reactions producing deuterated isotopologues could lead to a significant increase of the N(CH 2 DOH)/N(CH 3 OH) ratio.Hence, we decided to decrease the activation energy for those reactions by 200 K.We only varied the activation energy of specific addition reactions (marked with an asterisk in the overview in B).The abstraction reactions were left untouched.
Undertaking this variation produces a slightly higher CH 2 DOH column density (see Figure C.1).For early times (t=1 × 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 × 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 × 10 5 yr) and late (t = 1 × 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 Figure C.1).Consequently, the N(CH 2 DOH)/N(CH 3 OH) profile (see Figure C.2) experiences an increase as well.For early time steps the N(CH 2 DOH)/N(CH 3 OH) ratio reaches high values of up to 16.6 at t=1.15 × 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 × 10 6 yr), hitting the area of uncertainty, but presenting only a small upward shift compared to the 1D-3 model.

1D-10: Variation of the Formation Enthalpies
The incorporation of the more sophisticated reactive desorption mechanism required to expand the list of tabulated formation enthalpies H form , which are necessary to compute the reaction enthalpies ∆H.The complete list of formation enthalpies (and binding energies) of species involved in exothermic surface reactions is shown in Appendix A. Since experimental data for deuterated isotopologues is 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 Web-Book4 .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 have a run a model where we adopt the same values for non-deuterated and deuterated isotopologues.The column density profiles for CH 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 to 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 Figure 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 Figure 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 2 DOH)/N(CH 3 OH) ratios become very similar.4.2.2.1D-12: Variation of the Cosmic-Ray Ionisation Rate UV photons are not able to penetrate the inner, denser parts of molecular clouds with visual extinctions A V ≥ 1, as they are already efficiently absorbed by the outer layers of the cloud.Therefore, cosmic rays take their place as the main ionising agents in the central parts, constituting the start of the ion-molecule chemistry.The penetrating cosmic rays will ionise molecular hydrogen to form 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: The direction of the reaction from left to right is strongly favoured, due to the lower zero-point energy of H 2 D + , for temperatures below 30 K, that 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 + 3 is mostly frozen out on dust grains in the very inner part of the pre-stellar core.These convenient conditions promote a very efficient deuteration process in the inner parts of the core, which also quickly translates the high D/H ratios to more complex molecules.
While we have used the canonical value of the cosmic-ray ionisation rate per hydrogen molecule ζ(H) = 1.3 × 10 −17 s −1 for most of the models presented here, for the 1D-12 model we have 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 × 10 −17 s −1 in the innermost cell to ζ(H) = 4.89 × 10 −17 s −1 at the outer boundary of the core.
The 1D-12 model exhibits higher CH 3 OH and CH 2 DOH column densities as compared to the 1D-3 model (see Figure 8).Especially for early times (t = 1 × 10 5 yr), the CH 3 OH column densities are a factor of 13 higher than in the 1D-3 model.However, the difference to the 1D-3 model is decreasing over time: at intermediate times (t = 5 × 10 5 yr) it is approximately a factor of two and even lower at late times (t = 1 × 10 6 yr).The shape of the column density profile (see Figure 9) is almost identical between the two models, with small deviations at early times.For the increase of the singly deuterated methanol, we see a timedelayed behaviour compared to the non-deuterated isotopologue.CH 2 DOH is amplified by a factor of three at early times (t = 1 × 10 5 yr), but it never gets to the values observed for CH 3 OH.Nevertheless, the amplification for CH 2 DOH overtakes the one for CH 3 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 2 DOH)/N(CH 3 OH) ratio of ≈ 0.74 is reached at late times (t = 1 × 10 6 yr).

1D-13: Variation of the Cosmic-Ray Desorption Mechanism
A frequently adopted model for the cosmic-ray induced desorption (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-70MeV nucleon −1 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: 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 e.g.: f (a, T max ) = 1×10 −5 s/3.16×10 13 s = 3.16×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 CR 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 have chosen the CR 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×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 cosmic-ray induced desorption.The CH 3 OH column density profile (see Figure C.1) is slightly decreased, especially for the earlier time steps , while the CH 2 DOH profile is not significantly changed, resulting in a little higher N(CH 2 DOH)/N(CH 3 OH) ratio (see Figure C.2).It reaches its highest value of 0.088 at t = 1 × 10 5 yr in the very inner centre.However, the decline is again much steeper than is observed.

Conclusion
We presented several models for the prediction of column densities and deuterium fractions of methanol and its deuterated isotopologues in pre-stellar cores.As a comparison to observed quantities, we use single-dish observations of H 2 CO and CH 3 OH and some of their deuterated isotopologues towards the prestellar core L1544 conducted and analysed by Chacón-Tanarro et al. (2019).
All introduced models use a novel treatment of reactive desorption of molecules from the surface of interstellar dust grains.The treatment is based on experimental justification (Minissale et al. 2016b) and derives an individual reactive desorption efficiency for every species, depending on the forming chemical reaction(s), and on the type of underlying surface.
Our fiducial model serves as a comparison to the results of the models V17 (MONACO) and S16 (pyRate), presented in Chacón-Tanarro et al. (2019).It includes thermal diffusion (diffusion-to-binding energy ratio E d /E b = 0.55) as well as the diffusion of hydrogen and deuterium atoms by quantum tunneling.The chemical network does not comprise abstraction reactions for the methanol reaction scheme.We estimate a best fit time, which coincides with the occurrence of the methanol peak at 3.0 × 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 3 OH and CH 2 DOH are still underestimated.However, instead of 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 2 DOH)/N(CH 3 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 vertical direction and as is done in the V17 model, has an increasing effect on 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) showed on the one hand the necessity to employ an increased rate of surface diffusion and on the other hand to disregard abstraction reactions from the reaction scheme, in order to reproduce the observed order of magnitude for the methanol abundances.Therefore, we have 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 by quantum tunneling.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 several orders of magnitudes below the observed ones.However, we can not decide based solely on our models which enhanced diffusion process -fast thermal hopping or tunneling diffusion -matches the observed N(CH 2 DOH)/N(CH 3 OH) ratio profiles better.Both produce D/H ratios of a similar level, which are in good agreement with the observed values.For the best fit time, the two options have nearly identical profiles.
Also, we have 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 only including 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 have explored other modifications to our model, that we were suspecting to have an effect on the N(CH 2 DOH)/N(CH 3 OH) ratio: -two-phase model: higher CH 3 OH and CH 2 DOH column densities; deuterium fraction is severely underestimated by more than one order of magnitude -reactive desorption applied only to reactions with one reaction product: higher CH 3 OH and CH 2 DOH column densities, resulting in a slightly better agreement with the observations -location-dependent cosmic-ray ionisation rate: higher CH 3 OH and CH 2 DOH column densities only differing with 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 models: decrease of the activation energy by 200 K leading to deuterated isotopologues compared 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 of 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 to 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 mechanism as the Eley-Rideal mechanism, which is not treated by many chemical codes, or non-diffusive chemistry could play an important role for 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 by quantum tunneling.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.

Appendix A: Table of 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, on which we applied the newly developed reactive desorption mechanism laid out in Section 2.1.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).

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

Figure 5 Fig. 2 .
Fig.2.Physical model developed inKeto & Caselli (2010) yielding static radial profiles of the H 2 density n(H 2 ) (blue, in logarithmic scale), the visual extinction A V (red), the gas temperature T gas (orange) and the dust temperature T dust (green).

Fig. 3 .
Fig. 3. Gas phase abundance profiles of H 2 CO, CH 3 OH and CH 2 DOH for the 1D-4 model for four different time steps ranging from 10 5 yr (left) to 10 6 yr (right).The best fit time is t = 3 × 10 5 yr.

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

Fig. 6 .
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 5 yr for V17, t = 3.0 × 10 4 yr and t = 3.0 × 10 5 yr for S16 and t = 3.0 × 10 5 yr for 1D-4.The observed profiles are depicted in black (errors as grey-shaded areas).The CO column densities were obtained by using observations of C 17 O and the average isotopic ratios of 16 O/ 18 O = 557 and 18 O/ 17 O = 3.6 in the local ISM(Wilson et al. (1999)).
(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) Table

Fig. 7 .
Fig. 7. Modelled ratio between singly-deuterated methanol (CH 2 DOH) and non-deuterated methanol (CH 3 OH) for several models at four different time steps ranging from 10 5 yr to 10 6 yr.The best fit time is t = 3 × 10 5 yr.We show two models with slow thermal hopping (E d /E b = 0.55) (blue lines), two models with slow thermal hopping and tunneling diffusion of H and D (red lines) and two models with fast thermal hopping (E d /E b = 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).
3 OH and CH 2 DOH of this model are shown in Figure C.1 and the deuterium fraction profiles are shown in Figure C.2.The effect on both the CH 3 OH and CH 2 DOH column densities and on the N(CH 2 DOH)/N(CH 3 OH) ratio is vanishingly small.The difference with respect to the reference model, 1D-3, ranges around 4% for the best fit time (t = 3.0 × 10 5 yr.) 4.2.Physical variation 4.2.1.1D-11: Variation of the Gas Temperature

Table 1 .
Initial chemical abundances with respect to n H .

Table 2 .
Overview of the various models investigated in this work.
Table A.1.Tabulated formation enthalpies H form and binding energies E b .