NEW EXTENDED DEUTERIUM FRACTIONATION MODEL: ASSESSMENT AT DENSE ISM CONDITIONS AND SENSITIVITY ANALYSIS

In this work we are primarily concerned with providing a new extended deuterium network and assessing its capacity to model chemistry under both static cold and warm conditions in the ISM. We do not consider evolutionary models such as those pertaining to low-mass and high-mass star formation separately, although the static conditions we consider arise from evolutionary processes. The evolution of the ISM begins from fragmentation of turbulent, mainly atomic clouds with kinetic temperatures up to ∼100 K and densities of ∼10–100 cm⁻³. The denser clumps evolve into starless molecular cores (Bergin & Tafalla 2007) with temperatures of 8–15 K and densities of ∼10⁴–10⁶ cm⁻³ (Snow & McCall 2006; André et al. 2009; Launhardt et al. 2010; Nielbock et al. 2012). Some of these gravitationally bound cores may begin contracting, first isothermally, and then with increasing internal densities and temperatures. Then a central hydrostatic object forms, which starts heating up the surrounding gas. Protostars with a mass greater than 8 solar masses are generally referred to as high-mass protostellar objects (HMPOs). The collapsing envelope material can then reach temperatures of several hundred Kelvin closer to the central star, and peak densities of ∼10⁸ cm⁻³, conditions which define a "hot core" or, for low-mass protostars, a hot "corino" (van Dishoeck 2006). In this paper we concentrate on the evolutionary stages ranging from a cold molecular cloud to the warm envelopes of protostars. We choose a wide parameter space with a grid of 1000 points covering temperatures between 5 and 150 K and densities of 10³–10¹⁰ cm⁻³, and assume the standard ISM dust and a fixed A_V = 10 mag, meaning that the photochemistry is only driven by secondary UV photons. Assuming a fixed A_V reduces the problem to two dimensions, which is easier to analyze and visualize.

We have utilized the gas–grain chemical model "ALCHEMIC" developed by Semenov et al. (2010), where a detailed description of the code and its performance is presented. The code is optimized for modeling the time-dependent evolution of large chemical networks, including both gas-phase and surface species. In this paper we added a large set of reactions involving deuterated species. A few features of the "ALCHEMIC" model are summarized below.

The self-shielding of H₂ from photodissociation was calculated using Equation (37) from Draine & Bertoldi (1996). The shielding of CO by dust grains, H₂, and its self-shielding was calculated using the precomputed table of Lee et al. (1996, Table 11). We consider cosmic-ray particles (CRP) as the only external ionizing source, using a CRP ionization rate for H₂, ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹ (Herbst & Klemperer 1973), appropriate for molecular cloud environments and which has been utilized in several previous studies (such as Wakelam et al. 2006; Vasyunin et al. 2008; Druard & Wakelam 2012). The gas–grain interactions include dissociative recombination and neutralization of ions on charged grains, sticking of neutral species and electrons to uniformly sized 0.1 μm dust grains with a sticking coefficient of 1 and release of ices by thermal, CRP-, and UV-induced desorption, such that at high temperatures the surface population will be low as thermal desorption takes over. We do not allow H₂ and its isotopologues to stick to grains. We assume a UV photodesorption yield of 10⁻³ (e.g., Öberg et al. 2009a, 2009b). With our fixed visual extinction, the photon field derives from secondary electron excitation of molecular hydrogen followed by fluorescence.

We assume that each 0.1 μm spherical silicate grain provides ≈1.88 × 10⁶ sites (Biham et al. 2001) for surface recombination that proceeds solely through the classical Langmuir–Hinshelwood mechanism (e.g., Hasegawa et al. 1992). The grain surface topology, the presence of high- and low-energy binding sites, grain sizes and shapes are all separate parameters that may severely impact the chemistry. An accurate study of this impact will require a detailed treatment of the microscopic physics of molecules on various solid surfaces, which is far beyond the scope of the present study. For further reading we recommend papers by Perets & Biham (2006), Cuppen et al. (2009), and Vasyunin & Herbst (2013a, 2013b), where some of these issues have already been addressed.

Upon a surface recombination, we assume there is a 5% probability for the products to leave the grain due to the conversion of some of the exothermicity of reaction into breaking the surface-adsorbate bond (Garrod et al. 2007). We do not find significant differences (less than a factor of two) in D/H ratios and abundances of essential species, such as H₃⁺, water, and ammonia from varying this probability between 1% and 10%. However, we found a significant variation in ice abundances of formaldehyde and methanol of up to a factor of 6 at lower temperatures (≲ 25 K) when considering higher desorption probabilities ≳5%. Interestingly, ice abundances increase with the desorption probability and we find that this is due to a much more efficient formation of formaldehyde and methanol at intermediate times. Due to more intense gas–grain interactions precursor species of H₂CO and CH₃OH are able to form more readily in the gas phase and later stick to grains. Consequently, formaldehyde and methanol are formed faster via surface processes. Following experimental studies on the formation of molecular hydrogen on dust grains by Katz et al. (1999), we adopt the standard rate equation approach to the surface and ice chemistry without quantum-mechanical tunneling through the potential walls of the surface sites. We also do not consider competition kinetics between activation and diffusive barriers (Garrod & Pauly 2011).

A typical run, with relative and absolute tolerances of 10⁻⁵ and 10⁻²⁵, utilizing the original gas–grain network without deuterium chemistry (∼7000 reactions, ∼700 species) takes 1–5 s for 1 Myr of evolution with a Xeon 3.0 GHz CPU. With our new, almost ten-fold larger deuterium network, the same run takes approximately an order of magnitude longer to calculate. The linear dependence of the CPU time versus species number in the model is due to the advanced numerical scheme implemented in the ALCHEMIC code, which generates and uses sparse Jacobi matrices.

As input data, reaction rate coefficients and physical properties need to be specified, as do initial abundances. We have chosen to implement two different initial abundance sets and calculate the chemical evolution with the new deuterium network for 1 Myr.

For the first set, hereafter referred to as the "Primordial" model, we utilized the "low metals" abundances of Graedel et al. (1982) and Lee et al. (1998). Initially all deuterium is located in HD, with D/H =1.5 × 10⁻⁵ (Stancil et al. 1998; Linsky 2003); see Table 1. The abundances in the second set, the "Evolution" model, were calculated with our deuterium chemistry model, assuming a TMC1-like environment: T = 10 K and n_H = 10⁴ cm⁻³, at t = 1 Myr. Under such conditions elemental deuterium from HD is efficiently redistributed to other molecules, leading to their high initial D/H fractionation. These final abundances at 1 Myr are used as input in the Evolution model (see Table 2). The Evolution model serves as a simple example of a two-stage chemical model with physical conditions that can change dramatically at 1 Myr, unless the evolutionary model is run strictly under TMC-1 conditions.

