Chemical fractionation of deuterium in the protosolar nebula

TL;DR

This study models deuterium chemistry in a protostellar envelope, revealing how ice processes influence D/H ratios and explaining observed super-deuteration in certain molecules.

Contribution

It introduces a comprehensive astrochemical model including bulk-ice chemistry to analyze deuterium fractionation in star-forming regions.

Findings

Reproduces high D/H ratios in formaldehyde and dimethyl ether.

Shows ice chemistry hampers super-deuteration in most species.

Finds D2O/HDO ratio exceeds HDO/H2O due to gas-phase processes.

Abstract

Understanding gas-grain chemistry of deuterium in star-forming objects may help to explain their history and present state. We aim to clarify how processes in ices affect the deuterium fractionation. In this regard, we investigate a Solar-mass protostellar envelope using an astrochemical rate-equation model that considers bulk-ice chem- istry. The results show a general agreement with the molecular D/H abundance ratios observed in low-mass protostars. The simultaneous processes of ice accumulation and rapid synthesis of HD on grain surfaces in the prestellar core hampers the deuteration of icy species. The observed very high D/H ratios exceeding 10 per cent, i.e., super- deuteration, are reproduced for formaldehyde and dimethyl ether, but not for other species in the protostellar envelope phase. Chemical transformations in bulk ice lower D/H ratios of icy species and do not help explaining the super-deuteration. In the protostellar phase, the D2O/HDO abundance ratio was calculated to be higher than the HDO/H2O ratio owing to gas-phase chemistry. Species that undergo evaporation from ices have high molecular D/H ratio and a high gas-phase abundance.

Tables

Table 1. Table 1: Density at the center of the cloud core, radius r 1 subscript 𝑟 1 r_{1} of the core, and core mass at the start, at starbirth time, and at the end of the simulation.

$t$, Myr	$n_H$, cm-3	$r_1$, AU	$M_{\mathrm{core}}$, $M_{\odot }$
-1.5	$2.3\times {10}^3$	$1.0\times {10}^4$	0.035
0	$1.0\times {10}^7$	$1.0\times {10}^4$	1.18
0.5	$1.0\times {10}^7$	$1.0\times {10}^3$	—

Table 2. Table 2: Initial abundances of chemical species with respect to H nuclei.

Species	Abundance
H2	0.5
HD	1.00E-05
He	9.00E-02
C	1.40E-04
N	7.50E-05
O	3.20E-04
Na	2.25E-09
Mg	1.09E-08
Si	9.74E-09
P	2.16E-10
S	9.14E-08
Cl	1.00E-09
Fe	2.74E-09

Table 3. Table 3: Summary of the energies characterizing the mobility of icy species on grains in the model.

Surface    $|$			Bulk ice (mantle sublayers)
Energy	Notation	Value	Energy	Notation	Value
Desorption (adsorption)	$E_D$	. . . 111Data from the network made available by Albertsson et al. (2013).	Absorption	$E_B$	3.0$E_D$222Kalvāns (2015a, b)
Binding (diffusion)	$E_{b,s}$	$0.35E_D$333Garrod & Pauly (2011)	Binding (diffusion)	$E_{b,m}$	0.35$E_B$

Table 4. Table 4: Observational data on deuterium fractionation ratios for gas-phase chemical species in low-mass protostars and agreement with the model.

Isotopologues	Assumed $T$	Abundance	References	RD
	range, K	ratio, per cent		agreement?444Agreement between observations and calculations. If no, the closest calculated value is indicated.
HDO/H2O	100	0.3–18	(1, 2)	yes
HDO/H2O	$>$100	0.01–8	(2, 3, 4)	yes
D2O/HDO	$>$100	0.7–1.7	(4)	yes
D2O/H2O	$>$100	1.4E-3–2.9E-3	(4)	yes
NH2D/NH3	100	4–33	(5, 6, 7)	yes
NHD2/NH3	100	2.6–3.0	(5)	1.9
ND3/NH3	. . .	0.1	(8)	yes
N2D+/N2H+	100	0.2–36	(9, 10)	yes
DCO+/HCO+	. . .	0.39–14	(5, 6, 11, 12, 13)	yes
DCN/HCN	$>$25	0.47–19.5	(6, 11, 12, 13, 14)	yes
DNC/HNC	$>$25	1.5–6.1	(12, 13)	yes
HDS/H2S	$>$20	5–15	(12)	yes
HDCO/H2CO	100	1.5–45	(12, 13, 14, 15, 16)	yes
D2CO/HDCO	100	3–180	(15, 16)	yes
D2CO/H2CO	100	1–76	(5, 10, 15, 16, 17)	yes
CH2DOH/CH3OH	$>$30	10–94	(16, 18, 19)	17555See Sect. 3.3.
CH3OD/CH3OH	$>$30	0.8–6.9	(16, 18, 19)	yes
CH3OD/CH2DOH	$>$30	4.8–11.4	(16)	26
CHD2OH/CH3OH	$>$30	7–39	(16, 18, 19)	1.7
CD3OH/CH3OH	$>$50	0.2–2.8	(19)	0.07
DCOOCH3/HCOOCH3	$>$100	$\approx 15$	(20)	4.6
CH2DOCH3/CH3OCH3	. . .	$\approx 15$	(21)	yes
C2D/C2H	. . .	6–27	(12, 13)	yes
C3D/C3H	. . .	2.7–6.2	(22)	yes
C3HD/C3H2	. . .	4.8–9.4	(22)	yes
C4D/C4H	. . .	1.3–2.3	(22)	yes
C4HD/C4H2	. . .	1.3–5.1	(22)	yes
DC3N/HC3N	. . .	2.0–5.1	(22, 23)	yes
DC5N/HC5N	. . .	2.2–4.0	(22)	yes

