© ESO, 2016

1. Introduction

The first stages of the star formation process are dense starless and self-gravitating cores, i.e. the so-called pre-stellar cores (Ward-Thompson et al. 1999, Crapsi et al. 2005). Pre-stellar cores in nearby star-forming regions are typically cold (~10 K), dense (10⁴−10⁷ cm^-3), and quiescent (thermal pressure dominates over turbulent motions; e.g. Benson & Myers 1989, Fuller & Myers 1992, Lada et al. 2008, Keto & Caselli 2008).

Chemical differentiation takes place in pre-stellar cores (see e.g. Bergin & Tafalla 2007, di Francesco et al. 2007, Caselli 2011, for reviews). While CO is the second most abundant molecule in the interstellar medium, it tends to freeze onto dust grains in the dense, cold conditions at the centres of pre-stellar cores (e.g. Caselli et al. 1999, Bacmann et al. 2002). The level of CO depletion is usually measured as f_d = X_ref(CO) /N(CO)·N(H₂), where X_ref(CO) is the reference value of the CO fractional abundance, typically between 1 and 2 × 10^-4 (e.g. Frerking et al. 1982, Lacy et al. 1994). The typical value of the CO depletion factor in pre-stellar cores is 5−20 (Crapsi et al. 2005, Christie et al. 2012).

In such cold and dense gas, deuterated species are preferentially formed (e.g. Caselli & Ceccarelli 2012). H₂D⁺ is responsible for the enhancement of the deuterium fraction in most molecular species and is formed by the deuteron-proton exchange reaction H+HD⇌H₂D⁺+H₂+230 K (Millar et al. 1989). The deuteron-proton exchange reaction is exothermic and does not proceed from right to left at temperatures lower than 30 K if most of the H₂ molecules are in para form (e.g. Pagani et al. 1992). H₂D⁺ then reacts with other species to form deuterated ions via H₂D⁺+A→AD⁺+H₂, where A can be CO, N₂, and other neutral species (Herbst & Klemperer 1973, Dalgarno & Lepp 1984). When CO and other abundant neutral species that destroy H and H₂D⁺ are severely frozen onto dust grain surfaces, the deuterium fraction becomes significant. For example, the deuterium fraction in pre-stellar cores is 5−50%, while the elemental abundance of deuterium is ~1.5 × 10^-5 with respect to hydrogen atoms within 1 kpc of the Sun (Linsky et al. 2006, Caselli 2011).

In particular, the deuterium fraction of N₂H⁺ has been used to identify the earliest phases of star formation as the N₂H⁺ deuterium fraction peaks at the pre-stellar phase and towards the youngest protostars (Crapsi et al. 2005, Emprechtinger et al. 2009, Friesen et al. 2013, Fontani et al. 2014). The deuterium fraction in N₂H⁺ is usually given as R_D = N(N₂D⁺) /N(N₂H⁺), where N(i) is the column density of species i.

L1688 is a nearby, 119 ± 6 pc distant (Lombardi et al. 2008) low-mass star-forming region within the Ophiuchus molecular cloud complex. The multiple star-forming regions within L1688 contain more than 60 dense cores and 200 young stellar objects in different evolutionary stages (Motte et al. 1998, André et al. 2007, Simpson et al. 2008, Pattle et al. 2015, Dunham et al. 2015). L1688 is divided into ten regions (A−I; see Fig. 1) with different environmental properties. For instance, while the gas temperature is relatively constant within each region, it varies significantly from one region to another (≃10−17 K; Friesen et al. 2009).

The deuterium fraction across the entirety of L1688 has not been systematically studied yet. The deuterium fraction has only been measured for a few regions, for example, the B2 region has an average R_D ~ 3% (Friesen et al. 2010). The CO depletion across the whole of Ophiuchus has been found to be relatively low compared to the other Gould Belt star-forming regions with an average value less than 10 (Gurney et al. 2008, Christie et al. 2012).

In this paper, we present observations of N₂D⁺(1−0), N₂D⁺(2−1), N₂H⁺(1−0), C¹⁷O(1−0), and C¹⁷O(2−1) towards 40 cores to measure the deuterium fraction and CO depletion factor across the entire L1688 region. In Sect. 2, details regarding the observations are presented. Section 3 describes the results of hyperfine structure fitting as well as deuterium fraction and CO depletion calculations. In Sect. 4, we discuss the results and their relation to possible environmental effects. The conclusions are given in Sect. 5.

2. Observations and data reduction

Fig. 1 850 μm continuum emission of the L1688 region mapped by the Submillimeter Common-User Bolometer Array (SCUBA; Di Francesco et al. 2008), the beam size is 229. Contour levels go from 0.2 Jy beam-1 in steps of 0.2 Jy beam-1 (3σ). The dense cores studied here are indicated with blue filled squares and young stellar objects with open circles. The (0,0) offsets correspond to the J2000 equatorial position α = 16h00m00s, δ = −24°00′00′′. The fits file used to produce this map is available at http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/data/pub/JCMTSL/scuba_F_850umemi, file name scuba_F_353d1_16d8_850um.emi.fits