As a first step toward creating a consistent network with deuterium fractionation, we undertook a thorough search in the literature for updates to the reaction rates of the original non-deuterated network. We utilized the latest osu.2009 gas-phase chemical network and incorporated all essential updates as of 2012 December, adopted from Horn et al. (2004), Chabot et al. (2010), Hamberg et al. (2010a), Wakelam et al. (2010b), Laas et al. (2011), as well as those reported in the KInetic Database for Astrochemistry (KIDA).⁶ Further, to allow for the synthesis of the few complex molecules in our network such as methanol (CH₃OH), methyl formate (HCOOCH₃), and dimethyl ether (CH₃OCH₃), an extended list of surface reactions and photodissociation of ices was adopted from Garrod & Herbst (2006). Several tens of gas-phase photoreaction rates were updated using the new calculations of van Dishoeck et al. (2006).⁷

Next, we applied a cloning routine to this updated network (as described in Rodgers & Millar 1996), and added all additional primal isotope exchange reactions for H₃⁺ as well as CH$_3^+$ and C₂H$_{2}^+$ from Roberts & Millar (2000b), Gerlich et al. (2002), Roberts et al. (2004), and Roueff et al. (2005). In this cloning routine all reactions bearing hydrogen atoms are considered to have deuterated analogs, and "cloned" accordingly (assuming the same rate coefficient if no laboratory data are available). In cases where the position of the deuterium atom is ambiguous, we apply a statistical branching approach. In the resulting network we do not yet distinguish between the ortho/para states of molecules, and leave this for a separate paper.

A typical example of the outcome of the cloning procedure is presented for the reaction between C⁺ and CH₃:

$$\begin{equation} {\rm C^+ + CH_3} \rightarrow {\rm C_2H^+ + H_2} \Rightarrow \left\lbrace\begin{array}{@{}l}\rm C^+ + CHD_2 \rightarrow {\rm C_2D^+ + HD} \\ {\rm C^+ + CHD_2} \rightarrow {\rm C_2H^+ + D_2} \end{array} \right. \end{equation} \tag{ 1 }$$

The single ion–molecule reaction of C⁺ with CH₃ is cloned into two separate channels for CHD₂. Moreover, the branching of these two new channels is not equal. To visualize this, we label the two deuterium atoms in CHD₂ D_a and D_b. For the first reaction, which forms C₂D⁺, D_a can be placed on either product and D_b can be placed on the other, hence we have two possibilities: C₂D_a⁺ + HD_b or C₂D_b⁺ + HD_a. For the second channel, which forms C₂H⁺, both deuterons have to be placed on D₂ and we only have one possibility. This analysis assumes that the deuterons on D₂ are indistinguishable, which is in agreement with the Pauli exclusion principle. Alternatively, we could initially assume that they are distinguishable, but because half of the D₂ rotational-nuclear spin states are missing, the simple argument about 2/3 and 1/3 branching ratios remains valid.

To limit the size of the network we have restricted the cloning process to avoid any -OH endgroups. Observations of deuterated species suggest that fractionation of species with -OD endgroups is less important in low-mass protostars, but may still be important for high-mass protostars. For example, Parise et al. (2006) conducted a survey of deuterated formaldehyde and methanol in a sample of seven low-mass class 0 protostars, and found CH₃OD/CH₂DOH ≲ 0.1. A hypothesis of rapid conversion of CH₃OD into CH₃OH in the gas-phase due to protonation reactions that would affect only species for which deuterium is bound to the electronegative oxygen has been suggested by Charnley et al. (1997) and Osamura et al. (2004). We conducted a small study using a version of our deuterium network where -OH endgroups were cloned and found no significant changes in the resulting time-dependent molecular abundances.

Added and updated reactions are found in Tables 12 and 13. The resulting chemical network consists of ∼55,000 reactions connected by ≳ 1900 species, to our knowledge the most extended network for deuterium chemistry to date.⁸

Given the large size of the network with uncertainties in the adopted rate coefficients, reaction barriers, and branching ratios, it is educational to estimate how these uncertainties propagate in time-dependent modeled abundances. Before performing a detailed sensitivity analysis, as in our study of disk chemistry uncertainties (Vasyunin et al. 2008), it is of interest to characterize the influence of the reaction rate updates on the calculated abundances and the D/H ratios. This may help us to highlight the significance of recent laboratory astrochemistry activities, both for deuterated and un-deuterated species, in providing more accurate astrochemical data to the community.

First, we studied the effects of introducing deuterium chemistry into our model on abundances of un-deuterated species by comparing abundances throughout the parameter space to a non-deuterated version of the network. We found that species with relative abundance >10⁻²⁵ show mean values in abundance variations between the two networks within a factor 0.95–1.05. Since we did not find any particularly large variations in abundances for H-bearing species, we conclude that the results from our updated analysis are a pure effect of updated reaction rates and not caused by the additional pathways created by the cloning routine.

In order to separate the effect that recent updates have had on abundances, we generated an additional network by cloning an outdated network restricted to the reaction rate updates up until 2005. We then studied the impact of updated reaction rate coefficients by comparing the calculated time-dependent abundances between the "old" chemical network and the "new" network in the two-dimensional parameter space discussed in Section 2.1. In addition to the D/H abundance ratios, we will emphasize the differences in these ratios between the models, which we calculated by dividing the respective D/H ratios in the updated 2012 network by those from the outdated 2005 network, and will denote this ratio as R(D/H). The results have been obtained with the Primordial model only. In this comparison, we have excluded minor species with relative abundances below 10⁻²⁵. It should be noted that the R(D/H) ratios may remain unchanged when absolute abundances of species and their isotopologues in the updated and outdated networks increase in unison.

We list in Table 3 the arithmetic mean value calculated over the parameter grid (T = 5–150 K, n_H = 10³–10¹⁰ cm⁻³) of R(D/H) ratios for all species with fractional abundances ⩾10⁻²⁵ for which the mean value of R(D/H) ratios have changed by more than a factor of five. Among the listed species, we find light hydrocarbons (e.g., CHD₃, CH₂D₂), ions (CH₄D⁺, HD₂O⁺), and simple organic molecules (e.g., DCOOH, D₂CO), as well as key molecules such as doubly deuterated water and ammonia. Multi-deuterated species appear to be more affected by the updates than their singly deuterated analogs, as there are more intermediate pathways involved in their chemistry, as is most evident by comparing HD₂O⁺ and H₂DO⁺ in Table 3.