Equations

$$\mathrm{C}{\mathrm{H}}_2\mathrm{OD}+h\nu ⟶\mathrm{C}{\mathrm{H}}_2+\mathrm{OD},$$

$$\mathrm{C}{\mathrm{H}}_2\mathrm{OD}+h\nu ⟶\mathrm{CHD}+\mathrm{OH}.$$

$$b=\frac{N_{\mathrm{ice}}}{N_s}\times b_m,$$

$$k_{\mathrm{pd}}=\frac{\pi a^2{F}_{\mathrm{ph}}{Y}_{\mathrm{ph}}}{N_s},$$

Full text

††pagerange: Chemical fractionation of deuterium in the protosolar nebula–References††pubyear: 201X

1 Introduction

Deuterium fractionation $R_{\rm D}$ or D/H = [XD]/[XH] represents the enhancement of D content in molecules (the notation [X] stands for the relative abundance of chemical species X). $R_{\rm D}$ for molecules in Solar-system bodies and star-forming regions is often found to be much higher than the cosmic D/H abundance ratio of $\approx 10^{-5}$. To understand chemical processes in prestellar cores, protostellar envelopes, and protoplanetary discs, it is crucial to ‘decipher’ the information contained in the measured $R_{\rm D}$ values of molecules in different objects.

Deuterium fractionation in star-forming regions has recently become an increasingly active research field due to important observational, experimental, and theoretical advances. Overviews are available in, e.g., Albertsson et al. (2013), Taquet, Charnley & Sipilä (2014), Awad et al. (2014), and Ceccarelli et al. (2014). Specific attention has been paid to deuterium fractionation of water by Taquet et al. (2013b), Furuya et al. (2013), Albertsson, Semenov & Henning (2014), van Dishoeck, Herbst & Neufeld (2013), and Wakelam et al. (2014).

The modeling work by Taquet et al. (2013b) and Kalvāns & Shmeld (2013) indicates that deuterium-rich ice molecules mostly are concentrated on the outer surface of ice mantles residing on interstellar or circumstellar grains. To put it another way, ice molecules that are formed early in the evolution of a contracting cloud core (e.g., H2O) typically have low $R_{\rm D}$, while species that are synthesized in the later stages of prestellar cores typically, and are abundant on the outer surface, have a high D content (e.g., CH3OH). Moreover, the molecules in bulk ice are chemically uncoupled from the chemical processes on the surface and in the gas. The average $R_{\rm D}$ for bulk molecules is substantially lower than previously thought.

The different deuteration of surface and bulk-ice molecules can only be reproduced with astrochemical models that consider bulk-ice species as a separate phase. Recent deuterium chemistry models that consider this are those of Taquet, Charnley & Sipilä (2014) and Furuya et al. (2015). The former considers circumstellar envelopes in a physically detailed way and confirms that D/H is high for surface and low for bulk ice species. The monolayer structure of this model (Taquet, Ceccarelli & Kahane, 2012) allowed these authors determining the detailed composition of ice as a function of ice depth, measured in monolayers.

The above papers consider the icy mantles as rigid structures, where molecule diffusion and reactions occur only in the surface monolayer. However, photoprocessing of subsurface ice can be a major route for chemical synthesis. With the help of numerical simulations, such processes have been studied by, e.g., Kalvāns & Shmeld (2010); Garrod (2013), and Chang & Herbst (2014). Chemical processing of bulk ice may introduce changes in the D/H ratio for icy species (Kalvāns & Shmeld, 2013). Therefore, the specific aim of the present study is to investigate the deuterium fractionation of circumstellar molecules by using a model that considers active subsurface ice chemistry. For this purpose, a model was employed that considers time-dependent physical conditions in an infalling protostellar envelope (Sect. 2.1), ice depth-dependent composition of mantle layers on grains, and an extensive reaction network of deuterated molecules (Sect. 2.2). This model adds to the understanding of the processes that govern deuterium fractionation in gas, surface, and mantle in general (Sect. 3.2), and for specific compounds (Sect. 3.3). The outcome of these studies is summarized in Sect. 4.

2 Methods

Our numerical simulation is based on the model ‘Alchemic-Venta’ (Kalvāns, 2015b, c). A short description of the model is presented below; it includes a detailed description for items that differ from Kalvāns (2015b).

2.1 Physical model

Following the approach of Garrod & Herbst (2006), we adopt a two-stage hot-core model in order to represent a contracting prestellar molecular cloud core (Phase 1, integration time $t\leq 0$) and the subsequent heating of a circumstellar nebula (Phase 2, $t>0$). With this we assume that the protostar is formed at $t=0$.

