EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 584 (December 2015) A&A, 584 (2015) A124 Full HTML
Free Access
Issue
A&A
Volume 584, December 2015
Article Number A124
Number of page(s) 18
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201527050
Published online 04 December 2015

© ESO, 2015

1. Introduction

It is well established that water is predominantly formed via grain surface reactions, and is already abundant (10-4 with respect to hydrogen nuclei) in molecular clouds with line of sight visual extinction ≳3 mag (Whittet 1993, 2003). One of the major unresolved questions in the field of astrochemistry is whether, and how much, water in the solar system originated from the parent molecular cloud (e.g., recent reviews by Ceccarelli et al. 2014; van Dishoeck et al. 2014). The HDO / H2O ratio can be a useful tool to probe the processing of water during star- and planet-formation. We can trace water evolution by observing and comparing the HDO / H2O ratio in objects at different evolutionary stages. The upper limits of the HDO / H2O ice ratio in embedded protostars have been derived as (2−10) × 10-3 depending on the source (Dartois et al. 2003; Parise et al. 2003), which are much larger than the elemental abundance of deuterium (1.5 × 10-5) in the local interstellar medium (ISM, Linsky 2003). Recent interferometric observations have revealed the deuteration ratios of water vapor in the inner hot regions (T> 100 K) of the embedded protostellar sources, where water ice is sublimated (Jørgensen & van Dishoeck 2010; Coutens et al. 2013, 2014a; Persson et al. 2013, 2014; Taquet et al. 2013a). It has been found that the HDO / H2O ratio is as high as 10-3, which is similar to that in Oort cloud comets (e.g., Villanueva et al. 2009; Bockelée-Morvan et al. 2012). This similarity may imply that some of the cometary water originated from the embedded protostellar phase or earlier, while some processing may have been at work in the disk, taking the variation of the HDO / H2O ratio in the solar system objects into consideration (see e.g., recent theoretical work by Willacy & Woods 2009; Visser et al. 2011; Yang et al. 2013; Furuya et al. 2013; Albertsson et al. 2014; Cleeves et al. 2014).

On the other hand, to the best of our knowledge, there has been no detection of HDO ice in molecular clouds and prestellar cores. This limits our understanding of water evolution from molecular clouds to the embedded sources. Instead, there have been numerous theoretical studies on the HDO / H2O ice ratio in cold (~10 K) and dense (>104 cm-3) conditions (e.g., recent work by Taquet et al. 2013b; Sipilä et al. 2013; Lee & Bergin 2015). When the primary reservoir of deuterium is HD, deuterium fractionation is driven by isotope exchange reactions, such as (1)The reaction in the backward direction is endothermic, and thus prohibited at low temperatures of ≲20 K. The enrichment of deuterium in e.g., H is transfered to both gaseous and icy molecules through sequential chemical reactions (Millar et al. 1989; Roberts et al. 2004). However, the presence of ortho-H2 (o - H2) can change the situation. It can suppress the overall deuteration processes even at T ≲ 20 K, since the internal energy of o - H2, which is 170.5 K above that of para-H2 (p - H2), helps to overcome the endothermicity of the exchange reactions in the backward direction (Pagani et al. 1992; Gerlich et al. 2002; Walmsley et al. 2004). CO and electrons also play a role in deuterium chemistry; they destroy H2D+ to suppress further fractionation (Roberts et al. 2002). Figure 1 shows the threshold abundances of CO and electrons with respect to H2 as functions of the ortho-to-para nuclear spin ratio of H2 (OPR(H2)); below the lines for CO and electrons in this figure, Reaction (1) in the backward direction dominates the destruction (see Appendix A for more details). For example, the abundances of CO and electrons in dense cloud conditions are typically <10-4 and ~10-8−10-7, respectively. Then, even a small fraction of o - H2, ~10-3, can affect the deuterium chemistry at low temperatures.

thumbnail Fig. 1

Threshold abundances of CO and electrons to be the dominant reactant of H2D+. If CO and electron abundances with respect to H2 are below the lines, Reaction (1) in the backward direction dominates over H2D+ destruction by CO and electrons. See Appendix A for more details.

In astrochemical simulations, it has been generally assumed that hydrogen (and deuterium) is already locked in H2 (and HD) at the beginning of the simulations, and the focus is placed on the subsequent molecular evolution after H2 (and HD) formation. In this approach, the initial OPR(H2) is treated as a free parameter, since the OPR(H2) in the ISM is not well-constrained; it requires detailed physical and chemical modeling to estimate the OPR(H2) from observations of molecules other than H2 in cold clouds/cores (e.g., Pagani et al. 2009, 2011; Le Gal et al. 2014). The resultant deuteration ratio of icy molecules in the simulations strongly depends on the initial OPR(H2), since the time it takes for the OPR(H2) to reach steady state is longer (several Myr, Flower et al. 2006) than the timescale of ice formation (~the freeze-out timescale) in dense cloud core conditions.

In this paper, we investigate the evolution of water deuteration and the OPR(H2) during the chemical evolution from H i-dominated clouds through the formation of denser molecular clouds. One of the plausible scenarios of molecular cloud formation is that diffuse H i gas is swept up by global accretion flows, such as supernova explosions, via expanding H ii regions, or via colliding turbulent flows in the diffuse ISM (e.g., Hennebelle & Pérault 1999; Hartmann et al. 2001; Inoue & Inutsuka 2012; Inutsuka et al. 2015). A recent review of molecular cloud formation can be found in Dobbs et al. (2014). To simulate this process, we revisit a one-dimensional shock model developed by Bergin et al. (2004, hereafter B04) and Hassel et al. (2010, hereafter H10) with post-processing gas-ice chemistry calculations. This allows us to study OPR(H2) and water deuteration during the formation and early evolution of molecular clouds without making an arbitrary assumption of the initial molecular abundances, including the initial OPR(H2). B04 studied H2 and CO formation, while H10 focused on the formation of icy molecules without deuterium and nuclear spin-state chemistry.

This paper is organized as follows: we describe our physical model and chemical model in Sects. 2 and 3, respectively. In Sect. 4 we present the evolution of the OPR(H2) and water ice deuteration in our fiducial model. Parameter dependences are discussed in Sect. 5, and comparisons between our and previous studies are made in Sect. 6. Our findings are summarized in Sect. 7. In the Appendix, we derive some analytical formulae for the OPR(H2), its effect on deuterium chemistry, and deuterium fractionation of water ice under irradiation conditions, which may help readers to understand the numerical results and the dependencies on adopted parameters.

2. Physical model

To simulate the formation and physical evolution of a molecular cloud, we use the one-dimensional shock model developed by B04 and H10. Here we briefly outline the model, while more details can be found in the original papers.

The model describes the evolution of post-shock materials in a plane-parallel configuration. In order to obtain density and gas temperature evolution, the conservation laws of mass and momentum are solved with the energy equation, considering time-dependent cooling/heating rates and simplified chemistry (B04). The cosmic-ray ionization rate of H2 (ξH2) is set to be 1.3 × 10-17 s-1. Interstellar radiation is assumed to be incident from one side only, adopting the Draine field (Draine 1978). As time proceeds, the column density of post-shock materials increases, which assists molecular formation by absorbing interstellar radiation. The column density at a given time t after passing through the shock front is NH = n0v0t, where n0 and v0 are preshock H i gas density and velocity of the accretion flow, respectively. NH is converted into AV using the relation AV/NH = 5 × 10-22 mag/cm-2. In this work we choose the model with n0 = 10 cm-3 and v0 = 15 km s-1 (model 4 in H10). The preshock gas temperature is assumed to be 40 K. The preshock density and temperature agree with observed properties of H i gas in the ISM (Heiles & Troland 2003). The magnetic field in H i clouds is of the order of μG (Heiles & Troland 2005). In the scenario of molecular cloud formation via accretion of H i gas, molecular clouds are formed in the regions where the angle between the shock normal and the mean magnetic field is sufficiently low; otherwise the magnetic pressure prevents the accumulation of the gas (Inoue & Inutsuka 2009, 2012). In our physical model, the magnetic field component parallel to the shock front is set to be 0.01 μG, assuming that the shock normal and the mean magnetic field are almost parallel. In this setting, the formation timescale of molecular clouds with AV = 1 mag is ~4 Myr (AV/1 mag)(n0/10 cm-3)-1 (v0/15 km s-1)-1, which is comparable to or less than the estimated lifetime of molecular clouds (Hartmann et al. 2001; Kawamura et al. 2009). The impact of varying n0 and v0 is briefly discussed in Sect. 5.4.

Using density, gas temperature, and AV profiles, the dust temperature is calculated by equating the heating rate with the cooling rate. Heating processes include heating by interstellar radiation and collisions with gas molecules, while cooling processes include thermal emission and sublimation of ices (H10). Figure 2 shows the evolution of gas density, gas temperature (Tgas), and dust temperature (Tdust) adopted in this work.

In the physical model, the effect of self-gravity is not considered. Following the discussion by Hartmann et al. (2001), internal pressure due to the self gravity of a sheet-like cloud exceeds ram pressure due to the accretion flow when AV ≳ 4.5 mag (Pe/ 3.2 × 105kB)0.5, where Pe and kB are the ram pressure of an accretion flow and the Boltzmann constant, respectively. We use the shock model until AV reaches 3 mag, so that our simulation is in the regime where the ram pressure overwhelms the pressure due to self gravity.

thumbnail Fig. 2

Physical evolution after passing through the shock front. We adopt the same shock model as Model 4 of B04 and H10; i.e., the preshock density, temperature and velocity are 10 cm-3, 40 K, and 15 km s-1 respectively. The inset shows the zoom-in view of the earliest evolution at AV ≤ 0.1 mag, which corresponds to t ≤ 0.3 Myr.

3. Chemical model

To simulate the gas-ice chemistry during the formation and evolution of the cloud, we use a rate equation method adopting a three-phase model, which consists of gas, a chemically active icy surface, and an inactive inert ice mantle (Hasegawa & Herbst 1993b). In three-phase models, it is often assumed that chemical processes occur only in the uppermost monolayer (e.g., Hasegawa & Herbst 1993b; Garrod & Pauly 2011). However, the exact number of chemically active ice surface layers (Nact) remains uncertain. Laboratory experiments have shown that atomic hydrogen can penetrate into the top several monolayers of ice, and reactions can occur in the ice as well as on the surface (Fuchs et al. 2009; Ioppolo et al. 2010). Vasyunin & Herbst (2013a) find in their Monte-Carlo simulations that a model with Nact = 4 can reproduce the experimental thermal desorption profiles of mixed ices better than a model with Nact = 1. We set Nact to be 4 in our fiducial model, so that the top four monolayers are chemically active and have uniform chemical composition. In the rest of this paper, we refer to the outermost Nact layers as the active (surface) ice layers, and all of the layers including the active ice layers as the bulk ice mantle.

As chemical processes, we take into account gas-phase reactions, interactions between gas and (icy) grain surfaces, and surface reactions. For photodissociation and photoionization in the gas phase, we consider the self-shielding of H2, HD, CO, C i, and N2 and mutual shielding by H2 of HD, CO, C i, and N2 (Draine & Bertoldi 1998; Kamp & Bertoldi 2000; Visser et al. 2009; Wolcott-Green & Haaiman 2011; Li et al. 2013). Grain surface reactions are treated using the modified-rate equation method of Garrod (2008, method A). It can take into account the competition between surface processes to improve the accuracy of the original rate equation method in the stochastic regime. Our chemical reaction network is originally based on the gas-ice reaction network of Garrod & Herbst (2006). When the gas temperature exceeds 100 K, we use the reaction sets in Harada et al. (2010, 2012) instead of those in Garrod & Herbst (2006) for the neutral-neutral reactions and neutral-ion reactions; the network of Harada et al. (2010, 2012) was designed for modeling high-temperature gas-phase chemistry (100−800 K) and includes reactions with a high potential energy barrier. Our network has been extended to include mono, doubly, and triply deuterated species (Aikawa et al. 2012; Furuya et al. 2013) and nuclear spin states of H2, H, and their isotopologues (Hincelin et al. 2014; Coutens et al. 2014b). The rate coefficients for the H2 + H system, including the isotopologues HD, D2, H2D+, D2H+, and D with their nuclear spin states, are taken from Hugo et al. (2009). In the rest of this sections, we describe the treatment of interaction between gas and (icy) grain surface, surface reactions (Sect. 3.1), and the initial abundance (Sect. 3.2).

3.1. Gas-grain interaction and surface chemistry

