© ESO, 2010

1. Introduction

Although deuterium is a trace element (D/H  ~ 10^-5), the study of deuterium fractionation in molecules is a very vivid field of research. The main astrochemical interest derives from the understanding of fractionation processes ultimately being able to help to understand the formation of the related hydrogenated species. As an example, the first detection of a doubly-deuterated molecule, D₂CO, by Turner (1990) was interpreted as the hint that grain chemistry may be at work for synthesizing formaldehyde. However, the impact of deuterium fractionation studies rapidly outpassed the mere astrochemical interest. Since deuterated molecules are formed primarily at low temperatures and high densities, they became valuable probes for studying, e.g., the kinematics of cold cloud cores harboring the earliest stages of star formation (van der Tak et al. 2005).

The basic process transferring deuterium from the HD molecule to other molecules is through ion-molecule reactions involving ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, ${\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}$, and C₂${\mathrm{H}}_{\mathrm{2}}^{\mathrm{+}}$ (see Parise et al. 2009, and references therein). The route involving ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ is the dominant one at very low temperatures (T < 20 K), whereas the importance of ${\mathrm{CH}}_{\mathrm{3}}^{\mathrm{+}}$ takes over at slightly higher temperatures (20 < T < 40 K) owing to its higher endothermicity (Turner 2001; Roueff et al. 2007; Parise et al. 2009). H₂D⁺ is thus a very good probe of very cold objects, and has been extensively used in recent years to probe gas in which CO is strongly depleted, i.e., frozen out on dust grain surfaces (e.g. Caselli et al. 2003). The deuteration level is also a decisive argument in the study of the evolutionary state of early phases of star formation. The deuterium fractionation is believed to decrease once the central protostar starts to heat its envelope (e.g., Emprechtinger et al. 2009). For instance, partly based on molecular deuteration, Belloche et al. (2006) conclude that the dense core Cha-MMS1 in the Chamaeleon I molecular cloud is likely at the stage of the first hydrostatic core when H₂ has not been dissociated yet.

The study of several isotopologues of methanol toward the low-mass protostar IRAS16293-2422 (Parise et al. 2002, 2004) showed that grain chemistry models require a very high atomic D/H ratio accreting on the grains to explain the high CHD₂OH and CD₃OH abundances. At the time, gas-phase models could not reproduce such a high atomic D/H ratio, but a decisive step forward was made by including D₂H⁺ and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$ in the reaction networks (Roberts et al. 2003). Since then, D₂H⁺ has been searched for by observing of its ground para state line. The frequency of this line was initially published as 691.660440 GHz (± 19 kHz (1σ), Hirao & Amano 2003) and was later revised to 691.660483 GHz (±20 kHz (1σ), Amano & Hirao 2005). These astronomical searches are extremely difficult because of the high frequency of the molecular line, for which excellent instruments and weather conditions are required. Using the Caltech Submillimeter Observatory 10.4 m telescope (CSO), Vastel et al. (2004) detected a  ~3.3σ line (in peak) toward the Ophiuchus dense core 16293E (4.4σ in integrated intensity). This detection is the only published one to date, and its low signal-to-noise ratio, as well as the velocity shift between H₂D⁺ and D₂H⁺ (0.23 km s^-1 after remeasurement of o-H₂D⁺ and p-D₂H⁺ frequencies, Amano & Hirao 2005), certainly require a confirmation of the detection by improving the observation sensitivity and frequency calibration precision (which was  ~0.1 km s^-1 in the observations of Vastel et al. 2004). We present here the first secure detection of D₂H⁺ toward the prestellar core H-MM1 in L1688, the main molecular cloud in the Ophiuchus star-forming region.

2. Observations

2.1. The source

Fig. 1 The L1688 region, as mapped by SCUBA (850 μm), from the SCUBA legacy catalogs (see text). The contours are in steps of 3σ (=0.18 Jy/beam). The beamsize is 14′′ (HPBW).

The L1688 cloud (the main cloud in the Ophiuchus dark cloud) has been mapped both in molecular tracers (e.g. Loren et al. 1990) and in continuum emission (e.g. at 1.3 mm by Motte et al. 1998). These early mapping studies revealed different areas of high column density, named Oph-A to Oph-G. Using the Submillimeter Common User Bolometer Array (SCUBA) at the 15 m James-Clerck-Maxwell Telescope (JCMT), Johnstone et al. (2004) covered a larger region, leading to detection of two new dense cores, which they named Oph-H and Oph-I, following the same nomenclature. Oph-H contains a single continuum peak, H-MM1. The distance of the Ophiuchus region has recently been accurately determined to 120 pc (Loinard et al. 2008).