The simulation is performed for a single point (0D model) at the centre of a spherical core with radius $r_{1}$. The core was assumed to have a Plummer-like density profile with radius $r_{0}=1.3\times 10^{5}$ AU for the central density plateau and power-law parameter $\eta=5.0$. The core is embedded in a surrounding cloud with a fixed column density $N_{\rm out}=1\times 10^{21}$ cm*-2*. In Phase 1, the core undergoes free-fall gravitational collapse from an initial central density of $n_{\rm H,0}=2.3\times 10^{3}$ cm*-3*, which increases to $10^{7}$ cm*-3* over a time-scale of 1.5 Myr. The central plateau radius $r_{0}$ decreases proportionally to $n_{\rm H}^{-1/2}$ (Keto & Caselli, 2010; Taquet, Charnley & Sipilä, 2014) until it reaches $2\times 10^{3}$ AU, representing a centrally highly concentrated core.

These quantities were used to calculate the total column density $N_{\rm H}$ in the center of the spherical cloud core according to the Eq. (1) of Kalvāns (2015b). The obtained value of $N_{\rm H}$ is then used to calculate the interstellar extinction $A_{\rm V}$ (mag), for a line of sight through the center of the core. This $A_{\rm V}$ value is twice as high as that from cloud edge to its center. We use a conventional $N_{\rm H}/A_{\rm V}$ ratio of $2.0\times 10^{21}$ (see Valencic & Smith, 2015, and references therein). For calculating $A_{\rm V}$, we consider $N_{\rm out}$ as well as the matter inside $r_{1}$. Our simulation starts with $A_{\rm V}=1.34$ mag, corresponding to about 0.67 mag in the simulation of Furuya et al. (2015).

For the cold core phase, it was assumed that the cloud is externally heated, i.e., its temperature is governed by interstellar extinction. Dust temperature $T_{\rm dust}$ was calculated as a function of $A_{\rm V}$ with the method outlined by Garrod & Pauly (2011). The relation between gas temperature $T_{\rm gas}$ and $A_{\rm V}$ was derived from the data of the 3D collapsing cloud models by Hocuk et al. (2016). The difference between the two temperatures becomes less than 0.5 K when $A_{\rm V}>17$ mag and they equalize when $A_{\rm V}>89$ mag (remember that with $A_{\rm V}$ here we mean extinction along the whole line of sight).

In the second stage, the modeled gas parcel undergoes heating up to 200 K by energy influx from the newly formed protostar according to the T2 profile of Garrod & Herbst (2006). The warm-up time-scale is related to the lifetime of the dense envelope around the protostar. Infalling envelopes are thought to exist around protostars of Classes 0 and I, and around objects with flat spectral energy distribution (Stahler & Palla, 2005; Hartmann, 2009). Evans et al. (2009) estimate that the median lifetimes for these objects are 0.1–0.16 Myr, 0.44–0.54 Myr, and 0.35–0.40 Myr, respectively. Maury et al. (2011) estimate a median lifetime of 0.04–0.09 Myr for Class 0 objects. This means that the combined median lifetime of the envelope could lie in the range between 0.83 and 1.10 Myr although there are significant uncertainties (Schnee et al., 2012; Carney et al., 2016). The observations, which we use for comparing results (Table 4), have targeted both Class 0 and Class I objects.

For the purpose of the present model, we adopt a warm-up time-scale of 0.5 Myr, approximately half of the entire lifetime of a typical collapsing envelope. Volatile icy species, such as CO and CH4, evaporate near the end of the Class 0 phase, and the icy mantles are completely evaporated 0.4 Myr after starbirth, well into the Class I phase.

The gas parcel in consideration was assumed to be fully shielded from the protostellar radiation. Meanwhile, it was assumed that the radius of the core linearly decreased during the warm-up phase from the cloud core value of $r_{1}=10^{4}$ AU to $10^{3}$ AU, which is closer to the characteristic size of a protostellar nebula (e.g. Visser, Doty & van Dishoeck, 2011). Such a simple isopycnic (constant density) contraction was introduced to account for the mass loss in the circumstellar envelope, as it falls into the star and the protoplanetary disc. The resulting change in $A_{\rm V}$ has nearly no effect on the results. Fig. 1 shows the evolution of physical conditions in the modeled cloud core.

Table 1 summarizes the macrophysical parameters at key moments of the simulation. Because $n_{\rm H}$, $r_{1}$, and $r_{0}$ vary with time, the core mass $M_{\rm core}$ encompassed within $r_{1}$ is not constant.

2.2 Chemical network

The model adopts the publicly available full deuterium reaction network from Albertsson et al. (2013). This network was generated from the osu.2009 reaction database and does not discern between H and D atoms that are attached to different heavy atoms in a single molecule. In other words, H and D are interchangeable as they freely undergo intramolecular diffusion in reactions. For example, the cosmic-ray-induced photodissociation of the hydroxymethyl radical $\rm CH_{2}OD$ is included in the database as two reactions:

[TABLE]

[TABLE]

The second product set requires an intramolecular exchange of H and D atoms before the molecule is split in two. This behavior is not always possible, especially, when low energies and large molecules are involved. However, in the reaction set, intramolecular mixing of H and D is assumed to occur for all reactions, regardless of their type, reactants and products.