Figure 1 shows the L1688 region mapped in 850 μm dust continuum emission (Di Francesco et al. 2008). Forty dense cores, revealed by Motte et al. (1998) with 1.3 mm dust emission mapping, were selected for observation with the IRAM 30 m telescope and are shown with filled blue squares. The names and positions of the cores are given in Table A.1. The majority of the dense cores are starless with the exception of: VLA1623, well-known class 0 source (see e.g. André et al. 1993); E-MM3, young stellar object (YSO) with an edge-on circumstellar disk (Brandner et al. 2000); B1-MM4, B1B2-MM2, B2-MM8, B2-MM10, where the protostar position does not coincide with the millimetre dust peak from Motte et al. (1998) but it is within a half beam of our N₂D⁺(1−0) observations. Figure 1 also shows the positions of YSOs embedded in the cloud as open circles (Motte et al. 1998, Simpson et al. 2008, Dunham et al. 2015). The molecular line observations were performed with the IRAM 30 m telescope in June 1998, July 2000, and December 2004. The following transitions were observed: N₂D⁺(1−0), N₂D⁺(2−1), N₂D⁺(3−2), N₂H⁺(1−0), N₂H⁺(3−2), C¹⁷O(1−0), and C¹⁷O(2−1). The observations were obtained with the AB receiver and the VESPA backend. Typical system temperatures for the (1−0) transition observations were 100−200 K for N₂H⁺−N₂D⁺ and 200−360 K for the C¹⁷O line. The spectral resolution for the N₂H⁺, N₂D⁺, and C¹⁷O (1−0) lines varied from 6.5 to 40 kHz and the angular resolutions were 22, 32.1, and 26.6 arcsec for C¹⁷O(1−0), N₂D⁺(1−0), and N₂H⁺(1−0), respectively (see Table 1). The spectra were taken using the position switching (datasets 051-00 and 188-97) and frequency switching (dataset 066-04 with a frequency throw of 7.8 MHz) modes. In Table 1, the dates of the observing runs are given and each run denoted with a dataset number. Part of N₂H⁺(1−0) spectra from datasets 051-00 and 188-97 are already presented in André et al. (2007).

The data reduction was performed with the CLASS package¹. For each source, there were several spectra of the same line. We adjusted these spectra to have the same central frequency and summed together to improve the sensitivity. The integration time for different lines and objects varies from 4 to 30 minutes. We converted the intensity scale to the main-beam temperature scale according to the beam efficiency values given in Table 1.

The N₂D⁺(1−0), N₂D⁺(2−1), N₂H⁺(1−0), C¹⁷O(1−0), and C¹⁷O(2−1) lines have hyperfine splitting with 15, 40, 15, 3, and 9 components, respectively. As such, the spectra were analysed using the standard CLASS hyperfine structure (hfs) fitting method. The routine computes line profiles with the assumptions of Gaussian velocity distribution and equal excitation temperatures for all hyperfine components. The rest frequencies of the main components, the velocity offsets, and the relative intensities of the hyperfine components of the lines were taken from Frerking & Langer (1981), Pagani et al. (2009), and Dore, L. (priv. comm.). The N₂D⁺(3−2) and N₂H⁺(3−2) spectra have very poor baselines and reconstruction of the signal is not possible, so these data are not considered hereafter.

All spectra were initially fit assuming one velocity component. The hfs fitting routine returns both the rms of the baseline and the region with the spectral line. In case the rms of the spectral line region was greater than the rms of the baseline by a factor of 1.5, we redid the fit with an additional velocity component. This was repeated until the two rms agreed within a factor of 1.5. The largest number of velocity components we needed was three.

3. Results

3.1. Spectra

Fig. 2 N2H+(1−0), N2D+(1−0), and C17O(1−0) spectra towards the observed dense cores, labelled in each panel.

Fig. 2 continued.

The spectra of the N₂D⁺(1−0), N₂H⁺(1−0), and C¹⁷O(1−0) lines are shown in Fig. 2. The N₂H⁺(1−0), C¹⁷O(1−0), and C¹⁷O(2−1) emission was detected towards all 40 observed cores. N₂D⁺(1−0) emission was detected towards 23 out of 33 observed cores and N₂D⁺(2−1) emission was detected towards 25 out of 32 observed cores and the A-MM4 core only had a (2−1) detection. As the aim of the study is the measurement of deuterium fractions and their comparison to CO depletion factors, we focus on the (1−0) transitions for the remainder of the paper, since they have the most similar beam sizes and excitation conditions. Towards five cores (B1-MM3, B1-MM4, B2-MM2, B2-MM8, and F-MM2), the N₂H⁺(1−0) line shows two velocity components. The C¹⁷O(1−0) line towards all of the objects in regions C (except C-Ne and C-MM3) and E, and one object in region A (SM1N), shows two or three velocity components. The N₂D⁺(1−0) line shows two velocity components towards one core, B1-MM4. The results of the hfs fits are given in Tables A.2−A.6.

Fig. 3 VLSR of N2D+(1−0) in the left panel and C17O(1−0) in the right panel as a function of the VLSR of the N2H+(1−0) line. Black dotted lines are the lines of equal VLSR. The protostellar cores are shown with black open circles.

The centroid velocities, V_LSR, are determined from the hfs fitting and vary across L1688 from 3.3 to 4.6 km s^-1 with a velocity generally increasing from region A to F. Figure 3 shows the centroid velocities for the N₂H⁺(1−0), N₂D⁺(1−0), and C¹⁷O(1−0) lines. If two velocity components are detected in N₂H⁺(1−0) and N₂D⁺(1−0) (e.g. towards B1-MM4), both components are plotted in the left panel of Fig. 3. If only N₂H⁺(1−0) shows two components, we instead plot an average of the two components if the single N₂D⁺(1−0) line component appears to be a blend of multiple components (cores B2-MM2, B2-MM8, and F-MM2). We plot only the velocity of the closest N₂H⁺(1−0) component if the N₂D⁺(1−0) component does not appear to be blended with the second N₂D⁺(1−0) component presumably not detected above the noise level (core B1-MM3). The three cases are illustrated in Fig. 4, however, only a small fraction of the points (4/25) shown in Fig. 3 are for locations with different numbers of components detected in the two different tracers. For the right panel of Fig. 3, we plot each N₂H⁺(1−0) component against the closest C¹⁷O(1−0) component.

Fig. 4 Isolated hyperfine components of N2H+ and N2D+ (1−0) towards cores where two velocity components are found. The spectra are centred at the frequency of the isolated component. The dotted lines show the centroid velocities of the N2H+(1−0) components.

