THE D/H RATIO OF WATER ICE AT LOW TEMPERATURES

Jeong-Eun Lee; Edwin A. Bergin

doi:10.1088/0004-637X/799/1/104

1. INTRODUCTION

The fractionation of deuterium is an important chemical tracer of the physical history of the dense interstellar medium. The initiating reaction for the chemistry ( ${\rm H_3^+ + HD} \rightleftharpoons {\rm H_2D^+ + H_2}$ ) is exothermic in the forward direction, requiring additional energy of ∼232 K to activate the backward channel. As such, large deuterium enrichments in gaseous species have been isolated at low (T_gas ∼ 10 K) temperatures in dense (n $_{H_2}$ >10⁴ cm⁻³) star-forming gas (e.g., DCO⁺, DCN; Millar et al. 1989; Roueff & Gerin 2003). A key facet of the gas-phase chemistry is that many of the deuterium fractionation reactions produce an excess of deuterium atoms relative to hydrogen atoms. Due to high mobility, atoms are key players in the surface chemistry of molecular ices. This leads to the prediction that molecular ices formed at low temperatures will be enriched in deuterium (Tielens 1983). Indeed high levels of D-fractionation of water have long been observed to be present in warm (T ≳ 100 K) gas near young massive stars (Jacq et al. 1990). At these temperatures, the D/H ratio of gaseous molecules should be commensurate with the interstellar D/H ratio of ∼3 × 10⁻⁵ (Linsky et al. 2006) and the D-enriched water is therefore believed to originate from ices that formed at cold (T_gas = T_dust ∼ 10–30 K) temperatures prior to the formation of the star.

Historically, the measured water D/H enrichment in massive hot cores ranged from ∼20–100 above the galactic value (Jacq et al. 1990; Gensheimer et al. 1996; Helmich et al. 1996). The enriched D/H ratio of water has also been detected in the warm inner envelopes of low-mass young stellar objects (YSOs; Stark et al. 2004; Parise et al. 2005). The study of water in the interstellar medium has recently been given a large boost from the successful launch and operations of the Herschel Space Observatory (Pilbratt et al. 2010). Herschel has observed water vapor in numerous star-forming regions (van Dishoeck et al. 2011), which can be supplemented with ground-based or Herschel observations of HDO to derive a D/H ratio. Some of the Herschel results are finding high D/H ratios of water. For instance, Bergin et al. (2010) detected HD¹⁸O in Orion KL which implies a water D/H ratio of ∼0.003 (Neill et al. 2013). Toward some low-mass protostars, in both resolved and unresolved data, high levels of enrichments are found with HDO/H₂O > 0.003–0.05 (Liu et al. 2011; Coutens et al. 2012; Taquet et al. 2013). However, in some instances, lower values are found; for example, Jørgensen & van Dishoeck (2010) set an upper limit of HDO/H₂O ∼ 6 × 10⁻⁴ in one source while Persson et al. (2013) estimate ∼9 × 10⁻⁴ in yet another. Thus, a range of enrichment levels are now being inferred to have been present in cold water ice. This implied diversity calls for a theoretical exploration of how deuterium is implanted into ices.

Over the past few years there have been two advances in our theoretical understanding of gas-phase deuterium fractionation in the interstellar medium that have an impact on the composition of both gas and ices. First, we have recognized that the deuterium fractionation does not stop with H₂D⁺ as the sole initiating ion in the chemistry. Rather, successive reactions with HD can form additional deuterated isotopologues of H $_3^+$ :

$\begin{equation*} \rm {H_2D^+ \rightleftharpoons ^{HD}_{H_2} D_2H^+ \rightleftharpoons ^{HD}_{H_2} D_3^+} \end{equation*}$

(each is exothermic in the forward direction with a barrier for the reverse channel; Vastel et al. 2004; Roberts et al. 2003). This has important implications for the primary charge carriers and also the products of deuterium chemistry, which can tend to multiply deuterated forms (Roberts et al. 2003). Because HD and H₂ do not freeze in appreciable abundance on grain surfaces, some of these molecules (e.g., H₂D⁺ and D₂H⁺) have been isolated as tracers of the gas physics in the heavy-element, freeze-out-dominated dense centers of pre-stellar cores (Flower et al. 2004; van der Tak et al. 2005).

Second, it has been noted that deuterium fractionation reactions have an innate dependence on the ortho-to-para ratio (opr) of H₂ (Gerlich et al. 2002) and also on the relative abundances of the high- and low-energy spin isomers of each of the reactants (H₂D⁺, D₂H⁺, D $_3^+$ , H $_3^+$ ; Flower et al. 2004; Hugo et al. 2009). Thus, the presence of even a small amount of ortho-H₂, which has excess internal energy, can aid in overcoming the barrier for the reverse channel in the above reactions. This will reduce fractionation at the start of the chain (Walmsley et al. 2004, hereafter WFP04), while the relative spin state abundances of the other reactants can also change the resulting distribution of fractionated products.

In our work, we wish to explore the deuterium fractionation that becomes implanted within ices in the context of a cold and dense core and its dependence on the H₂ opr, the gas temperature, and the gas density. According to Flower et al. (2006), the opr for key species (H₂D⁺, D₂H⁺, D $_3^+$ , H $_3^+$ , H₂) varies during the core formation process (i.e., as the density grows). Furthermore, there is an additional dependence on the dynamical evolution (e.g., free-fall or steady-state contraction). In this work, we do not adopt any dynamical model for the core evolution. We rather consider the physical conditions (density and temperature) as free, but fixed, parameters for a given model. In this sense, a higher density can be considered as representative for later core evolutionary stages. We also consider the H₂ opr as a free parameter since it is unlikely in the steady state or equilibrium state (LTE) unless the timescale is long enough (see Figure 1 of Flower et al. 2006). Finally, in this study, the chemical evolution is calculated at a given physical condition with pre-calculated oprs.