When gaseous neutral species collide with a dust grain, they can stick to the surface. The radius of a dust grain is set to be 0.1 μm with the dust-to-gas mass ratio of 0.01. The sticking probabilities for H, D, H2, HD, and D2 are calculated as functions of gas and dust temperatures (Hollenbach & McKee 1979), while those for heavier species are set to be unity. Interactions between ions and dust grains are calculated in the same way as Furuya et al. (2012); the collisional rates between ions and grains are calculated considering Coulomb focusing (Draine & Sutin 1987), and the products and branching ratios are assumed to be the same as the corresponding electron recombination in the gas phase.

Adsorbed species can desorb thermally into the gas phase again. The rate of thermal desorption depends exponentially on the binding energy of absorbed species on a surface. The binding energy of a species depends on the molecular composition of a surface. For example, the binding energy of H2 on an H2O substrate is several hundred K (Hornekaer et al. 2005), while that of H2 on an H2 substrate is much smaller (23 K, Cuppen & Herbst 2007). In this work we use the binding energies on a water ice substrate, and modify them depending on the surface coverage of H2 (θH2, Garrod & Pauly 2011): (2)where is the binding energy of species A on a water substrate taken from Garrod & Herbst (2006). and are treated as free parameters, and we set them to be 550 K in our fiducial model. Without the modification, a three phase model can give unphysical solutions; in typical cold dark cloud conditions with K, H2 is depleted from the gas phase and thick pure H2 ice mantles form (Garrod & Pauly 2011). The desorption energy of atomic deuterium is set to be 21 K higher than that of atomic hydrogen ( K), following Caselli et al. (2002). For other deuterated species, we use the same desorption energies as for normal species.

For non-thermal desorption processes, we consider stochastic heating by cosmic-rays (Hasegawa & Herbst 1993a), photodesorption (Westley et al. 1995; Öberg et al. 2007), and reactive desorption (Garrod et al. 2007; Dulieu et al. 2013; Vasyunin & Herbst 2013b). We assume that roughly 1% of products formed by exothermic reactions are desorbed following the method of Garrod et al. (2007), except for the reactions OH + H → H2O, and O + O → O2. Recently Dulieu et al. (2013) found that desorption probabilities of these reactions are very high, >90% and 60%, on a silicate substrate, while the probabilities are much lower on an ice substrate. Minissale & Dulieu (2014) found that the desorption probability of O + O → O2 becomes smaller with the increasing ice coverage. Motivated by these experiments, we use the coverage dependent desorption probabilities for these two reactions and their deuterated analogues: (3)where θice is the surface coverage of the ice. Pbare is the desorption probability on a silicate substrate taken from Dulieu et al. (2013), while Pice is that on a ice substrate, calculated by the method of Garrod et al. (2007).

Photochemistry of ices is an essential process to prevent formation of thick ice mantles in gas with low extinction (e.g., Tielens 2005; Cuppen & Herbst 2007; Hollenbach et al. 2009). According to molecular dynamics (MD) simulations, there are several possible outcomes after a UV photon dissociates water ice; (i) the photofragments are trapped on the surface; (ii) either of the fragments is desorbed into the gas phase; (iii) the fragments recombine and the product is either trapped on the surface or desorbed into the gas phase, etc. (Andersson et al. 2006; Andersson & van Dishoeck 2008). The total photodissociation rates (cm-3 s-1) of water ice are calculated as follows (Furuya et al. 2013; Taquet et al. 2013b): where ngr is the number density of dust grains, a is the radius of dust grains, fabs,i is the fraction of the incident photons absorbed by species i (i.e., water ice here), and Nlayer is the number of monolayers of the bulk ice mantle. FUV is the flux of the FUV radiation field, assuming 2 × 108 photon cm-2 s-1 for interstellar radiation and 3 ×103 photon cm-2 s-1 for cosmic-ray induced UV radiation. Pabs,i is the absorption probability of the incident FUV photon per monolayer. We calculated Pabs,i by convolving wavelength dependent photoabsorption cross sections of water ice (Mason et al. 2006) with the emission spectrum of the interstellar radiation field (Draine 1978) and with that of the cosmic-ray-induced radiation field (Gredel et al. 1989). When we calculate fabs,i, we consider the absorption by up to four outermost monolayers regardless of Nact. The parameter γi is for the attenuation of interstellar radiation field by dust grains, and we apply the value for photodissociation in the gas phase. Recently, Arasa et al. (2015) performed MD simulations to study the possible outcomes of H2O, HDO, and D2O photodissociation in H2O ice and their probabilities per dissociation (bj) as functions of depth into the ice. With their results, the rate of each outcome is calculated by . The total photodesorption yield of water ice (desorbed as OX or X2O, X is H or D) per incident photon for thick pure water ice is ~3× 10-4 in our model.