We present in Fig. 1 a part of the SCUBA map of the Ophiuchus region, as available from the SCUBA Legacy catalogs¹ (data: scuba_F_353d1_16d8_850um.emi.fits, Di Francesco et al. 2008), showing the relative position of all dense cores. In this paper, we present observations toward the H-MM1 core, at coordinates α₂₀₀₀ = 16^h27^m58.3^s and δ₂₀₀₀ = −24°33′42.0′′.

H-MM1 corresponds to the Bolo25 source of the 1.1 mm BOLOCAM (CSO) map observed by Young et al. (2006). The FWHM beam size of these observations is 31′′. Assuming a temperature of 10 K and a dust opacity²κ_ν = 0.0114 cm² g^-1 at 1.12 mm, these authors infer a total mass of 0.24 M_⊙ in 40′′ and measure a source FWHM size of 47′′  ×   56′′. They measure a peak visual extinction, A_V, of 23 mag, and a mean density of 5 × 10⁵ cm^-3.

Khanzadyan et al. (2004) and Stanke et al. (2006) carried out an unbiased search for protostellar signatures in the ρ-Oph cloud. The closest source they identify in their continuum maps is MMS38 (at 9′′ from our position), which they classify as starless. No optical/NIR counterpart (2MASS) or sign of outflow is found in this source.

The region was also mapped with the Spitzer Space telescope, in the frame of the Cores to Disks legacy project (c2d). We downloaded the corresponding data from the Spitzer website³. No source is detected in any of the Spitzer bands, especially at 70 μm, which is a good tracer of embedded protostars (see Dunham et al. 2008). H-MM1 thus appears to be a dense prestellar core.

2.2. APEX line observations

Using the APEX2a single-pixel and CHAMP⁺ array receivers at the Atacama Pathfinder Experiment telescope (APEX), we observed H₂D⁺ (372 GHz) and D₂H⁺ (692 GHz) toward H-MM1 at the position given in Sect. 2.1. In each observation we used the position-switch mode, with the off-position at (0′′, 300′′) equatorial offsets. Calibration was done by measuring regularly the sky brightness and cold and hot loads in the cabin.

2.2.1. H₂D⁺

Using the APEX2a receiver (Risacher et al. 2006), we observed the ground-state o-H₂D⁺ line toward H-MM1. The frequency of this line is 372.421385 GHz (Amano & Hirao 2005). Observations took place on June 29, 2007, under very good weather conditions (precipitable water vapor PWV = 0.24 mm). The system temperature was 250 K. The backend, a Fast Fourier Transform Spectrometer (FFTS), was used with 8192 channels, yielding a 0.0983 km s^-1 channel separation. The effective frequency resolution of the FFTS in this configuration is 159 kHz (Klein et al. 2006, and B. Klein, priv. comm.), i.e. 0.13 km s^-1 at 372 GHz.

The beam efficiency of APEX at 372 GHz was measured as 0.73 for compact sources (as measured on Mars and Jupiter, which were respectively 8′′ and 32−35′′ at that time, Güsten et al. 2006). As the D₂H⁺ emission seems extended on a scale comparable to the footprint of CHAMP⁺ (~40′′, see Fig. 4), H₂D⁺ is likely to also be extended on that scale. We therefore adopt the T_mb scale for our analysis, using the beam efficiency of 0.73. The forward efficiency of the telescope is 0.95 at this frequency.

We reached an rms noise level of 72 mK (${\mathit{T}}_{\mathit{a}}^{\mathrm{\ast }}$ scale) at the initial resolution (0.13 km s^-1).

2.2.2. D₂H⁺

Using the CHAMP⁺ array receiver on the APEX telescope, we targeted the prestellar core H-MM1. The observations took place in June and August 2009, under good to very good weather conditions (most of the time with PWV < 0.45 mm). The typical system temperature in the central pixel of the array varied between 800 and 1400 K. Figure 2 presents the overlay of the CHAMP⁺ footprint on the 850 μm continuum emission. Each CHAMP⁺ pixel was tuned to the CO(6-5) frequency (691.4730763 GHz), which is separated by 187.4 MHz from the D₂H⁺ frequency. We adopt here the most recent measurement of Amano & Hirao (2005) of the D₂H⁺ line frequency: 691.660483 GHz. The beam size is 9′′. The CHAMP⁺ AFFTS was used with a channel separation of 183 kHz (0.079 km s^-1). The effective resolution is 212 kHz, i.e. 0.092 km s^-1.

The beam efficiency at 692 GHz was measured as 0.52 on Jupiter (angular size 47.3′′ at that time). Because the D₂H⁺ emission is observed to be extended (see below) on a typical size of 40′′, we adopt the beam efficiency of 0.52. The forward efficiency is 0.95 at this frequency.

