# Formation of Complex Organic Molecules in Cold Interstellar Environments
through non-diffusive grain-surface and ice-mantle chemistry

Mihwa Jin  Department of Chemistry, University of Virginia, Charlottesville,
VA, 22904, USA  Robin T. Garrod  Department of Chemistry, University of
Virginia, Charlottesville, VA, 22904, USA  Department of Astronomy, University
of Virginia, Charlottesville, VA, 22904, USA

###### Abstract

A prevailing theory for the interstellar production of complex organic
molecules (COMs) involves formation on warm dust-grain surfaces, via the
diffusion and reaction of radicals produced through grain-surface
photodissociation of stable molecules. However, some gas-phase O-bearing COMs,
notably acetaldehyde (CH3CHO), methyl formate (CH3OCHO), and dimethyl ether
(CH3OCH3), are now observed at very low temperatures, challenging the warm
scenario. Here, we introduce a selection of new non-diffusive mechanisms into
an astrochemical model, to account for the failure of the standard diffusive
picture and to provide a more generalized scenario of COM formation on
interstellar grains. New generic rate formulations are provided for cases
where: (i) radicals are formed by reactions occurring close to another
reactant, producing an immediate follow-on reaction; (ii) radicals are formed
in an excited state, allowing them to overcome activation barriers to react
with nearby stable molecules; (iii) radicals are formed through photo-
dissociation close to a reaction partner, followed by immediate reaction. Each
process occurs without the diffusion of large radicals. The new mechanisms
significantly enhance cold COM abundances, successfully reproducing key
observational results for prestellar core L1544 (catalog ). H-abstraction from
grain-surface COMs, followed by recombination, plays a crucial role in
amplifying chemical desorption into the gas phase. The UV-induced chemistry
produces significant COM abundances in the bulk ices, which are retained on
the grains and may persist to later stages. O2 is also formed strongly in the
mantle though photolysis, suggesting cometary O2 could indeed be interstellar.

Astrochemistry (75) – Interstellar dust processes (838) – Star formation
(1569) – Molecule formation (2076)

††journal: ApJS

##  1 Introduction

Complex organic molecules (COMs), usually defined as carbon-bearing molecules
with six or more atoms, have been detected within the interstellar medium and
in various protoplanetary environments (Blake et al., 1987; Fayolle et al.,
2015; Bottinelli et al., 2010; Arce et al., 2008; Öberg et al., 2010; Bacmann
et al., 2012). COMs synthesized at the early stages of star formation are
suggested to have been a starting point for the organic materials that went on
to seed the nascent solar system (Herbst & van Dishoeck, 2009). While the
degree to which the interstellar synthesis of COMs may contribute to pre-
biotic/biotic chemistry on Earth is still a matter of debate, many recent
studies have shed light on the possible interstellar/protostellar origins of
chemical complexity. For instance, the sugar-like molecule glycolaldehyde
(CH2(OH)CHO) has been detected toward the class 0 protostellar binary source
IRAS16293-2422 (Jørgensen et al., 2012), as well as the Galactic Centre source
Sgr B2(N) and other hot cores (Hollis et al., 2000; Beltrán et al., 2009). A
related molecule, ethylene glycol ($\ce{HOCH$_{2}$CH$_{2}$OH}$), has also been
found toward low-mass protostars (Maury et al., 2014; Jørgensen et al., 2016),
having first been detected in SgrB2(N) (Hollis et al., 2002). The simplest
amino acid, glycine ($\ce{NH2CH2COOH}$), has not been detected in the
interstellar medium, but numerous amino acids have been found in meteorites
(Kvenvolden et al., 1970), and glycine has been detected in two comets: Wild 2
(Elsila et al., 2009) and Churyumov-Gerasimenko (Altwegg et al., 2016). In the
laboratory, (Meinert et al., 2016) observed the generation of numerous sugar
molecules, including the aldopentose ribose, through photochemical and thermal
processing of interstellar ice analogs initially composed of
$\ce{H2O,CH3OH,andNH3}$.

One of the prevailing theories explaining the formation of interstellar COMs
is based on radical chemistry occurring on dust-grain particles, following the
formation of ice mantles on their surfaces. In dark clouds and prestellar
cores, where temperatures may be as low as around 10 K, hydrogen atoms are the
most mobile species on the grain. Their high mobility allows them to find
reaction partners on the surface easily. Thus, the accretion of other atoms
and simple molecules, such as O, C, N and CO, onto the grain surfaces often
leads to reactions with atomic H, producing an ice mantle composed of
commonly-observed stable hydrides like water, methane, and ammonia, through
repetitive hydrogenation. In dense regions, almost all of the CO produced in
the gas phase may ultimately freeze out onto the dust grains (e.g. Caselli et
al., 1999), producing an outer ice surface that is rich in CO. Atomic H may
also react with CO to produce methanol (CH3OH); although there are activation
energy barriers that prevent all of this CO being hydrogenated, a substantial
fraction of the ice – up to as much as around 30% with respect to water toward
some protostars, but generally of the order of 5% (Boogert et al., 2015) – is
found observationally to be composed of $\ce{CH3OH}$.

Grain-surface methanol, and indeed other molecules, may be broken down into
radicals through UV-induced photolysis, caused either by external photons or
by the secondary UV field induced by cosmic-ray collisions with gas-phase H2;
however, in the standard picture used in most astrochemical models, these
radicals are still immobile on the grain surfaces at low temperatures,
hindering their ability to react with anything other than mobile H atoms. As a
star-forming core evolves to be heated to temperatures greater than $\sim$20
K, radicals such as CH3 and HCO become more mobile, allowing them to meet via
diffusion and quickly react to form larger COMs (Garrod & Herbst, 2006; Garrod
et al., 2008). At this stage, most of the COMs produced in this way are unable
to desorb effectively into the gas phase where they may be directly detected;
only non-thermal mechanisms such as chemical desorption (Garrod et al., 2007)
and photo-desorption (Öberg et al., 2009a, b) are active, with the former
expected to yield only around 1% of reaction products to the gas phase, while
the latter would be weaker still in such regions where visual extinction is
high. Only when the dust-grains in the core reach a high enough temperature
(typically on the order of $\sim$100 K), through protostellar heating, can
large COMs sublimate efficiently from the grains. This process is reviewed in
more detail by Öberg et al. (2010).

According to this scenario, it would be unlikely that COMs other than methanol
would be effectively synthesized under the low-temperature conditions of dark
clouds and prestellar cores, due to the immobility of heavy species on the
grains, while methanol would still be formed as the result of the diffusion
and reaction of H alone. The detection of COMs in these cold environments
therefore presents a significant challenge to current astrochemical models
based on diffusive grain-surface chemistry. In particular, three O-bearing
COMs, acetaldehyde (AA; CH3CHO), methyl formate (MF; CH3OCHO), and dimethyl
ether (DME; CH3OCH3), have been detected in prestellar cores such as L1689B,
L1544 (catalog ) and B1-b (Bacmann et al., 2012; Vastel et al., 2014;
Cernicharo et al., 2012) as well as in the cold outer envelopes of protostars
(Öberg et al., 2010; Bergner et al., 2017).

Many chemical mechanisms have been proposed to explain the presence of COMs in
cold cores: The Eley-Rideal process combined with complex-induced reactions on
the grains (Ruaud et al., 2015); non-diffusive reactions (Chang & Herbst,
2016); oxygen insertion reactions (Bergner et al., 2017); cosmic-ray induced
chemistry (Shingledecker et al., 2018); and gas-phase formation through
radical reactions (Balucani et al., 2015). In particular, the latter gas-phase
formation scheme is noteworthy in that it suggests a chemical link between
gas-phase methyl formate and dimethyl ether via the radical $\ce{CH3OCH2}$,
which would be formed through the abstraction of a hydrogen atom from DME by
atoms such as F and Cl. The dimethyl ether itself would form through the
efficient radiative association of the radicals CH3 and CH3O in the gas phase.

Although many of the above-mentioned studies provide plausible routes to one
or more complex molecules, some suffer from poorly constrained parameters due
to the lack of experimental data or detailed computational chemistry studies
(e.g. the radiative association reaction to produce DME), or require somewhat
extreme chemical conditions. For example, the fractional abundance of methoxy
in the model of Balucani et al. (2015) is more than one order of magnitude
larger than the observed value, making it unclear whether dimethyl ether can
be formed solely in the gas phase in cold cores (Vasyunin & Herbst, 2013;
Chang & Herbst, 2016). Also, the chemical networks for F and Cl are relatively
sparse, introducing further uncertainties; the models of Vasyunin et al.
(2017) suggest that OH is instead a more important agent of H-abstraction from
dimethyl ether.

Cosmic ray-induced radiolysis models by Shingledecker et al. (2018) appear
promising, although it is unclear whether they can account for large
abundances of COMs on their own, or whether they require radicals also to be
produced by UV photons. For example, those authors suggest an average rate of
solid-phase water dissociation by cosmic rays of around $6\times 10^{-16}$ s-1
at the canonical cosmic-ray ionization rate. However, the typically-assumed
rate used in astrochemical models for the separate process of UV photo-
dissociation of gas-phase water by cosmic ray-induced photons is around
$1.3\times 10^{-14}$ s-1. Even if the rates of UV photo-dissociation in the
solid phase are around three times lower than those of gas-phase molecules
(Kalvāns, 2018), the influence of cosmic ray-induced photons in producing
radicals in ice mantles would, on average, appear to be greater than that of
direct cosmic-ray impingement by close to an order of magnitude. Other
calculations for CR-induced water dissociation (Garrod, 2019) would suggest an
even greater discrepancy. The formation of COMs through direct cosmic-ray
impacts into icy grains may therefore rely on the presence of pre-existing
radicals produced by CR-induced UV.

In order for chemical modeling studies to be directly comparable with
observations, they should ideally seek to reproduce the spatio-physical
conditions of the target source. In pre-stellar cores such as L1544 (catalog
), strong gradients in physical parameters such as dust temperature, density,
and visual extinction may sample a range of chemical regimes. The recent
detection of a chemically-active outer shell around the core center of L1544
(catalog ) sheds light on the complexity of its chemical structure (Bizzocchi
et al., 2014). Methanol column density in this source appears to peak at a
radius around 4000 au distant from the source center. Abundances of other COMs
seem to be enhanced at this position also (Jiménez-Serra et al., 2016).

The modeling study by Vasyunin et al. (2017) is notable in its use of an
explicit spatio-physical model to perform chemical modeling of L1544 (catalog
), successfully reproducing a similar feature to the observed methanol peak.
Those authors used a surface-specific treatment of chemical desorption from
grains in which a pure water-ice surface would allow only minimal desorption
of newly-formed molecules, while a more CO-rich surface would allow molecules
including methanol to desorb more easily. However, the peak fractional
abundances of methanol produced in this model appear high ($\sim$$10^{-7}$
with respect to total hydrogen) compared with the usual values observed in
cold sources ($\sim$$10^{-9}$). The efficiency of the chemical desorption of
methanol (which is responsible for its gas-phase abundance) and many other
species on water ice is poorly constrained by experiment, with an upper limit
of 8% for methanol formed through H addition (Minissale et al., 2016c); thus,
there is room for the efficiency of this mechanism to be varied in models to
produce the correct observed quantity. However, the Vasyunin et al. (2017)
models rely on gas-phase methanol as the feedstock for the production of
other, larger COMs. Lower peak abundances of gas-phase methanol would almost
certainly render the gas-phase routes investigated by those authors too weak
to account for the abundances of DME, MF and AA. Furthermore, using gas-phase
reactions between CH and CH3OH as a mechanism to form CH3CHO, they erroneously
extrapolated experimental reaction rates determined at room temperature
(Johnson et al., 2000) down to 10 K; this produces impossibly large rate
values that are around 3 orders of magnitude greater than the collisional
rates for neutrals at such temperatures. It is unclear how much influence
their alternative route for AA production, the radiative association of CH3
and HCO radicals, would have on that molecule’s peak gas-phase abundance;
however, in their network they assume a rate that is on the order of the
collisional rate (i.e. high efficiency), based on high-pressure experiments
that would in fact be more representative of collisional (rather than
radiative) de-excitation, and which therefore say nothing certain about
radiative association. Until detailed calculations for this, or indeed for the
radiative association of CH3O with CH3, are available, those rates may be
considered to be somewhat optimistic.

Recent laboratory results have added new impetus to understanding the behavior
of grain-surface chemistry at very low temperatures. Fedoseev et al. (2015,
2017) considered the co-deposition of CO and H at low temperatures, resulting
in significant production of complex organics such as glycolaldehyde, ethylene
glycol, and glycerol. They postulate that reactions between adjacent HCO
radicals, and/or HCO and CO, may produce such species with minimal thermal
diffusion required – a result not predicted by standard gas-grain chemical
models, due to the lack of diffusion involved. So-called microscopic Monte
Carlo models, unlike the more usual rate-equation models, are capable of
simulating the relative positions of all atoms and molecules on a grain
surface at each moment, allowing them automatically to account for this non-
diffusive process. Such models seem to indicate that COM production through
non-diffusive mechanisms is plausible (Fedoseev et al., 2015; Chang & Herbst,
2016; Dulieu et al., 2019, Ioppolo et al. 2020, subm., Garrod et al. 2020, in
prep.).

It would appear, therefore, that – regardless of other mechanisms such as gas-
phase processes or radiolysis – the standard rate description of grain-surface
reactions is insufficient to treat all temperature regimes; there are
situations in which reactants may be produced and rapidly meet (and react),
either without diffusion of radicals or with some minimal amount of diffusion
that does not obey the more general rate treatment.

In this study, we present a relatively simple formulation for non-diffusive
chemistry, for use in standard gas-grain chemical models, that allows a newly-
formed reaction product to react further with some other chemical species that
happens to be in close proximity to the product(s) of the first reaction. Due
to the instantaneous nature of this process, we refer to it here as a “three-
body reaction mechanism” (3-B). We also consider a similar, related mechanism
in which the new product has sufficient excitation energy to allow it to
overcome the activation-energy barrier to its reaction with some nearby
species; here, specifically with CO or H2CO. This mechanism is referred to
here as the “three-body excited formation mechanism” (3-BEF).

A functionally-similar mechanism was included in the recent model of the
solid-phase chemistry of cometary nuclei by Garrod (2019), but rather with the
initiating process being the production of radicals by UV-induced photo-
dissociation. In preliminary versions of that model that did not include that
mechanism, it was found that the photo-dissociation of bulk ice molecules at
temperatures 5–10 K was capable of producing implausibly high abundances of
reactive radicals, which were unable to react due to the lack of bulk thermal
diffusion at those temperatures. The Garrod (2019) model includes a new
reaction process whereby a newly-formed photo-product may react immediately
with a nearby reaction partner in the ice. This mechanism is included in the
present model also, which we label as the “photodissociation-induced reaction
mechanism” (PDI).

