© ESO, 2011

1. Introduction

Most, if not all, infrared dark clouds (IRDCs; Pérault et al. 1996; Egan et al. 1998) studied so far are fragmented into clumpy structures. Some of the clumps within IRDCs that have been studied with high-resolution interferometric observations are found to contain still denser cores, indicating that the clumps have further fragmented into smaller pieces (Rathborne et al. 2007, 2008; Zhang et al. 2009; Beuther & Henning 2009; Wang et al. 2011)¹.

Massive clumps in IRDCs are promising targets for studying the earliest stages of high-mass star- and star-cluster formation (e.g., Kauffmann & Pillai 2010, and references therein). (Sub)millimetre dust continuum maps of IRDCs have shown that the clumps within them have typical radii, masses, beam-averaged H₂ column densities, and volume-averaged H₂ number densities of  ~0.1−0.5 pc,  ~10² − 10³   M_⊙,  ≳ 10²² cm^-2, and  ~10⁴−10⁵ cm^-3, respectively (e.g., Rathborne et al. 2006; Parsons et al. 2009; Vasyunina et al. 2009; Rygl et al. 2010; Miettinen & Harju 2010). Molecular spectral-line observations towards these clumps have shown the typical gas kinetic temperature to lie in the range T_kin ≈ 10 − 20 K (Carey et al. 1998; Teyssier et al. 2002; Sridharan et al. 2005; Pillai et al. 2006b; Sakai et al. 2008, hereafter SSK08; Zhang et al. 2011; Devine et al. 2011; Ragan et al. 2011). Spectral energy distributions of clumps within IRDCs yield dust temperature values of T_dust ≈ 10 − 50 K (e.g., Rathborne et al. 2010; Henning et al. 2010). The star-formation process within some of these sources manifests itself through embedded infrared emission (e.g., Chambers et al. 2009; Ragan et al. 2009; Henning et al. 2010), outflows (Beuther & Sridharan 2007; Fallscheer et al. 2009; Sanhueza et al. 2010), and H₂O and CH₃OH maser emission (Pillai et al. 2006a; Ellingsen 2006; Wang et al. 2006). Hot cores (Beuther et al. 2005a, 2006; Rathborne et al. 2007, 2008) and UC HII regions (Battersby et al. 2010) have also been discovered within IRDCs. These are clear signs of high-mass star-formation activity. On the other hand, some of the massive clumps are dark at mid-infared (MIR) and do not show any other evidence of ongoing star formation (e.g., G024.61-00.33 MM2, Rygl et al. 2010; and G028.23-00.19 GLM1, Battersby et al. 2010). These sources are candidates for being the high-mass prestellar or precluster objects and are promising targets for the studies of the initial conditions of massive star-formation.

In this paper we aim to further constrain the physical and chemical characteristics of massive clumps associated with IRDCs. We study the degrees of CO depletion and deuterium fractionation, and the fractional ionisation in a sample of seven clumps. So far, only a few studies have been carried out to investigate molecular deuteration in IRDCs (Pillai et al. 2007; 2010; Fontani et al. 2011; Pillai et al. 2011). The level of deuteration has the potential to be used as an evolutionary indicator (e.g., Crapsi et al. 2005; Emprechtinger et al. 2009; Fontani et al. 2011) and as such it is useful to study it in IRDCs. The ionisation degree, x(e), and the cosmic-ray ionisation rate of H₂, ζ_H2, are also essential knowledge towards an understanding of the initial conditions of star formation (e.g., Caselli et al. 1998). Fractional ionisation determines the coupling between the gas and magnetic field, and can thus play a role in the dynamical evolution of a star-forming object. In addition, the ionisation fraction has an influence on the gas-phase chemical reactions, and thus assessing the ionisation structure of a source helps to better understand its overall chemistry. Because estimates of x(e) and ζ_H2 in massive clumps within IRDCs have not been reported so far, constraining these parameters through observations seems particularly worthwhile.

This paper is organised as follows. The observations, data-reduction procedures, and the source sample are described in Sect. 2. The direct observational results are presented in Sect. 3. In Sect. 4 we describe the analysis and present the results of the physical and chemical properties of the clumps. In Sect. 5 we discuss our results, and in Sect. 6 we summarise the main conclusions of this study.

2. Source selection, observations, and data reduction

2.1. Source selection

For this study we selected six clumps from the previous studies by Beuther et al. (2002a), Sridharan et al. (2005), and Rathborne et al. (2006), and one clump, ISOSS J18364-0221 SMM1², previously studied by Birkmann et al. (2006) and Hennemann et al. (2009). We selected sources with a high mass inside a typical clump radius of  ≲0.5 pc (revised masses from the above reference studies are  ~60−360   M_⊙; see Sect. 4.1.3), and which are likely to represent different evolutionary stages of (high-mass) star formation: two of the clumps are clearly associated with Spitzer 8-μm point-sources and five appear dark in the Spitzer 8-μm images. Moreover, three of the clumps are associated with H₂O and/or CH₃OH masers. As will be discussed later, the youngest sources in our sample are likely to be G015.05 MM1 and G015.31 MM3, whereas I18102 MM1 appears to be in the most advanced stage of evolution. In this way, we were aiming to investigate if there are any variations or evolutionary trends in the clump physical and chemical properties among the sample.

To gather the source sample, one criterion was that the source distance is not very long (revised values are  ~2.5−3.5 kpc; see Sect. 4.1.1 and references therein) to achieve a reasonable spatial resolution with the single-dish telescope. For example, 20″−30″ corresponds to  ~0.24−0.51 pc at the cloud distances, i.e., the spatial resolution is comparable to the clump sizes. We selected sources for which dust continuum data at (sub)mm wavelengths are available to derive the H₂ column density needed in the derivation of molecular fractional abundances. Also, to derive reliable physical parameters of the sources, they were supposed to have known gas kinetic temperature or dust temperature. The source list is given in Table 1 and the sources are discussed in more detail in Appendix B.

2.2. APEX molecular-line observations

In this subsection we describe the molecular-line observations carried out with the Atacama Pathfinder Experiment (APEX) 12-m telescope at Llano de Chajnantor (Chile). The observed lines, their spectroscopic properties, and observational parameters are listed in Table 2.

2.2.1. C¹⁷O(2−1)

The C¹⁷O(2−1) observations towards all the sources except J18364 SMM1 were carried out on 25 March 2008 with APEX during the observational campaign for the Science Verification (SV) of the Swedish Heterodyne Facility Instrument (SHFI or SHeFI; Belitsky et al. 2007; Vassilev et al. 2008a). The source J18364 SMM1 was observed on 26 April 2009. As frontend we used APEX-1 of the SHFI (Vassilev et al. 2008b). APEX-1 operates in single-sideband (SSB) mode using sideband separation mixers, and it has a sideband rejection ratio  >0 dB. The centre frequency for the image band is given by ν_image = ν_cen ± 12 GHz, where ν_cen is the passband centre frequency, and the  ±  refers to the image band corresponding to the upper or lower sideband (USB or LSB), respectively. The backend was the Fast Fourier Transfrom Spectrometer (FFTS; Klein et al. 2006) with a 1 GHz bandwidth divided into 16 384 spectral channels during SV and 8192 channels during observations towards J18364 SMM1.

The observations were performed in the wobbler-switching (WS) mode with a 150″ azimuthal throw (symmetric offsets) and a chopping rate of 0.5 Hz. The telescope pointing and focusing were checked by continuum scans on the planet Jupiter and the pointing was found to be accurate to  ~3″. Calibration was made by the chopper-wheel technique, and the intensity scale given by the system is ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$, the antenna temperature corrected for atmospheric attenuation.

The spectra were reduced using the CLASS90 programme of the IRAM’s GILDAS software package³. The individual spectra were averaged and the resulting spectra were Hanning-smoothed to increase the signal-to-noise (S/N) ratio of the data. A first- or third-order polynomial was applied to correct the baseline in the resulting spectra. In one case (G015.05 MM1) a polynomial baseline of order six had to be applied because there appeared to be some emission in the wobbler off-beams, which resulted in a very poor baseline. The resulting 1σ rms noise levels are 40–60 mK at the smoothed resolutions. The final velocity resolution of the smoothed spectra is 0.16 km s^-1, except for the spectrum towards J18364 SMM1, where it is 0.32 km s^-1.

The J = 2 − 1 transition of C¹⁷O contains nine hyperfine (hf) components. We fitted the hf structure of the transition using “method hfs” of CLASS90 to derive the LSR velocity (v_LSR) of the emission, and FWHM linewidth (Δv). The hf-line fitting can also be used to derive the line optical thickness, τ. However, in all spectra the hf components are blended and thus the optical thickness cannot be reliably determined. For the rest frequencies of the hf components we used the values from Ladd et al. (1998, Table 6). The adopted central frequency, 224 714.199 MHz, is that of the J_F = 2_9/2 → 1_7/2 hf component, which has a relative intensity of ${\mathit{R}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}$.

2.2.2. H¹³CO⁺(3−2) and DCO⁺(3−2)

The H¹³CO⁺(3−2) and DCO⁺(3−2) observations were carried out on 9–11 April and 3 June 2010. As frontend we used SHFI/APEX-1 and the backend was the FFTS with a 1 GHz bandwidth divided into 8192 channels. The observations were performed in the WS mode with similar settings as for the C¹⁷O(2−1) observations. The telescope pointing and focusing were checked by observing the planets Saturn and Jupiter and the pointing was again found to be accurate to  ~3″.

The individual spectra were averaged, smoothed, and a first- or third-order polynomial was applied to correct the baseline. The resulting 1σ rms noise values are 20–80 mK for H¹³CO⁺(3−2) and 30–40 mK for DCO⁺(3−2). The final velocity resolutions of the smoothed H¹³CO⁺(3−2) and DCO⁺(3−2) spectra are 0.28 and 0.34 km s^-1, respectively. As shown in the last column of Table 2, the on-source integration times were very different between different sources. This is particularly the case for H¹³CO⁺(3−2), where t_int varies by a factor of  ~11. This explains the considerable variation (factor of 4) in the noise values of this transition.

The spectral-line transitions of H¹³CO⁺ and DCO⁺ also have an hf structure (e.g., Schmid-Burgk et al. 2004; Caselli & Dore 2005). The J = 3 − 2 transition we observed is split up into eight and five hf components in the case of H¹³CO⁺ and DCO⁺, respectively. To fit the hf structure, we used the rest frequencies of the H¹³CO⁺(3−2) and DCO⁺(3−2) lines from the Cologne Database for Molecular Spectroscpy (CDMS; Müller et al. 2005)⁴. The adopted central frequency of H¹³CO⁺(3−2), 260 255.35250 MHz, is that of the J_F1,F = 3_7/2,4 → 2_5/2,3 hf component, which has a relative intensity of R_i ≃ 0.322. We note that this frequency is 13.5 kHz higher than the value determined by Schmid-Burgk et al. (2004; Table 3). In the case of DCO⁺(3−2) the central frequency used, 216 112.57790 MHz, is that of the J_F = 3₄ → 3₃ hf component; its relative intensity is ${\mathit{R}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{7}}$. Caselli & Dore (2005; their Table 5) determined a 2.1 kHz higher frequency for DCO⁺(3−2). The frequency intervals of the hf components for H¹³CO⁺(3−2) and DCO⁺(3−2) are very small. Therefore, the lines overlap significantly, which causes the hf structure to be heavily blended.

2.2.3. N₂H⁺(3−2) and N₂D⁺(3−2)

The N₂H⁺(3−2) observations towards three of our target positions (I18102 MM1, G015.05 MM1, and G015.31 MM3) were carried out on 29 April 2009. The N₂D⁺(3−2) observations were performed during 16–17 and 19–20 October 2008, 26 and 28 April 2009, and 7 June 2010. As frontend we used SHFI/APEX-2 and APEX-1 for N₂H⁺(3−2) and N₂D⁺(3−2), respectively. The backend was the 1 GHz FFTS with 8 192 channels. The observations were performed in the WS mode with the similar settings as described above. The pointing and focusing were made by continuum scans on the planets Jupiter, Mars, and Neptune. The pointing was found to be accurate to  ~3″−5″.

The data were reduced as described above, and a first- or third-order polynomial baseline correction resulted in 1σ rms noise values of 40–100 mK and 20–50 mK in the case of N₂H⁺ and N₂D⁺, respectively. The final velocity resolutions of the smoothed N₂H⁺ and N₂D⁺ spectra are 0.26 and 0.32 km s^-1, respectively. Again, the on-source integration times were very different for different sources, particularly for N₂D⁺(3−2) (see Col. 10 of Table 2).

The J = 3−2 transitions of both N₂H⁺ and N₂D⁺ contain 38 hf components. The hf lines were fitted using the rest frequencies from Pagani et al. (2009b; Tables 4 and 10). The adopted central frequencies of N₂H⁺(3−2) and N₂D⁺(3−2), 279 511.832 and 231 321.912 MHz, are those of the J_F1F = 3₄₅ → 2₃₄ hf component, which has a relative intensity of ${\mathit{R}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{11}}{\mathrm{63}}$. In these cases the hf components are blended as well and consequently the value of τ cannot be reliably determined through hf fitting.

2.3. Archival dust continuum data

Submillimetre and/or millimetre dust continuum maps of our sources are published in the papers by Beuther et al. (2002a; MAMBO 1.2 mm), Williams et al. (2004; SCUBA 450 & 850 μm), Rathborne et al. (2006; MAMBO-II 1.2 mm), and Birkmann et al. (2006; SCUBA 450 & 850 μm). For the sources I18102 MM1, I18182 MM2, and J18364 SMM1 we retrieved 850 μm dust continuum data from the JCMT/SCUBA archive (Di Francesco et al. 2008)⁵. No archival SCUBA data were available for the rest of our sources. Thus, for the sources G015.05 MM1, G015.31 MM3, and I18223 MM3 we extracted Bolocam 1.1-mm maps from images produced by the Bolocam Galactic Plane Survey (BGPS; see, e.g., Bally et al. 2010, and references therein)⁶. For I18151 MM2, no archival (sub)mm data were found and thus we used the MAMBO 1.2-mm map of this source published in the paper by Beuther et al. (2002a; kindly provided by H. Beuther).