In Section 2, we will present our model of the deuterium chemistry, while Section 3 presents the results showing variable levels of water-ice deuterium fractionation can be created depending on the H₂ opr and the gas temperature as well as the H₂ density. This can significantly alter the HD abundance and the composition of major charge carriers in the dense core center. In Section 4, we discuss the implications of our result.

2. MODEL

We adopt the chemical code used in Bergin et al. (1995), Bergin & Langer (1997), and Lee et al. (2004). The code deals with the basic gas chemistry (photodissociation, photoionization, and cosmic ray ionization) and gas-grain interactions (adsorption and desorption). The binding energies of H, O, OH, and H₂O are adopted from Cuppen & Herbst (2007), and they are 650, 800, 3500, and 5640 K, respectively, for the water-coated grain surfaces. For the case of bare silicate grain surfaces, these binding energies are divided by 1.47 (Lee et al. 2004). The sticking coefficient is assumed to be 1.0 for all species including H and D. According to Chaabouni et al. (2012), the sticking coefficients of H and D are 0.86 and 0.95, respectively, to a bare silicate grain surface. For water-ice-coated grains, the sticking coefficients must be greater. Therefore, we assume unity for the sticking coefficients of H and D. Three desorption mechanisms are included: thermal evaporation, cosmic-ray-induced heating, and direct photodesorption. We have assumed the H₂ cosmic-ray ionization rate (ζ) of 6 × 10⁻¹⁷ s⁻¹, which results in the equilibrium number density of atomic hydrogen of 2.5 cm⁻³ in our model; these numbers are consistent with the result of Goldsmith & Li (2005). We have updated the code to include the ortho-to-para-dependent deuterium chemistry in the gas and the surface chemistry for the formation of H₂ and water ice as described in the subsections below. Our chemical network includes ∼170 species (including grain counterparts and deuterium counterparts) and ∼1200 reactions.

Motivated by Hassel et al. (2010, hereafter HHB10), we assume the core is evolving from the lower density molecular cloud with most hydrogen as molecular, carbon in CO, and some initial water ice mantle present. In reality, a dense core that we model in this work must have evolved from lower density material. Chemical abundances in a forming core also change with time and density. However, most studies of dense core chemical evolution assume atomic initial conditions (e.g., C⁺, N, O, Fe⁺, etc.) with pre-existing H₂. HHB10 modeled the chemical evolution as the molecular cloud (the progenitor of the core) forms in a postshock region of atomic gas and suggested that their final chemical abundances could be a reliable initial condition when the chemical evolution is further calculated in dense cores. We adopt their results of Model 1 to set the initial abundances for our fiducial model. The most important initial condition in our work is the oxygen abundances in the atomic gas and the water ice because a relative fraction of atomic oxygen in the gas will determine the final abundance of newly formed (deuterated) water ice. This effect will be presented in Section 3.1.

First, we assume that all carbon at the start of our calculations is in the form of gaseous CO with an abundance of 2.8 × 10⁻⁴ (Lacy et al. 1994; Lee et al. 2003). We also assume that half of total oxygen is in CO. The remaining half of oxygen is found as water ice or gas phase atomic oxygen. Based on Model 1 of HHB10, the abundance ratio between atomic oxygen and water ice is about 2. Therefore, the abundances of atomic oxygen and water ice for our fiducial model are 1.9 × 10⁻⁴ and 0.9 × 10⁻⁴, respectively. We assume the initial atomic hydrogen abundance of 1 × 10⁻⁵ because the H density is about 1 cm⁻³ in molecular clouds (Chang et al. 2007) and n(H₂) =10⁵ cm⁻³ in our fiducial model. (However, the initial H abundance does not affect the final results of our model.) The initial abundances for our fiducial model are listed in Table 1.

Table 1. Initial Abundances

Species	Abundance (Relative to H₂)
H	1.00 × 10⁻⁵
He	0.18
He⁺	1.33 × 10⁻¹⁰
CO	2.80 × 10⁻⁴
HD	2.80 × 10⁻⁵
H₂O (ice)	9.00 × 10⁻⁵
O	1.90 × 10⁻⁴
N	4.28 × 10⁻⁵
SI⁺	3.40 × 10⁻⁹
Mg⁺	2.20 × 10⁻⁹
S⁺	5.60 × 10⁻⁷
Fe⁺	3.40 × 10⁻⁹

Download table as: ASCII Typeset image

2.1. Ortho-to-para-dependent Deuterium Chemistry in the Gas

The approach adopted here is to assume that the oprs of key elements (H $_3^+$ , H₂D⁺, D₂H⁺, and ${\rm D_3^+}$ ) are in equilibrium at the gas temperature. However, the H₂ opr is a free parameter. We then adopt the state to state rate coefficients for each reaction calculated by Hugo et al. (2009) modified by the spin state fractional amount for both reactants. For example, we take the following reaction:

$\begin{equation} {\rm o\!-\!H_2D^+ + o\!-\!H_2 \rightarrow p\!-\!H_3^+ + HD}. \end{equation} \tag{ 1 }$

Hugo et al. (2009) determine a rate of k₂ = 4.67 × 10⁻¹¹ exp(0.82/T) cm³ s⁻¹. In our network, we calculate the equilibrium value of the opr of H₂D⁺ for a given temperature while the H₂ opr is provided as an input. Then, we multiply the associated reaction rate by these fractions. In our example above, the modified rate is $k_2^{\prime } = k_2 \times f_{\rm o-H_2D^+} \times f_{\rm o-H_2}$ . We do not track the oprs of the products because we a priori assume the ratios are in equilibrium (or fixed for H₂). As a result, the rate coefficient for the reaction of o − H₂D⁺ + o − H₂ becomes 1.78 × 10⁻¹⁷ at 10 K with the H₂ opr of 7 × 10⁻⁵, regardless of the spin type of the product (i.e., o–H $_{3}^+$ and p–H $_{3}^+$ are not considered separately in the product). The rate coefficients ( $k_2^{\prime }$ ) derived from the state-dependent rates of Hugo et al. (2009) at 10, 20, and 30 K are listed in Table 2.

