© ESO, 2013

1. Introduction

Starless cores, which are sites of potential low-mass star formation, are condensations of dense and cold gas. Owing to the high density and low temperature attributed to these objects, it is expected that species containing heavy elements eventually deplete onto grain surfaces because of their relatively high binding energies. Indeed, depletion of several chemical species, such as CO (e.g. Willacy et al. 1998; Caselli et al. 1999) and CS (e.g. Tafalla et al. 2002; Pon et al. 2009), has been observed toward these objects.

For a long time (since the early 1970s), gas-phase chemical models were the main means of studying chemical evolution in starless cores. However, even though it has been known for a long time that H₂ must be formed on grain surfaces (Gould & Salpeter 1963; Hollenbach & Salpeter 1971) and consequently that surface chemistry should play an important role in the chemical evolution of these objects, it is only fairly recently (during the last 10 years or so) that models including the interaction between gas phase and grain surface chemistry have really taken over as the main means of numerical studies of the chemistry in starless cores (Aikawa et al. 2003, 2005; Garrod et al. 2007).

One of the main goals of simulating chemistry in starless cores is not only to explain the observed properties (such as line emission radiation from various species) of these cores, but also to search for new tracer species for the varying physical conditions. Particularly in the very centers of starless cores where density is high, species containing heavy elements are expected to be almost totally depleted onto grain surfaces, and it has been suggested that the light deuterium-bearing species H₂D⁺ and D₂H⁺ could serve as the main tracers of these conditions (Walmsley et al. 2004; Flower et al. 2004), after observations of the ground-state rotational transition of ortho-H₂D⁺ revealed an unexpected strong line toward the prestellar core L1544 (Caselli et al. 2003). However, in the so-called complete depletion model of Walmsley et al. (2004) and Flower et al. (2004), in which no heavy elements are present in the gas phase, a description of grain surface chemistry is not included, and it has remained somewhat unclear whether the abundances of H₂D⁺ and D₂H⁺ would be affected if one were to consider surface chemistry as well.

In this paper, we address this issue by studying the chemistry in physical conditions appropriate to the centers of starless cores using a gas-grain chemical model. To this end, we have constructed new chemical reaction sets for both gas phase and grain surface chemistry that include extensive descriptions of deuterium and spin state chemistry. Our goal is not only to study how the abundances of light deuterated species are affected when grain surface chemistry is taken into account, but also to study spin state abundance ratios and compare these against the predictions of the complete depletion model.

The paper is structured as follows. In Sect. 2, we introduce the new chemical reaction sets and the assumptions made in constructing them. We also discuss our chemical and physical models in general. Section 3 presents the results of our model calculations. In Sect. 4 we discuss our results and in Sect. 5 we present our conclusions.

2. Model

We have constructed new chemical reaction sets for both gas phase and grain surface chemistry. The new gas phase reaction set is based on the publicly available Ohio State University osu_03_2008¹ (hereafter OSU) reaction set, which was expanded to include the spin states (i.e. ortho and para forms) of the light hydrogen-bearing species ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$, H₂, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ (here, this process is referred to as ortho/para separation). In addition, species with up to four atoms were deuterated using a reaction cloning process. A new surface reaction set was also produced by performing a similar ortho/para separation and deuteration of the surface reaction set of Sipilä (2012), itself based on the surface chemistry network of Semenov et al. (2010); the network includes parameters (e.g. branching ratios, activation energies if applicable) for both ionization (by either external or cosmic ray-induced photons) processes and reactions via the Langmuir-Hinshelwood mechanism (see Sect. 2.2).

We present below the physical and chemical models used in this study and then discuss the ortho/para separation and deuteration process in detail.

2.1. Physical model

In this paper, we present radial abundance profiles of various chemical species. These profiles are produced identically to the method discussed in detail in Sipilä (2012): we have chosen a modified Bonnor-Ebert sphere (e.g. Keto & Field 2005; Sipilä et al. 2011; Sipilä 2012) as the physical core model and combine chemical and radiative transfer calculations to produce self-consistently calculated radial profiles for both chemical abundances and temperatures (T_dust ≠ T_gas generally). We assume the same cooling species as in Sipilä (2012), i.e. ¹²CO (and its isotopes ¹³CO and C¹⁸O), C, O and O₂. Here we concentrate on comparing our new chemistry model with the complete depletion model, and for this purpose we choose the “model A” core of Sipilä (2012), which is a modified Bonnor-Ebert sphere of mass 0.25 M_⊙. The density and temperature profiles of the model core are plotted in Fig. 1. The core has a high average density (n_H ~ 8 × 10⁵ cm^-3) and a low average temperature (T ~ 8 K, depending slightly on core age), so this core should be appropriate for comparing our new results with the complete depletion model. The gas temperature at the core edge is different for the two core ages because at t = 10⁵ years CO, the main coolant, is not yet fully depleted at the core edge (see Sipilä 2012). We note that the gas temperature at the core edge is slightly higher at t = 10⁶ years than in Sipilä (2012; Fig. 2). This is due to a slightly lower CO abundance at the core edge brought about by the inclusion of spin-state and deuterium chemistry in the present model that slightly modifies the chemistry of undeuterated species compared to the model of Sipilä (2012).

Fig. 1 Density (green solid line) and temperature (black lines) profiles of the model core. The solid and dashed black lines correspond to the gas temperature at t = 105 and t = 106 years, respectively (see also Sipilä 2012). The dotted line represents the dust temperature.

We assume that the model core is heavily shielded against external UV radiation – we set A_V = 10 mag at the core edge. The cosmic ray ionization rate is set to ζ = 1.3 × 10^-17 s^-1 and we assume spherical grains with radius a_g = 0.1 μm. The assumed initial element abundances with respect to total atomic hydrogen density are given in Table 1. The initial H₂ ortho/para ratio is (arbitrarily) set to 1.0 × 10^-3 (we return to this issue in Sect. 4.2); the initial HD abundance is taken from Sipilä et al. (2010; hereafter S10) and the rest of the initial abundances are taken from Semenov et al. (2010).

2.2. Chemical model

The chemical code used in this study is essentially the same as the one discussed in Sipilä (2012), with some minor bug fixes. Gas phase and grain surface reactions are solved simultaneously using the rate equation method, and gas-grain interaction occurs through adsorption and desorption processes. The rate coefficient for adsorption is ${\mathit{k}}_{\mathit{i}}^{\mathrm{ads}}\mathrm{=}{\mathit{\upsilon }}_{\mathit{i}}\hspace{0.17em}\mathit{S\sigma }$, where ${\mathit{\upsilon }}_{\mathit{i}}\mathrm{=}\sqrt{\mathrm{8}{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{gas}}\mathit{/}\mathit{\pi }{\mathit{m}}_{\mathit{i}}}$ is the thermal speed of species i, S is the sticking coefficient (set to unity for all species) and $\mathit{\sigma }\mathrm{=}\mathit{\pi }{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}$ is the grain cross section. Desorption occurs mainly² through transient heating of the grains by cosmic rays, with rate coefficient k_CR = f   k_TD(70   K), where k_TD is the thermal desorption rate coefficient (Hasegawa & Herbst 1993). Photodesorption is not included in the present model; we return briefly to this issue in Sect. 4.5.

Species adsorbed onto grain surfaces are assumed to be physisorbed and react via the Langmuir-Hinshelwood mechanism. Following the formalism of Hasegawa et al. (1992), the rate coefficient for surface reactions is given by ${\mathit{k}}_{\mathit{ij}}\mathrm{=}\mathit{\alpha }\hspace{0.17em}{\mathit{\kappa }}_{\mathit{ij}}\mathrm{(}{\mathit{R}}_{\mathit{i}}^{\mathrm{diff}}\mathrm{+}{{\mathit{R}}_{\mathit{j}}^{\mathrm{diff}}}^{\mathrm{)}}\mathit{/}{\mathit{n}}_{\mathrm{d}}$(1)where α is the branching ratio in the case that the reaction has multiple product channels, κ_ij is the reaction probability, ${\mathit{R}}_{\mathit{i}}^{\mathrm{diff}}$ is the diffusion rate of species i on the grain surface, and n_d is the density of grains. In addition to being destroyed in two-body reactions, species on the surface can be photodissociated either by external UV photons (k^phot = αexp   ( − γ/A_V)) or by cosmic ray induced UV photons (k^CRphot = αζ), although the former process is negligible in the high-extinction environment assumed here.