In the present paper the (sub)mm dust continuum data are used in particular to determine the H₂ column densities towards the line-observation positions. These are used to derive the fractional abundances of the observed molecules.

2.4. Spitzer infrared data

We also used Spitzer IRAC 8-μm and MIPS 24-μm images of the sources to investigate the clump associations with protostellar activity. The 8-μm images are from the Spitzer Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE; Benjamin et al. 2003)⁷. The 24-μm images (not shown in the paper) were taken from the Spitzer MIPS GALactic plane survey (MIPSGAL; Carey et al. 2009)⁸. In Fig. 1 we show the Spitzer/IRAC 8-μm images of the clumps overlaid with contours of (sub)mm dust continuum emission. We note that better quality Spitzer data of the sources are available and already published (see references in Appendix B). However, in the present work we did not use Spitzer data to derive the physical properties of the sources; the data are used to illustrate the source appearance at MIR wavelengths.

3. Observational results

3.1. Spectra

The Hanning-smoothed spectra are shown in Fig. 2. Note that all data are presented in units of ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$. The C¹⁷O(2−1) line is clearly detected towards all sources, but the hf structure of the line is completely blended in all cases except for G015.31 MM3, where it is partially resolved. The hf structure is not resolved because the linewidths are wider than the separation in velocity of individual hf components.

The H¹³CO⁺(3−2) line is also clearly detected towards all target positions, but towards G015.31 MM3 the line is very weak (only  ~3σ). The “double peak”-like profiles in some of the smoothed H¹³CO⁺(3−2) spectra can be attributed to the noise in most cases. In the case of J18364 SMM1, the self-absorption dip is a distinctive feature also in the unsmoothed spectrum, and thus likely to be real. DCO⁺(3−2) emission is seen towards all sources, but the line is very weak in I18102 MM1 and G015.31 MM3 (detected only at  ~3σ).

N₂H⁺(3−2) line-observations were performed only towards three of our sources, and the line is quite strong in all cases. The N₂H⁺(3−2) line of I18102 MM1 exhibits red asymmetry with a central dip near the velocity of the strongest hf component; this feature is also clearly visible in the unsmoothed spectrum. The N₂D⁺(3−2) line is detected towards all sources except G015.31 MM3, I18182 MM2, and J18364 SMM1, where the spectra have a low S/N ratio.

3.1.1. Other line detections

All observed sources show additional spectral lines in the observed frequency bands, mostly in the spectra covering the DCO⁺(3−2) and N₂D⁺(3−2) lines. The line identification was made with Weeds, which is an extension of CLASS (Maret et al. 2011), and applying the Jet Propulsion Laboratory⁹ (JPL; Pickett et al. 1998) and CDMS spectroscopic databases. We used the LTE modelling application of Weeds to check if all predicted lines of a candidate molecule were present in the observed spectrum (see Sect. 4.2). In some cases, we were able to reject some line candidates on the basis of non-detection of other transitions expected at nearby frequencies.

I18102 MM1 and I18151 MM2 show relatively line-rich spectra at  ~231 GHz and  ~216 GHz, in the USB and LSB, respectively. Two additional lines at  ~231 GHz (USB) are also visible in the spectrum of I18182 MM2 (see Fig. 3). The only unambiguous line identification in the spectrum of I18102 MM1 is ¹³CS(5−4). The second-strongest line (${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{=}\mathrm{0.12}$ K) can be assigned to C¹⁸O(2−1) at  ~219.56 GHz seen in the image side band (marked with “image”). The strongest line in the spectrum (${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{=}\mathrm{0.13}$ K) might be caused by CNCHO(J_Ka,Kc = 23_0,23 − 22_0,22) at  ~219.53 GHz arising from the image band, or it could represent an additional C¹⁸O(2−1) velocity component along the line of sight. The latter is likely to be true because of the relatively strong line intensity. The line at  ~231.06 GHz can be assigned to either OCS(19 − 18) or CH₃NH₂-E(J_Ka,Kc = 7_2,5 − 7_1,5) because these two transitions are blended. However, it could also be caused by yet another velocity component of C¹⁸O(2−1) seen towards I18102 MM1. From the spectrum towards I18151 MM2 we identified the C₂D(N = 3 − 2) line separated into two hf groups. This transition is split into 13 hf components, of which the detected two groups contain 10 components together. The hf group at  ~216.43 GHz is blended with CH₂CHC¹⁵N(J_Ka,Kc = 25_2,24 − 25_0,25). The J_Ka,Kc = 3_3,0 − 2_2,1 transition of ortho-c-C₃H₂ is seen towards I18151 MM2, but the line at  ~216.45 GHz is unidentified (marked with “U”). The CH₃COCH₃-EA(J_Ka,Kc = 15_9,7 − 14_8,6) line in the USB is possibly detected towards I18182 MM2. However, it is blended with C¹⁸O(2−1) entering via the image band. The relatively strong intensity of ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{=}\mathrm{0.17}$ K suggests the line is caused by C¹⁸O(2−1). The weaker nearby line could be related to an additional C¹⁸O(2−1) velocity component.

The J_Ka,Kc = 3_3,0 − 2_2,1 transition of ortho-c-C₃H₂ at  ~216.3 GHz in the LSB is also detected towards I18102 MM1, I18182 MM2, and I18223 MM3 (Fig. 4). A thermal (v = 0) SiO(6 − 5) line at  ~260.5 GHz in the USB is seen towards I18102 MM1, I18151 MM2, and I18223 MM3 (Fig. 5). In particular, the SiO(6 − 5) line of I18151 MM2 exhibits very broad wing emission. Finally, as in the cases of I18102 MM1 and I18182 MM2, we detected the J = 2 − 1 transition of C¹⁸O at  ~219.56 GHz in the image side band towards G015.05 MM1, I18151 MM2, G015.31 MM3, I18223 MM3, and J18364 SMM1 (Fig. 6). The intensities of these C¹⁸O lines, ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{~}\mathrm{0.09}\mathrm{-}\mathrm{0.31}$ K, are mostly comparable to those possibly detected towards I18102 MM1 and I18182 MM2. Note that the lines “leaking” from the rejected image sideband are heavily attenuated by the sideband filter. Thus, we cannot establish the correct intensity scale for the image-band lines.

In Table 3 we list the observed extra transitions, their rest frequencies (observed frequency for the U-line), and upper-state energies. The rest frequencies and upper-state energies were assigned using the JPL, CDMS, and Splatalogue¹⁰ (Remijan et al. 2007) spectroscopic databases. The additional lines/species are discussed in more detail in Appendix A.

3.2. Spectral line parameters

The spectral line parameters are given in Table 4. In this table we give the radial velocity (v_LSR), FWHM linewidth (Δv), peak intensity (${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$), integrated intensity (${}^{\mathrm{\int }}\mathit{T}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{d}\mathit{v}$), peak optical depth of the line (τ₀), and excitation temperature (T_ex). The values of v_LSR and Δv for C¹⁷O, H¹³CO⁺, DCO⁺, N₂H⁺, N₂D⁺, and C₂D were derived through fitting the hf structure. For the other lines these parameters were derived by fitting a single Gaussian to the line profile. The values of ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$ were also determined by fitting the lines with a single Gaussian. The velocity-integrated intensities ${}^{\mathrm{\int }}\mathit{T}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{d}\mathit{v}$ were calculated over the velocity range given in Col. (6) of Table 4. The uncertainties reported in v_LSR, Δv, ${\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}$, and ${}^{\mathrm{\int }}\mathit{T}_{\mathrm{A}}^{\mathrm{\ast }}\mathrm{d}\mathit{v}$ represent the formal 1σ errors determined by the fitting routine. We used RADEX¹¹ (van der Tak et al. 2007) to determine the values of τ₀ and T_ex for the lines of C¹⁷O, H¹³CO⁺, DCO⁺, N₂H⁺, c-C₃H₂, OCS, ¹³CS, and SiO. RADEX modelling is described in more detail in Sect. 4.2, where we discuss the determination of molecular column densities. For N₂D⁺(3−2) we used as T_ex the value derived for N₂H⁺(3−2) or, when not possible, we adopted T_ex [N₂H⁺(1 − 0)]  from SSK08. We note that the T_ex values obtained for the J = 3 − 2 transition of N₂H⁺ are comparable to those of J = 1 − 0 from SSK08.

4. Analysis and results

4.1. Revision of clump properties

In this subsection, we present recalculations of several clump properties presented previously in the literature (see Tables 1 and 5).

4.1.1. Kinematic distances

Sakai et al. (2008) calculated the near kinematic distance for all of our clumps except J18364 SMM1 (see their Table 4). They used the LSR velocity of N₂H⁺(1−0) and the rotation curve of Clemens (1985) with the standard rotation parameters (Θ₀, R₀) = (220   km   s^-1, 8.5 kpc), where Θ₀ is the circular orbital speed of the Sun, and R₀ is the galactocentric distance of the Sun. Because N₂H⁺ is a tracer of high-density gas, the radial velocity derived from N₂H⁺ lines is suitable to determine the source kinematic distance. The hf structure of N₂H⁺(1−0) was also resolved (to some degree) towards many of the sources by SSK08, and thus the hf-structure fitting is expected to yield a reliable centroid velocity. Towards the three sources for which we detected N₂H⁺(3−2) lines, the obtained values of v_LSR are comparable to those from N₂H⁺(1 − 0). We recalculated the distances from SSK08 using the recent rotation curve of Reid et al. (2009), which is based on measurements of trigonometric parallaxes and proper motions of masers in high-mass star-forming regions. The best-fit rotation parameters of Reid et al. (2009) are (Θ₀, R₀) = (254   km   s^-1, 8.4 kpc). The resulting galactocentric distances and near kinematic distances (R_GC and d) are given in Cols. (6) and (7) of Table 1. The revised distances differ at most by 0.2 kpc from those reported by SSK08.

The source J18364 SMM1 is an exception. Its CO(1 − 0) radial velocity of  ~33 km s^-1 corresponds to a kinematic distance of  ~2.2 kpc according to the Brand & Blitz (1993) rotation curve (see Birkmann et al. 2006). We recomputed this distance from the LSR velocity of C¹⁷O(2−1) (34.8 km s^-1) because it is tracing higher density gas than the main CO isotopologue (similar radial velocity was obtained for the other detected transitions). The revised distance is  ~2.5 kpc according to the Reid et al. (2009) rotation curve; the previous value of 2.2 kpc is within the errors.

We note that for the sources I18102 MM1, I18151 MM2, I18182 MM2, and I18223 MM3 there are kinematic distance estimates (2.6, 3.0, 4.5, and 3.7 kpc, respectively) based on the CS(2−1) velocity and the rotation curve of Brand & Blitz (1993) (see Sridharan et al. 2002). For both G015.05 MM1 and G015.31 MM3, Rathborne et al. (2006) used the kinematic distance of 3.2 kpc derived from the ¹³CO(1 − 0) velocity and the rotation curve of Clemens (1985). Our distances differ by 0.1–1.0 kpc from the above values.

4.1.2. Gas and dust temperatures

Sridharan et al. (2005) determined the NH₃ rotational temperature, T_rot(NH₃), for three of our sources (I18151 MM2, I18182 MM2, and I18223 MM3) at the 40″ angular resolution. Sakai et al. (2008) have determined T_rot(NH₃) for all of our sources except J18364 SMM1 (see their Table 10). The values of T_rot(NH₃) from the above two studies are otherwise similar except in the case of I18223 MM3. Sridharan et al. reported the value 32.7 K for this source, whereas SSK08 obtained a temperature lower by a factor of two ($\mathrm{16.}{\mathrm{2}}_{-0.9}^{\mathrm{+}\mathrm{1.0}}$ K) at an angular resolution about 1.8 times poorer (73″). The T_rot(NH₃) value for J18364 SMM1 at the 40″ resolution is 10.75 K (Krause 2003). We converted T_rot(NH₃) from Krause (2003) and SSK08 into the gas kinetic temperature, T_kin, using the relationship given by Tafalla et al. (2004, Appendix B). The uncertainty in T_kin was calculated by propagating the larger of the two T_rot-errors given by SSK08.

For three of our sources, G015.05 MM1, I18223 MM3, and J18364 SMM1, dust temperature estimates have been made by Rathborne et al. (2010), Beuther et al. (2010), and Birkmann et al. (2006), respectively. The determined dust temperature values of G015.05 MM1, 11.0–36.0 K, bracket the T_kin value of 17.2 K. A dust temperature of I18223 MM3 was also derived by Beuther & Steinacker (2007) from the source SED using the Spitzer/MIPS data at 24 and 70 μm, and MAMBO 1.2 mm and PdBI 3.2 mm data. The value 15 K they obtained for the cold part of the spectrum is three K lower than the value 18 K derived recently by Beuther et al. (2010). The latter also used the Herschel/PACS (70, 100, and 160 μm) and SPIRE (250, 350, and 500 μm) data, and SCUBA 850 μm data to construct the source SED (but not the 3.2-mm flux density). The dust temperature of I18223 MM3 is comparable to its gas temperature of 18.7 K. By utilising the far-infrared and submm flux density ratios, Birkmann et al. (2006) deduced the value ${\mathit{T}}_{\mathrm{dust}}\mathrm{=}\mathrm{16.}{\mathrm{5}}_{-3.0}^{\mathrm{+}\mathrm{6.0}}$ K for J18364 SMM1, which is higher than the gas temperature 11.4 K.

The gas and dust temperatures of the clumps are given in Cols. (2) and (3) of Table 5.

4.1.3. Clump masses, radii, and H₂ column and number densities

We calculated the masses and beam-averaged H₂ column densities of the clumps using the formulas $\mathit{M}\mathrm{=}\frac{{\mathit{S}}_{\mathit{\nu }}{\mathit{d}}^{\mathrm{2}}}{{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{R}}_{\mathrm{d}}}\mathit{,}$(1)$\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}}{{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{m}}_{\mathrm{H}}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{R}}_{\mathrm{d}}}\mathrm{·}$(2)In the above formulas, ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}$ and S_ν are the peak surface brightness and integrated flux density of the (sub)mm dust emission; note that ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{=}{\mathit{S}}_{\mathit{\nu }}^{\mathrm{peak}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{beam}}$, where ${\mathit{S}}_{\mathit{\nu }}^{\mathrm{peak}}$ is the peak flux density and Ω_beam is the solid angle of the telescope beam. B_ν(T_dust) is the Planck function for a dust temperature T_dust, κ_ν is the dust opacity per unit dust mass, R_d ≡  ⟨ M_dust/M_gas ⟩  is the average dust-to-gas mass ratio, μ_H2 is the mean molecular weight per H₂ molecule (2.8 for the He/H abundance ratio of 0.1), and m_H is the mass of a hydrogen atom.