For 80% of the cores, the V_LSR of the N₂H⁺(1−0) and N₂D⁺(1−0) lines are within 0.05 km s^-1 of each other. The largest centroid velocity difference is only 0.21 km s^-1. The C¹⁷O(1−0) and N₂H⁺(1−0) V_LSR can differ significantly with discrepancies up to 1 km s^-1, and 51% of the cores have velocities that differ by over 0.1 km s^-1. This suggests that N₂H⁺(1−0) and N₂D⁺(1−0) trace roughly the same gas, while C¹⁷O(1−0) likely traces more extended gas, which is expected considering the widespread distribution of CO in molecular clouds and its freeze-out in dense cold regions.

Fig. 5 Line widths, Δv, of N2D+(1−0) in the left panel and C17O(1−0) in the right panel in comparison with Δv of N2H+(1−0). Black dotted lines are the lines of equal Δv. Different colors show different parts of the cloud. The protostellar cores are shown with black open circles.

The line widths (full width at half maximum; FWHM hereinafter Δv) of N₂H⁺, N₂D⁺, and C¹⁷O (1−0) are shown on Fig. 5. As before, when N₂H⁺ and N₂D⁺ show two velocity components, we plot them separately (B1-MM4). In the case where two velocity components are present in N₂H⁺ and only one is seen in N₂D⁺, we either take the Δv of N₂H⁺ of the component with the closest V_LSR if it appears that one component of N₂D⁺ is missing because of the noise level (B1-MM3), or we take Δv = Δv₁/ 2 + Δv₂/ 2 + | V_LSR1 − V_LSR2| if the N₂D⁺ line appears to be a blend of two components (B2-MM2, B2-MM8, F-MM2). The three cases are shown in Fig. 4. For the comparison of the N₂H⁺ and C¹⁷O line widths, we took the Δv of the component with the closest V_LSR.

The N₂H⁺ and N₂D⁺ (1−0) line widths range from 0.2 to 0.7 km s^-1 except for the blended N₂D⁺ line at B2-MM8 (1.2 km s^-1; see left panel of Fig. 5). While the line widths of the N₂H⁺ and N₂D⁺ (1−0) lines are similar for most of the cores, 84% of the cores are within 0.1 km s^-1 of each other; the C¹⁷O(1−0) line widths are overall larger than those of N₂H⁺(1−0). The median width difference is 0.38 km s^-1 and the median line width ratio is 1.9. Similar to what was found with the line centroids, this suggests that the N₂H⁺ and N₂D⁺ trace the same gas, while C¹⁷O traces different, in particular, more turbulent, gas.

3.2. Non-thermal motions

Fig. 6 Ratio of non-thermal components to thermal components of the N2D+(1−0), N2H+(1−0), and C17O(1−0) lines. The dashed line shows the ratio equal to 1.

Figure 6 presents the ratio of non-thermal components Δv_NT of the N₂D⁺(1−0), N₂H⁺(1−0), and C¹⁷O(1−0) lines and thermal line widths of a mean particle, Δv_T. The non-thermal components are derived from the observed line widths Δv_obs via (1)where k is Boltzmann’s constant, T_k is the kinetic temperature, and m_obs is the mass of the observed molecule (Myers et al. 1991). To measure the non-thermal component, we use the kinetic temperature determined by Friesen et al. (2009) from ammonia observations. For those cores that were not observed in Friesen et al. (2009), we use dust temperatures determined by Pattle et al. (2015), assuming that the dust and gas temperatures are equivalent (assumption valid at volume densities above 10⁴ cm^-3; Goldsmith 2001). For most of the cores where both the gas and dust temperatures have been measured, the two values are indeed similar; the only exception to this is the B2 region, where the dust temperature is a few degrees lower than the gas temperature. This may be due to the effect of protostellar feedback, where shocks produced by outflows entraining the dense gas can heat the gas but not the dust (e.g. Draine 1980). Neither dust nor gas temperatures have ever been estimated for the I-MM1 core, so we adopt 11 K, which is the same as H-MM1, as these two cores have similar characteristics: both cores are relatively isolated and far away from the main source of irradiation and heating (see Sect. 4.5). The kinetic temperatures for all cores are given in Table A.7.

The thermal line widths, Δv_T, for a mean particle with mass 2.33 amu are 0.44−0.63 km s^-1 for typical temperatures of 10−20 K across L1688. The majority of C¹⁷O(1−0) lines are supersonic (78%), while most of N₂H⁺(1−0) and N₂D⁺(1−0) lines are subsonic (80% and 75%). The non-thermal to thermal line width ratio can be as high as 1.5 for N₂H⁺(1−0), 2.5 for N₂D⁺(1−0), and as high as 5 for C¹⁷O(1−0). The most turbulent region is Oph-B2 and the most quiescent regions are Oph-B1, C, E, H, and I. Oph-A contains turbulent as well as relatively quiescent cores. For most cores, N₂D⁺(1−0) and N₂H⁺(1−0) have similar non-thermal components and are narrower than the C¹⁷O(1−0) line.

3.3. Column densities and deuterium fractions

The hfs fits provide values needed to measure the excitation temperature (T_ex) and its error. These values are: the total optical depth, i.e. the sum of the optical depths of the various hyperfine components (τ); the quantity labelled T_ant × τ (see below); the FWHM, Δv of the line; and the centroid velocity relative to the local standard of rest (V_LSR). In case of optically thick lines, T_ant × τ is the total optical depth times the difference between the Rayleigh-Jeans equivalent excitation and background temperatures, while for optically thin lines it is the main beam temperature (T_mb). The T_ex can be calculated as (2)where h is the Planck constant, k is the Boltzmann constant, ν is the frequency of the observed transition, T_bg is the cosmic background temperature (2.7 K), J_ν(T_bg) is the equivalent Rayleigh-Jeans background temperature, and J_ν(T) is the function (3)The calculated excitation temperature depends on the value of τ. In the case of weak lines or low S/N, τ cannot be determined properly and the error of τ (Δτ) is high. In all cases where τ/ Δτ ≤ 3, we consider the lines to be optically thin and fix τ = 0.1 (the minimum opacity value) in CLASS. In this case of optically thin conditions, for N₂H⁺(1−0), the excitation temperature value is assumed to be the average T_ex found for optically thick N₂H⁺(1−0) lines, while for N₂D⁺(1−0) we adopt the (measured or assumed) N₂H⁺(1−0) excitation temperature towards the same dense core.