Table 2. Rate Coefficients

Reactions	Rate Coefficients
Reactions	10 K	20 K	30 K
p- ${\rm HD}_2^+$ + p-D₂ → D $_{3}^+$ + HD	2.61E-15	2.17E-12	1.63E-11
p- ${\rm HD}_2^+$ + o-D₂ → D $_{3}^+$ + HD	1.18E-11	1.32E-10	2.37E-10
o- ${\rm HD}_2^+$ + p-D₂ → D $_{3}^+$ + HD	2.65E-13	1.73E-11	6.01E-11
o- ${\rm HD}_2^+$ + o-D₂ → D $_{3}^+$ + HD	1.13E-09	9.97E-10	8.25E-10
m-D $_{3}^+$ + HD → ${\rm HD}_2^+$ + D₂	5.32E-17	9.95E-14	1.07E-12
o-D $_{3}^+$ + HD → ${\rm HD}_2^+$ + D₂	7.46E-17	1.35E-13	1.69E-12
o-D $_{3}^+$ + p-H₂ → H₂D⁺ + D₂	2.23E-25	5.52E-18	1.48E-15
o-D $_{3}^+$ + p-H₂ → ${\rm HD}_2^+$ + HD	8.77E-20	1.14E-14	5.61E-13
o-D $_{3}^+$ + o-H₂ → H₂D⁺ + D₂	7.67E-26	7.70E-19	9.03E-16
o-D $_{3}^+$ + o-H₂ → ${\rm HD}_2^+$ + HD	2.23E-16	1.27E-13	5.86E-12
p-H₂D⁺ + p-D₂ → ${\rm HD}_2^+$ + HD	3.28E-13	2.09E-11	5.94E-11
p-H₂D⁺ + p-D₂ → D $_{3}^+$ + H₂	4.81E-14	2.87E-12	8.02E-12
p-H₂D⁺ + o-D₂ → ${\rm HD}_2^+$ + HD	9.75E-10	8.43E-10	5.74E-10
p-H₂D⁺ + o-D₂ → D $_{3}^+$ + H₂	4.05E-10	3.37E-10	2.26E-10
o-H₂D⁺ + p-D₂ → ${\rm HD}_2^+$ + HD	6.19E-16	3.18E-12	3.55E-11
o-H₂D⁺ + p-D₂ → D $_{3}^+$ + H₂	7.36E-17	3.82E-13	4.28E-12
o-H₂D⁺ + o-D₂ → ${\rm HD}_2^+$ + HD	2.06E-12	1.46E-10	3.93E-10
o-H₂D⁺ + o-D₂ → D $_{3}^+$ + H₂	5.14E-13	3.37E-11	8.78E-11
p- ${\rm HD}_2^+$ + HD → H₂D⁺ + D₂	3.97E-16	8.49E-14	5.11E-13
p- ${\rm HD}_2^+$ + HD → D $_{3}^+$ + H₂	5.36E-12	6.27E-11	1.20E-10
o- ${\rm HD}_2^+$ + HD → H₂D⁺ + D₂	4.28E-16	8.84E-14	5.43E-13
o- ${\rm HD}_2^+$ + HD → D $_{3}^+$ + H₂	7.47E-10	6.58E-10	5.75E-10
m-D $_{3}^+$ + p-H₂ → H₂D⁺ + D₂	1.68E-25	4.27E-18	1.03E-15
m-D $_{3}^+$ + p-H₂ → ${\rm HD}_2^+$ + HD	3.90E-20	4.96E-15	2.04E-13
m-D $_{3}^+$ + o-H₂ → H₂D⁺ + D₂	4.92E-26	5.30E-19	5.56E-16
m-D $_{3}^+$ + o-H₂ → ${\rm HD}_2^+$ + HD	1.48E-16	8.25E-14	3.23E-12
o- ${\rm HD}_2^+$ + o-H₂ → H $_{3}^+$ + D₂	2.56E-23	4.05E-18	1.21E-15
o- ${\rm HD}_2^+$ + o-H₂ → H₂D⁺ + HD	1.19E-15	6.50E-14	1.62E-12
o- ${\rm HD}_2^+$ + p-H₂ → H $_{3}^+$ + D₂	1.19E-25	1.20E-18	2.31E-16
o- ${\rm HD}_2^+$ + p-H₂ → H₂D⁺ + HD	1.29E-18	1.91E-14	4.35E-13
p- ${\rm HD}_2^+$ + o-H₂ → H $_{3}^+$ + D₂	1.29E-27	4.62E-20	7.28E-17
p- ${\rm HD}_2^+$ + o-H₂ → H₂D⁺ + HD	4.62E-17	2.24E-14	1.11E-12
p- ${\rm HD}_2^+$ + p-H₂ → H $_{3}^+$ + D₂	7.83E-28	4.65E-19	3.21E-16
p- ${\rm HD}_2^+$ + p-H₂ → H₂D⁺ + HD	3.71E-18	4.12E-14	7.87E-13
o-H₂D⁺ + HD → H $_{3}^+$ + D₂	2.25E-17	1.28E-13	1.71E-12
o-H₂D⁺ + HD → ${\rm HD}_2^+$ + H₂	2.03E-12	1.41E-10	3.97E-10
p-H₂D⁺ + HD → H $_{3}^+$ + D₂	3.91E-18	5.33E-15	4.59E-14
p-H₂D⁺ + HD → ${\rm HD}_2^+$ + H₂	1.35E-09	1.20E-09	8.78E-10
p-H $_{3}^+$ + HD → H₂D⁺ + H₂	1.40E-09	1.15E-09	1.00E-09
o-H $_{3}^+$ + HD → H₂D⁺ + H₂	1.04E-10	4.17E-10	5.96E-10
p-H₂D⁺ + p-H₂ → H $_{3}^+$ + HD	3.76E-20	2.64E-15	7.98E-14
p-H₂D⁺ + o-H₂ → H $_{3}^+$ + HD	7.96E-17	1.39E-13	7.78E-12
o-H₂D⁺ + p-H₂ → H $_{3}^+$ + HD	2.04E-19	1.74E-14	5.36E-13
o-H₂D⁺ + o-H₂ → H $_{3}^+$ + HD	1.78E-17	3.96E-14	2.01E-12
p-H $_{3}^+$ + p-D₂ → H₂D⁺ + HD	2.25E-12	1.25E-10	4.08E-10
p-H $_{3}^+$ + p-D₂ → ${\rm HD}_2^+$ + H₂	1.77E-13	1.01E-11	3.32E-11
p-H $_{3}^+$ + o-D₂ → H₂D⁺ + HD	4.91E-10	3.80E-10	3.00E-10
p-H $_{3}^+$ + o-D₂ → ${\rm HD}_2^+$ + H₂	9.65E-10	7.26E-10	5.68E-10
o-H $_{3}^+$ + p-D₂ → H₂D⁺ + HD	1.54E-14	4.45E-12	2.50E-11
o-H $_{3}^+$ + p-D₂ → ${\rm HD}_2^+$ + H₂	1.46E-14	4.20E-12	2.36E-11
o-H $_{3}^+$ + o-D₂ → H₂D⁺ + HD	3.05E-11	1.35E-10	1.89E-10
o-H $_{3}^+$ + o-D₂ → ${\rm HD}_2^+$ + H₂	7.85E-11	2.94E-10	3.89E-10