The absorption probability of the incident FUV photon for the other icy species are calculated as , where is the unattenuated photodissociation rate in the gas phase. We also assume that all the photofragments are trapped in the active surface ice layers for simplicity. Instead, the photodesorption rates are calculated as (6)where Yi is the photodesorption yield per incident FUV photon for thick ice (Nlayer ≥ 4). The term min(Nlayer/ 4, 1) is for thin ice (Nlayer< 4), assuming that the yield is linear to ice thickness. We use the yields derived from experimental work for CO, CO2, O2, and N2 pure ices (Öberg et al. 2009; Fayolle et al. 2011, 2013). Bertin et al. (2012) found that photodesorption of CO is less efficient on a water substrate than on a CO substrate by more than a factor of 10. Infrared observations suggest that water ice is the main constituent of interstellar ice in the early stage of ice formation, while CO is the dominant constituent of the outer layers of the ice (e.g., Öberg et al. 2011). Considering the importance of CO abundance in the gas phase on deuterium chemistry, we treat the photodesorption yield of CO as a function of the surface coverage of CO: (7)where YCO - CO is the photodesorption yield for pure CO ice taken from Fayolle et al. (2011, 10-2, while YCO - H2O is the photodesorption yield for CO interacting with H2O. We assume YCO - H2O = 3 × 10-4. In this way, the CO freeze-out is self-limited in our model; with increasing the surface coverage of CO, the photodesorption yield of CO increases and the CO freeze-out is delayed. For species for which no data is available in the literature, we set Yi to be 10-3.

Surface reactions are assumed to occur by the Langmuir-Hinshelwood mechanism between physisorbed species; adsorbed species diffuse by thermal hopping and react with each other when they meet (e.g., Hasegawa et al. 1992). The reaction rate coefficient (s-1) of surface reaction, A + B → AB, is given by where khop(A) is the surface diffusion rate of species A by thermal hopping from one binding site to another, and NS is the number of binding sites on the surface per monolayer (~2 × 106). The energy barrier against hopping is set to be half of the desorption energy given in Eq. (2). The pre-exponential factor νA is the vibrational frequency of species A in the binding site and is evaluated by using the harmonic oscillator strength (Hasegawa et al. 1992). We do not allow the surface diffusion through quantum tunneling; laboratory experiments have shown that the thermal hopping is the dominant mechanism of the surface diffusion of atomic hydrogen on cold amorphous water ice (Watanabe et al. 2010; Hama et al. 2012). The efficiency factor εact is for reactions with the activation energy barriers, otherwise it is set to unity. When the reactants are in the same binding site, there would be the two possible outcomes; a reaction occurs overcoming the activation energy barrier or either one of the reactants hop to another binding site before the reaction occurs (Tielens & Hagen 1982, “reaction-diffusion competition”). This is naturally included in microscopic Monte Carlo simulations (e.g., Chang & Herbst 2014). Considering the competition, εact is evaluated as follows (Tielens & Hagen 1982; Awad et al. 2005; Chang et al. 2007): (10)where νA + B is the collisional frequency of the reactants in the binding site and κA + B is the probability to overcome the activation energy barrier per collision. Assuming the collisional frequency can be approximated by the vibrational frequency in the binding site, we set νA + B to be the largest of the vibrational frequency of either species A or B (Garrod & Pauly 2011). The barriers are assumed to be overcome either thermally or via quantum tunneling, whichever is faster. Transmission probabilities (probabilities of tunneling through the activation energy barriers) of reactions relevant for water ice are calculated assuming a rectangular potential barrier with a width of 1 Å. We treat the use of the reaction-diffusion competition as a free parameter in order to study its impact, because most published rate equation models do not consider the reaction-diffusion competition. When we do not consider the competition, Eq. (10) is replaced by as the analogy of gas phase reactions (e.g., Hasegawa et al. 1992). Figure 3 shows the efficiency factor for reaction OH + H2→ H2O + H, , with and without the competition in our models as functions of dust temperature. The activation energy barrier of the reaction is set to be 2100 K (Atkinson et al. 2004), and thus it is overcome via quantum tunneling at low temperatures (see Oba et al. 2012). When the competition is considered, is extremely higher than the transmission probability. Note that, when the competition is considered, decreases with increasing dust temperature and with decreasing Edes(H2), reflecting a increase of khop(H2).

We set the branching ratio for the H2 formation on a surface (o - H2):(p - H2) to be 3:1 as experimentally demonstrated by Watanabe et al. (2010). They also find that the spin conversion occurs on amorphous water ice in a short timescale (<1 h), while the mechanism remains uncertain. We do not consider spin conversion on a surface in this work.

thumbnail Fig. 3

Efficiency factor for reaction OH + H2→ H2O + H, , without (yellow) and with the reaction-diffusion competition in our models as functions of dust temperature. For the case with the competition, with different Edes(H2) are shown by different colors.

3.2. Initial abundances and parameters

The amount of material in the ISM that is available for gas and ice chemistry remains largely uncertain (e.g., Sofia 2004). B04 and H10 used the elemental abundances based on the observations of the ζ Oph diffuse cloud. We refer to elements heavier than oxygen as heavy metals in this paper. Graedel et al. (1982) found that when they adopt much lower abundances of heavy metals than the ζ Oph values, their gas-phase astrochemical model better reproduces the observation of various molecules in dense molecular clouds. The elemental abundances similar to the ζ Oph diffuse cloud are referred to “high metal (HM)”, while the reduced abundances are referred to “low metal (LM)” (Graedel et al. 1982). Following Graedel et al., the LM abundances have been widely used in astrochemical models of dense molecular clouds. The success of the LM abundances implies that the depletion of heavy metals from the gas phase occurs during the evolution from diffuse clouds to denser clouds, which is supported by some observational evidence (Joseph et al. 1986). Plausible mechanisms of the depletion are thought to be entrapment in ices (Andersson et al. 2013; Kama et al. 2015) and incorporation into refractory dust grains (Keller et al. 2002). The former process is naturally included in our simulations, while the latter is not.

To study the impact of different initial abundances, we consider both the HM abundances and the LM abundances, which are summarized in Table 1. All the elements are initially assumed to be in the form of either neutral atoms or atomic ions, depending on their ionization energy. We use the LM abundances in our fiducial model.

In Table 2, we summarize free chemical parameters we consider in this paper. The impact of the parameters is discussed in Sect. 5.

Table 1

Initial abundances with respect to hydrogen nuclei.

thumbnail Fig. 4

Fractional abundances of selected species with respect to hydrogen nuclei in the fiducial model as functions of visual extinction, or time. The solid lines represent species in the bulk ice mantle, while the dashed lines represent gaseous species. The H2 and H atom abundances are multiplied by a factor of 10-3.

Table 2

Summary of free chemical parameters.

4. Results from the fiducial model

4.1. Overview of molecular formation

Figure 4 shows the fractional abundances of selected species as functions of visual extinction. The horizontal axis is also read as time after passing through the shock front (see the labels at the top). Note that AV in our model is likely different from the line of site visual extinction estimated from observations; in reality, molecular clouds can be highly structured and UV photons can penetrate into the clouds along a direction where extinction is lower than the line of site visual extinction, while in our model interstellar radiation is assumed to be incident from one side only following the plane of the shock.

In the shock front, the gas is heated up to 7000 K. As the gas is cooled via the line emission of O i and C ii and the gas density increases from 10 cm-3 to 104 cm-3, the conversion of atomic hydrogen to molecular hydrogen occurs via grain surface reactions, supported by the self-shielding of UV radiation. Following H2 formation, other self-shielding species with smaller abundances, CO and HD, become abundant (B04). Another self-shielding species N2 becomes moderately abundant later (AV ~ 1.5 mag), though N2 is not the primary reservoir of nitrogen during the simulation, in agreement with previous studies (Maret et al. 2006; Daranlot et al. 2012). At AV ≲ 1 mag the primary reservoir is atomic nitrogen, while it is ammonia ice in the later times.

When the formation rate of water ice exceeds the photodesorption rate, dust grains start to be covered with water ice. Thick ice mantles are already formed at AV = 0.5 mag, and the water ice abundance reaches 10-4 at AV ~ 1.5 mag. The abundance of atomic oxygen in the gas phase is already less than 10-5 at AV ~ 1.5 mag, which corresponds to ~5% of elemental oxygen. The main surface formation pathways of water ice in our model are as follows: While the OH + H2 reaction in the gas phase has an activation energy barrier of 2100 K (Atkinson et al. 2004), Oba et al. (2012) experimentally demonstrated that Reaction (13) occurs on a cold substrate of 10 K via quantum tunneling. We used the barrier height of 2100 K to calculate the transmission probability of Reaction (13), assuming a rectangular potential barrier with a width of 1 Å. Taquet et al. (2013b) calculated the transmission probability taking into account the shape of the barrier with the Eckart model. The probability calculated by Taquet et al. (2013b) is 4 × 10-7, which is larger than that in our model by four orders of magnitude. It is likely therefore that we underestimate the rate coefficient of Reaction (13).

In our fiducial model, Reaction (12) is the most efficient to form water ice at AV ≲ 0.2 mag at which Tdust ≳ 14 K and n(H) /n(H2) ≳ 10-2, where n(i) is the number density of species i. As AV increases, the dust temperature and H/H2 abundance ratio decrease. A reduced dust temperature leads to a longer duration time of H2 in a binding site, and thus to a higher reaction efficiency (see Fig. 3); the relative contribution of Reaction (13) to water ice formation increases with increasing AV. Reaction (13) then dominates over Reaction (12) at AV ≳ 0.2 mag. Cuppen & Herbst (2007) found that Reaction (12) is the main water ice formation pathway in diffuse and translucent cloud conditions using microscopic Monte-Carlo simulations, while they found Reaction (13) to be the more important in colder and denser conditions. Our result is qualitatively consistent with their findings. We have confirmed that in our model without the reaction-diffusion competition, Reaction (12) always dominates over Reaction (13) during the simulation. As discussed in Sect. 4.3, water ice deuteration heavily depends on which formation path, Reaction (12) or (13), is more effective.

Water ice can also be formed via the sequential hydrogenation of O2 (Miyauchi et al. 2008; Ioppolo et al. 2008). This pathway is not important in our model, because the accretion rate of atomic oxygen onto the icy grain surface is lower than that of atomic hydrogen by two orders of magnitude during the simulation. Adsorbed atomic oxygen likely reacts with atomic hydrogen before another atomic oxygen accretes onto the surface.

The O2 abundance in the gas phase is as low as 10-9 during the simulation. Most of the oxygen is locked in either water ice or CO before AV reaches 1.5 mag. O2 is efficiently destroyed by photodissociation at AV< 1.5 mag, while the amount of atomic oxygen in the gas phase that is available for O2 formation is limited at AV> 1.5 mag. The model result is essentially the same as the scenario proposed by Bergin et al. (2000) to explain the low upper limits on the O2 abundance (<10-7) observationally derived in nearby clouds (Goldsmith et al. 2000; Pagani et al. 2003).

Figure 5 shows the number of icy layers in total and the fractional composition in the icy surface layers as functions of AV. The total number of ice layers formed during the simulation is ~90 monolayers. At AV ≲ 1.5 mag (Nlayer ≲ 60), the surface layers of ice mantle is predominantly composed of water, while in the later stage, other molecules, such as CO, become as abundant as water.

The ice composition at AV = 3 mag in our model is H2O:CO:CO2:CH4:NH3=100:38:1:15:15, which is compared with the observed median ice composition toward low-mass (high-mass) protostellar sources, 100:29(13):29(13):5(2):5(5) (Öberg et al. 2011). The largest discrepancy is in CO2. The CO2/H2O ratio in our model is about one order of magnitude lower than the observed ratio. We may overestimate the barrier against hopping for CO on the ice surface (Garrod & Pauly 2011). Since CO2 ice formation (e.g., CO + OH → CO2 + H) competes with water ice formation, our model may overestimate water ice formation. Our model overestimates the CH4/H2O and NH3/H2O ratios by a factor of three or more. Current gas-ice chemical models overestimate these two ratios in general (Garrod & Pauly 2011; Vasyunin & Herbst 2013a; Chang & Herbst 2014). The reason remains unclear.

thumbnail Fig. 5

The number of icy layers in total (top) and fractional compositions in the active surface ice layers (bottom) as functions of visual extinction, or time.

4.2. Ortho-to-para nuclear spin ratio of H2

OPR(H2) is determined by the competition between H2 formation on grain surfaces and ortho-para spin conversion in the gas phase through proton exchange reactions with H+ and H (Gerlich 1990; Le Bourlot 1991; Honvault et al. 2011).

If the timescales of these reactions are shorter than that of temporal variation of physical conditions, the OPR(H2) reaches the steady state value, where bo is the branching ratio to form o - H2 for H2 formation on grain surfaces (0.75, Watanabe et al. 2010). τo → p and τp → o are the spin conversion timescale between o - H2 and p - H2 through the proton exchange reactions, while τH2 is the formation timescale of H2 on grain surfaces. In the case with τH2 ≪ τo → p ≲ τp → o (or β2 ≫ 1), OPRst(H2) is ~bo/ (1−bo) = 3, i.e., the statistical value. In another extreme case, τH2 ≫ τo → p, OPRst(H2) is ~(β1 + boβ2), i.e., larger than the thermalized value β1. The derivation of Eq. (14) can be found in Appendix B (see also Le Bourlot 1991).

thumbnail Fig. 6

OPR(H2) in the fiducial model as a function of visual extinction. The dashed line shows the steady-state value given by Eq. (14) using molecular abundances in the simulations, while the dotted line shows the low-temperature thermalized value, 9exp(−170.5 /Tgas).

Figure 6 shows the OPR(H2) as a function of AV in the numerical simulation. The thermalized value (β1 = 9exp(−170.5 /Tgas)) and the steady-state value given by Eq. (14) are also shown. Since our physical model is time-dependent, the steady-state value can be calculated at given time. The expression of the thermalized value is valid at Tgas ≲ 100 K, where the population of higher rotational levels (J ≥ 2) is negligible. When most hydrogen is locked in molecular hydrogen, the OPR(H2) is already much less than the statistical value of three via ortho-para spin conversion in the gas phase, while it is much higher than the thermalized value. The OPR(H2) decreases with increasing AV, approaching the steady state value, and reaches ~5 × 10-4 at AV = 3 mag. Note that the steady state value itself depends on AV, reflecting local physical and chemical conditions. While the OPR(H2) evolves in a non-equilibrium manner, the difference between the numerical and steady-state values is at most a factor of five during the simulation. More details of the OPR(H2) evolution are discussed below.

In the H/H2 transition region at ≲0.1 mag, ortho-para spin conversion occurs predominantly through the following reactions: (17)Assuming that the gas temperature is 20 K and the H2 formation rate is half of the accretion rate of atomic hydrogen onto dust grains, we can evaluate OPRst(H2) in the H/H2 transition region by following (see Appendix C for the derivation): (18)where x(H2) and xgr are the abundances of H2 and dust grains, respectively, with respect to hydrogen nuclei, and ξH2/nH is the cosmic-ray ionization rate of H2 divided by the number density of hydrogen nuclei. Since β2< 1 is equivalent to τo → p<τH2, OPR(H2) is almost in the steady state in the H/H2 transition region as seen in Fig. 6. The equation reproduces the trend in the simulations that OPR(H2) decreases with increasing H2 abundance, and reproduces the numerical value of ~0.2 when the conversion almost finishes.

In this work, we adopt a cosmic-ray ionization rate of 1.3 × 10-17 s-1, which is appropriate for dense cloud cores (Caselli et al. 1998) and was used in B04. However, that is an order of magnitude lower than the values now thought to be appropriate for interstellar gas with low extinction (Le Petit et al. 2004; Indriolo & McCall 2012). If we adopt the higher value, which would be more appropriate for the H/H2 transition region, the OPR(H2) would be as low as 10-2 when the conversion of hydrogen into H2 almost finishes.

After the H/H2 transition, the abundance of atomic hydrogen continues to decrease toward the value n(H) ~ 1 cm-3 (ξH2/ 10-17 s-1) (e.g., Tielens 2005). The H+ abundance stays nearly constant at AV ≳ 0.5 mag, while H abundance increases with increasing AV, because of CO freeze-out and the drop of the ionization degree; so that the reaction of H with o - H2 dominates over Reaction (17) in the forward direction at AV ≳ 1.5 mag. In summary, as AV increases, the abundance of atomic hydrogen decreases leading to an increase of τH2, while the abundances of the light ions stay constant or increase leading to an decrease of τo → p. The different behaviors of a light ions and atomic hydrogen lead to a decrease of the OPR(H2) with increasing AV.

4.3. Water ice deuteration

Figure 7 shows the deuteration ratio of selected species as functions of AV. The H2D+/H and DH/H2D+ ratios increase with increasing AV, simply reflecting the decrease of the OPR(H2) and the freeze-out of CO. Although H is highly enriched in deuterium, non-deuterated H is still more abundant than its isotopologues. The atomic D/H and icy HDO/H2O ratios show more complicated behaviors. The HDO / H2O ratio in the bulk ice mantle preserves the past physical and chemical conditions which the materials experienced, while the HDO / H2O ratio in the active ice surface layers reflects local physical and chemical conditions. Interestingly, the latter ratio is less than the atomic D/H ratio by an order of magnitude except for the very early stage, where the two ratios are similar (see below for the reason). At AV = 3 mag, the HDO / H2O ratio in the bulk ice mantle is ~2 × 10-4, while that in the active surface ice layers is ~3 × 10-3. The HDO / H2O ratio in the bulk ice mantle is so low that methane and ammonia are the major deuterium reservoirs in the ice; their D/H ratios are ~10-2, while their abundances with respect to water are ≳0.1. Compared to the median composition of interstellar ices, our model overestimates the abundances of methane and ammonia with respect to water by a factor of three or more as noted in Sect. 4.1. Even considering this discrepancy, water is not the dominant deuterium reservoir in ice in our fiducial model. We note that most deuterium is in the gas phase as HD or atomic deuterium during the simulation, and only 4% of deuterium is in ice at most.

thumbnail Fig. 7

Deuteration of selected species and OPR(H2) in the fiducial model as functions of visual extinction, or time. The solid lines represent icy species in the active surface ice layers, while the dashed lines represent gaseous species. The HDO / H2O ratio in the bulk ice mantle is shown by the thick solid line.

4.3.1. The role of the ortho-to-para ratio of H2

Figure 8 shows the normalized destruction rates of H2D+ via the reactions H2D+ + CO, H2D+ + e, H2D+ + H2, and H2D+ + HD as functions of visual extinction. It shows that the main destroyer of H2D+ is H2 except at the lowest visual extinction during the simulation. Nevertheless, the effect of the OPR(H2) on ice deuteration is minor in the fiducial model; we confirmed that the difference in the HDO / H2O ratio in the bulk ice at AV = 3 mag between the models with and without nuclear spin chemistry is less than a factor of two. This is mainly because the dominant production pathway of atomic deuterium is photodissociation of HD at AV ≲ 1.5 mag, during which water ice is mainly formed. The contribution of D2 photodissociation is minor because of its low abundance. Only after the majority of elemental oxygen is locked up in water and CO does the dominant production pathway of atomic deuterium become recombination of deuterated ions with electrons. In addition, at AV ≳ 1.5 mag, the destruction rate of H2D+ by H2 is comparable to that by CO.

thumbnail Fig. 8

Normalized destruction rate of H2D+ by the reactions with CO and HD, recombination with electrons, and the backward reaction with H2 in the fiducial model as functions of AV. The normalized destruction rates by H2 assuming the thermalized ortho to para ratio of H2 are also shown for comparison.

4.3.2. The role of water ice photodissociation

In gas with low extinction, water ice is efficiently photodissociated by interstellar UV photons followed by reformation via surface reactions. The photodesorption yield of water (desorbed as OX or X2O, where X is H or D) per photodissociation is only a few % in the uppermost monolayer of the ice mantle, and the yield is lower in the deeper layers (Andersson et al. 2006). When photodesorption regulates the growth of ice mantle, icy molecules on the surface experience the dissociation-reformation cycle many times before being buried in the inert ice mantle. We find that this photodissociation-reformation cycle can cause the much lower HDO / H2O ratio in the active surface ice layers (fHDO, s) than the atomic D/H ratio (fD, s) as discussed below. In our fiducial case, fHDO, s is ≲10-3 during the simulation, which ensures that the resultant HDO / H2O ratio in the bulk ice is <10-3, although fD, s reaches >10-2 at the maximum.

Figure 9 summarizes the important reactions for H2O ice and HDO ice in our models. Photodissociation of HDO ice has two branches, (OD + H) and (OH + D) with a branching ratio of 2:1 (Koning et al. 2013). If most OH ice photofragments are converted to H2O ice, HDO ice photodissociation leads to the removal of deuterium from the water ice chemical network on the timescale of photodissociation. The efficiency of the deuterium removal depends on the probability for the OH ice photofragment to be converted back to HDO ice (pOH → HDO). pOH → HDO depends on the main formation pathway of H2O, OH + H2→ H2O + H or OH + H → H2O. The conversion of OH into HDO would mainly occur via OH + D regardless of the main formation pathway of H2O, which is confirmed in our models. There are two reasons: (i) OH + HD → HDO + H is much less efficient than OH + HD → H2O + D owing to the mass dependency of tunneling reactions (Oba et al. 2012), and (ii) the HD/H2 ratio is usually much less than the atomic D/H ratio. Then, pOH → HDO becomes smaller, i.e., the deuterium removal becomes more efficient, with an increase in the contribution of the OH+H2 pathway to the H2O formation.

thumbnail Fig. 9

Important reaction routes to form H2O and HDO on the (icy) grain surface in our models. The gray area indicates the cycle of photodissociation and reformation of water ice, which is discussed in the text.

thumbnail Fig. 10

Steady-state fHDO, s/fD, s ratio in the active ice layers as a function of the ratio of the reaction rate between OH + H2→ H2O + H and OH + H → H2O at three temperatures given by Eq. (D.10). See Appendix D for more details.

We can analytically evaluate the steady-state value of fHDO, s with respect to fD, s established through the photodissociation-reformation cycle (see Appendix D for details). Figure 10 shows the steady-state value of fHDO, s as a function of Γ ≡ ROH + H2/ROH + H, where RA + B is the rate of the surface reaction A + B. It demonstrates that when Γ ≫ 1, fHDO, s can be smaller than fD, s by more than one order of magnitude. As mentioned in Sect. 4.1, in our fiducial model, the OH+H pathway dominates over the OH+H2 pathway at AV ≲ 0.2 mag, while the OH+H2 pathway is dominant at later times. Reflecting the change of the main formation pathway, fHDO, s is similar to fD, s at AV ≲ 0.2 mag, while fHDO, s is much smaller than fD, s at later times.

The analogous mechanism is not at work for methane and ammonia ices in our model, because the reactions CH3 + H2 → CH4 + H and NH2 + H2 → NH3 + H are much less efficient than the hydrogenation of CH3 and NH2 via reactions with atomic hydrogen; the reactions with H2 have high activation energy barriers of ~6000 K (Kurylo & Timmons 1969; Aannestad 1973; Hasegawa & Herbst 1993a). Then, methane and ammonia ices are much more enriched in deuterium compared to water ice in our fiducial model.

4.4. Water vapor deuteration

The abundance of water vapor lies in the range of (1−10) × 10-9 during most of the simulation, as shown in Fig. 4. The HDO / H2O ratio in the gas phase is shown in Fig. 7. At AV ≲ 0.5 mag, HDO / H2O in the gas phase and that in the active surface ice layers are very similar; both H2O and HDO vapors are formed by photodesorption and destroyed by reaction with C+ to form HCO+ or DCO+, followed by dissociative electronic recombination to form CO. At AV ≳ 0.5 mag, HDO / H2O ratio in the gas phase is much higher than that in the active surface ice layers; the H2O vapor abundance is mostly determined by the balance between photodoesorption and photodissociation in the gas phase at AV ≲ 2.5 mag (Hollenbach et al. 2009), while HDO is mainly formed via the sequential ion-neutral reactions, starting from O + H2D+ → OD+ + H2, and destroyed by photodissociation. For example, at AV = 1.5 mag, the gas phase formation rate of HDO is faster than the photodesorption rate of HDO ice by a factor of 30, which corresponds to the difference between the HDO / H2O ratio in the gas phase and that in the active ice layers. At AV ≳ 2.5 mag, both H2O and HDO are formed via gas phase reactions and destroyed via photodissociation. Then, there is no relation between the ratios in the gas and in the active ice layers. The contribution of the photodesorption decreases as the UV field is attenuated and the fraction of water in the active ice layers decreases as shown in Fig. 5.

thumbnail Fig. 11

Fraction of elemental sulfur in the form of S+ (dashed lines), neutral atom S (dotted lines), and S-bearing molecules (the sum of CS, SO, and H2S; solid lines) in the gas phase as a function of AV in the fiducial model (blue lines) and model HM (red lines).

5. Parameter dependence

5.1. Elemental abundances of heavy metals

Figure 11 shows the fraction of elemental sulfur in gaseous species as functions of AV in the model with the high metal abundances (model HM) and in the fiducial model. The same physical model shown in Fig. 2 is used in model HM. Note that in the fiducial model 99.4% of elemental sulfur is assumed to be not available for gas and ice chemistry (Table 1). In model HM sulfur starts to be frozen out onto dust grains after the S+/S transition. When AV reaches 3 mag, 85% of elemental sulfur is locked in ice predominantly as HS and H2S, while 15% of sulfur remains in the gas phase. Then model HM is much richer in sulfur in the gas phase than the fiducial model. Figure 12 compares the abundance of electrons in model HM with that in the fiducial model. The higher abundance of S+ in model HM keeps the degree of ionization higher than the fiducial model, except for the very early stage, where the dominant positive charge carrier is C+.

thumbnail Fig. 12

HCS+/CS abundance ratio (solid lines) and electron abundance (dashed lines). Red lines indicate the model with the high metal abundances, while blue lines indicate the fiducial model with the low metal abundances. See Table 1 and Appendix E.

thumbnail Fig. 13

OPR(H2) (solid lines) and atomic D/H ratio (dotted lines) in the gas phase as functions of visual extinction. Red lines indicate the high metal case, while blue lines indicate the low metal case. The dashed lines indicate the steady-state values of the OPR(H2) given by Eq. (14).

Figure 13 shows OPR(H2) in model HM as a function of AV. The OPR(H2) in model HM is much higher at AV ≳ 0.5 mag than the fiducial model; the main destroyers of o - H2, H+ and H, are much less abundant in model HM because of the abundant gaseous sulfur. H+ is efficiently destroyed by the charge transfer to S atoms and to CS, which have lower ionization energies than atomic hydrogen. The higher degree of ionization leads to the reduced abundance of H through dissociative electron recombination. Reflecting the higher OPR(H2), the atomic D/H ratio is lower in model HM than the fiducial model especially at AV ≳ 1.5 mag.

Which model yields the better prediction for the OPR(H2) evolution in the ISM? In Appendix E, we compare observations of sulfur-bearing molecules toward diffuse/translucent/molecular clouds in the literature with our model results. We focus on the total abundance of selected sulfur-bearing molecules in the gas phase (CS, SO, and H2S) and the HCS+/CS abundance ratio. The former can probe the total abundance of elemental sulfur in the gas phase, while the latter can probe the ionization degree of the gas, which is related to the S+ abundance. We find that the fiducial model better reproduces the sulfur observations than model HM, and conclude that the fiducial model gives better predictions for the OPR(H2) evolution in the cold ISM.

5.2. The number of active surface ice layers

The parameter Nact was set to be four in the fiducial model. In the model with Nact = 1, the timescale in which surface icy molecules are buried in the inactive inert ice mantles is shorter than in the fiducial model. This limits the number of the photodissociation-reformation cycles so as to decrease the HDO/H2O ratio in the active ice layers. However, the effect on the resultant HDO/H2O ratio in the bulk ice mantle is minor (within a factor of two).

5.3. Desorption energies and reaction-diffusion competition

As discussed in Sect. 4.3 and shown in Fig. 10, water ice deuteration depends on the main formation pathway of water ice. The critical parameter is Γ = ROH + H2/ROH + H. The exact value of Γ at given physical conditions is highly uncertain at present; it strongly depends on chemical parameters and the treatment of reactions with activation energy barriers. In this subsection, we explore the impact of desorption energies of H2 and the reaction-diffusion competition on water ice deuteration. For simplicity, we assume that the desorption energies of atomic and molecular hydrogen are the same with a fixed hopping-to-desorption energy ratio of 0.5. In this case, Γ is reduced to be , where ⟨ A ⟩ is the population of species A in the active surface ice layers. Then, a smaller desorption energy leads to a smaller Γ (see Fig. 3). If the reaction-diffusion competition is turned off, Γ becomes smaller.

Figure 14 shows the impact of on the HDO / H2O ratio in the bulk ice mantle. Let us define the time when the main formation stage of water ice finishes as the time when the abundance of atomic oxygen becomes less than 1% of the elemental oxygen available for gas and ice chemistry. This happens at AV = 2−2.5 mag in our models. At that time, the HDO / H2O ratios of the bulk ice mantle are ~10-3, ~10-4, and ~10-5 in the models with K, 550 K, and 650 K, respectively. It is clear that the HDO / H2O ratio is sensitive to , while the fraction of elemental oxygen locked in water ice is not sensitive to . The models without the reaction-diffusion competition are also shown in Fig. 14. In this case, we do not see a strong dependence of the HDO / H2O ratio on , since Γ < 1 regardless of in the range of 450 K−650 K. In other words, Reaction (12) dominates over Reaction (13), so that the HDO / H2O ratio is more easily enhanced by a high atomic D/H ratio (see Fig. 10). When the water ice formation is finished, the HDO / H2O ratio in the bulk ice mantle is 6 × 10-3, which is higher than the models with the reaction-diffusion competition.

The observationally derived upper limits of the HDO / H2O ice ratio in the cold outer envelope of embedded protostars are (2−10) × 10-3 (e.g., Dartois et al. 2003), while there is still no observational constraint for the HDO / H2O ice ratio in molecular clouds. Let us assume that the HDO / H2O ice ratio in molecular clouds is less than 2 × 10-3. In Fig. 14, we can see that this low HDO / H2O ice ratio is reproduced by models with the reaction-diffusion competition, in which Reaction (13) contribute significantly to the water formation, while the ratio is reproduced only in a limited parameter space in the models without competition.

In addition, the present work suggests the possibility that water ice formed in molecular clouds is not so enriched in deuterium, compared to the water vapor observed in the vicinity of low-mass protostars where water ice is sublimated (~10-3, e.g., Persson et al. 2013). If this is the case, the enrichment of deuterium in water ice should mostly occur in the later prestellar core or/and protostellar phases, where interstellar UV radiation is heavily attenuated, CO is frozen out, and the OPR(H2) is lower than in molecular clouds. Note that even if water ice formed in molecular clouds is not enriched in deuterium compared to the elemental D/H ratio in the local ISM (1.5 × 10-5), an additional formation on small amount of water ice with a high HDO / H2O ratio in prestellar/protostellar core, e.g., x(H2O ice) = 5 × 10-6 with HDO / H2O = 2 × 10-2, can explain the observationally derived HDO / H2O ratio in the vicinity of the low-mass protostars (Furuya et al., in prep.).

thumbnail Fig. 14

HDO / H2O ratio in the bulk ice in the two grids (with and without the reaction-diffusion competition) of models as functions of AV. For the case without the competition (yellow), the solid line represents the model with Nact = 4, Edes(H2) = 550 K, and the low-metal abundances. The areas represent the models including variations in Nact and Edes(H2), while the light-colored areas represent the models including variations in Nact, Edes(H2), and the initial abundances. For the case with the competition, the models with different Edes(H2) are shown by different colors. The solid lines represent the models with Nact = 4. The areas represent models including the variations in Nact, while the light-colored areas represent the models including the variations in Nact and the initial abundances. The vertical dashed and dotted lines indicate the regions where the fraction of elemental oxygen in the atomic form becomes less than 10% and 1%, respectively, in the models represented by the solid lines.

thumbnail Fig. 15

Results from a weaker shock model plotted vs visual extinction to be compared with analogous figures for our fiducial model. a) Physical evolution after passing through the shock front (see Fig. 2). b) Fractional abundances of selected species with respect to hydrogen nuclei (see Fig. 4). c) The OPR(H2) (see Fig. 6). d) Deuteration ratios of selected species (see Fig. 7).