Garrod & Pauly (2011) used a similar mechanism to the 3-body mechanisms to
explain the formation of $\ce{CO2}$ ice at low temperatures. In their
treatment, the production of an OH radical via the reaction of H and O atoms
in proximity to a CO molecule could allow the immediate formation of CO2
(overcoming a modest activation energy barrier). Their models successfully
reproduced the observed behavior of CO, $\ce{CO2}$, and water ice in dense
clouds, and showed that such non-diffusive processes could be handled within a
standard gas-grain chemical model. More recently, Chang & Herbst (2016)
implemented a similar process, which they called a “chain reaction mechanism”
in their microscopic Monte Carlo simulation, achieving abundances of gas-phase
COMs high enough to reproduce the observational values toward cold cores (at a
temperature of $\sim$10 K), using a chemical desorption efficiency of 10$\%$
per reaction. Dulieu et al. (2019), seeking to explain the surface production
of NH2CHO in laboratory experiments involving H2CO, NO and H deposition,
introduced a non-diffusive reaction treatment for a single reaction (see also
Sec. 2.4).

Finally, a simple treatment for the Eley-Rideal process is included in the
present model, in which an atom or molecule from the gas-phase is accreted
directly onto a grain-surface reaction partner, resulting in immediate
reaction (mediated by an activation-energy barrier, where appropriate). Such
processes have been included in similar models before (e.g. Ruaud et al.,
2015), but are included here for completeness in the consideration of all
mechanisms by which reactants may instantly be brought together without a
mediating diffusion mechanism.

The formulations presented here for the above processes also allow, through a
repetitive application of the main non-diffusive reaction process (i.e. the
three-body mechanism), for the products of each of those processes themselves
to be involved in further non-diffusive reaction events (in cases where such
processes are allowed by the reaction network). Thus, for example, an Eley-
Rideal reaction may be followed by an immediate secondary non-diffusive
reaction. The importance of such repetitive processes will diminish with each
iteration.

All of the above non-diffusive reaction mechanisms are considered in the
model, with a particular emphasis on the production of the O-bearing COMs that
are now detected in the gas-phase in cold prestellar cores. The formulations
corresponding to each of the new mechanisms presented here are functionally
similar to each other, but quite different from the standard diffusive
reaction formula used in typical astrochemical models. However, they are fully
compatible with the usual treatment and may be used in tandem with it.

With the introduction of the new mechanisms, we run multi-point chemical
models of prestellar core L1544 (catalog ) to test their effectiveness in an
appropriate environment. A spectroscopic radiative-transfer model is
implemented here as a means to evaluate the observable column densities of
molecules of interest, allowing the direct comparison of the model results
with observations.

This paper is structured as follows. The chemical model and the newly-
implemented mechanisms are described in § 2. The results of the models are
explored in § 3, with discussion in § 4. Conclusions are summarized in § 5.

##  2 Chemical model

We use the astrochemical kinetic code MAGICKAL to study new grain-surface/ice-
mantle mechanisms that may effectively form complex organic molecules in cold
environments. The model uses a three-phase approach, based on that described
by Garrod & Pauly (2011) and Garrod (2013), in which the coupled gas-phase,
grain/ice-surface and bulk-ice chemistry are simulated through the solution of
a set of rate equations. Initial chemical abundances used in the model are
shown in Table 1, with elemental values based on those used by Garrod (2013).
The initial H/H2 abundances are chosen to agree approximately with the steady-
state values appropriate to our initial physical conditions, as determined by
the chemical model. Bulk diffusion is treated as described by Garrod et al.
(2017). Although the bulk ice is technically chemically active in this model,
at the low temperatures employed in this work, diffusive reactions in general
are negligibly slow, excluding processes involving H or H2 diffusion. However,
the addition of non-diffusive reactions to the model increases significantly
the degree of chemical activity within the bulk.

Table 1: Initial elemental and chemical abundances

Species, i |  $n(i)/n_{\textrm{H}}^{a}$ aaBoogert et al. (2015); values correspond to low-mass YSOs.
---|---
$\ce{H}$ | 5.0(-4)
$\ce{H2}$ | 0.49975
$\ce{He}$ | 0.09
$\ce{C}$ | 1.4(-4)
$\ce{N}$ | 7.5(-5)
$\ce{O}$ | 3.2(-4)
$\ce{S}$ | 8.0(-8)
$\ce{Na}$ | 2.0(-8)
$\ce{Mg}$ | 7.0(-9)
$\ce{Si}$ | 8.0(-9)
$\ce{P}$ | 3.0(-9)
$\ce{Cl}$ | 4.0(-9)
$\ce{Fe}$ | 3.0(-9)

Note. — ${}^{a}A(B)=A^{B}$

The model uses the modified-rate treatment for grain-surface chemistry
presented by (Garrod, 2008), which allows the stochastic behavior of the
surface chemistry to be approximated; the back-diffusion treatment of Willis &
Garrod (2017) is also used. Surface diffusion barriers ($E_{\mathrm{dif}}$)
are related to desorption (i.e. binding) energies ($E_{\mathrm{des}}$) such
that $E_{\mathrm{dif}}$=0.35$E_{\mathrm{des}}$ for all molecular species, with
bulk diffusion barriers taking values twice as high (Garrod, 2013). However,
the recent study by Minissale et al. (2016) estimated surface diffusion
barriers for atomic species such as N and O to be
$E_{\mathrm{dif}}$=0.55$E_{\mathrm{des}}$. We adopt a similar value
$E_{\mathrm{dif}}/E_{\mathrm{des}}$=0.6 for all atomic species; the impact of
this parameter is discussed in §4.3. These basic surface diffusion and
desorption parameters are also adjusted according to the time-dependent
abundance of H2 in the surface layer, following the method of Garrod & Pauly
(2011). All diffusion is assumed to be thermal, with no tunneling component
(see also Sec. 4.3).

Chemical desorption, whereby grain-surface reactions allow some fraction of
their products to desorb into the gas phase, is treated using the RRK
formulation of Garrod et al. (2007) with an efficiency factor
$a_{\mathrm{RRK}}=0.01$.

The grain-surface/ice-mantle photodissociation rates used in MAGICKAL are
based on the equivalent gas-phase rates, in the absence of other evidence
(e.g. Garrod et al., 2008), and they likewise assume the same product
branching ratios. Following the work of Kalvāns (2018), we adopt
photodissociation rates on the grain-surfaces and in the ice mantles that are
a factor of 3 smaller than those used for the gas phase. Photodissociation may
be caused either by external UV, or by the induced UV field caused by cosmic-
ray collisions with H2 molecules, and both sources of dissociation are
included in the model.

Methanol in cold clouds is mainly formed on the grain surfaces through ongoing
hydrogenation of CO (Fuchs et al., 2009). The methanol production network used
in the present network follows from that implemented by Garrod (2013), and
includes not only forward conversion of CO to methanol but also the backward
reactions of each intermediate species with H atoms.

The overall chemical network used here is based on that of Garrod et al.
(2017), with a few exceptions. In particular, a new chemical species, CH3OCH2,
has been added along with a set of associated gas-phase and grain-surface
reactions/processes listed in Tables 11 (gas-phase) and 12 (solid-phase). This
radical is a key precursor of DME in our new treatments (see §2.5). Its
inclusion also allows the addition of a grain-surface H-abstraction reaction
from DME, making this species consistent with methyl formate and acetaldehyde.

An additional reaction was included in the surface network, corresponding to
H-abstraction from methane (CH4) by an O atom ($E_{A}$=4380 K; Herron & Huie,
1973), as a means to ensure the fullest treatment for the CH3 radical in the
network.

The final change to the network is the adjustment of the products of
$\ce{CH3OH}$ photo-desorption to be $\ce{CH3}$+$\ce{OH}$ rather than
$\ce{CH3OH}$, roughly in line with the recent experimental study by Bertin et
al. (2016).

Each of the generic non-diffusive mechanisms that we include in the model (as
described below) is allowed to operate on the full network of grain-surface
and ice-mantle reactions that are already included for the regular diffusive
mechanism; reactive desorption, where appropriate, is also allowed to follow
on from each of these, in the case of surface reactions. The full model
therefore includes around 1600 surface and 1100 bulk-ice reaction processes
for each new generic mechanism included in the model (excluding 3-BEF, see
Section 2.5). All the grain-surface and bulk-ice reactions allowed in the
network are presented in machine-readable format in Table 12. The rate
formulations for diffusive and non-diffusive reaction mechanisms used in the
model are described below.

###  2.1 New chemical mechanisms

The standard formulation, as per e.g. Hasegawa et al. (1992), for the
treatment of a diffusive grain-surface chemical reaction (also known as the
Langmuir-Hinshelwood or L-H mechanism) between species $A$ and $B$ is based
on: the hopping rates of the two reactants, $k_{\mathrm{hop}}(A)$ and
$k_{\mathrm{hop}}(B)$; the abundances of both species on the grains, $N(A)$
and $N(B)$, which here are expressed as the average number of atoms or
molecules of that species present on an individual grain; the total number of
binding sites on the grain surface, $N_{S}$, often assumed to be on the order
of 1 million for canonically-sized grains; and an efficiency related to the
activation energy barrier (if any), $f_{\mathrm{act}}(AB)$, which takes a
value between zero and unity. Thus, the total rate of production (s-1) may be
expressed in the following form (which is arranged in such a way as to
demonstrate its correspondence with the non-diffusive mechanisms discussed
later):

$R_{AB}=f_{\mathrm{act}}(AB)\,\left[k_{\mathrm{hop}}(A)\,N(A)\right]\frac{N(B)}{N_{S}}\,+f_{\mathrm{act}}(AB)\,\left[k_{\mathrm{hop}}(B)\,N(B)\right]\frac{N(A)}{N_{S}}$
(1)

In the first term, the expression within square brackets corresponds to the
total rate at which particles of species $A$ may hop to an adjacent surface
binding site. The ratio $N(B)/N_{S}$ gives the probability for each such hop
to result in a meeting with a particle of species $B$. Multiplying these by
the reaction efficiency gives the reaction rate associated solely with the
diffusion of species $A$. The reaction rate associated with diffusion solely
of species $B$ is given by the second term. The total reaction rate is
commonly expressed more succinctly thus:

$\displaystyle R_{AB}=k_{AB}\,N(A)\,N(B)$ (2) $\displaystyle
k_{AB}=f_{\mathrm{act}}(AB)\,(k_{\mathrm{hop}}(A)+k_{\mathrm{hop}}(B))/N_{S}$
(3)

which provides a more standard-looking second-order reaction rate. The rate
coefficient, $k_{AB}$, may be further adjusted to take account of random walk,
in which a reactant that has yet to meet a reaction partner may re-visit
previous, unsuccessful sites. This effect typically reduces the overall
reaction rate by no more than a factor of a few (e.g. Charnley, 2005; Lohmar &
Krug, 2006; Lohmar et al., 2009; Willis & Garrod, 2017).

The individual hopping rate for some species $i$ is assumed in this model to
be a purely thermal mechanism, given by

$k_{\mathrm{hop}}(i)=\nu(i)\,\exp{\left(\frac{-E_{\mathrm{dif}}(i)}{T}\right)}$
(4)

where $\nu(i)$ is the characteristic vibrational frequency of species $i$ and
$E_{\mathrm{dif}}(i)$ is the barrier against diffusion in units of K.

The reaction efficiency factor for a reaction between species $A$ and $B$
considers the case where, if there is an activation energy barrier, the
diffusion of either species away from the other may compete with the reaction
process itself (see e.g. Garrod & Pauly, 2011), thus:

$f_{\mathrm{act}}(AB)=\frac{\nu_{AB}\,\kappa_{AB}}{\nu_{AB}\,\kappa_{AB}+k_{\mathrm{hop}}(A)+k_{\mathrm{hop}}(B)}$
(5)

where $\nu_{AB}$ is taken as the faster of either $\nu(A)$ or $\nu(B)$, and
$\kappa_{AB}$ is a Boltzmann factor or tunneling efficiency for the reaction
(see Hasegawa et al., 1992). The denominator represents the total rate at
which an event may occur when species A and B are in a position to react.

In order to formulate rates for non-diffusive reaction processes of whatever
kind, the total rate must again be decomposed into its constituent parts,
which, unlike in Eqs. (2) and (3), cannot generally be re-combined.

The generic form that we adopt for such processes is:

$R_{AB}=f_{\mathrm{act}}(AB)\,R_{\mathrm{comp}}(A)\,\frac{N(B)}{N_{S}}\,+f_{\mathrm{act}}(AB)\,R_{\mathrm{comp}}(B)\,\frac{N(A)}{N_{S}}$
(6)

where $R_{\mathrm{comp}}(i)$ is labeled the “completion rate” for the
reaction, corresponding specifically to the “appearance” of species $i$. The
determination of $R_{\mathrm{comp}}(i)$ values is explained in more detail in
the following subsections for each of the specific reaction mechanisms
considered. The above form is essentially the same as that given by Garrod
(2019) for photodissociation-induced reactions. The correspondence of Eq. (1)
(for diffusive reactions) with Eq. (6) is clear; the latter may be considered
a more general description of a surface reaction rate, which can be applied to
both diffusive and non-diffusive processes, according to the chosen form of
$R_{\mathrm{comp}}(i)$. The regular diffusive mechanism would use
$R_{\mathrm{comp}}(i)=k_{\mathrm{hop}}(i)\,N(i)$.

While the general form given in Eq. (6) is set up to describe grain-surface
processes, it may easily be adapted for bulk-ice processes by substituting
$N_{S}$ for $N_{M}$, the total number of particles in the ice mantle, with
$N(i)$ now representing the number of atoms/molecules of species $i$ present
in the mantle. In this case, mantle-specific diffusion rates should be used.

The several non-diffusive processes incorporated into the chemical model,
based on Eq. (6), are described below. Table 2 indicates which specific new
mechanisms are included in each of the model setups tested in the present
study.

###  2.2 Eley-Rideal reactions

The Eley-Rideal (E-R) reaction process occurs when some atom or molecule that
is adsorbing/accreting from the gas phase onto the grain surface immediately
encounters its grain-surface reaction partner as it adsorbs. Ruaud et al.
(2015) considered a more intricate treatment than we use here, in which an
adsorbing carbon atom could enter into a bound complex with a surface
molecule. Here we adopt a more generalized treatment in which we do not
differentiate between the binding properties of local surface species.

The rates for the Eley-Rideal process can easily be represented using Eq. (6).
For reactions that have no activation-energy barrier (and for which
$f_{\mathrm{act}}(AB)$ is therefore close to unity), this is achieved by
setting $R_{\mathrm{comp}}(i)=R_{\mathrm{acc}}(i)$, the total accretion rate
of species $i$ from the gas phase.

In the interests of completeness, it is necessary also to consider how to
treat the kinetics of E-R reactions that have at least some modest activation-
energy barrier. To this end, one may consider a hypothetical case where oxygen
atoms are slowly accreting onto an otherwise pure CO surface. For purposes of
illustration, it is initially assumed here that the surface-diffusion rates of
both O and CO are negligible, with the result that
$f_{\mathrm{act}}(\mathrm{O+CO})$=1.