Download table as: ASCII Typeset image

We estimate the oprs or meta-to-ortho ratios using the energy levels given in the CDMS database (Müller et al. 2005) and the spin statistics given in Table II of Hugo et al. (2009). The ratio is then calculated by

$\begin{equation} \frac{\rm ortho}{\rm para} = \frac{2I_{\rm o} + 1}{2I_{\rm p} + 1} \times \frac{\sum (2J_{\rm o} + 1)\; {\rm exp}(-E_{\rm o}/kT)}{\sum (2J_{\rm p} + 1)\; {\rm exp}(-E_{\rm p}/kT)}. \end{equation} \tag{ 2 }$

We list the oprs and meta-to-ortho ratios at 10, 20, and 30 K at Table 3.

Table 3. Spin Isotopomer Ratios

Species	10 K	20 K	30 K
H₂ o/p	7.00 × 10⁻⁵	1.79 × 10⁻³	3.06 × 10⁻²
D₂ o/p	3.62 × 10³	49.1	11.7
H $_3^+$ o/p	7.53 × 10⁻²	0.388	0.663
H₂D⁺ o/p	1.85 × 10⁻³	0.152	0.599
D₂H⁺ o/p	97.7	7.47	3.43
D $_3^+$ o/m	2.30 × 10⁻²	0.244	0.581

Download table as: ASCII Typeset image

In tracking the various ratios, the opr of H₂ remains the most important as it, along with temperature, must be below certain values to allow fractionation to proceed. Although it is a free parameter in our model, we also derive the equilibrium value of opr-H₂, adopting the usual relation of opr-H₂ = 9 × exp(−170.5/T) (Le Bourlot 1991) for below 80 K. Furthermore, WFP04 have shown that the H₂ opr does not reach the equilibrium value (opr-H₂ =3.5 × 10⁻⁷) at 10 K. Rather, over a wide range of densities, the value reached is instead opr-H₂ ∼7 × 10⁻⁵. We therefore have fixed the lower bound of the H₂ opr to the latter value. Flower et al. (2004) explored the temperature dependence of the opr for other key molecular species (H $_3^+$ , H₂D⁺, D₂H⁺, ${\rm D_3^+}$ ) and find that above ∼15 K the equilibrium ratio is a good approximation. For gas temperatures below ∼15 K, the ratios are generally above the equilibrium value. Outside of H₂, which is the most important player, we do not account for this effect.

2.2. Grain Surface Chemistry

In our calculation, two types of grain surface chemistry (formation of H₂ HD, and D₂; formation of OH, OD, H₂O, and HDO ices) are included. For the surface chemistry, we adopt the method used in Fogel et al. (2011), which followed the treatment of Hollenbach et al. (2009). For the formation of H₂ and water, an H or O atom must be frozen on the grain surface before another atom finds it. Therefore, the reaction rate coefficients for the surface chemistry vary with time depending on the abundance ratio between the ice species and grain. For example, if the abundance of H atoms on grains (H ice) is greater than the grain abundance, the reaction rate is reduced by the ratio of the grain abundance to the abundance of H ice because a gaseous H atom must collide with a grain to initiate the H₂ formation reaction. We scale the rate coefficient of each reaction in Fogel et al. (2011) by the mass ratio between the deuterium atom and the hydrogen atom to treat the formation of HD and deuterated water ice. According to Kristensen et al. (2011), the adsorption energies for D₂ and HD of 72.0 and 68.7 meV. The ratio of these binding energies is only 1.05, which is very small considering the mass ratio (1.3) between D₂ and HD. The mass added to the total mass of a molecule by D is a small fraction, so we assume the same binding energies for both deuterated and normal species, except for D and H, which have the binding energy ratio of 1.2 (Perets et al. 2007).

3. RESULTS

3.1. Deuterium Fractionation at T_gas = 10 K and Low H₂ opr

