A&A 492, 703-718 (2008)
DOI: 10.1051/0004-6361:20079009

Survey of ortho-H2D+ (11,0-11,1) in dense cloud cores

P. Caselli1,2 - C. Vastel3,4 - C. Ceccarelli5 - F. F. S. van der Tak6,7 - A. Crapsi8 - A. Bacmann5,9

1 - School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK
2 - INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
3 - Université de Toulouse, UPS, CESR, 9 avenue du colonel Roche, BP 44346, 31028 Toulouse Cedex 04, France
4 - CNRS, UMR 5187, 31028 Toulouse, France
5 - Laboratoire d'Astrophysique, Observatoire de Grenoble, BP 53, 38041 Grenoble Cedex 9, France
6 - Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
7 - National Institute for Space Research (SRON), Postbus 800, 9700 AV Groningen, The Netherlands
8 - Leiden Observatory, PO Box 9513, 2300 RA Leiden, The Netherlands
9 - Université Bordeaux 1, CNRS, OASU, UMR 5804, 33270 Floirac, France

Received 6 November 2007 / Accepted 16 September 2008

Aims. We present a survey of the ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) line toward a sample of 10 starless cores and 6 protostellar cores, carried out at the Caltech Submillimeter Observatory. The high diagnostic power of this line is revealed for the study of the chemistry, and the evolutionary and dynamical status of low-mass dense cores.
Methods. The derived ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities (N(ortho- ${\rm H}_{2}{\rm D}^{+}$)) are compared with predictions from simple chemical models of centrally concentrated cloud cores.
Results. The line is detected in 7 starless cores and in 4 protostellar cores. N(ortho- ${\rm H}_{2}{\rm D}^{+}$) ranges between 2 and 40 $\times $ 1012 cm-2 in starless cores and between 2 and 9 $\times $ 1012 cm-2 in protostellar cores. The brightest lines are detected toward the densest and most centrally concentrated starless cores, where the CO depletion factor and the deuterium fractionation are also largest. The large scatter observed in plots of N(ortho- ${\rm H}_{2}{\rm D}^{+}$) vs. the observed deuterium fractionation and vs. the CO depletion factor is likely to be due to variations in the ortho-to-para (o/p) ratio of ${\rm H}_{2}{\rm D}^{+}$ from >0.5 for $T_{\rm kin}$ < 10 K gas in pre-stellar cores to $\simeq$0.03 (consistent with $T_{\rm kin}$ $\simeq$ 15 K for protostellar cores). The two Ophiuchus cores in our sample also require a relatively low o/p ratio ($\simeq$0.3). Other parameters, including the cosmic-ray ionization rate, the CO depletion factor (or, more in general, the depletion factor of neutral species), the volume density, the fraction of dust grains and PAHs also largely affect the ortho- ${\rm H}_{2}{\rm D}^{+}$ abundance. In particular, gas temperatures above 15 K, low CO depletion factors and large abundance of negatively charged small dust grains or PAHs drastically reduce the deuterium fractionations to values inconsistent with those observed toward pre-stellar and protostellar cores. The most deuterated and ${\rm H}_{2}{\rm D}^{+}$-rich objects (L 429, L 1544, L 694-2 and L 183) are reproduced by chemical models of centrally concentrated (central densties $\simeq$106 cm-3) cores with chemical ages between 104 and 106 yr. Upper limits of the para- ${\rm H}_{3}{\rm O}^{+}$(11--21+) and para- ${\rm D}_{2}{\rm H}^{+}$ (11,0-10,1) lines are also given. The upper limit to the para- ${\rm H}_{3}{\rm O}^{+}$ fractional abundance is $\simeq$10-8 and we find an upper limit to the para- ${\rm D}_{2}{\rm H}^{+}$/ortho- ${\rm H}_{2}{\rm D}^{+}$ column density ratio equal to 1, consistent with chemical model predictions of high density (2 $\times $ 106 cm-3) and low temperature ( $T_{\rm kin}$ < 10 K) clouds.
Conclusions. Our results point out the need for better determinations of temperature and density profiles in dense cores as well as for observations of para- ${\rm H}_{2}{\rm D}^{+}$.

Key words: astrochemistry - stars: formation - ISM: clouds - ISM: molecules - radio lines: ISM - submillimeter

1 Introduction

In the past decade, astrochemistry has become more and more crucial in understanding the structure and evolution of star forming regions. There are no doubts that stars like our Sun form in gas and dust condensations within molecular clouds and that the process of star formation can only be understood by means of detailed observations of the dust (a good probe of the most abundant and elusive molecule, ${\rm H}_2$) and molecular lines (unique tools to study kinematics and chemical composition).

Millimeter and submillimeter continuum dust emission observations (see Ward-Thompson et al. 2007, for a detailed review of this topic) are very good probes of the density structure of dense cores, although uncertainties are still present regarding the dust opacity and temperature, both likely to change within centrally concentrated objects (but such variations have so far been hard to quantify observationally, see e.g. Bianchi et al. 2003; Pagani et al. 2004,2003). Stellar counts in the near-infrared provide an alternative way of measuring the dust (and ${\rm H}_2$) column and cloud structure (Lada et al. 1994), independent of any variation in dust properties, but they cannot probe regions with extinctions above about 40-50 mag (Alves et al. 1998), i.e. the central zone of very dense cores, such as L 1544, where $A_{\rm V} \simeq$ 100 mag within 11 ${\hbox{$^{\prime\prime}$ }}$ (Ward-Thompson et al. 1999). It appears that many starless cores can be approximated as Bonnor-Ebert spheres (Bonnor 1956; Ebert 1955), with values of the central densities ranging from about 105  ${{\rm cm}}^{-3}$ (as in the case of B 68; Alves et al. 2001) to 106  ${{\rm cm}}^{-3}$ (for e.g., L 1544, L 183, and L 694-2; Ward-Thompson et al. 1999; Pagani et al. 2003; Harvey et al. 2003). At the lower end of the central density range, dense cores appear to be isothermal, with gas temperatures close to 10 K (Tafalla et al. 2004; Galli et al. 2002), whereas higher density cores have clear evidence of temperature drops in the central few thousand AU, with dust temperatures approaching about 7 K (Crapsi et al. 2007; Pagani et al. 2004; Evans et al. 2001; Schnee & Goodman 2005; Pagani et al. 2003,2007; see also Bergin & Tafalla 2007, for a comprehensive review on starless cores).

Keto & Field (2005) have proposed that the ``shallower'' cores (such as B 68) are in approximate equilibrium and will not evolve to form protostars, whereas the centrally concentrated ones (such as L 1544) are unstable cores that are proceeding toward gravitational collapse and the formation of protostars. Indeed, this is in agreement with the findings of Lada et al. (2002), who claim that B 68 is oscillating around an equilibrium state, and those of Caselli et al. (2002a) and van der Tak et al. (2005), who studied the kinematic structure of L 1544 and found that it is consistent with contraction in the core nucleus (or central contraction).

To trace the gas properties, ${\rm NH}_{3}$ and ${\rm N}_2{\rm H}^+$ have been extensively used for several years (Benson & Myers 1989) and they seem to trace quite similar conditions, having comparable morphologies and line widths (Tafalla et al. 2002; Benson et al. 1998; Tafalla et al. 2004; Caselli et al. 2002c), despite the two orders of magnitude difference in the critical densities of the most frequently observed transitions ( ${\rm NH}_{3}$(1, 1) and ${\rm N}_2{\rm H}^+$(1-0), see Pagani et al. 2007). These two species are among the few that are left in the gas phase at volume densities above $\sim$105  ${{\rm cm}}^{-3}$. CO, CS and, in general, all the carbon bearing species so far observed (with the exception of CN; Hily-Blant et al. 2008) are heavily affected by freeze-out in the central parts of dense starless cores (Tafalla et al. 2006,2002; Pagani et al. 2005; Tafalla et al. 2004; Bacmann et al. 2003; Bergin et al. 2002; Crapsi et al. 2005; Caselli et al. 2002b,1999; Kuiper et al. 1996; Crapsi et al. 2004; Willacy et al. 1998; Pagani et al. 2007; Bergin et al. 2001; Bacmann et al. 2002). The freeze-out of neutral species (in particular of CO, Roberts & Millar 2000b; Dalgarno & Lepp 1984; Roberts & Millar 2000a) boosts the deuterium fractionation in species such as ${\rm N}_2{\rm H}^+$, ${\rm NH}_{3}$, H2CO, and ${\rm HCO}^+$(see e.g. Butner et al. 1995; Bacmann et al. 2003; Gerin et al. 2006; Crapsi et al. 2005; Tiné et al. 2000; Lis et al. 2006; Caselli et al. 2002b). In star forming cores, in particular in the direction of Class 0 sources, the first protostellar stage (e.g. André et al. 2000), immediately after the pre-stellar phase, the deuterium fractionation is found to be very large (Lis et al. 2002b,a; Loinard et al. 2002; Marcelino et al. 2005; Parise et al. 2006; Vastel et al. 2003; Ceccarelli et al. 1998; Parise et al. 2002; Crapsi et al. 2004; Emprechtinger et al. 2008; Parise et al. 2004; van der Tak et al. 2002; Ceccarelli et al. 2007). This is thought to be the signature of a recent event in which the star forming cloud core experienced the low temperature and high density conditions typical of the most centrally concentrated starless cores. Some deuterated molecules (e.g. deuterated ammonia) are formed in the gas phase (and stored on dust surfaces) whereas others (such as deuterated methanol and formaldehyde) are likely formed onto dust surfaces and then partially released to the gas phase upon formation (Garrod et al. 2007,2006). The interaction with the newly born protostar can (i) heat dust grains, leading to mantle evaporation (as in Hot Cores and Corinos; Bottinelli et al. 2004a,2007; Turner 1990; Cazaux et al. 2003; Bottinelli et al. 2004b); and/or (ii) ``erode'' dust mantles via sputtering in shocks produced by the associated energetic outflows (e.g. Lis et al. 2002a).

Questions that are still open include: (1)  ${\rm N}_2{\rm H}^+$, ${\rm NH}_{3}$ and their deuterated forms do trace the inner portions of centrally concentrated cores on the verge of star formation? Although Bergin et al. (2002) and Pagani et al. (2005, 2007) found evidence of depletion in the center of B 68 and L 183, respectively, there are no signs of freeze-out for NH3 in L 1544 (Crapsi et al. 2007) and for both nitrogen bearing species in L 1517B and L 1498 (Tafalla et al. 2004). At densities above $\sim$106  ${{\rm cm}}^{-3}$, the freeze-out time scale is quite short ($\sim$1000 yr) and all heavy species are expected to condense onto grain mantles. Moreover, recent laboratory measurements clearly show that ${\rm N}_2$, the parent species of both ${\rm NH}_{3}$ and ${\rm N}_2{\rm H}^+$ should freeze-out at the same rate as CO (having similar binding energies and sticking probabilities; Öberg et al. 2005; Bisschop et al. 2006). (2) For how long is the high degree of deuterium fractionation observed in starless cores maintained after the formation of a protostellar object?

