A&A 424, 905-917 (2004)
DOI: 10.1051/0004-6361:20040441

The chemistry of multiply deuterated species in cold, dense interstellar cores[*]

H. Roberts 1 - E. Herbst 2 - T. J. Millar 3

1 - Department of Physics, The Ohio State University, Columbus, OH 43210, USA
2 - Departments of Physics, Astronomy and Chemistry, The Ohio State University, Columbus, OH 43210, USA
3 - Department of Physics, UMIST, PO Box 88, Manchester M60 1QD, UK

Received 15 March 2004 / Accepted 7 June 2004

We present new models of interstellar deuterium chemistry that include more multiply deuterated species than previously considered. We first describe the updates that have been made, and compare the results given by two different underlying networks at the low temperatures and high densities characteristic of prestellar cores. We then use the models to explain recent observations of CO depletion and D2CO/H2CO ratios in prestellar cores. We find limited agreement between our models and observations for constant density scenarios, but for two of these cores we demonstrate reasonable agreement when we take the density structure of the core into account.

Key words: ISM: abundances - ISM: molecules - ISM: individual objects: $\rho$ Oph D, L1544

1 Introduction

Several new multiply deuterated species have been detected over the past two years, including ND3 (van der Tak et al. 2002; Lis et al. 2002), CHD2OH, CD3OH (Parise et al. 2002, 2004), D2S (Vastel et al. 2003) and HD2+ (Vastel et al. 2004). In addition, deuterated species in some regions have been observed to have abundances >10% of their un-deuterated analogues, a remarkable fact, given the low underlying abundance of deuterium ($\sim$10-5 with respect to H). These include D2CO (Loinard et al. 2001, 2002; Ceccarelli 2002; Ceccarelli et al. 2002), CH2DOH (Parise et al. 2002), NH2D (Saito et al. 2000; Hatchell 2003) and H2D+ (Caselli et al. 2003). The initial motivation for this study was, therefore, to add more multiply deuterated species to our gas-phase plus accretion chemical model (Roberts & Millar 2000a,b; Roberts et al. 2002) in order to explain these new observations.

Such large enhancements in the abundances of deuterium-bearing molecules can either be due to gas-phase or to grain-surface fractionation. Grain-surface reactions are undoubtedly important in producing saturated species such as methanol, water, ammonia, and hydrogen sulphide. Water ice is observed to be abundant and ubiquitous throughout the ISM (Whittet 2002) and enhanced abundances of gas-phase NH3, CH3OH, H2CO and H2S (among others) are observed in warmer regions around protostars where grain mantles have evaporated (e.g., Cesaroni et al. 1994; Hatchell et al. 1998; Buckle & Fuller 2002).

The most recent observational and theoretical evidence suggests that the deuterium fractionation in star-forming regions is set by gas-phase and grain-surface reactions during the cold, dense pre-protostellar phase. Species are believed to form on grain surfaces via H atom addition to CO, N, O and S. Deuterium fractionation on grains, therefore, comes from the relative amounts of atomic D and H which are accreting from the gas. The observations of deuterated methanol (Parise et al. 2002, 2004) and D2S (Vastel et al. 2003) require that the gas-phase atomic D/H ratio at the time the molecules formed was $\geq$0.1, a result which is now predicted by our most recent models of deuterium chemistry in prestellar cores (Roberts et al. 2003).

It is the freeze-out of molecules from the gas onto grain surfaces, particularly CO, that causes these enhancements in deuterium fractionation in the gas-phase. This effect was predicted by Brown & Millar (1989), and has been confirmed by the detection of large [DCO+]/[HCO+] ratios in L1544 in clumps in which CO is significantly depleted (Caselli et al. 1999), and by the correlation between CO depletion and D2CO fractionation in several prestellar cores (Bacmann et al. 2002, 2003). High NH2D fractionation has been observed close to the dust peaks of low-mass protostellar sources (Hatchell 2003) under conditions that support a gas-phase origin for the molecules.

Additionally, CO depletions have now been observed in many different sources (e.g., Caselli et al. 1999; Bergin et al. 2002; Redman et al. 2002; Bacmann et al. 2002; Tafalla et al. 2002; Pagani et al. 2002), so that heavy CO depletion may be more common than has been previously assumed.

2 The models

There are two main chemical reaction schemes used by the community: "RATE99'', the latest release of the UMIST database for astrochemistry (Le Teuff et al. 2000), and the so-called "New Standard Model'' (NSM), from the Ohio State University. The most recent release of the latter, known as osu2003 (Smith et al. 2004), is available at www.ohio-state.edu/~eric/research.html.  Here we utilize the previous version (Terzieva & Herbst 1998). Aikawa et al. (2003) have shown that the RATE99 and NSM schemes can give different predictions for interstellar chemistry, particularly for sulphur-bearing species. Their models of the collapsing prestellar core, L1544, include mono-deuterated analogues for each reaction set, but the fact that grain-surface chemistry and gravitational collapse are also included makes it difficult to pin down exactly how the different rate coefficients affect the chemistry. In this paper we consider the gas-phase deuterium fractionation predicted by each model for both purely gas-phase and gas-phase plus accretion models, and discuss the key rates that cause some results to differ.

We have taken a sub-set of each scheme, which includes only the elements H, He, C, N, O, Si, Fe, and S, and species with 6 or fewer carbon atoms. In the discussion that follows, we use "RATE99'' to refer to the reaction set based on the UMIST database, and "NSM'' to refer to the set based on the NSM. We also use the notation [X] or "fractional abundance'' to refer to the abundance of species X with respect to H2, and "fractionation of XD'', "molecular D/H ratio'' or "[XD]/[XH]'' to refer to the abundance of a deuterated species with respect to its totally hydrogenated analogue.

 New branching fractions for the dissociative recombination of N2H+; viz.,