Observations of gas-phase molecular depletion, which is dominated by freeze-out onto grain surfaces, have suggested that perhaps all heavy elements are frozen on grain surfaces (see Bergin & Tafalla 2007 for a more complete discussion). Under this circumstance, the only remaining gas-phase species with a dipole moment that are present are HD, H₂D⁺, and D₂H⁺. Of these, HD will not emit appreciably at 10 K, thus observations and several models focused on the question of formation and presence of the two ions in the centers of dense cores (see WFP04). Below we will compare our results to WFP04 in order to ascertain the similarities/differences between WFP04 and our (similar) chemical model and to isolate the reason of the similarities/differences.

Figure 1(a) shows the equilibrium abundances of the key ions as functions of density at the same physical condition as the standard model of WFP04 (a_g = 0.1 μm, T =10 K, ζ = 3 × 10⁻¹⁷ s⁻¹) with the complete depletion of heavy elements. The level of ionization is comparable between the two models, but we see major differences in the most abundant D-bearing ion. In our models ${\rm D_3^+}$ is the most abundant ion for all densities, while WFP04 find a changing distribution switching from H⁺ (low density) to D $_3^{+}$ (high density) as seen in Figure 2 of WFP04. Identical differences are found in comparison to the work of Flower et al. (2004) by Sipilä et al. (2010). They show that in adopting state-to-state coefficients of Hugo et al. (2009) there is a ∼factor of 5 difference in key fractionation reactions (see Sipilä et al. 2010). Our results are not entirely identical with Sipilä et al. (2010). For example, at n(H₂) = 10⁵ cm⁻³, H $_3^+$ is the main ion in our network compared to H⁺ in theirs. This is likely due to the fact that we treat our formation of H₂ differently by correcting for the possibility that there may be less than one hydrogen per grain (Fogel et al. 2011). Indeed, H⁺ becomes the main ion if we turn off our H₂ formation reaction. In addition, Sipilä et al. (2013) showed that H $_3^+$ becomes the main ion if the surface chemistry is explicitly included to their calculations (see Figure 3 of Sipilä et al. 2013).

In Figure 1(b), the effect of the surface chemistry of water ice on the gas-phase deuterium fractionation is presented. In this figure, we plot the abundances at a timescale of 10⁷ yr for all densities. The timescale for the chemistry to reach equilibrium varies with the density, i.e., a longer timescale for a lower density. However, the timescale in our density range is much shorter than 10⁷ yr. Therefore, our models reach similar conditions as WFP04 at 10⁷ yr. Earlier studies suggested the D-bearing ions as main tracers of dense and cold regions (Roberts et al. 2004; Aikawa et al. 2005). However, as shown in Figure 1(b), the abundances of deuterated ions drop at high densities when the surface chemistry of water ice formation is included. This is due to a decline in the abundance of HD as fractionation progresses and D atoms are locked in ices. In our calculations for n(H₂) ⩾10⁶ cm⁻³, the depletion of deuterated ions is somewhat exaggerated because our model does not consider other surface reactions that can produce HD (Sipilä et al. 2013).

Figure 2(a) illustrates the effect of HD depletion, which has been also characterized in Sipilä et al. (2013). Deuterium atoms, which are produced in the gas by the dissociative recombination of the deuterated ions, become frozen onto grain surfaces to form deuterated water ice, resulting in the reduction of the HD abundance in the gas. Therefore, the formation of water ice on grain surfaces significantly affects the gas-phase deuterium fractionation and changes the main charge carriers. The sharp increase of gas H and D abundances occurs around 10⁵ yr because the H–H₂ chemistry reaches initial equilibrium at the time. Before equilibrium, at low temperatures (T <20 K), the atomic hydrogen is quickly frozen on grain surfaces to form water ice as soon as it is dissociated from other species.

Figure 2(b) shows the evolution in the abundance of (deuterated) water ice and the D/H ratio in water ice for identical conditions as adopted for the results show in Figure 2(a). Both H₂O and HDO abundances increase with time, with a steeper gradient in HDO resulting in the increase of the D/H ratio with time. The steeper gradient leading to HDO ice formation is a direct result of the rise in the D/H ratio in the atomic pool (Figure 2(a)). The calculated maximum D/H ratio of water ice is about 0.05, which should be considered as the upper limit, because we only consider the water surface chemistry. In a more expansive model, the deuterium atoms could be incorporated into other deuterated ices (Tielens 1983); however, water ice still dominates.

Sipilä et al. (2013) also used the rates of Hugo et al. (2009) and incorporated the surface chemistry explicitly into their calculation, so it makes a good comparison with our work. We adopted their initial abundances to calculate the chemical evolution in the model with n(H₂) = 10⁵ cm⁻³ and T_gas = 10 K, which are similar to the conditions of the outermost shell in their model core. We also used their cosmic ray ionization rate (ζ = 1.3 × 10⁻¹⁷ s⁻¹), but our constant oprs and meta-to-ortho ratios at 10 K are adopted for the calculation.

In the calculation of Sipilä et al. (2013), the oprs vary with time; however, except for H₂, the variation in opr is within an order of magnitude. This does not produce a notable change in the resulting abundances of ${\rm H_3^+}$ and its deuterated ions as well as the H₂O/HDO ices. Figure 3 can be compared with Figure 4 (right panel) of Sipilä et al. (2013) directly. The overall trends are very similar although the abundances reach the steady state after ∼10⁶ yr because only water ice formation is included in our calculation. The biggest difference between this test and our fiducial model is the D/H ratio of water; in this test, it reaches 0.13, about three times higher than the D/H ratio in the fiducial model. This is because all oxygen in Sipilä et al. (2013) is initially atomic and available to eventually be incorporated into HDO and H₂O on grain surfaces.

**Figure 3.** Abundances (relative to H₂) of selected species as a function of time for comparison with Sipilä et al. (2013). The physical conditions of the model are T_gas = 10 K and n(H₂) = 10⁵ cm⁻³, and the same initial abundances and cosmic ray ionization rate as used in Sipilä et al. (2013) were adopted.
Download figure:
Standard image High-resolution image