The reaction O + CO $\rightarrow$ CO2 has an activation energy on the order of
1000 K; although the reaction will not be instant, in the absence of all other
competing processes it should nevertheless occur on some finite timescale.
Thus, the total timescale for the complete E-R process for an individual
accreting O atom encountering and reacting with a surface CO molecule would be
the sum of: (i) the accretion timescale of the oxygen atom onto the surface
and (ii) the lifetime against its subsequent reaction with CO, i.e.
$1/R_{\mathrm{acc}}(\mathrm{O})+1/(\nu_{\mathrm{O+CO}}\,\kappa_{\mathrm{O+CO}})$.
This would provide a total completion rate associated with O accretion (to be
employed in Eq. 6) of:

$R_{\mathrm{comp}}(\mathrm{O})=\frac{1}{1/R_{\mathrm{acc}}(\mathrm{O})+1/(\nu_{\mathrm{O+CO}}\,\kappa_{\mathrm{O+CO}})}$

The completion rate $R_{\mathrm{comp}}(\mathrm{O})$ should be viewed as the
rate at which the reaction process occurs successfully from the point of view
of an individual accreting O atom, taking into account all sequential steps in
the completion of the reaction process. Note that, in the full description,
the probability of encountering a CO molecule on the surface,
$N(\mathrm{CO})/N_{S}$, and the reaction efficiency,
$f_{\mathrm{act}}(\mathrm{O+CO})$, should both remain outside of the formula
for $R_{\mathrm{comp}}(\mathrm{O})$, as per Eq. (6). Neither of these values
affects the actual timescale over which an individual O atom successfully
accretes and reacts with a surface CO molecule; rather, they affect the
probability that a single such event is successful.

This expression for $R_{\mathrm{comp}}(\mathrm{O})$ could result in one of two
important outcomes, depending on the relative rates of accretion and reaction.
If accretion of O is very slow, and therefore reaction is comparatively fast,
then $R_{\mathrm{comp}}(\mathrm{O})\simeq R_{\mathrm{acc}}(\mathrm{O})$. Since
$N($CO$)/N_{S}\simeq 1$, this means that the total Eley-Rideal production rate
would initially be $R_{\mathrm{O+CO}}\simeq R_{\mathrm{acc}}(\mathrm{O})$. In
other words, the overall production rate of CO2 is only limited by the rate of
O accretion onto the surface, which is as one would expect for this case.

However, if reaction is slower than or comparable to the initiating accretion
process, each accretion of O would be followed by some significant lag-time
between accretion and reaction, which must be accounted for in the overall
rate; the incorporation of the above expression for
$R_{\mathrm{comp}}(\mathrm{O})$ into Eq. (6) indeed does this. Without this
expression, and instead using the value
$R_{\mathrm{comp}}(\mathrm{O})=R_{\mathrm{acc}}(\mathrm{O})$, the rate of
conversion of O and CO into CO2 would incorrectly be set to the accretion rate
of O. The correct formulation gives a total E-R reaction rate that is less
than the total accretion rate, allowing the build-up of O on the surface.

The final adjustment to the barrier-mediated E-R treatment comes into play
when one or other surface diffusion rate is non-negligible. If diffusion of
(say) O is indeed fast compared to reaction, then the reaction efficiency,
$f_{\mathrm{act}}(\mathrm{O+CO})$, becomes small, which reduces the total rate
of the reaction as per Eq. (6). However, the completion rate
$R_{\mathrm{comp}}(\mathrm{O})$ must also be adjusted, to correspond only to
the instances in which the O+CO reaction is actually successful. Successful
reactions would have to occur before the diffusive separation of the two
reactants could render the process unsuccessful, so the reaction timescale
would become shorter, even though the reaction probability (i.e.
$f_{\mathrm{act}}$) were reduced. For this reason, diffusion rates must also
be considered when formulating $R_{\mathrm{comp}}(\mathrm{O})$. Using a more
general description for reactants $A$ and $B$, the average lifetime against
some event occurring (including reaction itself), once the reactants are in a
position to react, may be described more fully by the expression:

$t_{\mathrm{AB}}=1/(\nu_{AB}\,\kappa_{AB}+k_{\mathrm{hop}}(A)+k_{\mathrm{hop}}(B))$
(7)

which can then be used in the general definitions:

$\displaystyle R_{\mathrm{comp}}(A)=\frac{1}{1/R_{\mathrm{app}}(A)+t_{AB}}$
(8) $\displaystyle
R_{\mathrm{comp}}(B)=\frac{1}{1/R_{\mathrm{app}}(B)+t_{AB}}$ (9)

where $R_{\mathrm{app}}(i)$ is the “appearance rate” of species $i$, which in
the case of the E-R mechanism is simply $R_{\mathrm{acc}}(i)$.

It should be noted that once diffusion becomes significant, a model even as
simple as the one used above to describe pure Eley-Rideal reaction processes
would be incomplete; the standard diffusive reactions described by Eq. (1)
must also be considered (as an entirely separate process) in such a model, to
handle the occasions where accreting atoms (e.g. O) do not immediately react
with their reaction partners (e.g. CO) before they diffuse away to another
binding site, where they may also have the ability to react. In this case, the
Eley-Rideal expressions would depend much less strongly on the time-lag effect
described above, meaning that $R_{\mathrm{comp}}$ and $R_{\mathrm{acc}}$ would
be similar in cases where diffusion of either reactant were relatively fast.
In practical application to astrochemical models, for non-diffusive reactions
whose reactants have slow or negligible diffusion rates, other processes could
also act to interfere with the reaction; for example, the UV-induced
dissociation of one or other reactant might occur on a shorter timescale than
a very slow reaction, or a hydrogen atom might arrive to react with one of
other reactant, before the reaction in question could occur. Competition from
processes such as these would prevent very slow reactions (i.e. those with
large activation barriers) from becoming important, even where diffusion of
the reactants were negligible. A yet more complete treatment of reaction
competition would include rates for these processes in Eqs. (5) and (7).

In our chemical model MAGICKAL, Eqs. (6)–(9) are used to set up Eley-Rideal
versions of all allowed grain-surface reactions in the network. Because the
E-R process is exclusively a surface process, no such processes in the ice
mantle are included. Note that, when incorporating the Eley-Rideal mechanism
into a model, no modification of the accretion (adsorption) rates themselves
is required, since the Eley-Rideal mechanism does not replace any part of the
adsorption rate. Rather, the E-R mechanism occurs immediately after
adsorption, and therefore acts as a sink on the surface populations of the
reactants, even though its rate is driven by the rate of arrival from the gas
phase of one or other reactant.

Equations with the general form of Eqs. (6)–(9) are used also to formularize
the remaining non-diffusive reaction mechanisms described below, where
$R_{\mathrm{app}}(i)$ is the only quantity to vary between processes. These
formulations can be used equally well for processes with or without
activation-energy barriers.

While Eq. (6) is still valid for the regular diffusive reaction mechanism with
completion rates of $R_{\mathrm{comp}}(i)=k_{\mathrm{hop}}(i)\,N(i)$, no
adjustment following Eqs. (7)–(9) should be used, nor would be needed. The
formulation required to model any lag-time for diffusive reactions is
different from that of non-diffusive processes (because $N(i)/N_{S}$ and
$f_{\mathrm{act}}(AB)$ cannot remain outside the $R_{\mathrm{comp}}(i)$
expression), but there are no circumstances in which such a lag-time would be
significant.

###  2.3 Photodissociation-induced reactions

Garrod (2019) suggested that the omission of non-diffusive, photodissociation-
induced reactions from models of interstellar ice chemistry may result in the
photolytic production of COMs being severely underestimated. Past models of
chemistry in star-forming regions (e.g. Garrod et al., 2008, 2017) have
allowed photodissociation to contribute to the production of COMs in the
surface and bulk-ice phases in only an indirect way, mediated by thermal
diffusion. That is, photodissociation of various molecules produces radicals,
which are separately allowed to react through the standard diffusive
mechanism. Thus, at very low temperatures, no significant COM production is
seen via radical-radical recombination, as diffusion of radicals is minimal.
However, the presence of radicals in or upon the ice means that in some
fraction of photodissociation events, the products may sometimes be formed
with other reactive radicals already present nearby. In this case, the
immediate products of photodissociation could react with the pre-existing
radicals either without diffusion, or following some short-ranged, non-thermal
diffusion process (possibly enabled by the excitation of the dissociation
products).

Eqs. (6)–(9) can again be used to describe this process, with an appropriate
choice for $R_{\mathrm{app}}(i)$, which is simply the total rate of production
of photo-product $i$ caused by all possible photodissociation processes:

$R_{\mathrm{app}}(i)=\sum_{\mathrm{all}\,j}\,R_{j}(i)\\\ $ (10)

where $R_{j}(i)$ is the production rate of $i$ via an individual
photodissociation process $j$. For the radical CH3, for example, this would
include the photodissociation of CH3OH, CH4, and various larger molecules
containing a methyl group.

If one were to consider, for example, the production of dimethyl ether through
this mechanism, an important reaction would be CH3 \+ CH3O $\rightarrow$
CH3OCH3, which is usually assumed to be barrierless. For this reaction, in Eq.
(6), species $A$ = CH3 and species $B$ = CH3O; the appearance rate of CH3
would be as described above. The main contribution to the appearance rate of
CH3O would likely be the photodissociation process CH3OH + h$\nu$
$\rightarrow$ CH3O + H. The formulation used for the dimethyl ether-producing
reaction simply states that some fraction of CH3 produced by photodissociation
of various molecules in the ice will immediately meet a CH3O radical that it
can react with, and vice versa.

Reactions affected by this PD-induced mechanism need not only be radical-
radical recombination reactions; the production, via photodissociation, of
atomic H in close proximity to CO, for example, could enhance the rate of the
reaction H + CO $\rightarrow$ HCO, which has an activation-energy barrier. The
treatment of barrier-mediated reactions in the generic Eqs. (7)–(9) is used
again for this purpose.

This treatment does not take into account any explicit consideration of
excitation of the photo-products, which could also enhance reaction rates (as
per e.g. Shingledecker et al., 2018, in the case of cosmic-ray induced
dissociation). It is also implicitly assumed that the rates of
photodissociation used in the network represent the rates at which
dissociation occurs without immediate recombination of those same photo-
products.

It is trivial to adapt the equations used for surface reactions to deal
instead with ice-mantle related processes, and this is indeed implemented in
the simulations presented here.

###  2.4 Three-body reactions

The laboratory results of authors such as Fedoseev et al. (2015, 2017), in
which H, CO and/or other species are deposited onto a cold surface, indicate
that surface reactions between radicals of low mobility may produce COMs, even
at low temperatures and without any energetic processing. The suggested
explanation is that pairs of radical species such as HCO may, on occasion, be
formed in close proximity to each other, allowing them to react either
immediately or after a very small number of surface hops. The HCO radicals
themselves would initially be formed through a more typical diffusive
(Langmuir-Hinshelwood) process or through an Eley-Rideal process, via the
barrier-mediated reaction of H and CO. Fedoseev et al. (2015) suggest that
reactions of HCO with CO may also be active, which would require no diffusion
of HCO at all, if the HCO itself is formed through the reaction of atomic H
and CO on top of a CO-rich surface.

In a similar vein, Garrod & Pauly (2011) found, using chemical kinetics
modeling, that the interstellar abundance of solid-phase CO2 could be
explained by the reaction H + O $\rightarrow$ OH occurring on a CO-rich dust-
grain ice surface. This allows the newly-produced OH to react rapidly with the
CO without any intervening thermal diffusion. They introduced into their
models a new reaction rate specifically for this process that was functionally
similar to Eq. (6).

Here, we use Eqs. (6)–(9) to calculate rates for what may be termed three-body
reactions, which include the above examples. This approach is extended to all
grain-surface reactions in the network, with a similar treatment for bulk-ice
processes. To do this, another dedicated expression for the appearance rate
$R_{\mathrm{app}}(i)$ to be used in Eqs. (8) and (9) must be constructed
specifically for three-body reactions. Eq. (10) can again be used, this time
where $R_{j}(i)$ is the production rate of $i$ (as determined using Eq. 1)
resulting from any diffusive (Langmuir-Hinshelwood) reaction, or for any non-
diffusive Eley-Rideal or photodissociation-induced reaction, whose rates are
described above. Thus, $R_{\mathrm{app}}(i)$ includes the production rates of
$i$ for all reactive mechanisms $j$ that could lead to a subsequent reaction.
From a technical point of view, the rates of all such reaction processes must
therefore be calculated in advance of the calculations for any three-body
reactions.

Using the example of the process considered by Garrod & Pauly (2011), the
reaction under consideration as a three-body process would be OH + CO
$\rightarrow$ CO2; thus, in Eq. (6), $A$ = OH and $B$ = CO. There are several
reactions in our network that could produce OH, but the main one is indeed H +
O $\rightarrow$ OH. The sum of the production rates of OH from all of these
reactions would comprise $R_{\mathrm{app}}($OH$)$. The appearance rate for CO
would also be constructed from the CO production rates of all reactions
leading to its formation.

In this way, CO2 could be formed via a three-body reaction process in which,
for example, an H and an O atom diffuse on a surface until they happen to meet
in a binding site where CO is in close proximity, they react to form OH in the
presence of the CO, and the OH and CO then subsequently react with no further
diffusion required. Alternatively, an oxygen atom could be situated in contact
with a CO molecule when an incoming H atom from the gas phase initiates an
Eley-Rideal process, leading to OH formation, followed by reaction with the
CO. Or, a CH4 molecule in close proximity to a CO molecule and an O atom could
be dissociated to H and CH3, with the H quickly reacting with the O atom to
produce OH, which would then react with CO. The prescription above would allow
many such scenarios to be included in the overall production rate of CO2,
including others relating to the formation of a CO molecule in close proximity
to an OH radical. The adoption of this generalized process means that the
special-case prescription for the OH + CO reaction introduced by Garrod &
Pauly (2011) is no longer required.

In the kind of chemical system considered by Fedoseev et al. (2015, 2017), in
which H and CO are deposited onto a surface, complex molecules could be built
up via three-body reactions between HCO radicals, initiated either by E-R or
L-H production of HCO.

Note that the new treatment does not explicitly differentiate between the case
where the newly-formed reactant is immediately in contact with the next
reaction partner, and the case where it has sufficient excess energy to allow
it to undergo a thermal hop in order to find its next reaction partner. It is
in fact highly probable that the products of exothermic reactions (which
includes virtually every surface reaction included in the network) would have
sufficient energy to allow some degree of non-thermal hopping immediately
following formation. The possibility of such energy also allowing barrier-
mediated three-body reactions to occur more rapidly is considered in the next
subsection.

To go yet a stage further, one may imagine a scenario in which the products of
three-body reactions themselves could also be involved in subsequent non-
diffusive three-body reactions. This possibility is also included in our
model, using the same equations as before, with appearance rates defined by:

$R_{\mathrm{app}}(i)=\sum_{\mathrm{all}\,j}\,R_{j,\mathrm{3B}}(i)\\\ $ (11)

