A&A 460, 709-720 (2006)
DOI: 10.1051/0004-6361:20066105

Searching for massive pre-stellar cores through observations of N2H+ and N2D+[*],[*]

F. Fontani1 - P. Caselli2 - A. Crapsi3 - R. Cesaroni 2 - S. Molinari4 - L. Testi2 - J. Brand1

1 - INAF, Istituto di Radioastronomia, CNR, Via Gobetti 101, 40129 Bologna, Italy
2 - INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
3 - Leiden Observatory, Postbus 9513, 2300 RA Leiden, The Netherlands
4 - INAF, Istituto di Fisica dello Spazio Interpalenatrio, Via Fosso del Cavaliere, 00133 Roma, Italy

Received 25 July 2006 / Accepted 27 August 2006

Aims. We have measured the deuterium fractionation and the CO depletion factor (ratio between expected and observed CO abundance) in a sample of high-mass protostellar candidates, in order to understand whether the earliest evolutionary stages of high-mass stars have chemical characteristics similar to those of low-mass ones. It has been found that low-mass starless cores on the verge of star formation have large values both of the column density ratio $N({\rm N_2D^+})/N({\rm N_2H^+})$ and of the CO depletion factor.
Methods. With the IRAM-30 m telescope and the JCMT we have observed two rotational lines of N2H+ and N2D+, the (2-1) line of C17O and DCO+, and the sub-millimeter continuum towards a sample of 10 high-mass protostellar candidates.
Results. We have detected N2D+ emission in 7 of the 10 sources of our sample, and found an average value $N({\rm N_2D^+})/N({\rm N_2H^+})\sim 0.015$. This value is $\sim$3 orders of magnitude larger than the interstellar D/H ratio, indicating the presence of cold and dense gas, in which the physical-chemical conditions are similar to those observed in low-mass pre-stellar cores. The integrated CO depletion factors show that in the majority of the sources the expected CO abundances are larger than the observed values, with a median ratio of 3.2.
Conclusions. In principle, the cold gas that generates the N2D+ emission can be the remnant of the massive molecular core in which the high-mass (proto-)star was born, not yet heated up by the central object. If so, our results indicate that the chemical properties of the clouds in which high-mass stars are born are similar to their low-mass counterparts. Alternatively, this cold gas could be located in one (or more) starless core (cores) near the protostellar object. Due to the poor angular resolution of our data, we cannot distinguish between the two scenarios.

Key words: stars: formation - ISM: molecules

1 Introduction

The initial conditions of the star formation process are still poorly understood. In recent years, studies of starless low-mass cores have begun to unveil the physical and chemical features that lead to the formation of low-mass stars (Kuiper et al. 1996; Caselli et al. 1999; Evans et al. 2001; Caselli et al. 2002a,b; Bergin et al. 2002; Tafalla et al. 200220042006; Harvey et al. 2003; Crapsi et al. 2005). On the other hand, the characterisation of the earliest stages of the formation of high-mass stars is more difficult than for low-mass objects, given their shorter evolutionary timescales, larger distances, and strong interaction with their environments.

Various authors (Molinari et al. 2000, 2002; Sridharan et al. 2002; Beuther et al. 2002; Fontani et al. 2005) have performed extensive studies aimed at the identification of precursors of ultracompact (UC) H II regions, i.e. very young (<105 yr), massive (M>8 $M_\odot$) objects which have not yet ionised the surrounding medium. In particular, from CS and mm continuum observations, it has been noted that in the earliest stages prior to the onset of massive star formation, the radial distribution of the intensity is quite flat, resembling the structure of starless cores in more quiescent and less massive molecular clouds (Beuther et al. 2002). More evolved objects show more centrally-peaked density structures, with power-law indices once again resembling those found in low-mass cores (i.e. $n_{\rm H_2}$ $\propto r^{-1.6}$; Shu 1977; Motte & André 2001). These results strongly suggest that one should apply to the high-mass regime the investigative techniques successful in the study of low-mass star forming cores.

It has been found that when the starless core is on the verge of dynamical collapse, most of the C-bearing species, including CS and CO, are frozen onto dust grains. CS observations of low-mass cores clearly show that this molecule avoids the central high density core nucleus, where the mm continuum peaks (Tafalla et al. 2002). The morphology of the C17O (1-0) and (2-1) line emission of IRAS 23385+6053, a reliable example of a precursor of an UC H II region (Fontani et al. 2004a), suggests that CO is depleted in the high-density nucleus traced by the continuum emission. On the other hand, species such as NH3, N2H+ and N2D+ are good tracers of the gaseous counterpart of the dust continuum emission (Tafalla et al. 2002; Caselli et al. 2002b; Crapsi et al. 2004, 2005), suggesting that they are not significantly affected by freeze-out (see Bergin et al. 2001 and Crapsi et al. 2004 for some evidence of N2H+ depletion toward the center of B68 and L1521F, respectively). This is probably due to the fact that the parent species N2, unlike CO, maintains a large gas phase abundance at densities around 106 cm-3, despite the recent laboratory findings that N2 and CO binding energies and sticking coefficients are quite similar (Öberg et al. 2005; Bisschop et al. 2006). The authors mentioned above also found a relation between the column density ratio $D_{\rm frac}=N(\mbox{N$_{2}$ D$^{+}$ })/N(\mbox{N$_{2}$ H$^{+}$ })$ and the core evolution: $D_{\rm frac}$ is predicted to increase when the core evolves towards the onset of star formation, and then it is expected to drop when the young stellar object starts to heat up its surroundings. Summarising, pre-stellar cores or cores associated with very young massive stars should show strong emission in the N2H+ and N2D+ lines, with large values of $D_{\rm frac}$, and little or no emission of molecular species such as CO and CS, if they evolve like their low-mass counterparts.

Although Sridharan et al. (2005) and Beltrán et al. (2006) have detected several candidate high-mass starless cores in the neighbourhood of some massive protostar candidate, a sample of well-identified high-mass pre-stellar cores is lacking. At present, the objects that are believed to be the closest to the earliest stages of the high-mass star formation process are those of two samples of high-mass protostellar candidates, selected from the IRAS-PSC by Molinari et al. (1996) and Sridharan et al. (2002), and then investigated by various authors (Molinari et al. 19982000; Brand et al. 2001; Fontani et al. 2004a,b; Beuther et al. 2002; Williams et al. 2004; Fuller et al. 2005). These studies have shown that a large fraction of these sources are indeed very young massive objects. Also, assuming that high-mass stars form in clusters (Kurtz et al. 2000), we may expect to find other cold and dense cores located close to our target sources. For this reason, we have observed with the IRAM-30 m telescope and SCUBA at JCMT a selected sample of candidate precursors of UC H II regions in N2H+, N2D+, C17O and sub-mm continuum. The sources have been selected from the two samples mentioned above on the basis of observational features indicative of very little evolution: they are associated with massive ( $M\sim 10^{2}~ M_{\odot}$) and compact (diameters less than 0.1 pc) molecular cores, in which the kinetic temperatures are around $\sim$20 K, and show faint or no emission in the radio continuum. Also, they have distances less than 5 kpc.

Section 2 describes the observations and the data reduction, while the results are presented in Sect. 3 and discussed in Sect. 4. A summary of the main findings is given in Sect. 5.

2 Observations and data reduction

2.1 Molecular lines

All molecular tracers were observed with the IRAM-30 m telescope. In Table 1 we give the molecular transitions observed (Col. 1), the line rest frequencies (Col. 2), the telescope half-power beam width (HPBW, Col. 3), the channel spacing (Col. 4) and the total bandwidth (Col. 5) of the spectrometer used. The observations were made in wobbler-switching mode. Pointing was checked every hour. The data were calibrated with the chopper wheel technique (see Kutner & Ulich 1981), with a calibration uncertainty of ${\sim}20\%$.

Table 1: Observed transitions.
molecular frequency HPBW $\Delta v^{\dagger}$ Bandwidth $^{\dagger}$
transition (GHz) ( $^{\prime\prime}$) (km s-1) (km s-1)
N2H+ (1-0) 93.173772 26 0.031/3.218 113/824
N2H+ (3-2) 279.511863 9 0.042/1.073 75/275
N2D+ (2-1) 154.217137 16 0.038/1.944 68/498
N2D+ (3-2) 231.321966 10 0.051/1.296 91/332
C17O (2-1) 224.714385 11 0.052/1.334 93/342
DCO+ (2-1) 144.077319 16 0.041/2.081 73/533

$^{\dagger}$ The two values refer to the two spectrometers used.