We reached an rms noise level of 34 mK (${\mathit{T}}_{\mathit{a}}^{\mathrm{\ast }}$ scale) at the initial resolution (0.09 km s^-1).

2.3. Continuum data

We retrieved the 850 μm map of the Ophiuchus region from the SCUBA Legacy catalogs (data: scuba_F_353d1_16d8_850um.emi.fits, Di Francesco et al. 2008). The region covering the L1688 cloud is shown in Fig. 1, and was discussed by Jørgensen et al. (2008). The beam of SCUBA at this frequency is 14′′. The local rms in the H-MM1 region of the map is found to be 60 mJy/beam.

Fig. 2 850 μm continuum emission of the H-MM1 core. The map is centered on the position of H-MM1 (see text). The contours are in steps of 3σ (=0.18 Jy/beam). The blue numbers show the different pixels of the CHAMP+ low-frequency array. The blue circles show the D2H+ beams. In red, the H2D+ beam.

Fig. 3 H2D+ (upper panel) and D2H+ (lower panel) toward the central position of H-MM1.

3. Results

The main result of this work is the first unambiguous detection of p-D₂H⁺ toward a prestellar core. Figure 3 presents the H₂D⁺ and D₂H⁺ emission at the (0′′, 0′′) position. The agreement of the LSR velocities, v_lsr, for both lines is remarkable (see Table 1) and further strengthens the assignment of the lines. Figure 4 presents the D₂H⁺ observation on each of the CHAMP⁺ pixels. The noise level is different from pixel to pixel, owing to the different performances of the off-axis pixels. The line is clearly detected in at least two adjacent pixels (pixels 2 and 4). There is also a tentative detection in pixel 3, with lower intensity. The D₂H⁺ emission is thus extended on the scale of the CHAMP⁺ footprint (~40′′), corresponding to a physical size of  ~4800 AU at the distance of 120 pc. The observational results are listed in Table 1.

Detection of p-D₂H⁺ in H-MM1 does not imply that it is a peculiar core. None of the searches for p-D₂H⁺ published in Caselli et al. (2008) towards five prestellar cores have reached upper limits that are deep enough to discard the detection of the line with the same peak temperature in those cores.

3.1. Linewidths

The measured linewidths are listed in Table 1. Because the lines are not very broad compared to the frequency resolution of the FFTS, we also computed deconvolved linewidths to correct for instrumental broadening (Col. 7).

The expected thermal linewidth for a molecule of mass m is $\begin{array}{ccc}\mathrm{\Delta }{\mathrm{v}}_{\mathrm{therm}\mathit{,}\mathrm{FWHM}}\mathrm{=}\mathrm{2}\mathrm{\times }\sqrt{\mathrm{2}\hspace{0.17em}\ln \mathrm{2}}\mathrm{\times }\sqrt{\frac{\mathit{kT}}{\mathit{m}}}\mathrm{·}& & \end{array}$For H₂D⁺, $\mathrm{\Delta }{\mathrm{v}}_{\mathrm{therm}}\mathrm{=}\mathrm{0.107}\mathrm{\times }\sqrt{\mathit{T}}$ km s^-1. For D₂H⁺, $\mathrm{\Delta }{\mathrm{v}}_{\mathrm{therm}}\mathrm{=}\mathrm{0.0955}\mathrm{\times }\sqrt{\mathit{T}}$ km s^-1, with T in Kelvin.

The small linewidth of D₂H⁺ on the (0′′, 0′′) position constrains the kinetic temperature to T < 13 K. For this upper limit value of the temperature, the rather large H₂D⁺ linewidth points to nonthermal velocity broadening of 0.3 km s^-1. The different nonthermal velocity contributions to both lines might be a sign that the two molecules do not trace exactly the same gas. The larger nonthermal velocity broadening of o-H₂D⁺ may imply that o-H₂D⁺ is present in a more diffuse and turbulent gas than p-D₂H⁺, as is also supported by chemical models, which predict that the p-D₂H⁺/ o-H₂D⁺ ratio raises with density (see Fig. 9). In this case, the p-D₂H⁺/ o-H₂D⁺ ratio that we derive in the densest gas will likely be a lower limit.

3.2. Correlation with continuum emission

The continuum emission at the position of the individual pixels of CHAMP⁺ is listed in Table 1. The positions where p-D₂H⁺ is detected correspond to the positions where the continuum emission is the strongest. The lower tentative detection in pixel 3 is also consistent with the continuum being weaker at this position. The correlation between the p-D₂H⁺ and the continuum emission is quite good. This implies that the p-D₂H⁺ emission comes from the densest regions of the core, in agreement with the narrow width of the line.

3.3. Column densities

3.3.1. LTE assumption