The ${\mathit{I}}_{\mathit{\nu }}^{\mathrm{dust}}$ values were determined from the (sub)mm maps shown as contours in Fig. 1. The angular resolutions of the maps are different: the beam size of the SCUBA 850 μm, Bolocam 1.1 mm, and MAMBO 1.2 mm data are 144, 31″, and 11″, respectively. The values of S_ν were taken from the literature as follows: for I18102 MM1, I18151 MM2, I18182 MM2, and I18223 MM3 we used the MAMBO 1.2-mm flux densities from Beuther et al. (2002a); for G015.05 MM1 and G015.31 MM3 we used the MAMBO-II 1.2-mm flux densities from Rathborne et al. (2006); for J18364 SMM1 we used the SCUBA 850-μm flux density from Birkmann et al. (2006). These flux densities are given in Col. (4) of Table 5.

As the dust temperature of the clumps we used the gas kinetic temperatures listed in Col. (2) of Table 5 and assumed that T_dust = T_kin. We extrapolated the values of κ_ν from the Ossenkopf & Henning (1994) model for dust grains with thin ice mantles, coagulated for 10⁵ yr at a gas density of n_H = n(H) + 2n(H₂) ≃ 2n(H₂) = 10⁶ cm^-3 [their Table 1, Col. (6)]. This is expected to be a reasonable dust model for the sources within cold and dense IRDCs. At the wavelengths considered in the present work, these values are κ_{850   μm} = 0.197 m² kg^-1, κ_1.1   mm = 0.121 m² kg^-1, and κ_1.2   mm = 0.106 m² kg^-1. The canonical value 1/100 was adopted for R_d. We note that the dust opacities are likely to be uncertain by a factor of  ≳ 2 (Ossenkopf & Henning 1994; Motte & André 2001). For comparison, Williams et al. (2004), Enoch et al. (2006), and Rygl et al. (2010) used the same dust model as we did and interpolated the values of κ_ν at 850 μm, 1.1 mm, and 1.2 mm to be 0.154, 0.114, and 0.1 m² kg^-1, respectively. These are slightly lower than our values, and the difference is likely to be caused by different extrapolation methods (e.g., log-interpolation) and/or different dust emissivity index, β, used by the authors to determine the value of κ_ν ∝ ν^β.

In the papers by Beuther et al. (2002a) and Rathborne et al. (2006), the reported clump sizes refer to their FWHM sizes (diameters) resulting from two-dimensional Gaussian fits. However, the integrated flux densities used to calculate the masses refer to larger clump areas, A. Accordingly, to calculate properly the volume-average H₂ number density,  ⟨ n(H₂) ⟩ , one way is to use the so-called effective radius of the clump, ${\mathit{R}}_{\mathrm{eff}}\mathrm{=}\sqrt{\mathit{A}\mathit{/}\mathit{\pi }}$ (see, e.g., Miettinen & Harju 2010). An alternative way to calculate  ⟨ n(H₂) ⟩  is to compute the amount of mass within the FWHM contour and use the FWHM radius. For sources with a Gaussian shape, the mass within the FWHM contour is a fraction ln   2 ≃ 0.693 of the total mass (see Kauffmann & Pillai 2010). Following Kauffmann & Pillai (2010), we reduce the clump masses by the above factor and use half the FWHM size as the effective radius of the clump. The density  ⟨ n(H₂) ⟩  was then calculated using the formula $\mathrm{⟨}\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\frac{\mathrm{⟨}\mathit{\rho }\mathrm{⟩}}{{\mathit{\mu }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{m}}_{\mathrm{H}}}\mathit{,}$(3)where  ⟨ ρ ⟩  = M/${\mathrm{(}}^{}\mathrm{4}\mathit{/}\mathrm{3}\mathit{\pi }{{\mathit{R}}_{\mathrm{eff}}^{\mathrm{3}}}^{\mathrm{)}}$ is the mass density, and M refers to the reduced clump mass as described above. For the source J18364 SMM1 we scaled the effective radius  ≈ 0.2 pc reported by Birkmann et al. (2006) using the revised distance (R_eff = 0.23 pc), and used the total mass within the clump area.

The results of the above calculations are presented in Cols. (5)–(8) of Table 5. The uncertainties in the derived quantities were propagated from the uncertainties in d and T_kin, but we note that the uncertainty in (sub)mm dust opacity is the largest source of error in M and N(H₂).

The derived values of M, N(H₂), and  ⟨ n(H₂) ⟩  are mostly within a factor of  ≲ 2 of those derived by Beuther et al. (2002a; where the revised equations by Beuther et al. 2005b, are employed), Rathborne et al. (2006), Birkmann et al. (2006), and Hennemann et al. (2009). We note that besides the different source distances, the dust parameters used by the authors in the above reference studies were different than here. Beuther et al. (2002a, 2005b) used T_dust values (35–50 K) resulting from the cold part of the source SEDs; the SEDs were constructed from the IRAS and mm data of the main core in the source region, which explains the rather high dust temperatures. In addition, Beuther et al. (2002a, 2005b) used the grain radius, mass density, and emissivity index values of a = 0.1 μm, ρ_dust = 3 g cm^-3, and β = 2, respectively, in the calculation of mass and column density. Following Hildebrand (1983), these values correspond to κ_1.2   mm ≈ 0.5 m² kg^-1. Rathborne et al. (2006) used the values T_dust = 15 K and κ_1.2   mm = 0.1 m² kg^-1. For the source J18364 SMM1, Birkmann et al. (2006) used ${\mathit{T}}_{\mathrm{dust}}\mathrm{=}\mathrm{16.}{\mathrm{5}}_{-3.0}^{\mathrm{+}\mathrm{6.0}}$ K instead of T_kin, and they had κ_{850   μm} = 0.180 m² kg^-1, which is slightly smaller than ours even though the adopted dust model was the same.

4.2. Molecular column densities and fractional abundances

The beam-averaged column densities of C¹⁷O, H¹³CO⁺, DCO⁺, N₂H⁺, SiO, OCS, ¹³CS, and o-c-C₃H₂ were derived using a one-dimensional spherically symmetric non-LTE radiative transfer code called RADEX (see Sect. 3.2). RADEX uses the method of mean escape probability for an isothermal and homogeneous medium. The molecular data files (collisional rates) used in the RADEX excitation analysis were taken from the LAMDA database (Schöier et al. 2005). The C¹⁷O, H¹³CO⁺, DCO⁺, and N₂H⁺ transitions are treated as a hypothetical unsplit transition. The input parameters in the off-line mode of RADEX are the gas kinetic temperature, H₂ number density, and the width (FWHM) and intensity of the spectral line. We used the values of T_kin and  ⟨n(H₂)⟩  listed in Table 5. However, we multiplied the densities by 1.2 (He/H₂ = 0.2) to take the collisions with He into account (see Sect. 4.1 of the RADEX manual¹²). As the input line intensity we used the main-beam brightness temperature ${\mathit{T}}_{\mathrm{MB}}\mathrm{=}{\mathit{T}}_{\mathrm{A}}^{\mathrm{\ast }}\mathit{/}{\mathit{\eta }}_{\mathrm{MB}}$. When the source is resolved, T_MB is equal to the Rayleigh-Jeans equivalent radiation temperature, T_R. The simulations aim to reproduce the observed line intensity and yield the values of the line peak optical thickness (τ₀), excitation temperature (T_ex), and the total column density of the molecule (N_tot). We varied T_kin and  ⟨ n(H₂) ⟩  according to their errors to estimate the uncertainties associated with τ₀, T_ex, and N_tot. The lines appear to be optically thin in most cases. The optical thickness is τ₀ ≳ 1 for all N₂H⁺ lines, for H¹³CO⁺ towards I18102 MM1 and J18364 SMM1, and for c-C₃H₂ of I18182 MM2. The T_ex values for C¹⁷O are close to T_kin, indicating that the lines are nearly thermalised. Also, T_ex(H¹³CO⁺) is found to be equal to T_ex(DCO⁺) within the errors. We note that the column density determination is in some cases hampered by self-absorbed line profiles.

For the rest of the observed molecules there are no molecular data files available in the LAMDA database. The line optical thicknesses and total beam-averaged column densities of these molecules were determined through LTE modelling with CLASS/Weeds. The input parameters for a Weeds model are N_tot, T_ex, source size (θ_s), linewidth (FWHM), and offset from the reference-channel velocity. The linewidth is directly determined from the observed line profile, so there are basically three free parameters left (N_tot, T_ex, θ_s). Some of the model parameters may be degenerate, and cannot be determined independently (Schilke et al. 2006; Maret et al. 2011). The source size is degenerate with excitation temperature for completely optically thick lines (τ ≫ 1), and with column density if the lines are completely optically thin (τ ≪ 1). We assumed that the source fills the telescope beam, i.e., that the source solid angle, ${\mathrm{\Omega }}_{\mathrm{s}}\mathrm{\propto }{\mathit{\theta }}_{\mathrm{s}}^{\mathrm{2}}$, equals the beam solid angle, ${\mathrm{\Omega }}_{\mathrm{A}}\mathrm{\propto }{\mathit{\theta }}_{\mathrm{HPBW}}^{\mathrm{2}}$, or θ_s = θ_HPBW. When the beam filling factor is unity, the line brightness temperature is T_B ≃ T_MB (see Eqs. (1) and (2) in Maret et al. 2011). For N₂D⁺, we used as T_ex the values obtained for N₂H⁺ from RADEX simulations or from SSK08 (Sect. 3.2). For the rest of the transitions we adopted the value T_ex = E_u/k_B for linear rotors, and T_ex = 2/3 × E_u/k_B for asymmetric top rotors. The adopted rotational excitation temperatures give a lower limit to the total beam-averaged column densities (Hatchell et al. 1998). The input N_tot was then varied until a reasonable fit to the line was obtained (see Fig. 7). We note that we detected only one transition per species for the complex organics. We therefore cannot apply a rotation diagram to derive rotational temperatures and molecular column densities (e.g., Goldsmith & Langer 1999).

We also determined the HCO⁺ column density from the column density of H¹³CO⁺. For this calculation, it was assumed that the carbon-isotope ratio  [¹²C] / [¹³C]  depends on R_GC according to the relationship given by Wilson & Rood (1994): $\frac{{\mathrm{[}}^{\mathrm{12}}\mathrm{C}\mathrm{]}}{{\mathrm{[}}^{\mathrm{13}}\mathrm{C}\mathrm{]}}\mathrm{=}\mathrm{7.5}\mathrm{\times }{\mathit{R}}_{\mathrm{GC}}\mathrm{\left[}\mathrm{kpc}\mathrm{\right]}\mathrm{+}\mathrm{7.6.}$(4)For the R_GC values considered here, 5.1−6.3 kpc, the above ratio lies in the range  ~46−55¹³.

For I18151 MM2, I18182 MM2, and I18223 MM3 we do not have N₂H⁺ data. For these sources we computed the N₂H⁺ column density from the J = 1−0 line parameters (T_ex and τ) from SSK08 (see, e.g., Eq. (10) in Miettinen et al. 2010). We obtain about 1.1–1.2 times higher N₂H⁺ column densities compared to SSK08, who used the optically thin approximation (see Sect. 5.1 for a more detailed discussion).

We calculated the fractional abundances of the molecules by dividing the molecular column density by the H₂ column density: x(mol) = N(mol)/N(H₂). For this purpose, the values of N(H₂) were derived from the (sub)mm dust continuum maps smoothed to the corresponding resolution of the line observations. With a resolution of 31″, the Bolocam 1.1-mm data could not be smoothed to correspond the resolution of the line observations. In these cases we used the original Bolocam data; the 31″ resolution is in most cases comparable to that of the line observations (), and therefore we do not expect this to be a significant source of error. The derived column densities and abundances are listed in Tables 6 and 7. The abundance errors were derived by propagating the errors in N(mol) and N(H₂).

4.3. CO depletion and deuterium fractionation