For optically thick transitions, the column density (N_tot) is given by (4)where λ is the wavelength of the observed transition, A_ul is the Einstein coefficient of the u → l transition, g_l and g_u are the statistical weights of the lower and upper levels, Q_rot is the partition function, and E_l is the energy of the lower level (Caselli et al. 2002b). For linear rotors, g_l and g_u are determined by g_J = 2J + 1, where J is the rotational quantum number. The partition function of linear molecules (such as N₂H⁺ and CO) is given by (5)where E_J = J(J + 1)hB and B is the rotational constant. For rotational transitions with hyperfine structure, τ refers to the total optical depth (given by the sum of the peak optical depths of all the hyperfine components) and Δv to the intrinsic line width. The error on N_tot is given by propagating the errors on Δv, τ, and T_ex in Eq. (4).

For optically thin lines (6)where W is the integrated intensity of the line (7)for a Gaussian line (Caselli et al. 2002b).

In case of non-detection of N₂D⁺, upper limits on the N₂D⁺ column density have been derived based on the 3σ uncertainty (3 σ_W) of the integrated intensity, with (8)where N_ch is the mean number of channels covering the velocity range of all the detected lines and Δv_res is the velocity resolution.

The column densities of N₂H⁺, N₂D⁺, and C¹⁷O are given in Tables A.2−A.6. The N₂D⁺ column densities derived from the (1−0) transition in most cores are larger than those derived from the (2−1) transition on average by only 10%. C¹⁷O column densities calculated with (1−0) lines in most cores are smaller than the C¹⁷O column densities calculated with (2−1) lines (see Sect. 3.4 for details). Figure 7 shows the column densities of N₂H⁺ and N₂D⁺. Where multiple components are detected, we plot the sum of the column densities.

Fig. 7 Column densities of N2D+ and N2H+ with lines of constant deuterium fraction (RD). The black line shows RD = 0.5 and the dotted line shows RD = 0.1. The protostellar cores are depicted with black open circles.

The deuterium fraction is defined as the ratio of column densities, R_D = N_tot(N₂D⁺) /N_tot(N₂H⁺), and has been measured for all cores where the N₂D⁺(1−0) line is detected. To calculate the deuterium fraction for cores where two velocity components are detected, we take the sum of the column densities derived from the two components. L1688 overall exhibits high levels of deuterium fractions with a large spread of values between different cores (R_D = 2−43%; see Fig. 8). Such a high level of deuteration (over 20%) was previously found towards other dense cores in different star-forming regions (e.g. Crapsi et al. 2005, Pagani et al. 2007, Emprechtinger et al. 2009, Fontani et al. 2011, 2014, Miettinen et al. 2012, Friesen et al. 2013). The deuterium fraction across the B2 region was previously studied by Friesen et al. (2010), who found slightly lower deuterium fractions (1−10%) than we found (3−18%). The column densities of N₂D⁺ (0.5−6 × 10¹¹ cm^-2) and N₂H⁺ (4−10 × 10¹² cm^-2) these authors obtain are also smaller to the column densities that we found (6−33 × 10¹¹ and 9−27 × 10¹² cm^-2). This difference could be because of different transitions used to calculate column densities, while using similar excitation temperatures derived from the N₂H⁺(1−0) line. Friesen et al. (2010) used N₂D⁺(3−2) and N₂H⁺(4−3) lines with 11′′ and 13′′ HBPW, while we use the (1−0) transition for both species (32′′ and 27′′). The factor of 2 difference in deuterium fraction could also arise from dense core coordinates that are not coincident.

3.4. CO depletion factor

Fig. 8 Deuterium fraction across L1688. The black points show the cores with no N2D+(1−0) detection. The coordinate offsets correspond to the J2000 equatorial position α = 16h24m16s, δ = −24°00′00′′.

In cold, dense, quiescent gas, CO freezes out onto dust grains and the level of this depletion is commonly expressed as a CO depletion factor, f_d, calculated as (9)where X_ref(CO) is the reference abundance and X(CO) is the observed abundance.

The reference abundance of CO in the local ISM has been found to be between 1 and 2 × 10^-4 (Wannier 1980, Frerking et al. 1982, Lacy et al. 1994). We use X_ref(C¹⁶O) = 2 × 10^-4 (Frerking et al. 1982), X(C¹⁸O) /X(C¹⁷O) = 4.11 (Wouterloot et al. 2005), and X(C¹⁶O) /X(C¹⁸O) = 560 (Wilson & Rood 1994) such that (10)To calculate the f_d in case of multiple velocity components, we consider the sum of the column densities of the individual components (as N(H₂) is derived from the millimetre dust continuum emission, which does not contain kinematic information). Since millimetre dust continuum emission is generally optically thin, the molecular hydrogen column density can be derived from the continuum flux density (11)

where S_λ is the flux in a single beam, Ω is the main beam solid angle, μ_H2 = 2.8 is the mean molecular weight per H₂ molecule (Kauffmann et al. 2008), m_H is the mass of atomic hydrogen, κ_λ is the dust opacity per unit mass column density at a given wavelength (κ_{850 μm} = 0.01 cm² g^-1; Johnstone et al. 2000), and B_λ(T_dust) is the Planck function for a dust temperature T_dust (Motte et al. 1998).