5.4. Ram pressure due to the accretion flow

In our fiducial physical model, the density of pre-shock materials and the velocity of accretion flow are assumed to be n0 = 10 cm-3 and v0 = 15 km s-1, respectively. The resultant post-shock material has the relatively high density of ~104 cm-3, compared to the average density of giant molecular clouds derived from observations (Bergin & Tafala 2007, and references therein). In this subsection, we briefly present a model with a smaller ram pressure, n0 = 3 cm-3 and v0 = 10 km s-1 (model 3 in H10), in which the density of post-shock materials is ~103 cm-3. The goal is to check if our primary results presented in Sect. 4 do not change qualitatively. In this model, the internal pressure due to self gravity overwhelms the ram pressure at AV> 1.6 mag (cf. Hartmann et al. 2001). Nevertheless, we use the shock model until AV reaches 3 mag.

Figure 15 summarizes the results from the model with the smaller ram pressure. We confirmed that most of our qualitative chemical results presented earlier in this paper do not change. The main differences are (i) molecular formation requires a larger shielding column density because of the lower gas density (B04; H10); (ii) the OPR(H2) is lower when most hydrogen is locked in H2 reflecting higher ξH2/nH (Eq. (18)); and (iii) the OPR(H2) is almost at steady-state during the simulation because of the longer timescale of the physical evolution.