There are also several species affected by the abundances that do no show any variance in D/H ratios, i.e., both deuterated and undeuterated species are similarly affected by updates. The abundances of CH₂D⁺ and CHD₂⁺ provide good examples of such behavior. These species show a coherent increase (within a factor of 1.1 between un-deuterated species and isotopologues) in their gas-phase abundances at 1 Myr and at high temperatures (≳ 100 K) and densities (n_H ≳ 10⁷ cm⁻³). As a result, their R(D/H) values remain close to unity. We identified the coherent increase in abundance as originating from an update taken from KIDA in the rate coefficient for the slow radiative association reaction forming CH₅⁺ via CH₃⁺ colliding with H₂. The rate coefficient of this reaction was lowered by almost two orders of magnitude, from 1.30 × 10⁻¹⁴ cm³ s⁻¹ to 4.10 × 10⁻¹⁶ cm³ s⁻¹ (at room temperature; see Wakelam et al. 2010b). We note that the older value was based on a misinterpretation of the original literature, which used 300 K in the formula for the rate coefficient but was intended only for temperatures up to 50 K (Herbst 1985). The same D/H ratio variation is transferred to HD₂O⁺ and H₂DO⁺ through the ion–neutral reaction with CH₅⁺ and its isotopologues reacting with free oxygen atoms. The abundances of the CH₄ isotopologues derive from the dissociative recombination of CH₅⁺ (and its isotopologues) as well as from ion–molecule reactions with CO, so CH₄ is directly affected by the updated reaction rate. It then transfers its D/H ratio variation into C₂H₂D₂ through CH₂D₂ reacting with CH. Another route involves the intermediary reaction between CH₂D₂ and S⁺ to form HD₂CS⁺, which later dissociatively recombines into D₂CS, which in turn reacts with CH₂D₂.

On the other hand, the deuterated analog of this radiative association reaction does not occur; instead CH$_3^+$ + HD produces CH₂D⁺ and H₂, and the corresponding rate constant is the same in the outdated and new networks (Millar et al. 1989). This slows down production of the key ion, CH$_5^+$, while deuterated isotopologues of CH$_3^+$ and CH$_5^+$ are produced with almost the same rate, consequently affecting D/H ratios of the gas-phase species listed in Table 3. It does not affect abundant organic species such as methanol and formaldehyde however, as they are mainly formed by surface hydrogenation of CO ice.

We find that there is a particularly large variance in abundances for the two sulfur-bearing species, C₂S and C₂S⁺, which show an increase in abundance by a factor of 187 and 26, respectively, compared with the non-deuterated network. We are not concentrating on the chemistry of sulfur-bearing molecules in this study because their chemistry is still poorly understood and often restricted to a few pathways. But the additional pathways that the cloning routine generates has a stronger effect on these two species as pathways reducing their abundances proceed much slower than their formation pathways.

We studied the impact of uncertainties in reaction rate coefficients on the resulting chemical abundances, and how they propagate throughout the chemical evolution. Two separate environments were chosen for this study, representing dark clouds (T = 10 K, n_H = 10⁴ cm⁻³; Nielbock et al. 2012) and lukewarm infrared dark clouds (IRDCs; T = 25 K, n_H = 10⁵ cm⁻³; Vasyunina et al. 2012). In these and all subsequent runs mentioned in this paper, we use the Primordial initial conditions unless stated to the contrary. We chose to concentrate on the uncertainties of the rate coefficients of the gas-phase reactions. Including variations in surface chemistry rates is a tricky problem as these depend on surface mobility, binding and diffusion energies of reactants, and properties of the surface itself (porosity, irregularities, etc.). Recently, Taquet et al. (2012) have studied the importance of deuterium fractionation on dust surfaces, using a multi-layered ice mantle model and quantified some of the associated errors in the calculated abundances due to the uncertainties in the surface chemistry, so we do not repeat such a study here.

Our analysis is based on the same method as employed by Vasyunin et al. (2004, 2008). We performed computations for a large set of models, using identical physical conditions and initial abundances, and the same chemical network but with randomly varied rate coefficients within their uncertainty limits. The rate uncertainties were taken from the KIDA database. Most of the reactions with deuterated species were created by the cloning procedure, and hence have unknown uncertainties. For these reactions we used the high standard error value in KIDA,⁹ with a normal logarithmic error distribution of a factor of two.

With this method, we generated 10,000 networks with the new rate coefficients k(T) randomly distributed as follows:

$$\begin{equation} k(T) = {\rm exp} ({\rm ln}\, k_{0}(T) + {\rm ln}\, F \times N[0,1]), \end{equation} \tag{ 2 }$$

where k₀(T) is the measured, calculated, or estimated rate coefficient at temperature T, F is the statistical distribution of the uncertainty, and N[0, 1] is a random value drawn from a standard Gaussian distribution with mean μ = 0 and variance 1. The time needed for a full run consisting of 10,000 networks for a specified temperature and density takes approximately one day of computational time on a Xeon 3.0 GHz CPU (with relative and absolute tolerances of 10⁻⁴ and 10⁻²⁰, respectively).

The huge size of our new deuterium chemical network makes it a challenging task to find a very precise correlation between the rate uncertainties and uncertainties in the molecular abundances of a particular species. Since for many deuterated species the set of primal pathways easily exceeds several reactions, the relative contribution of each individual reaction to the final uncertainty is likely to be small, ≲ 10%. To isolate the most problematic reactions for several key observed species, we conducted a sensitivity analysis using the same cross-correlation method as implemented by Vasyunin et al. (2008). We selected a handful of molecules, including their isotopologues and isomers, viz., H₃⁺, H₂O, HCN, HCO⁺, and CH₃OH, and for each of these species we calculated time-dependent linear correlation coefficients between the abundances and the rate coefficients for all the 10,000 network realizations and for each of our 30 logarithmically taken time steps. The linear correlation coefficients c_L(i, j, t) at specific density, i, and temperature, j, points and at time t are calculated by:

$$\begin{equation} c_{L} (i,j,t) = \frac{\Sigma _{l} \left(x_l^s (i,j,t) - \overline{x^s (i,j,t)}\right) \left(\alpha _l^j - \overline{\alpha ^j}\right) }{ \Sigma _l \left(x_l^s (i,j,t) - \overline{x^s (i,j,t)}\right)^2 \Sigma _l \left(\alpha _l^j - \overline{\alpha ^j}\right)^2 } \end{equation} \tag{ 3 }$$

with $x_s^l (i,j,t)$ being the molecular abundance for species s and iteration l, and $\overline{x_s^l (i,j,t)}$ and $\overline{\alpha ^j}$ signifying the standard (mean) abundances and rate coefficients for species s, respectively. Because key reactions can vary through time evolution, we calculate and use cumulative correlation coefficients in our results, for which we integrated the absolute values of time-dependent linear correlation coefficients over the 30 logarithmic taken time steps taken over 1 Myr. In our results in the next section, we restrict discussion to the cumulative correlation coefficients, designated as c.

We have determined that 10,000 realizations of the network may not be adequate for results of our sensitivity analysis to fully converge. For correlation coefficients c < 0.1 we sometimes still see small variations when we compare our results with a model containing only 9000 realizations. Therefore, we also ran a separate set of simulations with 20,000 realizations for a dark cloud environment, but found the same result with only minor deviations for the reactions with low correlation coefficients, c < 0.1. All correlation coefficients should stop fluctuating with size as soon as the number of the network realizations exceeds the number of reactions in the network, which is ∼50,000. Running the sensitivity analysis code with so many realizations, however, would be prohibitively time consuming. Therefore, below we consider reactions with c > 0.1 in our discussion.

Table 4 lists the subset of reactions with cumulative correlation coefficients c > 0.1, which are the most problematic for the chemical evolution of the following species (and their isotopologues and isomers): H₃⁺, HCO⁺, HCN, H₂O, CH₃OH, H₃O⁺, CH₃⁺, C₂H₂⁺, and CO. Because there are several additional key reactions with 0.05 < c < 0.1, we list them extensively in Table 11, including also reactions with c > 0.05 for the same set of species.