The 850 μm dust continuum emission flux measurements from the SCUBA survey (with the beam size of the convolved map of 229; Di Francesco et al. 2008) have been used to calculate N(H₂), adopting the dust temperature from Pattle et al. (2015) when available and the gas temperature from Friesen et al. (2009) in the other cases. We adopted the dust temperature estimated by Motte et al. (1998) for those cores not studied in the above mentioned papers. For Oph-I, where no dust or gas temperature has been measured, we assumed 11 K, which is the same as in Oph-H as Oph-I has similar characteristics to Oph-H, as already mentioned. The dust temperatures for all cores are given in Table A.7.

The depletion factor of CO is generally low in L1688, ranging from 0.2 to 2 in the B1, B1B2, C, E, and F regions, and from to 0.7 to 7.3 in the A, B2, H, and I regions. This result is consistent with previous large scale (e.g. Christie et al. 2012) and small scale (e.g. Bacmann et al. 2002, Gurney et al. 2008) studies. However, larger CO depletion factors are found for the A region compared to the work of Gurney et al. (2008): 1 to 6 instead of 1.5 to 4.5. This small discrepancy could be the result of slightly different pointings or from the use of different CO isotopologues and transitions. In particular, our C¹⁷O(2−1) observations tend to produce column densities larger than N_tot(C¹⁷O(1−0)) by an average factor of 1.5 in Oph-A and 1.1 in B1, B1B2, B2, C, E, and F regions (and thus depletion factors would be lower by 1.5 and 1.1), suggesting that temperature and density gradients along the line of sight may be present, which slightly affects the derived depletion factor depending on the CO transition used. To calculate column densities of C¹⁷O, LTE is assumed with the kinetic temperature equal to the dust temperature.

Christie et al. (2012) suggest that the low depletion factor of L1688 could be due to an unusual dust grain size distribution with a population of very large dust grains and very small spinning dust grains, which reduces the surface area available for freeze-out. However, L1688 is a complex region with active star formation and externally irradiated. Below, we discuss possible causes of the general low CO depletion factors and the significant variation in R_D and f_d found across L1688.

4. Discussion

4.1. Deuterium fraction

We have found a large range of deuterium fractions across L1688, from a minimum of 2% to a maximum of 43%. Previous studies of deuterium fraction did not show such a big spread: for example, 5−25% in Taurus (Crapsi et al. 2005), 1−10% in Ophiuchus B (Friesen et al. 2010), and 3−25% in Perseus (Emprechtinger et al. 2009, Friesen et al. 2013). The large values of deuterium fractions found towards some of the cores in L1688 in our more extensive survey may indicate the presence of centrally concentrated pre-stellar cores on the verge of star formation. The Ophiuchus molecular cloud is known to be denser on average than other nearby star-forming regions (e.g. Lada et al. 2013). The higher average densities, together with the generally higher molecular cloud temperatures (Pattle et al. 2015, Liseau et al. 2015), imply larger pressures that may accelerate the formation of denser cores and the rate of star formation (e.g. Kennicutt & Evans 2012) compared to the other nearby molecular cloud complexes. These high R_D values are found towards E-MM2d, C-Ne, I-MM1, H-MM1, B1-MM1, B1-MM3, and B1B2-MM1, which do not prominently appear in the 850 μm map and they are all relatively isolated structures (I and H) or in between bright sub-millimetre clumps (E, C-N, B1, B1B2). They are probably recently formed cold and dense structures on the verge of star formation.

4.2. Deuterium fraction and CO depletion

Deuterium fraction and CO depletion factor are expected to correlate because CO is one of the main destruction partners of H and its deuterated forms (Dalgarno & Lepp 1984). This correlation has been presented in several theoretical works (Crapsi et al. 2005, Caselli et al. 2008, Kong et al. 2015) and confirmed with observations (Crapsi et al. 2005, Emprechtinger et al. 2009, Friesen et al. 2013). Figure 9 shows the deuterium fraction as a function of CO depletion factor in L1688. The dotted curve on the figure shows the prediction from simple modelling by Crapsi et al. (2005), shifted to take into account the different reference CO abundance adopted here. We emphasize that the f_d values go below 1, which suggests that either our adopted X_ref(CO) has been underestimated by a factor of a few (2−3, considering the lowest f_d value in Table A.7); or the dust opacity, κ, which depends on the evolutionary stage and the properties of dust grains (see Henning et al. 1995), has been overestimated by a factor of 1.5−2; or the dust temperature has been overestimated by a few (1−3) K. However, to enable a comparison with recent literature work (where X_ref(CO) = 2 × 10^-4; see e.g. Hernandez et al. 2011, Fontani et al. 2012), and considering the factor of 2 uncertainty associated with X_ref(CO) (see also Miotello et al. 2014), we did not modify X_ref(CO), although the calculated f_d values may be underestimated by a factor of a few. The important message here is to see if previously found trends are reproduced and/or if L1688 hosts dense cores with a larger variety of chemical/physical properties than found in other nearby star-forming regions.

From Fig. 9, it is evident that two classes of cores are present in L1688. One group includes the A, B2, and I cores, which show a correlation between R_D and f_d similar to that found by Crapsi et al. (2005). The other group contains the B1, B1B2, C, E, F, and H cores, which completely deviate from the Crapsi et al. (2005) correlation. This group of cores shows R_D = 12−43% and f_d = 0.2−4.4. From the parameter space exploration of Caselli et al. (2008) and Kong et al. (2015), large values of R_D (>0.02) cannot be achieved in standard conditions if little CO freeze-out is present in the same gas traced by N₂D⁺. How do we reconcile theory with observations? One possibility is that the dense and cold regions responsible for the bright N₂D⁺ lines are smaller in size than the IRAM-30 m beam at 3 mm. In this case, the CO-depleted zone would be too diluted to be clearly detected within the larger scale CO-emitting region. Indeed, 2 of 6 cores with R_D> 20% and f_d< 4.4 (B1-MM3 and B1B2-MM1) are estimated to be 1300−1800 AU in size, which is less then 2618 AU corresponding to 22′′ at 119 pc, one unresolved (B1-MM1) (Motte et al. 1998), one (H-MM1) has no size estimate (outside of the mapped area in Motte et al. 1998). Higher angular resolution observations of dust continuum and molecular lines are needed to prove this point.