The explicit expressions of κ_ij and ${\mathit{R}}_{\mathit{i}}^{\mathrm{diff}}$ depend on whether quantum tunneling is assumed to occur. In case of thermal diffusion ${\mathit{R}}_{\mathit{i}}^{\mathrm{diff}}\mathrm{=}\frac{{\mathit{\nu }}_{\mathit{i}}}{{\mathit{N}}_{\mathrm{s}}}\exp \hspace{0.17em}\mathrm{(}\mathrm{-}{\mathit{E}}_{\mathit{i}}^{\mathrm{diff}}\mathit{/}{\mathit{T}}_{\mathrm{dust}}\mathrm{)}\mathit{,}$(2)where ${\mathit{\nu }}_{\mathit{i}}\mathrm{=}\sqrt{\mathrm{2}{\mathit{n}}_{\mathrm{s}}{\mathit{k}}_{\mathrm{B}}{\mathit{E}}_{\mathit{i}}^{\mathrm{b}}\mathit{/}{\mathit{\pi }}^{\mathrm{2}}{\mathit{m}}_{\mathit{i}}}$; ${\mathit{N}}_{\mathrm{s}}\mathrm{=}{\mathit{n}}_{\mathrm{s}}\hspace{0.17em}\mathrm{4}\mathit{\pi }{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}$ is the number of binding sites on the grain and ${\mathit{E}}_{\mathit{i}}^{\mathrm{b}}$ and ${\mathit{E}}_{\mathit{i}}^{\mathrm{diff}}$ are the binding and diffusion energies, respectively. As in Sipilä (2012), we assume n_s = 1.5 × 10¹⁵ cm^-2 for the surface density of binding sites and ${\mathit{E}}_{\mathit{i}}^{\mathrm{diff}}\mathrm{=}\mathrm{0.77}\hspace{0.17em}{\mathit{E}}_{\mathit{i}}^{\mathrm{b}}$ for the relation between diffusion and binding energy. If a reaction is exothermic and has no activation energy, κ_ij = 1, but for exothermic reactions with activation energy, κ_ij = exp   ( − E_a/T_dust) where E_a is the activation energy of the reaction. If quantum tunneling of light species (H and D) is assumed (see also Sect. 4.3), the reaction probability and diffusion rate are – in the reactions involving these species – replaced by ${\mathit{\kappa }}_{\mathit{ij}}\mathrm{=}\exp \mathrm{[}\mathrm{-}\mathrm{2}\hspace{0.17em}\mathrm{(}\mathit{a}\mathit{/}ħ\mathrm{)}\hspace{0.17em}\mathrm{(}\mathrm{2}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{k}}_{\mathrm{B}}{\mathit{E}}_{\mathit{a}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{]}$(3)and ${\mathit{R}}_{\mathit{i}}^{\mathrm{diff}\mathit{,}\mathrm{q}}\mathrm{=}\frac{{\mathit{\nu }}_{\mathit{i}}}{{\mathit{N}}_{\mathrm{s}}}\exp \left[\mathrm{-}\mathrm{2}\mathrm{(}\mathit{a}\mathit{/}{ħ}^{\mathrm{)}}\hspace{0.17em}{\mathrm{(}\mathrm{2}\hspace{0.17em}\mathit{m}\hspace{0.17em}{\mathit{k}}_{\mathrm{B}}{{\mathit{E}}_{\mathit{i}}^{\mathrm{diff}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\right]\mathit{,}$(4)respectively (Hasegawa et al. 1992), where μ is the reduced mass of the reactants. We assume a = 1   Å for the barrier width. We note that the above assumption of a rectangular barrier is not valid for endothermic reactions. Some of the reactions with activation barriers included in the surface reaction set may be endothermic. However, because reactions with high barriers turn out not to influence our results significantly, we have chosen to omit this possible problem with the tunneling probability.

In contrast to Sipilä (2012), we assume 500 K for the binding energy of H₂, which is in the range of typical values assumed by other authors, usually somewhere between 430 K (Garrod & Pauly 2011) and 600 K (Cazaux et al. 2010). In Sipilä (2012), the binding energy of H₂ was arbitrarily decreased to 100 K to avoid the problem of producing unphysical amounts of surface H₂. (This issue has been discussed by, e.g., Garrod & Pauly 2011.) As we now consider a binding energy value of 500 K, the surface H₂ abundance is typically high in our models; however, due to the low temperature (see below) considered here, a high surface H₂ abundance should not present a problem since the surface reactions involving H₂ have high activation energies (and H₂ is not allowed to tunnel).

As in Sipilä (2012), we assume a binding energy of 1390 K for atomic oxygen (Bergeron et al. 2008; Cazaux et al. 2011). Apart from H₂ and O, the binding energies for undeuterated species are adopted from Garrod & Herbst (2006). Following the usual approach in the literature (e.g. Cazaux et al. 2010; Taquet et al. 2013b), we assume that for each deuterated species, the binding energy is equal to that of the corresponding undeuterated species. We discuss the binding energies further in Sect. 4.4.

2.3. Ortho/para separation

To add the spin states of ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$, H₂, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ into the OSU reaction set, we first extracted the reactions involving these species from the OSU database. A Python script was written to analyze each reaction and to add the spin states using a predetermined set of rules.

Most reactions are separated according to spin selection rules. As an example, consider the reaction ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{O}\stackrel{\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}$(5)with rate coefficient k. This reaction has three possible pathways, $\begin{array}{ccc}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}& \mathrm{+}& \mathrm{O}\stackrel{\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\\ {\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}& \mathrm{+}& \mathrm{O}\stackrel{\frac{\mathrm{1}}{\mathrm{2}}\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\\ {\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}& \mathrm{+}& \mathrm{O}\stackrel{\frac{\mathrm{1}}{\mathrm{2}}\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathit{,}\end{array}$(6)and the branching ratios can be calculated with, e.g., the method of Oka (2004; see also Appendix A). We note that alternative approaches to the separation exist in the literature – for example, Flower et al. (2006) assumed branching ratios of 1:2 for the separation of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}$ in reactions analogous to (6), while Pagani et al. (2009) assumed 1:1, as is done here.

The above branching ratios are not applied to all reaction types. One example of this is charge transfer reactions where we have assumed that the ortho/para forms are conserved in the reaction. Also, whenever H₂ is formed in a reaction containing only reactants other than ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$, H₂, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, such as in the reaction $\mathrm{N}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{-}\mathrm{\to }\mathrm{HN}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}$(7)we assume that H₂ is formed in its para state (see Appendix A). This is a simplifying assumption, ensuring that we do not have to follow the spin states of every species with multiple protons. There are several special cases besides reaction (7) – they are listed in Appendix A along with our assumptions.

There are approximately 1000 reactions containing ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$, H₂, or ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ in the OSU reaction file. Performing the ortho/para separation results in  ~ 1700 reactions, so that the amount of new reactions brought about by the separation process is not significant (with respect to the total size of the reaction network, which is about 11600 reactions after ortho/para separation and deuteration; see Sect. 2.4). A similar separation process was applied to the surface reaction set; in this case, the separation process adds only about 50 new reactions.

Fig. 2 Radial abundances (with respect to n(H2)) of protons, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and electrons (indicated in the figure). Panel a) is a reproduction of the middle panel of Fig. 6 in Sipilä et al. (2010) (showing abundances at t = 106 years), calculated with the version of the chemical code used in that paper (see text). The two middle panels use the same physical model as in panel a), but chemical evolution has been calculated with the present version of the chemical code (see text) with quantum tunneling turned either on (panel b)) or off (panel c)). Panel d) plots again the same physical situation but calculated with the new reaction network presented in this paper, with initial heavy element abundances set to zero.

2.4. Deuteration process