{\rm ~N_2H^+ + e^- }\longrightarrow \left\{ \begin{array}{l}...
...% \\ {\rm ~NH +~H} \hspace{7mm} 65 \pm5\%, \end{array}\right .
\end{displaymath} (1)

recently measured using CRYRING (Geppert et al. 2004), have been incorporated. Previous models have always contained the assumption that 100% of the recombinations of N2H+ with electrons produce N2 and H, since molecular nitrogen has a triple bond. In gas-phase chemical models, this new branching fraction does not have a significant effect, since destruction of N2H+ by CO is as important as destruction by electrons, and because the NH can reform N2. When the interaction with grains is taken into account, however, CO is depleted by freezing onto grains and reaction (1) becomes more important. In Roberts et al. (2003) the assumption that N2 desorbs efficiently from grains at 10 K (Tafalla et al. 2002; Caselli et al. 2002; Bergin & Langer 1997) keeps nitrogenic species abundant in the gas phase for $\sim$106 yr, whereas CO, CS, O, etc. freeze out after a few times 104 yr. Adopting the updated branching fraction in reaction (1) causes the nitrogenic species to deplete from the gas phase on a much shorter time-scale. This is discussed further in Sect. 2.3.1.

2.1 Deuterating the reaction sets

Instructions for producing a deuterated reaction set have been published in previous papers (e.g., Millar et al. 1989; Rodgers & Millar 1996). We now have 340 gas-phase species: 125 of which are singly deuterated, 30 doubly deuterated, 17 triply deuterated, 4 with 4 deuterium atoms, and 1 species (CD3OD2+) with 5.  Deuterating RATE99 results in a set of 10 389 reactions, while deuterating the OSU scheme gives 10 799 reactions.

We have not included all the possible deuterated isotopomers of every species because the time taken to run a pseudo-time dependent model depends critically on the number of species. Instead, we considered the expected importance of each deuterated analogue in the overall chemistry. For example, for CH5+ we included only its singly, doubly, and triply deuterated analogues, disregarding CHD4+ and CD5+. The justification is that CH5+ is primarily destroyed by dissociative recombination with electrons, and so is important in forming methane (and in our model we include only the CH3D and CH2D2 isotopomers) and CH3 (which, in turn, forms formaldehyde). Detailed investigation shows that inclusion of the added isotopomers of CH5+ has little effect on the isotopomers of methyl or those of methane currently included in the models.

Fractionation reactions for mono-deuterated species are listed in Roberts & Millar (2000a). In these reactions, D and H atoms exchange places. The fractionation reactions that we have adopted for multiply deuterated species are listed in Table 1. For the analogues of H3+, the reaction rates are taken from those measured at 80 K (Giles et al. 1992).  No temperature dependence is assumed for the forward reactions, and a simplified purely exponential temperature dependence is assumed for the backwards (endothermic processes). Even in thermal systems, this temperature dependence is likely to be accurate only at the lowest temperatures (e.g. 10 K), where the rotational partition functions of reactants and products are close to unity. Non-thermal effects for the backwards reactions, similar to that considered by Gerlich et al. (2002) for H2D+ + H2, have not yet been derived. That there is some temperature dependence for the forward reactions is probable, because the rate coefficients determined by Gerlich et al. (2002) and Gerlich & Schlemmer (2002) at 10 K for reactions involving H3+, H2D+, CH3+, C2H2+, and C2H2D+ with HD tend to be lower than the Giles et al. values by up to a factor of five. It seemed better to us to maintain consistency rather than to lower the values of a subset of the fractionation reactions. Our view might change if theoretical justification for the low rates could be provided, in which case we would be able to lower all of them suitably. We have therefore not adopted the Gerlich et al. values for the forward rate coefficients of selected fractionation reactions. The effect of using the lower forward values in a model with mono-deuterated species has been studied by Roberts et al. (2002) and Gerlich et al. (2002) and found to reduce fractionation. Utilizing experimental rates for the backward reactions at very low temperatures without extreme care (see Gerlich et al. 2002) to remove excitation effects is unwise at this stage.

Table 1: Fractionation reactions added to the list given in Roberts & Millar (2000a). The forward rate is $k_{\rm f}$; the reverse rate, $k_{\rm r}$, is given by $k_{\rm f}$ e $^{-\gamma /{\rm T}}$.

In addition to fractionation reactions, there are many reactions in our model where deuterated species undergo non-exchange processes. Where the reactions of the deuterated species have more product channels than their un-deuterated analogues, we have assumed statistical branching ratios (Millar et al. 1989), except for the cases where experimental measurements indicate otherwise; e.g., for the recombination of the deuterated isotopomers of H2D+ (see Table 2) and for HD2O+ (Jensen et al. 2000) with electrons, and for the special cases described below. In the case of dissociative recombination reactions, this is a cautious view, since one could justify an alternative approach in which the isotope effects detected in the few experiments hold in a larger number of cases, which would have the effect of enhancing deuterium fractionation. This question has been studied by Le Petit & Roueff (2003).

We have included all possible deuterated isotopomers of methanol for the first time, adopting the assumption that the structure of each is preserved during reactions, an assumption meaning that deuterium atoms that are part of a methyl group on a reactant will also be in the methyl group on the product (see e.g., Millar et al. 1989; Rodgers & Millar 1996). Suggested by Charnley et al. (1997) to explain the observations of CH3OD and CH2DOH in the Orion Compact Ridge, this hypothesis is supported by recent quantum chemical calculations (Osamura et al. 2004) and by recent observations of deuterated methanol in the protostellar source IRAS 16293-2422 (Parise et al. 2002; 2004). Those isotopomers of methanol without a D on the oxygen atom (CH2DOH, CHD2OH and CD3OH) are highly fractionated, whereas CH3OD is not. The simplest explanation for this difference is that the CH3OD is preferentially destroyed in warm gas by protonation to CH3OHD+ followed by dissociative recombination to CH3OH, a result which arises naturally from the assumption that the methyl group is unaffected by such reactions.

Table 2: Rate coefficients ( $k_{d{\rm r}} = \alpha (T/300)^{\beta }$) for the dissociative recombination of H3+ and deuterated analogues.

Table 3: Rate coefficients ( $k~=~\alpha\left({T}/300\right)^{\beta}$) for radiative association reactions.

The rate coefficients and branching ratios that we have adopted for the dissociative recombination of H3+ and analogues with electrons are listed in Table 2. The rate coefficients for radiative association reactions of deuterated species, which were derived in Millar et al. (1989), have been updated and extended for multiply deuterated species, and are listed in Table 3. The assumption is that the deuterated species will associate more rapidly, due to their higher mass and lower vibrational frequencies, with enhancements based on the three-body measurements of Smith et al. (1982a). Selected radiative association rate coefficients in the table are set to zero; in these cases there is an exothermic exchange reaction that competes with the radiative association channel, drastically reducing it. Although the actual rate coefficients are likely to be small and non-zero, setting them to zero is a dramatic means of eliminating the process from our models.

A special word is in order for the association reactions involving methyl ion and water. Luca et al. (2002) have recently shown that the reaction rate coefficient for the undeuterated analogues at a temperature of $50 \pm 30$ K is significantly lower than utilized in our model. We are currently working with this group to obtain a new suitable value that includes temperature dependence. The results will not affect this paper since the gas-phase synthesis of methanol, heretofore not very efficient anyway, will become noticeably less so. It is likely that any methanol seen in the gas phase was produced on the surfaces of dust particles by hydrogenation of CO (Stantcheva & Herbst 2003) and desorbed into the gas.

2.2 The gas-phase models

\end{figure} Figure 1: Fractional abundance of selected species over time; T=10 K, n(H2)=104 cm-3. The filled symbols are the RATE99 model, the hollow symbols are the NSM model.
Open with DEXTER

Figure 1 compares the evolution of selected fractional abundances for non-deuterium bearing species from the RATE99 and NSM models. We use the "standard'' interstellar dark cloud parameters: $T_{\rm kin}=10$ K, n(H2)=104 cm-3, and assume no freeze-out of species onto grains. The gas-phase material is initially atomic (except for H2 and HD), with the elemental abundances listed in Table 4. We adopt a cosmic ray-ionisation rate of  $1.3 \times 10^{-17}$ s-1, and a visual extinction of 10 mag.

One of the major differences between the models is that the rates for ion - polar-molecule reactions in the NSM model have a temperature dependence, typically a factor of $\sqrt{300/{T}}$ (Herbst & Leung 1986), whereas the RATE99 reactions do not. In practice, this factor causes the destruction of polar neutral species (e.g., H2O, H2CO, CH3OH, NH3 and H2S), by molecular and atomic ions (e.g., H3+ and HCO+, C+ and H+) to be $\sim$5 times faster in the NSM model at 10 K. Hence, the steady-state abundances of these neutral species are lower than in the RATE99 model. The present state of measurements concerning ion-polar molecule reactions at low temperatures indicates that an inverse temperature dependence pertains to most but not all of the few reactions studied. The NSM/osu2003 networks build on this small data base and theoretical justification whereas the RATE99 network is based on a more cautious viewpoint; both choices are in our view defensible.

Table 4: The initial abundances, relative to total H, used in the models.

The time dependence of the abundances of the species may also vary significantly between the models. In the NSM model, C radiatively associates with C3 to form C4 which either reacts with O to form the un-reactive CO, or re-forms C3. The C + C3 association is not included in the RATE99model; however, so more C atoms react with H2 to form CH2 and with H3+ and HCO+ to form CH+. In the RATE99 model, therefore, carbon remains in a reactive form in the gas-phase for a longer time, leading to the peak in H2CO and CH3OH abundances at  $2\times 10^5$ yr (top two panels of Fig. 1).

Disagreement between the models also occurs for some of the sulphur-bearing species, especially at early times (bottom panel of Fig. 1). The first sulphur-bearing molecule to become abundant is CS. The destruction of CS by He+, which re-forms atomic sulphur, is $\sim$50 times faster in the NSM model at 10 K. This means that significant amounts of SO (and SO2) have formed by 104 yr in the NSM model, while in the RATE99 model it takes a few times 105 yr for them to become abundant.

Table 5: A comparison of steady-state molecular D/H ratios from gas-phase models. T=10 K, n(H2)=104 cm-3.

It has long been appreciated that one of the advantages of looking at molecular D/H ratios, rather than fractional abundances, is that they will be less sensitive to any variations in the underlying rate coefficients than will the absolute molecular abundances. This is borne out by Table 5, which lists deuterium fractionation for selected species at 10 K. The agreement between the models is generally good, well within a factor of 2 for all species except DCS+ and D2S. The branching ratios for the dissociative recombination of H3S+ and analogues differ between the two reaction schemes. In the RATE99 model 50% of the recombination reactions produce H2S and its analogues, while in the NSM model, the fraction is 75%. Thus, the conversion of H2S to H2DS+ to HDS, and so on to D2S, is more efficient in the NSM model. In the NSM model, HCS+ and DCS+ are primarily destroyed by O atoms (to form OCS+), since the rate for the reaction of HCS+ with O is 100 times larger in the NSM than in the RATE99 model. In the RATE99 model HCS+ and DCS+ are primarily destroyed by electrons, forming CS, which then reacts with H3+ and analogues to reform HCS+ and DCS+. As the bottom panel of Fig. 1 shows, the steady state abundance of CS in the NSM model is significantly lower than in the RATE99 model. Thus, in the NSM model significant amounts of HCS+ and DCS+ are made via the reaction of C+ with SO, forming CS+, followed by the reaction of CS+ with H2 and HD. As the HD/H2 ratio is much lower than the H2D+/H3+ ratio, this leads to a lower DCS+/HCS+ ratio in the NSM model.

Initially, all the deuterium is assumed to be contained in HD, with an abundance relative to H2 of  $3.2 \times 10^{-5}$ (Wright et al. 1999; Hébrard et al. 2002). Table 5 shows that in a gas-phase model, at steady state, this interstellar reservoir of deuterium is not significantly depleted. The abundances of the mono-deuterated species, relative to their un-deuterated analogues do become enhanced by factors of 1000-4000 over this underlying D/H ratio, but the fractionation in the multiply deuterated species is limited to 260 times the underlying D/H ratio (and is typically much lower). These results have been presented in previous models (Millar et al. 1989; Roberts & Millar 2000a,b), and suggest that the large abundances of multiply deuterated species that have now been observed cannot be produced by a purely gas-phase chemistry.

2.3 The accretion model

Models of deuterium chemistry based on the UMIST RATE99 and RATE95 networks, which include freeze-out of species onto grains, have previously been presented for H2 densities of 104 cm-3 (Roberts & Millar 2000a,b) and $3\times10^6$ cm-3 (Roberts et al. 2003). In the following discussion we present model results for an H2 density of 106 cm-3; all other parameters are the same as for the gas-phase models, unless otherwise stated. The models have been updated from those described in previous papers in a manner discussed below.

2.3.1 Updates from previous work

In Roberts & Millar (2000a,b), it was assumed that all species heavier than H2D+ eventually freeze onto the grains and remain there. We have now extended the model to include surface analogues of the neutral species, so we can consider desorption of surface species. Currently we have included only thermal desorption, using the method described in Hasegawa & Herbst (1993). We have taken desorption energies from them and from Aikawa et al. (1997). We now assume that molecular ions landing on the grains break up into neutral fragments, with the same branching fractions for the products as we use for dissociative recombination with electrons in the gas-phase, and that atomic ions landing on the grains are neutralised. The neutral products can return to the gas via evaporation. Consideration of ion-grain collisions is important at high densities, where the fractional electron abundance is low.

In our previous comparison with L1544 (Roberts et al. 2003), we made the assumption that N2 does not freeze out. This assumption is supported by observations of prestellar cores that show systematic differentiation between CO and CS, which disappear at the centres of the cores, and N2H+ and NH3, which remain in the gas (Tafalla et al. 2002; Bergin et al. 2002). In our new models we, therefore, follow Bergin & Langer (1997) and Aikawa et al. (1997, 2003) in lowering the binding energy of N2 on the surface so that it remains in the gas-phase even at 10 K.

In Roberts et al. (2003), we found that the inclusion of HD2+ and D3+ can give rise to enormous fractionation in the accretion models, so that eventually the reservoir of deuterium in the gas-phase, HD, can be used up. We can avoid the total loss of deuterium from the gas by assuming that all the atomic deuterium which freezes out is immediately returned to the gas as HD. This assumption is not physical, however, since grain-surface chemistry will also be occurring, and a proportion of the accreting H and D will react with heavier surface species. For example, CO, O, N and S will be converted to methanol, water, ammonia, hydrogen sulphide and their deuterated analogues.

A simple calculation using diffusion rates of various surface species, with a gas-phase atomic [D]/[H] ratio of 0.3 suggests that $\sim$90% of the accreting D will form HD, $\sim$5% will form D2 and $\sim$5% will react with heavier species. Rate equations may not be the best way of modelling surface chemistry though, especially when the average abundances of species per grain are $\la$1 (Tielens & Hagen 1982; Charnley et al. 1997; Caselli et al. 2002a), and so we have also used the model of Stantcheva & Herbst (2003) to see what proportion of atomic D is predicted to form HD and D2 in a stochastic calculation. Their model indicates that 53% of the accreting D will form HD, and 18% will form D2 (T. Stantcheva; private communication), so this branching is what we have incorporated into our model. We note, however, that the other gas-phase D/H ratios are not particularly sensitive to this assumption. In any case, at late times, the loss of deuterium from the gas-phase is much less catastrophic than in Roberts et al. (2003) because dissociative recombination of molecular ions on the grain surfaces also produces significant amounts of H, D, H2, HD and D2 which quickly evaporate.

As discussed above, reaction (1) causes the nitrogenic species to deplete from the gas phase on a much shorter timescale than we previously found, since, instead of 100% of the recombinations of N2H+ with electrons forming N2, which does not freeze out, 65% now form NH, which does. In order to keep the nitrogenic species in the gas-phase for significantly longer times than the C-, O-, and S-bearing species, we assume that the N2H+ and N2D+ that land on grains are immediately returned to the gas as N2 and H or D, rather than using the branching ratios from reaction (1). This assumption is justifiable because there is a finite binding energy of electrons to grains that must be subtracted from the exothermicity of dissociative recombination reactions.

The abundances of NH and ND are also affected by the adoption of these new branching ratios, since the dissociative recombination of N2H+ and N2D+ now becomes an important formation channel for them. In previous accretion models at 10 K and n(H2)=104 cm-3, at the point where the N2D+ fractionation peaks, ND/NH is 0.08, and [ND] is  $3\times10^{-11}$. With the branching ratios from above, however, ND/NH reaches 0.34 and [ND $] \sim 10^{-8}$. At a density of 106 cm-3, the effect of the new branching ratios is to increase the peak ND/NH ratio from 1.2 to 7.3. Since ND has a sub-mm spectrum (Saito & Goto 1993), it is observable.

2.3.2 A comparison of the models

\end{figure} Figure 2: Fractionation of selected species vs. time from the accretion models; T=10 K, n(H2)=106 cm-3. Left: RATE99 model; right: NSM model.
Open with DEXTER

Figure 2 shows molecular D/H ratios vs time from the accretion models, while Fig. 3 shows fractional abundances, using $T_{\rm kin}=10$ K and n(H2)=106 cm-3. Table 6 lists the fractionation of selected species at two different times for each model: at 104 yr, when all the gaseous species are still relatively abundant ([CO $] \sim 6{-}9\times 10^{-7}$), and at  $6\times10^4$ yr, when the C- and O-bearing species have mostly frozen out.

The first thing to notice is that the fractionation of D3+ in the NSM model remains similar to that of H2D+ and HD2+ until $\sim$ $3\times 10^4$ yr have passed, while in the RATE99 model the D3+ fractionation has become very large ($\sim$10) by 104 yr. The reason for this discrepancy is twofold. Firstly, the rates for destruction of H3+ and analogues by many neutral species, including NH3, H2O, OH, CH and NH, are faster in the NSM model than the RATE99 model at 10 K (due to the $\sqrt{300/{T}}$ factor, see above). This means that the abundances of the neutral species have to be more depleted in the NSM model than the RATE99 model to affect the fractionation in the primary ions by the same amount. The second reason is that N2 stays in the gas for a longer time in the NSM model, suppressing the fractionation (this is illustrated in Fig. 3). In the NSM model N2H+ reacts fairly rapidly with NH, NH2 and NH3, reforming N2. In the RATE99 model, though, the rates for those reactions are slower, so that NH, NH2 and NH3 are more likely to collide with the grains and be lost from the gas, while the N2H+ is more likely to recombine with electrons to form NH and N, which then also freeze out.

This difference means that the carbon, oxygen, and sulphur bearing species (e.g. DCO+, DCS+, HDCO and D2CO) do not become as highly fractionated in the NSM model, since the CO and O freeze out before the H3+ fractionation gets very large. It also causes a problem in forming large abundances of deuterated formaldehyde and methanol on the grain surfaces. Since grain-surface chemistry occurs quickly, relative to the gas-phase, the fractionation in species formed on the surface is expected to reflect the atomic [D]/[H] ratio in the gas at the time the species are accreting. For the time period over which the CO and O freeze out in the NSM model (1000-40 000 yr), the [D]/[H] ratio reaches 0.1, which is lower than the observations of methanol in IRAS 16293 require (Parise et al. 2002, 2004). The RATE99 model fares better here: over the time period when CO and O are freezing out, [D]/[H] is 0.2 after 104 yr and rises to 0.76.

On the lower panel of Fig. 2, when the fractionation of DCO+, N2D+, DCN and NH2D drop abruptly to "zero'', this actually means that one or both of the species' abundances has dropped to <10-20, so that taking a ratio of them would not be meaningful. Water, OH and OD, on the other hand, never completely freeze out and so have a high steady-state fractionation. This is because a small proportion of the atomic oxygen which desorbs from the grains reacts with D3+ and analogues, rather than immediately freezing out again, to form OD+ and OH+. Successive reactions with H2, and deuteration by D3+, give H3O+ and analogues, which recombine with electrons to produce OH, OD, and water. Thermal desorption of O is very slow at 10 K, however, so only a very small amount of water and OH remains in the gas-phase (<10-16 with respect to H2).

We have included the species HDS and D2S in Table 6 even though hydrogen sulphide is not produced efficiently in the gas-phase, and is considered a tracer of surface chemistry (e.g., Minh et al. 1991; Hatchell et al. 1999). Vastel et al. (2003) detected D2S in B1 and at the DCO+ peak of NGC 1333 IRAS 4A, with an abundance $\sim$10% of HDS. Based on observations towards L1689N (Vastel et al. 2003; Wakelam et al. 2004), the H2S and HDS abundances appear to be $\sim$a few times 10-10, whereas the fractional abundances of H2S and HDS in our models never reach 10-11. Vastel et al. concluded that the HDS and D2S they observed probably formed on grain surfaces during the cold, dense, pre-protostellar phase from S, H, and D atoms accreting, with a D/H ratio of $\sim$0.2-0.3.

\end{figure} Figure 3: Fractional abundance of selected species vs. time from the accretion models; T=10 K, n(H2)=106 cm-3. Left: RATE99 model; rightNSM model.
Open with DEXTER

Table 6: A comparison of fractionation from the accretion models. T=10 K, n(H2)=106 cm-3.

3 Comparisons with observations of prestellar cores

One of the strongest pieces of evidence linking freeze-out with high deuteration is the study by Bacmann et al. (2003) of CO depletion and D2CO fractionation in pre-stellar cores. These authors found D2CO/H2CO ratios ranging between 0.01 and 0.1, with the largest D2CO ratios from the sources which have the highest CO depletion.

In Fig. 4 we have updated their Fig. 2, so that instead of comparing with the model predictions derived from Roberts & Millar (2000b) at 104 cm-3, we compare D2CO/H2CO ratios vs. CO depletion from the RATE99 model at different densities. The best fit to the observations occurs for densities of  $5 \times 10^4{-}10^5$ cm-3, although depletion factors somewhat higher than those observed are still required to explain most of the D2CO/H2CO ratios. For n(H2)=105 cm-3, D2CO/H2CO reaches 0.04 after  $7\times 10^4$ yr of accretion, whereas, at  $5\times 10^4$ cm-3, it takes almost  $2\times 10^5$ yr.

For n(H2)=104 cm-3 the D2CO fractionation never exceeds 2%, since the overall fractionation is being limited by the electron abundance. As the density increases, the electron abundance is lowered (Brown & Millar 1989; Roberts & Millar 2000a), so that the peak deuterium fractionation is increased. For n(H $_2)=5\times10^4$ cm-3 the peak D2CO/H2CO ratio is 0.09, while for n(H2)=105 and 106 cm-3, it is 0.2 and 0.5, respectively. For the higher density models, this is not evident from Fig. 4 since at higher densities the freeze-out of CO onto grains happens more quickly, and so the peak D2CO/H2CO ratios occur when CO is depleted by factors >60. Figure 5, however, shows how the CO abundance and D2CO fractionation evolve over the whole time period when CO is freezing out at densities of 105 and 106 cm-3. The peak D2CO/H2CO ratio is two and a half times higher at 106 cm-3 than at 105 cm-3, but, since molecules freeze out more quickly at higher densities, the CO abundance is significantly lower. When the D2CO fractionation has reached $\sim$7%, for example, CO is depleted from its canonical value (10-4 with respect to H2) by a factor of 120 at a density of 105 cm-3, but by a factor of 440 when the density is 106 cm-3.

As discussed in Sect. 2.3.2, the NSM model requires even heavier depletions to enhance the molecular D/H ratios, so the agreement would be worse. We note, however, that the abundances of D2CO from Bacmann et al. (2003), a few times 10 -12-10-11, are reproduced (and even exceeded) by both models at times >1000 yr, until the formaldehyde freezes out. Thus, the problem is not that we do not make enough D2CO, rather that we produce too much H2CO.

The transfer of fractionation between the primary ions and formaldehyde is not particularly efficient. Deuterons are transferred from H2D+, HD2+ and D3+ to formaldehyde, but the branching ratios which determine the products when H2DCO+ and analogues then recombine with electrons do not favour the formation of deuterated formaldehyde. For H3CO+ recombining with electrons only one-third of the reactions in the models produce H2CO and H; the other two-thirds break up into CO + H + H2 and HCO + 2H.

Even so, our models may still be overestimating the proportion of HDCO and D2CO formed in this way, since the lowest energy isomer of protonated formaldehyde is actually "H2COH+'', rather than "H3CO+'' (Ohishi et al. 1996). If deuterons transferred from H2D+, HD2+ and D3+ to H2CO and HDCO, form H2COD+ and HDCOD+, respectively, it may not be possible for these ions to recombine to HDCO and D2CO.

Of course, our single density-single temperature model is not really representative of these cores. Bacmann et al. (2002) calculate an overall density and depletion factor for each prestellar core, but, in reality, we expect these quantities to increase towards the core centres. For example, from their C17O measurements, Bacmann et al. find CO to be depleted by factors of 14 in L1544 and 4.5 in L1698B. However, Lee et al. (2003) have observed C18O in L1544, finding that within 0.075 pc ($\sim$15 000 AU) of the centre, the depletion factor is at least 25, while Redman et al. (2002) have also observed C17O towards L1689B and estimate that over 90% of the CO has frozen out within the central 5000 AU. It could be, therefore, that the high D2CO/H2CO ratios observed by Bacmann et al. originate in the inner regions of the prestellar cores, where the CO is heavily depleted, while the bulk of the CO emission is coming from a lower density, less depleted, envelope.

\end{figure} Figure 4: Model predictions for D2CO fractionation as a function of CO depletion at 10 K, from the RATE99 model at various densities, compared with the observations of Bacmann et al. (2002, 2003) towards prestellar cores.
Open with DEXTER

\end{figure} Figure 5: Model predictions for the rise in D2CO fractionation at 10 K (right-hand axis) as CO freezes out (left-hand axis). The solid lines are for n(H2)=106 cm-3; the dashed lines are for n(H2)=105 cm-3.
Open with DEXTER

\end{figure} Figure 6: Physical parameters of prestellar cores as a function of radius. Top: L1544 (Tafalla et al. 2002); bottom: Oph D (Motte et al. 1998; Bergin et al. 2002). The dark lines represent the best fit to the observations, while the lighter `stairs' show the model fits we have adopted.
Open with DEXTER

Figure 6 shows density profiles for L1544 and for core D in $\rho$ Ophiuchi (hereafter Oph D) derived from the observations of Tafalla et al. (2002) and Motte et al. (1998), respectively. L1544 is probably the best studied of the prestellar cores, while Oph D is the source in which Bacmann et al. (2003) observed the highest D2CO fractionation ($10 \pm4$%). Using this information, we consider a set of models, which cover the range of physical conditions seen in these sources. The temperatures and densities of the model fits we have made are also shown in Fig. 6, so L1544 is divided into 5 "slabs'', Oph D into 6. A model is run for each slab, then any molecular abundance from each is multiplied by the path length and they are all summed to get a prediction for the column density through the centre of each core. The outer envelope of L1544 (T=17 K; n(H $_2 )=2\times10^4$ cm-3) is assumed to extend from 11 500 AU to 26 500 AU, while the envelope of Oph D (T=17 K; n(H $_2)=3\times10^4$ cm-3) extends from 14 200 to 24 200 AU. The cores are considered to be sufficiently embedded that photo-reactions are unimportant.

3.1 L1544

\end{figure} Figure 7: Fractional abundances vs. radius for selected molecules in L1544, using the density and temperature profiles from Fig. 6.
Open with DEXTER

We have calculated column densities for a core with a structure similar to L1544, as shown in Fig. 6. We initially ran a model using the RATE99 model under the physical conditions of the outer envelope for 105 yr. We then used the abundances it returned as inputs into the higher density models. Figure 7 shows the abundances we predict as a function of radius after a few times 104 yr, the time at which we find the best fit to the observations. In the less dense envelope, adsorption of species onto grains is slower, and at 17-20 K many species can desorb from the grains so CO remains abundant ( $8\times 10^{-5}$ with respect to H2). It is depleted more than 100 times in the inner region between 3000 and 7000 AU where the D2CO abundance peaks, and almost completely frozen out within $\sim$2000 AU of the centre. Atomic oxygen is also completely depleted in the core centre, whereas N2 and NH3 remain relatively abundant. The ion D3+ becomes more abundant than its analogues at $\sim$6000 AU, and is actually the dominant ion at the core centre.

\end{figure} Figure 8: Molecular D/H ratios vs. radius for selected molecules in L1544, using the density and temperature profiles from Fig. 6.
Open with DEXTER

Figure 8 shows how the molecular D/H ratios vary across the core. As we would expect, the fractionation is low in the warmer, less dense, envelope, and increases as we move towards the core centre. The atomic D/H ratio exceeds 0.1 at a radius of $\sim$7000 AU, so it is within this region that we would expect the most highly fractionated ices to be forming.

Table 7 lists fractional abundances in each slab and the predicted column densities through the centre of the core. We include data for those species which have been observed in L1544, as well as the deuterated species which the models predict to have the highest column densities.

The total column density of D2CO is $3.1\times 10^{12}$ cm-2, $\sim$3 times higher than was observed by Bacmann et al. (2003), while N(CO)/N(H $_2) = 4.6\times10^{-6}$, depleted by a factor of $\sim$20 from the canonical value, compared to a depletion factor of 14 derived by Bacmann et al. (2002). These numbers, along with our overall D2CO/H2CO ratio of 0.03, are now within the uncertainties of the observations, and the agreement is significantly better than the single density models in Fig. 4, where a D2CO/H2CO ratio of 0.03 required a CO depletion of at least 35. Also, the transitions of H2CO and D2CO observed by Bacmann et al. have critical densities $\geq$105 cm-3 (Aikawa et al. in prep.), so their observations may not be sensitive to the formaldehyde in the outer layers of the cloud. Figure 6 shows that the core density exceeds 105 cm-3 within $\sim$7000 AU of the centre. In the outer layers of the cloud, the abundance of D2CO is low, so neglecting them does not significantly affect the total column density. Almost 50% of the total H2CO, however, is contained in the layers of the cloud with n(H2)<105 cm-3. If this component is not detected by the observations, then our prediction for the D2CO/H2CO ratio (based on the central 3 slabs of the model) is $\sim$0.05.

The fractional abundance and column density of H2D+ we predict in L1544 are  $5\times 10^{-10}$ and  $5\times 10^{13}$ cm-2, respectively, reproducing the values observed by Caselli et al. (2003), and we predict that HD2+ will have a similar abundance. Our prediction for N(N2H+) is $\sim$3 times lower than the observation, while our predicted N2D+/ N2H+ and DCO+/HCO+ ratios are 1.2 and 0.2, respectively, $\sim$5 times higher than the observation towards the dust peak of L1544 (Caselli et al. 2002b).

Tafalla et al. (2002) observed several starless cores, including L1544, finding that CO and CS are heavily depleted at the core centres, that N2H+ has a constant fractional abundance across the core, and that the NH3 abundance increases towards the core centre. Figure 7 shows that, in our models, N2H+ and NH3 have fairly constant fractional abundances at radii >6000 AU, but NH3 is 34 times less at 3000 AU, while N2H+ is reduced 70 times. Our estimates for [N2H+] and [NH3] at the core centre are an order of magnitude lower than the observed values of  $7.5 \times 10^{-11}$ and  $4\times10^{-9}$, respectively (Tafalla et al. 2002).

Aikawa et al. (2001, 2003) have also presented models of prestellar cores similar to L1544, in which the dynamics of the cloud collapse are modelled to give the evolution of the density profile over time. The chemistry is based on the "NSM'' set of reaction rates, rather than the RATE99set adopted here. Aikawa et al. (2003) also include grain-surface reactions and find that the formation of molecular nitrogen on the grain surfaces, followed by its desorption into the gas, enhances the abundance of N2H+ and NH3 at the core centre. In our model of L1544, when n(H $_2)\sim 10^6$ cm-3 the fractional abundance of atomic nitrogen on the grain surfaces reaches 10-5 after only 500 yr, so we would also probably get better agreement with the observations of N2H+ and NH3 towards the centre of L1544 by including grain surface chemistry in the models.

Although deuterated ammonia has not been observed in L1544, it is worth noting how our model results compare to the detections which have been made. The highest NH2D/NH3, NHD2/NH3 and ND3/NH3 ratios which have been observed to data are 0.33 (in HH211; Hatchell 2003), 0.03 (in 16293 E; Loinard et al. 2001) and 0.001 (in NGC 1333 and Barnard 1; van der Tak et al. 2002; Lis et al. 2002), respectively, which are remarkable, given the low underlying abundance of deuterium, but for NHD2 and ND3 our predicted ratios in Table 7 exceed these by 6-20 times. Ammonia fractionation becomes very large due to its high proton affinity, the highest of any common interstellar molecule, which means that NH3 is efficiently converted to NH3D+ and then to NH2D, which forms NH2D2+, and so on until the ammonia is completely deuterated. Desorption of some N2 from the grain surfaces, discussed above, would act to reduce the overall fractionation to some degree (since it destroys H3+). Also, the sources mentioned above appear to be protostellar in nature, with associated outflows which could remove material from the grains and/or higher temperatures which would suppress deuterium fractionation.

3.2 Oph D

\end{figure} Figure 9: Fractional abundances ( top) and molecular D/H ratios ( bottom) vs. radius for selected molecules in Oph D, using the density and temperature profiles from Fig. 6.
Open with DEXTER

The central density of Oph D is not as high as L1544, so, as Fig. 9 shows, freeze-out at the core centre after a few times 104 yr is not so advanced. Again, CO and O are abundant in the outer envelope, but are depleted by factors of 200-300 within $\sim$12 000 AU. The N2, NH3 and N2H+ abundances remain approximately constant outside of 6000 AU, though within this region their abundances are starting to decrease.

Table 8 lists fractional abundances and molecular D/H ratios. The overall column density of D2CO is  $4.2 \times10^{12}$ cm-2, slightly higher than that observed by Bacmann et al. (2003), while the average D2CO/H2CO ratio is 0.05, half of the observed value, but within the uncertainty of the observations. The overall CO column density we calculate is  $4\times
10^{17}$ cm-2, making the fractional abundance of CO  $3.6 \times 10^{-6}$, depleted by a factor $\sim$28. Again, this agreement is better than for the single density models, where a depletion factor of at least 150 is required to give a D2CO/H2CO ratio $\geq$0.05. Considering only the regions of the cloud where n(H2) exceeds 105 cm-3, the column density of H2CO is reduced to  $4\times10^{13}$ cm-2 and the calculated D2CO/H2CO ratio becomes 0.086.

4 Conclusions

From our comparison of the reaction sets, we conclude that the RATE99 reaction set better reproduces the high observed fractionation of D2CO. For grain-surface chemistry to produce high abundances of deuterated methanol and hydrogen sulphide, a high atomic D/H ratio at the time that CO and S are freezing out is required; again, the RATE99 reaction set better reproduces the needed value. The NSM scheme is more successful, however, at keeping N2D+ and NH3 in the gas after the CO has frozen out. Both reaction networks could probably be improved in this respect by including the formation of molecular nitrogen on grain surfaces.

Our single density-single temperature models are too simplistic to explain the fractionation in prestellar cores fully. Introducing a density gradient across the core allows us to reconcile, to some degree, the observed high D2CO/H2CO ratios and relatively low CO depletions observed by Bacmann et al. (2002, 2003), but do not explain all the observations made towards L1544. In particular, the predicted fractionations of DCO+ and N2D+ are too high. Injecting N2, which has formed on the grains, and/or including other desorption mechanisms (such as cosmic-ray desorption) which return other material to the gas at 10 K would suppress the DCO+/HCO+ and N2D+/N2H+ ratios, but this might also reduce the D2CO/H2CO and atomic D/H ratios too much.

So far we have considered the effects of changing the temperature and density, since these, along with the underlying elemental abundances, are the parameters that most strongly affect deuterium fractionation. A change in the initial abundances of heavy elements by a common factor would only affect the time scale of the accretion models although lowering these abundances in the gas-phase models would indeed enhance fractionation. Roberts & Millar (2000a) investigated changing the C/O ratio and found that it didn't have a noticeable effect on the deuterium fractionation. The rate for cosmic ray ionisation will also have some effect on the deuterium fractionation, since an enhanced cosmic ray ionisation rate increases the electron abundance, and so reduces fractionation. In our density profile for Oph D, we have assumed that the dense core merges into the ambient cloud, and that the extinction is high enough that photo-reactions driven by external photons are unimportant in the outer layers of the core. Motte et al. (1998) find that the dust emission does merge with the ambient cloud in the east-west direction, but that there is a sharp edge to the southwest. If the visual extinction were lower in the outer layers of the cloud the abundances of more complex species such as H2CO might be lowered, possibly helping the overall D2CO fractionation, but in order to model this we would need to use a PDR code.

To improve these models further, we need more experimental data. For accurate desorption rates, the binding energies for species onto grains are required. In the gas-phase, we need to confirm the rates for the fractionation reactions at 10 K, and determine the branching ratios for dissociative recombination of more deuterated ions. In evaluating the two major gas-phase networks, it would be helpful to determine the temperature dependence of more ion - polar neutral reactions.

The authors thank Ted Bergin for providing data on the structure of Oph D, and for useful discussions. H.R. and E.H. acknowledge the support of the National Science Foundation (US) for their research in astrochemistry. Astrophysics at UMIST is supported by PPARC.



Online Material

Table 7: Fractional abundances in each model slab (with respect to H2), as well as total column densities and overall deuterium fractionation predicted for L1544.

Table 8: Fractional abundances in each model slab (with respect to H2), as well as total column densities and overall deuterium fractionation for Oph D.

Copyright ESO 2004