TY - JOUR
AU - Hocuk,, S.
AB - Abstract Understanding the 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 chemistry. 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 pre-stellar core hamper 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 a high molecular D/H ratio and a high gas-phase abundance. astrochemistry, molecular processes, stars: formation, ISM: molecules 1 INTRODUCTION Deuterium fractionation RD 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). RD 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 ≈10−5. To understand chemical processes in pre-stellar cores, protostellar envelopes and protoplanetary discs, it is crucial to ‘decipher’ the information contained in the measured RD 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, for example, 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 modelling 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 RD, while species that are synthesized in the later stages of pre-stellar 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 RD 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 et al. (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 (ML) structure of this model (Taquet, Ceccarelli & Kahane 2012) allowed these authors to determine the detailed composition of ice as a function of ice depth, measured in MLs. The above papers consider the icy mantles as rigid structures, where molecule diffusion and reactions occur only in the surface ML. 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, for example, 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 this 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 (Section 2.1), ice depth-dependent composition of mantle layers on grains and an extensive reaction network of deuterated molecules (Section 2.2). This model adds to the understanding of the processes that govern deuterium fractionation in the gas, surface and mantle in general (Section 3.2), and for specific compounds (Section 3.3). The outcome of these studies is summarized in Section 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 pre-stellar molecular cloud core (Phase 1, integration time t ≤ 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 r1. The core was assumed to have a Plummer-like density profile with radius r0 = 1.3 × 105 au for the central density plateau and power-law parameter η = 5.0. The core is embedded in a surrounding cloud with a fixed column density Nout = 1 × 1021 cm−2. In Phase 1, the core undergoes a free-fall gravitational collapse from an initial central density of nH, 0 = 2.3 × 103 cm−3, which increases to 107 cm−3 over a time-scale of 1.5 Myr. The central plateau radius r0 decreases proportionally to |$n_{\rm H}^{-1/2}$| (Keto & Caselli 2010; Taquet et al. 2014) until it reaches 2 × 103 au, representing a centrally highly concentrated core. These quantities were used to calculate the total column density NH in the centre of the spherical cloud core according to the equation (1) of Kalvāns (2015b). The obtained value of NH is then used to calculate the interstellar extinction AV (mag), for a line of sight through the centre of the core. This AV value is twice as high as that from the cloud edge to its centre. We use a conventional NH/AV ratio of 2.0 × 1021 (see Valencic & Smith 2015, and references therein). For calculating AV, we consider Nout as well as the matter inside r1. Our simulation starts with AV = 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 Tdust was calculated as a function of AV with the method outlined by Garrod & Pauly (2011). The relation between gas temperature Tgas and AV 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 AV > 17 mag, and they equalize when AV > 89 mag (remember that with AV here we mean extinction along the whole line of sight). In the second stage, the modelled 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 a 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, 0.44–0.54 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 of 0.83–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 r1 = 104–103 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 AV has nearly no effect on the results. Fig. 1 shows the evolution of physical conditions in the modelled cloud core. Figure 1. Open in new tabDownload slide Evolution of temperature, density and interstellar extinction. Density, gas temperature and dust temperature are plotted for the modelled centre part of the cloud core. AV corresponds for a line of sight passing through the centre. The protostar is formed at time t = 0. Figure 1. Open in new tabDownload slide Evolution of temperature, density and interstellar extinction. Density, gas temperature and dust temperature are plotted for the modelled centre part of the cloud core. AV corresponds for a line of sight passing through the centre. The protostar is formed at time t = 0. Table 1 summarizes the macrophysical parameters at key moments of the simulation. Because nH, r1 and r0 vary with time, the core mass Mcore encompassed within r1 is not constant. Table 1. Density at the centre of the cloud core, radius r1 of the core and core mass at the start, at starbirth time, and at the end of the simulation. t (Myr) . nH (cm−3) . r1 (au) . Mcore (M⊙) . −1.5 2.3 × 103 1.0 × 104 0.035 0 1.0 × 107 1.0 × 104 1.18 0.5 1.0 × 107 1.0 × 103 – t (Myr) . nH (cm−3) . r1 (au) . Mcore (M⊙) . −1.5 2.3 × 103 1.0 × 104 0.035 0 1.0 × 107 1.0 × 104 1.18 0.5 1.0 × 107 1.0 × 103 – Open in new tab Table 1. Density at the centre of the cloud core, radius r1 of the core and core mass at the start, at starbirth time, and at the end of the simulation. t (Myr) . nH (cm−3) . r1 (au) . Mcore (M⊙) . −1.5 2.3 × 103 1.0 × 104 0.035 0 1.0 × 107 1.0 × 104 1.18 0.5 1.0 × 107 1.0 × 103 – t (Myr) . nH (cm−3) . r1 (au) . Mcore (M⊙) . −1.5 2.3 × 103 1.0 × 104 0.035 0 1.0 × 107 1.0 × 104 1.18 0.5 1.0 × 107 1.0 × 103 – Open in new tab 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 data base 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 CH2OD is included in the data base as two reactions: \begin{equation} \mathrm{CH_2OD} + h\nu \longrightarrow \mathrm{CH_2} + \mathrm{OD}, \end{equation} (1) \begin{equation*} \mathrm{CH_2OD} + h\nu \longrightarrow \mathrm{CHD} + \mathrm{OH}. \end{equation*} The second product set requires an intramolecular exchange of H and D atoms before the molecule is split in two. This behaviour is not always possible, especially when low energies and large molecules are involved. However, in the reaction set, the 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 the 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 RD 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, ≈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 the 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 were 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 ED 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). ED 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 1600 K assumed by Kalvāns (2015b). Another parameter, which may affect the RD of interstellar and circumstellar molecules, is the ortho/para ratio of the H2 molecule (Flower, Pineau des Forêts & Walmsley 2004; Walmsley, Flower & Pineau des Forêts 2004; Flower, Pineau Des Forêts & Walmsley 2006a,b). 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 a thermal equilibrium value (Pagani et al. 2009; Sipilä, Caselli & Harju 2013, 2015). This happens even before the formation of the dense core. In this 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 AV > 4 mag and the ortho/para ratio is around 10−3, according to the models of Sipilä et al. (2013) and Furuya et al. (2015). Secondly, 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 pre-stellar 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 modelled (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 a 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; Chehrouri et al. 2011; Sugimoto & Fukutani 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 extents of zones with a 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 this study. 2.3 Chemical model Table 2 shows the initial abundances of chemical species. The adopted abundance of deuterium nuclei is 1.0 × 10−5 relative to H nuclei, consistent with findings that the overall cosmic D/H is lower in regions with higher H column densities and lower temperatures (Linsky et al. 2006; Prodanović, Steigman & Fields 2010). Abundances for major elements were taken from Garrod et al. (2008), and those for elements heavier than oxygen were taken 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). 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 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 Open in new tab 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 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 Open in new tab 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-shielding and mutual shielding for D2 in the same manner as for HD. The latter aspect has a 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 μm is the radius of the grain nucleus and b is the ice thickness. The latter is a time-dependent quantity, expressed as \begin{equation} b=\frac{N_{\mathrm{ice}}}{N_{\rm s}}\times b_{\rm m}, \end{equation} (2) where Nice is the total number of ice molecules per grain, Ns = 1.5 × 106 is the number of adsorption sites on the grain surface and also per ML and bm = 3.5 × 10−8 cm is the assumed ML 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 ED. These correspond to bulk absorption energies of EB = 3ED for species in the sublayers of the icy mantle (Kalvāns 2015b). Diffusion is governed by binding energies Eb, s = xED and Eb, m = xEB 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. This study does not consider quantum tunnelling 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. 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) ED – a Absorption EB 3.0EDb Binding (diffusion) Eb, s 0.35EDc Binding (diffusion) Eb, m 0.35EB Surface . . Bulk ice (mantle sublayers) . Energy . Notation . Value . . Energy . Notation . Value . Desorption (adsorption) ED – a Absorption EB 3.0EDb Binding (diffusion) Eb, s 0.35EDc Binding (diffusion) Eb, m 0.35EB aData from the network made available by Albertsson et al. (2013). bKalvāns (2015a,b). cGarrod & Pauly (2011). Open in new tab 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) ED – a Absorption EB 3.0EDb Binding (diffusion) Eb, s 0.35EDc Binding (diffusion) Eb, m 0.35EB Surface . . Bulk ice (mantle sublayers) . Energy . Notation . Value . . Energy . Notation . Value . Desorption (adsorption) ED – a Absorption EB 3.0EDb Binding (diffusion) Eb, s 0.35EDc Binding (diffusion) Eb, m 0.35EB aData from the network made available by Albertsson et al. (2013). bKalvāns (2015a,b). cGarrod & Pauly (2011). Open in new tab 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 ED is lower than the energy required for bulk diffusion, Eb, m, the 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 behaviour 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 \begin{equation} k_{\mathrm{pd}}=\frac{\pi a^2 F_{\mathrm{ph}} Y_{\mathrm{ph}}}{N_{\rm s}}, \end{equation} (3) where Fph (s−1 cm−2) is the photon flux (1.0 × 108 cm−2 s−1 for interstellar photons at AV = 0 mag and 4875 cm−2 s−1 for cosmic-ray-induced photons). Photodesorption yield Yph was taken to be 3 × 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). Equation (3) is physically more adequate than the approach of Kalvāns (2015c) because the ratio between Ns and grain cross-section πa2 is constant for spherical grains of different sizes. In other words, equation (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 × 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 modelled 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, the 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 ML is formed at AV = 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 et al. 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 a rapid core collapse. General ice chemistry and processes with a 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 ML of adsorbed species on the grain surface, the freeze-out proceeds at an increasing rate, thanks to the growing nH and decreasing Tdust. Species with high ED (e.g. water and ammonia) become depleted at t ≈ −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. Figure 2. Open in new tabDownload slide Calculated ice thickness and abundance in MLs for major species in the sublayers. The numbers in parentheses indicate the number of the respective sublayer, while ‘Other’ stands for all other icy species in that sublayer. Note that in the model, molecules within a sublayer are fully mixed. Figure 2. Open in new tabDownload slide Calculated ice thickness and abundance in MLs for major species in the sublayers. The numbers in parentheses indicate the number of the respective sublayer, while ‘Other’ stands for all other icy species in that sublayer. Note that in the model, molecules within a sublayer are fully mixed. 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 the ice surface. In total, about 10 per cent of molecules are chemically transformed into bulk ice. The subsurface (and other) reactions are driven by interstellar photons in the translucent cloud core (up to AV ≈ 7 mag; t ≤ −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 Eb, m is higher than the desorption energy ED, 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 the surface and subsurface mantles (sublayers 1, 2 and 3). The limited time-span for the RD 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 the average HDO/H2O in ice does not exceed 0.2 per cent at any time when a considerable ice mass (≥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 the HDO/H2O abundance ratio in ice is 2 per cent (Dartois et al. 2003; Parise et al. 2003). RD for water ice (bulk) depends on the availability of atomic D for reactions on grain surfaces during ice accumulation in Phase 1. Figure 3. Open in new tabDownload slide Calculated evolution of subsurface bulk-ice (black) and mantle surface (grey) D/H ratios for important hydrogenated species. The bulk ice contains >90 per cent of these molecules until evaporation. Figure 3. Open in new tabDownload slide Calculated evolution of subsurface bulk-ice (black) and mantle surface (grey) D/H ratios for important hydrogenated species. The bulk ice contains >90 per cent of these molecules until evaporation. Deuterium fractionation of surface molecules in the translucent core is affected by two competing effects. First, the photodissociation of HD and D2 occurs at a higher rate than that of H2 because of self-shielding effects. This ensures an elevated atomic D/H ratio, facilitating deuteration (especially visible in the early stages in Figs 3 and 4; see also Section 2.2). Secondly, 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 RD. Both effects were initially found by Furuya et al. (2015). Figure 4. Open in new tabDownload slide Calculated atomic D and HD abundances, and the gas-phase ratio for atomic D with respect to atomic H. The evolution of the total relative abundance of all icy species is also shown. One ice ML on the grain surface corresponds to a total ice species abundance of 2 × 10−5. Figure 4. Open in new tabDownload slide Calculated atomic D and HD abundances, and the gas-phase ratio for atomic D with respect to atomic H. The evolution of the total relative abundance of all icy species is also shown. One ice ML on the grain surface corresponds to a total ice species abundance of 2 × 10−5. The dip in the ice deuteration curves, visible for all species depicted in Fig. 3 at t ≈ −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 go 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 et al. 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 RD, 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 RD values are our main result, we also present calculated abundances. This allows evaluating whether the RD values are truly relevant. We compare the numerical simulation results against observed gas-phase molecular RD 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 et al. 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 RD 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 RD in Table 4 is only authors’ estimate to help classify literature data for the purposes of this paper. 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?a . 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) 17b 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 ≈15 (20) 4.6 CH2DOCH3/CH3OCH3 – ≈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 Isotopologues . Assumed T . Abundance . References . RD . . range (K) . ratio (per cent) . . agreement?a . 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) 17b 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 ≈15 (20) 4.6 CH2DOCH3/CH3OCH3 – ≈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 References: (1) Liu et al. (2011); Persson, Jørgensen & van Dishoeck (2013); (2) Coutens et al. (2013); (3) Taquet et al. (2013a); Persson et al. (2014); (4) Coutens et al. (2014); (5) Loinard et al. (2001); (6) Shah & Wootten (2001); (7) Hatchell (2003); (8) van der Tak et al. (2002); (9) Emprechtinger et al. (2009); Tobin et al. (2013); (10) Roberts & Millar (2007); (11) Jørgensen, Schöier & van Dishoeck (2004); (12) van Dishoeck et al. (1995); (13) Schöier et al. (2002); (14) Roberts et al. (2002); (15) Loinard et al. (2000); (16) Parise et al. (2006); (17) Ceccarelli et al. (2001); (18) Parise et al. (2002); (19) Parise et al. (2004); (20) Demyk et al. (2010); (21) Richard et al. (2013); (22) Sakai et al. (2009); (23) Cordiner et al. (2012). aAgreement between observations and calculations. If no, the closest calculated value is indicated. bSee Section 3.3. Open in new tab 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?a . 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) 17b 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 ≈15 (20) 4.6 CH2DOCH3/CH3OCH3 – ≈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 Isotopologues . Assumed T . Abundance . References . RD . . range (K) . ratio (per cent) . . agreement?a . 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) 17b 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 ≈15 (20) 4.6 CH2DOCH3/CH3OCH3 – ≈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 References: (1) Liu et al. (2011); Persson, Jørgensen & van Dishoeck (2013); (2) Coutens et al. (2013); (3) Taquet et al. (2013a); Persson et al. (2014); (4) Coutens et al. (2014); (5) Loinard et al. (2001); (6) Shah & Wootten (2001); (7) Hatchell (2003); (8) van der Tak et al. (2002); (9) Emprechtinger et al. (2009); Tobin et al. (2013); (10) Roberts & Millar (2007); (11) Jørgensen, Schöier & van Dishoeck (2004); (12) van Dishoeck et al. (1995); (13) Schöier et al. (2002); (14) Roberts et al. (2002); (15) Loinard et al. (2000); (16) Parise et al. (2006); (17) Ceccarelli et al. (2001); (18) Parise et al. (2002); (19) Parise et al. (2004); (20) Demyk et al. (2010); (21) Richard et al. (2013); (22) Sakai et al. (2009); (23) Cordiner et al. (2012). aAgreement between observations and calculations. If no, the closest calculated value is indicated. bSee Section 3.3. Open in new tab Fig. 5 shows that for gas-phase water and ammonia, there are several peaks in the deuteriumfractionation 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. Figure 5. Open in new tabDownload slide Evolution of calculated relative abundances [X] (left-hand panels) and RD ratios (right-hand panels) for water and ammonia in the circumstellar envelope. Solid and dashed lines in the left-hand panels indicate gas and ice-phase species, respectively. Horizontal dotted lines in the right-hand panels indicate observational constraints (upper and lower RD limits) from Table 4. Figure 5. Open in new tabDownload slide Evolution of calculated relative abundances [X] (left-hand panels) and RD ratios (right-hand panels) for water and ammonia in the circumstellar envelope. Solid and dashed lines in the left-hand panels indicate gas and ice-phase species, respectively. Horizontal dotted lines in the right-hand panels indicate observational constraints (upper and lower RD limits) from Table 4. The second peak at 0.14 Myr (22 K) coincides with the evaporation of solid N2 and CO. These species reappear in the gas phase and are partially converted into 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 the 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 into CH|$_3^+$| by radicals and ions in the gas. Most of these are generated by the 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 RD peak, induced by methane evaporation, is about 50 kyr long and ends only when most of CH4 has been converted into CO. This third RD 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 enables additional gas-phase deuterium fractionation reaction chains. The final gas-phase RD 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, RD 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 RD peaks can be discerned for most of the other species. Appendix Fig. A1 (available in electronic form) 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 the modelling of deuterium chemistry in circumstellar envelopes (Aikawa et al. 2012; Awad et al. 2014; Taquet et al. 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 ≈10−11 at the beginning of Phase 2, thanks to the reactions proposed by Roberts et al. (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 RD 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 RD for methanol and formaldehyde deuterated isotopologues. Data for other COMs are visualized in Appendix Fig. A2 (available in electronic form). For CHD2OH, CD3OH and DCOOCH3, the high observed RD values are not reached. In the case of CH2DOH, the lowest RD 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 magnitude lower than its maximum calculated abundance of 1.4 × 10−7 during evaporation at 124 K. At these later stages, the RD for CH2DOH is only about 0.1–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. Figure 6. Open in new tabDownload slide Evolution of calculated relative abundances (left-hand panels) and RD ratios (right-hand panels) for methanol and formaldehyde, as in Fig. 5. The abundance of CH3OD (not shown for clarity) is, on average, about half that of CH2DOH. Figure 6. Open in new tabDownload slide Evolution of calculated relative abundances (left-hand panels) and RD ratios (right-hand panels) for methanol and formaldehyde, as in Fig. 5. The abundance of CH3OD (not shown for clarity) is, on average, about half that of CH2DOH. 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 2015b,c). 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 us 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 × 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 et al. 2014). Similarly to the case of D2O, the high efficiency of formaldehyde deuteration is due 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 Appendix Fig. A3 (available in electronic form)]. 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), which can be formed in substantial amounts in subsurface ice layers rich in CO, N2 or NH3 (Appendix Fig. A1; see also Kalvāns 2015b). 4 CONCLUSIONS The scope and results of the simulation are broadly similar to those reported by Taquet et al. (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 RD ratios for some COMs, contrary to a suggestion by Taquet et al. (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 pre-stellar Phase 1, atomic D is largely lost to the 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 RD. The sequential evaporation of icy molecules in a warming protostellar envelope (Phase 2) affects 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 a high RD level, which, in turn, is transferred to other species in reactions involving chemical radicals and ions. This results in species’ RD that is rapidly changing in the modelled temperature range of a protostellar envelope. The evolution of molecular RD 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, RD can be strongly affected by the evaporation of other species that take part in hydrogenation reaction chains (Section 3.3). We are grateful to the anonymous referees for many valuable comments and suggestions that made this paper better. JK, IS and JRK thank the Ventspils City Council for its persistent support. This research has made use of NASA's Astrophysics Data System. REFERENCES Aikawa Y. , Herbst E., 1999 , ApJ , 526 , 314 Crossref Search ADS Aikawa Y. , Wakelam V., Hersant F., Garrod R. T., Herbst E., 2012 , ApJ , 760 , 40 Crossref Search ADS Albertsson T. , Semenov D. A., Vasyunin A. I., Henning T., Herbst E., 2013 , ApJS , 207 , 27 Crossref Search ADS Albertsson T. , Semenov D., Henning T., 2014 , ApJ , 784 , 39 Crossref Search ADS Andersson S. , van Dishoeck E. F., 2008 , A&A , 491 , 907 Crossref Search ADS Arasa C. , Koning J., Kroes G.-J., Walsh C., van Dishoeck E. F., 2015 , A&A , 575 , A121 Crossref Search ADS Awad Z. , Viti S., Bayet E., Caselli P., 2014 , MNRAS , 443 , 275 Crossref Search ADS Carney M. T. , Yıldız U. A., Mottram J. C., van Dishoeck E. F., Ramchandani J., Jørgensen J. K., 2016 , A&A , 586 , A44 Crossref Search ADS Caselli P. , Stantcheva T., Shalabiea O., Shematovich V. I., Herbst E., 2002 , Planet. Space Sci. , 50 , 1257 Crossref Search ADS Cazaux S. , Spaans M., 2004 , ApJ , 611 , 40 Crossref Search ADS Ceccarelli C. , Loinard L., Castets A., Tielens A. G. G. M., Caux E., 2000 , A&A , 357 , L9 Ceccarelli C. , Loinard L., Castets A., Tielens A. G. G. M., Caux E., Lefloch B., Vastel C., 2001 , A&A , 372 , 998 Crossref Search ADS Ceccarelli C. , Caselli P., Bockelée-Morvan D., Mousis O., Pizzarello S., Robert F., Semenov D., 2014 , Beuther H., Klessen R. S., Dullemond C. P., Henning T., , Protostars and Planets VI , Univ. Arizona Press , Tucson, AZ , 859 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chang Q. , Herbst E., 2014 , ApJ , 787 , 135 Crossref Search ADS Chehrouri M. et al. , 2011 , Phys. Chem. Chem. Phys. , 13 , 2172 Crossref Search ADS PubMed Cordiner M. A. , Charnley S. B., Wirström E. S., Smith R. G., 2012 , ApJ , 744 , 131 Crossref Search ADS Coutens A. et al. , 2013 , A&A , 560 , A39 Crossref Search ADS Coutens A. , Jørgensen J. K., Persson M. V., van Dishoeck E. F., Vastel C., Taquet V., 2014 , ApJ , 792 , L5 Crossref Search ADS Dartois E. , Thi W.-F., Geballe T. R., Deboffle D., d'Hendecourt L., van Dishoeck E., 2003 , A&A , 399 , 1009 Crossref Search ADS Demyk K. , Bottinelli S., Caux E., Vastel C., Ceccarelli C., Kahane C., Castets A., 2010 , A&A , 517 , A17 Crossref Search ADS Emprechtinger M. , Caselli P., Volgenau N. H., Stutzki J., Wiedner M. C., 2009 , A&A , 493 , 89 Crossref Search ADS Evans II N. J. et al. , 2009 , ApJS , 181 , 321 Crossref Search ADS Fayolle E. C. , Öberg K. I., Cuppen H. M., Visser R., Linnartz H., 2011 , A&A , 529 , A74 Crossref Search ADS Flower D. R. , Pineau des Forêts G., Walmsley C. M., 2004 , A&A , 427 , 887 Crossref Search ADS Flower D. R. , Pineau Des Forêts G., Walmsley C. M., 2006a , A&A , 449 , 621 Crossref Search ADS Flower D. R. , Pineau Des Forêts G., Walmsley C. M., 2006b , A&A , 456 , 215 Crossref Search ADS Furuya K. , Aikawa Y., Nomura H., Hersant F., Wakelam V., 2013 , ApJ , 779 , 11 Crossref Search ADS Furuya K. , Aikawa Y., Hincelin U., Hassel G. E., Bergin E. A., Vasyunin A. I., Herbst E., 2015 , A&A , 584 , A124 Crossref Search ADS Garrod R. T. , 2013 , ApJ , 765 , 60 Crossref Search ADS Garrod R. T. , Herbst E., 2006 , A&A , 457 , 927 Crossref Search ADS Garrod R. T. , Pauly T., 2011 , ApJ , 735 , 15 Crossref Search ADS Garrod R. T. , Weaver S. L. W., Herbst E., 2008 , ApJ , 682 , 283 Crossref Search ADS Hama T. , Kuwahata K., Watanabe N., Kouchi A., Kimura Y., Chigai T., Pirronello V., 2012 , ApJ , 757 , 185 Crossref Search ADS Hartmann L. , 2009 , Accretion Processes in Star Formation , 2nd edn, Cambridge Univ. Press , Cambridge Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hatchell J. , 2003 , A&A , 403 , L25 Crossref Search ADS Heine T. , Zhechkov L., Seifert G., 2004 , Phys. Chem. Chem. Phys. , 6 , 980 Crossref Search ADS Hocuk S. , Cazaux S., 2015 , A&A , 576 , A49 Crossref Search ADS Hocuk S. , Cazaux S., Spaans M., Caselli P., 2016 , MNRAS , 456 , 2586 Crossref Search ADS Jørgensen J. K. , Schöier F. L., van Dishoeck E. F., 2004 , A&A , 416 , 603 Crossref Search ADS Kalvaāns J. , 2015a , A&A , 573 , A38 Crossref Search ADS Kalvāns J. , 2015b , ApJ , 806 , 196 Crossref Search ADS Kalvāns J. , 2015c , ApJ , 803 , 52 Crossref Search ADS Kalvāns J. , Shmeld I., 2010 , A&A , 521 , A37 Crossref Search ADS Kalvāns J. , Shmeld I., 2013 , A&A , 554 , A111 Crossref Search ADS Katz N. , Furman I., Biham O., Pirronello V., Vidali G., 1999 , ApJ , 522 , 305 Crossref Search ADS Keto E. , Caselli P., 2010 , MNRAS , 402 , 1625 Crossref Search ADS Kristensen L. E. , Amiaud L., Fillion J.-H., Dulieu F., Lemaire J.-L., 2011 , A&A , 527 , A44 Crossref Search ADS Le Bourlot J. , 2000 , A&A , 360 , 656 Le Petit F. , Roueff E., Le Bourlot J., 2002 , A&A , 390 , 369 Crossref Search ADS Linsky J. L. et al. , 2006 , ApJ , 647 , 1106 Crossref Search ADS Liu F.-C. , Parise B., Kristensen L., Visser R., van Dishoeck E. F., Güsten R., 2011 , A&A , 527 , A19 Crossref Search ADS Loinard L. , Castets A., Ceccarelli C., Tielens A. G. G. M., Faure A., Caux E., Duvert G., 2000 , A&A , 359 , 1169 Loinard L. , Castets A., Ceccarelli C., Caux E., Tielens A. G. G. M., 2001 , ApJ , 552 , L163 Crossref Search ADS Maret S. et al. , 2004 , A&A , 416 , 577 Crossref Search ADS Martín-Doménech R. , Muñoz Caro G. M., Bueno J., Goesmann F., 2014 , A&A , 564 , A8 Crossref Search ADS Maury A. J. , André P., Men'shchikov A., Könyves V., Bontemps S., 2011 , A&A , 535 , A77 Crossref Search ADS Millar T. J. , Bennett A., Herbst E., 1989 , ApJ , 340 , 906 Crossref Search ADS Mispelaer F. et al. , 2013 , A&A , 555 , A13 Crossref Search ADS Nishimura Y. , Sakai N., Watanabe Y., Sakai T., Hirota T., Yamamoto S., 2013 , Kawabe R., Kuno N., Yamamoto S., , ASP Conf. Ser. Vol. 476, New Trends in Radio Astronomy in the ALMA Era: The 30th Anniversary of Nobeyama Radio Observatory , Astron. Soc. Pac. , San Francisco , 343 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pagani L. et al. , 2009 , A&A , 494 , 623 Crossref Search ADS Parise B. et al. , 2002 , A&A , 393 , L49 Crossref Search ADS Parise B. , Simon T., Caux E., Dartois E., Ceccarelli C., Rayner J., Tielens A. G. G. M., 2003 , A&A , 410 , 897 Crossref Search ADS Parise B. , Castets A., Herbst E., Caux E., Ceccarelli C., Mukhopadhyay I., Tielens A. G. G. M., 2004 , A&A , 416 , 159 Crossref Search ADS Parise B. , Ceccarelli C., Tielens A. G. G. M., Castets A., Caux E., Lefloch B., Maret S., 2006 , A&A , 453 , 949 Crossref Search ADS Persson M. V. , Jørgensen J. K., van Dishoeck E. F., 2013 , A&A , 549 , L3 Crossref Search ADS Persson M. V. , Jørgensen J. K., van Dishoeck E. F., Harsono D., 2014 , A&A , 563 , A74 Crossref Search ADS Prodanović T. , Steigman G., Fields B. D., 2010 , MNRAS , 406 , 1108 Richard C. et al. , 2013 , A&A , 552 , A117 Crossref Search ADS Roberts H. , Millar T. J., 2000 , A&A , 361 , 388 Roberts H. , Millar T. J., 2007 , A&A , 471 , 849 Crossref Search ADS Roberts H. , Fuller G. A., Millar T. J., Hatchell J., Buckle J. V., 2002 , A&A , 381 , 1026 Crossref Search ADS Roberts H. , Herbst E., Millar T. J., 2004 , A&A , 424 , 905 Crossref Search ADS Sakai N. , Sakai T., Hirota T., Yamamoto S., 2009 , ApJ , 702 , 1025 Crossref Search ADS Schnee S. , Di Francesco J., Enoch M., Friesen R., Johnstone D., Sadavoy S., 2012 , ApJ , 745 , 18 Crossref Search ADS Schöier F. L. , Jørgensen J. K., van Dishoeck E. F., Blake G. A., 2002 , A&A , 390 , 1001 Crossref Search ADS Shah R. Y. , Wootten A., 2001 , ApJ , 554 , 933 Crossref Search ADS Sipilä O. , Caselli P., Harju J., 2013 , A&A , 554 , A92 Crossref Search ADS Sipilä O. , Caselli P., Harju J., 2015 , A&A , 578 , A55 Crossref Search ADS Stahler S. W. , Palla F., 2005 , The Formation of Stars , Wiley , New York Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sugimoto T. , Fukutani K., 2011 , Nature Phys. , 7 , 307 Crossref Search ADS Taquet V. , Ceccarelli C., Kahane C., 2012 , A&A , 538 , A42 Crossref Search ADS Taquet V. , López-Sepulcre A., Ceccarelli C., Neri R., Kahane C., Coutens A., Vastel C., 2013a , ApJ , 768 , L29 Crossref Search ADS Taquet V. , Peters P. S., Kahane C., Ceccarelli C., López-Sepulcre A., Toubin C., Duflot D., Wiesenfeld L., 2013b , A&A , 550 , A127 Crossref Search ADS Taquet V. , Charnley S. B., Sipilä O., 2014 , ApJ , 791 , 1 Crossref Search ADS Tobin J. J. et al. , 2013 , ApJ , 765 , 18 Crossref Search ADS Valencic L. A. , Smith R. K., 2015 , ApJ , 809 , 66 Crossref Search ADS van der Tak F. F. S. , Schilke P., Müller H. S. P., Lis D. C., Phillips T. G., Gerin M., Roueff E., 2002 , A&A , 388 , L53 Crossref Search ADS van Dishoeck E. F. , Blake G. A., Jansen D. J., Groesbeck T. D., 1995 , ApJ , 447 , 760 Crossref Search ADS van Dishoeck E. F. , Herbst E., Neufeld D. A., 2013 , Chem. Rev. , 113 , 9043 Crossref Search ADS PubMed Visser R. , Doty S. D., van Dishoeck E. F., 2011 , A&A , 534 , A132 Crossref Search ADS Wakelam V. , Vastel C., Aikawa Y., Coutens A., Bottinelli S., Caux E., 2014 , MNRAS , 445 , 2854 Crossref Search ADS Walmsley C. M. , Flower D. R., Pineau des Forêts G., 2004 , A&A , 418 , 1035 Crossref Search ADS Watanabe N. , Kimura Y., Kouchi A., Chigai T., Hama T., Pirronello V., 2010 , ApJ , 714 , L233 Crossref Search ADS Whittet D. C. B. , Gerakines P. A., Hough J. H., Shenoy S. S., 2001 , ApJ , 547 , 872 Crossref Search ADS Willacy K. , Woods P. M., 2009 , ApJ , 703 , 479 Crossref Search ADS Wolcott-Green J. , Haiman Z., 2011 , MNRAS , 412 , 2603 Crossref Search ADS Öberg K. I. , Fayolle E. C., Cuppen H. M., van Dishoeck E. F., Linnartz H., 2009c , A&A , 505 , 183 Crossref Search ADS APPENDIX A: CALCULATED ABUNDANCES AND FRACTIONATION RATIOS FOR DEUTERATED SPECIES OBSERVED AROUND YOUNG STELLAR OBJECTS Figure A1. Open in new tabDownload slide Evolution of calculated relative abundances [X] (left) and RD ratios (right) for hydrogen isotopologues of observed gaseous (in)organic species in the circumstellar envelope. Solid and dashed lines in the left-hand panels indicate gas and ice-phase species, respectively. Horizontal dotted lines in the right-hand panels indicate observational constraints (upper and lower RD limits) from Table 4. Figure A1. Open in new tabDownload slide Evolution of calculated relative abundances [X] (left) and RD ratios (right) for hydrogen isotopologues of observed gaseous (in)organic species in the circumstellar envelope. Solid and dashed lines in the left-hand panels indicate gas and ice-phase species, respectively. Horizontal dotted lines in the right-hand panels indicate observational constraints (upper and lower RD limits) from Table 4. Figure A2. Open in new tabDownload slide Evolution of calculated relative abundances (left) and RD ratios (right) for methane, methyl formate, and dimethyl ether as in Fig. A1. Figure A2. Open in new tabDownload slide Evolution of calculated relative abundances (left) and RD ratios (right) for methane, methyl formate, and dimethyl ether as in Fig. A1. Figure A3. Open in new tabDownload slide Evolution of calculated relative abundances (left) and RD ratios (right) for observed carbon-chain species and acetylene as in Fig. A1. Figure A3. Open in new tabDownload slide Evolution of calculated relative abundances (left) and RD ratios (right) for observed carbon-chain species and acetylene as in Fig. A1. © 2017 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Chemical fractionation of deuterium in the protosolar nebula
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stx174
DA - 2017-05-21
UR - https://www.deepdyve.com/lp/oxford-university-press/chemical-fractionation-of-deuterium-in-the-protosolar-nebula-nQAXoPC0x5
SP - 1763
VL - 467
IS - 2
DP - DeepDyve
ER -