As can clearly be seen, ion–neutral processes dominate Table 4, accompanied by a few dissociative recombination and neutral–neutral reactions as well as cosmic ray ionization of the two critical species: H₂ and He. The last process may require more detailed description in astrochemical models, such as recently presented in Rimmer et al. (2012) and Glassgold et al. (2012), so we assigned a relatively large uncertainty of a factor of two for this group of processes. Approximately half of the bimolecular reactions are connected to the chemical evolution of water and light hydrocarbons in the gas. Abundances of species mostly produced on grains, such as methanol, will not be strongly affected by the uncertainties. However, a small fraction of these species is still present in the gas in the center of dense cores, for which uncertainties in the gas-phase chemistry may become important. Many of the deuterated reactions in the table are produced by our cloning procedure, so their error coefficients are only an approximation, and can, in fact, be larger than estimated. Also, isomerization reactions for HOC⁺ and DOC⁺ with reaction rate uncertainties of a factor of two possess strong correlation coefficients (0.3–0.4).

We find it clear that fractionation channels of the H₃⁺ and CH₃⁺ isotopologues require further study, as do reactions involving the isotopologues of H₃⁺ reacting with CO, water, OH, and their isotopologues, forming the initial steps toward more complex molecules. Reactions with H₃⁺ and H₂D⁺ are both well represented in the list and initiate the ion–molecule chemistry, while HD₂⁺ and D₃⁺ are often not abundant enough to have a significant effect in our models. H₃⁺ and H₂D⁺ react with CO to form the isotopologues of HCO⁺, with OH and OD to form ionized water (H₂O⁺, HDO⁺) as well as the water isotopologues, which strongly affect water abundances and D/H ratios. We also see many other interconnecting reactions among our set of key species. Several dissociative recombination reactions, which proceed very rapidly, show strong correlations. While their reaction rates can be accurately determined (Florescu-Mitchell & Mitchell 2006) the products and branching ratios of these reactions are not precisely known (see, e.g., Hamberg et al. 2010b; Geppert et al. 2006).

After calculating time-dependent abundances for all the realizations of the chemical network with varied reaction rate coefficients, and for each considered physical model, we fit Gaussians to the resulting abundance distributions at 1 Myr for all species. We carefully checked that such an approximation could be applied to the abundance distributions and found it to be the case for almost all species. In Figure 1, we show examples of abundance distributions, with their Gaussian fits, for H₂D⁺ and DCOOH in dark cloud and warm IRDC environments, respectively, both showing good fits. From the Gaussian fits, we then determined the FWHM (2.35σ) of these distributions and used the 1σ values to quantify the spread in abundances, which we henceforth refer to as the abundance uncertainties. We also applied the same procedure for the calculated D/H ratios, and show examples of the resulting D/H ratio distributions in Figure 2 for water and formic acid in the same two regions. We find good fits to all these distributions with estimated 1σ values of factors of 1.9 and 2.9 for the abundance distribution of water and formic acid, respectively, while the values for D/H ratios are a factor of 1.4 and 3.2, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 3, we plot the 1σ abundance uncertainties for deuterated species with up to three D-atoms at 1 Myr as a function of the number of atoms for both environments. There are two major trends visible in this plot. First, the abundance uncertainties are in general lower in the case of the IRDCs models compared with cold dark clouds.

Zoom In Zoom Out Reset image size

At such a low temperature (10–20 K) many reactions with barriers cannot proceed, lowering the overall chemical complexity and thus the cumulative rate uncertainties. One would expect uncertainties to be the lowest for dark clouds, but we suggest that as D/H ratios and abundances of deuterated species are also higher in colder environments, more reactions can occur to increase uncertainties. Second, there is a strong trend of increasing abundance uncertainties with the number of atoms in species. This is obvious as the more atoms a species has, the more reaction pathways lead to its production and destruction from initial composition, and thus the higher is the accumulating effect of their uncertainties on modeled abundances (see also the discussion in Vasyunin et al. 2008).

In general, using the 1σ confidence level, the abundances and column densities of species made of ≲ 3 atoms (e.g., CO, HCO⁺, DCO⁺) are uncertain by factors 1.5–5, those for species made of 4–7 atoms are uncertain by a factor of 1.5–7, and those for more complex species made of >7 atoms are uncertain by a factor of 2–10. The uncertainties for D- and H-bearing species are very similar in dark cloud environments. In warm IRDCs the typical uncertainties of larger H-bearing species (>4 atoms) are approximately a factor of two lower compared to D-bearing species, as de-fractionation begins at these elevated temperatures (25 K). Our estimates for the abundance uncertainties for deuterated species are comparable to the abundance uncertainties of un-deuterated species in protoplanetary disks (Vasyunin et al. 2008), as well as diffuse and dark dense clouds (e.g., Vasyunin et al. 2004; Wakelam et al. 2010a). Chemically simple species containing Mg, Na, and Si tend to have high uncertainties reaching up to three orders of magnitude in abundances because their chemical pathways remain poorly investigated. Also large molecular species, whether they are rather abundant hydrocarbons (C_nH_m with n, m ≳ 4), with fractional abundances up to 10⁻⁹–10⁻⁷, or complex and less abundant organic species (e.g., methyl formate, dimethyl ether), have large error bars in the computed abundances. For this latter group, the uncertainties can reach more than one order of magnitude.

In addition to the abundance uncertainties, we also determine uncertainties in the resulting D/H ratios, as shown in Figure 4. Overall, the uncertainties in D/H ratios are lower when compared with the uncertainties of the corresponding H- and D-bearing isotopologues. This is because abundances of the individual isotopologues are often affected by the rate uncertainties in the same way, given that a majority of the deuterium fractionation processes are cloned, thus inheriting the rate of the "ancestor" reaction. We find the same trends as for the abundance uncertainties, but only a hint of increasing uncertainties with number of atoms. Uncertainties for D/H ratios are generally about half of one order of magnitude, but may vary between a factor of 2 and 10. As in the case of the abundance uncertainties, we find the largest D/H uncertainties for large hydrocarbons (C_nH_m, with n, m ≳ 4), complex organics, and species containing Mg, Na, and Si.

Zoom In Zoom Out Reset image size

The question remains how the uncertainties compare between deuterated and un-deuterated species. To illustrate the overall relative uncertainties between deuterated species and their un-deuterated analogs, we plot in Figure 5 the ratios of abundance uncertainties of up to triply deuterated species and their un-deuterated analogs as a function of the number of atoms. Note that these relative uncertainties are not the same as the uncertainties in D/H ratios; the former can be labeled as u(D)/u(H), while the latter can be labeled as u(D/H), where u stands for uncertainty. Hence, the u(D)/u(H) ratio allows us to compare the relative errors between deuterated and un-deuterated species.

Zoom In Zoom Out Reset image size