After the ortho/para separation process was carried out, we produced deuterated versions of both the gas phase and surface reaction sets. In practice, deuterons are substituted in place of protons in each reaction and combinatorial arguments are invoked (as in the case of ortho/para separation) to work out branching ratios, where applicable. For example, deuterating the reaction $\mathrm{N}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\stackrel{\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{N}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{OH}$(8)results in nine new reactions including, for example, the reactions $\begin{array}{ccc}\mathrm{N}{\mathrm{H}}^{\mathrm{+}}& \mathrm{+}& \mathrm{HDO}\stackrel{\frac{\mathrm{2}}{\mathrm{3}}\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{NH}{\mathrm{D}}^{\mathrm{+}}\mathrm{+}\mathrm{OH}\\ \mathrm{N}{\mathrm{H}}^{\mathrm{+}}& \mathrm{+}& \mathrm{HDO}\stackrel{\frac{\mathrm{1}}{\mathrm{3}}\mathit{k}}{\mathrm{-}\mathrm{\to }}\mathrm{N}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{OD}\mathit{,}\end{array}$(9)where complete scrambling, i.e. a statistical approach, is assumed. Our deuteration routine is similar to that presented in Rodgers & Millar (1996) – recently an analogous method has also been applied to the osu_2009 network by Albertsson et al. (2011). As pointed out by Rodgers & Millar (1996), the complete scrambling assumption is not appropriate for all reactions, particularly for those involving more complex molecules with OH endgroups, for example. However, as we are interested in the chemistry of relatively simple species in the present paper, we have chosen a rather relaxed approach toward this problem and assume that the complete scrambling assumption is generally valid. Of course, there are a number of special cases; for example, in charge transfer reactions no atoms are interchanged, and the deuteration routine is designed to make sure this does not occur in the deuterated version of the reaction either. Other important special cases are deuterated analogues of reactions such as ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\mathrm{-}\mathrm{\to }{\mathrm{H}}_{\mathrm{3}}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{.}$(10)We have assumed that reactions such as this proceed by ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ donating a proton and thus, for example, the reaction D₂H⁺ + H₂O cannot result in HD₂O⁺ + H₂.

Deuterated reactions involving only the lightest species (up to helium) are taken from previous works. The deuteration chemistry of the light species in the context of a complete depletion model (no heavy elements in the gas phase, and no surface chemistry except for the formation of H₂, HD and D₂) was discussed originally by Walmsley et al. (2004) and Flower et al. (2004), and subsequently by many authors. The goal of the present paper is to study how the inclusion of heavy elements and a description of gas-grain interaction including surface chemistry affects the gas-phase chemistry of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologues; we have therefore incorporated the reaction set used in S10 into the expanded OSU set so that the final reaction set is consistent with the complete depletion model as far as the chemistry of light species is concerned. In practice, we replaced the appropriate reactions in the new reaction set with their S10 counterparts and added any reactions from S10 that did not exist in the new model³. This allows for a maximally consistent comparison between the two cases.

We have similarly added to the final (deuterated and spin state-separated) surface reaction set the reactions and the associated activation energies included in the models of Cazaux et al. (2010, 2011).

As a final step, all reactions in the final deuterated gas phase and surface reactions sets were checked for elemental balance, i.e., that both sides of all reactions contain the same amount of elements. The final gas phase reaction set, including deuterium and spin state chemistry, contains about 11 600 reactions. The final surface reaction set contains about 1350 reactions.

Fig. 3 Radial abundances (with respect to nH) of protons, electrons, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs (indicated in the figure). Panel a) corresponds to the complete depletion model at t = 106 years (using the reaction set of S10), panels b) and c) correspond to the present model at t = 105 years and t = 106 years, respectively.

3. Results

3.1. Comparison of the new chemical code against the S10 code

A previous (gas-phase chemistry only) version of our chemical code was used in S10 to study deuterium chemistry in the complete depletion limit; to test the new reaction set discussed here, and the chemical code itself, we reproduced some of the results of S10 using both the old and the new versions of the chemical code (Fig. 2). Panel (a) is a reproduction of the middle panel in Fig. 6 in S10 (abundances are shown at t = 10⁶ years), calculated using the old version of the chemical code and adopting the same density and temperature profiles, the same reaction set and assuming the same values for the physical parameters as in S10. It should be noted that because only gas-phase chemistry was considered in S10, the grain-surface formation mechanisms for H₂, HD and D₂ were hard-coded into the program and no surface species as such were included. The new version of the code, however, includes surface chemistry as described above and in Sipilä (2012); to compare these two approaches, we plot in panels (b) and (c) the results of calculations otherwise identical to those of panel (a), but calculated using the new code with quantum tunneling either included (panel b) or excluded (panel c). It can be seen that without quantum tunneling, the results of S10 are not reproduced. The reason for this is the very low temperature (~6.5 K in the center) of the model core; without quantum tunneling, surface production of HD is inefficient at this temperature and eventually leads to a drop in gas phase HD abundance. This in turn increases the ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and H₂D⁺ abundances and decreases the D₂H⁺ abundance, while the ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ abundance is relatively unaltered. The reasons for these effects are discussed in, e.g., Flower et al. (2004, 2006), and Pagani et al. (2009). With quantum tunneling included, the new code gives virtually identical results compared to the old code because of the rapid conversion of surface H and D to HD. It should be noted that the creation rate of HD on grain surfaces depends critically on the assumed temperature (e.g. Cazaux & Spaans 2009); already at  ~8 K, the difference between the tunneling and no tunneling cases is very small as demonstrated by comparing results at radii  ≳ 2000 AU in panels (b) and (c) of Fig. 2, where the temperature is  ~7.5 K.

Fig. 4 Abundances (with respect to nH) of selected species (indicated in the figure) as functions of time in the innermost shell (R ~ 240 AU, nH ~ 1.9 × 106 cm-3; left panel) and the outermost shell (R ~ 3600 AU, nH ~ 2.2 × 105 cm-3; right panel) of the model core.

In panel (d) we plot the same physical situation again as in panels (a − c), using the new reaction set described above but with initial heavy element abundances set to zero (and quantum tunneling included). As would be expected since the light element chemistry originates mainly from the S10 set (see Sect. 2.4), the abundances are practically identical to those presented in panels (a) and (b). Even though the deuteration process described above introduces some reactions involving only the light species that do not exist in the S10 set, it turns out that these reactions (such as some reactions involving the H ⁻  and D ⁻  anions) are rather insignificant (at least for this set of physical conditions).

In what follows, all modeling results correspond to the new ortho/para separated and deuterated reaction set with quantum tunneling switched on, unless otherwise stated.

3.2. Comparison with the complete depletion model

In the complete depletion model, it is assumed that “heavy” species (i.e. those containing elements heavier than He) are frozen onto grain surfaces at high gas densities, allowing the deuteration chemistry of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ to proceed unhindered. The deuteration chain of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ is highly dependent on the abundance of HD (e.g. Walmsley et al. 2004), which is produced on grain surfaces and is thus dependent on the available amount of atomic H and D. In the complete depletion model, grain-surface atomic H and D are spent only in the reactions creating H₂, HD, or D₂; in essence, it is assumed that heavy species on grain surfaces are locked in unreactive species, such as water or methanol. However, in a model where heavy elements are not totally depleted, the large amount of reacting species on the surface may reduce the amount of HD produced because of the additional reaction pathways for H and D (compared to the complete depletion model), and thus affect the abundances of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs in the gas phase.