To estimate the amount of CO depletion in the clumps, we calculated the CO depletion factor, f_D, following the analysis presented in the paper by Fontani et al. (2006). If x(CO)_can is the “canonical” (undepleted) abundance, and x(CO)_obs is the observed CO abundance, f_D is given by ${\mathit{f}}_{\mathrm{D}}\mathrm{=}\frac{\mathit{x}\mathrm{(}\mathrm{CO}{\mathrm{)}}_{\mathrm{can}}}{\mathit{x}\mathrm{(}\mathrm{CO}{\mathrm{)}}_{\mathrm{obs}}}\mathrm{·}$(5)The “canonical” CO abundance at the galactocentric distance R_GC was calculated using the relationship (Eq. (7) in Fontani et al. 2006) $\mathit{x}\mathrm{(}\mathrm{CO}{\mathrm{)}}_{\mathrm{can}}\mathrm{=}\mathrm{9.5}\mathrm{\times }{\mathrm{10}}^{-5}{\mathrm{e}}^{\mathrm{1.105}\mathrm{-}\mathrm{0.13}{\mathit{R}}_{\mathrm{GC}}\mathrm{\left[}\mathrm{kpc}\mathrm{\right]}}\mathit{.}$(6)This relationship results from using the value R₀ = 8.5 kpc, whereas our R_GC values are computed using R₀ = 8.4 kpc. This small discrepancy is negligible, however. At R_GC = 8.5 kpc, the above relationship gives the standard value 9.5 × 10^-5 for the abundance of the main CO isotopologue in the solar neighbourhood (Frerking et al. 1982). To calculate the “canonical” C¹⁷O abundance we take into account that the oxygen-isotopic ratio  [¹⁶O] / [¹⁸O]  depends on R_GC according to the relationship (Wilson & Rood 1994) $\frac{{\mathrm{[}}^{\mathrm{16}}\mathrm{O}\mathrm{]}}{{\mathrm{[}}^{\mathrm{18}}\mathrm{O}\mathrm{]}}\mathrm{=}\mathrm{58.8}\mathrm{\times }{\mathit{R}}_{\mathrm{GC}}\mathrm{\left[}\mathrm{kpc}\mathrm{\right]}\mathrm{+}\mathrm{37.1.}$(7)For the R_GC values of our clumps the above ratio ranges from about 337 to 407.5. When this relationship is combined with the  [¹⁸O] / [¹⁷O]  ratio, for which we use the standard value 3.52 (Frerking et al. 1982), the value of x(C¹⁷O)_can can be calculated as $\begin{array}{ccc}\mathit{x}\mathrm{(}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}{\mathrm{)}}_{\mathrm{can}}& \mathrm{=}& \frac{\mathit{x}\mathrm{(}\mathrm{CO}{\mathrm{)}}_{\mathrm{can}}}{{\mathrm{[}}^{\mathrm{18}}\mathrm{O}\mathrm{]}\mathit{/}{\mathrm{[}}^{\mathrm{17}}\mathrm{O}\mathrm{]}\mathrm{\times }{\mathrm{[}}^{\mathrm{16}}\mathrm{O}\mathrm{]}\mathit{/}{\mathrm{[}}^{\mathrm{18}}\mathrm{O}\mathrm{]}}\\ & \mathrm{=}& \frac{\mathit{x}\mathrm{(}\mathrm{CO}{\mathrm{)}}_{\mathrm{can}}}{\mathrm{3.52}\mathrm{\times }\mathrm{(}\mathrm{58.8}\mathrm{\times }{\mathit{R}}_{\mathrm{GC}}\mathrm{[}\mathrm{kpc}\mathrm{]}\mathrm{+}\mathrm{37.1}\mathrm{)}}\mathrm{·}\end{array}$(8)The depletion factor f_D is then calculated from f_D = x(C¹⁷O)_can/x(C¹⁷O)_obs, and the results, f_D = 0.6    ±    0.1 − 2.7 ± 1.3, are listed in Col. (2) of Table 9. The  ± -error quoted was calculated by propagating the uncertainty in x(C¹⁷O)_obs. The low values of f_D indicate that CO is not significantly depleted, if at all, in our clumps.

The degree of deuterium fractionation in HCO⁺ and N₂H⁺ was calculated by dividing the column density of the deuterated isotopologue by its normal hydrogen-bearing form as R_D(HCO⁺) ≡ N(DCO⁺)/N(HCO⁺) and R_D(N₂H⁺) ≡ N(N₂D⁺)/N(N₂H⁺). The values of R_D(HCO⁺) and R_D(N₂H⁺) are in the range 0.0002 ± 0.0001 − 0.014 ± 0.003 and 0.002 ± 0.001 − 0.028 ± 0.005, respectively; the error in R_D was derived from the errors in the corresponding column densities (see Cols. (3) and (4) of Table 9).

4.4. Ionisation degree

We have derived fractional abundances of different ionic species and their different isotopologues. In principle, a lower limit to the ionisation degree can be obtained by simply summing up these abundances (e.g., Caselli et al. 2002; Miettinenet al. 2009): $\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathit{>}\mathit{x}\mathrm{\left(}\mathrm{HC}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}{\mathrm{H}}^{\mathrm{13}}\mathrm{C}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}\mathrm{DC}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}{\mathrm{N}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathrm{+}\mathit{x}\mathrm{\left(}{\mathrm{N}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{\right)}\mathit{.}$(9)This is based on the gas quasi-neutrality: the electron abundance equals the difference between the total abundances of the cations and the anions. The resulting values, x(e) >  ∑ _ix(ions)_i, are listed in Col. (5) of Table 9. These lower limits are in the range 7 × 10^-9 − 1.2 × 10^-7.

The fractional ionisation in dense molecular clouds can also be inferred by employing ion-molecule chemical reaction schemes with the observed molecular abundances. All reactions considered here and notes on the associated rate coefficients are given in Table 8.

At first, we derived another estimate for a lower limit to x(e) through the ionisation balance determined by H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$, HCO⁺, N₂H⁺, and electrons following Qi et al. (2003). The H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and N₂H⁺ are mainly destroyed by CO; at steady state, reactions 3, 11, and 15 in Table 8 lead to the following equation for the lower limit to x(e) (cf. Eq. (6) in Qi et al. 2003): $\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathrm{\ge }\frac{{\mathit{k}}_{\mathrm{15}}\mathit{x}\mathrm{\left(}{\mathrm{N}}_{\mathrm{2}}{\mathrm{H}}^{\mathrm{+}}\mathrm{\right)}\mathit{x}\mathrm{\left(}\mathrm{CO}\mathrm{\right)}}{{\mathit{\beta }}_{\mathrm{11}}\mathit{x}\mathrm{\left(}\mathrm{HC}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}}\mathrm{·}$(10)The rate coefficients k₁₅ and β₁₁ were adopted from the UMIST database¹⁴ (Woodall et al. 2007). The derived values are listed in Col. (6) of Table 9. For I18151 MM2 and I18223 MM2 the x(e) values derived from Eq. (10) are similar to those computed from Eq. (9). For I18102 MM1, G015.31 MM3, and I18182 MM2 the summed abundance of different ionic species is clearly higher (by factors  ~ 3.5 − 16.3) than the lower limit to x(e) resulting from Eq. (10). In the case of G015.05 MM1, on the other hand, the summed ionic abundance is about six times lower. These discrepancies are not surprising because the chemical scheme behind Eq. (10) is certainly oversimplified.

Next, we determined the degree of ionisation by utilising the abundance ratios R_D(HCO⁺) ≡  [DCO⁺] / [HCO⁺]  and R_H ≡  [HCO⁺] / [CO] . The first studies of fractional ionisation based on the above ratios were carried out more than three decades ago (e.g., Guélin et al. 1977; Watson et al. 1978; Wootten et al. 1979). A similar analysis was subsequently applied in the papers by Caselli et al. (1998), Williams et al. (1998), Bergin et al. (1999), Anderson et al. (1999), and Caselli (2002).

We note that the following analysis includes several assumptions: i) HCO⁺ is mainly produced in the reaction between H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and CO (reaction 3 in Table 8); ii) all deuteration is caused by the reaction between H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and HD (reaction 2); iii) the presence of atomic deuterium, which could (slightly) increase the deuteration degree, is ignored; iv) ionic species are destroyed mainly by electrons, the most important neutrals (CO and O), and negatively charged dust grains; v) except CO and O, we neglect the contribution of some other neutral species, such as N₂, O₂, and H₂O, in the destruction of H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and H₂D⁺; and vi) we ignore the effects of refractory metals (Anderson et al. 1999; Caselli et al. 2002). Concerning part v), the abundances of neutrals, such as N₂ and H₂O, are poorly known and/or low. For instance, the first results from Herschel have shown that H₂O abundance is relatively low in high-mass star-forming regions (~10^-10 − 10^-8; van der Tak et al. 2010; Marseille et al. 2010a; Chavarría et al. 2010). Bergin et al. (1999) varied the nitrogen abundance in their chemical model and found that it does not affect the electron abundance. Also, the destructive reaction with O₂ is very slow (k = 9.3 × 10^-10 cm³ s^-1; UMIST).

By writing steady-state equations for the abundances of H₂D⁺ (reactions 2 and 7–10), DCO⁺ (reactions 7 and 13–14), and HCO⁺ (reactions 3, 7, and 11–12), it can be shown that ${\mathit{R}}_{\mathrm{D}}\mathrm{\left(}\mathrm{HC}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{{\mathit{k}}_{\mathrm{+}}\mathit{x}\mathrm{\left(}\mathrm{HD}\mathrm{\right)}}{{\mathit{k}}_{\mathrm{7}}\mathit{x}\mathrm{\left(}\mathrm{CO}\mathrm{\right)}\mathrm{+}{\mathit{\beta }}_{\mathrm{8}}\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{9}}\mathit{x}\mathrm{\left(}\mathrm{O}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{10}}\mathit{x}\mathrm{\left(}\mathrm{g}\mathrm{\right)}}\mathrm{·}$(11)Solving x(e) from the above formula yields $\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{8}}}\left[\frac{{\mathit{k}}_{\mathrm{+}}\mathit{x}\mathrm{\left(}\mathrm{HD}\mathrm{\right)}}{\mathrm{3}{\mathit{R}}_{\mathrm{D}}\mathrm{\left(}\mathrm{HC}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}}\mathrm{-}{\mathit{k}}_{\mathrm{7}}\mathit{x}\mathrm{\left(}\mathrm{CO}\mathrm{\right)}\mathrm{-}{\mathit{k}}_{\mathrm{9}}\mathit{x}\mathrm{\left(}\mathrm{O}\mathrm{\right)}\mathrm{-}{\mathit{k}}_{\mathrm{10}}\mathit{x}\mathrm{\left(}\mathrm{g}\mathrm{\right)}\right]\mathit{.}$(12)This represents the upper limit to x(e). We have assumed that k₇ = 1/3 × k₃, because H₂D⁺ can transfer a proton to CO producing HCO⁺; the rate for this is twice as high as for the channel producing DCO⁺. For the reaction rate of the deuteration reaction H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}\mathrm{+}\mathrm{HD}$ (k ₊ ) we used the recent results by Hugo et al. (2009, their Table VIII). For example, at 10 K the deuteration of H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ proceeds about 4.4 times faster according to Hugo et al. (2009) than suggested by laboratory measurements of Gerlich et al. (2002). Such a high rate was also used in the studies cited above, i.e., k ₊  = 1.5 × 10^-9 cm³ s^-1 (e.g., Caselli et al. 1998). We take the HD abundance to be twice the elemental D/H-ratio, i.e., x(HD) = 2 ×  [D] / [H]  ~ 3 × 10^-5 (e.g., Linsky et al. 2006; Prodanović et al. 2010). Because CO is not found to be depleted in our clumps, we assume that this is also the case for atomic oxygen and use the “standard” abundance relative to H₂ of x(O) = 3.52 × 10^-4, i.e., comparable to x(CO) (see Caselli et al. 1998; Caselli 2002). For the grain abundance, x(g), we use the value 2.64 × 10^-12, which is based on the values a = 0.1 μm, ρ_grain = 3 g cm^-3, and R_d = 1/100 (see, e.g., Eq. (15) in Pagani et al. 2009a). The derived upper limits to x(e) are shown in Col. (7) of Table 9. The values lie in the range  ~ 2 × 10^-6 −     ~ 3 × 10^-4, i.e., much higher than the lower limits estimated above.