where $R_{j,\mathrm{3B}}(i)$ is the production rate of $i$ caused by the
three-body reaction labelled $j$. Although these appearance rates will usually
be lower than those used in the first round of three-body reactions, the
second three-body reaction could be the most important for certain species if
they have no more dominant production mechanism. In the present models, we
allow a total of three rounds of three-body reactions to take place. Although
this could in theory be increased to any arbitrary number of rounds, the
influence of those processes rapidly diminishes beyond the second round.

As with the photodissociation-induced reactions, a similar method is employed
also for reactions in the bulk ice. In this case, the appearance rates of
reactants in the first round of three-body reactions would generally all be
products of photodissociation-induced reactions, as the Eley-Rideal process is
exclusively a surface mechanism, while the thermal diffusion of all species in
the bulk – excluding arguably H and H2 – would be very slow at the
temperatures considered in the simulations presented here.

Finally, we note that a method for treating what we label three-body reactions
was recently employed by Dulieu et al. (2019), for a single reaction between
H2CO and newly-formed H2NO. Those authors constructed a separate chemical
species to represent H2NO that is formed in contact with H2CO on a surface,
then included a special reaction in their network that occurs at a rate equal
to the vibrational collision frequency of the two contiguous species. Although
apparently different from our approach, such a treatment should also provide
the correct result; this is because, assuming that there are no competing
processes with the reaction in question, the abundance of the special H2NO
species is determined solely by its formation rate and its one destruction
rate. In that case, both the specific abundance of the special H2NO species
and its reaction rate coefficient are cancelled out in the overall rate
calculation, giving a total production rate for the reaction that is equal to
the rate that we employ in the present work. In such a case, therefore, the
chosen rate coefficient becomes immaterial to the result. It is presumably
possible to set up a large network of such reactions for newly-formed species;
however, the requirement to include a new chemical species for each reactant
pair would likely make this method prohibitive for large networks.

####  2.4.1 Specific reactions

Although the full model includes a range of three-body (3-B) processes capable
of producing acetaldehyde (CH3CHO), methyl formate (CH3OCHO) and dimethyl
ether (CH3OCH3), the dominant mechanisms for each (based on model results) are
presented below.

For acetaldehyde, the most important three-body mechanism is made up of a pair
of sequential two-body processes as follows:

$\displaystyle\rm H~{}+~{}CH_{2}$ $\displaystyle\rightarrow$ $\displaystyle\rm
CH_{3}$ (12) $\displaystyle\rm CH_{3}~{}+~{}HCO$ $\displaystyle\rightarrow$
$\displaystyle\rm CH_{3}CHO$ (13)

The most important sequential mechanisms for the other two COMs are:

$\displaystyle\rm H~{}+~{}CH_{2}$ $\displaystyle\rightarrow$ $\displaystyle\rm
CH_{3}$ (14) $\displaystyle\rm CH_{3}~{}+~{}CH_{3}O$
$\displaystyle\rightarrow$ $\displaystyle\rm CH_{3}OCH_{3}$ (15)
$\displaystyle\rm H~{}+~{}CO$ $\displaystyle\rightarrow$ $\displaystyle\rm
HCO$ (16) $\displaystyle\rm HCO~{}+~{}CH_{3}O$ $\displaystyle\rightarrow$
$\displaystyle\rm CH_{3}OCHO$ (17)

Each of these reaction pairs involves the addition of radicals in the second
step, and two of them involve the addition of atomic H to a radical in the
first step. The production of the COMs through these mechanisms should
therefore have a strong dependence on the instantaneous abundances of short-
lived reactive radicals.

The full network used in the models includes three-body versions of all the
reactions used for regular diffusive chemistry, for all surface and mantle
species.

###  2.5 Excited three-body reactions

Figure 1: Illustration of alternative mechanisms for acetaldehyde formation:
(a) the regular diffusive grain-surface reaction between radicals CH3 and HCO,
and (b) a postulated three-body excited-formation mechanism involving H, CH2
and CO, followed by a regular diffusive reaction between the radical product,
CH3CO, and another H atom. In case (a), reaction is slow at low temperatures.
In case (b), a rapid initiating reaction between H and CH2, with exothermicity
$\Delta H_{f}^{0}$ = 4.80 eV, provides enough energy to the product CH3 that
it may immediately overcome the barrier to reaction with a neighboring CO
molecule. This produces a precursor to acetaldehyde, CH3CO, that can easily be
hydrogenated by a mobile H atom to form CH3CHO.

Besides the three-body reaction process described in §2.4, we also consider a
mechanism whereby the initiating reaction produces a product that is
sufficiently excited that it is able to overcome the activation energy barrier
to a subsequent reaction. This is of particular interest if it may allow a
reaction with either CO or H2CO – both abundant surface species – that would
result in the production of a precursor to an important O-bearing COM. In this
picture, the energy of formation released by a reaction is held in the
vibrational excitation of the product species. That excited species can then
immediately react with a contiguous reaction partner.

Figure 1 shows the formation of CH3CHO via this three-body excited formation
(3-BEF) mechanism as an example. The original reaction network included the
direct association of $\ce{CH3}$ and HCO, mediated by radical diffusion, as
the main formation process for surface CH3CHO. The chance to form the CH3CHO
in that purely diffusive model is small, because it would require immobile
heavy radicals to meet at low temperature. The new 3-body process described
above, as well as the photodissociation-induced and Eley-Rideal processes,
would allow this reaction to occur non-diffusively. However, the excited
production of CH3 could also allow reaction with abundant surface CO. In the
first step, an H atom meets and then reacts with a $\ce{CH2}$ radical that is
adjacent to a CO molecule. This reaction is exothermic by 4.80 eV
($\sim$55,700 K), sufficient to overcome the barrier to the CH3 \+ CO reaction
(nominally 2,870 K, see below). Once this follow-up reaction has occurred, the
product CH3CO, which is a precursor to acetaldehyde, can easily be converted
into a stable species via hydrogenation by another H atom. The entire process
is described as follows:

$\displaystyle\rm H~{}+~{}CH_{2}$ $\displaystyle\rightarrow$ $\displaystyle\rm
CH_{3}^{*}~{}$ $\displaystyle\rm CH_{3}^{*}~{}+~{}CO$
$\displaystyle\xrightarrow{\text{E${}_{\mathrm{A}}$}}$ $\displaystyle\rm
CH_{3}CO~{}$ (18) $\displaystyle\rm H~{}+~{}CH_{3}CO$
$\displaystyle\rightarrow$ $\displaystyle\rm CH_{3}CHO~{}$

where an asterisk indicates an excited species. Similar reactions for the
production of $\ce{CH3OCH3}$ and $\ce{CH3OCHO}$ through the 3-BEF process are
as follows:

$\displaystyle\rm H~{}+~{}CH_{2}$ $\displaystyle\rightarrow$ $\displaystyle\rm
CH_{3}^{*}$ $\displaystyle\rm CH_{3}^{*}~{}+~{}H_{2}CO$
$\displaystyle\xrightarrow{\text{E${}_{\mathrm{A}}$}}$ $\displaystyle\rm
CH_{3}OCH_{2}$ (19) $\displaystyle\rm H~{}+~{}CH_{3}OCH_{2}$
$\displaystyle\rightarrow$ $\displaystyle\rm CH_{3}OCH_{3}$

$\displaystyle\rm H~{}+~{}H_{2}CO$ $\displaystyle\rightarrow$
$\displaystyle\rm CH_{3}O^{*}$ $\displaystyle\rm CH_{3}O^{*}~{}+~{}CO$
$\displaystyle\xrightarrow{\text{E${}_{\mathrm{A}}$}}$ $\displaystyle\rm
CH_{3}OCO$ (20) $\displaystyle\rm H~{}+~{}CH_{3}OCO$
$\displaystyle\rightarrow$ $\displaystyle\rm CH_{3}OCHO$

The 3-BEF process technically concerns only the first two reactions out of the
three, in each case; the final hydrogen-addition step most typically occurs
through the usual Langmuir-Hinshelwood mechanism that is already included in
the model, although non-diffusive mechanisms may also act to add the final H
atom.

Due to the more complicated requirement to consider the energy of formation in
each case, the three new 3-BEF processes shown above were individually coded
into the model, rather than constructing a generic mechanism. For this reason,
the 3-BEF mechanism is included only in the first round of three-body
processes. The production rate of the standard diffusive process for the
initiating reaction in each case is responsible for the entire value of
$R_{\mathrm{app}}(i)$, and only one term is required in Eq. (6). Crucially,
the reaction efficiency for the second reaction in the process (i.e. the
reaction whose rate is actually being calculated with the 3-BEF method) is
initially set to unity, to signify that the activation energy barrier is
immediately overcome.

Unfortunately, the activation energies of the above reactions between the
radicals and CO or H2CO are not well constrained. The chemical network of
Garrod (2013) included the CH3 \+ CO and CH3O + CO reactions, adopting a
generic activation energy barrier of 2,500 K, based on the approximate value
for the equivalent reactions of atomic H with CO and H2CO. A reaction between
CH3 and H2CO was also present in that network, with products CH4 and HCO, and
$E_{A}$=4440 K; this reaction is retained here in addition to the new pathway.

For the present network, we calculate an approximate activation energy of
2,870 K for the CH3 \+ CO reaction using the Evans-Polanyi (E-P) relation
(e.g. Dean and Bozzelli, 2000); this would be well below the energy produced
by the initiating reaction ($\sim$55,700 K). Due to the lack of comparable
reactions for a similar E-P estimate for the activation energy of the CH3 \+
H2CO reaction, the same value of 2,870 K might be assumed, placing it also
comfortably less than the energy produced by the H + CH2 $\rightarrow$ CH3
reaction. Few determinations exist for the activation energy of the CH3O + CO
reaction, although an experiment places it at 3,967 K (Huynh & Violi, 2008,
for temperatures 300-2500 K). This also is less than the energy produced by
the initiating reaction, H + H2CO $\rightarrow$ CH3O ($\sim$10,200 K). In any
case, the activation energies involved in each of the three reactions
mentioned here are sufficiently large that they should be of no importance
without the inclusion of the 3-BEF mechanism to provide the energy required,
while the 3-BEF mechanism itself is assumed to go at maximum efficiency.
However, the latter assumption may not necessarily be accurate, depending on
the form of the energy released by the reaction, and whether there is any
substantial loss prior to reaction actually occurring (see §3.3).

Table 2: Non-diffusive mechanisms included in model setups  Model |  Eley-Rideal |  Photodissociation-Induced |  Three-body |  Three-body Excited Formation |  Adjusted 3-BEF efficiency for MF
---|---|---|---|---|---
Control |  |  |  |  |
E-R | $\checkmark$ |  |  |  |
PD-Induced |  | $\checkmark$ |  |  |
3-B |  |  | $\checkmark$ |  |
3-B+3-BEF |  |  | $\checkmark$ | $\checkmark$ |
3-BEF Best |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$

###  2.6 Physical conditions

MAGICKAL is a single-point model, but a spatially-dependent picture of the
chemistry of L1544 (catalog ) can be achieved by running a set of models with
different physical conditions at specific positions within the prestellar
core. Recently, Chacón-Tanarro et al. (2019) determined the parameterized
density and temperature structure of L1544 (catalog ) as follows, considering
the optical properties of dust grains as a function of radius:

$T(r)\,\mathrm{[K]}=12-\frac{12-6.9}{1+\left(\frac{r}{28.07^{\prime\prime}}\right)^{1.7}}$
(21) $n_{\mathrm{H}_{2}}(r)\,\mathrm{[cm}^{-3}\mathrm{]}=\frac{1.6\times
10^{6}}{1+\left(\frac{r}{17.3^{\prime\prime}}\right)^{2.6}}$ (22)

where $r$ is measured in arcseconds. Based on this density structure, we
determine 15 densities at which the chemical models are to be run, ranging
logarithmically from the minimum of $n_{\textrm{H}}=4.4\times 10^{4}~{}{\rm
cm^{-3}}$ ($\sim$11,000 AU) to the maximum of $n_{\textrm{H}}=3.2\times
10^{6}~{}{\rm cm^{-3}}$ (core center). An additional eight density points are
then placed to achieve better spatial resolution toward the core center (where
the density profile is relatively flat). The appropriate temperature for each
point is then chosen from the profile, based on density/radius.

In order to take account of the gradual collapse of the gas into this final
density profile, the density used for each chemical model in the set is
independently evolved using a simple modified free-fall collapse treatment.
(The radial position of each model point is thus not explicitly considered
during this evolution). Each point begins with a gas density of
$n_{\textrm{H,0}}=3\times 10^{3}~{}{\rm cm^{-3}}$, with an initial
$\ce{H}/\ce{H2}$ ratio of $5\times 10^{-4}$. The density evolution stops once
each model reaches its specified final density, resulting in a marginally
different evolutionary time for each density point. The collapse treatment is
based on that used by Rawlings et al. (1992). Magnetic fields can play an
important role in the equilibrium of dense cores, significantly slowing down
the collapse process. Estimating an accurate collapse timescale is
challenging, although the ratio of the ambipolar diffusion time
($\tau_{\textrm{ap}}$) and the freefall time ($\tau_{\textrm{ff}}$) is
typically assumed to be $\tau_{\textrm{ap}}/\tau_{\textrm{ff}}\thicksim 10$
(see, e.g., Hennebelle & Inutsuka, 2019). Thus, the magnetic retardation
factor for collapse, B, is adopted here to control the collapse timescale.
This parameter takes a value between 0 (static) and 1 (free-fall) and is
technically density-dependent. In our model, this value is set for simplicity
to 0.3 for all density points, which results in a collapse timescale
approximately 3 times longer than the free-fall timescale. A time of a little
over $3\times 10^{6}$ year is therefore required to reach the final density at
each point, although much of this time is spent under relatively low-density
conditions as the collapse gradually ramps up.

The density evolution for each model is accompanied by increasing visual
extinction, which evolves according to the expression
$A_{\textrm{V}}=A_{\textrm{V,0}}(n_{H}/n_{H,0})^{2/3}$ (Garrod & Pauly, 2011).
The initial extinction values are set such that the values at the end of the
chemical model runs correspond to the linear integration of the density
profile, converted to visual extinction using the relationship
$N_{H}=1.6\times 10^{21}A_{\textrm{V}}$. An additional background visual
extinction of 2 is added for all positions and times, under the assumption
that L1544 (catalog ) is embedded in a molecular cloud (e.g. Vasyunin et al.,
2017). In contrast to density and visual extinction, the temperature is held
steady throughout the chemical evolution for each density point, with the same
value adopted for both the gas and the dust. Temperatures range from
approximately 8 to 14 K depending on radius, which is consistent with the
observational features (Crapsi et al., 2007).

##  3 Results

Figure 2:  Time evolution of the gas-phase abundances of the three O-bearing
COMs, at the core center, for models using each of the new mechanisms. The
abundance from the control model is denoted as a black dashed line. Diamonds
indicate the abundances at the end of each model run. The black dot-dash lines
and the gray shaded regions indicate the observational abundances and their
uncertainties, obtained from Jiménez-Serra et al. (2016) for L1544 (catalog ).
The gas density is indicated by the right-hand vertical axis and the similarly
colored line.  Figure 3:  Time evolution of the abundances in the dust-grain
ice mantles of the three O-bearing COMs at the core center, for models using
each of the new mechanisms. The abundance from the control model is denoted as
a black dashed line. Diamonds indicate the abundances at the end of each model
run. The gas density is included for reference, indicated by the right-hand
vertical axis and the similarly colored line.