The critical densities for the o-H₂D⁺ and p-D₂H⁺ transitions are 1.3 × 10⁵ and 5.6 × 10⁵ cm^-3 respectively (Hugo et al. 2009). Although the density of the dense core may be close to or only marginally higher than these critical densities, we start with an LTE analysis as a first approximation. We consider a non-LTE approach in the following section.

We consider the ortho and para species of each molecule as a different species, which do not interconvert via collisions with H₂. We thus compute separately the partition function of ortho and para species from the transitions given in the Cologne Database of Molecular Spectroscopy (CDMS, Müller et al. 2005). The partition functions are computed as $\begin{array}{ccc}\mathit{Q}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\sum _{\mathrm{levels}}{\mathit{g}}_{\mathit{i}}{\mathit{g}}_{\mathit{r}}{\mathrm{e}}^{\mathrm{(}\mathrm{-}\mathit{E}\mathit{/}\mathit{kT}\mathrm{)}}& & \end{array}$where g_i is the spin degeneracy (1 and 3 for para- and ortho- H₂D⁺ respectively, and 1 and 2 for para- and ortho- D₂H⁺) and g_r the rotational degeneracy (g_r = 2J + 1). The energies of the levels were all expressed relative to the ground level 0₀₀ of each molecule. These partition functions are presented in Table 2.

Fig. 4 D2H+ observations. Each panel shows a pixel of CHAMP+ (as displayed in Fig. 2). The vertical line shows the vlsr of the source as measured with the o-H2D+ line (4.23 km s-1). The y-axis is in Kelvin on the ${\mathit{T}}_{\mathit{a}}^{\mathrm{\ast }}$ scale.

The column density in the upper level of the transition can be directly derived from the measure of the opacity $\begin{array}{ccc}{\mathit{N}}_{\mathrm{u}}\mathrm{=}\frac{\mathrm{1}}{\mathit{h}{\mathit{B}}_{\mathrm{ul}}}\frac{\mathrm{1}}{{\mathrm{e}}^{\mathit{h\nu }\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}}}\mathrm{-}\mathrm{1}}{\mathit{\tau }}_{\mathrm{peak}}\mathrm{\Delta v}\mathit{.}& & \end{array}$The opacity of the transition, τ, is given by $\begin{array}{ccc}{\mathit{T}}_{\mathrm{mb}}\mathrm{=}\mathrm{\left[}\mathit{J}\mathrm{\right(}{\mathit{T}}_{\mathrm{ex}}\mathrm{)}\mathrm{-}\mathit{J}\mathrm{(}{\mathit{T}}_{\mathrm{bg}}\mathrm{\left)}\mathrm{\right]}\mathrm{[}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{]}& & \end{array}$(1)where $\mathit{J}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\frac{\mathit{h\nu }\mathit{/}\mathit{k}}{\exp \mathrm{(}\mathit{h\nu }\mathit{/}\mathit{kT}\mathrm{)}\mathrm{-}\mathrm{1}}\mathrm{·}$

Assuming LTE (i.e. T_ex = T_kin), we can then derive the total column density of the species $\begin{array}{ccc}{\mathit{N}}_{\mathrm{tot}}\mathrm{=}\frac{\mathrm{8}\mathit{\pi }{\mathit{\nu }}^{\mathrm{3}}}{{\mathit{c}}^{\mathrm{3}}}\frac{\mathit{Q}\mathrm{\left(}{\mathit{T}}_{\mathrm{ex}}\mathrm{\right)}}{{\mathit{g}}_{\mathrm{u}}{\mathit{A}}_{\mathrm{ul}}}\frac{{\mathrm{e}}^{{\mathit{E}}_{\mathit{u}}\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}}}}{{\mathrm{e}}^{\mathit{h\nu }\mathit{/}\mathit{k}{\mathit{T}}_{\mathrm{ex}}}\mathrm{-}\mathrm{1}} {\mathit{\tau }}_{\mathrm{peak}}\mathrm{\Delta v}\mathit{.}& & \end{array}$We used the A_ul of Ramanlal & Tennyson (2004). For both transitions, they are 11% different from the value quoted in the CDMS database. We used the accurate frequencies derived by Amano & Hirao (2005).

We note that Eq. (1) implies that T_mb ≤ J(T_ex). This condition for the 692 GHz line of D₂H⁺ translates into J(T_ex) ≥ 0.25 K, i.e. T_ex ≥ 7 K. This condition, as well as the previous linewidth argument, allow restricting the kinetic temperature to the range 7–13 K.

The results are presented in Table 3. Depending on the assumed temperature, the p-D₂H⁺/ o-H₂D⁺ ratio varies in the range 0.86–4.4.

3.3.2. Out-of-equilibrium radiative transfer