To demonstrate this, we plot in Fig. 3 the radial abundances of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs using the reaction set of either S10 (panel a) or this work (panels b and c); also plotted are the protons and electrons. Panels (a) and (c) correspond to a core age of 10⁶ years and panel (b) corresponds to t = 10⁵ years. It can be immediately seen that the abundances change dramatically depending on the adopted reaction set and that in the present model the abundances are highly time-dependent even at advanced core ages (in the complete depletion model, the abundances reach steady-state across the core at about t = 2 × 10⁵ years). In the dense center of the core, where one would expect a high degree of deuteration of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ based on the complete depletion model (e.g. Flower et al. 2004; Pagani et al. 2009), the present model predicts a completely opposite situation at t = 10⁶ years – in this case, the degree of deuteration increases toward the edge of the core, i.e. toward lower density. This can be understood through the surface chemistry; we plot in Fig. 4 the abundances of HD, HDO^∗ (in this paper, an asterisk represents a surface species) and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs as a function of time in the innermost (left panel) and outermost (right panel) shells of the model core. In the dense center of the core (n_H ~ 10⁶ cm^-3), H^∗ and D^∗ react primarily with O^∗ to produce OH^∗ and OD^∗, which react again with H^∗ and D^∗ to produce H₂O ^∗  and HDO^∗ (these processes are efficient because oxygen is initially in atomic form; see Table 1 and Sect. 4.1). This process almost completely suppresses the formation of HD^∗ through the reaction of H^∗ with D^∗ (see also Sect. 4.3). Thus, since its abundance is not efficiently regenerated on grains, HD ultimately depletes from the gas phase. HD depletion coincides with heavy element depletion – this is because the ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ abundance increases with heavy element depletion and starts to destroy HD, which is not efficiently created on grain surfaces owing to the reasons mentioned above.

The HD abundance experiences a slight increase at about t = 3.5    ×    10⁵ years. This arises because OH^∗ and OD^∗ disappear from the grain surfaces as they are processed into H₂O ^∗  and HDO^∗; this frees up some H^∗ to react with HDS^∗ and HDCO^∗ to produce HS^∗ + HD^∗ and HCO^∗ + HD^∗, respectively (see also Sect. 4.3). The HD^∗ thus formed desorbs and momentarily regenerates the HD abundance.

At the lower density of the core edge (~10⁵ cm^-3), represented by the right panel of Fig. 4, HD depletes less than at very high density (~10⁶ cm^-3), implying a higher deuteration degree of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ at an advanced core age. The slower depletion of HD at the lower density is due to the overall depletion timescale being longer in these conditions, which suppresses deuterium chemistry (as seen in the behavior of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ deuteration at early times), and thus there is somewhat less D^∗ available to form OD^∗ (and HDO ^∗  through OH^∗ + D^∗). As a consequence, the timescale for deuterium incorporation into HDO^∗ is longer than at high density. Also in this case, total HD depletion is prevented by the reactions H^∗ + HDS^∗ and H^∗ + HDCO^∗.

Fig. 5 Nuclear spin state abundance ratios of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and of H2 (multiplied by 105). Panel a) corresponds to the complete depletion model at t = 106 years, panel b) to the present model at t = 105 years and panel c) to the present model at t = 106 years.

Figures 3 and 4 demonstrate that the deuteration degree of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ is highly dependent on time. Indeed, at t = 10⁵ years, the agreement of the present model with the complete depletion model is much better than at t = 10⁶ years. With the exception of ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$, for which we obtain much lower abundances than in the complete depletion model, the agreement of the two models is rather good in the core center at t = 10⁵ years, but less so at the core edge where depletion is not as efficient. The detailed behavior of the deuteration degree depends on many factors, including the binding energies of various species, but the main result of Figs. 3 and 4 is the depletion of HD which does not occur in the complete depletion model.

From an observational standpoint it is interesting that the radial abundances of H₂D⁺ and D₂H⁺ are similar to each other in the new model at t = 10⁵ years (see Fig. 3), although their abundances are clearly lower in the new model than in the complete depletion model at the core edge. This agrees with observational results that the two species are present with comparable abundances (Parise et al. 2011; Vastel et al. 2012).

It is also observed in Fig. 3 that the ionization degree decreases when one adds heavy elements to the model. This is because in the present model, electrons are mainly removed in dissociative recombination reactions with HCO⁺ and HCNH⁺; the rate coefficients of these reactions are about a factor of 3 higher than those of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated forms, leading to a lower ionization degree in a model with heavy elements when compared to the one with zero heavy element abundances.

Fig. 6 Nuclear spin state abundance ratios of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and of H2 (multiplied by 105) as functions of time in the innermost shell (R ~ 240 AU, nH ~ 1.9 × 106 cm-3; left panel) and the outermost shell (R ~ 3600 AU, nH ~ 2.2 × 105 cm-3; right panel) of the model core.

3.3. Ortho/para ratios

In Figs. 2–4, we have plotted the abundances of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs as sums of the abundances of their respective nuclear spin states. To illustrate how the ortho/para ratios of these species differ in the present model compared with the complete depletion model, we plot in Fig. 5 the ortho/para ratios of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, H₂D⁺, D₂H⁺ and H₂ (multiplied by 10⁵) and the meta/ortho ratio of ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ as predicted by both the complete depletion model and the new model presented here. The panels depict the same core ages as in Fig. 3. Figure 5 is supplemented by Fig. 6, where the spin state ratios are plotted as a function of time in the innermost and outermost shells of the model core. As could be expected based on the evolution of the HD abundance in the present model, the spin state abundance ratios of the various species, which depend on the abundances of HD and ortho-H₂ (Walmsley et al. 2004; Flower et al. 2004; Sipilä et al. 2010), are somewhat different in the two models, depending also on core age. The ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, H₂D⁺ and D₂H⁺ spin state ratios are similar (within a factor of a few) in both models, particularly at t = 10⁵ years, but there are larger differences in the behavior of the ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{(}\mathrm{m}\mathit{/}\mathrm{o}\mathrm{)}$ ratio and the H₂(o/p) ratio depending on the model.

One can see from Figs. 5 and 6 that the o/p ratios (and the spin temperatures, T_spin) of H₂, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, and H₂D⁺ decrease as the core evolves. The effect is particularly marked for H₂, with T_spin dropping below 15 K after 10⁶ years in both the innermost and the outermost shells of the model core. Also for D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$, a slight decrease in T_spin (i.e. an increase in the o/p and m/o ratios, respectively) can be seen at late times (after 10⁶ years) when compared with early times ( < 10⁴ years). The ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{(}\mathrm{m}\mathit{/}\mathrm{o}\mathrm{)}$ ratio settles at  ~ 1 which indicates a clearly super-thermal T_spin of 42 K, but the ratio experiences a peak at around 10⁵ years with a m/o ratio of  ~ 8 (T_spin ~ 16 K). For ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, H₂D⁺ and D₂H⁺, the spin temperatures settle between 15 and 30 K at late stages of core evolution.

The ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ (m/o) ratio is much lower in the core center at t = 10⁶ years in the present model than in the complete depletion model. At the high density of the core center, the ratio first begins to increase following heavy element depletion (initiating the deuteration chemistry), but then decreases again due to HD depletion – comparing Figs. 4 and 6 reveals that a decrease of the ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ (m/o) ratio indeed coincides with HD depletion. This is because meta-${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ is most efficiently created in a spin-state-conversion reaction of ortho-${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ with HD and in the reaction between D₂H⁺(o) and HD; therefore HD depletion leads to a decrease of the ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ (m/o) ratio.

The H₂ (o/p) ratio is very different in the present model than in the complete depletion model. The effect of HD depletion can be seen here as well; H₂(o) is most efficiently created in the gas phase in the reactions $\begin{array}{ccc}& & {\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{HD}\mathrm{-}\mathrm{\to }{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}\\ & & {\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{+}\mathrm{HD}\mathrm{-}\mathrm{\to }{\mathrm{D}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}\end{array}$\arraycolsep1.75ptand ${\mathrm{D}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{HD}\mathrm{-}\mathrm{\to }{\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}$(13)so that HD depletion directly affects H₂(o) creation. HD depletion is stronger and occurs earlier at the core center than at the edge, and this accounts for the differences seen between panels (b) and (c) in Fig. 5. Note that the gas phase production of H₂(o) is now the most important factor controlling the H₂(o/p) ratio because quantum tunneling is included; see Sect. 4.3.

Interestingly, the present model predicts very similar H₂D⁺ and D₂H⁺ (o/p) ratios at t = 10⁵ years compared with the complete depletion model even at the core edge, where the H₂D⁺ and D₂H⁺ (ortho+para) abundances differ from the complete depletion model (Fig. 3). Finally we note that the ortho/para ratios predicted by our models at t = 10⁶ years are for the most part compatible with recent modeling results of Vastel et al. (2012; their complete depletion case), although our chemical network is much larger and the present model includes a time-dependent treatment of depletion.

4. Discussion

4.1. Initial heavy element abundances

In this paper, we assume that the gas is initially atomic, with the exceptions of hydrogen being in molecular form and deuterium being locked in HD (see Table 1). In other words, we do not take into account any possible chemical history of the core, i.e. the possibility that some of the heavier elements could be locked in ices at the start of the calculation. This can be particularly important for HDO^∗ which is created from OH^∗ and OD^∗ and is thus dependent on the amount of available atomic oxygen. The assumed initial HD abundance corresponds roughly to the average gas phase D/H abundance ratio, (D/H)_gas, in the Local Bubble (Wood et al. 2004). Substantially lower values of (D/H)_gas have been determined towards lines of sight with large column densities. Linsky et al. (2006) find a correlation between deuterium depletion and the depletion of metals, and suggest that D can be incorporated in carbonaceous dust grains. The assumption about the initial gas phase deuterium abundance determines the overall degree of deuteration, but is not likely to affect the relative abundances of deuterated species.

Fig. 7 Abundances of HD and HDO∗ as a function of time at a density of nH = 2 × 106 cm-3 calculated for progenitor cloud (see text) ages of 104 (solid lines), 105 (dashed lines) or 106 (dotted lines) years.

Because the present model does not include a description of core collapse (from, e.g., a diffuse cloud into the dense configuration that is the present model core), we have investigated the effect of the initial abundances on HD depletion by first calculating chemical evolution in a low-density cloud with n_H = 2 × 10³ cm^-3, visual extinction A_V = 3 mag and temperature T_gas = T_dust = 10 K (which should be appropriate for these values of density and visual extinction, Juvela & Ysard 2011). The abundances of all species corresponding to three different times (10⁴, 10⁵ and 10⁶ years) were then extracted and used as initial abundances for the dense core. Figure 7 shows the result of such calculations in a single-point chemical model with n_H = 2 × 10⁶ cm^-3, temperature T_gas = T_dust = 7.5 K and A_V = 100 mag, corresponding roughly to the interior of the model core discussed in this paper. It is observed that the older the progenitor cloud, the longer it takes for deuterium to get locked into HDO^∗. This is because after 10⁶ years of progenitor cloud evolution, about a half of the initial oxygen abundance has been converted to H₂O ^∗ , which limits somewhat the surface reactions leading to the production of HDO ^∗  in the dense core phase (the HDO ^∗ /H₂O ^∗  ratio is only about 10^-3 at the end of the progenitor cloud evolution). Notably, very little HD depletion takes place in the diffuse cloud. We also tested the effect of more extreme initial abundances by locking all oxygen initially in H₂O ^∗  (50%) and CO ^∗  (50%) and all nitrogen in $\mathrm{N}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{\ast }}$ – in these calculations, HD still depletes (by about two orders of magnitude in  ~ 10⁶ years) because surface D gets locked into multiple singly deuterated surface species (NH₂D ^∗ , HDCO ^∗  and HDO ^∗ ).

Fig. 8 Spin temperatures of gas-phase (solid line) and grain-surface (dashed line) H2 as a function of the kinetic temperature for n(H2) = 105 cm-3 and AV = 10 mag.

While the details of the chemistry change when one starts from molecular initial abundances (with respect to the model with initially atomic heavy element abundances), our main conclusion (HD depletion) remains unaffected. Figure 7 demonstrates this for n_H = 2 × 10⁶ cm^-3, but we have also performed test calculations at lower densities to confirm that HD depletion is largely unaffected in those conditions as well.

Fig. 9 As Fig. 4, but in the calculations quantum tunneling is not included.

4.2. Initial H₂ ortho/para ratio

Dense cores condense from diffuse molecular clouds where the observationally derived ortho/para column density ratios of H₂ indicate spin temperatures in the range 50–70 K (H₂(o/p) ~ 0.3 − 0.8, Crabtree et al. 2011). Crabtree et al. (2011) showed that the H₂ spin temperatures in diffuse clouds correspond to the gas kinetic temperatures because the ortho/para ratio of H₂ is thermalized through collisions with H⁺. Also in dense clouds with elevated temperatures, the spin-state conversion reactions ${\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{\to }{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}$(14)and ${\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{\to }{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}$(15)govern the H₂(o/p) ratio. As discussed by Flower et al. (2006; their Sect. 3.1), the spin temperature of H₂ follows closely the kinetic temperature down to T_kin ~ 20 K. The H₂(o/p) ratio corresponding to this temperature is 1.8 × 10^-3. Below 20 K, the endothermic para-ortho conversion reaction with H⁺ (reaction 15) starts to lose importance, and H₂(o) production is dominated by formation on grains. The consequence is that with a decreasing kinetic temperature, the H₂ spin temperature rises above the kinetic temperature. In cold, dense cores (T_kin ~ 10 K), the competition between the more favored ortho production on grains and the efficient para-ortho conversion in reactions with H⁺ and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ (and its deuterated counterparts) in the gas phase leads to an H₂ spin temperature which lies a few K above the kinetic temperature. This point is illustrated in Fig. 8 at high values of gas density (n(H₂) = 10⁵ cm^-3) and visual extinction (A_V = 10 mag); the gas and dust H₂ spin temperatures lie clearly above the kinetic temperature (dotted line) at low values of T_kin.

Based on the ideas presented above we have chosen, somewhat arbitrarily, the initial value 1.0 × 10^-3 for the H₂(o/p) ratio in our chemistry model (see Table 1). This ratio corresponds to an H₂ spin temperature of 18.7 K. Although the inital ortho/para ratio depends on the thermal history of the core, the spin temperature should in any case lie somewhere between the current kinetic temperature and 20 K, because during core condensation and cooling, H₂ has probably been thermalized down to 20 K. Therefore, we think the initial value quoted above is more realistic than the thermal value at 10 K (3.6 × 10^-7) or the statistical ratio at high temperatures (3) which is the assumed ortho/para formation ratio on grains. The choice of the initial abundances has some effect during the first  ~ 10⁵ years in our simulations; at late times they are forgotten because of the high density of the model core.

Fig. 10 As Fig. 6, but in the calculations quantum tunneling is not included.

Fig. 11 Left: ratio of the HDO∗ abundance to the initial gas phase HD abundance at different densities (indicated in the figure) as a function of time. Right: HDO ∗ /H2O ∗  ratio at different densities as a function of time.

4.3. Quantum tunneling

Whether diffusion of species via quantum tunneling occurs on grain surfaces is still a subject of some debate in modern astrochemistry. There is some experimental evidence that tunneling does not occur (e.g. Katz et al. 1999) and on the other hand some modeling results indicate that tunneling should be taken into account (e.g. Cazaux & Tielens 2004). In this paper, we have included quantum tunneling so that our analysis would be maximally consistent with the complete depletion model (see also Sect. 3.1 and Fig. 2). However, it is also prudent to take into account the possibility that tunneling does not occur and to see whether HD depletion is affected by this assumption; with this in mind, we have carried out calculations similar to those presented in the earlier sections with quantum tunneling turned off. The results of these calculations are presented in Fig. 9, which should be compared with Fig. 4. It is observed that while the detailed behavior of the abundances as a function of time changes when quantum tunneling is turned off, the main result of this work – the ultimate depletion of HD at high density – is not affected.

The details of the (deuterium) chemistry change when quantum tunneling is turned off, which can be explained as follows. When one assumes quantum tunneling, many of the surface reaction pathways for H and D that would be otherwise unavailable due to high activation energy become accessible because tunneling increases the reaction probability and the effective rate coefficient by increasing the diffusion rate (see Sect. 2.2). This effectively reduces the abundances of H^∗ and D^∗ so that it becomes increasingly difficult to form ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ and HD^∗ through the basic reactions H^∗ + H^∗ and H^∗ + D^∗, respectively. Indeed, the main pathway for producing HD^∗ at advanced core ages is H^∗ + HDCO^∗ (as discussed in Sect. 3.2) – when tunneling is included. However, without tunneling, all of the channels for which E_a ≠ 0 are effectively closed off owing to the very low temperature, and HD^∗ is mainly produced in H^∗ + D^∗ after OH^∗ and OD^∗ are processed into H₂O ^∗  and HDO^∗4. This reaction is more efficient than H^∗ + HDCO^∗ in the tunneling case and thus gives a slightly higher HD^∗ yield, and for this reason HD^∗ depletes slightly slower in the model with no tunneling (also observed in the form of a somewhat more pronounced bump in the HD abundance at high density); the decreased HD depletion also slightly increases the abundances of the deuterated isotopologs of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ at late times.

The model without tunneling also yields a larger amount of H₂(o) at late times compared to the model with tunneling, which is demonstrated in Fig. 10 (to be compared with Fig. 6). The reason for this is again the increased H^∗ abundance that leads to efficient H₂(o) creation through the H^∗ + H^∗ reaction. When tunneling is included, H₂(o) is not efficiently created on grain surfaces because 1) the H^∗ + H^∗ reaction is suppressed; and 2) other surface reactions primarily create H₂(p) owing to the adopted ortho/para separation rules (see Appendix A). An observationally important consequence is the increased H₂D⁺(o/p) ratio at late times when tunneling is turned off.

To summarize, the inclusion or exclusion of quantum tunneling modifies the detailed behavior of the deuterium chemistry and particularly the spin state ratios, but the main result of this paper (eventual HD depletion at high density) is not affected. While our assumptions regarding the production of H₂(o) naturally influence its abundance in the model, more (experimental) evidence is needed to justify either including or excluding quantum tunneling on grain surfaces so that the detailed behavior of the spin state abundances at advanced core ages can be better constrained.

4.4. Binding energies

Experimental evidence suggests that the binding energies of light species on grain surfaces depend strongly on the properties of the ice (e.g. Hornekær et al. 2005; Amiaud et al. 2006). Because in a realistic situation the properties of the surface (porosity, content of the ice) change as a function of time, one should consider time-dependent binding energies for the various species instead of the constant values considered here (and in the majority of other works). Without an appropriate model, it is difficult to quantify whether HD depletion would be affected if this were done, but to study the influence of varying binding energies we have run some test calculations using a single-point chemical model corresponding to the conditions in the center of the model core, varying the binding energies of H and H₂ and their deuterated forms between 300 K and 600 K (which should cover the typical values on either non-porous or porous amorphous solid water; Taquet et al. 2013b). In all test calculations, HD eventually depletes from the gas phase; however, the depletion is somewhat less pronounced with increasing H and D binding energy when tunneling is not included (because a higher H or D binding energy translates to lower surface reaction rates so that HDO^∗ is not formed as efficiently). If tunneling is included, the variation in binding energy does not influence HD depletion in the range of binding energies studied (unless H₂ is allowed to tunnel, see Sect. 4.5).

In this work, we assume that the binding energy of each deuterated species is the same as that of the undeuterated analogue. This is an approximation since in principle, one would expect deuterated species to be somewhat more strongly bound to grain surfaces than their undeuterated counterparts due to their slightly higher mass (Tielens 1983; Perets et al. 2005; Kristensen et al. 2011). However, we do not expect taking the binding energy difference between non-deuterated and deuterated species into account to influence the main result (HD depletion) of this paper; if we assumed higher binding energies for deuterated species, this would probably only lead to more efficient HD depletion owing to decreased desorption efficiency.

Fig. 12 As Fig. 11, but assuming H2 tunneling (in addition to H and D).

4.5. The HDO^∗ abundance

Looking at Fig. 4, it seems that the conversion of HD into HDO ^∗  depends only very loosely on density. In these high-density cases, oxygen is finally locked mainly in H₂O ^∗  while deuterium is locked in HDO ^∗ . However, at lower densities, the recycling of deuterium is slower and the late-time abundance of HDO ^∗  is consequently lower than at high density. This is illustrated in the left panel of Fig. 11, where the ratio of HDO ^∗  to the initial gas phase HD abundance at different densities is plotted as a function of time. For these calculations, we have assumed T_gas = T_dust = 10 K and A_V = 10 mag. Evidently, the relative amount of deuterium locked in HDO ^∗  at low density is small even at advanced core ages. Because HDO ^∗  is the most abundant D-bearing species on grain surfaces in our model, it follows that HD is hardly depleted at all at low density.

In the right panel of Fig. 11, we plot the HDO ^∗ /H₂O ^∗  ratio at different densities as a function of time. This plot can be compared with other studies of water ice deuteration, such as the recent study by Taquet et al. (2013b). Although our model is in many respects different from theirs, we predict similar values for the HDO ^∗ /H₂O ^∗  ratio for typical core lifetimes (10⁵ − 10⁶ years; see Figs. 8 and 10 in Taquet et al. 2013b). The agreement is very good considering the differences between our model and theirs; for instance, Taquet et al. (2013b) consider fixed values of the H₂ ortho/para ratio, whereas in our model the ratio is calculated as a function of time. Furthermore, in our model the binding energies of the various species do not vary with the H₂ content on the surface (see Sect. 2.2), which should make some difference on the progression of the surface chemistry, although this is difficult to quantify without a direct comparison of the two models (see also Sect. 4.4). Finally, our model (unlike that of Taquet et al. 2013b) does not include a treatment of photodesorption. Consequently, the desorption efficiency is somewhat lower in our model which may translate, in particular, to slightly higher HDO ^∗  abundances in the present model than in the model of Taquet et al. (2013b). It should be noted that it is unclear if HD depletion occurs in a multilayer model, such as that of Taquet et al. (2013b). Based on our results, HD depletion occurs because deuterium is mostly locked in HDO ^∗ , i.e. the late-time HDO ^∗  abundance is similar to the initial HD abundance. In the model of Taquet et al. (2013b), the H₂O ^∗  and HDO ^∗  abundances are about an order of magnitude lower than in our model in similar physical conditions, so that the entire deuterium reservoir in HD cannot be transferred solely to HDO ^∗ . However, in a multilayer model deuterium might still be efficiently locked into multiple surface species (e.g. H₂O, H₂CO, NH₃, ...), which could influence the gas-phase HD abundance – this issue should be investigated.

Aikawa et al. (2012) have recently studied deuterium chemistry in the context of a collapse model, and they predict HDO ^∗ /H₂O ^∗  ratios of  ~ 0.01. Their model does not include spin-state chemistry, and thus gives an upper limit on the strength of deuterium fractionation (which is reduced by H₂(o); Flower et al. 2006; Taquet et al. 2013b). In the model of Aikawa et al. (2012), the sticking coefficient is set to 0.5, whereas we use a value of unity, which means that the relative amount of reactive species on the surface should be larger in our model at any given time. Finally, in the model of Aikawa et al. (2012), the model core spends a relatively short time in the prestellar phase (the protostar is born in 2.5    ×    10⁵ years) – comparing the HDO ^∗ /H₂O ^∗  ratio given by their model near the core edge where n_H ~ 10⁵ cm^-3 and T ~ 10 K (Figs. 1 and 3 in Aikawa et al. 2012) with our results (Fig. 11), the agreement is again rather good.

Cazaux et al. (2011) report modeling results yielding very low HDO ^∗ /H₂O ^∗  ratios at low temperatures (and a strong temperature dependence of the ratio). In the model of Cazaux et al. (2011), H₂ is allowed to tunnel, and this allows H₂O ^∗  to be efficiently produced in the ${\mathrm{O}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ reaction despite its high activation barrier (3000 K in Cazaux et al. 2011), which in their model leads to a low HDO ^∗ /H₂O ^∗  ratio at low temperatures. We performed test calculations including the reactions of Cazaux et al. (2011) in the chemical model and allowing H₂ to tunnel as well (in addition to H and D) – the results of these calculations are plotted in Fig. 12, which assumes the same physical parameters as in Fig. 11. Evidently, HDO ^∗  is produced less efficiently at increasing densities. This is because H₂ tunneling allows the barriers associated with the reactions involving H₂ to be more easily overcome, and the ${\mathrm{O}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ reaction (in particular) becomes important. Consequently, O ^∗  and OH ^∗  are used up fast and the HDO ^∗ /H₂O ^∗  ratio levels off before the gas-phase deuterium chemistry, releasing atomic D, takes over. Thus, we find low HDO ^∗ /H₂O ^∗  ratios ( < 10^-3) even at high density, which leads to little HD depletion.

We note that while this is an interesting result, more experimental data not only on tunneling on grain surfaces, but also on the ${\mathrm{O}}^{\mathrm{\ast }}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ reaction at low temperatures, is needed to investigate this issue. Oba et al. (2012) claim to find no evidence for this reaction to proceed at low temperature. However, because of the difficulty to carry out laboratory experiments with atomic oxygen, more experiments are needed to confirm the Oba et al. result. Finally, the results of Fig. 12 are also affected by the H₂ binding energy; for example for ${\mathit{E}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{b}}\mathrm{=}\mathrm{300}$ K, surface H₂ desorbs fast, which greatly decreases the rates of the reactions involving H₂, and we get results closer to Fig. 11. Therefore, a model including time-dependent binding energies, which possibly yields less ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{\ast }}$ than the static-binding-energy model considered here, could also yield high HDO ^∗  abundances at high density even when H₂ tunneling is included. This point warrants further investigation.

4.6. Observational constraints on the HDO^∗ abundance

The depletion of heavy elements and the accompanying increase in the atomic D/H abundance ratio in the gas phase is expected to result in high degrees of deuteration in water, formaldehyde, and methanol incorporated in the icy mantles of grains (Tielens 1983). Observations suggest, however, that deuteration is clearly more pronounced in H₂CO ^∗  and CH₃OH ^∗  than in H₂O ^∗  (see e.g. Taquet et al. 2013b, and the references therein). This has been explained by the fact that water ice is formed at an early stage when CO depletion is not marked. Observations of the O-H and O-D stretching bands in the infrared have provided upper limits of  ~ 1% for the solid HDO ^∗ /H₂O ^∗  ratio towards intermediate-mass and low-mass protostars (Parise et al. 2003; Dartois et al. 2003). The O-D band at 4.1 μm is weak and broad. Gálvez et al. (2011) estimate that the detection of HDO in amorphous ice is nearly impossible if the HDO ^∗ /H₂O ^∗  abundance ratio is less than a few percent.

The gas-phase HDO/H₂O ratios in hot corinos around solar type protostars, determined recently using interferometric observations, range from  ~ 10^-3 to a few percent (Jørgensen & van Dishoeck 2010; Persson et al. 2013; Visser et al. 2013; Taquet et al. 2013a). These ratios, which are likely to reflect the deuteration of water ice before desorption, conform with relatively low HDO ^∗ /H₂O ^∗  ratios inferred from infrared absorption measurements.

The HDO ^∗ /H₂O ^∗  column density ratio towards the center of our model core is  ~ 0.03 and  ~ 0.07 at t = 10⁵ yr and t = 10⁶ yr, respectively. These values are close to the high end of the range of gas-phase HDO/H₂O abundance ratios derived in hot corinos. In principle, the HDO ^∗  column density in the model core becomes high enough to be determined through infrared absorption. Adopting the parameters used by Dartois et al. (2003), we estimate that the peak optical thickness of the 4.1 μm O-D absorption band from amorphous HDO ice (FWHM ~ 0.2   μm) is about 0.02 and 0.18 at the times t = 10⁵ yr and t = 10⁶ yr, respectively. These values are given for an offset from the center corresponding to one half of the outer radius where the visual extinction through the core is A_V ~ 23 mag. The HDO ^∗ /H₂O ^∗  abundance ratios at the two times quoted are  ~ 0.02 and  ~ 0.06. Including a low-density envelope decreases these values and increases the extinction through the cloud, but determining the ratio may be feasible towards a background star in a situation corresponding to a late evolutionary stage of the model core.

5. Conclusions

We have studied deuterium chemistry in starless cores with a gas-grain chemical model utilizing new chemical reaction sets for both gas phase and grain surface chemistry that include the spin states of the light species (${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$, H₂, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$) and deuterated forms of species with up to four atoms. It was found that at the high densities (n_H ≥ 10⁵ cm^-3) and low temperatures ( ≤ 10 K) appropriate to the centers of starless cores, HD eventually depletes from the gas phase because of surface chemistry: deuterium is efficiently locked into grain-surface HDO. Efficient HD depletion, particularly in the very centers of the cores, means that the molecular ions previously thought to be tracers of the innermost regions of starless cores, H₂D⁺ and D₂H⁺, will eventually disappear from the gas phase. It should be noted that in addition to these two ions, other deuterated species dependent on the abundance of H₂D⁺, such as DCO⁺ and N₂D⁺, show similar behavior. The timescale of HD depletion increases if chemical interaction between the gas and the dust precedes the starless core phase. In this case, part of the oxygen reservoir is locked in grain-surface H₂O in the beginning of the starless core phase, which limits surface HDO production.

While the present model predicts that HD ultimately depletes from the gas phase, our results do not imply that species such as H₂D⁺ and D₂H⁺ would be completely unusable as tracers of cold, dense gas. In contrast, at later stages of core evolution, the H₂D⁺ and D₂H⁺ abundances increase strongly towards the core edge (i.e. lower density) and the distributions become extended, which agrees with observational evidence of extended distributions of H₂D⁺ (Vastel et al. 2006, in L1544; Pagani et al. 2009, in L183) and D₂H⁺ (Parise et al. 2011, in Oph H-MM1). In this sense, a high observed H₂D⁺ and/or D₂H⁺ abundance toward a starless or prestellar core could (loosely) constrain the core age, although the detailed behavior of the deuteration chemistry is uncertain since it depends on various parameters, including the initial H₂(o/p) ratio. We note that the evaluation of the initial elemental abundances and the initial H₂(o/p) ratio would require detailed modeling of the core condensation phase. Also, more laboratory work on the O + H₂ surface reaction (and on tunneling on grain surfaces in general) is required, as this reaction can significantly affect our results.

A natural extension of this work is to perform a large-scale study of deuterium chemistry including comparisons with observations. This will be the subject of an upcoming paper. Finally, we note that the timescale of forming grain-surface HDO, critical to the gas phase abundance of HD according to our results, depends on the surface abundances of OH and OD. Because our model does not include a multilayer treatment of surface chemistry, it may overestimate the abundance of radicals available on the surface (as pointed out by Taquet et al. 2012). This could have implications for HD depletion, especially if HD depletion proceeds through conversion into HDO ^∗  as in our models. While the present model does not allow for it, this issue should be investigated.

Acknowledgments

We thank the anonymous referee for helpful comments that improved the paper. We also thank T. Albertsson, S. Cazaux, K. Crabtree, and D. Semenov for helpful discussions. O.S. and J.H. acknowledge support from the Academy of Finland grant 132291.

References

Appendix A: Details of the ortho/para separation

In this appendix we further discuss the ortho/para separation of the OSU reaction set. As mentioned in Sect. 2, there are several reaction types where branching ratios based on spin selection rules, such as in reaction (6), are not used; we list these special cases here. For each reaction type, we describe the assumptions behind the adopted branching ratios (if any) and give an example reaction.

A.1. Charge transfer reactions

We assume that in charge-transferring reactions such as ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{C}}_{\mathrm{2}}\mathrm{H}\stackrel{{\mathit{k}}_{\mathrm{A}\mathrm{1}}}{\mathrm{-}\mathrm{\to }}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$(A.1)the spin state of ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$ (or ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$) is conserved (we require that total nuclear spin in conserved in the reaction). The above reaction then separates into two reactions: ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}{\mathrm{C}}_{\mathrm{2}}\mathrm{H}\stackrel{{\mathit{k}}_{\mathrm{A}\mathrm{1}}}{\mathrm{-}\mathrm{\to }}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}$(A.2)and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}{\mathrm{C}}_{\mathrm{2}}\mathrm{H}\stackrel{{\mathit{k}}_{\mathrm{A}\mathrm{1}}}{\mathrm{-}\mathrm{\to }}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathit{,}$(A.3)with the same rate coefficient k_A1.

A.2. Reactions of species other than H$\begin{array}{}\mathrm{+}\\ 2\end{array}$ or H$\begin{array}{}\mathrm{+}\\ 3\end{array}$ where H₂ is created

In these reactions, neither ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$ nor ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ appears as a reactant; thus, reactions such as $\mathrm{N}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\stackrel{{\mathit{k}}_{\mathrm{A}\mathrm{21}}}{\mathrm{-}\mathrm{\to }}\mathrm{HN}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$(A.4)and ${\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\stackrel{{\mathit{k}}_{\mathrm{A}\mathrm{22}}}{\mathrm{-}\mathrm{\to }}{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$(A.5)fall into this category. For these reactions, we have assumed that only para H₂ is formed (and that, accordingly, the rate coefficient remains unchanged). Production of ortho-H₂ requires an exothermicity of at least  ~ 170 K (the energy difference between ortho and para states of H₂). The exothermicity of each reaction can be calculated if the enthalpies of each species involved in the reaction are known – we have however not attempted such calculations here, and proceeded on the assumption that the required exothermicity is not reached in general.

A.3. Reactions involving H$\begin{array}{}\mathrm{+}\\ 3\end{array}$ where H₂ is created

The majority of the ion-molecule reactions of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ are of the same type as reaction (5) in the main text where ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ donates a proton and spin selection rules are applied. However, there are also reactions such as ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{MgH}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{3}}}{\mathrm{-}\mathrm{\to }}\mathrm{M}{\mathrm{g}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathit{,}$(A.6)where multiple H₂ molecules are formed. For these reactions, we assume the following separation rules: $\begin{array}{ccc}& & {\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}\mathrm{MgH}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{3}}}{\mathrm{-}\mathrm{\to }}\mathrm{M}{\mathrm{g}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathit{,}\\ & & {\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}\mathrm{MgH}\stackrel{\frac{\mathrm{1}}{\mathrm{2}}{\mathit{k}}_{\mathit{A}\mathrm{3}}}{\mathrm{-}\mathrm{\to }}\mathrm{M}{\mathrm{g}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathit{,}\end{array}$\arraycolsep1.75ptand ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}\mathrm{MgH}\stackrel{\frac{\mathrm{1}}{\mathrm{2}}{\mathit{k}}_{\mathit{A}\mathrm{3}}}{\mathrm{-}\mathrm{\to }}\mathrm{M}{\mathrm{g}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathit{.}$(A.9)This logic is also applied to the (few) similar reactions that involve ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$ and a molecule with several hydrogen atoms. In principle, reaction (A.6) can also produce H₂(o) + H₂(o) which can be shown by calculating the branching ratios with the method of Oka (2004; see also Wirström et al. 2012). However, branching ratios based solely on nuclear spin statistics might not be appropriate for reactions with low exothermicies, and a simplifying assumption about the branching ratios is utilized. We note that while reactions such as (A.6) are not among the most important reactions controlling the spin states of H₂ and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ in starless cores, this important issue warrants further investigation.

A.4. Reactions involving H$\begin{array}{}\mathrm{+}\\ 3\end{array}$ or H$\begin{array}{}\mathrm{+}\\ 2\end{array}$ with no H₂ produced

These are reactions of the same type as the reaction ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}\mathrm{O}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{4}}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{.}$(A.10)In these cases, ortho and para states are destroyed with equal rates, so that the above reaction separates into ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}\mathrm{O}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{4}}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{,}$(A.11)and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}\mathrm{O}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{4}}}{\mathrm{-}\mathrm{\to }}\mathrm{O}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathit{.}$(A.12)