Table 2: Source list and detection summary.
IRAS RA (J2000) Dec. (J2000) $V_{\rm LSR}$ $d_{\rm kin}$ N2H+ $^{\dagger}$ N2D+ C17O DCO+ (sub-)mm
Source (h m s) ($^{\circ}$ ' '') (km s-1) (kpc)          
IRAS 05345+31571 05 37 52.4 32 00 06 -18.4 1.8 Y (0,-5) Y Y Y Y
IRAS 05373+23491 05 40 24.5 23 50 55 2.3 1.2 Y (-8,0) Y N Y Y
IRAS 18144-17231 18 40 24.5 -17 22 13 47.3 4.3 Y (-7,0) N Y(?) N Y
IRAS 18511+01461 18 53 38.1 01 50 27 56.8 3.9 Y (0,-5) Y $^{\clubsuit}$ Y Y Y
IRAS 18517+04372 18 54 13.8 04 41 32 43.9 2.9 Y (12,0) N Y N Ya
IRAS 19092+08411 19 11 38.6 08 46 30 58.0 4.5 Y (8,-4) N Y N Y
IRAS 20126+41042 20 14 26.0 41 13 32 -3.8 1.7 Y (-5,0) Y Y Y(?) Ya,b
IRAS 20216+41072 20 23 23.8 41 17 40 -2.0 1.7 Y (-5,5) Y Y N Ya,b
IRAS 20343+41292 20 36 07.1 41 40 01 11.5 1.4 Y (-8,4) (13,0) Y Y Y(?) Ya,b
IRAS 22172+55491 22 19 09.0 56 04 59 -43.8 2.4 Y (-5,0) Y $^{\clubsuit}$ Y N Y
1 $V_{\rm LSR}$ and $d_{\rm kin}$ from NH3 observations (Molinari et al. 1996). 2 $V_{\rm LSR}$ and $d_{\rm kin}$ from NH3 observations (Sridharan et al. 2002). $^{\dagger}$ Between brackets, the offset (in arcsec) of the N2H+ (1-0) line emission peak with respect to the map center (coordinates in Cols. 2 and 3) is given. $^{\clubsuit}$ Detected in the N2D+ (2-1) transition only. (?) Marginal detection. a Observed by Beuther et al. (2002) at 1.2 mm. b Observed by Williams et al. (2004) at 850 $\mu $m and/or 450 $\mu $m.

Table 3: Sources observed during the first observing run only. Columns 6-9 give the peak of each line or the rms level of the spectrum, both in K.
IRAS RA (J2000) Dec. (J2000) $V_{\rm LSR}$ $d_{\rm kin}$ N2H+ (1-0) N2H+ (3-2) N2D+ (2-1) N2D+ (3-2)
Source (h m s) ($^{\circ}$ ' '') (km s-1) (kpc)        
IRAS 05490+26582 05 52 12.9 26 59 33 0.8 2.1 0.4 ${\leq} 0.9$ ${\leq} 0.19$ ${\leq} 0.3$
IRAS 18567+07001 18 59 13.6 07 04 47 29.4 2.16 0.5 ${\leq} 0.4$ ${\leq} 0.13$ ${\leq} 0.2$
IRAS 20106+35451 20 12 31.2 35 54 46 7.8 1.6 0.5 0.45 ${\leq} 0.05$ ${\leq} 0.05$
IRAS 20205+39482 20 22 21.9 39 58 05 -1.7 4.5 $\leq$0.03 ${\leq} 0.1$ ${\leq} 0.05$ ${\leq} 0.05$
IRAS 20278+35211 20 29 46.9 35 31 39 -4.5 5.0 0.3 ${\leq} 0.2$ ${\leq} 0.06$ ${\leq} 0.08$
IRAS 22134+58342 22 15 09.1 58 49 09 -18.3 2.6 1.2 1.6 ${\leq} 0.06$ ${\leq} 0.08$
IRAS 23139+59392 23 16 09.3 59 55 23 -44.7 4.8 0.9 1.5 ${\leq} 0.16$ ${\leq} 0.3$
IRAS 23151+59122 23 17 21.0 59 28 49 -54.4 5.7 ${\leq} 0.08$ ${\leq} 0.5$ ${\leq} 0.13$ ${\leq} 0.36$
IRAS 23545+65082 23 57 05.2 65 25 11 -18.4 0.8 ${\leq} 0.07$ ${\leq} 0.8$ ${\leq} 0.15$ ${\leq} 0.42$

1 $V_{\rm LSR}$ and $d_{\rm kin}$ from NH3 observations (Molinari et al. 1996). 2 $V_{\rm LSR}$ and $d_{\rm kin}$ from NH3 observations (Sridharan et al. 2002).

2.1.1 N2H+ and N2D+

A sample of 19 sources, selected from two sample of protostar candidates as explained in Sect. 1, were observed in July 2002, obtaining single-point spectra of the N2H+ (1-0) and (3-2) lines, and N2D+ (2-1) and (3-2) lines towards the position of the IRAS source. Then, according to the results obtained, in a second observing run carried out in July 2003, we mapped in on-the-fly mode the 10 brightest sources in both the N2H+ transitions, in order to have a clear identification of the emission peak position. Finally, spectra of the N2D+ (2-1) and (3-2) lines were obtained with a deep integration towards the N2H+ (1-0) line peak position. These sources are listed in Table 2, while those observed during the first run only are listed in Table 3.

The on-the-fly maps of the N2H+ (1-0) and (3-2) transitions were obtained in a $1\hbox{$^\prime$ }\times 1\hbox{$^\prime$ }$ field centered on the position listed in Table 2. The two transitions were mapped at the same time and registered with two spectrometers: one with low spectral resolution and large bandwidth and a second with higher spectral resolution and narrower bandwidth (see Table 1). The antenna temperature, $T_{\rm A}^{*}$, and the main beam brightness temperature, $T_{\rm MB}$ are related as $T_{\rm A}^{*}=T_{\rm MB}~\eta_{\rm MB}$, with $\eta_{\rm MB}=B_{\rm eff}/F_{\rm eff}=0.80$ and 0.48 for the (1-0) and (3-2) lines, respectively. The spectra of the N2D+ (2-1) and (3-2) lines were obtained simultaneously, again with two spectrometers with different spectral resolutions and bandwidths (see Table 1). The values of $\eta_{\rm MB}$ are 0.73 and 0.57 for the (2-1) and (3-2) transitions, respectively.

The N2H+ (1-0), (3-2), and N2D+ (2-1), (3-2) rotational transitions have hyperfine structure. To take this into account, we fitted the lines using METHOD HFS of the CLASS program, which is part of the GAG software developed at the IRAM and the Observatoire de Grenoble. This method assumes that all the hyperfine components have the same excitation temperature and width, and that their separation is fixed to the laboratory value. The method also provides an estimate of the total optical depth of the lines, based on the intensity ratio of the different hyperfine components. The frequency of the N2H+ (1-0) line given in Table 1 is that of the main component (Dore et al. 2004), while that of the (3-2) transition has been taken from Crapsi et al. (2005), and it refers to the $F_{1} F = 4\;5\rightarrow 3\;4$ hyperfine component, which has a relative intensity of $17.46\%$. Rest frequencies of the N2D+ lines given in Table 1 refer to the main hyperfine component (Dore et al. 2004).

2.1.2 C17O and DCO+

On-the-fly maps of the C17O (2-1) and DCO+ (2-1) transitions were obtained on July 2003. These two transitions and the N2H+ (1-0) and (3-2) lines were mapped at the same time and with two spectrometers. The data were reduced using the CLASS program. Like the N2H+ lines, the C17O (2-1) line has hyperfine structure. However, the spectra were too noisy to resolve the different components, so that the lines were fitted with a single Gaussian.

2.2 Sub-millimeter continuum

Continuum images at 850 $\mu $m of 5 of the 10 sources mapped in N2H+ were taken on 1998 October with SCUBA at the JCMT (Holland et al. 1998). The standard 64-points jiggle map observing mode was used, with a chopper throw of 2$^\prime$ in the SE direction. Telescope focus and pointing were checked using Uranus and the data were calibrated following standard recipes as in the SCUBA User Manual (SURF). One of the maps of this dataset, that of IRAS 22172+5549, was published by Fontani et al. (2004b).

3 Results

Parameters and detection statistics of the 10 sources observed more extensively (see Sect. 2.1.1) are listed in Table 2. Column 1 gives the IRAS name, and the equatorial (J2000) coordinates of the maps center are listed in Cols. 2 and 3. In Cols. 4 and 5 we give the source velocities and kinematic distances, respectively. Columns 6-10 give the following information: detection (Y) or non-detection (N) in N2H+, N2D+, C17O, DCO+ and millimeter or sub-millimeter continuum, respectively. In Col. 6 we also give the offset of the position of the N2H+ (1-0) line emission peak with respect to the map center. As explained in Sect. 2.1.1, the N2D+ spectra have been taken toward this position. For the lines mapped, we have considered as detected those sources showing emission above the 3$\sigma$ level in the map field. Main parameters of the 9 sources observed only in the first observing run are listed in Table 3: the equatorial coordinates in Cols. 2 and 3 represent the observed position, while source velocity, $V_{\rm LSR}$, and kinematic distance, $d_{\rm kin}$, are listed in Cols. 4 and 5, respectively. In Cols. 6-9 we also give the peak temperature of the lines detected, or the rms level of the spectrum for those undetected.

In the following, we will consider only the results obtained for the sources listed in Table 2. Since they have been previously observed in different tracers by other authors, for completeness in Table 4 we give the most important physical parameters reported in the literature: the dust temperatures $T_{\rm d}$, obtained from gray-body fits to the observed spectral energy distribution of the source, the gas kinetic temperature $T_{\rm kin}$, derived from NH3 observations (Molinari et al. 1996; Jijina et al. 1999) or CH3C2H observations (Brand et al. 2001), and H2 column- and volume densities computed from continuum observations.

Table 4: Main physical parameters of our sources found in literature: dust temperature ($T_{\rm d}$), H2 column and volume densities ( $N_{\rm H_2}$ and $n_{\rm H_2}$, respectively) and gas kinetic temperature ( $T_{\rm kin}$). The references are shown between brackets.
Source $T_{\rm d}$ $N_{\rm H_2}$ $n_{\rm H_2}$ $T_{\rm kin}$
  (K) ( $\times 10^{22}$ cm-2) ( $\times 10^{5}$ cm-3) (K)
05345+3157 - 8.6 (1) 1.5 (1) 17.0 (2)
05373+2349 27(1) 5.9 (1) 2.4 (1) 21.2 (3)
18144+3157 29(1) 17.7 (1) 1.8 (1) 23.6 (3)
18511+0146 35 (1) 6.2 (1) 0.4 (1) 27.3 (4)
18517+0437 38 (5) 68 (6) 5.4 (6) 29 (3)
19092+0841 29 (1) 12.4 (1) 1.6 (1) 40.4 (4)
20126+4104 62 (5) 52 (6) 11.4 (6) 23.0 (5)
20216+4107 46 (5) 18 (6) 3.1 (6) 21.0 (5)
20343+4129 44 (5) 22 (6) 3.9 (6) 18 (5)
22172+5549 30 (1) 4.8 (1) 0.8 (1) 17.5 (3)
(1) Molinari et al. (2000);
(2) Jijina et al. (1999);
(3) Molinari et al. (1996);
(4) Brand et al. (2001);
(5) Sridharan et al. (2002);
(6) Beuther et al. (2002).

3.1 N2H+ and N2D+ spectra at the N2H+ emission peak

In the following, we will discuss the data taken with the spectrometer with the highest spectral resolution. All sources listed in Table 2 have been detected both in the N2H+ (1-0) and (3-2) transitions. Five have also been detected in both N2D+ (2-1) and (3-2) lines, and two (IRAS 18511+0146 and IRAS 22172+5549) in the N2D+ (2-1) transition only. All spectra of the N2H+ (1-0) and N2D+ (2-1) lines are shown in the Appendix, in Figs. A.1 and A.2.

In Table 5 we give the N2H+ (1-0) and (3-2) line parameters of the spectra taken at the peak position of the (1-0) line emission: in Cols. 3-7 we list integrated intensity ( $\int T_{\rm MB}{\rm d}v$), peak velocity ( $V_{\rm LSR}$), FWHM, opacity of the main component ($\tau_{10}$), and excitation temperature ( $T_{\rm ex}$) of the N2H+ (1-0) line, respectively. Columns 9-12 show the same parameters for the N2H+ (3-2) line, for which we do not list $T_{\rm ex}$ because we will not use it in the following. The integrated intensities have been computed over the velocity range given in Cols. 2 and 8, while for the other parameters we have adopted the fitting procedure described in Sect. 2.1.1. This procedure cannot provide $T_{\rm ex}$ for optically thin lines. However, for the sources with opacity ${\geq} 0.4$, the values of $T_{\rm ex}$ are comparable to the kinetic temperature $T_{\rm kin}$ listed in Col. 5 of Table 4, within a factor $\sim$2 (with the exception of IRAS 19092+0841, for which the difference is a factor of 4). Therefore, for sources with optically thin lines, we decided to assume $\mbox{\mbox{$T_{\rm ex}$ }} = T_{\rm kin}$.

Most of the hyperfine components of the N2H+ lines are not well resolved. This is due to the fact that the line widths ( ${\sim}1 {-} 3$ km s-1) are larger than the separation in velocity of most of the different components. As an example, we show in Fig. 1 the spectra of IRAS 05345+3157 with the position of the hyperfine components. Regarding the assumption of equal excitation conditions made to treat the hyperfine structure of the N2H+ (1-0) line, Daniel et al. (2006) found with theoretical calculations that this assumption is not valid in clouds with densities in the range 104-106 cm-3 and line optical depths greater than 20. Since all of our sources have comparable densities (see Sects. 3.3 and 3.4) but the lines have much smaller opacities (see Table 5), the assumption of equal excitation conditions is considered to be valid.

The line parameters of the N2D+ transitions, derived using the same fitting procedure adopted for N2H+, are listed in Table 6. We give an upper limit for the integrated intensity of the lines for the undetected sources, using the observed rms and an average value of the line FWHM of 1.5 km s-1, both for the (2-1) and (3-2) transitions. With respect to Table 5, we do not give the line opacities and $T_{\rm ex}$ because in most of the spectra the uncertainties are comparable to or larger than the values obtained. The two sources with well-determined opacity of the main component of the N2D+ (2-1) line are IRAS 20126+4104 and IRAS 20216+4107, for which we obtain $\tau=0.13\pm 0.09$ and $0.96\pm 0.03$, respectively. For the remaining ones, i.e. IRAS 05345+3157, IRAS 05373+2349, IRAS 18511+0146, IRAS 20343+4139 and IRAS 22172+5549, we have fitted the lines forcing the optical depth to be 0.1. The observed LSR velocities and line widths are in good agreement with those derived from N2H+, indicating that the two molecular species trace the same material.

{\includegraphics[angle=0,width=9cm]{} }\end{figure} Figure 1: Spectra of IRAS 05345+3157 in N2H+ (1-0) and (3-2) lines ( left panel) and N2D+ (2-1) and (3-2) lines ( right panel). The lines under the N2H+ spectra indicate the position of the hyperfine components.

Table 5: N2H+ line parameters.
  N2H+ (1-0)       N2H+ (3-2)    

vel. range $\int T_{\rm MB}{\rm d}v$ $v_{\rm LSR}$ FWHM $\tau_{10}$ $T_{\rm ex}$   vel. range $\int T_{\rm MB}{\rm d}v$ $v_{\rm LSR}$ FWHM $\tau_{32}$
  (km s-1) (K km s-1) (km s-1) (km s-1)   (K)   (km s-1) (K km s-1) (km s-1) (km s-1)  
05345+3157 -28, -9 14.66 -17.97 1.59 0.5 $\pm$ 0.2 22.0   -23, -13 16.33 -18.03 1.75 4.3 $\pm$ 0.5
05373+2349 -8, 11 20.86 2.31 1.57 0.8 $\pm$ 0.3 20.14   -2, 7 26.13 2.26 1.47 3.7 $\pm$ 0.2
18144-1723 36, 60 38.44 48.01 3.47 0.49 $\pm$ 0.03 26.23   42, 53 31.41 48.42 3.01 $\pm$ 1
18511+0146 46, 68 26.76 57.12 2.86 0.1 -   52, 62 22.11 56.91 3.33 $\pm$ 1
18517+0437 33, 54 26.01 43.84 2.72 0.1 -   39, 48 28.38 43.71 3.21 0.7 $\pm$ 0.3
19092+0841 45, 70 13.20 57.64 3.76 0.5 $\pm$ 0.1 10   50, 63 13.13 57.07 5.84 0.1
20126+4104 -14, 7 36.15 -3.87 2.01 0.1 -   -8, 1 39.92 -3.36 1.57 7.82 $\pm$ 0.04
20216+4107 -12, 7 9.39 -1.66 1.30 1.2 $\pm$ 0.3 9.4   -6, 2 14.41 -1.89 0.99 8.6 $\pm$ 0.5
20343+4139a 0, 22 14.22 10.87 1.43 1.79 $\pm$ 0.07 9.2   7, 14 14.45 10.81 1.21 0.1
20343+4139b 0, 22 18.39 11.48 2.44 1.37 $\pm$ 0.07 9.3   7, 14 21.51 11.75 3.69 0.1
22172+5549 -54, -34 6.97 -43.59 2.04 0.1 -   -49, -39 14.65 -43.56 2.63 0.1

a Peak at (-8 $^{\prime\prime}$, 4 $^{\prime\prime}$). b Peak at (13 $^{\prime\prime}$, 0 $^{\prime\prime}$).

Table 6: N2D+ line parameters.
  N2D+ (2-1)   N2D+ (3-2)
source vel. range $\int T_{\rm MB}{\rm d}v$ $v_{\rm LSR}$ FWHM   vel. range $\int T_{\rm MB}{\rm d}v$ $v_{\rm LSR}$ FWHM
  (km s-1) (K km s-1) (km s-1) (km s-1)   (km s-1) (K km s-1) (km s-1) (km s-1)
05345+3157 -19, -16 0.33 -17.68 1.01   -20, -15 0.44 -17.70 1.57
05373+2349 0, 4 0.30 2.43 1.65   0, 7 0.43 2.38 1.78
18144-1723   $\leq$0.45         $\leq$0.68    
18511+0146 55, 60 0.23 56.90 0.99     $\leq$0.54    
18517+0437   $\leq$0.60         $\leq$0.42    
19092+0841   $\leq$0.60         $\leq$0.68    
20126+4104 -6, 0 0.49 -4.14 1.52   -6, -1 0.31 -3.86 1.05
20216+4107 -4, 1 0.31 -1.57 1.06   -4, 6 0.25 -11.85 0.72
20343+4139a 0,15 0.48 10.57 0.80   9, 14 0.79 10.66 1.21
20343+4139b 8,15 0.60 11.23 2.93   9, 15 0.84 11.62 2.17
22172+5549 -45, -38 0.23 -43.26 2.2     $\leq$0.14    
a Peak at (-8 $^{\prime\prime}$, 4 $^{\prime\prime}$). b Peak at (13 $^{\prime\prime}$, 0 $^{\prime\prime}$).

Table 7: Deconvolved source angular diameters (in arcseconds) derived from each tracer assuming a Gaussian source.
source N2H+(1-0) N2H+(3-2) C17O (sub-)mm continuum DCO+
05345+3157 31.8 20.8 7.9 21.4 34.4
05373+2349 34.9 21.4 - 14 29.6
18144-1723 23.8 22.2 - 14.8 -
18511+0146 41.3 29.9 41.3 29.3 23.6
18517+0437 25.2 18.3 28.7 29.0* -
19092+0841 - 15.1 17.6 11.7 -
20126+4104 32.0 22.7 24.7 18.0* -
20216+4107 23.6 14.9 14.1 22.7* -
20343+4129 - 21.8 32.5 26.0* -
22172+5549 18.1 - 24.0 15.7 -
* From 1.2 mm continuum (Beuther et al. 2002). They are the geometric mean of the major and minor axes obtained from two-dimensional Gaussian fits to the observed emission.
IRAS 20126+4104, IRAS 20216+4107 and IRAS 20343+4129 have also been observed at 850 $\mu $m by Williams et al. (2004) but they do not estimate the diameters.

3.2 Distribution of the integrated intensity

The maps of the integrated intensity of the lines (N2H+ (1-0) and (3-2), C17O (2-1) and DCO+ (2-1)) are shown in the Appendix, from Fig. A.3 to A.12. For the sources observed with SCUBA, the emission map of each line has been superimposed on the continuum map at 850 $\mu $m, while for the others the position of the (sub-)millimeter peaks found in literature are indicated. We have not shown the maps of the sources marginally detected in C17O and DCO+ (2-1) lines.

Table 8: Source properties derived from (sub-)mm continuum emission: deconvolved angular diameter ( $\theta _{\rm cont}$), integrated 850 $\mu $m continuum flux ($F_{\nu }$, for IRAS 18517+0437 we give the integrated 1.2 mm continuum flux), linear size (D), gas+dust mass ( $M_{\rm cont}$), visual extinction ($A_{\rm v}$) and H2 column and volume densities ( $N_{\rm H_2}$ and $n_{\rm H_2}$), derived from $M_{\rm cont}$ and $A_{\rm v}$.
source $\theta _{\rm cont}$ D $F_{\nu }$ $M_{\rm cont}$ $A_{\rm v}$ $N_{\rm H_2}^{\dagger}$ $N_{\rm H_2}^{\dagger\dagger}$ $n_{\rm H_2}^{\dagger}$ $n_{\rm H_2}^{\dagger\dagger}$
  (arcsec) (pc) (Jy) ($M_\odot$) (mag) ( $\times 10^{23}$ cm-2) ( $\times 10^{23}$ cm-2) ( $\times 10^{5}$ cm-3) ( $\times 10^{5}$ cm-3)
05345+3157 21.4 0.19 4.42 179 220 6.26 2.07 10.9 3.6
05373+2349 14 0.08 2.87 35 239 6.81 2.25 27.8 9.3
18144-1723 14.8 0.31 7.76 1118 496 14.1 4.67 14.8 4.9
18511+0146 29.3 0.55 10.23 959 137 3.89 1.29 2.3 0.8
18517+0437 29 0.41 7.2 $^{{\bf 1}}$ 336 87 2.47 0.8 2.0 0.6
19092+0841 11.7 0.25 5.57 426 283 8.05 2.7 10.3 3.4
20126+4104 18 0.15 21.9 $^{{\bf 2}}$ 504 982 28.0 9.2 61.2 20.2
20216+4107 22.7 0.19 6.1 $^{{\bf 2}}$ 160 196 5.58 1.84 9.7 3.2
20343+4129 26 0.18 19 $^{{\bf 2}}$ 426 586 16.7 5.51 30.7 10.3
22172+5549 15.7 0.18 3.47 239 307 8.73 2.88 15.5 5.2
1 From Beuther et al. (2002), derived from the 1.2 mm continuum.
2 From Williams et al. (2004). $^{\dagger}$ From $M_{\rm cont}$. $^{\dagger\dagger}$ From $A_{\rm v}$.

Figures A.3-A.12 show that in all sources mapped with SCUBA the distribution of the 850 $\mu $m emission and the N2H+ (1-0) and (3-2) integrated line emission are in good agreement. The emission contours at half of the maximum (FWHM) overlap fairly well. For almost all of the sources, the angular separation between the lines emission peaks and that of the sub-mm emission (or the 1.2 mm for IRAS 18517+0437) is much smaller than the beam size, with the exception of the N2H+ (3-2) line in IRAS 18144-1723 and IRAS 18517+0437, for which the peak of the map is offset by ${\sim}12$ $^{\prime\prime}$ from the dust peak. This result confirms that the N2H+ molecule and the dust continuum emission trace similar material, as already found by several authors who observed the N2H+ lines in both low-mass and high-mass star formation regions (Caselli et al. 2002a; Crapsi et al. 2005; Fuller et al. 2005).

The distribution of the C17O (2-1) line integrated emission follows that of N2H+ and the 850 $\mu $m continuum in IRAS 18511+0146, IRAS 18517+0437, IRAS 19092+0841, IRAS 20216+4107 and IRAS 20343+4129, while in IRAS 22172+5549 it has a different distribution and in the other sources the signal is faint and irregular. The integrated intensity of the DCO+ (1-0) line, clearly detected only toward IRAS 05345+3157, IRAS 05373+2349 and IRAS 18511+0146, presents significant differences from one source to another: in IRAS 18511+0146 (Fig. A.6) the DCO+ emission is in good agreement with that of the continuum; in IRAS 05345+3157 it is more extended (see Fig. A.3), while in IRAS 05373+2349 the two tracers are almost totally separated (Fig. A.4).

Assuming that the emission in each tracer has a Gaussian profile, we have derived the angular diameter ($\theta$) of the sources deconvolving the observed FWHM in that tracer with the appropriate Gaussian beam. The results are listed in Table 7. We could not compute the diameters of IRAS 19092+0841 and IRAS 20343+4129 in the N2H+ (1-0) line and of IRAS 22172+5549 in the N2H+ (3-2) line because these are not resolved. The diameters derived from the N2H+ (1-0) line range from 18.1 $^{\prime\prime}$ to 41.3 $^{\prime\prime}$, while those from the (3-2) transition are between 14.9 $^{\prime\prime}$ and 29.9 $^{\prime\prime}$, indicating that the higher excitation transition traces a more compact region. The values of $\theta$ derived from C17O are typically in between those deduced form the N2H+ (1-0) and (3-2) lines, with the exception of IRAS 05345+3157 and IRAS 18517+0437. However, $\theta$ has been derived under the hypothesis of a Gaussian source, which is only a rough assumption for the C17O emission maps. For the sources mapped with SCUBA by Williams et al. (2004), i.e. IRAS 20126+4104, IRAS 20216+4107 and IRAS 20343+4129, the authors do not give any estimate of the angular diameter. For this reason, in Table 7 we list the diameters estimated from 1.2 mm continuum maps, derived calculating the geometric mean of the major and minor axes of the sources obtained by Beuther et al. (2002) from two-dimensional Gaussian fits (see their Table 2).

The diameters derived from the 850 $\mu $m continuum emission are on average smaller than those derived from the 1.2 mm continuum by Beuther et al. (2002). Since the angular resolutions are comparable, we believe that this can be due to the different observation technique. As pointed out in Sect. 2.2, the 850 $\mu $m observations were carried out in "jiggle'' map observing mode, which provides maps with lower sensitivity than the 1.2 mm maps obtained by Beuther et al. (2002) in the dual-beam on-the-fly observing mode, and therefore can be less sensitive to extended structures.

3.3 Physical properties from dust emission

The physical parameters derived form the dust emission are presented in Table 8: in Cols. 2 and 3 we give the angular and linear diameters ( $\theta _{\rm cont}$ and D, respectively); in Cols. 4-6 integrated flux density ($F_{\nu }$), gas+dust mass ( $M_{\rm cont}$) and visual extinction ($A_{\rm v}$) are listed; Cols. 7 and 8 give the H2 column densities ( $N_{\rm H_2}$) derived in two different ways which will be explained in the following, and in Cols. 9 and 10 the corresponding H2 volume densities are given.

The linear diameters have been computed using the kinematic distances listed in Table 2, and are between $\sim$0.08 and $\sim$0.5 pc, typical of clumps hosting high-mass forming stars (Kurtz et al. 2000; Fontani et al. 2005). The flux densities, $F_{\nu }$, are obtained integrating the SCUBA maps by the 3$\sigma$ level in the maps. For IRAS 20126+4104, IRAS 20216+4107 and IRAS 20343+4129 we give the values listed in Williams et al. (2004). For IRAS 18517+0437, observed neither in this work nor by Williams et al. (2004) at 850 $\mu $m, we list the value obtained at 1.2 mm by Beuther et al. (2002).

From the continuum flux density, assuming constant gas-to-dust ratio, optically thin and isothermal conditions, the total gas+dust mass is given by:

 \begin{displaymath}M_{\rm cont}=\frac{F_{\nu} d^2}{k_{\nu}B_{\nu}(T)}\cdot
\end{displaymath} (1)

In Eq. (1), d is the distance, $k_{\nu}$ is the dust opacity coefficient, derived according to $k_{\nu}=k_{\nu_0}(\nu/\nu_0)^{\beta}$ (where $\nu_0=230$ GHz and $k_{\nu_0}=0.005$ cm2 g-1, which implies a gas-to-dust ratio of 100, Kramer et al. 1998), and $B_{\nu}(T)$ is the Planck function calculated at the dust temperature T. For this latter, we decided to use the values listed in Col. 5 of Table 4 (i.e. the gas temperature derived from NH3 or CH3C2H) instead of the dust temperature given in Col. 2 because this latter has been derived from the spectral energy distribution of the sources, which includes the IRAS measurements, not sensitive to temperatures lower than $\sim$30 K. Therefore, this value can overestimate the dust temperature throughout the envelope. We have assumed $\beta=1.5$, which is a plausible "average'' value for high-mass star formation regions (see e.g. Molinari et al. 2000; Hatchell et al. 2000; Williams et al. 2004).

The visual extinction, $A_{\rm v}$, has been derived from the continuum flux density using Eq. (2) of Kramer et al. (2003):

 \begin{displaymath}A_{\rm v}=\frac{1.086}{\kappa^{\prime}}\frac{F_{\rm 850~\mu m}}{B_{\rm 850~\mu m}(T)\Omega_{\rm s}},
\end{displaymath} (2)

where $F_{\rm 850~\mu m}$ is the integrated flux density at 850 $\mu $m, and $B_{\rm 850~\mu m}(T)$ is the Planck function calculated at the temperature T, for which we have taken the gas kinetic temperature (Col. 4 of Table 4). $\Omega_{\rm s}$ is the source solid angle and $\kappa'=\kappa_{850}/\kappa_{\rm V}$is the ratio of absorption and extinction coefficients. The values of $\Omega_{\rm s}$ have been computed from the angular diameters given in Table 7, and we have assumed $\kappa'=1.7\times 10^{-5}$, which is the average value found by Kramer et al. (2003) for $\beta=1.5$. For IRAS 18517+0437, for which $F_{\rm 850~\mu m}$ is not available, we scaled the flux at 1.2 mm to 850 $\mu $m assuming once again a dust opacity index $\beta=1.5$ (i.e. a spectral index $\alpha=3.5$).

Finally, we have determined the molecular hydrogen total column density following two approaches: (1) from $M_{\rm cont}$ and D assuming a spherical source; (2) from $A_{\rm v}$ using the relation given in Frerking et al. (1982):

\begin{displaymath}N_{\rm H_{2}}=0.94\times 10^{21}\times A_{\rm v}({\rm mag})~[{\rm cm^{-2}}].
\end{displaymath} (3)

Both approaches use $F_{\nu }$ but the latter has the advantage of not assuming any gas-to-dust ratio, which is typically affected by large uncertainties.

We find masses distributed between $\sim$30 to $\sim$900 $M_\odot$, visual extinctions from $\sim$90 to $\sim$103 mag, average column densities between $\sim$1023 and $\sim$1024 cm-2 and average volume densities of $\sim$10 5-106 cm-3. Mass, visual extinction and gas column density were previously estimated from the (sub-)mm continuum emission in six sources of our sample by Molinari et al. (2000), and in the remaining ones by Beuther et al. (2002). Even though both Molinari et al. and Beuther et al. used marginally different assumptions in computing the dust opacity, and a different approach to determine $N_{\rm H_2}$ and $A_{\rm v}$, their estimates are in good agreement with ours. Table 8 shows that the column- and volume density estimates derived from $A_{\rm v}$ is $\sim$3 times smaller than those obtained from $M_{\rm cont}$. We believe that this systematic difference is due to the different hypotheses made on the dust absorption, since all the other parameters used in Eqs. (1) and (2) are the same. However, given the large uncertainties associated with these measurements, we conclude that the two estimates of $N_{\rm H_2}$ are in good agreement, indicating that the assumed gas-to-dust ratio of 100 is a reasonable value for our sources.

3.4 N2H+ and N2D+ column densitites and deuterium fractionation

The deuterium fractionation, $D_{\rm frac}$, can be derived by determining the ratio of the column density of a hydrogen-bearing molecule and its deuterated isotopologue. In cold and dense clouds the molecular species used to compute $D_{\rm frac}$ can be affected by depletion. For example, HCO+ and H2CO and their deuterated counterparts are highly depleted in the high-density nucleus of low-mass starless cores (e.g. Carey et al. 1998). Therefore, the deuterium fractionation derived from these species is representative only of an outer, lower-density shell. On the other hand, it has been well established that N2H+ and N2D+ do not freeze-out even in the densest portions of low-mass molecular clouds (e.g. Caselli et al. 2002a,b; Crapsi et al. 2005). Therefore these species should give an estimate of the deuterium fractionation representative also of the densest regions of our sources.

We have derived the N2H+ and N2D+ column densities, $N({\rm N_2H^+})$ and $N({\rm N_2D^+})$, following the method outlined in the Appendix of Caselli et al. (2002b), which assumes a constant excitation temperature, $T_{\rm ex}$. The values of $T_{\rm ex}$ for the N2H+ (1-0) lines have been derived from the hyperfine fitting procedure described in Sect. 2.1.1. For the optically thin lines, for which the fitting procedure cannot provide $T_{\rm ex}$, we have assumed $T_{\rm ex}$= $T_{\rm kin}$ (Col. 5 of Table 4) as already pointed out in Sect. 3.1. For the N2D+ (2-1) line, for which we do not have good estimates of the optical depth, we have assumed the excitation temperature of the N2H+ (1-0) line.

The results obtained are listed in Table 9: the N2H+ column densities are of the order of a few ${\sim} 10^{13}$ cm-2, while the N2D+ column densities are a few ${\sim} 10^{11}$ cm-2. The values of $D_{\rm frac}$ are between 0.004 and 0.02, with an average value of ${\sim} 0.015$. This value is $\sim$3 orders of magnitude larger than the cosmic "average'' value of ${\sim} 10^{-5}$ (Oliveira et al. 2003), and close to that found by Crapsi et al. (2005) in their sample of low-mass starless cores, indicating that also the physical conditions of the gas responsible for the N2D+ emission are similar.

We have also computed the N2H+ abundance relative to H2 by dividing $N({\rm N_2H^+})$ by the molecular hydrogen column densities listed in Col. 8 of Table 8 (i.e. the estimates derived from $A_{\rm v}$). The average value is ${\sim} 2.05\times 10^{-10}$, which is in excellent agreement with the $2\times 10^{-10}$ found by Caselli et al. (2002c) in a sample of 18 low-mass starless cores.

Table 9: N2H+ and N2D+ total column densities ( $N({\rm N_2H^+})$ and $N({\rm N_2D^+})$), deuterium fractionation ( $D_{\rm frac}$), and N2H+ chemical abundance relative to H2( $X_{\rm N_2H^+}$). $X_{\rm N_2H^+}$ has been calculated from the molecular hydrogen column density given in Col. 8 of Table 8.
Source $N({\rm N_2H^+})$ $N({\rm N_2D^+})$ $D_{\rm frac}$ $X_{\rm N_2H^+}$
  ( $\times 10^{13}$ cm-2) ( $\times 10^{11}$ cm-2)   ( $\times 10^{-10}$)
05345+3157 3.14 $\pm$ 0.07 3.2 $\pm$ 0.9 0.010 $\pm$ 0.007 1.5
05373+2349 4.34 $\pm$ 0.06 2.8 $\pm$ 0.8 0.006 $\pm$ 0.001 1.9
18144+3157 9.75 $\pm$ 0.08 $\leq$4.8 $\leq$0.004 2.1
18511+0146 6.47 $\pm$ 0.01 2.5 $\pm$ 0.4 0.004 $\pm$ 0.001 5.0
18517+0437 6.48 $\pm$ 0.02 $\leq$6.7 $\leq$0.01 8.1
19092+0841 1.9 $\pm$ 0.5 $\leq$5.3 $\leq$0.028 0.7
20126+4104 7.62 $\pm$ 0.01 4.9 $\pm$ 0.7 0.006 $\pm$ 0.001 0.8
20216+4107 1.4 $\pm$ 0.2 2.8 $\pm$ 0.5 0.020 $\pm$ 0.007 0.8
20343+4129a 2.2 $\pm$ 0.1 4.4 $\pm$ 0.7 0.020 $\pm$ 0.005 0.4
20343+4129b 2.9 $\pm$ 0.1 5.4 $\pm$ 0.8 0.02 $\pm$ 0.01 0.5
22172+5549 1.226 $\pm$ 0.006 2.1 $\pm$ 0.4 0.017 $\pm$ 0.004 0.4
a Peak at (-8 $^{\prime\prime}$, 4 $^{\prime\prime}$). b Peak at (13 $^{\prime\prime}$, 0 $^{\prime\prime}$).

Table 10: H2 volume densities and N2H+ and N2D+ column densities, derived in the LVG approximation.
  N2H+   N2D+
Source $n_{\rm H_2}$ $N({\rm N_2H^+})$   $n_{\rm H_2}$ $N({\rm N_2D^+})$
  ( $\times 10^{6}$ cm-3) ( $\times 10^{13}$ cm-2)   ( $\times 10^{6}$ cm-3) ( $\times 10^{11}$ cm-2)
05345+3157 2.4 3.2   8.0 8.6
05373+2349 2.8 5.0   9.0 6.7
18144+3157 1.6 11   - -
18511+0146 1.0 3.9   - -
18517+0437 1.7 3.5   - -
19092+0841 1.3 7.2   - -
20126+4104 2.6 8.6   6.4 8.0
20216+4107 1.5 1.1   1.0 6.0
20343+4129a 0.4 4.0   3.0 2.9
20343+4129b 1.0 9.6   6.0 30
22172+5549 4.0 1.4   - -
a Peak at (-8 $^{\prime\prime}$, 4 $^{\prime\prime}$). b Peak at (13 $^{\prime\prime}$, 0 $^{\prime\prime}$).

The N2H+ and N2D+ column densities have been also computed using the LVG program described in the Appendix C of Crapsi et al. (2005), with which we have also estimated the $n_{\rm H_2}$ volume densities. The results are given in Table 10. We note that both the N2H+ and N2D+ column densities are in good agreement with the estimates given in Table 9. It is interesting to compare the H2 volume density estimates given in Tables 8 and 10. The values obtained in LVG approximation are on average ${\sim} 2$ times larger than the estimates given in Col. 9 of Table 8, and ${\sim} 7$ times larger than those listed in Col. 10 of Table 8. However, the sources with higher discrepancy are IRAS 18511+0146 and IRAS 18517+0437, for which the density estimates have large uncertainties: IRAS 18511+0146 is the most distant of our sources, while for IRAS 18517+0437 we have derived the flux at 850 $\mu $m from an extrapolation of the 1.2 mm flux. If we neglect these two sources, the average ratio between the density estimates from Col. 10 of Table 8 and those in Table 10 becomes ${\sim} 4$. Taking into account that the densities derived in the LVG approximation are affected by large uncertainties (even 1 order of magnitude), we conclude that for the sources with more accurate measurements the different density estimates are consistent within the uncertainties, confirming that the N2H+ molecule is a good H2 density tracer.

Crapsi et al. (2005) have found that the H2 volume densities derived from the LVG code are systematically lower than those obtained from dust, suggesting a possible partial depletion of N2H+ in the inner nucleus of their starless cores. However, this nucleus has a diameter of $\sim$2500 AU (see e.g. Caselli et al. 2005), whose contribution is negligile at a distance of some kiloparsecs with our angular resolution. Also, the H2 volume densities derived from dust assume an isothermal cloud, which may not be a good assumption for our sources.

Table 11: Parameters used to determine the CO depletion factor: source Galactocentric distance ( $D_{\rm GC}$), "expected'' C17O abundance ( $X_{\rm C^{17}O}^{\rm E}$), integrated intensity of C17O at the dust peak position ( $I({\rm C^{17}O}))$), beam filling factor ( $\eta _{\nu }$), C17O total column density ( $N({\rm C^{17}O})$), observed C17O abundance ( $X_{\rm C^{17}O}^{\rm O}$) and CO depletion factor ($f_{\rm D}$).
Source $D_{\rm GC}$ $X_{\rm C^{17}O}^{\rm E}$ $I({\rm C^{17}O}))$ $\eta _{\nu }$ $N({\rm C^{17}O})$ $X_{\rm C^{17}O}^{\rm O}$ $f_{\rm D}$
  (kpc) ( $\times 10^{-8}$) (K km s-1)   ( $\times 10^{15}$cm-2) ( $\times 10^{-8}$)  
05345+3157 6.38 8.62 4.10 0.34 6.04 2.9 3.0
05373+2349 6.96 7.38 $\leq$2.38 0.62 $\leq$1.99 $\leq$0.9 $\geq$8.2
18144+3157 8.85 4.62 6.78 0.64 5.61 1.2 3.8
18511+0146 11.5 2.56 12.84 0.93 7.73 6.0 0.4
18517+0437 7.85 6.0 8.11 0.88 7.39 9.2 0.7
19092+0841 5.64 10.6 9.46 0.72 9.03 3.3 3.2
20126+4104 4.67 14.3 5.75 0.83 3.65 0.4 35.8
20216+4107 7.47 6.47 7.13 0.62 5.93 3.2 2.0
20343+4129 6.88 7.54 3.86 0.91 2.17 0.4 18.8
22172+5549 9.42 4.05 2.39 0.83 1.45 0.5 8.1

3.5 CO depletion factor

From the C17O maps we give an estimate of the integrated CO depletion factor, $f_{\rm D}$, defined as the ratio between the "expected'' abundance of CO relative to H2, $X_{\rm CO}^{\rm E}$, and the "observed'' value  $X_{\rm CO}^{\rm O}$:

\begin{displaymath}f_{\rm D}=\frac{X_{\rm CO}^{\rm E}}{X_{\rm CO}^{\rm O}},
\end{displaymath} (4)

where $X_{\rm CO}^{\rm O}$ is the ratio between the observed CO column density and the observed H2 column density.

The C17O column density has been derived from the integrated line intensity at the dust peak, assuming optically thin conditions. This assumption is based on the results of several surveys of C17O performed in high-mass star formation regions with comparable mass and density (e.g. Hofner et al. 2000; Fontani et al. 2005). Under this assumption, and using the Raleigh-Jeans approximation (valid for our frequencies), one can demonstrate that the column density of the upper level, J, is related to the integrated line intensity $I({\rm C^{17}O})$ according to:

 \begin{displaymath}N_J=1.209\times 10^{3}~~ \frac{2J+1}{J}\frac{k}{\eta_{\nu} \mu^{2}\nu }I({\rm C^{17}O})
\end{displaymath} (5)

where k is the Boltzmann constant, $\nu$ the line rest frequency, $\eta _{\nu }$ the beam filling factor and $\mu $ the molecule's dipole moment (0.11 Debye for CO). To compute $\eta _{\nu }$, we have used the angular diameters listed in Table 7 derived from C17O for all sources but IRAS 05345+3157 and IRAS 18144+3157, for which we have taken the diameter of the sub-mm continuum. The total column density has been obtained from:

 \begin{displaymath}N_{\rm C^{17}O}=\frac{N_J}{g_J}Q(T_{\rm ex}){\rm exp}~\left(\frac{E_J}{kT_{\rm ex}}\right)
\end{displaymath} (6)

where gJ and EJ are statistical weight (=2J+1) and energy of the upper level, and $Q(T_{\rm ex})$ is the partition function at the temperature $T_{\rm ex}$. We have used as excitation temperature the kinetic temperatures listed in Col. 4 of Table 4. Then, we computed the observed C17O abundances by dividing $N_{\rm C^{17}O}$ by the H2 column densities listed in Col. 8 of Table 8, i.e. the values obtained from $A_{\rm v}$. We decided to use these latters rather than those derived from $M_{\rm cont}$ because they do not assume any gas-to-dust ratio.

The C17O "expected'' abundance has been obtained for each source taking into account the variation of carbon and oxygen abundances with the distance from the Galactic Center. Assuming the standard value of $9.5\times 10^{-5}$ for the abundance of the main CO isotopologue in the neighbourhood of the solar system (Frerking et al. 1982), we have computed the expected CO abundance at the Galactocentric distance ( $D_{\rm GC}$) of each source according to the relationship:

\begin{displaymath}X_{\rm CO}^{\rm E}=9.5\times 10^{-5}{\rm exp}~(1.105-(0.13~D_{\rm GC}({\rm kpc}))) ,
\end{displaymath} (7)

which has been derived according to the abundance gradients in the Galactic Disk for 12C/H and 16O/H listed in Table 1 of Wilson & Matteucci (1992), and assuming that the Sun has a distance of 8.5 kpc from the Galactic Center. Then, following Wilson & Rood (1994), we have assumed that the oxygen isotope ratio 16O/18O depends on $D_{\rm GC}$ according to the relationship 16O/ $^{18}{\rm O}=58.8\times D_{\rm GC}({\rm kpc})+37.1$; finally, taking the standard ratio 18O/ $^{17}\rm O=3.52$ (Frerking et al. 1982), we have computed the expected abundance of the C17O, $X_{\rm C^{17}O}^{\rm E}$, for each source according to:

 \begin{displaymath}X_{\rm C^{17}O}^{\rm E}=\frac{X_{\rm CO}^{\rm E}}{3.52~(58.8~D_{\rm GC}+37.1)}\cdot
\end{displaymath} (8)

The results are listed in Table 11: Cols. 2 and 3 give the source Galactocentric distance and the corresponding $X_{\rm C^{17}O}^{\rm E}$, respectively; Cols. 4-8 list the C17O (2-1) integrated line intensity ( $I({\rm C^{17}O})$), the beam filling factor ( $\eta _{\nu }$), the C17O column density ( $N({\rm C^{17}O})$), the observed C17O abundance ( $X_{\rm C^{17}O}^{\rm O}$) and the CO depletion factor ( $f_{\rm D}=X_{\rm C^{17}O}^{\rm E}/X_{\rm C^{17}O}^{\rm O}$).

{\includegraphics[angle=-90,width=8cm]{} }
\end{figure} Figure 2: Deuterium fractionation ( $D_{\rm frac}$) versus integrated CO depletion factor ($f_{\rm D}$) for our sources (filled circles) and the low-mass pre-stellar cores of Crapsi et al. (2005) (open squares). The arrows indicate the upper limits on $D_{\rm frac}$ (see Table 9) or the lower limit on $f_{\rm D}$ (see Table 11).

Column 7 of Table 11 shows that in 7 out of 9 sources in which we have detected C17O, $f_{\rm D}$ is smaller than 10, while in the remaining two sources, IRAS 20126+4104 and IRAS 20343+4129, $f_{\rm D}$ is much higher than 10. For IRAS 18511+0146 we obtain an unusual value of 0.4. We believe that this is due to the fact that this source is the most further away from the Galactic Center ( $D_{\rm GC}\sim 12$ kpc), and the expected CO abundance at such a distance could be different from that calculated from Eq. (8). These results indicate that in the majority of our sources the observed abundance is comparable to or marginally lower than the expected value. However, the discussion of these results requires three main comments: first, the values of the "canonical'' CO abundance measured by other authors in different objects vary by a factor of 2 (see e.g. Lacy et al. 1994; Alves et al. 1999). Second, the integrated CO depletion factor is an average value along the line of sight. Therefore, given the large distance to our targets and the poor angular resolution of our data, the effect of not-depleted gas associated with the more external gaseous envelope of the source can significantly affect the measured $f_{\rm D}$, which hence is to be taken as a lower limit. Finally, the angular resolution of our observations also allow us to derive only average values of $f_{\rm D}$ over the sources, which may have complex structure.

4 Discussion

In the following we discuss the parameters presented in the previous sections, and compare them with the findings from the literature.

4.1 Deuterium fractionation and CO depletion

Theoretical models predict that the deuterium fractionation is correlated to the amount of CO depletion (Caselli et al. 2002b; Aikawa et al. 2005). Observations of low-mass objects have partially confirmed this theoretical prediction. Bacmann et al. (2003) have found a good correlation between $f_{\rm D}$ and $D_{\rm frac}$ (derived from D2CO/H2CO) in 5 pre-stellar low-mass cores. More recently, Crapsi et al. (2005) have obtained $D_{\rm frac}$ using ${\rm N_2D^+/N_2H^+}$ in a sample of 31 low-mass starless cores, and they have shown that cores with higher $D_{\rm frac}$ have also higher $f_{\rm D}$.

In order to check if such a correlation also exists in high-mass objects, we have added to Fig. 5 of Crapsi et al. (2005), in which $D_{\rm frac}$ is plotted against the CO depletion factor, the values obtained from our high-mass objects. The result is shown in Fig. 2: all the objects of our sample have $D_{\rm frac}$ smaller than those of the Crapsi et al. sample, and nearly half of them also have smaller $f_{\rm D}$. The plot also indicates that the correlation found by Crapsi et al. (2005) and Bacmann et al. (2003) for low-mass sources cannot be extended to the sources of our sample. In fact, $D_{\rm frac}$ does not show any dependence on $f_{\rm D}$.

Neverthless, when discussing these results we have to keep in mind several caveats. First, the sources of Crapsi et al. (2005) are embedded in low-mass molecular clouds ${\sim} 100 {-} 200$ pc away from the Sun, therefore their estimates of the deuterium fractionation and the CO depletion factor, obtained with the same angular resolution, are less affected than ours by the contribution of the non-deuterated and non-depleted gas along the line of sight associated with the molecular envelope of the source. Additionally, we are comparing a sample of well known starless cores, most of which are in the pre-stellar phase, with a sample of high-mass protostellar candidates, in which the heating produced by the forming protostar may affect both deuteration and CO depletion. Observations of deuterated molecules in low-mass protostars have shown that the deuterium fractionation is similar to the values found in pre-stellar cores, as shown for example by Loinard et al. (2002), Hatchell (2003), and Parise et al. (2006), who have observed formaldehyde, ammonia and methanol towards several low-mass protostars. However, these values are relative to the gaseous envelope in which the low-mass protostars are born, which is thought to maintain the initial conditions of the star formation process, while the heating produced by a forming high-mass protostar is expected to dramatically push down the deuteration and the CO depletion of the molecular environment (see e.g. Turner 1990). Finally, as already mentioned in Sect. 3.5, the angular resolution of the maps allow us to derive only average values of $D_{\rm frac}$ and $f_{\rm D}$ over the sources, which in principle have complex morphology and may host objects in different evolutionary stages (see e.g. Kurtz et al. 2000). The origin of such cold gas in our sources is thus unclear. We further discuss this point in Sect. 4.3.

4.2 Effect of temperature on deuterium fractionation and CO depletion

The deuterium fractionation and the CO depletion factor in a molecular core are related to its physical properties. In particular, in sources with comparable gas density, $D_{\rm frac}$ and $f_{\rm D}$ should be larger in colder clouds. In this section, we will describe a simple chemical model which has been used to interpret the present observational results on the depletion factor, the deuterium fractionation and the gas temperature. Unlike low mass cores, the structure of massive cores cannot be derived in detail with current observations. The physical properties listed in Table 4 and 8 are average values within telescope beams whose angular sizes are comparable to source sizes (compare HPBW in Table 1 with $\theta _{\rm cont}$in Table 8). Therefore, any attempt of modeling massive cores as smooth objects with densities and temperature gradients similar to those found in low-mass cores is subject to large uncertainties, also considering that massive star forming regions are probably clumpy, so that dense and cold material may be confined in small regions with filling factors significantly smaller than unity.

For these reasons, we decided to use a chemical model simpler than the one used for low mass cores (see e.g. Crapsi et al. 2005). In particular, the detailed physical structure is neglected and the clouds are treated as homogeneous objects with density and temperature equal to the average values from Col. 10 of Table 8 and Col. 5 of Table 4, respectively. As already pointed out in Sect. 3.3, we used the gas temperature instead of the dust temperature to compare with our model predictions, given that the latter has been determined using also the IRAS measurements, which are not sensitive to temperatures lower than 30 K (unlike the data used to determine $T_{\rm gas}$), and it is thus expected to overestimate the dust temperature throughout the cloud. If clumps denser and colder than the average values are indeed present, our calculation is expected to underestimate the observed CO depletion factor and the deuterium fractionation in our sources, unless the filling factor is significantly less than unity.

Other simplifications include the determination of the electron fraction, which is now simply assumed to be given by (following McKee 1989):

$\displaystyle x(e)=1.3 \times 10^{-5} \times n({\rm H_2})^{-0.5} .$     (9)

Dust grains (all assumed to be negatively charged[*]) are distributed in size according to Mathis et al. (1977; hereafter MRN) (but with the minimum size value enlarged by a factor of 10, as in the case of low-mass cores; if we assume the MRN value, conclusions do not appreciately change in the range of physical parameters considered here). We also included in the model all the multiply deuterated forms of H3+ and their destruction via dissociative recombination with electrons, recombination onto negatively charged grains, and through reactions with CO, O and N2 has been taken into account. The freeze-out of CO, O and N2 is balanced by thermal (see e.g. Hasegawa et al. 1992) and non-thermal (cosmic-ray impulsive heating, Hasegawa & Herbst 1993) desorption. In analogy with similar work done on low-mass cores, the binding energies used are those for CO and N2 onto H2O (see Öberg et al. 2005), and the oxygen binding energy is fixed at 800 K (Tielens & Allamandola 1987). The cosmic-ray ionization rate used here is $3\times 10^{-17}$ s-1, an appropriate value for high-mass star forming regions (e.g. van der Tak & van Dishoeck 2000).

The model has been run for different temperature values ( $T_{\rm gas} =
T_{\rm dust}$ from 10 to 100 K) and volume density fixed at the average value of the objects in our sample ( $n({\rm H_2}) = 6\times 10^5~{\rm cm^{-3}}$). All models start with undepleted abundances of the neutral species and run for a time given by the free-fall time appropriate for the chosen density ( $t_{\rm ff} = 5\times10^4$ yr), which implies that chemical and dynamical timescales are assumed to be comparable.

Figure 3 shows the model predictions for $D_{\rm frac}$, i.e. the deuterium fractionation in species such as N2H+, and the CO depletion factor $f_{\rm D}$ versus the gas (dust) temperature. The thin curves show models where the Gerlich et al. (2002) rate coefficients for the proton-deuteron exchange reactions are used, whereas the dotted curve is for models with the (factor of  3) larger rate coefficients typically used in chemical models (e.g. Roberts et al. 2003, 2004; but see Walmsley et al. 2004; Flower et al. 2005). The recently measured rate coefficients still suffer from some uncertainties, especially due to the importance of the so-called "back reactions'' of the deuterated forms of H3+ with ortho-H2, which lower the deuteration depending on the unknown ortho/para H2 ratio.

{\includegraphics[angle=0,width=8cm]{} }\end{figure} Figure 3: Top left: deuterium fractionation versus kinetic temperaturederived from NH3 observations. The arrows indicate the deuterium fractionation upper limits. The curves represent the prediction of the theoretical model described in Sect. 4.2. The solid curve is model prediction when using the rate coefficients measured by Gerlich et al. (2002) for the proton-deuteron exchange reactions. The dotted curve is from the standard model, with (a factor of about 3) larger rate coefficients. Dashed curve is from the standard model, when 70% of CO molecules are trapped in H2O ice. Bottom: same as top for the integrated CO depletion factor. Top right: deuterium fractionation versus depletion factor for the same models and data in the previous two panels. Note that the presence of trapped CO does not significantly affect the deuterium fractionation.

The data points in Fig. 3 suggest that the model is in quite good agreement with observations for the majority of the objects, despite its simplicity. Considering the uncertainties in the rate coefficients and in the gas/dust temperature, the observed deuterium fractionation is well reproduced in the observed range of temperatures. In the case of the CO depletion factor, there are three objects that show too large $f_{\rm D}$ for their adopted temperatures: IRAS 20126+4107, IRAS 18144-1723, and IRAS 19092+0841. In the case of IRAS 20126+4107 and IRAS 19092+0841 we note that the excitation temperature derived from N2H+ (1-0) is indeed close to 10 K, suggesting that in these cases the deuterated gas may be confined in smaller and colder clumps, compared to the high mass star forming region where they are embedded. Therefore, lower dust/gas temperatures (close to 10 K) should be more appropriate.

However, for IRAS 18144-1723, the N2H+ (1-0) excitation temperature is quite large ( 26 K), so that some other mechanisms may be present to maintain a large fraction of CO frozen onto dust grain (unless the dust temperature is at least 5 K lower than the gas temperature). One possibility could be the trapping of CO molecules in H2O ice. Indeed, if a large fraction of CO is trapped in water ice, dust temperatures between 30 and 70 K are needed in order to release all the CO back in the gas phase (Collings et al. 2003; Viti et al. 2004). To test this possibility, we run models where a certain fraction of CO molecules were assumed to stick onto dust grains with binding energies appropriate for H2O ice (4820 K; Sandford & Allamandola 1990). In the case of IRAS 18144-1723, to reproduce the observed $f_{\rm D}$ value (3.8) at the measured temperature (23.6 K), 70% of solid CO needs to be trapped into water and this is what is shown in Fig. 3 by the dashed curves, which is also able to reproduce the $f_{\rm D}$ value observed in the warmest source (IRAS 19092+0841). However, the CO trapping does not affect the deuterium fractionation, which is limited by the gas temperature. We then conclude that a significant fraction of CO may be indeed trapped in H2O ice, although more accurate determination of the dust and gas temperature are needed to give more quantitative estimates.

4.3 Where are the pre-stellar cores?

The most important result of this work is the detection of N2D+ emission in 7 of our sources, with values of deuterium fractionation close to those found in low-mass starless cores. In this section we discuss this result, focusing attention on the origin and the nature of the cold and dense gas that gives rise to this emission.

While Crapsi et al. (2005) have established that in their objects the N2D+ emission arises from the densest core nucleus, the angular resolution of our data does not allow us to determine the accurate location of this emission. As already noted in Sects. 4.1 and 4.2, it is well known that high-mass stars typically form in clusters, so that the surroundings of a massive protostar often show complex morphology (see Molinari et al. 2002; Cesaroni et al. 2003; Fontani et al. 2004a) and may be fragmented in objects with different masses and in different evolutionary stages (see e.g. Kurtz et al. 2000; Fontani et al. 2004b), whose angular separation can be comparable to or even smaller than the angular resolution of the observations.

For this reason, in principle the observed cold and dense gas responsible for the N2D+ emission may be associated with the high-mass protostellar object or with another object located very close to it. In the first case, such gas is located in the most external shell of the parental cloud not yet heated up by the high-mass (proto-)star: a "record'' of the early cold phase. Relatively high values of deuterium fractionation have been found also in some hot cores (e.g. Oloffson 1984, Turner 1990), i.e. the molecular environment in which massive stars have been recently formed. Such a scenario would indicate that the molecular cloud in which low-mass and high-mass stars are born have similar chemical and physical conditions. In the second case, the N2D+ emission is due to one or more molecular cores in the pre-stellar phase located close to the central protostar. High-angular resolution observations of one of our targets, IRAS 05345+3157, have clearly shown that the molecular gas is fragmented into several cores (see Molinari et al. 2002), separated on average by $\sim$5 $^{\prime\prime}$-10 $^{\prime\prime}$. Such a separation is smaller than the angular resolution of our data. Therefore, the observed N2D+ emission might arise from some of these cores. Given the comparable distances, in principle this could be the case also for the other sources. To solve this problem, observations with higher angular resolution are required.

5 Conclusions

We have observed several rotational lines of N2H+, N2D+, C17O, DCO+ and sub-mm continuum emission with the IRAM-30 m telescope and the JCMT in 10 high-mass protostellar candidates. Our main goal is to measure the deuterium fractionation through the ratio $N({\rm N_2D^+})/N({\rm N_2H^+})$, and the CO depletion factor (ratio between "expected'' and observed CO abundance), in order to shed light on the chemical and physical properties of the molecular clouds in which high-mass stars are born. The main results of this study are the following:

Our results allow us to conclude that most of our sources host cold and dense gas with properties similar to those found in low-mass starless cores on the verge of the gravitational collapse. The problem to solve now is the location and the origin of this gas. We propose two scenarios: in the first one, the cold gas is distributed in an external shell not yet heated up by the high-mass protostellar object, a remnant of the parental massive starless core. In the second one, the cold gas is located in cold and dense cores close to the high-mass protostar but not associated with it. This scenario is supported by high-angular resolution observations of some of our objects, which have revealed their clumpy structure. Observations at higher angular resolutions will be very helpful in deciding which solution is correct.

We thank C. M. Walmsley for his valuable comments and suggestions. Antonio Crapsi was supported by a fellowship from the European Research Training Network "The Origin of Planetary Systems'' (PLANETS, contract number HPRN-CT-2002-00308) at Leiden Observatory.



6 Online Material

Appendix: spectra and maps

  \begin{figure}\par {\includegraphics[angle=0,width=14cm]{} }
\end{figure} Figure A.1: Spectra of the N2H+ (1-0) line ( left panels) and the N2D+ (2-1) line ( right panels) for IRAS 05345+3157, IRAS 05373+2349, IRAS 18144+3157, IRAS 18511+0146 and IRAS 18517+0437 ( from top to bottom). All spectra have been taken at the position of the N2H+ (1-0) line emission peak.

  \begin{figure}{\includegraphics[angle=0,width=14cm]{} }
\end{figure} Figure A.2: Same as Fig. A.1 for IRAS 19092+0841, IRAS 20126+4104, IRAS 20216+4107, IRAS 20343+4129 and IRAS 22172+5549. For IRAS 20343+4129, we show the spectra obtained towards position a (see Table 5).

  \begin{figure}{\includegraphics[angle=-90,width=16cm]{} }
\end{figure} Figure A.3: Left panel: plot of the N2H+ (1-0) ( top panel) and (3-2) ( bottom panel) integrated emission, superimposed on the 850 $\mu $m map (grey-scale), obtained with SCUBA towards IRAS 05345+3157. Right panel: same as left panel for the C17O (2-1) line ( top panel) and DCO+ (2-1) line ( bottom panel). The levels are: from 20$\%$ to 95$\%$ of the peak (step 15$\%$) for the SCUBA map; from 30$\%$ to 90 $\%$ of the peak (step 10$\%$) for the N2H+ lines; from 40$\%$ to 90$\%$ (step 10$\%$) for the other lines. The contours at half of the maximum are indicated by thick lines for each tracer.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.4: Same as Fig. A.3 for IRAS 05373+2349, which is undetected in C17O.

  \begin{figure}{\includegraphics[angle=0,width=7.5cm]{} }
\end{figure} Figure A.5: Same as Fig. A.3 for IRAS 18144+3157. We do not show the C17O map because it is too noisy, nor the DCO+ map because the source is not detected in this tracer.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.6: Same as Fig. A.3 for IRAS 18511+0146.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.7: Left panel: integrated N2H+ (1-0) ( top panel) and (3-2) ( bottom panel) line emission in IRAS 18517+0437. Right panel: integrated emission of the C17O (2-1) line. In each panel, the filled triangles at map center indicate the position of the 1.2 mm emission peak (Beuther et al. 2002).

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.8: Same as Fig. A.3 for IRAS 19092+0841.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.9: Same as Fig. A.7 for IRAS 20126+4104. The triangle corresponds to the position of the 850 $\mu $m emission peak (Williams et al. 2004).

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.10: Same as Fig. A.9 for IRAS 20216+4107.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.11: Same as Fig. A.9 for IRAS 20343+4129.

  \begin{figure}{\includegraphics[angle=-90,width=11cm]{} }
\end{figure} Figure A.12: Same as Fig. A.3 for IRAS 22172+5549.

Copyright ESO 2006