The time-evolution of the fractional abundances at the core-center position is
presented in Fig. 2 (gas phase) and Fig. 3 (solid phase), for each of the main
chemical model setups. In the control model, no new mechanisms are added. In
each of the other model setups, a single new mechanism is added to the control
model setup, except for model 3-B+3-BEF, in which it is assumed that the 3-BEF
mechanism could not occur without the 3-B mechanism also being active.

As seen in Fig. 2, every new mechanism introduced here, excluding E-R,
significantly increases the abundances of $\ce{CH3OCH3}$ and $\ce{CH3OCHO}$ in
the gas phase during core evolution, while $\ce{CH3CHO}$ is only substantially
increased via 3-BEF. However, it should be noted that the increased fractional
abundances rapidly drop as density increases toward the end-time of all the
models, mostly converging to the control-model values. This indicates that the
new mechanisms may hardly affect the gas-phase COMs at the core center, but
may be more effective at more distant radii (i.e. lower density regions); this
would nevertheless result in higher abundances toward the core-center position
when averaged over the line of sight to include lower-density gas.

The presence of the COMs in the gas phase following their formation on grain
surfaces is the result of chemical desorption. All surface reactions that form
a single product have a small possibility of returning that product to the gas
phase. The upper limit on the ejection probability per reaction is 1%.

Similar to the gas phase, every mechanism excluding E-R significantly
increases the solid-phase populations of the COMs (Fig. 3). Note that the
solid-phase population of $\ce{CH3CHO}$, whose gas-phase abundance is only
strongly increased by the 3-BEF mechanism, increases even in the 3-B (only)
and PD-Induced models. The 3-B and 3-B+3-BEF models converge to essentially
the same value at the end-time. Dimethyl ether in the mantle is produced in
similar quantities by each of the three effective mechanisms, while 3-BEF and
then 3-B are more important than PD-Induced formation in the case of methyl
formate. The E-R mechanism produces only marginal increases in mantle
abundances of acetaldehyde and methyl formate. The increase in dimethyl ether
production caused by E-R is around an order of magnitude throughout most of
the evolution, although this is dwarfed by the effects of the other
mechanisms.

Figure 4:  Radial distribution of the gas-phase abundances of the three
O-bearing COMs depending on the mechanisms. The black dotted line and grey
shaded region indicates the abundance from the observation and error,
respectively. The observational data and its error is referred from the
Jiménez-Serra et al. (2016). The observational error of the $\ce{CH3CHO}$ is
not provided in the reference. The abundance from the control model is denoted
as a black dashed line and a black dotted line, respectively. While the result
from the 3-B combined with the 3-BEF is denoted as a black solid line, a green
line represent the 3-B only model. The red and blue solid line indicates the
abundances from the E-R model and PD-Induced model, respectively.  Figure 5:
Radial distribution of the solid-phase abundances of the three O-bearing COMs
depending on the mechanisms. The black dotted line and grey shaded region
indicates the abundance from the observation and error, respectively. The
observational data and its error is referred from the Jiménez-Serra et al.
(2016). The observational error of the $\ce{CH3CHO}$ is not provided in the
reference. The abundance from the control model is denoted as a black dashed
line and a black dotted line, respectively. While the result from the 3-B
combined with the 3-BEF is denoted as a black solid line, a green line
represent the 3-B only model. The red and blue solid line indicates the
abundances from the E-R model and PD-Induced model, respectively.

Fig. 4 shows the radial distribution of gas-phase COM abundances using the
full radius–density–temperature profile model results; abundances shown
correspond to the end-time abundances, at which the final density profile is
achieved. The observational values that are also indicated in the figure are
for the core-center position; however, those observations correspond to a beam
of radius $\thicksim$1900 AU, and would also sample a range of physical
conditions along the line of sight – some caution should therefore be taken in
directly comparing them with the local fractional abundance values.

It may be seen that the general trend, even for the control model, is for COM
abundances to increase toward greater radii. The 3-B+3-BEF model produces
maximum molecular abundances for acetaldehyde and methyl formate similar to
the observational values. For the latter molecule, the modeled fractional
abundance exceeds the observational values at radii greater than around 2500
AU, although the absolute gas density begins to fall off at these positions,
so that they should contribute less to the total column density of the
molecule. At the largest radii modeled, the 3-B (only) model produces methyl
formate sufficient to match the observed abundance (although, again, perhaps
with little contribution to total column density). Acetaldehyde also reaches
its peak abundance at large radii, although it reaches a similar abundance at
smaller radii.

For each of the new models, the local fractional abundance of $\ce{CH3OCH3}$
is greatest at positions away from the core center, but a significant increase
in abundance is achieved at almost all positions for every model, versus the
control. However, the maximum value achieved (for the 3-B model) is still at
least two orders of magnitude lower than the observations, both in the inner
regions and at the outer edge. Curiously, for dimethyl ether, the most
effective model is the 3-B (only) model, whereas the 3-B+3-BEF model is the
most productive for the other two COMs.

In each of the 3-B, 3-B+3-BEF and PD-Induced models, acetaldehyde and dimethyl
ether abundances show a prominent peak feature at around 2000 AU. This feature
is also present for methyl formate in the PD-Induced and 3-B (only) models.
Observations of L1544 (catalog ) by Jiménez-Serra et al. (2016) show higher
fractional abundances of $\ce{CH3CHO}$ and $\ce{CH3OCH3}$ toward an off-center
position at $r\simeq$4000 AU, versus those at the core-center, with
$\ce{CH3OCHO}$ arguably showing similar behavior. Our 3-B, 3-B+3-BEF and PD-
Induced models all show this general behavior (methyl formate in the 3-B+3-BEF
model notwithstanding), albeit at a somewhat different radius from the
observations. The origin of this peak and its similarity to the observations
are discussed in more detail in § 4.2.

Figure 5 shows the radial distribution of the ice-mantle abundances at the
end-time of each model, plotted as a function of the water abundance in the
ice at each position. In contrast with the gas-phase, all the new mechanisms
but E-R significantly increase the solid-phase abundances of COMs at all
radii. This is partly because the ice mantle preserves the earlier surface
layers during the evolution of the prestellar core, when significant
enhancement of the COM abundances in the gas-phase is found (see Fig. 2),
which is itself caused by increased efficiency in the production of COMs on
the grain surfaces. However, the PD-Induced model permits substantial ongoing
processing of mantle material itself.

While COM production is not especially important in the control or E-R models,
the others attain substantial COM abundances in the ices, comparable with gas-
phase values observed in hot molecular cores. The maximum abundance achieved
by methyl formate in the 3-B+3-BEF model is close to 1% of water abundance at
the core center, i.e around $10^{-6}$ with respect to total hydrogen. This
value may thus be too high to agree with observations of hot cores/corinos, if
the abundances achieved in the prestellar stage should be preserved intact to
the later, warmer stages of evolution.

It is also noteworthy that the COM fraction in the ices is in general somewhat
greater at larger radii, although the dimethyl ether abundance is fairly
stable through the core in the 3-B and 3-B+3-BEF models, and methyl formate is
also stable throughout in the 3-B+3-BEF model.

Figure 6: Radial distribution of the ice abundances of the main ice
constituents. The abundance of CO and $\ce{CO2}$ is denoted as a black solid
and a black dashed line, respectively. A blue, green, and red solid line
represents the abundances of $\ce{CH4}$, $\ce{NH3}$, and $\ce{CH3OH}$,
respectively.

Figure 6 shows the final radial distribution of the main ice constituents as a
fraction of the local water ice abundance, for the control model and for the
3-B+3-BEF model. The absolute abundance of profile of water ice is fairly
constant across the profiles. The latter is taken as representative of all the
new models, due to their similarity, except for the E-R model, which is rather
similar to the control. Based on the absolute abundance profiles, column
densities for each ice species are derived by integrating along the line of
sight (without beam convolution); the resulting abundances with respect to
$\ce{H2O}$ ice column density are summarized in Table 3. Comparable
observational ice abundances are also shown, taken from Boogert et al. (2015),
who provided median values with respect to $\ce{H2O}$, along with the full
range of the observed abundances (from subscript to superscript value). Both
of our model setups produce a centrally-peaked distribution of $\ce{CH3OH}$
ice, while CO ice is approximately as abundant as $\ce{H2O}$ ice, especially
toward inner radii where the most extreme depletion occurs. With the new
mechanisms, a more gently-sloped distribution appears, and a better match with
observational abundances of CO and $\ce{CO2}$ is achieved. The other ice
components in the 3-B+3-BEF model are within the observational range as well.

Table 3: Abundances relative to the $\ce{H2O}$ ice column density.  Model | $\ce{CO}$ | $\ce{CO2}$ | $\ce{CH4}$ | $\ce{NH3}$ | $\ce{CH3OH}$
---|---|---|---|---|---
Observation aaBoogert et al. (2015); values correspond to low-mass YSOs. | $0.21^{0.85}_{(<0.03)}$ | $0.28^{0.50}_{0.12}$ | $0.05^{0.11}_{0.01}$ | $0.06^{0.10}_{0.03}$ | $0.06^{0.25}_{(<0.01)}$
Control | 0.83 | 0.07 | 0.19 | 0.09 | 0.07
3-B+3-BEF | 0.57 | 0.12 | 0.11 | 0.08 | 0.01

###  3.1 Column density analysis

Observational abundances may not accurately represent the true local
abundances within a source. This is because the observational intensities are
not only averaged over the line of sight, but are also affected by the
excitation characteristics of each observed species, and by the response of
the telescope beam. Considering this, it is indispensable to perform spectral
simulations for better comparison with the observations. The spectral model
used here simulates molecular lines (of COMs) that are expected to be
observable, and uses chemical abundances shown in Fig. 4 as the underlying
distribution. The 1-D chemical/physical model is treated as spherically
symmetric, so that molecular emission can be simulated along lines of sight
passing through the core at various offsets (including directly on-source),
assuming local thermodynamic equilibrium. Each line of sight passes through a
range of gas densities, temperatures and chemical abundances. The resulting
2-D simulated intensity maps for each frequency channel are then convolved
with a Gaussian telescope beam of appropriate size, dependent on frequency and
the telescope in question. (For a more detailed description of the spectral
model, see Garrod, 2013). The FWHM of the molecular lines is assumed here to
be 1 km/s, with a spectral resolution of 250 kHz, although the simulations are
quite insensitive to the precise choice of parameters.

The integrated intensities of the ensemble of molecular lines is used in a
rotational diagram analysis (Goldsmith & Langer, 1999) to obtain column
densities ($N_{\textrm{tot}}$) and rotational temperatures
($T_{\textrm{rot}}$) for each molecule. These quantities can then be compared
directly with those obtained from observations. Beam sizes were assumed to be
$\thicksim$28′′-31′′ between 79-87 GHz and $\thicksim$24′′-26′′ between 94-103
GHz, based on the size of the observing beam of the IRAM 30 m telescope. The
distance to the model prestellar core is assumed to be 140pc (Elias, 1978).

Table 4: List of targeted transitions with their spectroscopic properties

Molecule | Transition | Frequency (GHz) |  $E_{\textrm{up}}$ (K) |  $A_{\textrm{ij}}$ (s${}^{-}1$)
---|---|---|---|---
A-$\ce{CH3CHO}$ |  21,2-10,1 | 84.21976 | 4.96 | $2.4\times 10^{-6}$
E-$\ce{CH3CHO}$ |  50,5-40,4 | 95.94744 | 13.9 | $3.0\times 10^{-5}$
A-$\ce{CH3CHO}$ |  50,5-40,4 | 95.96346 | 13.8 | $3.0\times 10^{-5}$
E-$\ce{CH3CHO}$ |  51,4-41,3 | 98.86331 | 16.6 | $3.0\times 10^{-5}$
∗A-$\ce{CH3CHO}$ |  51,4-41,3 | 98.90094 | 16.5 | $3.0\times 10^{-5}$
AA-$\ce{CH3OCH3}$ |  32,1-31,2 | 84.63680 | 11.0 | $4.4\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  22,0-21,1 | 86.22872 | 8.36 | $3.5\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  22,1-21,2 | 89.70281 | 8.36 | $3.8\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  60,6-51,5 | 90.93754 | 19.0 | $5.7\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  32,2-31,3 | 91.47931 | 11.1 | $4.9\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  42,3-41,4 | 93.85964 | 14.7 | $5.6\times 10^{-6}$
AA-$\ce{CH3OCH3}$ |  52,4-51,5 | 96.85246 | 19.3 | $6.2\times 10^{-6}$
∗AA-$\ce{CH3OCH3}$ |  41,4-30,3 | 99.32600 | 10.2 | $8.8\times 10^{-6}$
A-$\ce{CH3OCHO}$ |  72,6-62,5 | 84.45475 | 19.0 | $8.0\times 10^{-6}$
A-$\ce{CH3OCHO}$ |  73,4-63,3 | 87.16129 | 22.6 | $7.8\times 10^{-6}$
A-$\ce{CH3OCHO}$ |  81,8-71,7 | 89.31664 | 20.1 | $1.0\times 10^{-5}$
A-$\ce{CH3OCHO}$ |  72,5-62,4 | 90.15647 | 19.7 | $9.8\times 10^{-6}$
A-$\ce{CH3OCHO}$ |  91,9-81,8 | 100.0805 | 24.9 | $1.5\times 10^{-5}$
A-$\ce{CH3OCHO}$ |  81,7-71,6 | 100.4907 | 22.8 | $1.5\times 10^{-5}$
A-$\ce{CH3OCHO}$ |  90,9-80,8 | 100.6834 | 24.9 | $1.5\times 10^{-5}$
∗A-$\ce{CH3OCHO}$ |  82,6-72,5 | 103.4787 | 24.6 | $1.5\times 10^{-5}$

Note. — Acetaldehyde line data from JPL catalogue based on the data set of
Bauder et al. (1976). Dimethyl ether line data from JPL catalogue based on the
data set of Lovas et al. (1979); Neustock et al. (1990). Methyl formate line
data from JPL catalogue based on the data set of Ilyushin et al. (2009);
Plummer et al. (1984). The representative molecular transition used for the
normalized convolved intensity analysis in § 3.4 is denoted with ∗.

The radiative transfer and rotational diagram analysis is performed toward the
on-source position, and toward two offset positions: (i) the peak of the COM
abundances (2000 AU), and (ii) the low-density outer-shell (9000 AU). By
considering these three positions, we may compare the modeled COM peaks with
the observational ones, and determine the dependence of the chemical reactions
on the local physical conditions in the prestellar core.