A majority of deuterated species show larger abundance uncertainties with respect to their un-deuterated analogs. There are two major reasons for this behavior. First, the majority of reactions with deuterated species originate from our cloning procedure, and thus have larger assumed uncertainties. Second, to produce a deuterium isotopologue of a molecule, additional chemical pathways (e.g., isotope exchange processes) are required, increasing the accumulation of rate uncertainties. For many hydrocarbons (C_nH_m, C_nH_m⁺ with m, n = 2, 3) the abundance uncertainties of their deuterated isotopologues are comparable to those of the main isotopologues (with ratios of ∼0.7–2). These hydrocarbons form through ions of hydrocarbons reacting with H₂ or smaller neutral hydrocarbons, such as CH₄ and C₂H₃, and their reactions originate purely from the cloning procedure.

For a limited number of species, the abundance uncertainties of their deuterated isotopologues are even smaller than for the main isotopologues (with ratios of ∼0.7–0.9; see Figure 5), e.g., C₂D⁺ and D⁻. These are simple radicals and ions produced by a limited set of reactions, with relatively well-known rate coefficients and thus small uncertainties. The relevant deuterium fractionation chemistry is also limited and has comparably low uncertainties (∼0.7–1.1). Note that the spread in abundance uncertainties ratios appears to decrease with increasing number of atoms in species, getting closer to unity. This effect in Figure 5 occurs because we plot only species with relative abundances exceeding 10⁻²⁵ (with respect to hydrogen), whose numbers decrease substantially with size. If we also add species with such low abundances, this feature disappears and the trend in the uncertainties between small and large molecules is similar.

For HD and D₂, we find that the ratio of their abundance uncertainties to that of H₂ can reach very large values, ∼1000 and ∼10,000, respectively. This effect occurs because abundances of H₂ are very well constrained (uncertainties are ∼10⁻⁵), while HD and D₂ have typical values of abundance uncertainties up to one order of magnitude.

Finally, in Figure 6 we plot the abundance uncertainties of deuterated species as a function of number of D-atoms, once again restricted to species with relative abundances >10⁻²⁵. We see an increase in uncertainty with number of D-atoms, as expected from the increasing number of reactions involved for subsequently adding more D-atoms. We also notice a wider spread in uncertainties for species with lower levels of deuteration because smaller species, such as OD, HDO, and DCO⁺, often have lower uncertainties because there are not as many steps involved in their formation compared to larger species. Larger multiply deuterated species cannot have these low uncertainties as there are too many steps involved in their formation. Once again we also find that dark clouds have a larger spread in uncertainties than warm IRDCs. Again we argue that this difference occurs because dark cloud environments overall have higher abundances of deuterated species which are processed through more reactions.

Zoom In Zoom Out Reset image size

In this section, we discuss general trends in the modeled D/H ratios in our two-dimensional parameter space (see Section 2.1). The computed D/H fractionation ratios at the final time of 1 Myr are shown in Figure 7 for the Primordial (left panels of each separate block) and Evolution model (right panels of each separate block), which were introduced in Section 2.3, as functions of density and temperature for the following key gaseous molecules: H₂DO⁺, H₂D⁺, HDO, D₂O, ND, ND₃, DCOOH, and DCN. For most species, the D/H ratios can reach high values of ≳ 10⁻³ at T ≲ 30–80 K. At higher temperatures, ≳ 100 K, the computed D/H ratios begin to approach the elemental ratio of ≈1.5 × 10⁻⁵. The higher D/H ratios of ≳ 0.1 have been observed for many species in the ISM, such as CH₂DOH, D₂CO (Ceccarelli 2002), D₂O (Butner et al. 2007), H₂D⁺ (Caselli et al. 2003), HDO (Liu et al. 2011), NH₂D (Hatchell 2003; Roueff et al. 2005), and NHD₂ (Roueff et al. 2005). For comparison of our model results with observations, see Section 4.1.

Zoom In Zoom Out Reset image size

We isolate species that are either mostly sensitive to kinetic temperature or initial abundances, the former referring to the standard temperature dependence for deuterium fractionation, and the latter referring to the difference between the initial abundances in the primordial and evolutionary models. We note that usually species for which D/H ratios depend more strongly on the initial abundances also show (a weaker) temperature dependence, because the release of CO and other radicals decreases the abundances of the H₃⁺ isotopologues which lowers the efficiency of transferring D atoms to other species. These two groups are listed in Table 5. The temperature-dependent species can be further divided into two subgroups by the temperature where D/H ratios decrease most sharply, and hence where deuterium fractionation becomes less pronounced: (1) at low temperatures, T ∼ 20–40 K (related to the fractionation via H₃⁺ isotopologues) and (2) at higher temperatures T ∼ 80 K (related to the fractionation via CH₂D⁺ and C₂HD⁺).

The left two columns of Figure 7 show species that demonstrate a strong dependence of D/H ratios on the initial abundances while the right two columns show species with a strong temperature gradient for D/H ratios and no distinct dependence on initial abundances. We define here a strong dependence on the initial abundance as an overall variation by a factor of five in D/H ratios at temperatures >100 K for the final time step. For the case of HDO in Figure 7, D/H ratios are approximately two orders of magnitude higher in the Evolution model compared to the Primordial model at temperatures >100 K and hence its D/H ratios are considered to depend strongly on the initial abundance. Frozen molecules show D/H distributions similar to their gas-phase analogs as we do not specifically consider selective substitution of H by D in surface species (yet we do consider production of deuterated molecules by surface chemistry). Therefore, we do not discuss the distribution of D/H for ices separately.

The placement of species in different groups can be explained by the relative pace of their chemical evolution. If steady-state abundances are reached by 1 Myr in the primordial stage, then they remain at steady state throughout the evolutionary model unless the physical conditions are changed. In Figure 8 we show the evolution of D/H ratios at three separate times of 10⁴, 10⁵, and 10⁶ yr, and at a density n_H = 10⁴ cm⁻³ as a function of temperature for H₂D⁺, CH₂D⁺ and HDO. We discuss the specific case for the evolution of HDO later and begin with the two ions (two top panels) which are highly reactive, with fast chemical timescales associated with ion–neutral reactions and dissociative recombinations. Thus, the chemical "memory" of the pristine state of the fast-evolving deuterated species is completely lost at the final considered time of 1 Myr, with no apparent differences between abundances computed with the two distinct initial abundance sets.

Zoom In Zoom Out Reset image size

We find that the majority of neutral species, which includes several H-, C-, O-, N- and S-bearing species, show a dependence on the initial abundance as their evolution is significantly slower compared to ions. These are predominantly neutral species unless they are strongly associated with one of the major ions (such as H₃⁺, CH₃⁺, HCO⁺ and their isotopologues). This is illustrated in Figure 8 for HDO (bottom panel), where we see differences in D/H ratios of ∼3 orders of magnitude between the two models at higher temperatures approaching 100 K. In the Primordial model we start with atomic abundances and the water D/H ratios are lower for ices than for gas vapor because. This is because the formation of key ingredients for HDO ice through OD + H or OH + D is slower gas-phase fractionation via ion–molecule processes. In contrast, the Evolution model starts with high abundances and D/H ratios of water, mostly frozen onto grains. At low temperatures <30 K water quickly reach the steady-state abundances while at higher temperatures, because gas-phase fractionation is not efficient, the D/H ratios of water vapor decreases with increasing temperature. Desorption plays an important role at temperatures above 20–30 K where, even in highly obscured regions, CRP-driven desorption occurs. At temperatures ≳ 100 K we see a drop in D/H ratios down to similar values of water ice because desorption replenishes the water vapor constantly and the water vapor also inherits the same D/H ratios as the ice (∼10⁻²).