The detection of strong ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) emission in the direction of L 1544 (Caselli et al. 2003), and the conclusion that ${\rm H}_{3}^{+}$, with its deuterated counterparts, is one of the most abundant molecular ions in core centers, have opened a new way to study the chemical evolution (Aikawa et al. 2005; Flower et al. 2005; Roberts et al. 2003; Flower et al. 2004,2006b; Walmsley et al. 2004; Roberts et al. 2004; Flower et al. 2006a) and the kinematics (van der Tak et al. 2005) of the central few thousand AU of starless cores. Thus, ${\rm H}_{2}{\rm D}^{+}$ is an important tool to understand the chemical and physical properties of the material out of which protoplanetary disks and ultimately planetary systems form.

In the gas phase at temperatures below $\sim$20 K, the deuterium fractionation is mostly regulated by the proton-deuteron exchange reaction:

$\displaystyle %
{\rm H}_{3}^{+} + {\rm HD} \rightarrow {\rm H}_{2}{\rm D}^{+} + {\rm H}_2 + \Delta {E},$     (1)

where $\Delta E$ (=230 K, Millar et al. 1989) prevents the reverse reaction from being fast in cold regions, unless a significant fraction of ${\rm H}_2$ is in the ortho form (Walmsley et al. 2004; Gerlich et al. 2002). Flower et al. (2006a) showed that values of the ortho-to-para (o/p) ${\rm H}_2$ ratio much higher than 0.03 are inconsistent with the observed high levels of deuteration of the gas (see their Fig. 6). The above reaction, together with the freeze-out of neutral species (which boosts the production rate of ${\rm H}_{2}{\rm D}^{+}$ compared to ${\rm H}_{3}^{+}$; e.g. Aikawa et al. 2001), allows the ${\rm H}_{2}{\rm D}^{+}$/ ${\rm H}_{3}^{+}$ ratio to overcome the cosmic D/H ratio by several orders of magnitude. Not only has ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) been found to be $\simeq$1 K in L 1544, but also ${\rm D}_{2}{\rm H}^{+}$ has been detected toward another starless condensation (16293E, in Ophiuchus; Vastel et al. 2004). The ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) line has been mapped in L 1544 (Vastel et al. 2006a), finding that the ${\rm H}_{2}{\rm D}^{+}$ emitting region has a radius of about 5000 AU, comparable to the size of the ${\rm N}_2{\rm D}^+$(2-1) map made in the same region by Caselli et al. (2002a).

In this paper we present new ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) observations, carried out with the Caltech Submillimeter Observatory (CSO) antenna, in the direction of 10 starless cores and 6 cores associated with very young protostellar objects. As it will be shown, the line has been detected in 7 of the 10 starless cores and in 4 out of 6 star forming cores. In Sect. 2 the observational details are given. Results and ortho- ${\rm H}_{2}{\rm D}^{+}$ spectra are shown in Sect. 3, together with a brief discussion of the upper limits of the para- ${\rm D}_{2}{\rm H}^{+}$( 11,0-10,1) lines (para- ${\rm H}_{3}{\rm O}^{+}$(11--21+) upper limits can be found in the on-line Appendix). ${\rm H}_{2}{\rm D}^{+}$ column densities are derived in Sect. 4. A chemical discussion, aimed at interpreting the observations, is given in Sect. 5 and conclusions are in Sect. 6.

2 Observations

2.1 Technical details

Observations of the ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) line ($\nu_0$ = 372.421385(10) GHz; Amano & Hirao 2005) were carried out with the Caltech Submillimeter Observatory (CSO) on Mauna Kea (Hawaii), between October 2002 and April 2005. The spectra were taken in wobbler switching mode, with a chop throw of 300 $^{\prime\prime}$. The backend used was an acousto-optical spectrometer (AOS) with 50 MHz bandwidth. The velocity resolution, as measured from a frequency comb scan, is 0.1  ${\rm km~s}^{-1}$. The beam efficiency ( $\eta_{\rm b}$) at $\nu$ = 372 GHz was measured on Saturn, Mars and Jupiter and is listed in Table 2. Measurements for extended sources were made for only a few sources (L 1544, L 183, NGC 1333-DCO+, B 1, NGC 2264G-VLA2) and were found to be $\sim$70%, compared to 60% for planets measurements. ${\rm H}_{2}{\rm D}^{+}$ is likely to be extended (e.g. L 1544: Vastel et al. 2006a; L 183: Vastel et al. 2006b). Consequently, we used the extended source beam efficiency whenever available. However the difference is not highly significant. At 372 GHz, the CSO 10.4-m antenna has a half power beamwidth of about 22 ${\hbox{$^{\prime\prime}$ }}$.

Similar setups were used for the para- ${\rm H}_{3}{\rm O}^{+}$(11--21+) line at 307.1924100 GHz (JPL catalogue to be found at http://spec.jpl.nasa.gov/), which was observed at CSO in October 2002 and June 2003. Table A.1 presents the beam efficiencies that were used for H3O+ data. Pointing was measured every two hours and found to be better than 3 $^{\prime\prime}$.

Observations of the para-D2H+ (11,0-10,1) line ($\nu$ = 691.660483(20) GHz; Amano & Hirao 2005) were carried out at CSO between April 2003 and April 2005 under very good weather conditions (225 GHz zenith opacity less than 0.065). We used the 50 MHz AOS with a spectral resolution better than 0.04 km s-1. The observations were performed using the wobbler with a chop throw between 150 ${\hbox{$^{\prime\prime}$ }}$ and 180 ${\hbox{$^{\prime\prime}$ }}$ according to its stability. The beam efficiency was carefully and regularly checked on Mars, Venus, Saturn and Jupiter, and found to be $\sim$40%. For more extended sources, the beam efficiency was measured to be $\sim$60%. This value has been adopted for 16293E, the only source where D2H+ has been detected, assuming that in this case the emission covers an area not significantly smaller than the beam. For the other sources, we use $\eta_{\rm b}$ = 40% (Table 4) and consider a factor of 1.5 uncertainty in the column density upper limit value due to the unknown source size. Pointing was monitored every 1.5 h and found to be better than 3 ${\hbox{$^{\prime\prime}$ }}$. At 692 GHz, the CSO 10.4-m antenna has a half power beamwidth of about 11 ${\hbox{$^{\prime\prime}$ }}$.

2.2 Source selection

The source list is given in Table 1, which reports the coordinates, the Local Standard of Rest velocity ( $V_{\rm LSR}$) at which we centered our spectra, and the distance to the source.

Table 1: Source sample.

The selection criteria for starless cores is similar to those described in Crapsi et al. (2005), where sources with bright continuum and ${\rm N}_2{\rm H}^+$ emission have been selected to include chemically evolved cores, where CO is significantly frozen onto dust grains and where ${\rm H}_{2}{\rm D}^{+}$ is thus expected to be more abundant. The sample consists of ``shallow'' cores, with central densities of $\sim$105  ${{\rm cm}}^{-3}$ (L 1498, TMC-2, L 1517B, B 68) and more centrally concentrated ones, with central densities of $\sim$106  ${{\rm cm}}^{-3}$ (TMC-1C, L 1544, L 183, Oph D, L 429 and L 694-2).

The star forming regions observed have been selected as being representative of the early phases of protostellar evolution, so that any detection of ${\rm H}_{2}{\rm D}^{+}$ can be compared with starless cores to see if any evolutionary trends appear. Among the protostellar cores we selected:

NGC 1333 ${\rm DCO}^+$, close to IRAS 4A, in Perseus, where large abundances of deuterated species have been observed (ND3, van der Tak et al. 2002; NH2D, Hatchell 2003; D2S, Vastel et al. 2003; ND2H, Roueff et al. 2005).

B 1, one of the highest column density cores in the Perseus Complex, with active low-mass star formation going on (e.g. Hirano et al. 1999) and with large deuterium fractionations, as shown by the detection of triply deuterated ammonia (ND3, Lis et al. 2002b,2006), doubly deuterated hydrogen sulfide (D2S, Vastel et al. 2003), and doubly deuterated thioformaldehyde (D2CS, Marcelino et al. 2005).

IRAM 04191, a very low luminosity Class 0 source in Taurus, driving a powerful outflow, but embedded in a dense core which appears to maintain many of the starless core characteristics (Belloche & André 2004).

L 1521F, initially selected as a starless core with a chemical and physical structure similar to L 1544 (but different kinematics; Crapsi et al. 2004), and recently found associated with a L < 0.07 $L_{\odot}$ protostellar object, thanks to the high sensitivity of the Spitzer Space Telescope (Bourke et al. 2006).

Ori B9, a massive dense core in Orion B, with peculiarly narrow molecular line widths and low gas temperature (Lada et al. 1991; Caselli & Myers 1995; Harju et al. 1993), thus an ideal target (among massive cores) to detect deuterated species.

NGC2264-VLA2, studied by, e.g., Ward-Thompson et al. (1995), who found the Class 0 driving source of the bipolar outflows, Girart et al. (2000), who determined the gas temperature and volume density of the surrounding core, and by Loinard et al. (2002), who observed D2CO and found an extremely large deuterium fractionation ([D2CO]/[H2CO] = 0.4, equivalent to an enrichment over the cosmic D/H ratio of more than 9 orders of magnitude).

\end{figure} Figure 1: ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) spectra toward the 16 sources of our sample. The units are main beam brightness temperature in K (y-axis, assuming a unity filling factor) and velocity in ${\rm km~s}^{-1}$ (x-axis). A star in the top right indicates dense cores associated with protostellar objects. The vertical dotted line is the $V_{\rm LSR}$ velocity measured with ${\rm N}_2{\rm H}^+$(1-0), whereas the vertical dashed line marks $V_{\rm LSR}$ as measured with the present ${\rm H}_{2}{\rm D}^{+}$ observations. Within the uncertainties, the two values are identical. Note that the strongest emission is present in the densest starless cores (L 1544, L 183, L 429, L 694-2) and in the star forming regions L 1521F, B 1, and NGC 2264G. Note also the large variation in linewidth among the various sources.
Open with DEXTER

3 Results

3.1 ortho-H2D+

The ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) spectra are shown in Fig. 1, and the results of Gaussian fits to the lines are listed in Table 2. One striking result from the figure and the table is the large variation in intensity (a factor of 5) and linewidths (a factor of 4), and not just between starless cores and protostellar cores. In Col. 5 of Table 2 we also report the non-thermal line width, defined as (see Myers et al. 1991):

$\displaystyle %
\Delta {v}_{\rm NT} = \sqrt{\Delta {v}_{\rm obs}^2 - \Delta {v}_{\rm T}^2},$     (2)