One strategy to apply this radiative transfer and rotational diagram technique
would be to simulate precisely the same molecular lines used in individual
observational datasets for L1544 (catalog ). However, since the present aim is
to determine a well-defined column density (and rotational temperature) based
on the models, with which observed column densities may be directly compared,
we instead choose a selection of lines that may plausibly be (or indeed have
been) detected toward cold sources, and which include a range of upper energy
levels. Emission lines of CH3CHO and CH3OCHO recently detected toward the cold
dark cloud B5 (Taquet et al., 2017) are chosen for this analysis. While Taquet
et al. (2017) detected a relatively large number of molecular lines for CH3CHO
and CH3OCHO, only four transitions of CH3OCH3 with a limited range of
$E_{\textrm{up}}$ (10–11 K) were detected by those authors. Our adoption of
only those lines could therefore cause substantial uncertainty in the
determination of $N_{\textrm{tot}}$(CH3OCH3). For this reason, we choose eight
bright (i.e. high Aij) AA-transitions of CH3OCH3 with $E_{\textrm{up}}$
ranging from 8–19 K, using the Splatalogue web tool
111http://www.cv.nrao.edu/php/splat. The spectroscopic data originate from the
JPL line list 222https://spec.jpl.nasa.gov (Bauder et al., 1976; Lovas et al.,
1979; Neustock et al., 1990; Ilyushin et al., 2009; Plummer et al., 1984); the
COM transitions considered in this analysis are listed in Table 4.

Table 5: The molecular column densities of the COMs at the core center

Model | $\ce{CH3CHO}$ | $\ce{CH3OCH3}$ | $\ce{CH3OCHO}$
---|---|---|---
| $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$
Observation aaJiménez-Serra et al. (2016) |  1.2$\times 10^{12}$ |  1.5$\times 10^{12}$ (4.0$\times 10^{11}$) |  4.4$\times 10^{12}$(8.0$\times 10^{12}$)
Control |  2.6$\times 10^{10}$ (4.3$\times 10^{7}$) |  1.2$\times 10^{10}$ (1.9$\times 10^{7}$) |  1.1$\times 10^{9}$ (9.9$\times 10^{6}$)
3-B |  3.8$\times 10^{10}$ (1.1$\times 10^{8}$) |  2.6$\times 10^{10}$ (1.8$\times 10^{8}$) |  4.7$\times 10^{11}$ (4.4$\times 10^{9}$)
3-B+3-BEF |  1.2$\times 10^{12}$ (4.1$\times 10^{9}$) |  2.6$\times 10^{9}$ (7.0$\times 10^{7}$) |  1.1$\times 10^{14}$ (1.3$\times 10^{12}$)
E-R |  2.7$\times 10^{10}$ (4.5$\times 10^{7}$) |  1.3$\times 10^{10}$ (1.1$\times 10^{7}$) |  1.1$\times 10^{9}$ (9.8$\times 10^{6}$)
PD-Induced |  4.0$\times 10^{10}$ (1.2$\times 10^{8}$) |  1.4$\times 10^{10}$ (6.0$\times 10^{7}$) |  1.3$\times 10^{9}$ (9.6$\times 10^{6}$)

Note. — Values in parentheses indicate observational or rotational diagram
line-fitting (model) errors.

Table 6: The molecular column densities of the COMs toward the off-center COM
peak

Model | $\ce{CH3CHO}$ | $\ce{CH3OCH3}$ | $\ce{CH3OCHO}$
---|---|---|---
| $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$
Observation aaJiménez-Serra et al. (2016) | 3.2$\times 10^{12}$ | 7.7$\times 10^{11}$ (3.2$\times 10^{11}$) | 2.3$\times 10^{12}$ (2.8$\times 10^{12}$)
Control | 3.3$\times 10^{10}$ (3.7$\times 10^{7}$) | 2.8$\times 10^{10}$ (3.6$\times 10^{7}$) | 1.4$\times 10^{9}$ (6.9$\times 10^{6}$)
3-B | 3.6$\times 10^{10}$ (4.2$\times 10^{7}$) | 4.2$\times 10^{10}$ (5.7$\times 10^{7}$) | 5.3$\times 10^{11}$ (2.1$\times 10^{9}$)
3-B+3-BEF | 1.1$\times 10^{12}$ (4.6$\times 10^{9}$) | 1.8$\times 10^{9}$ (1.2$\times 10^{7}$) | 1.3$\times 10^{14}$(1.3$\times 10^{12}$)
E-R | 3.3$\times 10^{10}$ (3.9$\times 10^{7}$) | 2.9$\times 10^{10}$ (2.9$\times 10^{7}$) | 1.4$\times 10^{9}$ (6.8$\times 10^{6}$)
PD-Induced | 3.3$\times 10^{10}$ (3.7$\times 10^{7}$) | 2.8$\times 10^{10}$ (3.6$\times 10^{7}$) | 1.5$\times 10^{9}$ (7.3$\times 10^{6}$)

Note. — The off-position for the observations is 4000 AU, corresponding to the
observational methanol peak of L1544 (catalog ) (Jiménez-Serra et al., 2016).
In the models, the fractional abundance peak occurs around 2000 AU. Values in
parentheses indicate observational or rotational diagram line-fitting (model)
errors.

Table 7: The molecular column densities of the COMs toward the outer-shell
(radius 9000 AU)

Model | $\ce{CH3CHO}$ | $\ce{CH3OCH3}$ | $\ce{CH3OCHO}$
---|---|---|---
| $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$
Control | 2.5$\times 10^{10}$ (3.8$\times 10^{6}$) | 3.9$\times 10^{10}$ (8.8$\times 10^{7}$) | 1.3$\times 10^{9}$ (1.4$\times 10^{7}$)
3-B | 3.0$\times 10^{10}$ (3.1$\times 10^{6}$) | 5.7$\times 10^{10}$ (4.4$\times 10^{7}$) | 3.7$\times 10^{11}$ (3.5$\times 10^{9}$)
3-B+3-BEF | 7.5$\times 10^{11}$ (2.0$\times 10^{8}$) | 1.5$\times 10^{9}$ (2.2$\times 10^{5}$) | 9.5$\times 10^{13}$ (6.0$\times 10^{11}$)
E-R | 2.6$\times 10^{10}$ (3.8$\times 10^{6}$) | 4.0$\times 10^{10}$ (8.7$\times 10^{7}$) | 1.3$\times 10^{9}$ (1.4$\times 10^{7}$)
PD-Induced | 2.5$\times 10^{10}$ (3.7$\times 10^{6}$) | 3.9$\times 10^{10}$ (8.8$\times 10^{7}$) | 1.5$\times 10^{9}$ (1.5$\times 10^{7}$)

Note. — Values in parentheses indicate rotational diagram line-fitting errors.

Figure 7: Column densities calculated from each model. Red, blue, and green
histogram represents the molecular column densities of AA, DME, and MF,
respectively. The errors of the modeled column densities are not presented
here as those values are relatively small ($<3\%$) compared to the bar size.
The observed value and its error bound are presented together with solid and
dashed horizontal lines respectively. For the COM peak position (lower-left
panel), the the observations was performed towards 4000 AU, corresponding to
the observational methanol peak of L1544 (catalog ) (Jiménez-Serra et al.,
2016), while the off-position for the model is 2000 AU, where fractional
abundance peak occurs.

###  3.2 Column densities of O-bearing COMs toward the core center

Tables 5 – 7 (core center, 2000 AU, 9000 AU) compare the molecular column
densities obtained from the RD analysis of different chemical models with
observational values from the literature; observational errors and rotational
diagram line-fitting error estimates are given in parentheses. Figure 7 shows
the molecular column densities for each model at three different positions as
histograms. The observed value (a solid horizontal line) and its error bounds
(dashed horizontal lines) are presented together for comparison. While every
chemical model introduced here significantly underproduces CH3OCH3, both the
3-B and 3-B+3-BEF models result in meaningful differences from the control
model for the other two COMs at the core center (Table 5). One thing to note
is that the COMs are more actively formed via 3-BEF than solely by the 3-B
mechanism. For example, while 3-B+3-BEF significantly increases the column
density of $\ce{CH3CHO}$ as well as $\ce{CH3OCHO}$, 3-B substantially
increases $\ce{CH3OCHO}$ only. Furthermore, even though both the 3-B and
3-B+3-BEF mechanisms enhance the $\ce{CH3OCHO}$ population significantly, the
increment is much higher in the 3-B+3-BEF model (3-B+3-BEF even substantially
overproduces $\ce{CH3OCHO}$ – see also § 3.3).

In either the 3-B or 3-BEF models, the key quantities through which the
production rates of COMs (on the grains or in the gas phase) may be understood
are the surface abundance of reactants, and the production (i.e. appearance)
rates of their reaction partners. The latter quantity is an explicit component
of the new expressions for non-diffusive processes, whereas in the regular L-H
formulation it does not appear. The higher formation rates in the 3-BEF model
can be explained by the fact that this mechanism involves the addition of
radicals to stable compounds (which are thus more abundant on the grain
surface) in the second step of the consecutive reaction chain (Eqs. 18-20),
while the 3-B process involves the addition of sparse radicals (Eqs. 13, 15,
and 17).

Figure 8: Time evolution of the fractional ice composition of the reactants on
the grain-surface related to the 3-B+3-BEF mechanisms forming CH3OCH3(red
lines) and CH3OCHO(black lines). Note that the abundances shown refer
specifically to species on the outer grain/ice surface and not within the ice
mantles. A fractional abundance of $\sim$ $1.3\times 10^{-12}$ corresponds to
one particle per grain.

The greater rate of formation of $\ce{CH3OCHO}$ over that of $\ce{CH3CHO}$ in
either the 3-B or 3-BEF model can also be understood in the same context. In
the 3-B model, the reactants $\ce{CH3}$ and $\ce{CH3O}$ are technically
competing with each other to form either $\ce{CH3CHO}$ or $\ce{CH3OCHO}$ by
reacting with HCO radicals on the grain surface. As seen in Fig. 8, the
fractional grain-surface abundance of $\ce{CH3O}$ (shown for the core-center
position) is much higher than that of $\ce{CH3}$. The production rate of CH3O
is also much greater than that of CH3, which is partly why its surface
abundance is higher. Similarly, in the 3-BEF model, $\ce{CH3}^{*}$ and
$\ce{CH3O}^{*}$ are competing with each other to react with CO on the grain
surface to form either $\ce{CH3CHO}$ or $\ce{CH3OCHO}$; CO is abundant, and
the appearance rates of $\ce{CH3}^{*}$ and $\ce{CH3O}^{*}$ directly determine
the formation rates of $\ce{CH3CHO}$ and $\ce{CH3OCHO}$.

Note that only a fraction of newly-formed methyl radicals can take part in the
formation of the COMs through the 3-BEF mechanism, because only the excited
methyl radicals formed via hydrogenation of $\ce{CH2}$ have this mechanism
available; abstraction of H from $\ce{CH4}$ by other H atoms is slightly
endothermic, so it should not produce $\ce{CH3}^{*}$.

Thus, the radical $\ce{CH3}$ acts as a bottleneck to the formation of the COMs
in the non-diffusive models. Also, the gradual depletion of C and related
hydrocarbons from the gas phase, while CO remains abundant, means that
production of CH3 cannot keep up with the production of CO-related radicals.
The reaction of CH3 with H or H2 from the gas phase to re-form CH4 also keeps
the average grain-surface CH3 abundance low. The production of HCO and CH3O
radicals continues to be effective as methanol builds up; while the net
direction of the CO chemistry is to convert it gradually to methanol, there
are backward reactions at every step, including H-abstraction from CH3OH, that
allow the intermediate radicals to maintain some level of surface coverage and
a sustained rate of production/appearance.

Given that the formation of $\ce{CH3OCH3}$ is related to $\ce{CH3}$ in both
the 3-B or 3-BEF models, the lower column density of $\ce{CH3OCH3}$ in those
models can be explained. CO and $\ce{H2CO}$ are competing to form either
$\ce{CH3CHO}$ or $\ce{CH3OCH3}$ by reacting with the excited methyl radicals
on the grain surface, but CO is much more abundant than $\ce{H2CO}$ (Fig. 8).
The small amount of excited methyl radicals on the surface are thus
preferentially consumed to form $\ce{CH3CHO}$.

The E-R mechanism does not make a substantial difference to the gas-phase
abundances versus the control (see Fig. 4 and Table 5). This is because the
E-R process requires high surface coverage of the reactive species on the
grains to be effective. This result is not exactly consistent with the results
from Ruaud et al. (2015). They find that the combination of E-R and their
complex-induced reaction mechanisms is efficient enough to reproduce the
observed COM abundances at temperatures as low as 10 K. Beyond the uncertainty
in the level of contribution of either mechanism, the different model
parameters of both studies should be noted: Ruaud et al. (2015) mainly focus
on the accretion of carbon atoms and assume a much higher binding energy
(3600–8400 K) than ours (800 K). This may cause higher concentration of
reactive species on the grain surface, allowing the E-R process to be
efficient.

The PD-induced reaction process is ineffective in increasing the population of
COMs in the gas phase at the core center. However, the PD-induced model
significantly increases (more than 2 orders of magnitude) the amount of COMs
in the ice mantles throughout the core’s evolution (see figure 5). Other
studies suggest indeed that the bulk ice is where the majority of physico-
chemical changes caused by radiation chemistry are likely to occur (Johnson,
1990; Spinks & Woods, 1990; Shingledecker et al., 2017). The enhanced
population of the COMs in the ice mantle does not actively affect the
population in the gas phase, because the COM products are preserved in the
mantle rather than diffusing to the grain-surface, which is directly coupled
to the gas phase. Even though this process does not make a prominent
difference in the gas-phase abundance for the prestellar core, it would
significantly affect the chemistry during the warm-up period of a protostellar
core in which accumulated mantle manterial is ejected from the grains.

###  3.3 Optimization of the 3-BEF model

As discussed in §3.2, the formation of methyl formate (CH3OCHO) through the
3-BEF mechanism is so efficient that CH3OCHO is significantly overproduced,
while this is not the case for acetaldehyde (CH3CHO). The 3-BEF mechanism as
described in §2.5 is assumed to proceed with 100 % efficiency; however, the
appropriate value in individual cases could be lower if the exothermic energy
available from the initiating reaction ($E_{\mathrm{reac}}$) is similar in
magnitude to the activation energy barrier ($E_{\mathrm{A}}$) of the
subsequent reaction. Assuming the energy is initially released into the
vibrational modes of the excited species, it may not be available in the
required mode for reaction with an adjacent species to occur before that
energy is lost to the surface, or indeed that the excited species diffuses
away entirely from its reaction partner. If the excited product has $s$
internal vibrational modes, the 3-BEF process would be expected to have
substantially sub-optimal efficiency in the case where
$E_{\mathrm{A}}>E_{\mathrm{reac}}/s$, while it would not occur at all in the
case where $E_{\mathrm{A}}>E_{\mathrm{reac}}$. The former condition would
appear to hold for the reactions shown in Eqs. (20), in which methyl formate
is produced; here, $s$(CH3O)=9, $E_{\mathrm{reac}}\simeq 10,200$ K, and
$E_{\mathrm{A}}$=3,967 K for the CH3O + CO $\rightarrow$ CH3OCO reaction
(Huynh & Violi, 2008).