H₂DO⁺ was the only ion we found showing a significant difference in the distribution between the two sets of initial abundances, which is related to the protonation of HDO by H₃⁺ and which thus behaves chemically similar to HDO (see above). Complex neutral species, which are not a direct outcome of fast dissociative recombination reactions, evolve more slowly via neutral–neutral and surface processes, and thus show a dependence on the adopted initial abundances. The high initial abundances of deuterated species accumulated at 10 K in the Evolution model enable rich deuterium chemistry during the entire time span of 1 Myr even at warm temperatures of T ≲ 50–100 K. As a result, some species show almost uniformly high D/H ratios of ≳ 10⁻³ at T = 10–150 K (see Figure 7, left column). On the other hand, the same species in the Primordial model (with only HD available initially) show a significant decrease in the final D/H ratios toward warmer temperatures (≳ 40–80 K). The effect is strong for saturated species that are at least partly formed on dust surfaces, and thus are sensitive to the choice of the initial abundances (and temperature), such as water or ammonia isotopologues. Also, most multiply deuterated species, whose formation involves several isotope exchange reactions, show strong dependence on the initial abundances.

Ions constitute the majority of species for which calculated D/H ratios are only dependent on temperature (see Figure 7, right two columns). Species such as the H₃⁺ isotopologues that are sensitive to freeze-out (or specifically, the freeze-out of CO), "daughter" molecules such as DCO⁺, as well as (ionized) light hydrocarbons related to CH₂D⁺ and C₂HD⁺, all show a strong temperature dependence. From the exact gas kinetic temperature, labeled the critical temperature, at which the D/H ratios start to approach the elemental value, we can separate these species into two subgroups. The deuterated species formed via low-temperature fractionation channels involving isotopologues of H₃⁺ have a critical temperature of 20–30 K, whereas other deuterated species synthesized via high-temperature fractionation channels involving CH₂D⁺ and C₂HD⁺ show a higher critical temperature of ∼80 K (Parise et al. 2009). Most multi-deuterated species belong exclusively to the low-temperature group as do the majority of ions, which inherit the temperature dependence from the H₃⁺ isotopologues via proton/deuteron transfer reactions. The best example of such species is H₂D⁺, shown in Figure 7. It shows a D/H turnover point at 20 K, after which the fractionation ratios decrease smoothly and reach levels of ∼10⁻⁵ at temperatures of $\mathord {\sim }$100 K.

There are only a few species that show a dependence on temperature where deuterium fractionation is reduced at T ∼ 80 K. Those are directly synthesized from deuterated light hydrocarbon ions such as CH₂D⁺, which retain high D/H ratios at elevated temperatures, e.g., DCN and CH₂D₂. HDCS is unique in showing no clear difference between the two models at least within our defined temperature range, as its chemical evolution depends almost equally on the two distinct deuterium fractionation routes (via H₃⁺ and light hydrocarbon isotopologues) via the formation of the ion H₂DCS⁺.

Also, we find that the overall sensitivity of the calculated D/H ratios to density is weak, except for species sensitive to the freeze-out of CO and other radicals, i.e., H₂D⁺, HD$_2^+$, and their direct dissociative recombination products. These species experience a rapid decrease in fractionation ratios toward lower densities of <10⁴ cm⁻³, at which CO cannot severely deplete onto the dust grains even after 1 Myr of evolution. However, if we would follow the evolution of our chemical model beyond 1 Myr, at some moment the CO depletion can become severe enough for low-density regions to allow high D/H ratios for the H₃⁺ isotopologues (if the adopted non-thermal desorption rate is not too high).

To understand in more detail the differences between the evolution of deuterated molecules, we have analyzed the chemistry of several key isotopologues. We considered four distinct astrophysical environments: densities of 10⁴ and 10⁸ cm⁻³, and temperatures of 10 and 80 K for both the Primordial and Evolution initial abundance sets. We list key formation and destruction pathways for the assorted species and their main reactants in Tables 14–18 in Appendix B.

In Tables 6–8 we list the most abundant, (potentially) observable deuterated species with ALMA in dark clouds, IRDCs, and HMPOs, respectively. In order to predict the observability of molecules in these different environments we use the Cologne Database of Molecular Spectroscopy (CDMS; Müller et al. 2001, 2005) to calculate line fluxes under local thermal equilibrium (LTE) conditions for various deuterated molecules with transitions observable by ALMA. For these estimates we assume the following parameters for the different environments: for dark clouds temperature 10 K, number density 10⁴ cm⁻³ and H₂ column density 10²² cm⁻² (Launhardt et al. 2013), for warm IRDCs 25 K, 10⁵ cm⁻³ and H₂ column density 10²⁴ cm⁻² (e.g., Rathborne et al. 2008), and for HMPOs 75 K, 10⁵ cm⁻³ and column density 10²⁴ cm⁻² (e.g., Beuther et al. 2007; Rodón et al. 2008). There are known differences in H₂ column densities between interferometry and single-dish observations (Vasyunina et al. 2009). While we concentrate here on the higher end of column densities for our calculations, the calculated line fluxes can be adopted to the lower column densities by simply dividing them by 10.

For dark clouds and warm IRDCs we implement the atomic initial abundances, while for HMPOs we implemented instead the Evolution model (with initially high D/H ratios). Line intensities for a few molecules (e.g., DCO⁺) can also be estimated using the non-LTE molecular radiative transfer tool RADEX (van der Tak et al. 2007). Because RADEX and our calculations using CDMS gave different results we scaled our CDMS calculations to match RADEX values for linear molecules (e.g., DCO⁺). We look for transitions in bands 3–9 for ALMA (84–720 GHz with gaps between bands; see science.nrao.edu/facilities/alma/observing) and consider a line sensitivity limit of 1 mK for all bands. In Tables 6–8 are only the three strongest transitions listed, while a full list is available upon request from the authors.

Amongst these listed species are tracers (through D/H ratios) of initial abundances (e.g., D₂CO, HDO, NH₂D, CH₂D⁺) and temperature (D₂O, DCO⁺, DCN, D₂S) as we have listed in Table 5. A few ions are observable with ALMA, such as DCO⁺, H₂DO⁺, CH₂D⁺. We find that tracers of ionization in the cold ISM regions with high depletion, H₂D⁺ and HD₂⁺, will be hard to detect with ALMA, as found by Chapillon et al. (2011) for protoplanetary disks. We note that all isotopologues of ammonia are easily observable in the all three ISM environments. The metastable doublet lines of NH₃ are used to constrain gas temperature (e.g., Ho & Townes 1983; Maret et al. 2009). We think that using the relative abundances of the minor NH₃ isotopologues one could also discern the previous temperature history of the environment. We determine that several sulfur-bearing species, such as HDCS, D₂CS, HDS, should be observable by ALMA, and as sulfur chemistry is not yet well understood, observations of these species could serve as proxies for future studies. Finally we note that the HDO lines are expected to be observable with ALMA, but it is not included in CDMS, and we could not calculate its line intensities assuming LTE. On the other hand, water has a complex level structure, with some of the lines that are masing and many that become highly optically thick, and the escape probability non-LTE method of RADEX is not capable of modeling its line intensities reliably. For estimation of the water line fluxes one has to perform a full line radiative transfer modeling for each individual object.