**Figure 4.** Effect of the initial water ice abundance. f(H₂O_ice) is the initial ratio between the water ice abundance and the total oxygen abundance excluding oxygen in CO (f(H₂O_ice) = X(H₂O_ice)/[X(H₂O_ice)+X(O_atom)]). In our fiducial model, f(H₂O_ice) = 0.3. The timescale for the plot is 3 × 10⁵ yr. The solid and dotted lines present the D/H ratio of water ice when the gas temperature is 10 K and 20 K, respectively. The density is the same as 10⁵ cm⁻³ for all models. If more oxygen is initially locked in water ice, the final D/H ratio of water ice becomes lower. Only for the model of f(H₂O_ice) = 0, i.e., with no initial water ice, do we adopt the binding energies of molecules to the bare silicate grain, which are lower than those to the water-ice-coated grain by a factor of 1.47.
Download figure:
Standard image High-resolution image

We tested different initial abundance ratios between water ice and gaseous atomic oxygen at two temperatures, 10 and 20 K. The equilibrium H₂ opr at each temperature was adopted. The result is presented in Figure 4. The D/H ratio of water ice at a timescale of 3 × 10⁵ yr depends on the relative disposition of water ice and atomic oxygen at the start of our calculation. For a higher initial abundance of water ice relative to atomic oxygen, we find a lower water D/H ratio on grain surfaces. This is because if more oxygen is initially locked in water ice, less oxygen is available to form H₂O and HDO on grain surfaces at later times, resulting in a reduced water D/H ratio. In addition, because of the active backward reactions of deuterium fractionation at the higher temperature of 20 K, the overall D/H ratio of water ice is lower than that at 10 K.

Recently, Taquet et al. (2013) developed a full deuterium chemical network, adopting explicit surface chemical reactions for the formation of various molecular ices as well as water ice, to show that other molecular ices such as CH₃OH and H₂CO can be more deuterated than the water ice. However, the water ice is the most abundant ice on grain surfaces, with an abundance higher than other ices by orders of magnitude in most conditions as shown in Taquet et al. (2013). Therefore, the biggest reservoir of deuterium in grain surfaces is water ice, which justifies our focus on water surface chemistry. Figure 5 shows the D/H ratio of water ice versus the density, which has been calculated with the same parameters used in Figure 10 of Taquet et al. (2013), i.e., the H₂ opr of 3 × 10⁻⁶ and A_v = 10 mag. For the calculation, we adopted our fiducial initial abundances. The results of our calculations are very consistent with those by Taquet et al. (2013), and the D/H ratio of water vapor observed in the hot corino of IRAS 16293 (Coutens et al. 2012) is well explained by our model with this H₂ opr. Dislaire et al. (2012) determined that the H₂ opr is 1 × 10⁻³ in the dense core of IRAS 16293. With this ratio we can also explain the observed D/H ratio water as presented dotted lines in Figure 5; if the temperature is 20 K, the density of IRAS 16293 must be higher than 10⁵ cm⁻³.

**Figure 5.** D/H ratio of water ice vs. the density at 3 × 10⁵ yr. The same parameters (H₂ *opr* = 3 × 10⁻⁶ and *A_v* = 10 mag) as used in Figure 10 of Taquet et al. (2013) have been adopted (solid lines). The blue and red lines indicate the results at 10 K and 20 K, respectively. The green box refers to the D/H ratio of water observed in the hot corino of IRAS 16293 (Coutens et al. 2012). The dotted lines are for H₂ *opr* = 1 × 10⁻³.
Download figure:
Standard image High-resolution image

3.2. Deuterium Fractionation for Variable Gas Temperature, H₂ opr, and Density

The D/H ratio of water ice is a function of the efficiency of deuterium fractionation in the gas, which in turn depends on the H₂ opr and the gas temperature. In the dense ISM, there exist a range of environments where the ambient dense star-forming material has gas temperatures of ∼10 K (e.g., Taurus), but other regions have higher values (e.g., Orion). Figure 6 shows the two-dimensional dependence for the D/H ratio of water ice and HD from T_gas = 10–30 K and for H₂ opr from ∼10⁻⁵ to the high temperature ratio of 3. In Figure 6, the fractionation ratio is calculated at the density of 10⁵ cm⁻³ with our fiducial initial abundances at the timescale of 3 × 10⁵ yr. In this calculation, the H₂ opr is set to be a constant as a free parameter, independent of temperature. The white dotted line in Figure 6 indicates the equilibrium H₂ opr calculated by the relation, 9 × exp(− 170.5/T). The equilibrium value becomes lower than 7 × 10⁻⁵ at T ≲ 15 K. In the case, we plot 7 × 10⁻⁵. The D/H ratios along the white lines are the values when H₂ opr reaches the equilibrium at given temperatures.

Clearly, in Figure 6 there is a strong dependence on the D/H ratio for water ice on both parameters. However, the D/H ratio is insensitive to the H₂ opr when the H₂ opr is smaller than ∼10⁻³–10⁻². Flower et al. (2006) showed that the H₂ opr decreases with collapse, while we simply explore the effects of fractionation for constant H₂ opr. However, the D/H ratio of water ice at a given temperature is not sensitive to the H₂ opr as long as it is smaller than 10⁻², which is true at all times in the Flower et al. collapse model. Moreover, the formation of the molecular cloud itself will precede the formation of the dense core. Thus, there is potentially ∼10⁶–10⁷ yr for the gas to evolve to a lower H₂ opr ratio (and setting the stage for subsequent D enrichments) prior to the condensation and collapse of a dense molecular core (e.g., Bergin et al. 2004; Clark et al. 2012).