By deriving a steady-state equation for the H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ abundance (see reactions 1–6 in Table 8) and applying it in the corresponding equation for HCO⁺, it can be shown that ${\mathit{R}}_{\mathrm{H}}\mathrm{=}\frac{\mathrm{[}{\mathit{\zeta }}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{/}\mathit{n}\mathrm{(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}{\mathit{k}}_{\mathrm{3}}}{\left[{\mathit{\beta }}_{\mathrm{4}}\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathit{x}\mathrm{\left(}\mathrm{CO}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{5}}\mathit{x}\mathrm{\left(}\mathrm{O}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{6}}\mathit{x}\mathrm{\left(}\mathrm{g}\mathrm{\right)}\right]\left[{\mathit{\beta }}_{\mathrm{11}}\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathrm{+}{\mathit{k}}_{\mathrm{12}}\mathit{x}\mathrm{\left(}\mathrm{g}\mathrm{\right)}\right]}\mathrm{·}$(13)After H$\begin{array}{}\mathrm{+}\\ \mathrm{2}\end{array}$ has formed via cosmic-ray ionisation of H₂ (reaction 1), it quickly reacts with H₂ to form H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ (Solomon & Werner 1971). Thus the H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ abundance is governed by the rate ζ_H2. When the fractional ionisation in the source is determined, Eq. (13) can be used to infer the cosmic-ray ionisation rate of H₂. To calculate ζ_H2, for each source we adopted as x(e) the summed abundance of ionic species. Using the derived upper limits to x(e) would result in unrealistically high values of ζ_H2. This suggests that the upper x(e) limits are clearly higher than the true values. The obtained results are shown in the last column of Table 9. In most cases, the ζ_H2 values lie in the range  ~1 × 10^-17−5 × 10^-16 s^-1. However, towards I18102 MM1 we obtain a very high rate of  ~1 × 10^-15 s^-1. Because I18102 MM1 is associated with a bright MIR point-source and shows no signs of CO depletion, it is probably a rapidly evolving source. Therefore, the steady-state assumption used in the above analysis may be invalid for I18102 MM1.

5. Discussion

In the following subsections we discuss our results and compare them with the results from previous studies.

5.1. Molecular column densities and abundances

Some of the molecular column densities and fractional abundances derived here have been determined in previous studies. Assuming LTE conditions and optically thin emission, SSK08 estimated N₂H⁺ column densities from the J = 1−0 transition for I18102 MM1, G015.05 MM1, and G015.31 MM3, which are clearly lower than our values (by factors  ~2−6). Sakai et al. (2008) assumed that T_ex(N₂H⁺) = T_rot(NH₃) even though the T_ex(N₂H⁺) values are considerably lower than T_rot(NH₃) (their Table 4). On the other hand, T_ex [N₂H⁺(1−0)]  from SSK08 are very similar to those from our RADEX simulations for the J = 3 − 2 transition. Also, the line optical thicknesses SSK08 list in their Table 4 indicate optically thick emission (τ ~ 2−11). This can explain the discrepancy in the derived column densities.

Marseille et al. (2008) modelled the N₂H⁺, HCO⁺, and H¹³CO⁺ emission of I18151 MM2, and obtained abundances that are comparable within the errors to those we have derived. This strengthens the reliability of our HCO⁺ abundances derived from H¹³CO⁺ by utilising the [¹²C]/[¹³C] abundance ratio.

Sakai et al. (2010, hereafter SSH10) derived H¹³CO⁺ column densities from the J = 1 − 0 transition for I18102 MM1, I18151 MM2, and I18223 MM3. Their results are otherwise similar to ours except that for I18102 MM1 the column densities differ by a factor of  ~4. Sakai et al. (2010) assumed optically thin emission and that T_ex(H¹³CO⁺) = T_rot(CH₃OH) ± 5 K = 16.7 ± 5 K ~ 3.6 × T_ex [H¹³CO⁺(3−2)] . These authors also derived SiO column densities from the J = 2 − 1 transition for the above three sources. We have obtained 36–86 times higher column densities for these sources from the J = 6 − 5 transition through RADEX calculations. The spatial resolution of SiO observations by SSH10 was 18″, and as T_ex(SiO) they used T_rot(CH₃OH) ± 5 K, or assumed that T_ex(SiO) = 20 K, which are over four times higher than T_ex [SiO(6 − 5)]  derived here. For comparison, we performed the LTE modelling of SiO(6−5) lines in Weeds, and found that column densities comparable to those obtained from RADEX are needed to explain the observed line intensities. Of course, similarly to SSH10, our RADEX and Weeds analyses of a single SiO transition include several (uncertain) assumptions. The large discrepancy in the derived N(SiO) values could also be due to the fact that the J = 2−1 and J = 6−5 lines, which have completely different upper-state energies (6.3 K and 43.8 K, respectively), originate in different parts within a clump.

In Fig. 8 we show stock charts of the derived molecular column densities and fractional abundances. For comparison, we plot column densities and abundances from several previous studies. The quoted reference studies deal with IRDCs (Ragan et al. 2006; SSK08; Gibson et al. 2009; Beuther & Henning 2009; Chen et al. 2010; SSH10; Vasyunina et al. 2011), high-mass young stellar objects (HMYSOs; Fontani et al. 2006; Thomas & Fuller 2008), and low-mass starless and protostellar cores (Butner et al. 1995). The latter study was included because we found no other reported DCO⁺ surveys towards IRDCs or HMYSOs. We also plot values from Fontani et al. (2011), who studied a sample of 27 high-mass cores (including starless cores, HMYSOs, and UC HII regions). The C¹⁸O abundances from Ragan et al. (2006) were converted into x(C¹⁷O) by using the abundance ratio [¹⁸O]/[¹⁷O]  = 3.52.

As can be seen in Fig. 8, the column densities and abundances we have derived are generally comparable to those observed in other sources. The most evident difference is the high SiO and ¹³CS column densities and abundances we have obtained compared to the IRDCs studied by SSH10 and Vasyunina et al. (2011).

Sakai et al. (2010) found that, in general, the N(SiO)/N(H¹³CO⁺) ratio is higher for the MSX dark sources than for the sources associated with MIR emission, and suggested that this ratio represents the fraction of the shocked gas. For the sources I18102 MM1, I18151 MM2, and I18223 MM3, SSH10 found the ratios $\mathit{N}\mathrm{\left(}\mathrm{SiO}\mathrm{\right)}\mathit{/}\mathit{N}\mathrm{\left(}{\mathrm{H}}^{\mathrm{13}}\mathrm{C}{\mathrm{O}}^{\mathrm{+}}\mathrm{\right)}\mathrm{=}\mathrm{2.}{\mathrm{9}}_{-1.2}^{\mathrm{+}\mathrm{1.9}}$, $\mathrm{2.}{\mathrm{4}}_{-0.9}^{\mathrm{+}\mathrm{1.5}}$, and $\mathrm{1.}{\mathrm{6}}_{-0.6}^{\mathrm{+}\mathrm{1.0}}$, respectively. Using our column densities, these ratios become very different: 57.1 ± 46.8, 220 ± 74.4, and 80 ± 20, respectively. The latter values agree better with the general trend found by SSH10 for a sample of 20 sources. These authors suggested that the SiO emission from the MIR dark objects originates in newly formed shocks, whereas the SiO emission from more evolved, MIR bright objects could originate in gas shocked earlier in time. This is supported by the fact that we see the narrowest SiO line towards I18102 MM1, presumably the most evolved source in our sample with detected SiO emission (see Sect. B.1 of Appendix B). Follow-up studies of extensive source samples are needed to test the statistical significance of the hypothesis by SSH10.

5.2. Depletion and deuteration

The CO depletion factors derived in this work are small, only in the range  ~0.6−2.7. Thus, CO molecules do not appear to be significantly depleted, if at all, in the studied clumps. This conforms to the fact that the clumps show signs of star-formation activity, such as outflows, which presumably release CO from the dust grains back into the gas phase (Appendix B). The presence of outflow-shocked gas is revealed by the SiO line wings (particularly towards I18151 MM2) and very high SiO abundances found in the present work. Also, the temperature in the envelopes of I18102 MM1 and I18151 MM2 is, likely owing to heating by embedded YSOs, slightly higher than the CO sublimation temperature of 20 K (e.g., Aikawa et al. 2008). Moreover, relatively high cosmic-ray ionisation rates towards the clumps were found. Thus, the cosmic-ray impulsive heating could also play a role in returning CO into the gas phase (Hasegawa & Herbst 1993).

However, the derived CO depletion factors represent only the average values along the line of sight within the 278 beam. Because the studied clumps lie at kiloparsec distances from us, the observations could be “contaminated” by non-depleted gas along the line of sight (Fontani et al. 2006). It is also uncertain whether the “canonical” CO abundance used to calculate f_D is exactly correct. Its value is known to vary by a factor of  ~2 between different star-forming regions (Lacy et al. 1994). Finally, the H₂ column densities derived from dust emission suffer from the uncertainty of dust opacity, which may be a factor of 2–3. These uncertainties are likely to explain the unsual values of f_D < 1.

For comparison, Fontani et al. (2006) determined the amount of CO depletion in HMYSOs. They found CO depletion factors in the range 0.4–35.8 with a median value of 3.2. Thomas & Fuller (2008) found depletion factors  ≲ 10 for their sample of HMYSOs and deduced that the sources are young, a few times 10⁵ yr.

At high angular resolution it is possible to identify more CO depleted regions within the clumps. Zhang et al. (2009) studied two massive clumps in the IRDC G28.34+0.06 at  ~1″ resolution with the SMA. They found CO depletion factors up to  ~10² − 10³ towards the two clumps, where the highest values were found in the cores within one of the clumps. Similarly, with the  ~3″ resolution PdBI observations, Beuther & Henning (2009) found a high CS depletion factor of  ~100 in IRDC 19175-4.

The degrees of deuteration in our clumps are  ~0.0002−0.014 in HCO⁺ and  ~0.002−0.028 in N₂H⁺. Compared to the average cosmic D/H ratio of  ~ 1.5 × 10^-5 (Sect. 4.4), the derived values are  ~10¹ − 10³ times higher. Relatively low deuteration degrees in some of our sources are consistent with the oberved low CO depletion factors. CO is the main destroyer of H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and H₂D⁺, and thus its presence can lower the deuterium fractionation (e.g., Caselli et al. 2008). For example, Swift (2009) estimated the fractional abundance of ortho-H₂D⁺ in two IRDCs of only  ~3−5 × 10^-13, which is up to three orders of magnitude lower than found in low-mass cores (e.g., Caselli et al. 2008; Friesen et al. 2010). From a theoretical point of view, if there are no depletion gradients in the source, one would expect to see the equality R_D(HCO⁺) = R_D(N₂H⁺) (Rodgers & Charnley 2001b). That we found somewhat lower values for R_D(HCO⁺) than R_D(N₂H⁺) could be related to differential depletion of molecules, and radial density gradients (Caselli et al. 2002). For high-mass star-forming clumps, the density profile is found to be of the form n(r) ∝ r ^{− (1.6 ± 0.5)} (e.g., Beuther et al. 2002a). CO is the parent species of HCO⁺, but it destroys the N₂H⁺ molecules. Deuteration proceeds most efficiently in regions where CO is mostly frozen onto dust grains. In the warmer envelope layers, where CO is not depleted, HCO⁺ can have a relatively high abundance, i.e., a lower value of R_D(HCO⁺) compared to that of N₂H⁺ (Emprechtinger et al. 2009). In general, the amount of molecular deuteration is expected to decrease during protostellar evolution because of internal heating of the surrounding envelope (e.g., Emprechtinger et al. 2009). Indeed, there is a decreasing trend in R_D(HCO⁺) with the gas kinetic temperature as shown in Fig. 9. When the temperature becomes  ≳ 20 K, reaction 2 in Table 8 can proceed in both directions, and thus the destruction rate of H₂D⁺ increases. Consequently, other deuterated molecules that form via reactions with H₂D⁺ start to decrease in abundance.

For comparison, Fontani et al. (2006) determined R_D(N₂H⁺) in HMYSOs and found the values in the range  ≲ 0.004 − 0.02, with an average value of  ~0.015. Roberts & Millar (2007) derived the upper limits to R_D(N₂H⁺) of  <0.002 and  <0.006 for the hot cores G34.26 and G75.78, respectively. Our source I18102 MM1 appears similar in this regard. Chen et al. (2010) found that R_D(N₂H⁺) = 0.017 − 0.052 in three cores within the IRDC G28.34+0.06. Moreover, they found that R_D(N₂H⁺) is lower in the more evolved stages of protostellar evolution. Recently, Fontani et al. (2011) found values of R_D(N₂H⁺) = 0.012 − 0.7, 0.017– ≤ 0.4, and 0.017–0.08 for their sample of high-mass starless cores, HMYSOs, and UC HII regions, respectively. The above results are quite similar to our findings. Another study of deuteration in IRDCs is that by Pillai et al. (2007), who derived NH₂D/NH₃ column density ratios in the range 0.005–0.386 for their sample of IRDC clumps. More recently, Pillai et al. (2011) found [NH₂D/NH₃] ratios in the range 0.06–0.37 towards IRDCs. The [NH₂D/NH₃] ratios obtained by Pillai et al. are even higher than those observed in low-mass dense cores (see Pillai et al. 2011, and references therein).

Again, high-resolution interferometric observations could reveal substructures within the clumps with a higher degree of deuteration. For example, Fontani et al. (2006) found an average value R_D(N₂H⁺) ~ 0.01 towards the high-mass star-forming region IRAS 05345+3157, but at high angular resolution (~3″) they resolved it into two N₂D⁺ condensations each with ten times higher deuteration degree [R_D(N₂H⁺) = 0.11; Fontani et al. 2008]. These high values are comparable to those found in low-mass starless cores [Crapsi et al. 2005; R_D(N₂H⁺) < 0.02 − 0.44] and Class 0 protostellar envelopes [Emprechtinger et al. 2009; R_D(N₂H⁺) < 0.029−0.271]. For their sample of low-mass dense cores, Caselli et al. (1998) found R_D(HCO⁺) to lie in the range 0.025–0.07. For a sample of more massive cores, Bergin et al. (1999) found that most sources have values of R_D(HCO⁺) < 0.035. This trend conforms to our results towards more massive clumps. It is possible that the cold phase during which the deuterium fractionation takes place is so short for HMYSOs, that very high deuteration degrees are not reached (Roberts & Millar 2007). On the other hand, the high [NH₂D/NH₃] ratios found by Pillai et al. (2007, 2011) could indicate a different production mechanism of deuterated ammonia compared to those of N₂H⁺ and HCO⁺, which are both purely produced in the gas phase.

One caveat to the deuteration analysis is that the FWHM linewidths of the detected N₂H⁺ transitions are significantly broader than those of N₂D⁺ [Table 4, Col. (4)]. Similarly, the linewidths of H¹³CO⁺(3−2) and DCO⁺(3−2) are quite different compared to each other towards some of the sources. The above transition pairs have similar critical densities, and the lines were observed at comparable spatial resolutions. Therefore, the difference between the linewidths suggests that the corresponding transitions may not be tracing the same gas component. In this picture, the calculated deuteration degrees should only be interpreted as an average value along the line of sight.

5.3. Ionisation

We have calculated the average ionisation degree along the line of sight towards the clumps by utilising three different methods. The summed abundances of the ionic species yield lower limits to x(e) in the range  ~0.7−12.3 × 10^-8. A very simple chemical scheme including HCO⁺ and N₂H⁺ yielded another estimates for the lower limits, about 0.3 − 5.6 × 10^-8. Thirdly, we used the R_D(HCO⁺) values to derive an upper limit to fractional ionisation, and found the range x(e) = 0.2−29.2 × 10^-5. To our knowledge, these are the first reported estimates of the ionisation degree in IRDCs made so far. As discussed below, the derived lower limits are likely to be more realistic estimates of x(e), whereas the upper limits appear to be very high. The x(e) values should be taken as rough estimates because of the simplicity of the analytical model we used. For example, our analysis is based on chemical equilibrium, which may not be valid (Caselli et al. 1998; Lintott & Rawlings 2006). Other factors that cause uncertainties are the possible variation of CO depletion and fractional molecular abundances along the line of sight. Electron abundance strongly depends on the depletion of CO because CO destroys ionic species such as H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and N₂H⁺. In this regard, one could think that the lower limit to x(e) derived by summing up the abundances of ionic species is the most reliable one. However, there is also a possibility that the line emission from different molecules originate in different layers of the clump. Indeed, the different linewidths between N₂H⁺(3−2) and N₂D⁺(3−2) point towards this possibility (Sect. 5.2). In this case, the physical meaning of the summed abundance of ionic species is unclear. Furthermore, that we have no information about the abundances of some important molecular ions, such as H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and H₃O⁺, nor the metal ions, hampers the determination of x(e). Metal ions play an important role in the ionisation balance in the outer layers of dense clouds (Caselli 2002). That the sources show low degrees of deuteration conforms to the fairly high ionisation levels of gas. However, detailed chemical modelling is needed to perform a more accurate analysis of the ionisation degree.

The cosmic-ray ionisation rates of H₂ we found are in the range ζ_H2 ~ 1 × 10^-17 − 1 × 10^-15 s^-1. Six of the seven sources have estimated values between  ~1 × 10^-17 and  ~5 × 10^-16 s^-1, but the highest rate of  ~10^-15 s^-1 is seen towards I18102 MM1. W emphasise that the ζ_H2 values were derived from the lower limits to x(e). However, as discussed above, we believe that these estimates are more accurate than those using the upper limits to x(e). Also, the applicability of the volume-averaged H₂ number densities used to compute ζ_H2 is uncertain because all detected molecular-ion transitions have a much higher critical density (Hezareh et al. 2008). The standard relation between the electron abundance and the H₂ number density is x(e) ~ 1.3 × 10^-5n(H₂) ^− 1/2 (McKee 1989). This is based on the pure cosmic-ray ionisation with the rate 1.3 × 10^-17 s^-1 and includes no depletion of heavy elements. By using the densities derived in the present work, 1.1 − 18.5 × 10⁴ cm^-3, the standard relation yields the values x(e) ~ 3−12 × 10^-8. These are comparable to the estimated lower limits to x(e).

How do our results compare to those found in other sources? Caselli et al. (1998) determined x(e) in a sample of 24 low-mass cores consisting of both starless and protostellar objects. Their analysis was based on observations of CO, HCO⁺, and DCO⁺, and the resulting values were in the range 10^-8 to 10^-6. They argued that the variation in x(e) among the sources is caused by variations in metal abundance and ζ_H2. Williams et al. (1998) used observations of C¹⁸O, H¹³CO⁺, and DCO⁺ to determine the values 10^-7.5 ≲ x(e) ≲ 10^-6.5 in a similar sample of low-mass cores as Caselli et al. (1998), but they used a slightly different analysis. Applying the same analysis as Williams et al. (1998), Bergin et al. (1999) found the ionisation levels of 10^-7.3 ≲ x(e) ≲ 10^-6.9 towards more massive cores (in Orion) than to those studied by Williams et al. The most massive sources were found to have the lowest electron abundances, x(e) < 10^-8. This conforms to the results by de Boisanger et al. (1996), who found that x(e) ~ 10^-8 in the massive star-forming regions NGC 2264 IRS1 and W3 IRS5. Hezareh et al. (2008) studied the high-mass star-forming region DR21(OH) in Cygnus X, and found that x(e) = 3.2 × 10^-8. These are consistent with our lower limits to x(e) in massive clumps within IRDCs. The very high upper limits to x(e) we obtained, particularly that in I18102 MM1 (~10^-4), resemble those found in PDRs, where carbon can provide most of the charge and x(e) ≃ x(C⁺) (Goicoechea et al. 2009).

Observational results presented in the literature suggest relatively large environmental variations in ζ_H2¹⁵. van der Tak & van Dishoeck (2000) used H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ and H¹³CO⁺ observations to constrain the ionisation rate towards HMYSO envelopes, and obtained the best-fit value ζ_H2 = (2.6 ± 1.8) × 10^-17 s^-1. We found comparable rates towards three of our sources. Caselli et al. (1998) inferred the values of ζ_H2 spanning a range of two orders of magnitude between 10^-18 − 10^-16 s^-1 in low-mass cores. Some of this variation could be caused by different cosmic-ray flux in the source regions. Williams et al. (1998) deduced a mean value of ζ_H2 = 5 × 10^-17 s^-1. On the other hand, H$\begin{array}{}\mathrm{+}\\ \mathrm{3}\end{array}$ observations towards the diffure ISM suggest higher values in the range ζ_H2 ~ 0.5 − 1.2 × 10^-15 s^-1 (McCall et al. 2003; Indriolo et al. 2007). The ionisation rate towards the Galactic Centre clouds is also found to be very high, ζ_H2 > 10^-15 s^-1 up to  ~ 10^-13 s^-1 (e.g., Oka et al. 2005; Yusef-Zadeh et al. 2007; Goto et al. 2008). However, these high rates can explain the high temperatures of  ~ 70 − 100 K in the Galactic Centre clouds, whereas the dense YSO envelopes are typically much colder. van der Tak et al. (2006) estimated that ζ_H2 ~ 4 × 10^-16 s^-1 in the Sgr B2 envelope, whereas Hezareh et al. (2008) found a lower value of ζ_H2 = 3.1 × 10^-18 s^-1 towards DR21(OH). Observations of the Horsehead Nebula by Goicoechea et al. (2009) could only be reproduced with ζ_H2 = 7.7 ± 4.6 × 10^-17 s^-1. Despite all the uncertainty factors in ζ_H2, it seems possible that some of our clumps could be exposed to somewhat higher cosmic-ray flux, but the observed clump temperatures,  ≃ 20 K, suggest that ζ_H2 is  ~ 10^-17 − 10^-16 s^-1 (see Bergin et al. 1999).

Because (most of) our sources are associated with high-mass star formation, which could increase the local UV radiation field, the role played by photoionisation should also be considered. Photoionisation is expected to become negligible when the extinction is A_v ≳ 4 mag (McKee 1989). In units of N(H₂) this corresponds to about  ≳ 4.9 × 10²¹ cm^-2 (Vuong et al. 2003). Therefore, the ionisation in the studied clumps is very likely dominated by cosmic ray particles and the ambient UV radiation field is not important in this regard. If the clumps are inhomogeneous, however, external UV radiation could penetrate more efficiently into them and enhance the ionisation level (Boissé 1990). Our sources also show high HCO⁺ abundances (~5 × 10^-9 − 9 × 10^-8), whereas in the presence of significant UV radiation one would expect to observe lower values of x(HCO⁺) because the photoproduced electrons raise the HCO⁺ recombination rate (Jansen et al. 1995). On the other hand, the HCO⁺ abundance could be enhanced by the presence of protostellar outflows where shock-generated UV radiation is present (Rawlings et al. 2000). This could explain the high HCO⁺ abundance in I18102 MM1 (Sect. B.1).

5.4. Coupling between the ions and neutrals

The knowledge of the gas ionisation degree enables one to estimate the level of coupling between the neutral gas component and the pervasive magnetic field. The coupling strength can be quantified through the wave-coupling number, which is defined as the ratio of the clump radius to the MHD cutoff wavelength, W ≡ R/λ_min (Myers & Khersonsky 1995). The wavelength λ_min represents the minimum wavelength for propagation of MHD waves, and is given by λ_min ≡ πv_Aτ_in, where v_A is the Alfvén speed, and τ_in is the ion-neutral collision time. When λ < λ_min, the frequency of the wave is higher than the collision frequency between the ions and neutrals (ν ≫ ν_in). In this case, the neutral fluid component is decoupled from the magnetic-field dynamics. If the clump is much larger than λ_min, i.e., W ≫ 1, the magnetic field is strongly coupled to the neutral component, and MHD waves can stabilise the clump.

We calculated the W number by utilising Eq. (6) in Williams et al. (1998):

$\mathit{W}\mathrm{=}\frac{\mathit{x}\mathrm{\left(}\mathrm{e}\mathrm{\right)}\mathit{n}\mathrm{(}{\mathrm{H}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{⟨}{\mathit{\sigma }}_{\mathrm{in}}\mathit{v}\mathrm{⟩}}{\mathit{\pi }}{\begin{array}{}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{5}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{p}\mathrm{)}}{\mathrm{4}\mathit{\pi \mu }{\mathit{m}}_{\mathrm{H}}\mathit{G}\mathrm{-}\mathrm{15}{\mathit{\sigma }}_{\mathrm{T}}^{\mathrm{2}}\mathit{/}\left[\mathit{n}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}{\mathit{R}}^{\mathrm{2}}\right]}\begin{array}{}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(14)where  ⟨ σ_inv ⟩  is the ion-neutral collision rate, p ∈  [0,   1]  is a parameter describing the level of turbulence (0 for a minimum amount of turbulent motions, 1 for a maximum level of turbulence), μ is the mean molecular weight per particle (taken to be 2.33), G is the gravitational constant, and ${\mathit{\sigma }}_{\mathrm{T}}\mathrm{=}\sqrt{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{kin}}\mathit{/}\mathit{\mu }{\mathit{m}}_{\mathrm{H}}}$ is the one-dimensional thermal velocity dispersion. Following Williams et al. (1998), we take  ⟨ σ_inv ⟩  = 1.5 × 10^-9 cm³ s^-1. We used as x(e) the summed abundance of ionic species, and consequently the values of W should be taken as lower limits. We assumed that the clumps are turbulent (p = 1), and clump radii were taken to be the effective radii listed in Col. (7) of Table 5. Note that the latter were also used to derive the H₂ number densities. The computed values of W and λ_min are given in Table 10. The derived W values are in the range  ~ 2 − 19, whereas the minimum wavelengths lie between 0.018 and 0.174 pc. Note that the clump radii range from 0.15 to 0.34 pc, and are mostly comparable to each other. Thus, the large variation of the W values cannot (solely) be explained by different source sizes. Because W represents the ratio of maximum to minimum wavelength, the maximum wavelength corresponds to the clump radius. For comparison, Bergin et al. (1999) deduced that massive cores in Orion (comparable in size to our clumps) are characterised by the value W = 20. The MHD wave power transmission at the lowest frequency of wave propagation (i.e., the longest wavelength) rapidly drops below unity when W ≲ 100 (Myers & Lazarian 1998; their Fig. 1). Therefore, the values W < 100 imply that the coupling is “marginal”, and that MHD waves have a quite narrow band of wavelengths to propagate above λ_min. However, for I18102 MM1 and J18364 SMM1 the allowed wavelengths of MHD waves could extend from the effective radii  ~ 0.2 − 0.3 pc down to λ_min ~ 0.02 − 0.04 pc. This is expected to lead to a clumpy structure of a medium (Myers & Khersonsky 1995), just like is found to be the case in J18364 SMM1, which is fragmented into two subcores (see Sect. B.6).

Another approach to investigate the coupling of the neutral gas and the magnetic field is to determine the ratio between the ambipolar diffusion (AD) timescale and free-fall timescale, τ_AD/τ_ff. In the process of AD, neutrals drift quasi-statically (under their own self-gravity) relative to ions and magnetic field (Mestel & Spitzer 1956). It can be shown that the ratio τ_AD/τ_ff is approximately given by W ~ τ_ff/τ_in (see Mouschovias 1991; Williams et al. 1998). The free-fall timescales of our clumps, ${\mathit{\tau }}_{\mathrm{ff}}\mathrm{=}\sqrt{\mathrm{3}\mathit{\pi }\mathit{/}\mathrm{32}\mathit{G}\mathrm{⟨}\mathit{\rho }\mathrm{⟩}}\mathrm{~}\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{-}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}$ yr, imply AD timescales in the range τ_AD ~ W × τ_ff ~ 2 × 10⁵ − 4 × 10⁶ yr [see the first values in Col. (4) of Table 10]. We note that τ_AD strongly depends on the ionic composition and can be even much longer than the above rough estimates. Assuming that HCO⁺ is the dominant ionic species, the AD timescale is given by τ_AD ~ 1.3 × 10¹⁴x(e) yr (see Eq. (5) in Walmsley et al. 2004). Using again the lower limits to x(e), we obtain timescales in the range τ_AD ~ 1 − 16 × 10⁶ yr, which are about four times longer than the above estimates [see the second values in Col. (4) of Table 10]. Therefore, on the scale probed by our single-dish observations, AD is not expected to be a main driver of clump evolution unless it occurs on timescales much longer than 10⁶ yr. However, the situation may be different at smaller scales where the density is higher, and therefore x(e) is lower and τ_AD is shorter.