where $\Delta {v}_{\rm obs}$ ($\equiv $ $\Delta {v}$ in Table 2) is the observed line width and $\Delta {v}_{\rm T}$ is the thermal linewidth of the observed molecule, calculated assuming the kinetic temperatures listed in Col. 3 of Table 3. The ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) non-thermal line width is on average two times larger than that derived from the ${\rm N}_2{\rm D}^+$(2-1) data of Crapsi et al. (2005) (using the same kinetic temperature), the only exception being B 68, where both the ${\rm N}_2{\rm D}^+$(2-1) and ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) lines are totally thermally broadened. The larger non-thermal line widths (and the moderate optical depths, see Sect. 4) indicate that the ortho- ${\rm H}_{2}{\rm D}^{+}$ lines are tracing a region within the cores with more prominent internal motions than the ${\rm N}_2{\rm D}^+$(2-1) lines, even in L 1544, where the emission of these two lines has the same morphology and extension (Vastel et al. 2006a).

Table 2: Gaussian fits of the ortho- ${\rm H}_{2}{\rm D}^{+}$(11,0-11,1) lines in the source sample.

Table 3: ortho- ${\rm H}_{2}{\rm D}^{+}$column densities.

In the category of ``shallow'' cores, we have 3 non detections among 4 objects. The only shallow core detected in our survey is L 1517B, which is in fact the most compact and centrally concentrated of its class and has one of the narrowest ${\rm H}_{2}{\rm D}^{+}$ lines observed so far (together with the denser core Oph D): 0.4  ${\rm km~s}^{-1}$, only 1.2 times larger than the ${\rm H}_{2}{\rm D}^{+}$ thermal linewidth at 10 K, as measured with ${\rm NH}_{3}$ observations (Tafalla et al. 2004). Observations carried out with the APEX telescope revealed a probable ortho- ${\rm H}_{2}{\rm D}^{+}$ emission in one of our non-detected shallow cores, B 68 (Hogerheijde et al. 2006). The detected line is indeed relatively faint ( $T_{\rm mb}$ = 0.2 K), very narrow (0.3  ${\rm km~s}^{-1}$, practically thermal), and consistent with our non-detection ( $T_{\rm rms}$ = 0.3 K, see Table 2). The observations of the other two objects in this group (L 1498, TMC-2) have better sensitivities ( $T_{\rm rms}$ $\simeq$ 0.13 K), but they may still hide the ${\rm H}_{2}{\rm D}^{+}$ line if its intensity is similar to the one in B 68.

Among the five most centrally concentrated objects in the starless core sample, four show strong ( $T_{\rm mb}$ = 0.7-0.9 K) ${\rm H}_{2}{\rm D}^{+}$ emission, whereas TMC-1C has $T_{\rm mb}$ $\simeq$ 0.4 K. The line widths span the range between 0.4  ${\rm km~s}^{-1}$ (for Oph D) and 0.7  ${\rm km~s}^{-1}$ (for L 429), possibly reflecting contraction motions in different stages of core evolution or large optical depths (although a simple analysis seems to discard this last hypothesis; see Sect. 4).

In the six (young) protostellar cores, the detection rate is quite high, with four ${\rm H}_{2}{\rm D}^{+}$ lines detected. L 1521F has a line shape similar to that in L 1544 ( $\Delta {v}$ = 0.5  ${\rm km~s}^{-1}$), but the brightness temperature is 1.7 times lower, in agreement with the two times lower deuterium fractionation observed (N( ${\rm N}_2{\rm D}^+$)/N( ${\rm N}_2{\rm H}^+$) = 0.1 and 0.2, in L 1521F and L 1544, respectively; Caselli et al. 2002a; Crapsi et al. 2004). B 1 and NGC 1333- ${\rm DCO}^+$ present the largest linewidth in the sample, suggesting that the active star formation is probably injecting energy in the form of non-thermal motions and turbulence. We are confident that the broad line in NGC 1333- ${\rm DCO}^+$ is not a baseline artifact, in particular because both the centroid velocity and the linewidth are coincident (within the errors) with those observed in D2S by Vastel et al. (2003) (but not with the ND3 line observed by van der Tak et al. 2002, which remains a puzzle). The ${\rm H}_{2}{\rm D}^{+}$ line in NGC 2264G is narrower than in B 1 and NGC 1333- ${\rm DCO}^+$ and more similar to L 1521F, probably indicating that the circumstellar environment is still quite pristine.

3.2 para-D2H+

The para- ${\rm D}_{2}{\rm H}^{+}$(11,0-10,1) line has been searched for in four sources (B 1, NGC 1333- ${\rm DCO}^+$, NGC 2264G-VLA2, and L 183). Table 4 lists the spectral resolution ( $\Delta {v}_{\rm res}$, Col. 2), the system temperature ( $T_{\rm sys}$, Col. 3) and the integration time ( $t_{\rm int}$, Col. 4) of the observations. The rms noise and the upper limits of the radiation temperature (or brightness temperature, assuming a unity filling factor) are in Cols. 5 and 6, respectively. From these data, we calculated the corresponding upper limits of the para- ${\rm D}_{2}{\rm H}^{+}$ column density in each source, applying the same method as for ortho- ${\rm H}_{2}{\rm D}^{+}$ (see Sect. 4) and assuming a critical density for the transition of 105  ${{\rm cm}}^{-3}$, a line width equal to that measured for the ${\rm H}_{2}{\rm D}^{+}$ line (Table 2), except for 16293E (for which we used the para- ${\rm D}_{2}{\rm H}^{+}$ line width observed by Vastel et al. 2003), the kinetic temperature and ${\rm H}_2$ volume densities listed in Table 3, and the following parameters for the para- ${\rm D}_{2}{\rm H}^{+}$(11,0-10,1) transition: frequency $\nu_0$ = 691.660483 GHz, Einstein coefficient for spontaneous emission $A_{\rm ul}$ = 4.55 $\times $ 10-4 s-1, rotational constants as given in Table 4 of Amano & Hirao (2005), lower state energy of 50.2 K, and degeneracy of the upper and lower levels equal to 9.

Table 4: p- ${\rm D}_{2}{\rm H}^{+}$ column density upper limits.

The last column of Table 4 shows the ratio between the upper-limit column densities of para- ${\rm D}_{2}{\rm H}^{+}$ and the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities listed in Table 3. The estimated values are well within those calculated by Flower et al. (2004, see their Fig. 7) for cloud cores with volume densities 2 $\times $ 106  ${{\rm cm}}^{-3}$ and temperature ranges between 10 and 15 K.

4 Analysis

4.1 Derivation of the average ortho-H2D+ column densities

In this section we estimate the average ortho- ${\rm H}_{2}{\rm D}^{+}$ column density in each source, assuming that the line is emitted in homogeneous spheres at the density and temperature quoted in the literature and reported in Table 3. Values of the ${\rm H}_2$ column densities are also reported in Table 3, Col. 5, and they are used to determine the fractional abundance of ortho- ${\rm H}_{2}{\rm D}^{+}$ (see Sect. 4.3). The N(${\rm H}_2$) values for L 1544, L 429 and L 694-2 have been determined by Crapsi et al. (2005) assuming constant temperature. However, a temperature gradient has been measured toward L 1544 (Crapsi et al. 2007) and assumed (because of the similar physical structure) in L 429 and L 694-2, so that N(${\rm H}_2$) needs to be modified. As found by Pagani et al. (2004) and Stamatellos et al. (2007), the temperature drop in L 183 and Oph D implies an increase in N(${\rm H}_2$) by a factor of about 1.4. Given that L 1544, L 429 and L 694-2 have structures similar to L 183 and Oph D, we simply multiplied the Crapsi et al. (2005) values by the same correction factor (1.4) to account for the temperature gradient, as reported in Table 3.

We further assume that the level structure of the ortho- ${\rm H}_{2}{\rm D}^{+}$ molecule is reduced to a two-level system. This is especially justified in starless cores, considering that the first excited state is 17.4 K above ground and the second one is 110 K. In star forming regions, we assume that the observed ${\rm H}_{2}{\rm D}^{+}$ line arises from gas with characteristics not significantly different from starless cores (which is probably true in L 1521F and NGC 2264G-VLA2, see Sect. 3). In the case of B 1 and NGC 1333- ${\rm DCO}^+$, the broad lines suggest that the embedded young stellar objects have probably increased the degree of turbulence in the region and may have locally altered the conditions where ${\rm H}_{2}{\rm D}^{+}$ emits. However, we should point out that where the gas temperature increases above 20 K, and/or where shocks are present, ${\rm H}_{2}{\rm D}^{+}$ should not survive for long, considering that in these conditions the backward reaction (1) can quickly proceed and that dust mantles can be either evaporated or sputtered back into the gas phase, with the consequence of increasing the CO abundance and thus the destruction rate of  ${\rm H}_{2}{\rm D}^{+}$.

To estimate the average ortho- ${\rm H}_{2}{\rm D}^{+}$ column density, we evaluate the excitation temperature $T_{\rm ex}$ and the line optical depth $\tau$ simultaneously, by solving iteratively the following equations (valid for $T_{\rm bg} \ll E_{\rm ul}/k$):

$\displaystyle %
T_{\rm ex} = \frac{T_{\rm kin}}{\frac{T_{\rm kin}}{E_{\rm ul}}
\ln \left( 1 + \frac{\beta n_{\rm cr}}{n_{\rm tot}} \right) +1}$     (3)

$\displaystyle %
\beta = \frac{0.75}{\tau}\left(1+\frac{{\rm e}^{-2\tau}}{\tau}+
\frac{({\rm e}^{-2\tau}-1)}{2\tau^2}\right)$     (4)

$\displaystyle %
\tau = - \ln \left[ 1 - \frac{T_{\rm mb}}{J_{\nu}(T_{\rm ex}) - J_{\nu}(T_{\rm bg})} \right]$     (5)

Eq. (3) is the definition of $T_{\rm ex}$ corrected for the line opacity: in practice, the critical density is reduced by the probability that the emitted photon can indeed escape absorption. We adopted the photon escape probability $\beta$ of Eq. (4), valid for a homogeneous sphere (Osterbrock 1989). Finally the line optical depth $\tau$ is derived from the observation via Eq. (5). In the above equations, $E_{\rm ul}$ is the energy of the transition ( $E_{\rm ul}/k_{\rm B}$ = 17.4 K, with $k_{\rm B}$ the Boltzmann constant), $n_{\rm tot}$ is the particle volume density (assumed to be 1.2 $\times $ n(${\rm H}_2$), to account for He), $T_{\rm kin}$ is the gas temperature, and $J_{\nu}(T_{\rm ex})$ and $J_{\nu}(T_{\rm bg})$ are the equivalent Rayleigh-Jeans excitation and background temperatures.