Over this range the ratio of HD/H₂ shows an inverse dependence (Figure 6, bottom). That is, the D/H ratio of molecular hydrogen increases with the gas temperature and H₂ opr, trending toward the cosmic value. When H₂ opr is <10⁻² and T_gas < 20 K, the abundance of HD exhibits depletions. The lowest HD/H₂ ratio from this calculation at n(H₂) = 10⁵ cm⁻³ is ∼7.8 × 10⁻⁶, which is lower than the comic ratio (∼3 × 10⁻⁵) by a factor of about 4. The fraction is even lower (∼2 × 10⁻⁷) by more than two orders of magnitude compared to the cosmic ratio at n(H₂) = 3 × 10⁵ cm⁻³.

Figure 7 presents the abundance of ${\rm H_3^+}$ and its deuterated isotopologues across the identical two-dimensional grid. ${\rm H_3^+}$ is always the most dominant ion. For a temperature lower than ∼20 K, the freeze-out of deuterium as ices produces low abundance ratios of deuterated ions relative to ${\rm H_3^+}$ . Furthermore, the ratios drop sharply at higher temperatures due to the activation of the backward reactions. In our fiducial model, the grain mantle is coated with water ice, so the binding energies of molecules are greater compared to those to the bare silicate grain. As a result, even at 30 K, most CO, which is a primary destroyer of H $_3^+$ and H₂D⁺ via the reactions of ${\rm H_3^+ + CO \rightarrow HCO^+ + H_2}$ and H₂D⁺ + CO → DCO⁺ + H₂, is still frozen on grain surfaces, leaving H $_3^+$ abundant. However, in the model with T_gas = 25 K with no initial water ice (i.e., with binding energies to the bare silicate grain), CO evaporates and the abundance of H $_3^+$ is lowered to 10⁻⁹.

Figure 8 shows the dependence of the D/H ratio of water ice and the H₂ gas on the temperature and density when the equilibrium H₂ opr, which is marked as the white dashed lines in Figures 6 and 7, is assumed at each given temperatures. In the densities greater than ∼3 × 10⁵ cm⁻³, the D/H ratio is only dependent on the gas temperature. The highest D/H ratio of water ice is almost 0.1 for densities greater than 3 × 10⁵ cm⁻³ and T ⩽ 15 K. As noted above, this value must be considered as a maximum ratio since we considered only the water ice formation and did not include other potential deuterium carriers (e.g., H₂CO, CH₃OH, and NH₃). However, water ice is still the most abundant ice and the major reservoir of deuterium in the icy grain mantle.

Figure 9 presents the abundance distribution of H $_3^+$ and its deuterated species in the same domain of temperature and density as used in Figure 8. H $_3^+$ is not sensitive to temperature as seen in Figure 9 and predominantly exhibits a density dependence due to the increase recombination rates. The abundance of the deuterated daughter products of H $_3^+$ are sensitive to density at temperatures lower than 20 K. At gas temperatures greater than 20 K, their abundances drop sharply and are almost constant. This is due to the activation of the backward reactions which slows down deuterium fractionation for T > 20 K. The effective formation of the deuterated water ice on grain surfaces at higher densities results in the significant depletion of HD in the gas (Figure 8), which in turn results in the heavy depletion of deuterated ions at the cold and densest regions (Figure 9). Figure 10 shows the abundance ratio between H₂D⁺ and H $_3^+$ , which presents a similar trend to the explained above; the abundance ratio changes very sharply around 3 × 10⁵ cm⁻³ at T ⩽ 20 K, but at T > 20 K, the abundance ratio drops with temperature and is insensitive to density.

**Figure 10.** D/H ratio (log (D/H)) of H $_3^+$ derived from Figures 9(a) and (b). The contours represent log (D/H) of −1.3, −1.8, −2.3, and −2.8, which are the same as used in Figure 9(a). The D/H ratio of H $_3^+$ becomes very low at the density greater than 3 × 10⁵ cm⁻³ and the temperature lower than 20 K although it is still higher than the cosmic value.
Download figure:
Standard image High-resolution image

**Figure 10.** D/H ratio (log (D/H)) of H $_3^+$ derived from Figures 9(a) and (b). The contours represent log (D/H) of −1.3, −1.8, −2.3, and −2.8, which are the same as used in Figure 9(a). The D/H ratio of H $_3^+$ becomes very low at the density greater than 3 × 10⁵ cm⁻³ and the temperature lower than 20 K although it is still higher than the cosmic value.
Download figure:
Standard image High-resolution image

4. DISCUSSION AND IMPLICATIONS

It is clear from numerous investigations discussed above that the deuterium chemistry has an innate dependence on the ortho to para ratio of H₂, the gas temperature, and the gas density. In this paper, we have attempted to provide an initial look at the interdependence of these key parameters in terms of the overall initiation of deuterium chemistry in the gas, but also the deuterium implantation into ices. We note that we are not alone in exploring this parameter space and we also refer the reader to Taquet et al. (2013) and Sipilä et al. (2013). Below we attempt to use observations to provide additional perspective.

4.1. Water Ice