The latter aspect decreases the extent of selective deuterium fractionation. For example, let us assume that CH2OH is produced via CH2+OH surface reaction. If CH2 is more enriched in deuterium than OH, then this information is lost in the synthesis, as well as when the molecule disintegrates via reactions (1). This means that $R_{\rm D}$ for the C–H bond is lowered and that of the O–H bond is increased by fictitious intramolecular diffusion. However, H and D are not fully interchangeable in the interstellar medium. This is manifested by the low CH3OD/CH2DOH abundance ratio, $\approx$10 per cent, as observed by Parise et al. (2006). With H and D interchangeable, such a low ratio cannot be reproduced.

To change this unsatisfying situation and achieve more realistic results in abundance ratios of complex organic molecule (COM) isotopologues, we revised the reaction network by hand. The products that involve intramolecular diffusion of H and D were removed for reactions involving CH2OH, CH3OH, HCOOCH3, HCOOH, CH3CHO, CH3OH+, CH3OH${}_{2}^{+}$, CH3O(H)CH${}_{3}^{+}$, and CH3CHO+. For the methyl formate cation, the species COOCH${}_{4}^{+}$ were replaced with HCOOCH${}_{3}^{+}$. We also add methoxy radical CH3O to the network. For the latter, necessary reaction data was adopted from the network employed by Garrod, Weaver & Herbst (2008). These changes do not fully eliminate the interchangeability for H and D in COMs.

Several additional changes were introduced in this network. The desorption energy $E_{D}$ for H2 was adopted to be 430 K (Cazaux & Spaans, 2004; Heine, Zhechkov & Seifert, 2004; Garrod & Herbst, 2006). For HD and D2, the values were assumed to be higher by 2 and 32 K, respectively (Kristensen et al., 2011). $E_{D}$ was taken to be 450 K for H and 471 K for D (Caselli et al., 2002; Aikawa et al., 2012). The network was updated with reactions 5–12, 14, and 16–19 from Kalvāns (2015b). The deuterium analogues for these reactions were included, with an exception of $\rm D+CO\longrightarrow DCO$. The activation barrier of this reaction (1400 K) in the network of Albertsson et al. (2013) is already lower than the 1600 K assumed by Kalvāns (2015b).

Another parameter, which may affect $R_{\rm D}$ of interstellar and circumstellar molecules, is the ortho/para ratio of the H2 molecule (Walmsley, Flower & Pineau des Forêts, 2004; Flower, Pineau des Forêts & Walmsley, 2004; Flower, Pineau Des Forêts & Walmsley, 2006b, a). Molecular hydrogen formed on the grains is released into the gas phase with the statistical ortho/para ratio of 3:1. The subsequent interaction with H+ then reduces this value towards thermal equilibrium value (Pagani et al., 2009; Sipilä, Caselli & Harju, 2013, 2015). This happens even before the formation of the dense core. In the present study, we do not consider the spin states of H2 or other species. This approach eases calculations of the nearly 80,000 molecular processes in the model and can be justified by the following discussion.

First, ice formation in the present model occurs in a dense medium, where the H2 ortho/para ratio is assumed to be low, and therefore has a limited practical effect. More than 90 per cent of ice molecules are deposited when $A_{\rm V}>4$ mag and ortho/para ratio is around $10^{-3}$, according to the models of Sipilä, Caselli & Harju (2013) and Furuya et al. (2015).

Second, the main source for a high abundance ratio of gas-phase neutral D atoms with respect to similar H atoms (atomic D/H) in the prestellar Phase 1 is HD photodissociation, as in Furuya et al. (2015). The photodissociation of HD occurs faster than that of H2 due to shielding effects. As a result, the atomic D/H in the modeled (translucent) core is higher than the cosmic D/H ratio and largely independent of binary gas-phase reactions involving the isotopologues of H${}_{3}^{+}$, where the H2 ortho/para ratio is important. The latter aspect means that H2 spin states do not have a decisive effect on ice deuteration in diffuse medium either.

A third point regarding the effect of H2 spin states on ice deuteration is that there is considerable evidence that the ortho/para ratio decreases when H2 interacts with interstellar ice analogues (Le Bourlot, 2000; Watanabe et al., 2010; Sugimoto & Fukutani, 2011; Chehrouri et al., 2011; Hama et al., 2012). As far as we know, this effect has not yet been considered in astrochemical models. Thanks to such interactions, the relaxation time of the H2 ortho/para ratio in cold cores is reduced, reducing also the spatial and temporal extent of zones with high ortho/para ratio. This adds to the uncertainties regarding the H2 ortho/para ratio, which limits the usefulness of the inclusion of H2 spin states in the present study.

2.3 Chemical model

Table 2 shows the initial abundances of chemical species. The adopted abundance of deuterium nuclei is $1.0\times 10^{-5}$ relative to H nuclei, consistent with findings that the overall cosmic D/H is lower in regions with higher H column density and lower temperatures (Linsky et al., 2006; Prodanović, Steigman & Fields, 2010). Abundances for major elements were taken from Garrod, Weaver & Herbst (2008), those for elements heavier than oxygen – from Aikawa & Herbst (1999). This means that we consider a medium where refractory elements are significantly depleted, while the ice forming species are not. This approach is similar to that of Furuya et al. (2015).