6. Summary and conclusions

We have carried out a molecular-line study of seven massive clumps associated with IRDCs. We used APEX to observe transitions of C¹⁷O, H¹³CO⁺, DCO⁺, N₂H⁺, and N₂D⁺. The principal aim of this study was to investigate the depletion, deuterium fractionation, and the degree of ionisation in the sources. Our main results are summarised as follows:

• 1.

  The CO molecules do not appear to be significantly depleted, ifat all, in the observed sources. The largest CO depletion factor, ~2.7, is found towards G015.31 MM3, which is dark at MIR. Inmany sources, CO is likely to be released from dust grains into thegas phase owing to shock evaporation.

• 2.

  The deuteration degree in HCO⁺ was found to range from 0.0002 to 0.014, whereas that in N₂H⁺ lies in the range 0.002–0.028. These are lower values than found in low-mass starless cores and protostellar envelopes, but still significantly higher than the cosmic D/H-ratio  ~10^-5. The degree of deuteration in HCO⁺ appears to decrease with increasing gas temperature, as expected theoretically. This likely reflects the evolutionary stage of the clump: at early stages the source is cold and as the evolutionary stage progresses the temperature rises. In the course of evolution the amount of deuterium fractionation decreases.

• 3.

  For the first time, we have estimated the level of ionisation in IRDCs. The ionisation degree is difficult to determine accurately, especially from an analytical basis like the one in the present work. However, a lower limit to electron abundance appears to be x(e) ≳ 10^-8−10^-7, which is comparable to those seen in other star-forming regions. The cosmic-ray ionisation rate of H₂ were estimated to vary from  ~10^-17 s^-1 up to  ~10^-15 s^-1 in one target position. Similar estimates have been found in a variety of Galactic environments.