A wide range of D/H ratios of water vapor have been reported in the hot gas near protostars from <6 × 10⁻⁴ to >0.01 (Jacq et al. 1990; Gensheimer et al. 1996; Helmich et al. 1996; Stark et al. 2004; Parise et al. 2005; Bergin et al. 2010; Jørgensen & van Dishoeck 2010; Liu et al. 2011; Coutens et al. 2012; Neill et al. 2013; Persson et al. 2013). These values are believed to trace the deuterium enrichment of evaporated water ice which formed during earlier colder phases. In addition, Coutens et al. (2012) argued that the water should form before the gravitational collapse of a protostar based on their calculation of nearly constant water D/H ratios (∼0.01) from the hot core to the low-density outer envelope of IRAS 16293-2422. Thus, one interpretation of the D/H ratios is that these variations reflect differences in the initial conditions, in particular the H₂ opr, gas temperature, and density (and also cosmic ray ionization rate), under which the ices form. Thus, to account for this range, one possibility is that the H₂ opr could vary from source to source or even along the line of sight, while assuming that the gas generally remains cold (T <15 K). It is generally believed that H₂ forms on grains with an opr of 3, which is then gradually thermalized toward the ratio at the local gas temperature via reactions with H⁺ and H $_3^+$ . The interconversion of ortho and para-H₂ is slow (requiring timescales >3 × 10⁶ yr for steady state; Flower et al. 2006). However, Pagani et al. (2009, 2011) demonstrate that the conversion is dependent on the local gas density with order of magnitude variations in the H₂ opr. Thus, differences in the density evolution during the core formation phase can result in variations in the opr of H₂. Clearly, this is complicated by the fact that the gas temperature is also not constant and can vary from source to source.

If we focus on the extremes the high enrichment levels >0.01, such as toward the low-mass protostar IRAS 16293 (Coutens et al. 2012), require both low H₂ opr <0.01 and low gas temperatures (<20 K) in our models. Moreover, the gas would need to stay within these conditions for ∼10⁵ yr, which appears plausible given current timescales for core formation (see discussion in Pagani et al. 2011). Low ratios, such as found toward the low-mass protostar NGC 1333 IRAS 4B (Jørgensen & van Dishoeck 2010), cannot be achieved unless the gas temperature during much of the evolution is >15 K while also requiring an H₂ opr > 0.1. Thus, our models, with a different set of parameters, can match the observed range. However, both of these sources are in low-mass star-forming regions where the ambient gas is cold with characteristic temperatures of ∼10–15 K (Bergin & Tafalla 2007). Moreover, it would require an H₂ opr that varies by over two orders of magnitude, despite the roughly similar conditions of ambient material. In this light, perhaps not all measured D/H ratios are intrinsic (i.e., set at birth) and some mechanisms could be operative to alter the ratio in the very hot gas near protostars (see discussion in Neill et al. 2013). However, as suggested by Jørgensen & van Dishoeck (2010), higher resolution observations can help in providing a more meaningful comparison between the various measurements.

We note that attempts to detect solid HDO in YSOs have been also made (Dartois et al. 2003; Parise et al. 2003) setting an upper limit of the D/H ratio of ≲ 0.02 in water ice. Considering the parameters associated with the surface chemistry such as binding energies, these limes are consistent with our predictions. The highest equilibrium D/H ratio of water ice in our fiducial model is ∼0.05 at 10 K and 10⁵ cm⁻³. This value is a firm upper limit because we have a limited gas-phase network that explores water ice formation. Thus, we do not consider hydrogenation of other ices such as H₂CO, CH₃OH, and NH₃.

4.2. Ions in the Cold Gas

As summarized in the review of dense cold cores by Bergin & Tafalla (2007), dense pre-stellar molecular cores often exhibit a high degree of central concentration. Furthermore, most of their mass resides at a low temperature of ∼10 K. Under these conditions, chemical models have previously predicted that H₂D⁺, D₂H⁺, and ${\rm D_3^+}$ would be excellent tracers of the dense core center where many other molecules will be frozen onto grain surfaces (Roberts et al. 2004; Aikawa et al. 2005). The dense protoplanetary disk midplane presents similar cold dense environment and Ceccarelli & Dominik (2005) predict that the deuterated forms of ${\rm H_3^+}$ will be the primary charge carriers. However, these studies considered only the depletion of molecules without dealing with actual surface chemistry. Because HD is not frozen on grain surfaces, in the condition of the high degree of molecular depletion, the deuterated ions would be the main charge carriers in the gas.

However, D atoms, produced by the dissociative recombination of deuterated ions, are also frozen on grains surfaces to eventually be incorporated into molecular ices. Although our study considered only water ice formation, it can catch the representative effect of the depletion of deuterium, leaving ${\rm H_3^+}$ as the primary species in the gas. If our results are correct, the H₂D⁺ abundance of 10⁻¹⁰ ∼ 10⁻⁹ derived by Caselli et al. (2003) might correspond to the region with the density of 1 to 3 × 10⁵ cm⁻³ rather than ∼10⁶ cm⁻³, where almost all deuterium is converted to the deuterated water ice as seen in Figure 1(b) (X(H₂D⁺) ⩽10⁻¹²). Thus, molecular ions like H₂D⁺,D₂H⁺, ${\rm D_3^+}$ would not trace core centers or the dense midplane of protoplanetary disks, but rather dense intermediate layers that are closer to the core or disk surface where HD can exist with higher abundance. According to Caselli et al. (2008), which surveyed dense starless/protostellar cloud cores in the o − H₂D⁺ (1_{1, 0}–1_{1, 1}) line, the column density of o − H₂D⁺ is not constrained solely by one parameter such as CO depletion factor or density. As presented in this study, the abundances of ions are functions of density, temperature, H₂ opr, and the gas-grain interaction.

Nevertheless, HD is significantly destroyed only at n(H₂) ⩾3 × 10⁵ cm⁻³ and T < 20 K (Figure 8). For instance, the abundance of H₂D⁺ is about 2 × 10⁻⁹ at n(H₂) = 10⁶ cm⁻³ and T = 20 K (Figure 9). Therefore, deuterium bearing ions will still trace the densest gas at T ⩾ 20 K.

We thank the anonymous referee, whose comments led to improvements in the paper. This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education of the Korean government (grant No. NRF-2012R1A1A2044689). This work was also supported by the National Science Foundation under grant No. 1008800 and by NASA under grant NN08AH23G from the Astrophysics Theory and Fundamental Physics and Origins of Solar System programs.

THE D/H RATIO OF WATER ICE AT LOW TEMPERATURES

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. MODEL