4.3. Deuterium fraction and non-thermal motions

Fig. 9 Deuterium fraction as a function of CO depletion factor. The protostellar cores are depicted with black open circles.

Fig. 10 Deuterium fraction as a function of the non-thermal component of the N2H+(1−0) line width (left) and the C17O(1−0) line width (right). The protostellar cores are depicted with black open circles.

To understand further the characteristics of the highly deuterated, but CO-rich cores in Fig. 9, we plot the deuterium fraction as a function of non-thermal line width of N₂H⁺(1−0) and C¹⁷O(1−0) in Fig. 10. This figure shows that these cores preferentially occupy the left area of the panels, indicating that on average they have narrower C¹⁷O(1−0) and especially N₂H⁺(1−0) lines. Thus, the highly deuterated cores are overall more quiescent than the rest of the sample, which is in agreement with their relatively isolated nature and maybe smaller size, as mentioned in the previous section. Also, these are the cores with relatively small differences between C¹⁷O(1−0) and N₂H⁺(1−0) LSR velocities (see Fig. 11), again suggesting quiescent conditions. The rest of the sample displays broader line widths, suggestive of faster internal motions (in case of gravitational contraction) or external stirring, e.g. due to proximity to active sites of star formation. Indeed, relatively large N₂H⁺(1−0) line widths have been found by Crapsi et al. (2005) towards some of the most evolved starless cores in their sample (L1544 and L429; see their Fig. 6). Line widths tend to increase towards the centre of L1544 (Caselli et al. 2002a) because of contraction motions. The quiescent and highly deuterated cores found in L1688 may then represent an earlier evolutionary stage, compared to L1544 and other contracting pre-stellar cores, where the core has just started to become centrally concentrated but contraction has not started yet (or it has not affected scales large enough to be detected with the current single-dish observations). High angular resolution observations are needed to investigate this conclusion.

4.4. R_D and f_d versus molecular hydrogen column density and temperature

The measured R_D and f_d values are plotted as a function of molecular hydrogen column density, N(H₂), and dust temperature, T_dust, in Fig. 12. No correlation is found for R_D vs. N(H₂), while f_d appears to increase with N(H₂) with different slopes depending on the region. In particular, f_d is increasing faster with N(H₂) in the B2 region compared to the A and C regions. For the other regions, it is hard to see any trend, which is probably because of the more limited range of N(H₂) values detected. One possible cause of the different slopes in the f_d−N(H₂) correlations is the different amount of external heating due to the proximity of Oph-A to HD 147889 (see next section and Fig. 13). This extra illumination maintains the dust grains in region A at a higher temperature compared to the other L1688 regions (see also Liseau et al. 2015), so larger column/volume densities are needed to reach dust temperatures low enough (<25 K) to allow CO molecules to freeze-out onto the dust grains. The C region, which instead follows the same trend as A region, is next closest to HD 147889 after the A region (see Fig. 13), and probably is also illuminated by the star, although the temperature (≤15 K) and column density (≤10²² cm^-2) of region C are as low as in the other regions.

The difference in dust temperatures among the various regions in L1688 also causes the scatter plot in the right panels of Fig. 12. Here we stress that the different amount of external illumination impinging on the different L1688 regions causes the well-known f_d vs. T_dust correlation (e.g. Kramer et al. 1999) to disappear. Oph-A is the only region in L1688 with dust temperatures larger than 15 K and significant CO freeze-out. Oph-B2 displays a sharp drop of f_d with increasing T_dust, which is expected given the exponential dependence on T_dust of the CO evaporation rate (see e.g. Hasegawa et al. 1992).

The deuterium fraction is also expected (and has been measured) to drop with dust temperatures above about 15−20 K (e.g. Emprechtinger et al. 2009, Caselli et al. 2008, Kong et al. 2015). However, this trend is not observed in L1688, as shown in the top left panel of Fig. 12. Once again, the non-uniform conditions among the various regions (in particular, the amount of external illumination, gas volume density, etc.) make it difficult to see a well-defined pattern. Even within the same regions, we do not notice any trend, which may be caused by dust temperatures (mainly derived from Herschel data) that not are representative of the cold regions within which the deuterium fractionation is taking place. More detailed and higher resolution data are needed to explore this possibility.

Fig. 11 Deuterium fraction as a function of the difference between the central velocities of the C17O(1−0) and N2H+(1−0) lines ΔVLSR. The protostellar cores are depicted with black open circles.

4.5. Distance to heating sources

It is well known that star formation in L1688 is affected by the OB-association Sco OB 2 (Pattle et al. 2015), which is located ~145 pc from the Sun (de Zeeuw et al. 1999), i.e. ~25 pc behind L1688. Sco OB 2 is a moving group of more than 120 stars, which mostly consists of B and A spectral-type stars. This group of stars occupies an area of about 15° diameter on the sky (de Zeeuw et al. 1999). We looked for a correlation between R_D, f_d, the distance to nearby stars (ρ Oph, HD 147889, V 2246 Oph, Oph S1), and the closest YSOs (see Fig. 13). ρ Oph is a multiple system of B2IV−B2V spectral type stars; it is a member of the Sco OB 2 association, located to the north of L1688 and is the most distant from L1688 among the four nearby stars. The cores in Oph-A, the nearest region to ρ Oph, show a tentative correlation between the CO depletion factor and the projected angular distance D to the ρ Oph system with f_d increasing with D (see Fig. 14). The cores in the others sub-regions do not show any correlation.