Self-shielding of HD, and its mutual shielding by H2 was added (Le Petit, Roueff & Le Bourlot, 2002; Wolcott-Green & Haiman, 2011), in addition to the shielding of gaseous and adsorbed H2, CO, and N2. For consistency, we include self and mutual shielding for D2 in the same manner as for HD. The latter aspect has little effect, because of the low abundance of D2.

Neutral gaseous species can be adsorbed on to interstellar grains with radius $a+b$, where $a=0.1\,\mu$m is the radius of grain nucleus and $b$ is ice thickness. The latter is a time-dependent quantity, expressed as

[TABLE]

where $N_{\rm ice}$ is the total number of ice molecules per grain, $N_{s}=1.5\times 10^{6}$ is the number of adsorption sites on the grain surface and also per monolayer (ML), and $b_{m}=3.5\times 10^{-8}$ cm is the assumed monolayer thickness.

The ice mantle was described as consisting of four layers – one surface and three subsurface layers (‘sublayers’). Chemical reactions may occur in all sublayers, in line with the approach developed by Kalvāns (2015c). Dissociation of icy species by interstellar and cosmic-ray-induced photons is possible; it is assumed that subsurface species are shielded by the layers above. Each ML has a radiation attenuation probability of 0.007 (Andersson & van Dishoeck, 2008).

We note that experiments show that the evaporation of icy mixtures occurs sequentially, i.e., through diffusion from subsurface ice layers (Martín-Doménech et al., 2014), although partial entrapment is possible (Fayolle et al., 2011). This supports our model with mobile bulk-ice molecules. Chemical photoprocessing of such molecules is important because it can significantly influence the abundances of minor icy species, such as COMs (e.g., Kalvāns & Shmeld, 2010; Garrod, 2013).

The rate of desorption of surface molecules is governed by their adsorption energies $E_{D}$. These correspond to bulk absorption energies of $E_{B}=3E_{D}$ for species in the sublayers of the icy mantle (Kalvāns, 2015b). Diffusion is governed by binding energies $E_{b,s}=xE_{D}$ and $E_{b,m}=xE_{B}$ for ice surface and mantle (sublayer) phases, respectively. The value of the parameter ’$x$’ was taken to be 0.35 (Garrod & Pauly, 2011). The various energies of the icy species are summarized in Table 3. The present study does not consider quantum tunneling for molecule diffusion on the surface and in the ice, which agrees with Katz et al. (1999) and Watanabe et al. (2010). Inter-sublayer diffusion for ice molecules was included (Kalvāns, 2015c), as it is essential in transferring molecules from bulk ice to the surface when the ice temperature rises in Phase 2.

Six desorption mechanisms for species in the surface layer were considered – desorption by interstellar photons, cosmic-ray-induced photodesorption, reactive desorption, indirect reactive desorption, evaporation, and desorption by cosmic-ray-induced whole-grain heating. The efficiency for indirect reactive desorption – or desorption by H+H surface reaction heat – was calibrated using the approach explained in detail in the previous works (Kalvāns, 2015a, c). Evaporation and photodesorption of species in subsurface layers was allowed if the number of the above MLs was lower than unity.

Because the molecular desorption energy $E_{D}$ is lower than the energy required for bulk diffusion, $E_{b,m}$, ice evaporation rate at rising temperatures is limited by diffusion, not thermal desorption. This is in an agreement with the experimental results by Öberg et al. (2009c) and Mispelaer et al. (2013). Such a behavior arises because in dense icy mixtures the (evaporating) surface becomes saturated with refractory species and the evaporation is hampered, unless the volatile molecules are able to diffuse efficiently from the bulk ice to the surface layer.

The rate coefficient (s*-1*) for desorption by interstellar and secondary photons was calculated according to

[TABLE]

where $F_{\rm ph}$ (s*-1cm-2*) is photon flux ($1.7\times 10^{8}$ cm*-2* s*-1* for interstellar photons at $A_{\rm V}=0$ mag and 4875 cm*-2* s*-1* for cosmic-ray-induced photons). Photodesorption yield $Y_{\rm ph}$ was taken to be $3\times 10^{-4}$ for interstellar photons (Arasa et al., 2015; Furuya et al., 2015) and a 1.5 times lower value for cosmic-ray-induced photons (Kalvāns, 2015c). Eq. 3 is physically more adequate than the approach of Kalvāns (2015c) because the ratio between $N_{s}$ and grain cross section $\pi a^{2}$ is constant for spherical grains of different sizes. In other words, Eq. (3) means that the actual photodesorption rate is independent of the grain size and ice thickness $b$. The calculated photodesorption rates are lower than in our previous studies. A standard cosmic-ray ionization rate of $1.3\times 10^{-17}$ s*-1* was used.

3 Results

In this section, we briefly mention the general processes governing ice chemistry. Then, we turn to the novelties of our deuterium chemistry model before focusing on the calculation results in the context of observations of circumstellar deuterated molecules. We remind that the simulation starts at $t=-1.5$ Myr, while the protostar is formed at 0 Myr, and the envelope is modeled up to 0.5 Myr, at which point it reaches a maximum temperature of 200 K. We try to attribute chemical processes to gas or solid (ice) phases although these phases are closely intertwined. For example, while the most efficient deuteration occurs in the gas-phase, depletion of heavy species and formation of H2 on grain surfaces are prerequisites for this process.