For deuterated ices we only list the most abundant species in Table 9 and do not try to predict their observability. Amongst the deuterated ices we find that several polyynes (C_nH₂, with n > 4) are abundant, especially at temperatures >10 K. These species have been observed in the solar system such as in Titan's atmosphere (e.g., Teanby et al. 2009), but should also be abundant in ices in interstellar space. Several organics such as deuterated formic acid and hydroxylamine (both singly and doubly deuterated) are abundant, even among the most abundant at higher temperatures (≳ 25 K). Also both ammonia and water (singly and doubly deuterated) show high abundances.

Uncertainties for the predicted abundances are factors of ∼1.5–5 for species made of ≲ 3–4 atoms and ≳ 1.5–10 for more complex ones. Apart from the estimated abundance and D/H ratio uncertainties due to the errors in the rate coefficients, we also have error bars associated with the exact estimation of physical properties of the observed environment (uncertainties in dust emissivity, temperature, poorly known dust-to-gas mass ratio, etc.) Also, many of the published observed D/H values are based on measurements made from several sources, for which derived physical properties can differ significantly.

We now discuss how our model results compare with assorted observations. A more extensive comparison of our model results to observations of deuterated species can be found in Table 19. van Dishoeck et al. (1995) observed several deuterated species toward the protobinary source IRAS 16293, namely, DCO⁺, DCN, C₂D, HDS, HDCO, NH₂D, detected in the different regions around the protostar. The first region is the warm and dense inner core (T ≳ 80 K, n_H ∼ 10⁷ cm⁻³) found to be rich with organic molecules, the second is the circumbinary envelope with T ∼ 40 K and n_H ∼ 10⁶–10⁷ cm⁻³, where molecules such as DCO⁺ and HDCO where found, and the third is the colder, low-density outer part of the envelope with T ∼ 10–20 K and n_H ∼ 10⁴–10⁵ cm⁻³ with radicals such as CN, C₂H, and C₃H₂. We find good agreement between our calculated values and those derived for all observed species. As an example, HDCO shows a D/H ratio of 0.13 in our Evolution model compared with the observed value of 0.14, and for DCN the D/H ratio is 0.027 compared with the observed ratio of 0.013. The worst agreement we find is for C₂D and DCO⁺ with a difference of a factor of ∼5 in the D/H ratios, which is still acceptable agreement. Uncertainties for these species range between factors of 2–5, with largest uncertainties for C₂D and DCN. Considering these uncertainties, our predicted D/H values are in agreement with the IRAS 16293 observations.

Caselli et al. (2003) detected ortho-H₂D⁺ toward the prestellar core L1544, with derived abundances of 7.2 × 10⁻¹⁰ and 3.2 × 10⁻¹⁰ at the peak and off-peak positions, respectively. For our Primordial model with a core density $n_{\rm {H_2}} = 10^{6}$ cm⁻³, a temperature of 7 K, and an appropriate equilibrium o/p ratio of 3:1 taken from Walmsley et al. (2004), we find reasonable agreement with calculated abundances of 0.7–4.8 × 10⁻¹⁰ in the core (peak position) and 2.3–4.2 × 10⁻¹⁰ at the off-peak position, with a lower density, assuming a density of 10⁴–10⁵ cm⁻³. We however find a higher abundance at the off-peak position; there are several possible reasons for the discrepancy such as incorrect treatment of the ortho–para species in the network or the lack of a detailed physical model. A detailed study is not within the scope of this paper but we note that abundances are within a factor of 2–3 of observed abundances.

Stark et al. (2004) observed H₂D⁺, DCO⁺, HCO⁺, HDO and H₂O toward the protobinary source IRAS 16293 (A and B) as well as the cold prestellar object IRAS 16293 E. They measured the H₂D⁺  abundance to be 2 × 10⁻⁹ in the cold, outer envelope with $n_{{\rm H}_2}$ = 10⁴–10⁵ cm⁻³ and T < 20 K, where our Evolution model predicts a similar abundance of ∼10⁻⁹. In the inner envelope the temperature is higher, depleting the H₃⁺ isotopologues as CO returns to the gas phase, and the abundance decreases to ∼10⁻¹² cm⁻³. The temperature in the inner envelope is not well constrained, but with a central density of ∼10⁶ cm⁻³ we find the best agreement at temperatures ≈30 K, with an abundance of 1.2 × 10⁻¹², and the agreement worsens if the temperature is increased, as more CO returns to the gas phase and the overall deuterium fractionation ceases. DCO⁺ and HCO⁺ were observed with abundances 2 × 10⁻¹¹ and <1 × 10⁻⁹, respectively, which agree with our model values from 17 to 26 K with DCO⁺ abundances 8.4 × 10⁻¹¹–8.2 × 10⁻¹¹ and HCO⁺ abundances 4.7 × 10⁻¹⁰–1.1 × 10⁻⁹, leading to a D/H ratio of 0.078–0.18. HDO and H₂O were observed with abundances 3 × 10⁻¹⁰ and 3 × 10⁻⁷–4 × 10⁻⁹ respectively, which agree reasonably well with our modeled abundances of 5.9 × 10⁻¹²–6.3 × 10⁻¹¹ and 9.9 × 10⁻¹⁰–1.7 × 10⁻⁸, respectively. Lastly, HDO, DCO⁺ and HCO⁺ were also observed in the prestellar core object IRAS 16293 E. From estimated temperatures 16–25 K and densities 1.1–1.6 × 10⁶ cm⁻³, we estimate DCO⁺ abundances to be 7.0 × 10⁻¹¹–2.9 × 10⁻⁹ and HCO⁺ abundances 8.9 × 10⁻¹¹–1.10 × 10⁻⁹, leading to a D/H ratio 0.26–0.79. This agrees well with the observed abundances of 5.0 × 10⁻¹¹ and 1.0 × 10⁻¹⁰ for DCO⁺ and HCO⁺ with a D/H ratio 0.5. We find the same agreement for HDO with an observed abundance 2 × 10⁻¹⁰ comparable to our modeled abundance range of 3.0–7.5 × 10⁻¹⁰.