A.5. Other reactions

In addition to the reaction types presented above, there are also a number of reactions that do not fall into the general categories discussed above. Whenever this is the case, we apply a custom logic unique to each reaction that is, as far as possible, consistent with the above general cases. For example, the reaction ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\stackrel{{\mathit{k}}_{\mathit{A}\mathrm{5}}}{\mathrm{-}\mathrm{\to }}\mathrm{H}{\mathrm{S}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}$(A.13)separates into the reactions ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{o}\mathrm{\right)}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\begin{array}{}\mathrm{-}\mathrm{\to }\end{array}\mathrm{H}{\mathrm{S}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{(}\mathrm{o}\mathit{/}\mathrm{p}\mathrm{)}\mathit{,}$(A.14)and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}\mathrm{\left(}\mathrm{p}\mathrm{\right)}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\begin{array}{}\mathrm{-}\mathrm{\to }\end{array}\mathrm{H}{\mathrm{S}}^{\mathrm{+}}\mathrm{+}\mathrm{H}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{(}\mathrm{o}\mathit{/}\mathrm{p}\mathrm{)}\mathit{,}$(A.15)with 1:1 branching ratios for both reactions. For the former reaction, the method of Oka (2004) yields an H₂ ortho/para ratio of 2:1 (and 1:1 for the latter reaction), when H₂S can exist in either ortho or para form. But since we do not know the spin state of H₂S, we assume that it exists totally in para form and make the above simplifying assumption.