• 4.

  Additional molecular species were detected towards some of the sources. In particular, ortho-c-C₃H₂ was detected in four sources, and SiO towards three sources. This is the first reported observation of c-C₃H₂ towards IRDCs. The SiO lines showed broad linewidths indicative of outflows. Also for the first time towards IRDCs, the deuterated ethynyl, C₂D, was detected in I18151-1208 MM2. We estimate the C₂D/C₂H column density ratio to be about 0.07, which agrees well with the results from other sources (Appendix A).

• 5.

  The ambipolar-diffusion timescale implied by the lower limits to x(e) is typically  ~several Myr. Therefore, on the scale probed by our observations, clump evolution cannot be mainly driven by ambipolar diffusion unless it occurs on timescales much longer than 10⁶ yr.

The present single-dish study cannot follow the chemical structure of the sources down to small scales where the evolution of the massive young stellar object(s) is taking place (e.g., the case of J18364 SMM1). Follow-up high-resolution studies are needed i) to resolve whether the clumps are fragmented into smaller subunits (as is known to be the case in J18364 SMM1); ii) to examine the kinematics (infall, outflow, relative motions) of the subfragments, and, importantly, their chemical partition (such as depletion and deuteration); and iii) to locate the detected complex molecules (envelope or hot core?). In this regard, follow-up observations with ALMA seem particularly worthwhile. Finally, detailed chemical models are needed to better understand the observed molecular abundances, and how these are related to the timescale of the high-mass star-formation process.

Online material

Appendix A: Other detected transitions

Besides the molecular species discussed in Sect. 5.1, the sources I18102 MM1, I18151 MM2, I18182 MM2, and I18223 MM3 show emission lines also from other species, some of which are complex organics. These include o-c-C₃H₂, C₂D, C¹⁸O (seen in the image band), and possibly CNCHO (in the image band), CH₃NH₂ (blended with OCS), CH₂CHC¹⁵N, and CH₃COCH₃. In this appendix we discuss each of these species (except C¹⁸O) and their derived column densities and abundances.

Cyclopropenylidene (c-C₃H₂).    We have detected the ortho form of cyclic-C₃H₂ in four sources, and found the column densities and abundances in the range  ~7.0 × 10¹² − 7.7 × 10¹³ cm^-2 and  ~1.9 × 10^-10 − 1.2 × 10^-8, respectively. The c-C₃H₂ molecule has been detected in several different Galactic sources, but to our knowledge, this is the first reported detection of this molecule towards massive clumps associated with IRDCs. For comparison, the beam-averaged c-C₃H₂ column densities and abundances of  ~4−9 × 10¹³ cm^-2 and 3−7 × 10^-11 have been found in Sgr B2 (Turner 1991; Nummelin et al. 2000). These column densities are comparable to those in our clumps, but the abundances appear lower in Sgr B2.

The major formation channel of c-C₃H₂ is ${\mathrm{C}}_{\mathrm{3}}{\mathrm{H}}^{\mathrm{+}}\mathrm{+}{\mathrm{H}}_{\mathrm{2}}\mathrm{\to }\mathrm{c}\mathrm{-}{\mathrm{C}}_{\mathrm{3}}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}$; $\mathrm{c}\mathrm{-}{\mathrm{C}}_{\mathrm{3}}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{c}\mathrm{-}{\mathrm{C}}_{\mathrm{3}}{\mathrm{H}}_{\mathrm{2}}$ (e.g., Park et al. 2006). The very low rotational excitation temperatures we derived, 3.1–4.8 K, suggest that c-C₃H₂ emission mainly comes from the cool and relatively low-density envelope (Turner 1991). Also, higher electron abundance in the outer layers can promote the formation of c-C₃H₂. This is consistent with the derived high fractional abundances.

Cyanoformaldehyde or formyl cyanide (CNCHO).    The CNCHO transition is possibly seen in the image sideband towards I18102 MM1. Because of attenuated intensity, we cannot determine the column density of this molecule. Remijan et al. (2008) determined the values N(CNCHO) = 1−17 × 10¹⁴ cm^-2 and x(CNCHO) = 0.7 − 11 × 10^-9 towards Sgr B2(N). They suggested that CNCHO is likely formed in a neutral-radical reaction of formaldehyde (H₂CO) and the cyanide (CN) radical.

Methylamine (CH₃NH₂).    Rotational levels of CH₃NH₂ undergo the so-called A- and E-type transitions due to the internal rotation of the CH₃-group. We have a possible detection of the E-type CH₃NH₂ transition, and we derived the values N(CH₃NH₂) = 5 × 10¹⁴ cm^-2 and x(CH₃NH₂) = 5.1 ± 0.6 × 10^-8. However, our CH₃NH₂ line is blended with OCS line, and therefore these values should be taken with caution. Beam-averaged column densities of N(CH₃NH₂) ~ 1 × 10¹⁴ − 8 × 10¹⁵ cm^-2 have been found in Sgr B2 (Turner 1991; Nummelin et al. 2000). The origin of CH₃NH₂, and nitrogen-bearing organics in general, could be in the evaporation of ice mantles, indicative of a hot-core chemistry (e.g. Rodgers & Charnley 2001a).

Carbonyl sulfide (OCS).    The high-J transition of OCS we have possibly detected (blended with CH₃NH₂) implies a column density and fractional abundance of 2.3 ± 0.8 × 10¹⁵ cm^-2 and 2.4 ± 0.9 × 10^-7, respectively. For example, Qin et al. (2010) detected OCS in the high-mass star-forming region G19.61-0.23. Assuming T_rot = E_u/k_B, they derived the beam-averaged column density and fractional abundance of 2.2 ± 0.1 × 10¹⁶ cm^-2 and 2.7 ± 0.1 × 10^-8, respectively. Note that Qin et al. achieved a much higher angular resolution with their SMA observations, and thus their high column density value could partly be caused by filtering out the extended envelope. Our coarser angular resolution probably causes a significant beam dilution. Sakai et al. (2010) derived OCS column densities of 1.2 − 5.5 × 10¹⁴ cm^-2 for their sample of 20 massive clumps associated with IRDCs, i.e., lower than found here towards I18102 MM1.

The OCS molecules form on grain surfaces through the addition of S atom to CO, or via the O atom addition to CS (e.g., Charnley et al. 2004). Solid-state OCS has been detected towards high-mass star-forming regions by, e.g, Gibb et al. (2004). Sakai et al. (2010) suggested that OCS is released from the grains into the gas phase through protostellar shocks. This could well be the case in I18102 MM1 (Sect. B.1).

Deuterated ethynyl (C₂D).    We have made the first detection of C₂D towards IRDCs. The column density and abundance estimated from the N = 3 − 2 transition are  ~1.5 × 10¹³ cm^-2 and 4.2 ± 0.5 × 10^-10, respectively. Sakai et al. (2010) detected the normal isotopologue C₂H(N = 1 − 0) towards I18151 MM2. The C₂H column density they derived,  ~2.2 × 10¹⁴ cm^-2, together with N(C₂D) derived by us, suggest a deuteration degree of  ~0.07 in C₂H. Vrtilek et al. (1985) detected the C₂D(N = 2 − 1) transition near the Orion-KL position, and obtained the column density  ~1.8 × 10¹³ cm^-2. Moreover, they derived the N(C₂D)/N(C₂H) ratio of 0.05. These are very similar to what we have found. Quite similarly, Parise et al. (2009) found upper limits of N(C₂D) < 2.5    ×    10¹³ cm^-2 and x(C₂D) < 2 × 10^-10 in a clump associated with the Orion Bar.

The formation of C₂D is believed to take place in the gas phase through the route $\mathrm{C}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{\to }\mathrm{C}{\mathrm{H}}_{\mathrm{2}}{\mathrm{D}}^{\mathrm{+}}\mathrm{\to }{\mathrm{C}}_{\mathrm{2}}\mathrm{D}$ (see, e.g., Parise et al. 2009). The N(C₂D)/N(C₂H) ratio we have obtained is comparable with those predicted by the low-metal abundance model by Roueff et al. (2007) at temperature around 30–40 K. We note that the C₂H abundance,  ~6 × 10^-9, calculated from the observed deuteration degree (0.07) is comparable to the values x(C₂H) = 2.5 × 10^-9 − 5.3 × 10^-8 recently found by Vasyunina et al. (2011) towards IRDCs.

Acrylonitrile or vinyl cyanide, ¹⁵N isotopologue (CH₂CHC¹⁵N).    The Weeds modelling suggests that there is a CH₂CHC¹⁵N transition blended with the hf group of the detected C₂D line. Thus, reliable column density estimate cannot be performed. Assuming that the detected line is completely due to CH₂CHC¹⁵N emission, we derive very high values of 2.5 × 10¹⁶ cm^-2 and  ~7 × 10^-7 for the column density and fractional abundance. For comparison, a column-density upper limit of 3 × 10¹⁵ cm^-2 towards Sgr B2(N) was derived by Müller et al. (2008).

If present, this molecule would indicate a hot-core chemistry. The main isotopologue C₂H₃CN is expected to form through gas-phase reactions after the ethyl cyanide (C₂H₅CN), forming on dust grains, evaporates into the gas phase (Caselli et al. 1993).

Acetone (CH₃COCH₃).    The first detection of CH₃COCH₃ in the interstellar medium was made by Combes et al. (1987) towards Sgr B2. They found the column density and fractional abundance of this molecule to be 5 × 10¹³ cm^-2 and 5 × 10^-11. These are significantly lower than we have estimated towards I18182 MM2 (2 × 10¹⁵ cm^-2 and  ~4 × 10^-7). Combes et al. (1987) found that the CH₃COCH₃ abundance is about 1/15 of its precursor molecule CH₃CHO (acetaldehyde), and suggested the formation route of acetone to be $\mathrm{C}{\mathrm{H}}_{\mathrm{3}}^{\mathrm{+}}\mathrm{+}\mathrm{C}{\mathrm{H}}_{\mathrm{3}}\mathrm{CHO}\mathrm{\to }\mathrm{(}\mathrm{C}{\mathrm{H}}_{\mathrm{3}}{\mathrm{)}}_{\mathrm{2}}\mathrm{CH}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}\mathit{h\nu }\mathrm{;}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathrm{C}{\mathrm{H}}_{\mathrm{3}}{\mathrm{)}}_{\mathrm{2}}\mathrm{CH}{\mathrm{O}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}}\mathrm{\to }\mathrm{C}{\mathrm{H}}_{\mathrm{3}}\mathrm{COC}{\mathrm{H}}_{\mathrm{3}}\mathrm{+}\mathrm{H}$. However, Herbst et al. (1990) showed that this radiative association reaction is likely too slow to be consistent with the observed abundance in Sgr B2. The chemistry behind the formation of acetone is not clear. It could be caused by some other gas-phase ion-molecule reactions, or caused by grain chemistry (Herbst et al. 1990).

Friedel et al. (2005) found that acetone in Orion BN/KL is concentrated towards the hot core. They derived the beam-averaged column densities of  ~2−8 × 10¹⁶ cm^-2. More recently, Goddi et al. (2009) found the column density N(CH₃COCH₃) = 5.5 × 10¹⁶ cm^-2 in Orion BN/KL, in agreement with the Friedel et al. results. The very high column densities and abundances of CH₃COCH₃ found in Orion BN/KL and in I18182 MM2 in the present work require some other reaction pathway(s) than the above radiative association reaction to be efficient. Grain surface and hot-core gas-phase chemistry may both play critical roles (Garrod et al. 2008).

The detection of complex molecules in clumps associated with IRDCs, indicating the presence of hot cores, supports the scenario that high-mass star formation can take place in these objects. For example, Rathborne et al. (2007, 2008) found that the clumps G024.33+00.11 MM1 and G034.43+00.24 MM1 are both likely to contain a hot molecular core. On the other hand, complex organics could also be ejected from grain mantles through shocks, and cosmic rays heating up the dust can also have some effect (e.g., Requena-Torres et al. 2008, and references therein).

Appendix B: Discussion on individual sources

B.1. IRDC 18102-1800 MM1

The clump I18102 MM1 is the warmest (21.3 K) and most massive (~36 M_⊙) source of our sample. It is associated with Spitzer point sources at 8 and 24 μm, and high-mass star formation is taking place within it as indicated by the presence of the 6.7 GHz Class II CH₃OH maser (Beuther et al. 2002b). Among our sample, the lowest degree of deuteration in both HCO⁺ and N₂H⁺ is found for I18102 MM1. On the other hand, the source shows the highest degree of ionisation.