3.1 Evolution of the ice layer

The first water ice monolayer is formed at $A_{\rm V}=3.0$ mag ($t=-0.32$ Myr). This is consistent with observations (e.g., Whittet et al., 2001) and some models (Garrod & Pauly, 2011; Taquet, Charnley & Sipilä, 2014; Hocuk & Cazaux, 2015). Ice appears later than in the model of Furuya et al. (2015). This is largely because the simulation of Furuya et al. starts at an earlier point in cloud evolution and takes longer before the onset of rapid core collapse.

General ice chemistry and processes with similar model for star-forming cores have been described before (Kalvāns, 2015a, c). Fig. 2 shows the calculated abundances for major ice species in specific sublayers. The evolution of total ice thickness can also be seen. After the appearance of the first monolayer of adsorbed species on grain surface, the freeze-out proceeds at an increasing rate thanks to the growing $n_{\rm H}$ and decreasing $T_{\rm dust}$. Species with high $E_{D}$ (e.g., water and ammonia) become depleted at $t\approx-0.05$ Myr. The depicted temporal evolution of the ice layer on grains provides the context for the discussion that follows. Five species – H2O, CO, CO2, NH3, and N2 – are major ice ingredients. Together they constitute 98 per cent of all ice molecules at $t=0$. In the model, CO and N2 are products of gas-phase chemistry, while $>90$ per cent of H2O, CO2, and NH3 molecules are formed via surface reactions.

Ice species are initially formed either in the gas phase, and then get adsorbed, or directly formed on the surface, before they are buried below ice surface. In total, about 10 per cent of molecules are chemically transformed in bulk ice. The subsurface (and other) reactions are driven by interstellar photons in the translucent cloud core (up to $A_{\rm V}\approx 7$ mag; $t\leq-0.16$) and cosmic-ray-induced photons in the dark core and the envelope. The subsequent heating of the protostellar envelope releases the icy species. Calculation results show that CO, CO2, and H2O evaporate at approximate temperatures of 24, 58, and $>$115 K, respectively. Because the molecular diffusion energy in bulk ice $E_{b,m}$ is higher than desorption energy $E_{D}$, species often diffuse to the surface and evaporate over an extended period of time.

3.2 Deuterium chemistry in ice

Fig. 3 shows the calculated deuterium fractionation of the most important hydrogenated ice species for surface and subsurface mantle (sublayers 1, 2, and 3). The limited time-span for $R_{\rm D}$ curve of bulk ice species indicates the period of existence for the subsurface layers. The figure illustrates the concentration of deuterated species on the surface. The HDO/H2O surface abundance ratio can reach 3 per cent, although average HDO/H2O in ice does not exceed 0.2 per cent at any time when a considerable ice mass ($\geq 1$ ML) is present. This is in agreement with the non-detection of solid HDO in the interstellar or circumstellar medium. The upper limit for HDO$/$H2O abundance ratio in ice is 2 per cent (Parise et al., 2003; Dartois et al., 2003). $R_{\rm D}$ for water ice (bulk) depends on the availability of atomic D for reactions on grain surfaces during ice accumulation in Phase 1.

Deuterium fractionation of surface molecules in the translucent core is affected by two competing effects. First, photodissociation of HD and D2 occurs at a higher rate than that of H2, because of self-shielding effects. This ensures elevated atomic D/H, facilitating deuteration (especially visible in the early stages in Figs. 3 and 4; see also Sect. 2.2). Second, species on the surface undergo repeated photodissociation and recombination (which occurs faster with H atoms), before being incorporated in the bulk ice. This results in a low molecular $R_{\rm D}$. Both effects were initially found by Furuya et al. (2015).

The dip in the ice deuteration curves, visible for all species depicted in Fig. 3 at $t\approx-0.10$ Myr, arises because the ice formation epoch partially overlaps with a rapid synthesis of molecular hydrogen (including HD) on grains. D atoms are consumed in surface reactions with atomic H faster than they are accreted from the gas phase. Because the majority of D atoms goes into HD, ice species formed during this period are poor in deuterium. This coincidence at $-0.10$ Myr is illustrated with Fig. 4.

Afterwards, the freeze-out removes all heavy hydrogenated species from the gas. The abundance of H and D atoms is low because hydrogen has been locked in molecules. Immediately after the rapid accretion phase, conditions for efficient deuterium fractionation are maintained by gas-phase chemistry. These deuteration processes, however, affect only the outer surface of ices.

3.3 Deuterium fractionation of gaseous species

Chemical models tend to produce higher-than-observed maximum gas-phase abundances for major ice species, such as water, that are evaporating in the vicinity of a protostar (e.g. Garrod & Herbst, 2006; Taquet, Charnley & Sipilä, 2014; Wakelam et al., 2014). This discrepancy probably arises because the observations often have a limited resolution and sample a whole line-of-sight in the nebula. Despite this, the observed molecular ratios, such as $R_{\rm D}$, can be reliable. The necessary condition is that most of the molecules of the observed species occur in regions with similar history and physical conditions. This can be verified by accurately evaluating the excitation temperatures of the different isotopologues (Nishimura et al., 2013). While calculated $R_{\rm D}$ values are our main result, we also present calculated abundances. This allows evaluating whether the $R_{\rm D}$ values are truly relevant.