Coutens et al. (2012) observed multiple lines of HDO and H₂¹⁸O toward IRAS 16293A with an estimated D/H ∼0.034 in the hot corino region and D/H ∼0.005 in the outer envelope, utilizing a standard isotopic ratio of H₂¹⁸O/H₂¹⁶O = 500. In order to reproduce the observed line emission, they added an outer absorbing layer with an H₂O column density of 1.23 × 10¹³ cm⁻². Depending on the exact choice of density and temperatures, our models give for the cold envelope (n_H ∼ 10⁵ cm⁻³, T ≲ 30 K) D/H ratios of ∼0.01–0.1, for both sets of initial abundances, in agreement with observations. For hot cores (n_H ∼ 10⁸ cm⁻³, T ∼150 K) the Primordial model estimates the D/H ratio to be 0.0001–0.001, while the Evolution model predicts that the D/H ratio is ∼0.0001–0.01. Although our evolutionary model is a poor solution and does not account for the gradual warm-up of the environment, our predicted D/H ratios from the Evolution model are in agreement with those estimated by Coutens et al. (2012), albeit only with upper limits of our estimates.

The radical OD was observed for the first time outside of the solar system by Parise et al. (2012) along the line of sight toward the low-mass protostar IRAS 16293A. They also observed HDO and found a high OD/HDO ratio of ∼10–100. Parise et al. (2012) compared their observations to the modeled values of OH/H₂O and found their calculated values to be too low. The agreement was slightly better when they implemented a simple evolutionary model with increasing temperature with time, but the result was still lower than observations, with the highest modeled values reaching 5.7. Studying the chemical evolution in our Primordial model for temperatures T < 30 K and densities n_H = 10⁴–10⁶ cm⁻³, we find that the large OD/HDO ratio is mainly due to the efficient reaction OH + D $\longrightarrow$ OD + H, as originally suggested by Millar et al. (1989). Via this reaction, the OD/HDO ratios can reach values approximately one order of magnitude larger than OH/H₂O. Furthermore, toward temperatures ∼30 K the OD/HDO ratios can even be as high as two orders of magnitude. Thus, our model is in agreement with the observed OD/HDO ratios from Parise et al. (2012), without the need for a warm-up phase.

Roberts & Millar (2000a, 2000b) have investigated the chemical evolution with deuterium fractionation for temperatures 10–100 K and densities 3 × 10³–3 × 10⁸ cm⁻³ on a less resolved grid, consisting of only 100 points. They used a time-dependent chemical gas-phase model based on the UMIST'95 database (Millar et al. 1997). Their resulting network consists of ∼300 species linked by >5000 reactions, but only includes singly deuterated species and limited surface chemistry for H₂ and HD. We compared the results between our models for a number of species, including DCO⁺, HDCO, DCN, DNC, and DC₅N, looking at the distribution of D/H fractionation ratios and time-dependent abundances at 10⁵ yr, and we found good agreement between our models. We also studied the molecular abundances under conditions typical of the TMC-1 environment in our Primordial model and under these conditions we found that the quantitative agreement in the D/H ratios is better than an order of magnitude for all species, with the worst agreement for NH₂D where the ratio between the two models is 0.14. The comparison for D/H values is shown in Table 10. The intrinsic uncertainty in the abundance of DC₅N as predicted in our sensitivity analysis is very large, ∼1–1.5 orders of magnitude, and is comparable to the difference between our and Roberts' & Millar's model. We note however that our modeled D/H ratios show a better agreement with the observations of DC₅N and HC₅N in TMC-1 by MacLeod et al. (1981) than calculated values by Roberts & Millar. Roberts & Millar (2000a) expanded their study to include doubly deuterated species, allowed species to freeze onto grains and looked at a different selection of species. We compared their predictions with our results for singly and doubly deuterated isotopologues of NH₃, H₂O, H₂CO, and found reasonably good agreement for all singly deuterated species. In our model, NH₂D shows enhanced D/H ratios (∼10⁻³–10⁻¹) up to temperatures of 30–40 K, while the enhanced D/H ratios in the model of Roberts & Millar only appear up to 20–30 K. For the doubly deuterated species D₂O, NHD₂, and D₂CO, we predict D/H ratios similar to Roberts & Millar up to temperatures of ∼50 K, with values ∼10⁻³–10⁻¹, while at larger temperatures our models diverge. Our model predicts a strong decrease in the respective D/H ratios to ∼10⁻⁵, while the D/H ratios of Roberts & Millar decrease more smoothly and do not reach the same value until at ∼100 K.

In the study of Roberts et al. (2004), the chemical evolution in a sample of prestellar cores using two subsets of the Rate'99 and osu.2003 chemical networks was compared. With the networks limited to include species with six or fewer carbon atoms and no surface chemistry, Roberts et al. (2004) used the chemical models to successfully explain observations of the CO depletion, and its relevance to the D₂CO and HDCO fractionation ratios. If we compare the calculated fractionation ratios between the steady state abundances of Roberts et al. (2004; see their Table 5) and our models for a TMC-1 environment, the D/H ratios for the majority of key species such as H₂D⁺, N₂D⁺, DCO⁺, and HDO agree reasonably well. However, we found significant discrepancies for D₂O, HD₂⁺, D₃⁺, and NHD₂. The reason why our D/H ratios for doubly deuterated species differ from those of Roberts et al. (2004) appears to be twofold. First, they have not considered surface chemistry and assumed that all atomic D that freezes out is immediately returned to the gas as HD. In contrast, in our model, the accreted deuterium atoms are incorporated in surface species and do not easily return to the gas phase. Second, they use only steady-state abundances while we use time-dependent abundances at 1 Myr.

Walmsley et al. (2004) studied steady-state chemistry in a completely depleted, low-mass prestellar core, with an emphasize on explaining observations of ortho-H₂D⁺ toward L1544 by Caselli et al. (2003). While our model does not yet include nuclear-spin state chemistry, we compared the calculated abundances of the H₃⁺ isotopologues at densities of n_H = 10⁵, 10⁶, 10⁷ cm⁻³, T = 10 K assuming a cosmic ray ionization rate ζ = 3 × 10⁻¹⁷ s⁻¹ (note that at such conditions para-H₂ will be the dominant form of molecular hydrogen and thus fractionation will proceed with a high efficiency as assumed in our model). For both the Primordial and Evolution models we found good overall agreement for the D/H ratios of the H₃⁺ isotopologues and the electron abundances. The only difference occurs at high densities of ∼10⁷ cm⁻³, where the HD₂⁺ and D₃⁺ abundances are one and two orders of magnitude lower in our model, respectively. We could not find out whether this difference increases at higher densities. The likely reason for such a discrepancy is the lack of surface chemistry and the assumption of complete freeze-out in the Walmsley et al. (2004) model. Even at such high densities and 10 K, the depletion is not complete in our model, so that the H₃⁺ isotopologues can still be destroyed by ion–molecule reactions in addition to dissociative recombination with electrons or negatively charged grains. Moreover, in our model, atomic D released upon dissociative recombination can stick to a grain and be incorporated in the surface molecules. At n_H ≳ 10⁷ cm⁻³ and after 1 Myr of evolution, a substantial fraction of the gas-phase reservoir of the elemental D can be chemically "transferred" to ices, unable to directly come back to the gas phase. Consequently, it will increase surface fractionation and abundances of deuterated ices.