The other three stars are pre-main sequence (PMS) or very young stars related to the Ophiuchus star-forming region. Oph S1 is a B4-K8 type binary star, where a fainter component is a T-Tauri star (Gagné et al. 2004), to the east of the A region; V 2246 Oph is a Herbig Ae/Be star on the west side of the A region; and HD 147889 is a young B2 star on the west side of L1688 (e.g. van den Ancker et al. 1997, Liseau et al. 1999). In all regions of L1688, f_d does not show any correlation with the distance to the brightest young and PMS stars of L1688.

The A region contains more YSOs than other regions and is the closest to the external heating sources. Oph-A has the smallest values of R_D and the highest values of f_d in L1688. Relatively low deuteration could be due to this proximity to nearby irradiating sources, while the large f_d values may be because this region is also the densest in L1688, thus, probably harbouring the densest cores. Here, the dust temperatures are probably low enough (<25 K) to allow molecular freeze-out to proceed at a higher rate due to the larger gas-dust collision rates (because of the overall larger densities), but not large enough to promote deuterium fractionation. This is likely because of an increase of the ortho-to-para H₂ ratio in warmer environments (see e.g. Flower et al. 2006). An alternative to the low values of R_D found in Oph-A could be that N₂D⁺ cores are compact and small compared to our beam, while N₂H⁺ is extended and abundant owing to the large average densities (see e.g. Friesen et al. 2014, Liseau et al. 2015), as also found in infrared dark clouds where massive stars and star clusters form (Henshaw et al. 2014). Indeed, Friesen et al. (2014) detected compact (a few hundred AU in size) dust continuum condensation and ortho-H₂D⁺ emission towards one of the cores embedded in Oph-A. However, the line widths of N₂H⁺(1−0) and N₂D⁺(1−0) lines are the same within the errors (see Fig. 5), so this alternative scenario may be harder to justify. Again, higher angular resolution observations of N₂D⁺ are needed to disentangle these two scenarios.

Fig. 12 Left upper panel: deuterium fraction depending on the molecular hydrogen column density, N(H2). Right upper panel: deuterium fraction as a function of dust temperature taken from Pattle et al. (2015), Friesen et al. (2009), and Motte et al. (1998). Left lower panel: CO depletion factor versus the molecular hydrogen column density. Right lower panel: CO depletion factor as a function of dust temperature. The protostellar cores are depicted with black open circles; the color coding is given in Fig. 11.

Fig. 13 Left upper panel: deuterium fraction depending on the angular distance to HD 147889. Right upper panel: deuterium fraction as a function of the angular distance to the closest YSO. The positions of YSOs were taken from Motte et al. (1998), Simpson et al. (2008), and Dunham et al. (2015). Left lower panel: CO-depletion factor depending on the angular distance to HD 147889. Right lower panel: CO depletion factor depending on the angular distance to the closest YSO. The color coding is given in Fig. 11.

No correlation between R_D and distance to any heating source (embedded or external) was found, although the largest R_D values are found at projected distances larger than 50 arcsec from embedded YSOs and greater than 20 arcsec from HD 147889. These are the already discussed large R_D – low f_d relatively isolated cores, which we discussed in Sect. 4.2. Moreover, the region closest to ρ Oph and V 2246 Oph, the sub-region Oph-A, has the lowest deuterium fractions.

Friesen et al. (2010) studied the deuterium fraction in the B2 region and found a correlation between R_D and distance to the nearest protostar. The deuterium fraction we measure in this region (0−18%, median ~12%) is systematically higher than that found by Friesen et al. (2010) (0−10%, mostly <4%). According to our study, there is no correlation between R_D and the distance to the closest YSO (see Fig. 13 right upper panel). This difference in the results of the studies (R_D− YSO correlation) might be because we studied fewer cores in this particular region (nine cores instead of full mapping of the B2 region, as done by Friesen et al. 2010).

Fig. 14 CO depletion factor as a function of angular distance to the ρ Oph system.

5. Conclusions

This paper presents single point observations of the ground state transitions of N₂H⁺, N₂D⁺, and C¹⁷O towards pre-stellar cores in L1688. We measure the deuterium fraction, R_D, and CO depletion factor, f_d, and study the correlation between these two parameters as well as their dependence on physical parameters, which vary across the cloud; these include dust temperature, molecular hydrogen column density, level of turbulence, projected distance to stellar sources that externally irradiate the cloud, and embedded sources that can internally stir and shock the gas.

We found that the L1688 cores show a large spread of deuterium fractions, 2−43%, and moderate CO depletion factors of up to 7. The reference value of the CO fractional abundance adopted here, 2 × 10^-4, however, may be underestimated by a factor of 2−3 or the dust opacity value used, 0.01 cm² g^-1, may be overestimated by a factor of 1.5−2.

The largest R_D values are found towards (B1, B1B2, C, E, H, I) dense cores, which present relatively quiescent (subsonic) motions as measured from the width of the high density tracer lines; thus, they are probably in an evolutionary stage just preceding the contraction towards protostellar birth as found in other more evolved pre-stellar cores such as L1544 in Taurus (Caselli et al. 2002a). These highly deuterated cores are also relatively isolated (they are typically found in between strong sub-millimetre dust continuum emission and far away from embedded protostars) and the CO freeze-out is low, in contrast with chemical model predictions. This dichotomy (large R_D and low f_d) can be understood if the deuterated gas is confined in a region smaller than the beam size, so that the CO-depleted region is too small to be revealed with observations of the widespread C¹⁷O(1−0) emitting gas. Higher angular resolution observations are needed to confirm this statement.

Except for this sub-group of highly deuterated cores, widths are generally supersonic in C¹⁷O(1−0) (78%) and subsonic in N₂H⁺(1−0) and N₂H⁺(1−0) (75% and 80%). The B2 region stands out here with supersonic widths in all tracers in the majority of cores. This is probably a combination of internal systematic motions (e.g. contraction) and external stirring (turbulence, interaction with outflows driven by embedded protostars).

The correlation between R_D and f_d already found in other studies of starless cores is maintained for a sub-group of cores (those in A, B2, and I, when plotted together). The highly deuterated cores show a significantly steeper rise of R_D with f_d, suggesting that the CO observations are not sensitive to the CO-depleted zone, as mentioned above.