6. Comparisons with previous theoretical studies

Water ice deuteration has been numerically studied over wide ranges of the physical conditions: nH = 104−106 cm-3, T = 10−20 K, and AV = a few-10 mag. Most of the previous studies used a pseudo-time-dependent chemical model, in which physical conditions are constant with time (e.g., Taquet et al. 2013b; Sipilä et al. 2013; Lee & Bergin 2015), while some investigated deuterium chemistry in gravitationally collapsing clouds/cores (e.g., Aikawa et al. 2012; Taquet et al. 2014). The general conclusions of the previous studies are that (i) the HDO / H2O ice ratio is controlled by the temporal evolution of the atomic D/H ratio (one exception is Cazaux et al. 2011, but see below) and (ii) the evolution of the atomic D/H ratio strongly depends on the initial OPR(H2), which is treated as a free parameter. The rationale for point (i) is that H2O and HDO are predominantly formed via sequential reactions of atomic hydrogen/deuterium with atomic oxygen in the previous models.

The present study has considered the chemistry during which cold and dense conditions are first achieved by the passage of a shock through precursor H i-dominated clouds. This approach allows us to study the evolution of OPR(H2) and water ice deuteration without arbitrary assumptions concerning the initial OPR(H2). Contrary to point (i), the present study has shown that the HDO / H2O ice ratio is not always controlled by the atomic D/H ratio. We have shown that deuterium can be removed from water ice chemistry on the timescale of water ice photodissociation, when Reaction (13) dominates over Reaction (12) or equivalently Γ > 1. As a consequence, the HDO / H2O ratio in the active surface ice layers can be much smaller than the atomic D/H ratio. On the other hand, when Γ < 1, the HDO / H2O ice ratio scales with the atomic D/H ratio as predicted by the previous studies. The former condition (Γ > 1) favors gas with high but not very high extinction, where the dust temperature is low and the gaseous H/H2 abundance ratio is low, while the latter condition favors gas with a low extinction. The onset of water ice mantle formation requires a threshold extinction, which depends on the UV field, the gas density, and the dust temperature (e.g., Tielens 2005). Therefore, the critical question here is what fraction of water ice is formed with Γ > 1 in molecular clouds. In our fiducial case, for example, the majority of water ice is formed with Γ > 1. However, the quantification depends on uncertain chemical parameters, such as the hopping timescale of H and H2 on a surface and the tunneling transmission probability of Reaction (13), as well as adopted physical models. We also note that when the reaction-diffusion competition is turned off, Γ is always less than unity in our models. To the best of our knowledge, no previous study focusing on ice deuteration considered the reaction-diffusion competition (except for Cazaux et al. 2011). This explains why the previous studies commonly claimed point (i).

Cazaux et al. (2011) studied the deuteration of water ice in a translucent cloud followed by gravitational collapse. They claimed that the HDO / H2O bulk ice ratio scales with the D/H2 ratio at Tdust ≲ 15 K, while it scales with the atomic D/H ratio at higher dust temperatures. They demonstrated that for Tdust ≲ 15 K, OH and OD are mainly formed by O + H2 → OH + H and O + D → OD, respectively, followed by hydrogenation to form H2O and HDO. At higher dust temperatures, OH and OD are mainly formed by the reaction of atomic hydrogen (deuterium) with atomic oxygen in their model. At Tdust ≲ 15 K, the HDO / H2O bulk ice ratio is ~10-6 at the end of their simulations, while it is ~10-3 at Tdust ≳ 15 K. The key assumption they made was that the surface reaction, O + H2 → OH + H, occurs efficiently through quantum tunneling. Since the reaction O + H2 → OH + H is endothermic by 960 K, the validity of the assumption has been questioned by several authors (Oba et al. 2012; Taquet et al. 2013a; Lamberts et al. 2014). Recently, Lamberts et al. (2014) concluded that the contribution of the O + H2 pathway to the total OH formation in the ISM is small (likely much less than 10%), based on their experiments and modeling.

7. Summary

We have investigated water deuteration and the OPR(H2) during the formation and early evolution of a molecular cloud, following the scenario that H i-dominated gas is swept and accumulated by global accretion flows to form molecular clouds. We used the one-dimensional shock model developed by Bergin et al. (2004) and Hassel et al. (2010) to follow the physical evolution of post-shock materials, combined with post-processing detailed gas-ice chemical simulations. Our findings are summarized as follows.

  • 1.

    In the H/H2 transition region, the OPR(H2) is already much less than the statistical value of three, because the timescale of spin conversion through the reaction of o - H2 with H+ is shorter than the timescale of H2 formation on the grain surface. When most hydrogen is locked in H2, the OPR(H2) is nearly at steady state with the value of ~0.1 in our fiducial model. The exact value depends on the ionization conditions and the gas temperature (see Appendix C and Fig. C.1). At later times, the OPR(H2) decreases in a non-equilibrium manner as the gas accumulates, reflecting local physical and chemical conditions. The evolution of the OPR(H2) depends on the abundances of elements heavier than oxygen, especially sulfur. A higher sulfur abundance leads to the decrease of the abundances of H+ and H, which leads to a longer ortho-para spin conversion timescale in the gas phase.

  • 2.

    The HDO / H2O ratio in the bulk ice is as low as 10-4 at the end of our fiducial simulation where most of oxygen is already locked up in molecules. The key mechanism to suppress water ice deuteration is the cycle of photodissociation and reformation of water ice on the icy surface, which removes deuterium from water ice chemistry at the timescale of photodissociation. The efficiency of the mechanism depends on which formation path, the barrierless reaction OH + H → H2O or the barrier-mediated reaction OH + H2→ H2O + H, is more effective, because the pathway including barrier-mediated reactions favors hydrogenation over deuteration. Depending on the contribution of the OH + H2 pathway to water ice formation, the resultant HDO / H2O ratio in the bulk ice ranges from 10-5 to 10-3. The contribution of the OH + H2 pathway to water ice formation depends strongly on whether or not we include the reaction-diffusion competition and hopping timescale of H2 on the surface. The OPR(H2) plays a minor role in water ice deuteration, because the production of atomic deuterium is dominated by photodissociation of HD at the main formation stage of water ice, rather than the electron dissociative recombination of deuterated ions.

  • 3.

    The above results suggest the possibility that water ice formed in molecular clouds is deuterium-poor, compared to the water vapor observed in the vicinity of protostars where water ice is sublimated. If this is the case, the enrichment of deuterium in water should mostly occur at somewhat later evolutionary stages of star formation, i.e., cold prestellar/protostellar cores, where interstellar UV radiation is heavily attenuated, CO is frozen out, and OPR(H2) is lower than in molecular clouds. When the protostellar core begins to warm up, the situation becomes more complex.

  • 4.

    The HDO / H2O ratio in the gas phase, the active surface ice layers, and in the bulk ice mantle can have different values. Even if the H2O vapor abundance is determined by the balance between photodesorption and photodissociation, it does not necessarily mean that the HDO / H2O ratio in the gas phase and in the surface ice layers are similar.

Acknowledgments