Rice-Ramsperger-Kassel (RRK) theory may be introduced to obtain a statistical
estimate of the efficiency. Using the same formulation that is employed to
determine the probability of chemical desorption in the model (Garrod et al.,
2007), the probability of a successful 3-BEF process would be:

$P=\left[1-\frac{E_{\textrm{A}}}{E_{\textrm{reac}}}\right]^{s-1}$ (23)

where $s$ now includes an additional vibrational mode representing the
reaction coordinate (i.e. $s$(CH3O)=10). For the reactions forming CH3OCHO,
the values provided above give a probability of 1.2%, while for the reactions
producing CH3CHO and CH3OCH3 the probability would be 73%. This shows that the
P(CH3OCHO) of 1 originally introduced in our 3-B+3-BEF model was too high,
explaining the overproduction of that species.

For the present models, in which only three 3-BEF processes are explicitly
considered, we empirically test a selection of efficiencies for the reaction
to form CH3OCHO ranging incrementally from 100% to 0.1% in factors of 10; the
other two 3-BEF reactions are assumed to operate at maximum efficiency. It is
found that a probability of 0.1% best reproduces the molecular column
densities from the observations.

The empirically-determined optimal efficiency is clearly lower than the simple
RRK treatment above would suggest. However, the latter does not include
competition between reaction and diffusion of the excited species, which could
account for at least a factor of a few, representing several diffusion
directions. Likewise, additional translational degrees of freedom of the
excited species could be considered in Eq. (23), rather than just one reaction
coordinate. We note also that these modifications would reduce the efficiency
of the other two 3-BEF reactions considered here, perhaps bringing them closer
to around 10%. The molecular dynamics study by Fredon et al. (2017) of
reaction-induced non-thermal diffusion would indeed suggest that translational
motion would be a necessary factor to consider in a detailed treatment of the
3-BEF process (although it should be noted that those authors assumed all of
the energy to be immediately released into translational modes, rather than
distributed also into internal vibration and/or rotation.

Figure 9 shows rotational diagrams obtained from LTE radiative transfer
calculations based on the 3-BEF Best model molecular profiles, with the beam
directed toward the core center. Table 8 compares the molecular column
densities towards the core center from this model with the observational
literature values. The errors (in parentheses) for modeled column densities
are derived from the standard deviation of linear regression fitting in
rotation diagrams. The three-body mechanisms introduced here are efficient
enough to reproduce the amount of CH3OCHO and CH3CHO in the prestellar core
when an appropriate efficiency for the 3-BEF mechanism is adopted.

Figure 9:  Rotation diagrams for the three O-bearing COMs at the core center.
The black dashed lines represent the fit.  Table 8: Column densities of COMs
in the 3-BEF Best model

Species | $\ce{CH3CHO}$ | $\ce{CH3OCH3}$ | $\ce{CH3OCHO}$
---|---|---|---
| $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$ | $(\textrm{cm}^{-}2)$
Observation (Core Center) aaJiménez-Serra et al. (2016) |  1.2 $\times 10^{12}$ |  1.5 $\times 10^{12}(2.0\times 10^{11})$ |  4.4 $\times 10^{12}(4.0\times 10^{12})$
Core Center |  1.7 $\times 10^{12}(5.7\times 10^{9})$ |  3.4 $\times 10^{10}(3.3\times 10^{8})$ |  4.0 $\times 10^{12}(2.7\times 10^{10})$
2000 AU |  1.4 $\times 10^{12}(5.4\times 10^{9})$ |  4.8 $\times 10^{10}(1.1\times 10^{8})$ |  4.8 $\times 10^{12}(1.1\times 10^{10})$
9000 AU |  8.5 $\times 10^{11}(3.3\times 10^{8})$ |  5.7 $\times 10^{10}(3.0\times 10^{7})$ |  4.0 $\times 10^{12}(3.8\times 10^{10})$

Note. — Values in parentheses indicate observational or rotational diagram
line-fitting (model) errors.

Figure 10: Comparison of chemical distribution of the 3-BEF Best (solid lines)
in the gas with those of the normal 3-B+3-BEF (dashed lines). Control model
results are also shown (dotted lines).

Figure 10 compares the chemical distribution of the 3-BEF Best results (solid
lines) with those of the control (dotted lines) and the normal 3-BEF (dashed
lines) models. While the amount of CH3OCHO is significantly reduced in the
3-BEF Best model compared to the normal 3-BEF in both gas- and solid-phase,
the population of CH3OCH3 increases (by roughly an order of magnitude in the
gas); the weakening of the 3-BEF mechanism for methyl formate production
leaves more of the CH3O radical available to participate in other reactions,
including the regular 3-B mechanism (CH3 \+ CH3O) that produces dimethyl
ether. A commensurate increase is seen in the column density values.

The adjustment to the efficiency of MF production through the 3-BEF process
also reduces the solid-phase abundance of that molecule with respect to water
back to more plausible values that are in line with the maximum typical values
observed in hotter sources (i.e. around 10-8 with respect to H2). The fraction
is higher beyond around 5,000 AU, but the total ice abundance at these
positions would also be somewhat lower.

###  3.4 $\ce{CO}$ hydrogenation and $\ce{CH3OH}$ abundances

The $\ce{CH3OH}$ map of Bizzocchi et al. (2014) shows a highly asymmetric non-
uniform ring surrounding the dust peak of L1544 (catalog ). This morphology is
consistent with central depletion and preferential release of methanol in the
region where CO starts to freeze out significantly. Jiménez-Serra et al.
(2016) shows that COMs are actively formed and already present in this
methanol peak.

Figure 11:  (a) Radial distribution of fractional abundance of $\ce{CH3OH}$
for each model. The abundance value from the observations toward the core
center is denoted a black dotted line. Gas density as a function of radius is
also indicated. (b) The normalized, convolved intensity of a representative
emission line for methanol and for each of the three COMs of interest, shown
as a function of the offset of the beam from the on-source position.

The upper panel of Fig. 11 shows the radial distribution of $\ce{CH3OH}$
fractional abundance for each of the chemical models; abundances are very
similar for all models at all positions. Methanol in the gas is mainly formed
as the result of the hydrogenation of grain-surface CO all the way to CH3OH,
followed by chemical desorption. The radial distribution of the fractional
abundance of gas-phase methanol has its peak well beyond where the
observations would suggest. However, it should be noted that the gas density
in these more distant regions drops off significantly, according to the
physical profile. The location of the peak in absolute abundance would provide
a better comparison directly with observations, although the best method is to
consider the column density structure of methanol explicitly.

The lower panel of Fig. 11 shows the normalized convolved intensity of a
representative emission line of methanol, as a function of the beam offset
from the center, using the radiative transfer model already described with the
3-BEF Best model data. Since the lines are optically thin and are well
represented by an LTE treatment (see Fig. 12 and Table 9), the line intensity
profile scales well with the column density profile along each line of sight.
The modeled methanol emission shows a peak near to 4000 AU as reported in the
observations, even though this feature is not so obvious, as the slope is
quite gentle. The same treatment is shown for the other three COMs of
interest. Methyl formate shows a fairly similar distribution of emission to
that of methanol, while the other two COMs show peaks at 2000 AU as seen in
the fractional abundances. The representative molecular transition used for
this analysis is denoted with an asterisk in Tables 4 and 9.

Figure 12: The rotation diagram for $\ce{CH3OH}$ . The black dashed lines
represent the fit.  Table 9: List of methanol transitions with their
spectroscopic properties

Molecule | Transition | Frequency (GHz) |  $E_{\textrm{up}}$ (K) |  $A_{\textrm{ij}}$ (s${}^{-}1$)
---|---|---|---|---
E-$\ce{CH3OH}$ |  5-1-40 | 84.52117 | 40.4 | $2.0\times 10^{-6}$
E-$\ce{CH3OH}$ |  2-1-1-1 | 96.73936 | 12.5 | $2.6\times 10^{-6}$
E-$\ce{CH3OH}$ |  20-10 | 96.74455 | 20.1 | $3.4\times 10^{-6}$
E-$\ce{CH3OH}$ |  30-20 | 145.09375 | 27.1 | $1.2\times 10^{-5}$
E-$\ce{CH3OH}$ |  3-1-2-1 | 145.09744 | 19.5 | $1.1\times 10^{-5}$
E-$\ce{CH3OH}$ |  10-1-1 | 157.27083 | 15.4 | $2.2\times 10^{-5}$
∗E-$\ce{CH3OH}$ |  20-2-1 | 157.27602 | 20.1 | $2.2\times 10^{-5}$

Note. — Methanol line data from JPL catalogue based on the data set of Xu et
al. (2008). The representative molecular transition used for the normalized
convolved intensity analysis in § 3.4 is denoted with ∗.

The full RD analysis is performed for methanol as for the other COMs. The
seven E-transition lines of $\ce{CH3OH}$ that were detected by Taquet et al.
(2017) are chosen for this (Table 9). A single fit to all lines provides a
column density 8.6$\times 10^{13}~{}\textrm{cm}^{-2}$ at the core center. This
value is roughly consistent with the observation (2.6$\times
10^{13}~{}\textrm{cm}^{-2}$, Bizzocchi et al., 2014). The precise value, as
with those of the other COMs, will be dependent on the fidelity of the
chemical desorption treatment used here.

##  4 discussion

Of the several new non-diffusive processes tested here, the Eley-Rideal
mechanism appears to have the least effect, due largely to the low surface
coverage of reactive species. Those reactive species that might benefit from
the spontaneous arrival of a reaction partner from the gas phase always
maintain low fractional surface coverage due to their reactivity with highly
diffusive surface species e.g. atomic H. Species that do build up a large
surface coverage, like CO, tend to have large barriers to reaction, so that
incoming species are more likely to diffuse away than to react spontaneously.
The importance of the E-R process to typical surface reactions is unlikely to
be substantial under any physical conditions as long as atomic H remains
mobile.

Photodissociation-induced reactions, in which the PD process acts
spontaneously to bring a reactive radical into contact with some other
species, has no significant influence on the gas-phase abundances of complex
organics, but has a strong effect on the COM content of the ice mantles. The
basic three-body process provides substantial improvement in the gas-phase
abundances of COMs, notably methyl formate and dimethyl ether, by allowing the
products of diffusive reactions (in some fraction of cases) to find a reaction
partner themselves without requiring further diffusion. However, the excited
formation mechanism, which allows the reaction of excited, newly-formed
radicals with stable species (in spite of activation energy barriers) has the
strongest effect, and is again most important for methyl formate and
acetaldehyde. An adjustment to the efficiency of these processes, based on the
available energy from the initiating reaction, appears to provide the best
match with observational column densities of those molecules.

It is important that the process that seems to reproduce most effectively the
gas-phase abundances of the COMs (3-BEF) is one that occurs on the grain/ice
surface itself, rather than deep within the mantle, allowing chemical
desorption to return some fraction of the product to the gas phase.

The details of the various mechanisms and their implications are discussed in
more detail below.

###  4.1 H-abstraction/recombination as an amplifier of chemical desorption

The models show substantial success in reproducing observed gas-phase column
densities, through molecular production mechanisms operating on the surfaces
of the icy dust grains. Consideration should therefore be given to the
efficiency of the desorption mechanism that releases surface molecules into
the gas phase. Although photo-desorption is included in all of the models
presented here (with the explicit assumption of fragmentation of methanol as
the result of this process), the most important ejection mechanism for grain-
surface COMs is chemical desorption. In these models, this occurs with a
maximum efficiency per reaction of 1%; this efficiency is further lowered
according to the RRK-based treatment described by Garrod et al. (2007).

Thus the formation of, for example, acetaldehyde, through Eqs. 18, culminating
in the addition of an H-atom to the CH3CO radical, may sometimes produce gas-
phase CH3CHO. However, the immediate desorption following its formation is not
the only factor in ejecting those molecules into the gas. The chemical
desorption effect is considerably amplified by the abstraction of H atoms from
existing surface COMs, followed rapidly by recombination of the resulting
radical with another H atom, inducing the ejection into the gas of some
fraction of the product molecules. In the case of methanol, for instance, once
it is formed on the grain surface through the repetitive addition of H to CO,
the abstraction of H from CH3OH by another H-atom allows it to be transformed
back to its precursor ($\ce{CH3O/CH2OH}$), providing additional chances for
chemical desorption – indeed, this process of addition and abstraction was
suggested by Minissale et al. (2016) as a mechanism by which the depletion of
CO from the gas phase could be slowed and its grain-surface conversion to
methanol delayed. Similar H-abstraction/addition processes are present for
each of the larger COMs of interest in our models.

To understand how significantly this process takes part in the overall
chemical desorption scheme, four additional test models were run for
conditions appropriate to the core center, turning off the H-abstraction
reaction for each molecule (the three larger COMs plus methanol). The local
fractional abundances of COMs from each test model are compared with the
control in Table 10. When the H-abstraction reaction of a specific COM is
turned off, the gas-phase abundance of that molecule decreases by $\thicksim$1
order of magnitude. Furthermore, when H abstraction from methanol is switched
off, it reduces the fractional abundance of other COMs such as $\ce{CH3OCHO}$
and $\ce{CH3OCH3}$, whose surface production is closely related to the
$\ce{CH3O}$ radical. The abstraction of H from methanol by other H atoms in
fact strongly favors the production of the CH2OH radical; the network employed
here uses surface reaction rates for these processes calculated by F. Goumans
and S. Andersson (see Garrod, 2013) based on harmonic quantum transition state
theory. However, as per the network of Garrod (2013), the recombination of
CH2OH with H is assumed to produce either methanol or H2CO+H2 with a branching
ratio of 1:1. The production of formaldehyde in this way can then lead to
reaction with H atoms again; this forward process strongly favors production
of the CH3O radical, thus influencing the production of DME and MF.

Table 10: Local fractional abundances of COMs at the core center when the H-abstraction reaction listed is switched off in the 3-BEF Best.  | All on | $\ce{CH3OH}+\ce{H}$ | $\ce{CH3CHO}+\ce{H}$ | $\ce{CH3OCHO}+\ce{H}$ | $\ce{CH3OCH3}+\ce{H}$
---|---|---|---|---|---
$\ce{CH3OH}$ | $2.9\times 10^{-11}$ | $5.4\times 10^{-12}$ | $2.9\times 10^{-11}$ | $2.8\times 10^{-11}$ | $2.9\times 10^{-11}$
$\ce{CH3CHO}$ | $4.5\times 10^{-13}$ | $4.5\times 10^{-13}$ | $3.8\times 10^{-14}$ | $4.5\times 10^{-13}$ | $4.5\times 10^{-13}$
$\ce{CH3OCHO}$ | $1.4\times 10^{-12}$ | $5.5\times 10^{-13}$ | $1.4\times 10^{-12}$ | $7.2\times 10^{-14}$ | $1.4\times 10^{-12}$
$\ce{CH3OCH3}$ | $2.6\times 10^{-14}$ | $8.4\times 10^{-15}$ | $2.6\times 10^{-14}$ | $2.6\times 10^{-14}$ | $1.4\times 10^{-15}$

###  4.2 COM distribution and COM peaks

As seen in Fig. 10 for the 3-BEF Best model, the COMs in the gas phase have
their lowest fractional abundances at the core center, gradually increasing
toward the outer shell of the prestellar core. This general feature is
observed regardless of model type (Fig. 4). Interestingly, a local fractional
abundance peak for COMs is found at around 2000 AU, especially for CH3CHO and
CH3OCH3. This result suggests at least qualitative agreement with the
observational result of Jiménez-Serra et al. (2016); those authors performed
deep observations of the COMs toward the low-density outer shell (4000 AU) as
well as the core center of L1544 (catalog ). While they observed higher
abundances for all three COMs at the outer position, the level of enhancement
for CH3OCHO was ambiguous, due to its large observational error.

The behavior seen in the models indicates that there are two possible peak
features (or two plausible causes for observed peaks) that could become
apparent in column densities or line intensities (e.g. lower panel of Fig. 11)
as opposed to fractional abundances. The first of these relates simply to the
increased fractional abundances of COMs at large radii, combined with the
drop-off in overall gas density at the greatest extents, producing a peak in
the absolute molecular abundances that manifests in the resulting column
density or line intensity profiles. This behavior is especially apparent for
methyl formate (which does not show the bump-like feature at around 2000 AU).
This peak seems to be in reasonably good agreement with the observational peak
position; Fig. 11 indicates peak line intensities around 4000–6000 AU. A major
cause of the lack of COMs in the gas phase at small radii (in terms of
fractional abundance, e.g. Fig. 10) is that most of the gas-phase material at
those locations has already accreted onto the grains and become locked into
the ice mantles by the end-time of the models; little CO exists in the gas
phase (on the order 10-7 with respect to total H), thus grain-surface
chemistry involving CO-related products is limited. At the greatest radii,
freeze-out is incomplete and CO chemistry is still active with the accretion
of new CO. Somewhat greater gas-phase abundances of atomic H at large radii,
caused by the density slope, also encourage H-abstraction from COMs on the
surfaces, followed by recombination and chemical desorption.

The local peak at 2000 AU in the fractional abundances of acetaldehyde and
dimethyl ether occurs particularly in models that use the 3-B and 3-BEF
processes, and manifests also in the resultant column density profiles of
those molecules (lower panel of Fig. 11). The gas density at the 2000 AU
position is at least three times higher than at the outer-peak region
(4000–6000 AU), so the inner peak tends to dominate over the outer in its
contribution to column densities, for the models/molecules in which that inner
peak occurs.

What is the origin of the inner peak at 2000 AU? It is related to the freeze-
out of gas-phase material through the core. It traces a position at which the
net rate of accretion of gas-phase material onto the grains is close to zero,
caused by the high degree of depletion that has already occurred for most
major gas-phase species other than hydrogen. For example, the gas-phase
abundance of CO reaches a local minimum at this position. At radii internal
and external to the 2000 AU peak region, there is slightly more gas-phase
material remaining to be accreted onto dust grains at the end of the model
runs, thus the ice mantles still continue (slowly) to grow. This local peak in
freeze-out is due to the combined density and temperature profiles used in the
models. The adsorption rates of neutrals scale with gas density, which is
greatest at the core center, but they also scale with the square root of the
gas temperature, which is greater at larger radii. The 2000 AU position is the
point where the two profiles combine to give the largest total adsorption
rate. The position of the maximum freeze-out position is thus strongly
dependent on the density and temperature profiles. Furthermore, given a
slightly longer model run time, the freeze-out peak would likely widen, as
other positions reached a state of near-zero net accretion onto grains.

The stronger production of COMs (acetaldehyde and dimether ether) around this
2000 AU position peak is a consequence of the changing freeze-out conditions
described above. Once the net rate of freeze-out reaches zero, it indeed
undergoes a reversal in which there is a small, net rate of loss of material
from the grains. This loss is caused by the desorption of molecular hydrogen
from the grain surface, which is slowly replenished by the gradual outward
diffusion of H2 molecules embedded deep in the ice mantles. This H2-loss
process occurs throughout all the model runs, but is of little importance
until the adsorption of non-volatile species diminishes, when gas-phase
species become depleted. Once this net loss of material from the grains starts
to occur, some molecules embedded in the upper layer of the ice mantles are
“uncovered”, becoming available for surface chemical processing. Most
importantly, this includes CH4, from which an H-atom may be chemically
abstracted through several mechanisms, increasing both the production rate of
CH3 and its surface abundance. This drives up the three-body production of
acetaldehyde and dimethyl ether (Eqs. 12–15), which are chemically desorbed
into the gas phase – either directly, or as the result of H-abstraction and
recombination on the surface.

The behavior of the inner peak in certain COMs should therefore be treated
with a degree of skepticism. Not only does its position depend on the
interplay of the observationally-determined physical profiles, but its
strength must be time-dependent. Furthermore, the ability of the chemical
model to treat accurately the return (“uncovering”) of mantle material to the
ice surface is limited by the use of only a single mantle phase, rather than
the consideration of distinct layers within the ice (cf. Taquet et al., 2014).
If most of the methane residing in the mantles is present mainly in the
deepest layers, the inner-peak effect described above would be overestimated
here. It is also the case that, even with this mechanism in play, the gas-
phase abundance of dimethyl ether is insufficient to reproduce observed column
densities in L1544 (catalog ) (although see § 4.6). If such a mechanism is
active, considering the uncertainty in its precise position (based on models),
it may not be easily distinguished from the outer peak at 4000+ AU.

The peak in methanol column density occurs at the outer peak position, caused
again by a peak in absolute abundance of that molecule. It is noteworthy that
in the present models, the local fractional abundance of methanol does not
need to exceed a value of a few 10-9 to be able to reproduce the observed
column density. Again, the strength of the methanol peak will be dependent on
the efficiency of chemical desorption for that molecule, which is not well
constrained through purely experimental means.

###  4.3 The effect of diffusion barriers

Figure 13: Comparison of chemical distribution of the COMs in the 3-BEF Best
models with different values of $E_{\mathrm{dif}}$:$E_{\mathrm{des}}$ for
atoms.  Figure 14: Comparison of chemical distribution of the reactants in the
3-BEF Best models with different values of
$E_{\mathrm{dif}}$:$E_{\mathrm{des}}$ for atoms.

In many astrochemical models including MAGICKAL, chemistry on the grains is
governed by the diffusion of surface species via thermal hopping (any non-
diffusive processes notwithstanding). The energy required for a particular
species to hop from one surface binding site to another is given by the
diffusion barrier $E_{\mathrm{dif}}$; this value is parameterized in the
chemical model as some fraction of the desorption energy, i.e.
$E_{\mathrm{dif}}$/$E_{\mathrm{des}}$. Even though this is a key parameter to
describe the mobility of species on grain surfaces, the exact value has not
historically been well constrained, broadly ranging from 0.3 to 0.8. In the
present models, this parameter was set to
$E_{\mathrm{dif}}$/$E_{\mathrm{des}}$=0.6 for all atomic species, which leans
toward a high value based on recent experimental estimates by Minissale et al.
(2016), who suggested 0.55 for atoms. In our past models (e.g. Garrod 2013),
atoms and molecules were assigned the same fractional barrier value of 0.35,
based on the optimum value for CO. All molecular species in the present models
retain the 0.35 value.

At the very low surface temperatures that are found in prestellar cores, the
diffusion of atoms in particular is of great importance. For this reason, test
models were also run using the previous fractional diffusion barrier of 0.35
for atoms. Figure 13 shows a comparison of COM abundances for the two cases,
shown for the end-time of the 3-BEF Best model run using the L1544 (catalog )
physical profiles as before. Using the higher
$E_{\mathrm{dif}}$/$E_{\mathrm{des}}$ ratio, the COMs typically show much
higher abundances at positions near the core center. This result is somewhat
contradictory to the expectations of Vasyunin et al. (2017), who suggest that
the $E_{\mathrm{dif}}$/$E_{\mathrm{des}}$ ratio would not play a crucial role
in cold environments, as diffusion of H and $\ce{H2}$ via tunneling is
dominant. In our model,while tunneling through chemical barriers is included,
surface diffusion via tunneling is not, as the barriers are assumed to be too
broad for tunneling to be effective. In this case, the higher diffusion
barrier for the atomic species means that the time taken for H atoms to reach
and react with surface radicals is increased. This consequently raises the
lifetimes of those radicals on the surface (see Fig. 14), which in turn
renders the non-thermal mechanisms explored here more effective, increasing
the production of COMs.

It should be noted that the higher $E_{\mathrm{dif}}$/$E_{\mathrm{des}}$ ratio
does not always result in a larger amount of COMs (or radicals) on the grain
surface. For example, the discrepancy in the COM abundances between the two
models decreases at large radii, and methyl formate and acetaldehyde here are
even a little more abundant in the case where atomic diffusion barriers are
lower, due to slightly more effective H-abstraction from methane to produce
CH3.

The variation of the $E_{\mathrm{dif}}$/$E_{\mathrm{des}}$ ratios thus has an
important effect on the chemical model results; the higher value for atomic
species, and for H in particular, reproduces COM abundances more effectively,
through the increase in radical lifetimes. Senevirathne et al. (2017)
calculated the distribution of binding energies and diffusion barriers for H
on an amorphous water surface, suggesting representative values for each;
although their H binding energy (661 K) is higher than the value used in our
models (450 K), their diffusion barrier (243 K) is close to the value we use
here (270 K) for the $E_{\mathrm{dif}}$/$E_{\mathrm{des}}$=0.6 models. We note
also that Senevirathne et al. (2017), based on their calculations of diffusion
rates using quantum transition-state theory, aver that tunneling (as opposed
to the thermal mechanism) is likely of limited importance under most
temperature conditions in dark clouds; our use of purely thermal diffusion
rates in MAGICKAL is thus broadly consistent with that work.

In other models, in which non-diffusive chemical processes were not included,
the variation of the H diffusion barrier might have less of an effect, as most
of the active chemistry in that case would involve only atomic H. The lifetime
of the radicals would therefore be of less relevance, since H would still be
the dominant reaction partner. In the present models, the mobility of atomic
hydrogen is a major determinant of the effectiveness of non-thermal processes
in producing complex species.

###  4.4 Gas-phase processes

Due perhaps to the generally low abundances of DME that our chemical models
provide, they do not appear to reproduce the correlation between CH3OCHO and
CH3OCH3 sometimes observed in various evolutionary stages of star-forming
regions (Jørgensen et al., 2011; Brouillet et al., 2013; Jaber et al., 2014).
As a means by which such a relationship might arise, Brouillet et al. (2013)
suggested protonated methanol $\ce{CH3OH2+}$ in the gas-phase as the common
precursor to form CH3OCHO and CH3OCH3 via reactions with $\ce{HCOOH}$ and
$\ce{CH3OH}$. As a test, the proposed reactions were incorporated into our
chemical network; however, they were too slow to be effective in producing
CH3OCHO and CH3OCH3 in our model, due to the low abundance of protonated
methanol in the gas.

Recently, potentially influential gas-phase reactions were proposed by Shannon
et al. (2013, 2014), who found that reactions of either OH or
C(${}^{3}\textrm{P}$) with methanol are efficient at low temperatures, due to
quantum tunneling:

$\displaystyle\ce{OH}~{}+~{}\ce{CH3OH}~{}$ $\displaystyle\rightarrow$
$\displaystyle~{}\ce{CH3O}~{}+~{}\ce{H2O}$
$\displaystyle\ce{C}(^{3}\textrm{P})~{}+~{}\ce{CH3OH}~{}$
$\displaystyle\rightarrow$ $\displaystyle~{}\ce{CH3}~{}+~{}\ce{HCO}$

The gas-phase methanol reactions could act not only as an efficient loss
process for gas-phase methanol, but also produce more radicals that would be
available as reactants to form other COMs when they accrete onto grain
surfaces, or directly in the gas phase if such processes are efficient.
Vasyunin & Herbst (2013) suggested a gas-phase radiative association reaction
between the radicals that are produced by the above mechanisms, to form DME:

$\displaystyle\ce{CH3O}~{}+~{}\ce{CH3}$ $\displaystyle\rightarrow$
$\displaystyle\ce{CH3OCH3}+h\nu$

To understand how significantly these reactions would affect the overall
formation of COMs in our chemical model, we ran a test model that included all
three (with rate coefficients on the order of $10^{-10}$ cm-3). However, gas-
phase methanol was still predominantly destroyed by ion-molecule reactions at
the core center. The contribution of the above neutral-neutral reactions to
the loss of gas-phase methanol was minor ($\thicksim 2\%$), hardly changing
the abundances of methanol and the three COMs, while the radiative association
reaction also showed minimal influence.

Balucani et al. (2015) proposed a gas-phase mechanism that would form methyl
formate from dimethyl ether through the radical $\ce{CH3OCH2}$. The dimethyl
ether itself would form through the efficient radiative association of the
radicals CH3 and CH3O :

$\displaystyle\ce{CH3OCH3}~{}+~{}\ce{(F,Cl)}\rightarrow\ce{CH3OCH2}+~{}\ce{(HF,HCl)}$
$\displaystyle\ce{CH3OCH2}~{}+~{}\ce{O}\rightarrow\ce{HCOOCH3}~{}+~{}\ce{H}$

Although our network does not include fluorine, the incorporation of the other
reactions into our model did not make a meaningful difference to the results,
because they involve a one-way process where CH3OCH3 is converted into
CH3OCHO. In our chemical model, neither the radiative association of the
radicals CH3 and CH3O nor any other processes were efficient enough to form
abundant CH3OCH3. As such, several key reactions concerning gas-phase
chemistry of COMs do not affect our chemical model significantly.

Thus, at least under the conditions tested in our physical model, we found no
efficient gas-phase mechanisms that could produce either DME or MF.

###  4.5 Other surface processes

The 3-body excited-formation mechanism included here is especially efficient
for the initiating reaction H + CH2 $\rightarrow$ CH3, which is highly
exothermic, but which also results in a small product, CH3, that has only a
limited number of vibrational modes in which the resulting energy may be
stored. The models suggest that when this is coupled with highly abundant CO
on the grain surface, the subsequent reaction between the two proceeds at a
sufficient pace to produce enough CH3CO (and thence CH3CHO) to be able to
explain the gas-phase abundance of the latter molecule (given an adequate
desorption mechanism). The production of CH3O via the hydrogenation of
formaldehyde is also exothermic, but not sufficiently so to allow the
subsequent reaction with abundant CO to proceed at high efficiency.
Nevertheless, this low-efficiency mechanism is capable of producing enough
CH3OCO (and thence CH3OCHO) to account for the presence of methyl formate in
the gas phase.

As noted in Section 3.3, a more detailed treatment of the 3-BEF mechanism
should include not only the energy partition between bonds, but also
translational degrees of freedom of the excited species. This would impact the
RRK calculation, but could also provide an alternative outcome to the process.
The RRK treatment as formulated in Section 3.3 assumes that the efficiency of
the process is determined solely by the competition between energy going into
the “reaction mode” and energy being lost to the surface. However, if