The ortho- ${\rm H}_{2}{\rm D}^{+}$ column density is then derived from $\tau$:

$\displaystyle %
N({\rm ortho{-}H_2D^+}) = \frac{8\pi\nu^3}{c^3}\frac{Q(T_{\rm e...
...\frac{{\rm e}^{E_u/T_{\rm ex}}}{{\rm e}^{E_u/T_{\rm ex}}-1}
\int \tau {\rm d}v.$     (6)

A key parameter is the critical density $n_{\rm cr}$ of the transition. Recent unpublished calculations by Hugo & Schlemmer (private communication) suggest a collisional coefficient with ortho and para ${\rm H}_2$ of $\sim$10-9 cm3 s-1 at 10 K, with a very shallow dependence on the temperature. The implied critical density is $\sim$105 cm-3, namely a factor of 10 lower than assumed in previous work (Black et al. 1990; van der Tak et al. 2005). The ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities, calculated both with the new and old critical density values, are reported in Table 3. The lines are either optically thin (e.g. in B 68) or marginally thick ($\simeq$1, in L 1544, L 183, L 694-2, L 429). From the table, it is clear that a factor of 10 variation in the collisional coefficient implies a change in the column density value of factors between 1.4 and 4.5, depending on the volume density and kinetic temperature. In the following analysis we use the column densities calculated with the 105 cm-3 critical density.

The four objects that show the largest values of N(ortho- ${\rm H}_{2}{\rm D}^{+}$) are among the most centrally concentrated cores in the sample: L 429, L 1544, L 694-2 and L 183. The two other dense cores, TMC-1C and Oph D have significantly lower ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities, and this may be related to different evolutionary stages. We will further discuss these issues in Sect. 5.

4.2 Evaluation of the errors

The estimates of the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities reported in Table 3 suffer from several sources of uncertainties. We already mentioned a basic source of uncertainty, that associated with the collisional coefficient of the transition. In addition to that, the densities and temperatures used to derive the excitation temperatures are also relatively uncertain, not only because of the uncertainty in deriving these values at the centers of the studied sources but also because the ${\rm H}_{2}{\rm D}^{+}$ line emission may originate in regions that are denser than the quoted average density gas due to the ${\rm H}_{2}{\rm D}^{+}$ abundance distribution. In this context, the errors associated with the rms of the observations reported in Table 2 are certainly the smallest in the error propagation chain. Although it is difficult to exactly quantify the error in the determination of the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities, we estimate here how reasonable changes in the gas temperature and density would affect the reported column densities. Increasing or decreasing the density by a factor of 2 results in decreasing/increasing the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities by less than 30%. However, a change in the kinetic temperatures of Table 3 by 1 K would change the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities by up to a factor of 2 in the coldest objects (because of the exponential in the level population equation).

In summary, considering also the beam efficiency variation between 0.4 and 0.7 (Table 2), the ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities reported in Table 3 are likely to be uncertain by about a factor of $\simeq$2.

4.3 Correlations

We have looked for possible correlations between the column density or fractional abundance of ortho- ${\rm H}_{2}{\rm D}^{+}$ and physical parameters such as the volume density, the ${\rm H}_2$ column density, the kinetic temperature and the non-thermal line width. No significant correlations have been found, with the exception of x(ortho- ${\rm H}_{2}{\rm D}^{+}$) ($\equiv $N(ortho- ${\rm H}_{2}{\rm D}^{+}$)/N(${\rm H}_2$), with N(${\rm H}_2$) from Table 3) vs. $T_{\rm kin}$, for which we find (see Fig. 2):

$\displaystyle %
{\rm Log}~ x({\rm ortho}{-}{\rm H}_{2}{\rm D}^{+}) = (-8.6 \pm 0.3) - (0.16 \pm 0.03) T_{\rm kin},$     (7)

with a correlation coefficient of -0.83. What is causing the ortho- ${\rm H}_{2}{\rm D}^{+}$ drop observed between 7 K and 15 K? Partially, this may be due to the fact that the warmest sources: (i) may have a smaller ortho- ${\rm H}_{2}{\rm D}^{+}$ emitting region than found in L 1544 (about 60 $^{\prime\prime}$; Vastel et al. 2006a), because of a drop of the ortho- ${\rm H}_{2}{\rm D}^{+}$ abundance (see Fig. 6 of Flower et al. 2004) and (ii) are all at distances >300 pc (with the exception of 16293E, which follows the trend). Thus, the ${\rm H}_{2}{\rm D}^{+}$ emission may be affected by beam dilution (not considered in this study). However, as we will see in the next section, variations in the CO depletion factor, linked to the gas density and, in turn, to the gas and dust temperature (at least in starless cores, larger temperatures are typically associated with lower gas densities, lower CO depletion factors, and lower deuterium-fractionations), can also (at least partially) produce the observed trend. In any case, a more detailed physical and chemical structure (such as that recently developed by Aikawa et al. 2008) is needed to investigate these points, although the lack of information on the extent and morphology of the ${\rm H}_{2}{\rm D}^{+}$ emission prevents us from a more quantitative analysis of the correlations between the gas traced by ${\rm H}_{2}{\rm D}^{+}$ and the physical properties of the selected cores.

\end{figure} Figure 2: ortho- ${\rm H}_{2}{\rm D}^{+}$ fractional abundances vs. the gas volume density (n(${\rm H}_2$), left panel) and gas temperature ( $T_{\rm kin}$, right panel). Empty symbols refer to starless cores, whereas filled symbols refer to cores associated with young stellar objects. Note that upside-down triangles are upper limits. The line in the right panel is the least square fits to the x(ortho- ${\rm H}_{2}{\rm D}^{+}$) vs. $T_{\rm kin}$ data, the only significant correlation that has been found (see text).
Open with DEXTER

In the following we will concentrate on molecular abundances and use both a simple chemical model applied to a homogeneous cloud and a slightly more detailed chemical-physical model of a centrally concentrated and spherically symmetric cloud to reproduce and interpret the observed variations of ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities, deuterium fractionation and CO depletion fraction in the selected cores.

Table 5: The forward rate coefficient is given by $\alpha (T/300)^{\beta }$. The reverse rate is given by $\alpha (T/300)^{\beta }{\rm e}^{-\gamma /T}$.

5 Chemical discussion

The chemistry of starless cores and their degree of deuteration has been investigated in detail by Roberts et al. (2003), Flower et al. (2004), Roberts et al. (2004), Walmsley et al. (2004), Aikawa et al. (2005), Flower et al. (2005), Flower et al. (2006a), Flower et al. (2006b). From these models it appears that there are several sources of uncertainty that can profoundly affect the chemistry: (i) surface chemistry (diffusion and reaction rates are still quite uncertain and highly dependent upon the poorly known surface of dust grains); (ii) freeze-out and desorption rates (binding energies may be changing throughout the core, due to changes in grain mantle composition, and nonthermal desorption processes are poorly known; see Garrod et al. 2007,2006); (iii) the dynamical evolution of dense cores is hard to constrain chemically, and different theoretical models of star formation predict significantly different time scales (e.g. Shu et al. 1987; Hartmann et al. 2001); (iv) the cosmic-ray ionization rate $\zeta $ is not well constrained, and our ignorance of the cosmic-ray energy spectrum (especially at low energies) prevents us from making quantitative estimates of the possible variations of $\zeta $ within dense cores (see Padoan & Scalo 2005, for a recent discussion on this point; and Dalgarno 2006, for a more general review); (v) the fraction of ortho-${\rm H}_2$, which affects the deuteration in the gas phase (the backward reaction (1) proceeds faster with ortho-${\rm H}_2$, because of its higher energy compared to para-${\rm H}_2$; Gerlich et al. 2002, hereafter GHR02; Walmsley et al. 2004); (vi) the dust grain size distribution (if dust grains coagulate in the densest regions of starless cores, the freeze-out rate diminishes, altering the gas phase chemical composition; see, e.g. Vastel et al. 2006a; Flower et al. 2006a); (vii) the abundance of polycyclic aromatic hydrocarbons (PAHs), unknown in dark clouds (where observational constraints are yet to be found), which may significantly affect the electron fraction (Lepp & Dalgarno 1988; Flower & Pineau des Forêts 2003; Tielens 2005); (viii) the fraction of atomic oxygen in the gas phase, which is thought to be low (mainly to explain the stringent SWAS upper limits on the water abundance; e.g. Bergin & Snell 2002), but observations of dark clouds with the ISO satellite appear to disprove this (Lis et al. 2001; Vastel et al. 2000; Caux et al. 1999), although the limited velocity resolution of ISO LWS is a severe limit on these results. In the following section, the effects of some of the above parameters on the deuterium fractionation will be presented for the simple case of a homogeneous cloud. In Sect. 5.2, a more detailed physical structure and a slightly more comprehensive chemical model will be considered to make an attempt at constraining some of the unknown parameters.

5.1 Simple theory: what affects deuterium fractionation

Ignoring for the moment the density and temperature structure of molecular cloud cores and any gas-dust interaction, which will be considered in Sect. 5.2, we show here simple relations between the ${\rm H}_{2}{\rm D}^{+}$/ ${\rm H}_{3}^{+}$ abundance ratio and parameters such as the gas kinetic temperature, the grain size distribution, the CO depletion factor and the cosmic-ray ionization rate. We use a simple chemical network which includes all the multiply deuterated forms of ${\rm H}_{3}^{+}$, formed following the reaction scheme listed in Table 5 and destroyed by CO, electrons (dissociative recombination) and negatively charged dust grains (recombination). The adopted rate coefficients are the same as in Table 1 of Ceccarelli & Dominik (2005), with the exception of the reaction of ${\rm H}_{3}^{+}$ (and deuterated isotopologues) with CO, and the reaction of ${\rm H}_{3}^{+}$ and electrons, for which we used the values listed in the latest release of the UMIST database (RATE06) available at http://www.udfa.net/. We note that the rate coefficient of the ${\rm H}_{3}^{+}$ + CO reaction in the UMIST database is temperature-independent, as expected for ion-molecule reactions where the neutral species has a small dipole moment (Herbst, private communication).

\end{figure} Figure 3: [ ${\rm H}_{2}{\rm D}^{+}$]/[ ${\rm H}_{3}^{+}$], [ ${\rm D}_{2}{\rm H}^{+}$]/[ ${\rm H}_{3}^{+}$], [ ${\rm D}_{3}^{+}$]/[ ${\rm H}_{3}^{+}$] and $R_{\rm D}$ abundance ratios as a function of the gas temperature ( $T_{\rm kin}$) in a dense cloud with uniform volume density n(${\rm H}_2$) = 105  ${{\rm cm}}^{-3}$. ( Top row) The abundance ratios are plotted against $T_{\rm kin}$ for different values of the depletion factor $f_{\rm D}$ (=1, 5, 10, 50 and 100), with fixed values of $a_{\rm min}$ = 50 Å and $\zeta $ = 3 $\times $ 10-17 s-1. The dotted curves are for $f_{\rm D}$ = 10 clouds with n(${\rm H}_2$) = 1 $\times $ 106  ${{\rm cm}}^{-3}$ (top dotted curve) and n(${\rm H}_2$) = 1 $\times $ 104  ${{\rm cm}}^{-3}$ (bottom dotted curve). Note the large increase of the [ ${\rm D}_{3}^{+}$]/[ ${\rm H}_{3}^{+}$] ratio for $f_{\rm D}$ = 50 and 100. ( Central panel) Abundance ratios vs. $T_{\rm kin}$ for different values of $a_{\rm min}$ and values of $f_{\rm D}$ and $\zeta $ fixed at 10 and 3 $\times $ 10-17 s-1, respectively. ( Bottom panel) Abundance ratios vs. $T_{\rm kin}$ for different values of the cosmic-ray ionization rate, $\zeta $. Here, $f_{\rm D}$ = 10 and $a_{\rm min}$ = 50 Å.
Open with DEXTER

The steady state equations for this simple system are:

    $\displaystyle \frac{x({\rm H}_{2}{\rm D}^{+})}{x({\rm H}_{3}^{+})} = \frac{k_1 ...
...HD}) k_{-3} k_3]}
{D_1 D_2 D_3 - x({\rm HD}) (k_{-3} k_3 D_1 - k_{-2} k_2 D_3)}$ (8)
    $\displaystyle \frac{x({\rm D}_{2}{\rm H}^{+})}{x({\rm H}_{2}{\rm D}^{+})} = \frac{k_2 x({\rm HD}) D_3}{D_2 D_3 - k_{-3} k_3 x({\rm HD})}$ (9)
    $\displaystyle \frac{x({\rm D}_{3}^{+})}{x({\rm D}_{2}{\rm H}^{+})} = \frac{k_3 x({\rm HD})}{k_{-3} + k_{\rm co} x({\rm CO}) + k_{\rm rec3} x(e) + k_g x(g)},$ (10)