We compare the numerical simulation results against observed gas-phase molecular $R_{\rm D}$ in low-mass protostars. It is assumed that the temporal evolution of calculated abundances may qualitatively represent a 1D spatial structure of the protostellar core (Garrod, Weaver & Herbst, 2008). This means that higher temperatures in the warm-up stage of the model represent regions closer to the protostar.

The observational data on abundance ratios for isotopologues of different species have been summarized in Table 4. The calculated abundance and $R_{\rm D}$ curves are shown graphically in figures including those in Appendix A (available in electronic form). Observations with only upper limits reported were not considered here. For some species (e.g., water), the temperature ranges for Table 4 were specified in the reference papers. For others, it was usually assumed that the excitation temperatures represent the lower limit of the kinetic temperature. Because the excitation temperatures can be very low, this condition does not always constrain the temperature range. The temperature range indicated for molecular $R_{\rm D}$ in Table 4 is only authors’ estimate to help classifying literature data for the purposes of the present paper.

Fig. 5 shows that for gas-phase water and ammonia, there are several peaks in the deuterium fractionation ratios. The first peak occurs a few kyr before starbirth, when the remaining gaseous heavy molecules are being completely depleted during freeze-out. Efficient deuterium enrichment occurs in the cold gas of the dense core (Millar, Bennett & Herbst, 1989; Roberts & Millar, 2000; Roberts, Herbst & Millar, 2004). While this peak corresponds to Phase 1 in our model, it might be relevant to the outer, colder parts of early circumstellar envelopes.

The second peak at 0.14 Myr (22 K) coincides with the evaporation of solid N2 and CO. These species re-appear in the gas phase and are partially converted to N2D+ and DCO+ by the products of cosmic-ray-induced dissociation and ionization of H2 and He. N2D+ and DCO+ then serve as deuteration agents. Meanwhile, the presence of abundant CO now removes the heavy isotopologues of H${}_{3}^{+}$ that were the main deuteration agents in the cold phase.

The third, smaller peak at 0.2 Myr is caused by evaporation of methane and its D isotopologues. Significant amounts of methane ice are in inner ice layers, near the grain nuclei and have to diffuse to the surface before evaporating. The total evaporation of methane takes almost 20 kyr. Methane can be transformed to CH${}_{3}^{+}$ by radicals and ions in the gas. Most of these are generated by cosmic-ray induced dissociation and ionization of H2 and He. CH${}_{3}^{+}$ isotopologues are affected by reactions that are selective with respect to H and D atoms (Table 13 of Albertsson et al., 2013). The high deuteration of the abundant methane and its associated ions and radicals is transferred to other species via hydrogen atom exchanges or dispatchment of free H and D atoms. Because of these factors, the $R_{\rm D}$ peak, induced by methane evaporation, is about 50 kyr long and ends only when most of CH4 has been converted to CO.

This third $R_{\rm D}$ peak is followed by a number of smaller peaks caused by the successive evaporation of other molecules, such as C2H2, C2H4, H2S, and CH3OCH3. The gas-phase processing of these species enable additional gas-phase deuterium fractionation reaction chains. The final gas-phase $R_{\rm D}$ peak at $t=0.34$ Myr (94 K) is related to the evaporation of water and ammonia themselves. The peak occurs before the bulk of these species has evaporated because deuterium-rich molecules are concentrated on the outer icy surface and evaporate first. After evaporation, $R_{\rm D}$ is mainly governed by gas-phase chemistry, much like the case of methane. With the evaporation of water, deuterated hydronium H3O+ replaces DCO+ as the main molecular deuteration agent. Similar $R_{\rm D}$ peaks can be discerned for most of the other species. Appendix Fig. 7 shows the calculated results for other (in)organic species observed in protostellar envelopes.

In the protostellar Phase 2, we obtain a higher gas-phase D2O/HDO abundance ratio than HDO/H2O for most of the time (cf. Fig. 5), which agrees with the recent observations by Coutens et al. (2014). We did not find similar results among recent papers on modeling of deuterium chemistry in circumstellar envelopes (Aikawa et al., 2012; Awad et al., 2014; Taquet, Charnley & Sipilä, 2014). The cause of the high D2O/HDO abundance ratio is gas-phase chemistry – this ratio is inherited from the high deuteration of the trihydrogen cation. Its four isotopologues H${}_{3}^{+}$, H2D+, HD${}_{2}^{+}$, and D${}_{3}^{+}$ have almost equal relative abundances of $\approx 10^{-11}$ at the beginning of Phase 2, thanks to the reactions proposed by Roberts, Herbst & Millar (2004). This deuterium enrichment is transferred to water and other molecules via ion-neutral reactions. Because a significant part of many species, including water, is formed on surfaces via the addition of atomic D, the $R_{\rm D}$ of these species never attains the high maximum values characteristic for purely gas-phase products, such as the formyl cation (DCO+/HCO${}^{+}=58$ per cent at 0.14 Myr) or diazenylium (N2D+/N2H${}^{+}=81$ per cent).