K.F. thanks Ewine F. van Dishoeck, Catherine Walsh, and Maria N. Drozdovskaya for stimulating discussions, Carina Arasa for providing useful data on photolysis, and Vianney Taquet for constructive comments on the manuscript. We thank the anonymous referee for helpful suggestions on the manuscript. K.F. is supported by the Research Fellowship from the Japan Society for the Promotion of Science (JSPS). Y.A. acknowledges the support by JSPS KAKENHI Grant Numbers 23103004 and 23540266. A.V. acknowledges the financial support of the European Research Council (ERC; project PALs 320620). E.H. acknowledges the support of the National Science Foundation for his program in astrochemistry, and support from the NASA Exobiology and Evolutionary Biology program under subcontract from Rensselaer Polytechnic Institute. Some kinetic data we used have been downloaded from the online database KIDA (Wakelam et al. 2012, http://kida.obs.u-bordeaux1.fr).

References

  1. Aannestad, P. A. 1973, ApJS, 25, 223 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aikawa, Y., & Herbst, E. 1999, ApJ, 526, 314 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  3. Aikawa, Y., Wakelam, V., Hersant, F., Garrod, R., & Herbst, E. 2012, ApJ, 760, 40 [NASA ADS] [CrossRef] [Google Scholar]
  4. Albertsson, T., Semenov, D., & Henning, T. 2014, ApJ, 784, 39 [NASA ADS] [CrossRef] [Google Scholar]
  5. Andersson, S., & van Dishoeck, E. F. 2008, A&A, 491, 907 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Andersson, S., Al Halabi, A., Kroes, G.-J., & van Dishoeck, E. F. 2006, J. Chem. Phys., 124, 4715 [Google Scholar]
  7. Anderson, D. E., Bergin, E. A., Maret, S., & Wakelam, V. 2013, ApJ, 779, 141 [NASA ADS] [CrossRef] [Google Scholar]
  8. Arasa, C,Koning, J., Kroes, G.-J., Walsh, C., & van Dishoeck, E. F. 2015, A&A, 575, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Atkinson, R., Baulch, D. L., Cox, R. A., et al. 2004, Atmos. Chem. Phys., 4, 1461 [NASA ADS] [CrossRef] [Google Scholar]
  10. Awad, Z., Chigai, T., Kimura, Y., Shalabiea, O. M., & Yamamoto, T. 2005, ApJ, 626, 262 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bergin, E. A., & Tafalla, M. 2007, ARA&A, 45, 339 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bergin, E. A., Melnick, G. J., Stauffer, J. R., et al. 2000, ApJ, 539, L129 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bergin, E. A., Hartmann, L.W., Raymond, J. C., & Ballesteros-Paredes, J. 2004, ApJ, 612, 921 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bertin, M., Fayolle, E. C., Romanzin, C., et al. 2012, Phys. Chem. Chem. Phys., 14, 9929 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bockelée-Morvan, D., Biver, N., Swinyard, B., et al. 2012, A&A, 544, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bouwman, J, Mattioda, A. L., Linnartz, H., & Allamandola, L. J. 2011, A&A, 525, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Brünken, S., Sipilä, O., Chambers, E. T., et al. 2014, Nature, 516, 219 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  18. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [Google Scholar]
  19. Caselli, P., Stantcheva, T., Shalabia, O., Shematovich, V. I., & Herbst, E. 2002, Planet. Space Sci., 50, 1257 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222 [NASA ADS] [CrossRef] [Google Scholar]
  21. Cazaux, S., Caselli, P., & Spaans, M. 2011, ApJ, 741, L34 [NASA ADS] [CrossRef] [Google Scholar]
  22. Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, Protostars & Planets VI, 859 [Google Scholar]
  23. Chang, Q., & Herbst, E. 2014, ApJ, 787, 135 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chang, Q., Cuppen, H. M., & Herbst, E. 2007, A&A, 469, 973 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Cleeves, L. I., Bergin, E. A., Alexander, C. M. O., et al. 2014, Science, 345, 1590 [NASA ADS] [CrossRef] [Google Scholar]
  26. Coutens, A., Vastel, C., & Caux, E., et al. 2012, A&A, 539, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Coutens, A., Vastel, C., Cabrit, S., et al. 2013, A&A, 560, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Coutens, A., Jørgensen, J. K., Persson, M., et al. 2014a, ApJ, 792, L5 [NASA ADS] [CrossRef] [Google Scholar]
  29. Coutens, A., Vastel, C., Hincelin, U., et al. 2014b, MNRAS, 445, 1299 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cuppen, H. M., & Herbst, E. 2007, ApJ, 668, 294 [Google Scholar]
  31. Dartois, E., Thi, W. F., Geballe, T. R., et al. 2003, A&A, 399, 1020 [Google Scholar]
  32. Dalgarno, A., Black, J. H., & Weisheit, J. C. 1973, ApJ, 14, L77 [Google Scholar]
  33. Daranlot, J., Hincelin, U., Bergeat, A., et al. 2012, PNAS, 109, 10233 [Google Scholar]
  34. Dickens, J. E., Irvine, W. M., Snell, R. L., et al. 2000, ApJ, 542, 870 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dobbs, C. L., Krumholz, M. R., Ballesteros-Paredes, J., et al. 2014, Protostars & Planets VI, 3 [Google Scholar]
  36. Draine, B. T. 1978, ApJS, 36, 595 [NASA ADS] [CrossRef] [Google Scholar]
  37. Draine, B. T., & Bertoldi, F. 1996, ApJ, 468, 269 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] [Google Scholar]
  39. Dulieu, F., Congiu, E., Noble, J., et al. 2013, Nat. Sci. Rep., 3, 1338 [NASA ADS] [CrossRef] [Google Scholar]
  40. Fayolle, E. C., Bertin, M., Romanzin, C., et al. 2011, ApJ, 739, L36 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fayolle, E. C., Bertin, M., Romanzin, C., et al. 2013, A&A, 556, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Flower, D. R., Pineau Des Forêts, G., Walmsley, C. M. 2006, A&A, 449, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Furuya, K., Aikawa, Y., Tomida, K., et al. 2012, ApJ, 758, 86 [NASA ADS] [CrossRef] [Google Scholar]
  45. Furuya, K., Aikawa, Y., Nomura, H., Hersant, F., & Wakelam, V. 2013, ApJ, 779, 11 [NASA ADS] [CrossRef] [Google Scholar]
  46. Garrod, R. T. 2008, A&A, 491, 239 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Garrod, R. T., & Herbst, E. 2006, A&A, 457, 927 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Garrod, R. T., & Pauly, T. 2011, ApJ, 735, 15 [NASA ADS] [CrossRef] [Google Scholar]
  49. Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Gerlich, D. J. 1990, Chem. Phys., 92, 2377 [Google Scholar]
  51. Gerlich, D., Herbst, E., & Roueff, E. 2002, Plan. Sp. Sci., 50, 1275 [NASA ADS] [CrossRef] [Google Scholar]
  52. Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJS, 151, 35 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  53. Goldsmith, P. F., Melnick, G. J., Bergin, E. A., et al. 2000, ApJ, 539, L123 [NASA ADS] [CrossRef] [Google Scholar]
  54. Graedel, T. E., Langer, W. D., & Frerking, M. A. 1982, ApJS, 48, 321 [NASA ADS] [CrossRef] [Google Scholar]
  55. Gredel, R., Lepp, S., Dalgrano, A., & Herbst, E. 1989, ApJ, 347, 289 [NASA ADS] [CrossRef] [Google Scholar]
  56. Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, ApJ, 757, 185 [NASA ADS] [CrossRef] [Google Scholar]
  57. Harada, N., Herbst, R., & Wakelam, V. 2010, ApJ, 721, 1570 [NASA ADS] [CrossRef] [Google Scholar]
  58. Harada, N., Herbst, R., & Wakelam, V. 2012, ApJ, 756, 104 [NASA ADS] [CrossRef] [Google Scholar]
  59. Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852 [NASA ADS] [CrossRef] [Google Scholar]
  60. Hasegawa, T. I., & Herbst, E. 1993a, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  61. Hasegawa, T. I., & Herbst, E. 1993b, MNRAS, 263, 589 [NASA ADS] [CrossRef] [Google Scholar]
  62. Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167 [NASA ADS] [CrossRef] [Google Scholar]
  63. Hassel, G. E., Herbst, E., & Bergin, E. A. 2010, A&A, 515, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Heiles, C., & Troland, T. H. 2003, ApJ, 586, 1067 [NASA ADS] [CrossRef] [Google Scholar]
  65. Heiles, C., & Troland, T. H. 2005, ApJ, 624, 773 [NASA ADS] [CrossRef] [Google Scholar]
  66. Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309 [NASA ADS] [Google Scholar]
  67. Hincelin, U., Herbst, E., Chang, Q., et al. 2014, 69th Int. Symp. Mol. Spectrosc., 9 [Google Scholar]
  68. Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [NASA ADS] [CrossRef] [Google Scholar]
  69. Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155 [NASA ADS] [CrossRef] [Google Scholar]
  70. Hollenbach, D., Kaufman, M. J., Bergin, E. A., & Melnick, G. J. 2009, ApJ, 690, 1497 [CrossRef] [Google Scholar]
  71. Honvault, P., Jorfi, M., González-Lezana, T., Faure, A., & Pagani, L. 2011, Phys. Rev. Lett., 107, 023201 [NASA ADS] [CrossRef] [Google Scholar]
  72. Hornekaer, L., Baurichter, A., Petrunin, V. V., et al. 2005, J. Chem. Phys., 122 [Google Scholar]
  73. Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  74. Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
  75. Inoue, T., & Inutsuka, S. 2009, ApJ, 704, 161 [NASA ADS] [CrossRef] [Google Scholar]
  76. Inoue, T., & Inutsuka, S. 2012, ApJ, 759, 35 [NASA ADS] [CrossRef] [Google Scholar]
  77. Inutsuka, S.-I., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Ioppolo, S., Cuppen, H. M., Romanzin, C., van Dishoeck, E. F., & Linnartz, H. 2008, ApJ, 686, 1474 [NASA ADS] [CrossRef] [Google Scholar]
  79. Ioppolo, S., Cuppen, H. M., Romanzin, C., van Dishoeck, E. F., & Linnartz, H. 2010, Phys. Chem. Chem. Phys., 12, 12065 [NASA ADS] [CrossRef] [Google Scholar]
  80. Iqbal, W., Acharyya, K., & Herbst, E. 2014, ApJ, 784, 139 [NASA ADS] [CrossRef] [Google Scholar]
  81. Jørgensen, J. K., & van Dishoeck, E. F. 2010, ApJ, 725, L172 [NASA ADS] [CrossRef] [Google Scholar]
  82. Joseph, C. L., Snow, Jr., T. P.,Seab, C. G., & Crutcher, R. M. 1986, ApJ, 309, 771 [NASA ADS] [CrossRef] [Google Scholar]
  83. Kama, M., Caux, E., Lopez-Sepulcre, A., et al. 2015, A&A, 574, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Kamp, I., & Bertoldi, F. 2000, A&A, 353, 276 [NASA ADS] [Google Scholar]
  85. Kawamura, A., Mizuno, Y., Minamidani, T., et al. 2009, ApJS, 184, 1 [NASA ADS] [CrossRef] [Google Scholar]
  86. Keller, L. P., Hony, S., Bradley, J. P., et al. 2002, Nature, 417, 148 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  87. Koning, J., Kroes, G. J., & Arasa, C. 2013, J. Chem. Phys., 138, 104701 [NASA ADS] [CrossRef] [Google Scholar]
  88. Kurylo, M. J., & Timmons, R. B. 1969, J. Chem. Phys., 50, 5076 [NASA ADS] [CrossRef] [Google Scholar]
  89. Lamberts, T., Cuppen, H. M., Fedoseev, G., et al. 2014, ApJ, 570, 57 [Google Scholar]
  90. Le Bourlot, J. 1991, A&A, 242, 235 [NASA ADS] [Google Scholar]
  91. Lee, J.-E., & Bergin, E. A. 2015, ApJ, 799, 104 [NASA ADS] [CrossRef] [Google Scholar]
  92. Le Gal, R., Hily-Blant, P., Faure, A., et al. 2014, A&A, 562, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Le Petit, F., Roueff, E., & Herbst, E. 2004, A&A, 417, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Li, X., Heays, A. N., Visser, R., et al. 2013, A&A, 555, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Linsky, J. L. 2003, Space Sci. Rev., 106, 49 [NASA ADS] [CrossRef] [Google Scholar]
  96. Liu, F.-C., Parise, B., Kristensen, L., et al. 2011, A&A, 527, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Lucas, R., & Liszt, H. S. 2002, A&A, 384, 1054 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Maret, S., Bergin, E. A., & Lada, C. J. 2006, Nature, 442, 425 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  99. Mason, N. J., Dawes, A., Holtom, P. D., et al. 2006, Faraday discuss., 133, 311 [NASA ADS] [CrossRef] [Google Scholar]
  100. Millar, T. J., Bennet, A., & Herbst, E. 1989, ApJ, 340, 906 [NASA ADS] [CrossRef] [Google Scholar]
  101. Minissale, M., & Dulieu, F. 2014, J. Chem. Phys., 140, 074705 [NASA ADS] [CrossRef] [Google Scholar]
  102. Miyauchi, N., Hidaka, H., Chigai, T., et al. 2008, Chem. Phys. Lett., 456, 27 [NASA ADS] [CrossRef] [Google Scholar]
  103. Montaigne, H., Geppert, W. D., Semaniak, J., et al. 2005, ApJ, 631, 653 [NASA ADS] [CrossRef] [Google Scholar]
  104. Oba, Y., Watanabe, N., Hama, T., et al. 2012, ApJ, 749, 67 [NASA ADS] [CrossRef] [Google Scholar]
  105. Öberg, K. I., Fuchs, G. W., Awad, Z., et al. 2007, ApJ, 662, L23 [NASA ADS] [CrossRef] [Google Scholar]
  106. Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 496, 281 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011, ApJ, 740, 109 [NASA ADS] [CrossRef] [Google Scholar]
  108. Ohishi, M., Irvine, W. M., & Kaifu, N. 1992, in Astrochemistry of Cosmic Phenomena, IAU Symp., 150, 171 [NASA ADS] [CrossRef] [Google Scholar]
  109. Pagani, L., Wannier, P. G., Frerking, M. A., et al. 1992, A&A, 258, 472 [NASA ADS] [Google Scholar]
  110. Pagani, L., Olofsson, A. O. H., Bergman, P., et al. 2003, A&A, 402, L77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Pagani, L., Vastel, C., Hugo, E., et al. 2009, A&A, 494, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Pagani, L., Roueff, E., & Lesaffre, P. 2011, A&A, 739, L35 [Google Scholar]
  113. Parise, B., Simon, T., Caux, E., et al. 2003, A&A, 410, 897 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Persson, M. V., Jørgensen, J. K., & van Dishoeck, E. F. 2013, A&A, 549, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Persson, M. V., Jørgensen, J. K., van Dishoeck, E., & Harsono, D. 2014, ApJ, 563, 74 [Google Scholar]
  116. Roberts, H., Herbst, E., & Millar, T. J. 2002, MNRAS, 336, 283 [NASA ADS] [CrossRef] [Google Scholar]
  117. Roberts, H., Herbst, E., & Millar, T. J. 2004, A&A, 424, 905 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Sofia, U. J. 2004, in Astrophysics of Dust, eds. A. N. Witt, G. C. Clayton, & B. T. Draine (San Francisco, CA: ASP), ASP Conf. Ser., 309, 393 [Google Scholar]
  119. Sipilä, O., Caselli, P., & Harju, J. 2013, A&A, 554, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Taquet, V., López-Sepulcre, A., Ceccarelli, C., et al. 2013a, ApJ, 768, L29 [NASA ADS] [CrossRef] [Google Scholar]
  121. Taquet, V., Peters, P., Kahane, C., et al. 2013b, A&A, 550, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Taquet, V., Charnley, S. B., Sipilä, O. 2014, ApJ, 791, 1 [NASA ADS] [CrossRef] [Google Scholar]
  123. Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge: Cambridge Univ. Press), 219 [Google Scholar]
  124. Tielens, A. G. G. M., & Hagen, W. 1982, A&A, 114, 245 [NASA ADS] [Google Scholar]
  125. Turner, B. E. 1996, ApJ, 461, 246 [NASA ADS] [CrossRef] [Google Scholar]
  126. van Dishoeck, E. F., Bergin, E. A., Lis, D. C., & Lunine, J. I. 2014, Protostars & Planets VI, 835 [Google Scholar]
  127. Vasyunin, A. I., & Herbst, E. 2013a, ApJ, 762, 86 [NASA ADS] [CrossRef] [Google Scholar]
  128. Vasyunin, A. I., & Herbst, E. 2013b, ApJ, 769, 34 [Google Scholar]
  129. Vidali, G., Ihm, G., Kim, H., & Cole, M. 1991, Surface Sci. Rep., 12, 135 [Google Scholar]
  130. Villanueva, G. L., Mumma, M. J., Bonev, B. P., et al. 2009, ApJ, 690, L5 [NASA ADS] [CrossRef] [Google Scholar]
  131. Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  132. Visser, R., Doty, S. D., & van Dishoeck, E. F. 2011, A&A, 534, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  133. Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371 [NASA ADS] [CrossRef] [Google Scholar]
  134. Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, ApJS, 199, 21 [NASA ADS] [CrossRef] [Google Scholar]
  135. Walmsley, C. M., Flower, D. R., & Pineau Des Forêts, G. 2004, A&A, 418, 1035 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Watanabe, N., Kimura, Y., Kouchi, A., et al. 2010, ApJ, 714, L233 [NASA ADS] [CrossRef] [Google Scholar]
  137. Westley, M. S., Baragiola, R. A., Johnson, R. E., & Baratta, G. A. 1995, Nature, 373, 405 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  138. Willacy, K., & Woods, P. M. 2009, ApJ, 703, 479 [NASA ADS] [CrossRef] [Google Scholar]
  139. Whittet, D. C. B. 1993, in Dust and Chemistry in Astronomy (Bristol and Philadelphia: Institute of Physics Publishing), 9 [Google Scholar]
  140. Whittet, D. C. B. 2003, Dust in the galactic environment (Bristol: Institute of Physics Publishing) [Google Scholar]
  141. Whittet, D. C. B., Shenoy, S. S., Bergin, E. A., et al. 2007, ApJ, 655, 332 [NASA ADS] [CrossRef] [Google Scholar]
  142. Whittet, D. C. B., Poteet, C. A., Chiar, J. E., et al. 2013, ApJ, 774, 102 [NASA ADS] [CrossRef] [Google Scholar]
  143. Wolcott-Green, J., & Haiman, Z. 2011, MNRAS, 412, 2603 [NASA ADS] [CrossRef] [Google Scholar]
  144. Yang, L., Ciesla, F. J., & Alexander, C. M. O. 2013, Icarus, 226, 256 [NASA ADS] [CrossRef] [Google Scholar]