\begin{eqnarray*}&& D_1 = k_{-1} + k_{\rm co} x({\rm CO}) + k_{\rm rec1} x(e) + ...
... k_{-3} + k_{\rm co} x({\rm CO}) + k_{\rm rec3} x(e) + k_g x(g).

In the above expressions, k1, k2, k3 and k-1, k-2, k-3 are the forward and backward rate coefficients relative to reactions of ${\rm H}_{3}^{+}$, ${\rm H}_{2}{\rm D}^{+}$, and ${\rm D}_{2}{\rm H}^{+}$, respectively, with HD; x(i) is the fractional abundance (w.r.t. ${\rm H}_2$) of species i; $x({\rm HD})$ = 3 $\times $ 10-5 is the fractional abundance of HD (assuming that [D]/[H] = 1.5 $\times $ 10-5, Oliveira et al. 2003, and that in molecular clouds the deuterium is mainly in the form of HD); $k_{\rm CO}$ is the rate coefficient of the destruction reaction of all ${\rm H}_{3}^{+}$ isotopologues with CO; $x({\rm CO})$ = $x_{\rm can}(\rm CO)/f_{\rm D}$, where $x_{\rm can}(\rm CO)$ = 9.5 $\times $ 10-5 is the canonical abundance of CO as measured by Frerking et al. (1982), and $f_{\rm D}$ is the CO depletion factor (1/$f_{\rm D}$ is the fraction of CO molecules left in the gas phase, see Sect. 5.2.2); $k_{\rm rec1}$, $k_{\rm rec2}$, $k_{\rm rec3}$ are the dissociative recombination rate coefficients for ${\rm H}_{2}{\rm D}^{+}$, ${\rm D}_{2}{\rm H}^{+}$, and ${\rm D}_{3}^{+}$, respectively. Following Draine & Sutin (1987; see also Crapsi et al. 2004), the rate coefficient for the recombination onto dust grains has the form:
$\displaystyle %
k_{\rm g} = 1.6\times 10^{-7} \left( \frac{a_{\rm min}}{10^{-8}...
...-4} \frac{T_{\rm kin}}{10~{\rm K}}
\frac{a_{\rm min}}{10^{-8}~{\rm cm}} \right)$     (11)