Thanks to the recent computation of all elastic and nonelastic collision rate coefficients by Hugo et al. (2009), we can derive the column density using an out-of-equilibrium radiative transfer model. We used the non-LTE radiative transfer code RADEX (van der Tak et al. 2007), in the isothermal sphere geometry, with the escape probability as a function of the opacity as follows: $\begin{array}{ccc}\mathit{\beta }\mathrm{=}\frac{\mathrm{1.5}}{\mathit{\tau }}\left[\mathrm{1}\mathrm{-}\frac{\mathrm{2}}{{\mathit{\tau }}^{\mathrm{2}}}\mathrm{+}\left(\frac{\mathrm{2}}{\mathit{\tau }}\mathrm{+}\frac{\mathrm{2}}{{\mathit{\tau }}^{\mathrm{2}}}\right){\mathit{e}}^{\mathrm{-}\mathit{\tau }}\right]\mathrm{·}& & \end{array}$Because the rates for collisions with o-H₂ and p-H₂ are quite different, the nature of the collider matters for the excitation. Our chemical models (presented in Sect. 4) predict that the steady-state ortho/para-H₂ is within a factor of two around 10^-4 for all physical conditions probed (Fig. 8). We checked that the excitation is unchanged (within numerical precision) when considering only para-H₂ or para-H₂ with traces of ortho-H₂ (o/p-H₂  ≤ 2 × 10^-4). In a first approximation we thus assumed collisions with para-H₂ only in the derivation of the column density of o-H₂D⁺and p-D₂H⁺. Pagani et al. (2009), however, showed that the o/p-H₂ ratio is unlikely in steady-state in prestellar cores (see also Sect. 4.2 and Fig. 7). For completeness, we checked the computed column densities in the case where the o/p-H₂ ratio is higher than the steady-state value. We find that the derived column densities start to deviate noticeably from the ones derived with only p-H₂ for values of the o/p-H₂ ratio on the order of unity.

The results with excitation with para-H₂ only are presented in Fig. 5a. The results with an o/p-H₂ ratio of 10^-2 are found to differ from Fig. 5a by less than 1%. The column densities derived within the LTE approximation are also shown on the figure. Even at a density of 10⁶ cm^-3, the LTE assumption is not accurate, especially for p-D₂H⁺, which has the highest critical density. One can indeed check that the excitation temperature does not reach the kinetic temperature in any of these cases (Fig. 5b).

Fig. 5 a) Column densities at the (0, 0) position, computed from the nonequilibrium radiative transfer assuming the beam filling factor for both molecules is 1, as a function of the physical conditions (kinetic temperature, density). The thick black lines show the column densities computed in the LTE assumption. b) Excitation temperature in the same conditions and color conventions as panel a). The black line represents Tex = Tkin.

Assuming o/p-H₂ = 3 (the highest possible ratio at thermodynamical equilibrium) results in column densities about 10–20% higher than in the case of excitation with only para-H₂, the effect being stronger at low densities, as expected (at high densities, the system is close to LTE and does not depend on the nature of the collider). The abundance ratio of the two molecules is even less sensitive to the o/p-H₂ ratio, as shown in Fig. 6. This figure presents the p-D₂H⁺/ o-H₂D⁺ ratio resulting from our radiative transfer analysis, for different densities and in the two cases o/p-H₂ = 0 and 3. In the following, we aim at investigating the implications of the p-D₂H⁺/ o-H₂D⁺ ratio with the help of a state-of-the art chemical model.

Fig. 6 p-D2H+/ o-H2D+abundance ratio resulting from the radiative transfer modeling of the observations, for different densities and o/p-H2 ratio.

4. Astrochemical modeling

4.1. Our model

We solve the chemistry in the gas phase using a code developed by our group and based on the DLSODES solver, which is part of the ODEPACK package (www.netlib.org). The code was benchmarked against the results of Pagani et al. (2009) and Sipilä et al. (2010). We use the same approach as Sipilä et al. (2010) in building the reaction file for the complete depletion case. We take into account all different spin states of the different molecules and ions (i.e., we consider separately ortho and para species of H₂, D₂, H₂D⁺, D₂H⁺, ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and ortho, para and meta species of ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$). We use the reaction rates between these species of Hugo et al. (2009). In addition we include a simple chemistry of CO and N₂, as in Pagani et al. (2009). We use the recombination rates for all ions of Pagani et al. (2009).