Online material

Appendix A: The condition for the nuclear spin ratio of H2 to slow down deuteration processes

Here we derive the condition determining how much o - H2 is needed to slow down deuterium fractionation driven by Reaction (1). It is well established that the destruction of H2D+ by electrons and CO suppresses the fractionation as well as by H2. H2D+ also reacts with HD to form the doubly deuterated species D2H+, leading to further fractionation. Considering the balance of formation and destruction, the steady-state H2D+/H ratio can be calculated as (A.1)where n(i) is the number density of species i, kCO is the rate coefficient of H2D+ with CO, ke is the rate coefficient of dissociative electron recombination of H2D+, kHD is the rate coefficient of H2D+ with HD, and k1f and k1b are the effective rate coefficients for Reaction (1) in the forward direction and in the backward direction, respectively. Since the forward Reaction (1) is exothermic, k1f does not depend on the OPR(H2). We can say that the presence of o - H2 slows down the fractionation when k1bn(H2) is greater than kCOn(CO), ken(e), and kHDn(HD).

The endothermicity of Reaction (1) in the backward direction depends on the nuclear spin states of H2 and H2D+. The following set of reactions can occur: The rate coefficients of the reactions can be found in Hugo et al. (2009), who calculated them with the assumption that the reaction proceeds by a scrambling mechanism in which all protons are equivalent; results are different if long-range hopping is the effective mechanism. The effective rate coefficient k1b can be defined as follows (Gerlich et al. 2002; Lee & Bergin 2015): (A.6)From Eq. (A.6), we can express k1b as a function of ortho-to-para ratios of H2 and H2D+ (OPR(H2D+)), and gas temperature.

The OPR(H2D+) in the dense ISM is mainly determined by the following reactions (Gerlich et al. 2002; Hugo et al. 2009): Then, assuming steady state, the OPR(H2D+) is given as a function of OPR(H2) and gas temperature: (A.11)From Eqs. (A.6) and (A.11), we can express k1b as a function of OPR(H2) and the gas temperature. In Fig. 1, we show threshold abundances of CO and electrons as functions of the OPR(H2), below which Reaction (1) in the backward reaction is more efficient than the destruction by CO and electrons. Roughly speaking, the backward reaction rate exceeds the destruction rate by CO when OPR(H2) ≳ 20n(CO) /n(H2) at ≤20 K, while it exceeds the recombination rate with electrons when OPR(H2) ≳ 3 × 103n(e) /n(H2) at ≤20 K. Assuming the canonical HD abundance with respect to H2, 3 × 10-5, the condition k1bn(H2) >kHDn(HD) corresponds to OPR(H2) ≳ 6 × 10-4 at ≤20 K.

Appendix B: Analytical treatment of H2 ortho-to-para ratio

In this Appendix, we derive an analytical solution of OPR(H2) when abundances of ionic species and the H2 formation rate on grain surfaces are given. Let us consider the following set of differential equations, which describe the temporal variations of the abundances of o - H2 and p - H2 in the gas and solid phases:

where F, W, D are the adsorption rate, desorption rate, and rate of H2 destruction via e.g., photodissociation, respectively. ⟨ A ⟩ in Eqs. (B.3) and (B.4) is the population of species A in the active surface ice layers, and thus ⟨ A ⟩ ngr is the number density of species A in the surface ice layers per unit volume of gas. RH2 is the formation rate of H2 on grain surfaces, while bo is the branching ratio to form o - H2. In Eqs. (B.1) and (B.2), we neglect H2 formation in the gas phase for simplicity. In Eqs. (B.3) and (B.4), we do not consider nuclear spin conversion on the surface, though it is straightforward to include the spin conversion on the surface in the following analysis.

Combining Eqs. (B.1)−(B.4), we get (B.5)where we used n(o - H2) + n(p - H2) = n(H2) and n(H2) ≫ ⟨ H2 ⟩ ngr. The latter should be valid in molecular gases; the binding energy of H2 on a H2 substrate is only 23 K, corresponding to an sublimation timescale of ~10-10 s at T = 10 K (Vidali et al. 1991; Cuppen & Herbst 2007), while the adsorption timescale is ~109/nH yr. The H2 formation timescale, τH2, was defined as n(H2) /RH2. We also assumed that the rate coefficients of reactions to destroy o - H2 and p - H2 are the same, i.e., D(o - H2) /D(p - H2) = n(o - H2) /n(p - H2). Equation (B.5) describes the time evolution of the OPR(H2). The terms τo → p, τp → o, and τH2 are time-dependent in general. It is straightforward to solve Eq. (B.5) in the astrochemical simulations without nuclear spin state chemistry, and one may obtain a reasonable approximation of the temporal variation of the OPR(H2). In order to solve Eq. (B.5) analytically, we consider τo → p, τp → o, and τH2 as constant. This assumption is valid when the ortho-to-para spin conversion time scale is longer than the chemical (formation and destruction) timescale of hydrogen and light ions. The solution is given as follows: where τopr gives the characteristic timescale of OPR(H2) evolution, and OPR0(H2) is the initial OPR(H2). The steady-state value of OPR(H2) (OPRst(H2) ≡ OPR(H2)(t → ∞)) is given in Eq. (14), which was derived by Le Bourlot (1991) in a different manner.

When OPRst(H2) ≪ OPR0(H2) ≪ 1, Eq. (B.6) can be simply rewritten as OPR(H2)(t) ≈ OPRst(H2) + OPR0(H2)exp(−t/τopr). Then, it takes a greater time than ln [ OPR0(H2) / OPRst(H2) ] τopr to reach the steady state value.

Appendix C: H2 ortho-to-para ratio when the conversion of hydrogen into H2 is almost complete

Evolution of molecular abundances have often been investigated via pseudo-time dependent models in which hydrogen is assumed to be in H2 at the beginning of the calculation. In such models, the molecular D/H ratio depends on the initial OPR(H2), which is treated as a free parameter (e.g., Flower et al. 2006). In this appendix we analytically derive the OPR(H2) when the conversion of hydrogen into H2 is almost complete.

During the H/H2 transition, ortho-para spin conversion occurs through Reaction (17). H+ is primarily formed via cosmic-ray/X-ray ionization of atomic hydrogen, while it is mainly destroyed via recombination with electrons, and charge transfer to other species, such as atomic deuterium and oxygen (e.g., Dalgarno et al. 1973). The former is valid when n(H) /n(H2) >ξH/ (bH+ξH2) ≈ 0.05, where ξH and ξH2 are the ionization rates of atomic hydrogen and H2, respectively. bH+ is the branching ratio to form H+ for the ionization of H2. At steady state, the number density of H+ can be given as (C.1)where (C.2)where k(X + Y) is the rate coefficient of reaction X + Y. We assumed that the number density of electrons is equal to that of carbon ions. The carbon ion is the dominant form of carbon in the H/H2 transition region, while oxygen and deuterium are predominantly in atomic form. With Eq. (C.1), we can evaluate β1 and β2 (Eqs. (15) and (16)) by the following: where x(i) is the abundance of species i with respect to hydrogen nuclei, nH is the number density of hydrogen nuclei, vth is the thermal velocity of atomic hydrogen, and k17f and k17b are the rate coefficients of Reaction (17) in the forward direction and the backward direction, respectively. We used ξH2 = 2.2ξH. We assumed that the H2 formation rate is half of the accretion rate of atomic hydrogen onto dust grains. By substituting Eqs. (C.3) and (C.4) into Eq. (14), we can get the OPR(H2) in the H/H2 transition regime under the steady state assumption as a function of ξH2/nH and gas temperature. In Fig. C.1, we show the OPR(H2) when the the conversion of hydrogen into H2 is (almost) complete, i.e., x(H2) = 0.5. Above the thick dashed line, where τo → p is smaller than τH2, the steady-state assumption is justified.

It is clear that the OPR(H2) is greater than unity only when H2 formation occurs at low ionization (ξH2/nH< 10-22 cm3 s-1) and/or warm conditions. In our fiducial model, this occurs at Tgas> 50 K, but the exact value depends on ξH2/nH. In those regions with such warm gas temperatures, the dust temperature would also be warm. At Tdust ≳ 20 K, the H2 formation rate may drop considerably, depending on characteristics of chemisorption sites on grain surfaces (e.g., Hollenbach & Salpeter 1971; Cazaux & Tielens 2004; Iqbal et al. 2014), which is not considered here.

thumbnail Fig. C.1

OPR(H2) when the conversion of hydrogen into H2 is almost complete under the steady-state assumption. Above the thick dashed line, where τo → p is smaller than τH2, the steady-state assumption is justified. Above the thin dashed line, where β1 is larger than β2, the OPR(H2) is similar to the low-temperature thermalized value of 9exp(−170.5 /Tgas). See Appendix C.

Appendix D: Scaling relation of the HDO/H2O ratio in the surface ice layers with the atomic D/H ratio