where $a_{\rm min}$ (=50 Å) is the minimum radius of dust grains in the Mathis et al. (1977) (MRN) size distribution ( $a_{\rm max}$ = 2.5 $\times $ 10-5 cm). Finally, the fractional abundance ($x_{\rm g}$) of dust grains in a MRN size distribution is given by:
$\displaystyle %
x_{\rm g} = 5.3 \times 10^{-6} \left( \frac{a_{\rm max}}{10^{-4...
... \right)^{-0.5} \left( \frac{a_{\rm min}}{10^{-8}~{\rm cm}} \right)^{-2.5}\cdot$     (12)

The electron fraction is calculated as in Caselli et al. (2002b) using a simplified version of the reaction scheme of Umebayashi & Nakano (1990) where we compute a generic abundance of molecular ions ``mH+'' assuming formation due to proton transfer with ${\rm H}_{3}^{+}$ and destruction by dissociative recombination on grain surfaces (using rates from Draine & Sutin 1987). ${\rm H}_{3}^{+}$ is formed as a consequence of cosmic-ray ionization of ${\rm H}_2$ and destroyed by proton transfer with CO and ${\rm N}_2$. Metals are also taken into account and their fractional abundance has been assumed to be 10-7 (following the initial abundances of ``low-metal'' models, appropriate for dark clouds; Lee et al. 1996).

The deuterium fractionation in species such as ${\rm HCO}^+$ or ${\rm N}_2{\rm H}^+$ ($R_{\rm D}$ $\equiv $ [ ${\rm DCO}^+$]/[ ${\rm HCO}^+$] or [ ${\rm N}_2{\rm D}^+$]/[ ${\rm N}_2{\rm H}^+$]) in this chemical scheme is simply given by:

$\displaystyle %
R_{\rm D} = \frac{1/3 x({\rm H}_{2}{\rm D}^{+}) + 2/3 x({\rm D}...
..._{3}^{+}) + 2/3 x({\rm H}_{2}{\rm D}^{+}) + 1/3 x({\rm D}_{2}{\rm H}^{+})}\cdot$     (13)

Figure 3 shows the [ ${\rm H}_{2}{\rm D}^{+}$]/[ ${\rm H}_{3}^{+}$], [ ${\rm D}_{2}{\rm H}^{+}$]/[ ${\rm H}_{3}^{+}$], [ ${\rm D}_{3}^{+}$]/[ ${\rm H}_{3}^{+}$] and $R_{\rm D}$ abundance ratios as a function of gas temperature, for different values of the depletion factor $f_{\rm D}$ (top panel), the minimum value of the dust grain radius in the grain-size distribution, $a_{\rm min}$ (central panel), and the cosmic-ray ionization rate, $\zeta $ (bottom panel). The first thing to note is the sharp drop in the deuterium fractionation at temperatures above $\simeq$15 K, when reaction (1) and the analogous ones for the formation of ${\rm D}_{2}{\rm H}^{+}$ and ${\rm D}_{3}^{+}$ also start to proceed backward. Between 5 and 10 K, the deuteration ratios increase because of the inverse temperature dependence of the ${\rm H}_{2}{\rm D}^{+}$ destruction rate coefficients $k_{\rm rec1}$, $k_{\rm rec2}$, $k_{\rm rec3}$ and kg (present in the denominator of Eqs. (8)-(10)). Note the high peaks in the $f_{\rm D}$ = 50 and 100 [ ${\rm D}_{3}^{+}$]/[ ${\rm H}_{3}^{+}$] curves. This is due to the fact that in regions where most of the neutrals are frozen (although here we only talk about CO, we can generalize this statement including in the CO-depletion factor all the neutral species which react with the ${\rm H}_{3}^{+}$ isotopologues), the deuterium fractionation proceeds rapidly and ${\rm H}_{3}^{+}$ is efficiently converted into ${\rm D}_{3}^{+}$, as already shown by Walmsley et al. (2004) in the case of dense cloud cores and by Ceccarelli & Dominik (2005) in the case of protoplanetary disks. The large abundances of ${\rm D}_{3}^{+}$ yield large $R_{\rm D}$ values ($\geq$1).

Figure 3 also shows that at temperatures below $\simeq$17 K, the deuterium fractionation is very much dependent upon the CO depletion factor (a well known phenomenon; e.g. Dalgarno & Lepp 1984), the gas volume density (see the dotted curves in the top figures), the fractional abundance of dust grains and the cosmic-ray ionization rate. In particular, a value of $a_{\rm min}$ = 5 $\times $ 10-8 cm (5 Å) corresponds to $x_{\rm g}$ = 2 $\times $ 10-7, which may be regarded as a possible value for the fractional abundance of PAHs (e.g. Lepp & Dalgarno 1988; Tielens 2005). Thus, if PAHs are abundant in dense cores and they are the main negative charge carriers, the deuterium fractionation is expected to be $\leq$0.1, because the four ${\rm H}_{3}^{+}$ isotopologues quickly recombine. The fact that molecules such as ${\rm N}_2{\rm H}^+$ and ${\rm NH}_{3}$ show large deuterium fractionations (or large $R_{\rm D}$ values) in the direction of pre-stellar cores and Class 0 sources ($R_{\rm D}$ > 0.1; see Fig. 4) thus suggests that (negatively charged) PAHs have abundances significantly below 10-7.

The bottom panel of Fig. 3 shows that the larger the cosmic-ray ionization rate the smaller the deuteration ratios. This is mainly due to the fact that a larger value of $\zeta $ implies a larger electron fraction, and a consequently larger dissociative recombination rate (see denominator of Eqs. (8)-(10)). Again, the large deuterium fractionation observed in pre-stellar cores and Class 0 objects can be used to put upper limits on $\zeta $ (see Dalgarno 2006).

The $R_{\rm D}$ values obtained in this analysis for typical parameters ((i) $f_{\rm D}$ $\simeq$ 10, as typically observed in pre-stellar cores; (ii)  $a_{\rm min}$ = 50 Å, as in the MRN distribution; and (iii) $\zeta $ $\simeq$$\times $ 10-17 s-1) reach $\simeq$0.5 for $T_{\rm kin}$ $\leq$ 15 K. This value may appear too large when compared to the deuterium fractionation measured in pre-stellar cores by Crapsi et al. (2005), but it is quite close to the [NH2D]/[ ${\rm NH}_{3}$] ratio found by (1) Crapsi et al. (2007) in the nucleus of the pre-stellar core L 1544 using interferometric observations and by (2) Pillai et al. (2007) in Infrared Dark Clouds. However, the results presented in this section apply to an ideal homogeneous cloud and are based on ``standard'' rate coefficients for the proton-deuteron exchange reactions ${\rm H}_{3}^{+}$, ${\rm H}_{2}{\rm D}^{+}$, ${\rm D}_{2}{\rm H}^{+}$ + HD. In fact, GHR02 have measured slower rates which, if adopted, lead to $R_{\rm D}$ values about a factor of 3 lower compared to those obtained in Fig. 3. This will be discussed in the next sub-section.

\end{figure} Figure 4: ortho- ${\rm H}_{2}{\rm D}^{+}$ column density as a function of the observed deuterium-fractionation $R_{\rm D}$ (see text for details). Filled circles are protostellar cores, empty circles are starless cores. Downward arrows denote upper limits in our estimate of the column density. The names of the cores associated with each mark in the figure are shown in the top panel. The bottom panel shows the same plot, where theoretical predictions from chemical models (see text) are superposed. Solid curves are for models that use the ``standard'' rate coefficients for the proton-deuteron exchange reactions, whereas dashed curves are for models adopting the smaller GHR02 rate coefficients. Each curve has six points (filled (empty) diamonds for the models using the ``standard'' (GHR02) rate coefficients), which correspond to model cores with different central density (from 4 $\times $ 104 cm-3 to 4 $\times $ 108 cm-3, from left to right). Different curves are for different ortho/para ${\rm H}_{2}{\rm D}^{+}$ ratios, from 0.03 to 2.0. The cosmic-ray ionization rate has been fixed to 1.3 $\times $ 10-17 s-1.
Open with DEXTER

5.2 Simple chemical-physical model

In this section, our estimates of ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities are correlated with the deuterium fractionation and the depletion factor, previously measured in the same objects, and simple chemical models are used to investigate the observed variations among the various sources. The model used is similar to that described in Vastel et al. (2006a) and first applied by Caselli et al. (2002b) in the case of L 1544, but with the deuterium fractionation chemistry and rate coefficients as described in the previous sub-section. We consider a spherical cloud with a given density and temperature profile where dust and gas are present. Initially (besides ${\rm H}_2$), CO and N2 are present in the gas phase with abundances of 9.5 $\times $ 10-5 (Frerking et al. 1982) and 3.75 $\times $ 10-5, respectively. The abundance of ${\rm N}_2$ assumes that 50% of the nitrogen is in the atomic form (but no nitrogen chemistry is considered, except for the ${\rm N}_2$ adsorption/desorption onto/from dust grains and the formation and destruction of ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$). This is a totally arbitrary choice, but the ${\rm N}_2$ abundance is extremely uncertain (see e.g. Stepnik et al. 2003; Flower et al. 2006b; Maret et al. 2006) and in any case, its variation in the gas phase only affects, in our simple model, the absolute abundance of ${\rm N}_2{\rm H}^+$, without significantly affecting the ${\rm N}_2{\rm D}^+$/ ${\rm N}_2{\rm H}^+$ column density ratio. Atomic oxygen has not been included in the chemistry in order to avoid one extra (uncertain) parameter in the model. Adding atomic oxygen to the chemical network lowers the deuterium fractionation if its binding energy (also not well constrained) is sufficiently low (see e.g. Caselli et al. 2002b, and their discussion in Sect. 3.2).

The dust grain distribution follows MRN and we assume a gas-to-dust mass ratio of 100. However, the size of the minimum dust radius, $a_{\rm min}$, has been increased by an order of magnitude, following recent (indirect) evidence of grain growth toward the center of dense cores (e.g. Flower et al. 2005,2006b; Vastel et al. 2006a; Bergin et al. 2006). The higher $a_{\rm min}$ adopted here (5 $\times $ 10-6 cm) lowers the freeze-out rate by a factor of 5 compared to the MRN value (by changing the surface area of dust grains), and it is the ``best-fit'' value for the L 1544 chemical model (see Vastel et al. 2006a; Caselli et al. 2002b). The freeze-out time scale of species i ( $R_{\rm freeze}(i)$) is given by:

                             $\displaystyle %
t_{\rm freeze}(i)$ = $\displaystyle \frac{1}{S ~ \Sigma ~ \langle v_{\rm th}(i) \rangle n_{\rm H}}$  
  = $\displaystyle \frac{2 \times 10^{4}~{\rm yr}}{S} \times
\left[ \left( \frac{a_{...
...-0.5} -
\left( \frac{a_{\rm max}}{10^{-5}~{\rm cm}} \right)^{-0.5} \right]^{-1}$  
    $\displaystyle \times~ A(i)^{0.5} \left( \frac{10~{\rm K}}{T_{\rm gas}} \right)^{0.5}
\left( \frac{10^5~{\rm cm^{-3}}}{n_{\rm H}} \right)$ (14)

where S is the sticking coefficient ($\simeq$1, as recently found by Bisschop et al. 2006), $\langle v_{\rm th}(i) \rangle$ is the average thermal velocity of species i, $n_{\rm H}$ is the total number density of H nuclei ($n_{\rm H}$ = n(H) + 2n(${\rm H}_2$)), A(i) is the atomic mass number of species i, and
$\displaystyle %
\Sigma = (4.88\times 10^{-25} + 4.66\times 10^{-25}) \times
(a_{\rm min}^{-0.5} - a_{\rm max}^{-0.5})$     (15)

is the grain surface area per H nucleon in the MRN distribution (see also Weingartner & Draine 2001). The electron fraction is calculated as in Sect. 5.1.

The binding energies for CO and ${\rm N}_2$ have been taken from Öberg et al. (2005), assuming that the mantle composition is a mixture of CO and H2O ice ($E_{\rm D}$(CO)/$k_{\rm B}$ = 1100 K and $E_{\rm D}$(${\rm N}_2$) = 0.9 $E_{\rm D}$(CO)). The cosmic-ray ionization rate used here is $\zeta $ = 1.3 $\times $ 10-17 s-1, but we have also changed it to explore the effects on the chemistry (see next subsections). Different density structures have been considered (see below) and the (gas = dust[*]) temperature profile is similar to the one found by Young et al. (2004) and parametrized so that:

$\displaystyle %
T_{\rm dust} \sim T_{\rm gas} \simeq 3 \times [8.7 - \log~(n({\rm H}_2)] ~~ {\rm K}.$     (16)

The minimum (maximum) allowed temperature is 4 K (14 K). We also consider models with temperature profiles increasing inwards, to simulate the heating of the embedded young stellar object. In this case, we assume a central temperature of 50 K at a distance of 80 AU, and a radial dependence r-0.6 for r > 80 AU. If the temperature drops below the one described in Eq. (16), the latter value is used.

CO and ${\rm N}_2$ can freeze out (with rates given by 1/ $t_{\rm freeze}$, see Eq. (14)) and return to the gas phase via thermal desorption or cosmic-ray impulsive heating (following Hasegawa et al. 1992; Hasegawa & Herbst 1993). The abundance of the molecular ions ( ${\rm N}_2{\rm H}^+$, ${\rm HCO}^+$, ${\rm H}_{3}^{+}$ and their deuterated isotopologues) are calculated in terms of the instantaneous abundances of neutral species (assumption based on the short time scale of ion chemistry, compared to the depletion time scale; see Caselli et al. 2002b, for details).

5.2.1 N(ortho-H2D+) vs. the observed R ${_{\rm D}}$

In Fig. 4, the column density of ortho- ${\rm H}_{2}{\rm D}^{+}$ is plotted as a function of the observed deuterium fractionation ratio ($R_{\rm D}$). $R_{\rm D}$ is equivalent to N( ${\rm N}_2{\rm D}^+$) /N( ${\rm N}_2{\rm H}^+$), in the case of starless cores (plus L 1521F), and this value has been taken from the survey of Crapsi et al. (2005). In star forming regions, the ${\rm N}_2{\rm D}^+$/ ${\rm N}_2{\rm H}^+$ column density ratio is not available, and other column density ratios have been used: (i)  $N({\rm NH_2D})/N({\rm NH_3})$ for NGC 1333- ${\rm DCO}^+$ (Hatchell 2003) and for B 1 (Roueff et al. 2005); (ii)  $\sqrt{N({\rm D_2CO})/N({\rm H_2CO})}$ for NGC 2264G VLA2 (Loinard et al. 2002). No deuterium fractionation estimates are available for IRAM 04191. Given that ${\rm NH}_{3}$ and ${\rm N}_2{\rm H}^+$ appear to trace similar zones of dense cores (e.g. Benson et al. 1998; Caselli et al. 2002c), and that both derive from the same parent species (${\rm N}_2$), one expects that the D-fractionation observed in the two species is also similar (and linked to the theoretical $R_{\rm D}$ in Eq. (13)). However, it is not obvious that formaldehyde is actually tracing the same region (indeed H2CO is centrally depleted in the two starless cores studied by Tafalla et al. 2006, unlike ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$). Figure 4 shows that the deuterium-fractionation in NGC 2264G is the largest one in the whole sample, probably suggesting that different deuteration mechanisms (beside the ${\rm H}_{3}^{+}$ fractionation) may be at work for H2CO. One possibility is that surface chemistry is needed to explain the observed amount of deuterated formaldehyde and methanol, as originally discussed by Charnley et al. (1997), Ceccarelli et al. (1998), and more recently by Parise et al. (2006).

In the top panel, each data point is labelled with the corresponding name, whereas the bottom panel shows the same data points with model curves superposed (see below). The first thing to note is that, on average, dense protostellar cores (filled symbols) have lower N(ortho- ${\rm H}_{2}{\rm D}^{+}$) values than starless cores, but on average they show quite large deuterium fractionations (especially in the case of NGC 2264G VLA2, where $R_{\rm D}$ is from measurements of doubly deuterated formaldehyde, as already mentioned). Another thing to note is that there is not any clear correlation between N(ortho- ${\rm H}_{2}{\rm D}^{+}$) and the observed $R_{\rm D}$. To investigate this unexpected result, we used the model described in Sect. 5.2 and simulate an evolutionary sequence, similar to what was done in Crapsi et al. (2005) in their Fig. 5.

We consider Bonnor-Ebert (BE) spheres with density structures analogous to the radial (cylindrical) density profile of the contracting disk-like cloud at different stages of evolution in the model of Ciolek & Basu (2000), namely those at times t = t1 (=2.27 Myr and central densities $n_{\rm c1}$(${\rm H}_2$) = 4 $\times $ 104  ${{\rm cm}}^{-3}$), t2 (=2.60 Myr, and $n_{\rm c2}$(${\rm H}_2$) = 4 $\times $ 105  ${{\rm cm}}^{-3}$), t3 (=2.66 Myr, and $n_{\rm c3}$(${\rm H}_2$) = 4 $\times $ 106  ${{\rm cm}}^{-3}$), t4 (=2.68 Myr, and $n_{\rm c4}$(${\rm H}_2$) = 4 $\times $ 107  ${{\rm cm}}^{-3}$), and t5 (=2.684 Myr, and $n_{\rm c5}$(${\rm H}_2$) = 4 $\times $ 108  ${{\rm cm}}^{-3}$). The BE density profile is reasonably well reproduced by the parametric formula (Tafalla et al. 2002):

$\displaystyle %
n(r) = \frac{n_{\rm c}({\rm H}_2)}{(1 + (r/r_0)^{\alpha})}$     (17)

where r0 = 13 000, 3000, 800, 300, and 80 AU for t1, t2, t3, t4, and t5, respectively, and $\alpha$ = 2. The five chemical models (i.e. the CO and N2 depletion chemistry in the five model clouds with different physical structure) have been run for a time interval given by (ti - t1), for i = 2-5, whereas the t1 model has been run for 2.27 Myr. We also consider a model cloud with the same density structure as the t5 model, but with a temperature profile resembling that of a centrally heated protostellar core (model t5a), with a central temperature of 50 K and a temperature gradient proportional to r-0.6 (as mentioned above). The abundance profiles of ${\rm H}_{2}{\rm D}^{+}$, ${\rm N}_2{\rm H}^+$ and ${\rm N}_2{\rm D}^+$ calculated by the models have been convolved with 22 $^{\prime\prime}$, 26 $^{\prime\prime}$, and 17 $^{\prime\prime}$ FWHM antenna beams, respectively, to simulate the observations and calculate the column densities toward the core center (as well as off peak).

\end{figure} Figure 5: ( Top panel) Same models as in the bottom panel of Fig. 4 but with a fixed value o/p- ${\rm H}_{2}{\rm D}^{+}$ (=0.5), and various cosmic-ray ionization rates (from 6 to 30 $\times $ 10-18 s-1). The thick dashed curve connects model results using GHR02 rates, $\zeta $ = 6 $\times $ 10-18 s-1 and o/p- ${\rm H}_{2}{\rm D}^{+}$ = 0.3. This is an attempt to reproduce the low ortho- ${\rm H}_{2}{\rm D}^{+}$ column density values observed toward the $R_{\rm D}$-rich protostellar cores and the two Ophiuchus pre-stellar cores Oph D and 16293E. ( Bottom panel) Plot showing predictions of models with density structures equal to the Ciolek & Basu (2000) cloud at times ti (i = 2-5) and corresponding central densities $n_{ci} \equiv n_i$ (see text) at different stages of chemical evolution (at t = 103, 104, 105, and 106 yr, from bottom left to top right). The o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio has been fixed at 0.5. The thick curves refer to the model of a protostellar envelope at different evolutionary times, where o/p- ${\rm H}_{2}{\rm D}^{+}$ = 0.1, to simulate a possible (slight) increase of the gas temperature due to the central heating source.
Open with DEXTER

The results of these models are the small diamonds in each of the curves in the bottom panel of Fig. 4, with t1 lying on the left-end and $t_{\rm 5a}$ on the right-end of the curve. Solid curves represent models with standard rate coefficients for the proton-deuteron exchange reactions, whereas dashed curves models use the about 3 times smaller GHR02 rates (see previous sub-section). The different curves correspond to models with different values of the o/p ratio of ${\rm H}_{2}{\rm D}^{+}$ (o/p- ${\rm H}_{2}{\rm D}^{+}$), from 0.03 (bottom curves) to 2.0 (top curves). As discussed by Flower et al. (2004), the o/p ratio is a sensitive function of the o/p H2 ratio and, ultimately, of the gas temperature (see their Fig. 6) and at $T_{\rm gas}$ < $\sim$15 K, it changes from $\simeq$0.03 to values probably larger than one (this last statement is valid if the curve in Fig. 6 of Flower et al. (2004) is simply extrapolated at temperatures lower than 9 K, the minimum value in the figure[*]). In all curves, the t5 models show a slight decrease in both the ${\rm H}_{2}{\rm D}^{+}$ column density and in the deuterium fractionation. In fact, at such high central densities: (i)  ${\rm H}_{2}{\rm D}^{+}$ is efficiently converted into ${\rm D}_{2}{\rm H}^{+}$ and D3+, thus lowering the total ${\rm H}_{2}{\rm D}^{+}$ column along the line of sight; (ii) ${\rm N}_2$ significantly freezes onto dust grains, so that the ${\rm N}_2{\rm D}^+$/ ${\rm N}_2{\rm H}^+$ column density ratio - our measure of the observed $R_{\rm D}$ - traces regions away from the center, where the density is lower and the temperature is higher. If a protostar is present, the ${\rm H}_{2}{\rm D}^{+}$ column density increases again because of the less efficient transformation of ${\rm H}_{3}^{+}$ into ${\rm D}_{3}^{+}$, while $R_{\rm D}$ decreases (see also Fig. 3) because of the less abundant ${\rm D}_{3}^{+}$. From the comparison between the data and the model, one is tempted to conclude that indeed variations in the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio (and ultimately in the gas temperature of the central few thousand AU of dense cores) can explain the observed scatter among the cores. It is interesting to note that the two pre-stellar cores with high values of $R_{\rm D}$ and relatively low ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities are both embedded in the Ophiuchus Molecular Cloud Complex: Oph D and 16293E. Chemical abundances observed in these two cores are consistent with a lower ${\rm H}_{2}{\rm D}^{+}$ o/p ratio, suggesting possible (slightly) larger kinetic temperatures.

Figure 5 shows two other attempts to interpret the data. The top panel considers exactly the same models as in Fig. 4 but now the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio is fixed at 0.5 (consistent with dense gas at 9 K; see Flower et al. 2004), whereas the free parameter is the cosmic-ray ionization rate $\zeta $, which is varied from 6 $\times $ 10-18 s-1(bottom dashed and solid curves) and 3 $\times $ 10-17 s-1 (top curves). We note that variations in the cosmic-ray ionization rate are known to exist in the Galaxy (van der Tak et al. 2006). The four pre-stellar cores with the largest N(ortho- ${\rm H}_{2}{\rm D}^{+}$) values (L 492, L 1544, L 694-2, and L 183) can all be reproduced by t2-t4 (t1-t2) models with the GHR02 (standard) rates and $\zeta $ = 1 $\times $ 10-17 s-1 (with L 429 (L 183) being the most (least) dynamically evolved). Significantly lower values of $\zeta $ (<6 $\times $ 10-18 s-1) appear to be required in the protostellar and Ophiuchus cores. However, an alternative way to reproduce these data, without requiring extremely low $\zeta $ values, is to lower the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio, as found before (the thick dashed curves represent models with o/p- ${\rm H}_{2}{\rm D}^{+}$ = 0.3 and $\zeta $ = 6 $\times $ 10-18 s-1). In B 68, only the low o/p- ${\rm H}_{2}{\rm D}^{+}$ model curve (at early evolutionary times) can reproduce the low ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities and deuterium fractionations, the lowest in the sample.

In the bottom panel of Fig. 5, we consider five clouds with structure as in the ti (i = 2, ..., 5) models (see Eq. (17)) and for each of them we follow the chemical evolution, checking the results after 103, 104, 105, and 106 yr. This allows us to explore how different chemical ages can change the gas composition and avoid the problem of fixing the (unknown) chemical evolution time as in the bottom panel of Fig. 4. As for the top panel, the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio has been fixed at 0.5, except for the thick curves representing the $t_{\rm 5a}$ protostellar cloud models, where o/p = 0.1, assuming that the gas has been (slightly) heated compared to the pre-stellar cores. The cosmic-ray ionization rate is fixed at 1.3 $\times $ 10-17 s-1. The four ortho- ${\rm H}_{2}{\rm D}^{+}$ - rich pre-stellar cores can be reproduced by t3 models, with (chemical) ages between 104 and 106 yr, when GHR02 rate coefficients are used. On the other hand, the protostellar cores and the two Ophiuchus cores are better matched by the more dynamically evolved (centrally concentrated) $t_{\rm 5a}$ and t5 models, respectively, after about 103-104 years of (chemical) evolution.

5.2.2 N(ortho-H2D+) vs. the observed $f_{\rm D}$

From the models described in the previous subsection, one can directly derive the abundance of CO within each cloud model as a function of cloud radius and, as before, $N({\rm CO})$is obtained, after integrating the CO number density along the line of sight and smoothing the results with a 22 $^{\prime\prime}$ beam width (to simulate observations carried out at the IRAM 30 m antenna, where most of the cores have been observed). From the model column density, the CO depletion factor, $f_{\rm D}$(CO) is easily calculated using the expression:

$\displaystyle %
f_{\rm D}({\rm CO}) = \frac{x_{\rm can}({\rm CO})}
{N({\rm CO})/N({\rm H}_2)},$     (18)

where $x_{\rm can}$(CO) is the ``canonical'' or ``undepleted'' abundance of CO (assumed here equal to 9.5 $\times $ 10-5, following Frerking et al. 1982, but known to be uncertain by a factor of about 2; see e.g. Lacy et al. 1994).

To compare our model predictions with the data, we collect from the literature the values of observed $f_{\rm D}$(CO) and plotted them versus N(ortho- ${\rm H}_{2}{\rm D}^{+}$) in Fig. 6. The majority of the $f_{\rm D}$(CO) data comes from Crapsi et al. (2005), except for: (i) TMC-1C ($f_{\rm D}$(CO) = 6.9, from Schnee et al. 2007a); (ii) L 183 ($f_{\rm D}$ = 5, from Pagani et al. 2005); (iii) NGC 1333- ${\rm DCO}^+$ ($f_{\rm D}$ = 12, from Jørgensen et al. 2002); (iv) B 1 ($f_{\rm D}$(CO) = 3.2, from Lis et al. 2002b); (v) IRAM 04191 ($f_{\rm D}$(CO) = 3.5 from Belloche & André 2004); (vi) OriB 9 ($f_{\rm D}$(CO) = 3.6, from Harju et al. 2006, for N(${\rm H}_2$); and Caselli & Myers 1995, for $N({\rm CO})$). Also, for L 1544, L 429, L 694-2 and Oph D, the new values of N(${\rm H}_2$), adopted to take into account the temperature structure (see Table 3), imply different values of $f_{\rm D}$ (larger by a factor of about 1.4, the ratio between the new and old N(${\rm H}_2$) values, as explained in Sect. 4.1) from those reported by Crapsi et al. (2005). The data and model results are shown in Fig. 6. In general, the presence of embedded young stellar objects appears to lower the ${\rm H}_{2}{\rm D}^{+}$ column density, without much affecting the amount of CO freeze-out, which is probably still large in the high density and cold protostellar envelopes. The possible (small) temperature increase caused by the central heating is thus not sufficient to significantly release CO back into the gas phase, but it can affect the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio (and, consequently, the ortho- ${\rm H}_{2}{\rm D}^{+}$ column density), as discussed in the previous sub-section.

\end{figure} Figure 6: N(ortho- ${\rm H}_{2}{\rm D}^{+}$) vs. the observed CO depletion factor, $f_{\rm D}$(CO). The same models ti (i=1, ..., 5) described for the previous figures are used here, varying the cosmic-ray ionization rate $\zeta $ ( top panel) and the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio ( bottom panel) within the range of Fig. 4 and 5. The thick dashed curve refers to (GHR02) models with o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio of 0.1 and $\zeta $ = 6 $\times $ 10-18 s-1. As found before, variations in the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio, plus differences in the physical structure, provide a way to reproduce the whole observed spread in the data. The labels in the top panel correspond to the source names: A $\equiv $ L 1498, B $\equiv $ TMC-2, C $\equiv $ TMC-1C, D $\equiv $ L 1517B, E $\equiv $ L 1544, F $\equiv $ L 183, G $\equiv $ Oph D, H $\equiv $ B 68, I $\equiv $ L 429, J $\equiv $ L 694-2, K $\equiv $ NGC 1333- ${\rm DCO}^+$, L $\equiv $ B 1, M $\equiv $ IRAM 04191, N $\equiv $ L 1521F, O $\equiv $ Ori B9.
Open with DEXTER

In the top panel of Fig. 6, the data are compared with the same models described in Fig. 5 (left panel; ti, i = 1, 2, ...5), where the cosmic-ray ionization rate is varied from 6 $\times $ 10-18 s-1 to 3 $\times $ 10-17 s-1. The o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio has been fixed at 0.5, except for the thick dashed curve, where o/p = 0.1, a value probably more appropriate for protostellar cores (see previous sections). Now, t2-t3 models (depending on the rate coefficient values adopted for the proton-deuteron exchange reactions), with $\zeta $ $\simeq$ 10-17 s-1, can explain the observations toward the ortho- ${\rm H}_{2}{\rm D}^{+}$ - rich pre-stellar cores, with the exception of L 183, where the low $f_{\rm D}$ value suggests a younger dynamical phase, consistent with what is found in the previous sub-section (we also note that L 183 and TMC-1C appear to have similar ages, but with different cosmic ray ionization rates and/or o/p- ${\rm H}_{2}{\rm D}^{+}$ ratios, see bottom panel of Fig. 6). Lower o/p- ${\rm H}_{2}{\rm D}^{+}$ ratios are needed to reproduce the protostellar and Ophiuchus cores, as found for the N(ortho- ${\rm H}_{2}{\rm D}^{+}$) vs. $R_{\rm D}$ relation.

The bottom panel of Fig. 6 shows the same set of data and models, but now $\zeta $ has been fixed at 1.3 $\times $ 10-17 s-1, whereas the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio has been changed as in Fig. 4. As already noted, L 1544, L 429, L 694-2, and L 183 data are best reproduced by large values of the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio, and the appropriate physical structure is that of t1-t3 models, similarly to the top panel (with L 183 being the least evolved). This implies cores slightly more evolved than what is found in Figs. 4 (bottom panel) and 5 (top panel), where (standard rate) models between t1 and t2 were preferred. The small discrepancy can be understood if the predicted CO depletion factor is too large compared to observations. Reasons for this could be: (i) we are missing an important desorption mechanism (besides the cosmic-ray impulsive heating and the thermal desorption, the latter not being efficient at the temperatures of these objects); (ii) unlike our spherical and isolated model cores, real cloud cores are embedded in molecular clouds where CO is quite abundant. Thus, the observed ``extra'' gaseous CO may be part of the undepleted material accreting onto the core from the surrounding molecular cloud (see also Schnee et al. 2007a, for a similar conclusion in the particular case of TMC-1C).

6 Conclusions

Low-mass dense cloud cores have been observed with the CSO antenna at the frequency of the ortho- ${\rm H}_{2}{\rm D}^{+}$(11,0-11,1) line. The main conclusions of this work are:

In starless cores, the line has been detected in 7 (out of 10) objects. The brightest lines ( $T_{\rm mb}$ = 0.7-0.9 K) are observed toward the densest and more dynamically evolved starless cores (L 1544, L 183, Oph D, L 429, and L 694-2), where the deuterium fractionations and the CO depletion factors are largest. Non-detections are found in cores less centrally concentrated and dense than the rest of the sample (L 1498, TMC-2, and B 68). In B 68, recent APEX observations have detected a faint ortho- ${\rm H}_{2}{\rm D}^{+}$ line with intensities consistent with our upper limit.

Significant differences in line widths are also observed among starless cores, with the narrowest lines ( $\Delta {v}$ $\simeq$ 0.4  ${\rm km~s}^{-1}$) found in TMC-1C, L 1517B, Oph D, and L 183, whereas L 1544, L 429, and L 694-2 show $\Delta {v}$ = 0.5-0.7  ${\rm km~s}^{-1}$.

The ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) line has been detected in 4 out of 6 protostellar cores, where we find $T_{\rm mb}$ $\leq$ 0.5 K (the largest value observed is toward L 1521F, which hosts a very low luminosity object recently detected by Spitzer). The broadest lines are observed toward the two cores in Perseus NGC 1333- ${\rm DCO}^+$ and B 1, where $\Delta {v}$ $\simeq$ ${\rm km~s}^{-1}$. The ortho- ${\rm H}_{2}{\rm D}^{+}$ (11,0-11,1) lines have broader non-thermal widths than ${\rm N}_2{\rm D}^+$(2-1) lines, even in L 1544, where the two transitions have similar extension and morphology.

The ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities range between 2 and 40 $\times $ 1012 cm-2 in starless cores and between 2 and 9 $\times $ 1012 cm-2 in protostellar cores. Thus, protostars in the earliest stages of their evolution appear to have already changed the chemical structure of their envelopes by lowering the ortho- ${\rm H}_{2}{\rm D}^{+}$ fractional abundance but not their deuterium fractionation. This is probably due to (a few degrees) variation of the gas temperature, which strongly affects the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio. A similar effect is also found in the two Ophiuchus cores, suggesting a (slightly) higher gas temperature than in the other (mainly Taurus and Perseus) cores.

Simple models suggest that variations in the gas kinetic temperature, CO depletion factor, grain size distribution, cosmic-ray ionization rate and volume densities can largely affect the ${\rm H}_{2}{\rm D}^{+}$ abundances relative to ${\rm H}_{3}^{+}$. In particular, gas temperatures above 15 K, low CO depletion factors and large abundance of negatively charged small dust grains or PAHs drastically reduce the deuterium fractionation to values inconsistent with those observed toward pre-stellar and protostellar cores.

Plots of the ortho- ${\rm H}_{2}{\rm D}^{+}$ column density vs. (i) the deuterium fractionation observed in species such as ${\rm N}_2{\rm H}^+$ and ${\rm NH}_{3}$; and (ii) the observed CO depletion factor, show a large scatter. The data can be reproduced by chemical models of centrally concentrated cores assuming variations in the o/p ratio of ${\rm H}_{2}{\rm D}^{+}$ (ultimately linked to kinetic temperature variations). Changes in the cosmic-ray ionization rate between 6 and 30 $\times $ 10-18 s-1 can also explain part of the scatter, but not those objects with large deuterium fractionations and low ortho- ${\rm H}_{2}{\rm D}^{+}$ column densities (such as the Ophiuchus pre-stellar cores and the protostellar ones), for which a decrease of the o/p- ${\rm H}_{2}{\rm D}^{+}$ ratio is required. The most deuterated and ${\rm H}_{2}{\rm D}^{+}$-rich objects are better reproduced by chemical models of centrally concentrated cores with $n_{\rm c}$ = a few $\times $ 106  ${{\rm cm}}^{-3}$ and chemical ages between 104 and 106 yr. Protostellar cores, plus the two Ophiuchus cores, are better matched by lower o/p- ${\rm H}_{2}{\rm D}^{+}$ ratios, more dynamically evolved (central densities $\ge$107 cm-3) models and chemical ages of $\simeq$103-104 years. To better constrain dynamical and chemical ages, the rate coefficient for the proton-deuteron exchange reaction needs to be well defined.

The para- ${\rm H}_{3}{\rm O}^{+}$(11--21+) upper limits are consistent with radiative transfer calculations if the fractional abundance of ${\rm H}_{3}{\rm O}^{+}$ is $\la$10-8. The para- ${\rm D}_{2}{\rm H}^{+}$(11,0-10,1) upper limits provide an upper limit of the para- ${\rm D}_{2}{\rm H}^{+}$ to ortho- ${\rm H}_{2}{\rm D}^{+}$ column density ratio ($\equiv $p/o). We find p/o < 1 in all the sources (except for L 183 and NGC 1333- ${\rm DCO}^+$, where p/o < 3), consistent with chemical model predictions of high density (2 $\times $ 106  ${{\rm cm}}^{-3}$) and low temperature ( $T_{\rm kin} < 10$ K) clouds (Flower et al. 2004).
More accurate determinations of temperature and density profiles, as well as observations of the para- ${\rm H}_{2}{\rm D}^{+}$(10,1-10,0) line at 1.37 THz, are sorely needed to place more stringent constraints on gas-grain chemical processes in dense cloud cores.

We thank the anonymous referee for his/her very detailed review, which greatly improved the paper. P.C. acknowledges support from the Italian Ministry of Research and University within a PRIN project. Part of this research has been supported by NSF grant AST 05-40882 to the CSO. The authors thank the staff of the CSO telescope for their support. We also thank E. Hugo and S. Schlemmer for providing their collisional rates prior to publication.

Appendix A: The search of para-H3O+

Table A.1: p- ${\rm H}_{3}{\rm O}^{+}$ column density upper limits.

The ${\rm H}_{3}{\rm O}^{+}$ line was observed to investigate the chemistry of oxygen in dense cores and, with the help of chemical models, to place some constraints on the oxygen abundance, which, analogously to CO, significantly affects the deuterium fractionation.

We searched for para- ${\rm H}_{3}{\rm O}^{+}$(1-1-2+1) in seven dense cores but only upper limits were measured. Table A.1 lists the results of this search, including the rms noise and the corresponding upper limits of the column density, which have been calculated in two different ways: (i) using the volume density and kinetic temperature values listed in Table 3 (N1, Col. 6); and (ii) assuming a volume density of 105  ${{\rm cm}}^{-3}$ and $T_{\rm kin}$ = 10 K (N2, Col. 8). In both cases, the RADEX[*] code has been used (van der Tak et al. 2007). The observed upper limits (see Table A.1) are compatible with para- ${\rm H}_{3}{\rm O}^{+}$ fractional abundance upper limits of 10-8 according to calculation made using the Ratran code (Hogerheijde & van der Tak 2000).



Copyright ESO 2008