The densest region in L1688, Oph-A, hosts dense cores with significant amount of CO freeze-out (f_d close to 7) but no corresponding large R_D values. This can be explained if the gas and dust temperatures are low enough (<25 K) to allow CO freeze-out but high enough to significantly increase the ortho-to-para ratio compared to cooler regions. Alternatively, the N₂D⁺ cores may be small and diluted within our beam, whereas the N₂H⁺ is abundant all across the region as a consequence of the large average densities.

The nearby ρ Oph star in the Sco OB2 association appears to affect the amount of CO freeze-out in Oph-A cores because f_d increases with projected distance to this star. No other regions appear to be chemically affected by their proximity to external stars or embedded young stellar objects, except for the B2 region, where f_d decreases with increasing dust temperature.

Our observations hint at the present of compact starless cores (smaller than those found towards less dense and cooler molecular cloud complexes, such as Taurus) with large deuterium fractions (12−43%) and small CO depletion factors (0.2−4.4). These cores have relatively low temperatures compared to their surroundings, so they do not appear as clear structures in dust continuum emission maps at 850 μm. Their compact nature may be a consequence of the overall higher densities and temperatures across L1688 compared to other nearby star-forming regions. High angular resolution observations are needed to test these predictions.

Acknowledgments

The authors thank J. Pineda for thoughtful comments, which helped to improve the manuscript. The authors acknowledge the financial support of the European Research Council (ERC; project PALs 320620); Andy Pon acknowledge that partial salary support was provided by a CITA National Fellowship.

References

Appendix A: Hyperfine splitting fit results

In Tables A.2–A.6, the quantity labelled T_ant × τ, is τ × (J_ν(T_ex) − J_ν(T_bg)) in the case of optically thick transition (τ> 0.1) or the main beam temperature (T_mb) in the case that the line is optically thin (with an adopted value of τ = 0.1). For the details, see Sect. 3.3.

All Tables

All Figures

Fig. 1 850 μm continuum emission of the L1688 region mapped by the Submillimeter Common-User Bolometer Array (SCUBA; Di Francesco et al. 2008), the beam size is 229. Contour levels go from 0.2 Jy beam-1 in steps of 0.2 Jy beam-1 (3σ). The dense cores studied here are indicated with blue filled squares and young stellar objects with open circles. The (0,0) offsets correspond to the J2000 equatorial position α = 16h00m00s, δ = −24°00′00′′. The fits file used to produce this map is available at http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/data/pub/JCMTSL/scuba_F_850umemi, file name scuba_F_353d1_16d8_850um.emi.fits

In the text

	Fig. 2 N2H+(1−0), N2D+(1−0), and C17O(1−0) spectra towards the observed dense cores, labelled in each panel.
In the text

	Fig. 2 continued.
In the text

	Fig. 3 VLSR of N2D+(1−0) in the left panel and C17O(1−0) in the right panel as a function of the VLSR of the N2H+(1−0) line. Black dotted lines are the lines of equal VLSR. The protostellar cores are shown with black open circles.
In the text

	Fig. 4 Isolated hyperfine components of N2H+ and N2D+ (1−0) towards cores where two velocity components are found. The spectra are centred at the frequency of the isolated component. The dotted lines show the centroid velocities of the N2H+(1−0) components.
In the text

	Fig. 5 Line widths, Δv, of N2D+(1−0) in the left panel and C17O(1−0) in the right panel in comparison with Δv of N2H+(1−0). Black dotted lines are the lines of equal Δv. Different colors show different parts of the cloud. The protostellar cores are shown with black open circles.
In the text

	Fig. 6 Ratio of non-thermal components to thermal components of the N2D+(1−0), N2H+(1−0), and C17O(1−0) lines. The dashed line shows the ratio equal to 1.
In the text

	Fig. 7 Column densities of N2D+ and N2H+ with lines of constant deuterium fraction (RD). The black line shows RD = 0.5 and the dotted line shows RD = 0.1. The protostellar cores are depicted with black open circles.
In the text

	Fig. 8 Deuterium fraction across L1688. The black points show the cores with no N2D+(1−0) detection. The coordinate offsets correspond to the J2000 equatorial position α = 16h24m16s, δ = −24°00′00′′.
In the text

	Fig. 9 Deuterium fraction as a function of CO depletion factor. The protostellar cores are depicted with black open circles.
In the text

	Fig. 10 Deuterium fraction as a function of the non-thermal component of the N2H+(1−0) line width (left) and the C17O(1−0) line width (right). The protostellar cores are depicted with black open circles.
In the text

	Fig. 11 Deuterium fraction as a function of the difference between the central velocities of the C17O(1−0) and N2H+(1−0) lines ΔVLSR. The protostellar cores are depicted with black open circles.
In the text

Fig. 12 Left upper panel: deuterium fraction depending on the molecular hydrogen column density, N(H2). Right upper panel: deuterium fraction as a function of dust temperature taken from Pattle et al. (2015), Friesen et al. (2009), and Motte et al. (1998). Left lower panel: CO depletion factor versus the molecular hydrogen column density. Right lower panel: CO depletion factor as a function of dust temperature. The protostellar cores are depicted with black open circles; the color coding is given in Fig. 11.

In the text

Fig. 13 Left upper panel: deuterium fraction depending on the angular distance to HD 147889. Right upper panel: deuterium fraction as a function of the angular distance to the closest YSO. The positions of YSOs were taken from Motte et al. (1998), Simpson et al. (2008), and Dunham et al. (2015). Left lower panel: CO-depletion factor depending on the angular distance to HD 147889. Right lower panel: CO depletion factor depending on the angular distance to the closest YSO. The color coding is given in Fig. 11.

In the text

	Fig. 14 CO depletion factor as a function of angular distance to the ρ Oph system.
In the text