Fuller et al. (2005) detected central dips and red asymmetries in the spectral lines HCO⁺(1 − 0), HCO⁺(4 − 3), N₂H⁺(1 − 0), and H₂CO(2_1,2 − 1_1,1). A similar line profile is seen in the J = 3 − 2 transition of N₂H⁺ in the present study. These profiles indicate the presence of expanding or outflowing gas (see, e.g., Park et al. 2000 and references therein). Beuther & Sridharan (2007) detected very broad (42.2 km s^-1 down to zero intensity) SiO(2−1) wings towards I18102 MM1, indicative of bipolar outflows. This source showed the second-broadest SiO line in the sample of Beuther & Sridharan (2007). Also, the SiO(2−1) line detected by SSH10 towards I18102 MM1 was very wide, Δv = 13 ± 1 km s^-1. The SiO(6−5) line detected in the present study has Δv ≃ 4.7 km s^-1, and the width down to zero intensity is  ~10 km s^-1. We have also detected the ¹³C isotopologue of CS in this source. Sakai et al. (2010) suggested that CS could originate in the shock evaporation of grain mantles or radiative heating.

Beuther & Sridharan (2007) detected CH₃CN(6_K − 5_K) and CH₃OH(5_K − 4_K) lines towards I18102 MM1, and derived the abundances of 3 × 10^-11 and 4 × 10^-11 for CH₃CN and CH₃OH, respectively. Sakai et al. (2008) detected CH₃OH(7_K − 6_K) and HC₃N(5−4) lines in this source, and SSH10 detected C₂H(N = 1 − 0) and a hint of CH₃OH(2₁ − 1₁) A⁻. We note that SSK08 and SSH10 did not detected the lines of CCS(4₃ − 3₂), SO(2₂ − 1₁), or OCS(8 − 7) in their surveys. The upper limit they derived for the OCS column density,  ≲ 2.7 × 10¹⁴ cm^-2, is about 8.5 times less than the N(OCS) value we obtained. The observational results cumulated so far indicate hot-core chemistry in I18102 MM1. This conforms to the fact that this clump is giving birth to a high-mass star(s), and possibly through disk accretion as indirectly suggested by the outflow signatures. The associated hot core is likely to be in its later stages of evolution because it is associated with a methanol maser (cf. Rathborne et al. 2008).

B.2. G015.05+00.07 MM1 and G015.31-00.16 MM3

G015.05 MM1 is the lowest mass (63 M_⊙) and G015.31 MM3 is the second-coldest (13.7 K) clump of our sample. Both clumps are dark in the Spitzer 8 and 24 μm images. The highest deuteration degree in N₂H⁺ (0.028) was found towards G015.05 MM1, whereas G015.31 MM3 shows the largest CO depletion factor (2.7) in our sample. G015.05 MM1 is associated with H₂O maser (Wang et al. 2006), indicative of star-formation activity. Rathborne et al. (2010) derived the following properties for G015.05 MM1 from the broadband SED: T_dust = 11 − 36 K, L = 15.5 − 362   L_⊙, and M = 35 − 158   M_⊙ (scaled to the revised distance 2.6 kpc). The values T_kin = 17.2 K and M = 63   M_⊙ derived in the present study lie at the low end of the above values of temperature and mass (at the derived clump densities, it is expected that T_kin ≃ T_dust, Goldsmith & Langer 1978).

Sakai et al. (2008) barely detected the HC₃N(5−4), CH₃OH(7_K − 6_K), and CCS(4₃ − 3₂) lines towards G015.05 MM1 and G015.31 MM3. Indeed, all the CH₃OH detected objects in the survey by SSK08 are associated with the Spitzer 24-μm sources.

Both G015.05 MM1 and G015.31 MM3 are likely to be in a very early stage of evolution. Moreover, both sources are massive enough to allow high-mass star formation. G015.31 MM3 could represent or host the so-called high-mass prestellar core, whereas G015.05 MM1 could be slightly more evolved with H₂O maser emission but still lacking IR emission at 8 and 24 μm. The small wave-coupling number of W ~ 2 for G015.31 MM3 suggests that “magnetic turbulence” is not able to fragment the clump into smaller pieces, strengthening the possibility that it is a massive prestellar “core”.

B.3. IRDC 18151-1208 MM2

The I18151 MM2 clump is dark in the Spitzer 8-μm image (24 μm not available). It is associated with H₂O (Beuther et al. 2002b) and Class I CH₃OH masers (Marseille et al. 2010b). This clump has the highest volume-averaged H₂ number density (~1.9 × 10⁵ cm^-3) among our sources.

Beuther & Sridharan (2007) detected the broadest SiO(2−1) wing emission (65 km s^-1 down to zero intensity) towards I18151 MM2 in their sample. Also, the J = 2 − 1 and 3 − 2 SiO lines detected recently by López-Sepulcre et al. (2011) are very broad, i.e., FWZP = 84.3 and 103.1 km s^-1, respectively. The SiO(6 − 5) line we detected is also very broad with the FWHM 46.8 km s^-1. The outflow activity within the clump was confirmed by Marseille et al. (2008) who, for the first time, found that I18151 MM2 is driving a CO outflow and hosts a mid-IR-quiet, possibly a Class 0-like HMYSO (cf. Motte et al. 2007). Marseille et al. (2008) modelled the dust continuum emission (SED) of I18151 MM2 and found that the bolometric luminosity, mass, and the mean temperature of the source are L_bol = 2190   L_⊙, M_gas = 373 ± 81   M_⊙, and  ⟨ T ⟩  = 19.4 ± 0.2 K (scaled to the revised distance 2.7 kpc). Within the errors, these mass and temperature values are comparable to the values derived in the present paper. Besides the Class I methanol maser tracing the molecular outflow, Marseille et al. (2010b) detected blue asymmetry in CH₃OH(5 _− 1,5 − 4_0,4) E, indicating infall motions. The 22 GHz H₂O maser is probably excited in the outflow shocked gas (cf. Furuya et al. 2011).

Beuther & Sridharan (2007) also detected CH₃CN(6_K − 5_K) and CH₃OH(5_K − 4_K) lines towards I18151 MM2, which is a sign of hot-core chemistry, and they derived the abundances of 8 × 10^-11 and 6 × 10^-10 for CH₃CN and CH₃OH, respectively. SSK08 detected CH₃OH(7_K − 6_K) lines towards I18151 MM2, but not CCS(4₃ − 3₂) or HC₃N(5−4). In their line survey, SSH10 detected C₂H(N = 1 − 0), but not SO(2₂ − 1₁), OCS(8 − 7), or CH₃OH(2₁ − 1₁) A⁻ lines. The C₂H column density they derived,  ~ 2.2 × 10¹⁴ cm^-2, together with the value N(C₂D) = 1.5 × 10¹³ cm^-2 derived by us, suggest a deuteration degree of  ~ 0.07 in C₂H. Marseille et al. (2008) concluded from their modelling of CS transitions that CS is depleted in I18151 MM2. We derived only a small CO depletion factor of  ~ 1.6 for this clump, and low degrees of deuteration in HCO⁺ and N₂H⁺, namely 0.3% and 1%, respectively. Also, this source shows the lowest lower limit to ionisation degree, only 7 × 10^-9.

B.4. IRDC 18182-1433 MM2

The filamentary clump I18182 MM2 is associated with Spitzer 8 and 24-μm sources. This source shows the highest degree of deuteration in HCO⁺ (0.014).

Beuther & Sridharan (2007) detected CH₃OH(5_K − 4_K) lines towards I18182 MM2, and derived the CH₃OH abundance of 2.1 × 10^-10. Sakai et al. (2008) did not detect CCS(4₃ − 3₂), CH₃OH(7_K − 6_K), or HC₃N(5−4) lines towards this source in their survey. Our tentative detection of the O-bearing species CH₃COCH₃ indicates hot-core chemistry within the clump. Moreover, O-bearing species are sign of the early stage of chemical evolution (e.g., Shiao et al. 2010). This conforms to the presumable young age of the clump as it is associated with IRDC. Interestingly, the nearby clump I18182 MM1, which is associated with the HMYSO IRAS 18182-1433, also appears to contain a hot core (Beuther et al. 2006).

B.5. IRDC 18223-1243 MM3

I18223 MM3, part of a long filamentary IRDC, is a high-mass clump harbouring an embedded accreting low- to intermediate-mass protostar that could evolve to a high-mass star at some point in the future (see Beuther & Steinacker 2007; Beuther et al. 2010). Fallscheer et al. (2009) detected a molecular outflow in this source and found evidence for a large rotating structure, or toroid, perpendicular to the outflow. Evidence for outflow activity in this source was already found by Beuther et al. (2005c) and Beuther & Sridharan (2007), who found that there are 4.5 μm emission features at the clump edge and that the spectral lines of CO, CS, and SiO show broad wing emission. The line profile of SiO(6 − 5) detected in the present study also indicates outflowing gas.

The clump shows 24 μm emission but is dark at the Spitzer IRAC wavelengths (3.6, 4.5, 5.8, and 8.0 μm). Beuther et al. (2010) derived an SED for this source between 24 μm and 1.2 mm, including the recent Herschel PACS and SPIRE data, and obtained the total luminosity of 539 L_⊙ (scaled to d = 3.5 kpc). They also found that the ratio between the total and submm luminosity (integrated longward of 400 μm) is only 11, suggesting that the source is very young.

Beuther & Sridharan (2007) detected CH₃CN(6_K − 5_K) and CH₃OH(5_K − 4_K) lines towards I18223 MM3, indicative of

hot-core chemistry, and derived the abundances of 9 × 10^-11 and 6.1 × 10^-10 for CH₃CN and CH₃OH, respectively. Sakai et al. (2008) detected HC₃N(5−4), only weak CH₃OH(7_K − 6_K) lines, and no CCS(4₃ − 3₂) lines towards this source. Sakai et al. (2010) did not detect the lines SO(2₂ − 1₁), OCS(8 − 7), or CH₃OH(2₁ − 1₁) A⁻ in their survey; however, they detected the C₂H(N = 1−0) line. Based on the column densities of different species (e.g., SiO and H¹³CO⁺), SSH10 suggested that I18223 MM3 is in early stage of evolution. This conforms to the fact that the second-highest value of R_D(N₂H⁺) in our sample (0.013) is found towards I18223 MM3.

B.6. ISOSS J18364-0221 SMM1

At the distance of  ~2.5 kpc, the clump J18364 SMM1 is the nearest source in our sample. It is also the coldest (11.4 K) clump of our sample. The clump is dark at 8 μm but it is associated with the 24-μm point source. The second-highest value of R_D(HCO⁺) in our sample, 0.012, is found towards this source. It also shows high cosmic-ray ionisation rate of H₂ (~5 × 10^-16 s^-1). Birkmann et al. (2006) studied this clump through J = 3 − 2 transitions of HCO⁺ and H¹³CO⁺ (beam size  ~ ). The former line showed blue asymmetric profile, indicating infall. They also found significant CO(2−1) line wings, indicating the presence of outflows. The H¹³CO⁺(3−2) line we observed shows an asymmetric profile with a central dip and slightly stronger red peak, contrary to that observed by Birkmann et al. (2006) in HCO⁺(3−2). This difference is probably caused by the larger beam size of our observations (24″), yielding a signature of outflowing gas motions.

More recently, J18364 SMM1 was studied in detail by Hennemann et al. (2009). They resolved this contracting clump in the interferometric mm continuum into two compact cores, named SMM1 North and South separated by (0.12 pc). Their positions are indicated in Fig. 1. The peak H₂ column densities and dust temperatures were found to be 2.7 × 10²³ cm^-2 and 15 K for SMM1 North and, 2.4 × 10²³ cm^-2 and 22 K for SMM1 South. Using the revised distance 2.5 kpc, the radius, mass, and luminosity of the northern core are 0.06 pc, 19 M_⊙, and 26 L_⊙, whereas for the southern core these values are 0.05 pc, 13 M_⊙ and 230 L_⊙. Thus, the cores within the clump have comparable sizes and masses but the southern one is associated with the Spitzer 24 and 70 μm sources and is more luminous. Hennemann et al. (2009) found that SMM1 South drives an energetic molecular outflow and that the core centre is supersonically turbulent. On the other hand, the IR-dark core SMM1 North shows lower levels of turbulence, but it also drives an outflow. Both the outflows from SMM1 North and South are quite collimated and their estimated ages are  < 10⁴ yr. The HCN(1 − 0) modelling results by Hennemann et al. (2009) showed that the spectrum of SMM1 South can be explained with a collapse of the core. They obtained an infall velocity of 0.14 km s^-1 and an estimated mass infall rate of  ~ 3.4 × 10^-5   M_⊙ yr^-1 (scaled to the revised distance).

In summary, the clump J18364 SMM1 is fragmented into two cores that both harbour protostellar seeds, possibly evolving into intermediate- to high-mass stars. As discussed in Sect. 5.4, MHD wave propagation could have played a role in fragmenting the parent clump. As the southern core harbours a 24 μm source, is highly turbulent in the central region, and shows jet features at large distance from the driving source, it appears to be more evolved than the northern core.

---

Acknowledgments

We are indebted to the referee for his/her valuable comments, which helped significantly in making the paper more readable. We thank the staff at the APEX telescope for performing the service mode observations presented in this paper, and the ESO archive team for their help with the data products. O.M. acknowledges support from the Academy of Finland grant 132291. Many thanks to Henrik Beuther for providing us the MAMBO 1.2-mm dust continuum map of IRDC 18151-1208. We would also like to thank the people who maintain the CDMS and JPL molecular spectroscopy databases, and the Splatalogue Database for Astronomical Spectroscopy. The present study makes use of the JCMT Archive project, which is a collaboration between the Canadian Astronomy Data Centre (CADC), Victoria, and the James Clerk Maxwell Telescope (JCMT), Hilo. This work also makes use of observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory (JPL), California Institute of Technology under a contract with NASA. Moreover, this research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Infrared Science Archive, which is operated by the JPL, California Institute of Technology, under contract with the NASA.

References

All Tables

All Figures