All Tables

All Figures

	Fig. 1 Density (green solid line) and temperature (black lines) profiles of the model core. The solid and dashed black lines correspond to the gas temperature at t = 105 and t = 106 years, respectively (see also Sipilä 2012). The dotted line represents the dust temperature.
In the text

Fig. 2 Radial abundances (with respect to n(H2)) of protons, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and electrons (indicated in the figure). Panel a) is a reproduction of the middle panel of Fig. 6 in Sipilä et al. (2010) (showing abundances at t = 106 years), calculated with the version of the chemical code used in that paper (see text). The two middle panels use the same physical model as in panel a), but chemical evolution has been calculated with the present version of the chemical code (see text) with quantum tunneling turned either on (panel b)) or off (panel c)). Panel d) plots again the same physical situation but calculated with the new reaction network presented in this paper, with initial heavy element abundances set to zero.

In the text

	Fig. 3 Radial abundances (with respect to nH) of protons, electrons, and ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs (indicated in the figure). Panel a) corresponds to the complete depletion model at t = 106 years (using the reaction set of S10), panels b) and c) correspond to the present model at t = 105 years and t = 106 years, respectively.
In the text

	Fig. 4 Abundances (with respect to nH) of selected species (indicated in the figure) as functions of time in the innermost shell (R ~ 240 AU, nH ~ 1.9 × 106 cm-3; left panel) and the outermost shell (R ~ 3600 AU, nH ~ 2.2 × 105 cm-3; right panel) of the model core.