Here we derive the scaling relation of the HDO / H2O ratio in the chemically active surface ice layers with the atomic D/H ratio, which is applicable under the UV irradiation conditions where photodesorption regulates the growth of ice mantles. Figure 9 shows the important reactions for H2O ice and HDO ice in our models, and they are considered in the following analysis. Let us consider the following set of differential equations which describe temporal variations of the abundances of H2O, HDO, OH and OD in the active surface ice layers: where kph and kphdes are the rate coefficients of photodissociation and photodesorption in the active surface ice layers, respectively, of H2O ice and HDO ice. bOH is the branching ratio to form OH + D for HDO ice photodissociation (~0.3, Koning et al. 2013). Combining Eqs. (D.1)−(D.4), the evolutionary equation of the HDO / H2O ratio in the active ice layers (fHDO, s = ⟨ HDO ⟩ / ⟨ H2O ⟩) can be written as follows: (D.5)where we used the inequalities ⟨ OH ⟩ ≪ ⟨ H2O ⟩, ⟨ OD ⟩ ≪ ⟨ HDO ⟩, and fHDO, s ≪ 1, which are verified by our numerical simulations. We defined fD, s = ⟨ D ⟩ / ⟨ H ⟩. We also used the relations, δhop ≡ khop(D) /khop(H) ~ kO + D/kO + H ~ kOH + D/kOH + H, which should be valid because of the much higher hopping rates of atomic H and D compared to those of atomic O and OH radical. The energy barrier difference against hopping between atomic D and H is ~10 K on amorphous water ice (Hama et al. 2012). When the right hand side of Eq. (D.5) is less (more) than zero at given fD, s, the ratio fHDO, s decreases (increases) on a timescale shorter than bOHkph. Since the photodissociation timescale in gas with low extinction is much shorter than the dynamical timescale (>1 Myr in the case of our cloud formation model), let us assume steady state, which corresponds to the minimum (or maximum) of fHDO, s/fD, s. Then we obtain (D.6)To produce water ice mantles under UV irradiation, the rate of OH formation through the hydrogenation of atomic oxygen (i.e., kO + H ⟨ O ⟩ ⟨ H ⟩) should be larger than the photodesorption rate of H2O (see Fig. 9). When the dynamical timescale is longer than the timescale of the OH formation, the chemistry evolves as the OH formation rate and the H2O photodesorption rate are almost balanced: (D.7)This relation is confirmed in our simulations. The population of OH can be evaluated from the following equation: (D.8)Then (D.9)where we used kph ⟨ H2O ⟩ ≫ kphdes ⟨ H2O ⟩ ≈ kO + H ⟨ O ⟩ ⟨ H ⟩. From Eqs. (D.6), (D.7), and (D.9), we get the scaling relation, (D.10)where Pphdes is kphdes/kph ~ 0.02 (Andersson et al. 2006; Arasa et al. 2015). The term δhopfD, s/ (1 + Γ) corresponds to pOH → HDO which is discussed in the main text. In the case with Γ ≪ 1 (or, in other words, OH + H is the dominant formation pathway of H2O ice), we get fHDO, s ≈ fD, s. In another extreme case Γ ≫ 1, we get fHDO, s ≈ (Γ-1 + 0.02)fD, s, i.e., fHDO, s is much smaller than fD, s. For example, Γ is ~104 at AV = 1 mag in our fiducial simulation, and it further increases with the increase of AV. The minimum of the fHDO, s/fD, s ratio is ~0.02 in our fiducial simulation, which is well reproduced by Eq. (D.10).

Appendix E: Comparisons with observations of sulfur-bearing species

As discussed in Sect. 5.1, the evolution of the OPR(H2) depends on the assumed heavy metal abundances. The goal of this appendix is to determine which case, the LM abundances (fiducial case) or the HM abundances, gives the better predictions on the OPR(H2) evolution in the ISM.

There are some estimates of OPR(H2) in cold dense clouds/cores from observations of molecules other than H2 in the literature. We do not use them for the constraint here, because how to estimate OPR(H2) from the observations is not well-established, and because the evolution of OPR(H2) in the ISM can be in the non-equilibrium manner and may vary among sources depending on their past physical evolution. Instead, we compare observations of sulfur-bearing molecules toward diffuse/translucent/molecular clouds in the literature with our model results. We focus on the total abundance of selected sulfur-bearing molecules in the gas phase (CS, SO, and H2S) and the HCS+/CS abundance ratio. The former can probe the total abundance of elemental sulfur in the gas phase, though the main form of gaseous sulfur is likely to be the neutral atom S or S+ especially in gas with low extinction. The latter, the HCS+/CS ratio, can probe the ionization degree of the gas, which is related to the S+ abundance (and abundances of the other heavy metal ions). The HCS+/CS ratio is anticorrelated with the electron abundance in our models, because CS is formed by dissociative electron recombination of HCS+ and destroyed by photodissociation. Although the HCS+/CS ratio also depends on the photodissociation timescale of CS, it is common between the fiducial model and model HM at given AV. The rate coefficient and the branching ratios for the electron recombination of HCS+ are taken from Montaigne et al. (2005).

Turner (1996) derived the molecular abundances of H2S, CS, and SO in translucent clouds with line of sight visual extinction of up to 5 mag. He found that the total abundance of these three species with respect to hydrogen nuclei is typically 10-8–10-7. In the dense molecular clouds, TMC-1 and L134N, with higher line of sight visual extinction than the samples in Turner (1996), the total abundance is 10-9–10-8 (Ohishi et al. 1992; Dickens et al. 2000). The lower total abundance in the dense molecular clouds implies that depletion of gaseous sulfur is significant in the dense clouds (Joseph et al. 1986). It is not obvious which AV in our model can be compared with the observations towards the translucent and dense molecular clouds. Towards TMC-1(CP) and L134N, there is evidence of CO freeze-out; the gaseous CO abundance with respect to hydrogen nuclei (4 × 10-5) is less than the canonical value of 10-4 and infrared absorption by CO ice is detected in the line of sight towards the vicinity of them (Whittet 2007, 2013). We assume that our results at AV> 2 mag, where the gaseous CO abundance decreases to less than ~4 × 10-5 because of the freeze-out, can be compared with the observations in TMC-1 and L134N, while the results at AV< 2 mag are compared with the observations in clouds with lower line of site extinction.

Figure 11 shows the total abundance of H2S, CS, and SO as a function of visual extinction in our models. Model HM better reproduces the total abundance of the S-bearing molecules observed in the translucent clouds, while the fiducial model better reproduces the total abundance observed in the dense molecular clouds. This result implies that again, the depletion of sulfur from the gas phase occurs during the formation and evolution of molecular clouds, and that model HM underestimates the degree of the sulfur depletion.

The HCS+/CS ratio has been derived toward diffuse clouds (0.08, Lucas & Liszt 2002), clouds with line of sight visual extinction of up to 5 mag (0.01−0.1, Turner 1996), and the dense molecular clouds, TMC-1 and L134N (0.06, Ohishi et al. 1992). Roughly speaking, the ratio is in the range of 0.01−0.1 in the diffuse and dense ISM. Figure 12 shows the HCS+/CS ratio and the electron abundance in the fiducial model and model HM. The fiducial model reproduces the HCS+/CS ratio much better than model HM, though the both models tend to underestimate the HCS+/CS ratio compared with observations, i.e., overestimate the ionization degree of the gas.

Considering the HCS+/CS ratio is reproduced by the fiducial model much better than by model HM, we conclude that the fiducial model gives a better prediction for the S-bearing species and thus for the OPR(H2). In particular, model HM seems to overpredict the S-bearing species in the gas phase. The non-success of the model HM probably means that a non-negligible fraction of sulfur is incorporated into dust grains (Keller et al. 2002) and/or the non-thermal desorption rates of icy sulfur are overestimated in the current model; the partitioning of elemental sulfur between gas and ice in our model depends on the assumed non-thermal desorption rates, especially chemisorption probabilities, which remain uncertain at the current stage. In our models, the dominant sulfur reservoirs in ices are HS and H2S, though there has been no detection of H2S ice in the ISM. HS ice is formed from H2S ice via the hydrogen abstraction reaction, H2S + H, with an activation energy barrier of 860 K (Hasegawa et al. 1992). HS ice is hydrogenated to form H2S ice again, desorbing ~1% of the product H2S through chemisorption. This loop is the main mechanism of the desorption of icy sulfur in our model as in Garrod et al. (2007). If we set the chemisorption probability for the hydrogenation of HS ice to be 0.1%, only 3% of elemental sulfur remains in the gas phase at AV = 3 mag in model HM, though gaseous sulfur is still more abundant than in the fiducial model.

All Tables

Table 1

Initial abundances with respect to hydrogen nuclei.

In the text
Table 2

Summary of free chemical parameters.

In the text

All Figures

thumbnail Fig. 1

Threshold abundances of CO and electrons to be the dominant reactant of H2D+. If CO and electron abundances with respect to H2 are below the lines, Reaction (1) in the backward direction dominates over H2D+ destruction by CO and electrons. See Appendix A for more details.
In the text
thumbnail Fig. 2

Physical evolution after passing through the shock front. We adopt the same shock model as Model 4 of B04 and H10; i.e., the preshock density, temperature and velocity are 10 cm-3, 40 K, and 15 km s-1 respectively. The inset shows the zoom-in view of the earliest evolution at AV ≤ 0.1 mag, which corresponds to t ≤ 0.3 Myr.
In the text
thumbnail Fig. 3

Efficiency factor for reaction OH + H2→ H2O + H, , without (yellow) and with the reaction-diffusion competition in our models as functions of dust temperature. For the case with the competition, with different Edes(H2) are shown by different colors.
In the text
thumbnail Fig. 4

Fractional abundances of selected species with respect to hydrogen nuclei in the fiducial model as functions of visual extinction, or time. The solid lines represent species in the bulk ice mantle, while the dashed lines represent gaseous species. The H2 and H atom abundances are multiplied by a factor of 10-3.
In the text
thumbnail Fig. 5

The number of icy layers in total (top) and fractional compositions in the active surface ice layers (bottom) as functions of visual extinction, or time.
In the text
thumbnail Fig. 6

OPR(H2) in the fiducial model as a function of visual extinction. The dashed line shows the steady-state value given by Eq. (14) using molecular abundances in the simulations, while the dotted line shows the low-temperature thermalized value, 9exp(−170.5 /Tgas).
In the text
thumbnail Fig. 7

Deuteration of selected species and OPR(H2) in the fiducial model as functions of visual extinction, or time. The solid lines represent icy species in the active surface ice layers, while the dashed lines represent gaseous species. The HDO / H2O ratio in the bulk ice mantle is shown by the thick solid line.
In the text
thumbnail Fig. 8

Normalized destruction rate of H2D+ by the reactions with CO and HD, recombination with electrons, and the backward reaction with H2 in the fiducial model as functions of AV. The normalized destruction rates by H2 assuming the thermalized ortho to para ratio of H2 are also shown for comparison.
In the text
thumbnail Fig. 9

Important reaction routes to form H2O and HDO on the (icy) grain surface in our models. The gray area indicates the cycle of photodissociation and reformation of water ice, which is discussed in the text.
In the text
thumbnail Fig. 10

Steady-state fHDO, s/fD, s ratio in the active ice layers as a function of the ratio of the reaction rate between OH + H2→ H2O + H and OH + H → H2O at three temperatures given by Eq. (D.10). See Appendix D for more details.
In the text
thumbnail Fig. 11

Fraction of elemental sulfur in the form of S+ (dashed lines), neutral atom S (dotted lines), and S-bearing molecules (the sum of CS, SO, and H2S; solid lines) in the gas phase as a function of AV in the fiducial model (blue lines) and model HM (red lines).
In the text
thumbnail Fig. 12

HCS+/CS abundance ratio (solid lines) and electron abundance (dashed lines). Red lines indicate the model with the high metal abundances, while blue lines indicate the fiducial model with the low metal abundances. See Table 1 and Appendix E.
In the text
thumbnail Fig. 13

OPR(H2) (solid lines) and atomic D/H ratio (dotted lines) in the gas phase as functions of visual extinction. Red lines indicate the high metal case, while blue lines indicate the low metal case. The dashed lines indicate the steady-state values of the OPR(H2) given by Eq. (14).
In the text
thumbnail Fig. 14

HDO / H2O ratio in the bulk ice in the two grids (with and without the reaction-diffusion competition) of models as functions of AV. For the case without the competition (yellow), the solid line represents the model with Nact = 4, Edes(H2) = 550 K, and the low-metal abundances. The areas represent the models including variations in Nact and Edes(H2), while the light-colored areas represent the models including variations in Nact, Edes(H2), and the initial abundances. For the case with the competition, the models with different Edes(H2) are shown by different colors. The solid lines represent the models with Nact = 4. The areas represent models including the variations in Nact, while the light-colored areas represent the models including the variations in Nact and the initial abundances. The vertical dashed and dotted lines indicate the regions where the fraction of elemental oxygen in the atomic form becomes less than 10% and 1%, respectively, in the models represented by the solid lines.
In the text
thumbnail Fig. 15

Results from a weaker shock model plotted vs visual extinction to be compared with analogous figures for our fiducial model. a) Physical evolution after passing through the shock front (see Fig. 2). b) Fractional abundances of selected species with respect to hydrogen nuclei (see Fig. 4). c) The OPR(H2) (see Fig. 6). d) Deuteration ratios of selected species (see Fig. 7).
In the text
thumbnail Fig. C.1

OPR(H2) when the conversion of hydrogen into H2 is almost complete under the steady-state assumption. Above the thick dashed line, where τo → p is smaller than τH2, the steady-state assumption is justified. Above the thin dashed line, where β1 is larger than β2, the OPR(H2) is similar to the low-temperature thermalized value of 9exp(−170.5 /Tgas). See Appendix C.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.