The pseudo-time-dependent chemical model only contains gas phase reactions, except that formation of H₂, HD, and D₂ on the grain surface is also included. We use the “large grain limit” (Lipshtat et al. 2004), with the further assumption that all the influx of H and D atoms is converted into H₂, HD, and D₂ molecules, and the diffusion rates of H and D atoms on the grain surface are approximately the same. With these assumptions, the formation rates of these three molecules can be written explicitly as $\begin{array}{ccc}\mathit{R}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{n}}_{\mathrm{G}}{\mathit{\sigma }}_{\mathrm{G}}\frac{\mathrm{(}{\mathit{V}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{H}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{H}}\mathrm{+}{\mathit{V}}_{\mathrm{D}}{\mathit{n}}_{\mathrm{D}}}\\ \mathit{R}\mathrm{\left(}\mathrm{HD}\mathrm{\right)}& \mathrm{=}& {\mathit{n}}_{\mathrm{G}}{\mathit{\sigma }}_{\mathrm{G}}\frac{{\mathit{V}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{H}}{\mathit{V}}_{\mathrm{D}}{\mathit{n}}_{\mathrm{D}}}{{\mathit{V}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{H}}\mathrm{+}{\mathit{V}}_{\mathrm{D}}{\mathit{n}}_{\mathrm{D}}}\\ \mathit{R}\mathrm{\left(}{\mathrm{D}}_{\mathrm{2}}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{n}}_{\mathrm{G}}{\mathit{\sigma }}_{\mathrm{G}}\frac{\mathrm{(}{\mathit{V}}_{\mathrm{D}}{\mathit{n}}_{\mathrm{D}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{V}}_{\mathrm{H}}{\mathit{n}}_{\mathrm{H}}\mathrm{+}{\mathit{V}}_{\mathrm{D}}{\mathit{n}}_{\mathrm{D}}}\mathrm{·}\end{array}$Here n_G and σ_G are the density and cross section of grain particles, V_H, n_H, V_D, and n_D are the thermal velocity and number density of H and D atoms. In this limit, the rather controversial, detailed formation mechanism of H₂ and the uncertain evaporation and diffusion energy barriers do not enter the equations. The ortho and para forms of H₂ and D₂ are produced according to the statistical weights (3:1 and 2:1).

The influence of the different parameters on the chemistry has been studied in great detail by Walmsley et al. (2004), Flower et al. (2004), Pagani et al. (2009), and Sipilä et al. (2010). Here, we adopt a cosmic ray ionization of 3 × 10^-17   s^-1, a dust-to-gas mass ratio of 1.3 × 10^-2, and a grain radius of 0.1 μm. These parameters are not varied, because of the lack of observational constraints at hand.

We simulate different levels of CO depletion by using different CO initial abundances. We run the models with a constant elemental C and O abundance (i.e. not allowing depletion on the grains during the evolution). In all our models, we assume that depletion affects CO and N₂ abundances at the same level.

4.2. Timescales

The free-fall timescale τ_ff of an isothermal sphere of n = 10⁶ cm^-3 is 3 × 10⁴ yrs. Magnetic fields can increase the dynamical timescale, then controlled by ambipolar diffusion in the early evolution of subcritical cores. Timescales of core evolution can be inferred empirically by comparing the number of prestellar and protostellar cores in observed regions (for a review, see Ward-Thompson et al. 2007,and references therein). These empirical timescales lie roughly in the range 2–5 × τ_ff.

It is instructive to compare the chemical timescales to this dynamical timescale. Figure 7 shows the time evolution of the o/p-H₂ ratio and p-D₂H⁺/o-H₂D⁺ ratio for constant physical conditions. The temperature was fixed to 12 K, and three different densities were considered. As already pointed out by, e.g., Pagani et al. (2009), the longest timescale in the chemistry is the conversion between the two forms of ortho and para H₂. The other ions adjust to the o/p ratio of H₂ typically in less than 10² yrs. The ortho-para interconversion of H₂ takes around 10⁴ yrs (typical timescale for decrease in the o/p-H₂ ratio by a factor 1/e, see Fig. 7).

The free-fall timescale varies as 1/$\sqrt{\mathit{n}}$ and the chemical timescale varies even slower with the density. The timescale for o/p-H₂ conversion is 6 × 10³ yrs at density n = 10⁶ cm^-3, and 9 × 10³ and 1.7 × 10⁴ yrs for densities of 10⁵ and 10⁴ cm^-3, respectively. At high densities, the chemical timescales are somewhat shorter than the free-fall timescales (thus roughly one order of magnitude smaller than the empirical core timescales), and the difference even increases at lower densities. This comparison between timescales justifies at first order considering the chemistry independently of the dynamical evolution, as already discussed by Walmsley et al. (2004).

Fig. 7 Time evolution of the chemistry for T = 12 K. The different colors stand for different densities: n = 104 (blue), 105 (green), 106 cm-3 (red). The thick lines represent the o/p H2 ratio. The thin lines represent the p-D2H+/o-H2D+ ratio. The different line styles stand for different CO abundances: complete depletion (full), xCO = 10-6 (dash), xCO = 10-5 (dash-dot).

Because the evolutionary stage of the H-MM1 core is not constrained, in the following we focus on the study of the chemical composition when the steady state is obtained. This choice avoids including the age of the core as an additional parameter, but does not invalidate the general applicability of our conclusions. As shown in Fig. 7, taking the chemical composition at long times overestimates the p-D₂H⁺/o-H₂D⁺ ratio in case the core is not yet in chemical steady state. This note will be important for the extrapolation of our conclusions in the next sections.

Fig. 8 Steady-state ortho-to-para ratio of H2.

Fig. 9 Chemical model predictions of the p-D2H+/o-H2D+ ratio (same convention as previous figures), at steady state. The thick decreasing lines show the ratio derived from the non-LTE analysis of the observations, in the assumption of different densities and different o/p-H2 ratios (same curves as presented in Fig. 6 but with the same linestyle for both o-p H2 ratios).

4.3. Physical conditions

We could in principle derive the density profile of the core from the continuum emission at 850 μm. However, this would require knowing the temperature profile, which is not known. The Gould Belt Survey (André et al. 2010) made with the Herschel Space Observatory should allow accurate temperature profiles to be derived. Until then, since we cannot constrain the density profile without ad-hoc assumptions on the temperature, we restrict ourselves to running chemical models at different densities, temperatures, and levels of CO depletion. Our goal is to investigate under which average conditions the model can reproduce the observations. A full physical model of the source will be presented in a forthcoming detailed study, after collecting new observational constraints.

Figure 9 shows the prediction of the model for the p-D₂H⁺/o-H₂D⁺ ratio. The thin lines are the model predictions, as a function of temperature, and for different densities and CO depletion factors. At fixed density, the p-D₂H⁺/o-H₂D⁺ ratio increases with the CO depletion level, because HD becomes more and more competitive over CO for reactions with ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$ and isotopologues. The ratio increases with temperature in this temperature range, as already explained in detail by Flower et al. (2004). The three thick lines are the ratios as calculated with RADEX from the observed line intensities (as a function of temperature and for different densities and o/p-H₂ ratios, as already shown in Fig. 6). We discuss in the following section the comparison between the chemical model and the observations.

5. Discussion

Figure 9 allows investigation of the effect of the temperature, density, and CO depletion level on the p-D₂H⁺ / o-H₂D⁺ratio, and inferring under which average conditions our model could reproduce the observed value. It should be noted that the observational ratio also depends on the assumed density, because of the non-LTE effects described in Sect. 3.3.2. Agreement between observations and model can only be reached for densities strictly higher than 10⁵ cm^-3. This is consistent with the relatively high measured average density of the core (5  ×  10⁵ cm^-3, Young et al. 2006). Even at the high density of 10⁶ cm^-3, the model can only reproduce the high observed p-D₂H⁺/o-H₂D⁺ ratio if T > 10 K and if the CO depletion level is substantial. The model points to the fact that the CO abundance should be  <10^-5. This corresponds to a CO depletion level of more than 10 (for T = 13 K). The required depletion is even more severe if the density is a bit lower or if the temperature is closer to 10 K. Complete depletion (i.e.  ≫ 100) is required at 10 K.

Bacmann et al. (2002) have studied the level of CO depletion toward a sample of seven prestellar cores. These cores are located at distances between 120 and 200 pc, and three of them are in the Ophiuchus molecular cloud. The cores have mean densities in the range (0.9–6)  ×  10⁵ cm^-3, similar to H-MM1. Based on IRAM observations of dust continuum and CO isotopologues (at angular resolution of 22′′), the authors find depletion factors at the continuum peak of the cores in the range 4.5–15.5, but this factor drops to 2−5 at offset positions where the continuum is 70% of the peak value. Although the study of Bacmann et al. (2002) was done at a 22′′ resolution (to be compared to our 9–14′′ beams), preventing us from doing an accurate comparison, the average depletion level we deduce in H-MM1 seems to be in the high range. Moreover, the D₂H⁺ emission is extended over 40′′. Should the H₂D⁺ also be extended on that scale, this would imply that the model requires a high CO depletion (≫ 10) over an extended region, contradicting the results of Bacmann et al. (2002). It is clear that the observation of the H₂D⁺ spatial distribution and the direct measurement of the level of CO depletion are required to confirm whether the astrochemical model is indeed unable to describe the roots of deuterium chemistry satisfactorily.

Following the discussion of the linewidths in Sect. 3.1, the level of the thick observational curves in Fig. 9 is in fact likely underestimated. Similarly, following the discussion in Sect. 4.2, the model predictions are overestimated if the core has not yet reached chemical equilibrium. These remarks show that the difficulties of the model might be even more severe than shown in Fig. 9.

The model also seems to exclude kinetic temperatures below 10 K, but several very cold prestellar cores are known (6–7 K, Crapsi et al. 2007; Pagani et al. 2007; Harju et al. 2008). The temperature in H-MM1 still needs to be determined, but this point could be another difficulty of the model.

Further tests of our understanding of the root of deuterium fractionation will depend on the observation of H₂D⁺ at the offset positions where D₂H⁺ was detected and on the independent measure of the kinetic temperature and CO depletion level, for which we have an ongoing APEX project.

6. Conclusions

We presented the first secure D₂H⁺ detection in the ISM, toward the Oph H-MM1 core. The centroid velocity of the line is consistent with the one of the o-H₂D⁺ line. The emission is extended over 40′′ at least. Using an out-of-equilibrium radiative transfer model, we inferred the column densities of both o-H₂D⁺ and p-D₂H⁺. We analyzed the chemistry with a state-of-the-art astrochemical model including all spin states of ${\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$, H₂D⁺, D₂H ⁺ , and ${\mathrm{D}}_{\mathrm{3}}^{\mathrm{+}}$. The model can reproduce the data only if the average density is high (n > few10⁵ cm^-3), the kinetic temperature is between 10 and 13 K, and the depletion level of CO is high (>10 at 13 K,  >100 at 10 K). Further observational studies may show whether these conditions indeed apply to H-MM1, and if the CO depletion level stays high on the large scales at which p-D₂H⁺ emits.

Should new observational studies confirm the failure of the model, new ways will have to be investigated. We list here several possibilities: include a more complete description of the chemistry (treating fully e.g. the CO chemistry and gas-grain interaction), recompute or measure the reaction and collision rates of all (spin separated) species, include a correct treatment of the excitation of the different molecules (at the moment, all adopted reactive and non-reactive collision rates describe the molecules in their rotational ground state) in the chemical model.

Acknowledgments

We thank the technical team who built the CHAMP⁺ receiver for APEX. This work is supported by the German Deutsche Forschungsgemeinschaft, DFG Emmy Noether project number PA1692/1-1.

References

All Tables

All Figures

	Fig. 1 The L1688 region, as mapped by SCUBA (850 μm), from the SCUBA legacy catalogs (see text). The contours are in steps of 3σ (=0.18 Jy/beam). The beamsize is 14′′ (HPBW).
In the text

	Fig. 2 850 μm continuum emission of the H-MM1 core. The map is centered on the position of H-MM1 (see text). The contours are in steps of 3σ (=0.18 Jy/beam). The blue numbers show the different pixels of the CHAMP+ low-frequency array. The blue circles show the D2H+ beams. In red, the H2D+ beam.
In the text

	Fig. 3 H2D+ (upper panel) and D2H+ (lower panel) toward the central position of H-MM1.
In the text

	Fig. 4 D2H+ observations. Each panel shows a pixel of CHAMP+ (as displayed in Fig. 2). The vertical line shows the vlsr of the source as measured with the o-H2D+ line (4.23 km s-1). The y-axis is in Kelvin on the ${\mathit{T}}_{\mathit{a}}^{\mathrm{\ast }}$ scale.
In the text

Fig. 5 a) Column densities at the (0, 0) position, computed from the nonequilibrium radiative transfer assuming the beam filling factor for both molecules is 1, as a function of the physical conditions (kinetic temperature, density). The thick black lines show the column densities computed in the LTE assumption. b) Excitation temperature in the same conditions and color conventions as panel a). The black line represents Tex = Tkin.

In the text

	Fig. 6 p-D2H+/ o-H2D+abundance ratio resulting from the radiative transfer modeling of the observations, for different densities and o/p-H2 ratio.
In the text

	Fig. 7 Time evolution of the chemistry for T = 12 K. The different colors stand for different densities: n = 104 (blue), 105 (green), 106 cm-3 (red). The thick lines represent the o/p H2 ratio. The thin lines represent the p-D2H+/o-H2D+ ratio. The different line styles stand for different CO abundances: complete depletion (full), xCO = 10-6 (dash), xCO = 10-5 (dash-dot).
In the text

	Fig. 8 Steady-state ortho-to-para ratio of H2.
In the text

	Fig. 9 Chemical model predictions of the p-D2H+/o-H2D+ ratio (same convention as previous figures), at steady state. The thick decreasing lines show the ratio derived from the non-LTE analysis of the observations, in the assumption of different densities and different o/p-H2 ratios (same curves as presented in Fig. 6 but with the same linestyle for both o-p H2 ratios).
In the text