Deuterated analogues of several COMs – methanol CH3OH, formaldehyde H2CO, dimethyl ether CH3OCH3, and methyl formate HCOOCH3 – have been observed near low-mass protostars (Table 4). Fig. 6 shows the calculated abundances and $R_{\rm D}$ for methanol and formaldehyde deuterated isotopologues. Data for other COMs are visualized in Fig. 8. For CHD2OH, CD3OH, and DCOOCH3 the high observed $R_{\rm D}$ values are not reached. In the case of CH2DOH, the lowest $R_{\rm D}$ margin from observations is 10 per cent. This number is nominally reached in the model at $t=0.26$ Myr and $T=60$ K. At this temperature, methanol does not evaporate yet and its gas-phase relative abundance is below $10^{-10}$. This is orders of magnitudes lower than its maximum calculated abundance of $1.4\times 10^{-7}$ during evaporation at 124 K. At these later stages, $R_{\rm D}$ for CH2DOH is only about 0.1 to 1 per cent. This discrepancy means that our model does not adequately reproduce the CH2DOH/CH3OH abundance ratio. A major reason for this disagreement is the synthesis of deuterium-poor methanol in bulk ice.

The synthesis of more than 90 per cent methanol CH3OH, dimethyl ether CH3OCH3, and methyl formate HCOOCH3 molecules occurs in the subsurface bulk ice layers (Kalvāns, 2015c, b). The latter two species have also gas-phase formation pathways that become effective when methanol and formaldehyde have evaporated (Garrod & Herbst, 2006). We use the low activation energies for CO hydrogenation from Kalvāns (2015b). This allows to reproduce the abundances of methanol and its daughter species. At the same time, formaldehyde H2CO is largely consumed in ice before evaporation and is unable to reach the high relative abundance of about $1\times 10^{-7}$ observed in protostellar envelopes (Ceccarelli et al., 2000, 2001; Schöier et al., 2002; Maret et al., 2004). Formaldehyde is produced in gas-phase to attain a relative abundance of nearly 10*-9* (see Garrod & Herbst, 2006).

The deuterium fractionation of formaldehyde H2CO is reproduced in our model, which has been a problem for some of the previous models (e.g. Aikawa et al., 2012; Taquet, Charnley & Sipilä, 2014). Similarly to the case of D2O, the high efficiency of formaldehyde deuteration is thanks to to gas-phase chemistry with the network of Albertsson et al. (2013).

The deuteration of carbon-chain species is reproduced correctly (cf. Table 4 and Fig. 9). The abundances and deuteration of some organic species become almost steady after the evaporation events, thanks to high-temperature ion-neutral reactions (e.g. Willacy & Woods, 2009; Awad et al., 2014). In our model, carbon chains are formed mainly in reactions on grain outer surfaces. Exceptions are cyanides HC3N and HC5N (also inorganic species HCN and HNC) that can be formed in substantial amounts in subsurface ice layers rich in CO, N2, or NH3 (Fig. 7; see also Kalvāns, 2015b).

4 Conclusions

The scope and results of the simulation are broadly similar to those reported by Taquet, Charnley & Sipilä (2014). A major difference was that we consider bulk-ice chemical processes. It was found that subsurface photoprocessing facilitates the synthesis of deuterium-poor species in bulk ice. Therefore, bulk ice chemistry does not help in explaining the very high observed $R_{\rm D}$ ratios for some COMs, contrary to a suggestion by Taquet, Charnley & Sipilä (2014). Solid deuterium-rich organic species are confined to ice surfaces and those abundant in the bulk ice have lower D/H.

With the cloud core contracting, darkening, and cooling in the prestellar Phase 1, atomic D is largely lost to synthesis of the HD molecule on grain surfaces. In our model, this process coincides with a period of rapid ice accumulation. Such competition for adsorbed D atoms results in that the majority of ice molecules are formed with relatively low $R_{\rm D}$.

The sequential evaporation of icy molecules in a warming protostellar envelope (Phase 2) affect the deuteration of many gaseous species. This happens either by the introduction of abundant deuteration agents in the gas (e.g., N2D+) or by providing abundant hydrogenated species (e.g., methane) that are processed by selective deuteration reactions to temporarily achieve high $R_{\rm D}$ level, which, in turn, is transferred to other species in reactions involving chemical radicals and ions. This results in species’ $R_{\rm D}$ that is rapidly changing in the modeled temperature range of a protostellar envelope.

The evolution of molecular $R_{\rm D}$ with increasing temperature consists of a series of peaks.In the protostellar envelope, a combination of high gas-phase abundance and high D/H ratio is observed for species that are undergoing evaporation. Additionally, $R_{\rm D}$ can be strongly affected by the evaporation of other species that take part in hydrogenation reaction chains (Sect. 3.3).

Acknowledgements

We are grateful to the anonymous referees for many valuable comments and suggestions that made this paper better. JK, IS, and JRK thank Ventspils City Council for its persistent support. This research has made use of NASA’s Astrophysics Data System.

Appendix A Calculated abundances and fractionation ratios for deuterated species observed around young stellar objects.