In the text

	Fig. 5 Nuclear spin state abundance ratios of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and of H2 (multiplied by 105). Panel a) corresponds to the complete depletion model at t = 106 years, panel b) to the present model at t = 105 years and panel c) to the present model at t = 106 years.
In the text

	Fig. 6 Nuclear spin state abundance ratios of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and its deuterated isotopologs and of H2 (multiplied by 105) as functions of time in the innermost shell (R ~ 240 AU, nH ~ 1.9 × 106 cm-3; left panel) and the outermost shell (R ~ 3600 AU, nH ~ 2.2 × 105 cm-3; right panel) of the model core.
In the text

	Fig. 7 Abundances of HD and HDO∗ as a function of time at a density of nH = 2 × 106 cm-3 calculated for progenitor cloud (see text) ages of 104 (solid lines), 105 (dashed lines) or 106 (dotted lines) years.
In the text

	Fig. 8 Spin temperatures of gas-phase (solid line) and grain-surface (dashed line) H2 as a function of the kinetic temperature for n(H2) = 105 cm-3 and AV = 10 mag.
In the text

	Fig. 9 As Fig. 4, but in the calculations quantum tunneling is not included.
In the text

	Fig. 10 As Fig. 6, but in the calculations quantum tunneling is not included.
In the text

	Fig. 11 Left: ratio of the HDO∗ abundance to the initial gas phase HD abundance at different densities (indicated in the figure) as a function of time. Right: HDO ∗ /H2O ∗  ratio at different densities as a function of time.
In the text

	Fig. 12 As Fig. 11, but assuming H2 tunneling (in addition to H and D).
In the text