© ESO, 2014

1. Introduction

The fundamental role of massive stars in shaping their direct environment and the galaxies they are in, makes the understanding of the process through which they form one of the major objectives yet to be achieved by astrophysicists. Progress in this sense requires a detailed observational knowledge of high-mass star forming regions from a physical and chemical point of view. Molecular line and dust continuum emission in the submm regime are among the main tools to study the first stages of massive star formation, where the sources are still embedded in the molecular gas.

In the low-mass regime, molecules such as CO and CS tend to freeze onto the dust grains in the densest and coldest parts of starless cores (e.g. Caselli et al. 1999; Kramer et al. 1999; Bergin & Tafalla 2007; Caselli 2011). This is observed as an abundance drop for these molecules (referred to as depletion) in the central parts of the core, identified by means of dust continuum emission (e.g. Tafalla et al. 2002, and references therein). The evolution of the abundance of different molecules on large scales can be reproduced by time-dependent models for gas-phase chemistry (e.g. Langer et al. 2000). Comparing observations of starless cores in the mm-continuum and in molecular lines from (among other species) CO isotopologues, Tafalla et al. (2002) find that depletion of CO can be substantial, with abundances with respect to H₂ in the central regions that are up to 1−2 orders of magnitude below the canonical ones. CO depletion is a temperature- and density-sensitive process; at low temperatures and high densities the depletion is higher, because under those conditions it is easier for the molecules to attach themselves to the grains. When protostars are formed in the core, the temperature increases, and at temperatures T ~ 20 − 25 K (commonly found in high-mass clumps, e.g. Wienen et al. 2012) the molecules evaporate from grains back into the gas phase, and the abundance returns to canonical levels (for CO ~ 10^-4 with respect to molecular hydrogen, in the solar neighbourhood; see e.g. Fontani et al. 2006). The depletion in a molecular core can thus vary substantially during the process of star formation, and can be used as an evolutionary indicator.

The Atacama Pathfinder Experiment (APEX¹, Güsten et al. 2006) is a 12 m submm telescope located on Chajnator plane in Chile. The APEX Telescope Large Area Survey of the Galaxy (ATLASGAL, Schuller et al. 2009) is the first complete survey of the inner Galactic plane, carried out in the submm continuum at 870 μm (345 GHz), and thus sensitive to very cold dust, potentially tracing the pristine condensations where the process of star formation still has to begin. Contreras et al. (2013) compiled a catalogue of compact sources (extracted with SExtractor, Bertin & Arnouts 1996) of the ATLASGAL survey, and Csengeri et al. (2014) filtered the extended emission and decomposed the remainder into Gaussian components using the GAUSSCLUMPS algorithm (Stutzki & Güsten 1990), optimised to extract the properties of embedded sources. The properties of a sample of massive star forming clumps in the ATLASGAL survey were studied by Urquhart et al. (2013a) and Urquhart et al. (2013b), searching for sources associated with methanol masers (from the methanol multi-beam survey, Green et al. 2009, 2010, 2012; Caswell et al. 2010, 2011) and (ultra-)compact Hii regions (from the CORNISH survey, Hoare et al. 2012; Purcell et al. 2013), respectively.

On the other hand, the starless-phase in the high-mass regime is elusive, due to its short duration (e.g. Motte et al. 2007; Russeil et al. 2010; Tackenberg et al. 2012). Candidates have been found looking for compact and massive molecular condensations without signs of active star formation (e.g. Chambers et al. 2009; Rygl et al. 2010, 2013; Butler & Tan 2012; Giannetti et al. 2013; Beuther et al. 2013). Recent Herschel observations have also been used to search for massive starless objects (Nguyen Luong et al. 2011). These clumps show lower temperatures than their actively star forming counterparts (10 − 15 K vs. 20 − 40 K). The latter are characterised by tracers of current star formation such as PAH emission, radio continuum and maser emission (see above).

A number of studies, usually targeting a limited sample of sources, have been carried out to investigate depletion in high-mass clumps (e.g. Zinchenko et al. 2009; Miettinen et al. 2011; Hernandez et al. 2011; Rygl et al. 2013; Liu et al. 2013), and very few address the variation with time of this essential parameter for the chemistry of the source (Fontani et al. 2012). Discordant evidence for depletion exists for high-mass objects, with both claims of significant CO freeze-out onto grains (Hernandez et al. 2011; Fontani et al. 2012; Rygl et al. 2013) and of canonical gas-phase abundances (Zinchenko et al. 2009; Miettinen et al. 2011). In this work we investigate the CO abundance in a large sample of massive clumps, selected from ATLASGAL, by means of its rarer (optically thin) isotopologues. The clumps have been selected to be in different evolutionary phases, so that we can also study changes in CO abundance in massive clumps during their evolution.

The paper is organised as follows. In Sect. 2 we describe the sample and its selection; in Sect. 3 we briefly summarise the observations. In Sect. 4 we derive the physical properties of the clumps, and the carbon and oxygen isotopic ratios across the inner Galaxy; Sect. 5 is dedicated to the discussion of these results, to the derivation of the CO-depletion factor (f_D = X^E/X^O, where X^O and X^E are the observed and expected abundances, respectively), and to the study of the physical and chemical changes occurring in the sources as they evolve. Finally, in Sect. 6 we draw some general conclusions and summarise our findings.

2. The sample

The ATLASGAL survey constitutes an excellent sample for selecting massive clumps. To select sources in various evolutionary phases, the clumps were extracted from the ATLASGAL database with the following criteria, each one defining a group of objects:

• The 32 brightest sources of the whole survey, ex-cluding the Central Molecular Zone (i.e. the centralfew × 10² pc) are called IRB group. These clumps are also detected in the IR;

• The 25 brightest sources, that are classified as a Massive Young Stellar Object (MYSO) in the Red MSX Sources survey (Urquhart et al. 2008), hereafter referred to as the RMS group;

• The 23 brightest objects dark at 8 μm constitute the D8 group;

• Finally, the 22 brightest sources that are dark at 24 μm, hereafter called D24.

Above, “brightest” refers to the submm peak flux at 870 μm. The members of each group were selected from the ATLASGAL catalogue after removing those of the preceding groups. Sources are classified as 24 μm or 8 μm dark if their average 24 μm or 8 μm flux density within the ATLASGAL beam is smaller than the average flux at the same wavelength in their vicinity.

This sample is termed the TOP100 sample of the ATLASGAL survey. Clumps in the first two groups of objects most likely contain high-mass stars and/or massive young stellar objects, while the other two are in an earlier phase of evolution (see Motte et al. 2007; Nguyen Luong et al. 2011), hosting very deeply embedded (and thus heavily extincted) (proto-)stars or being still starless. Separating the sample in these four categories allows us to study possible variations of depletion as a function of evolution.

This is the first of a series of papers concerning follow-up observations of the ATLASGAL survey.

3. Observations

Single-pointing observations were carried out for all sources in the TOP100 sample with APEX/FLASH (First-light APEX submillimeter heterodyne instrument; Heyminck et al. 2006) in C¹⁷O(3 − 2), on 14, 15, 24 June 2011 and 11−12 August 2011. Lower-excitation lines have been observed with APEX-1 for southern sources (¹³CO, C¹⁸O J = 2 − 1; observations carried out on 13−15 November 2008, and on 30 October and 1−2 November 2009), and with the EMIR (Eight mixer receiver) on the IRAM 30-m telescope for northern sources (C¹⁸O, ¹³CO, C¹⁷O, ¹³C¹⁸O J = 1 − 0; on 8−11 April 2011). The frequency of each transition and the beam size of the antenna for that molecular line are listed in Table 1. We have divided our sample into three subsamples, named S1, S2 and S3, according to the observational data we had available. Their positions and basic characteristics are listed in Tables A.1–A.3. Subsample S1 consists of 35 sources with EMIR follow-ups. In subsample S2 there are 44 sources, with APEX-1 follow-ups and subsample S3 (23 sources) were observed only with FLASH. Single Gaussians were fitted to the spectra for ¹³CO, C¹⁸O and ¹³C¹⁸O lines, while for C¹⁷O(3 − 2) and C¹⁷O(1 − 0) the appropriate hyperfine splitting was taken into account. The line parameters are listed in Tables A.4 to A.9. Figures from B.1 to B.5 show the spectra for all the observed sources; an example is given in Fig. 1, showing C¹⁷O(1 − 0), C¹⁸O(1 − 0), ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$ and C¹⁷O(3 − 2) for AGAL010.624-00.384. The velocity resolution and the rms noise of the spectra are typically between 0.3 − 0.5 km s^-1 and 0.05 − 0.10 K (in units of main beam temperature, T_MB), respectively.

4. Results

4.1. Distance

In order to derive the masses and to investigate the relative abundance of different isotopes as a function of Galactocentric distance, we searched the literature for a distance determination for each of the sources in our sample. We use direct maser parallax or spectrophotometric measurements where available (7 and 16 sources respectively), otherwise (79 sources) we use the kinematic distance (using the rotation curve of Brand & Blitz 1993), resolving the near-far ambiguity with Hi self-absorption data from Wienen et al. (in prep.), or from the literature. After this process, only two sources out of 102 still have a distance ambiguity, and for these we use the near kinematic distance. On the other hand, no such ambiguity is present for the Galactocentric distance. Our findings are summarised in Tables A.1−A.3. We choose to use the Brand & Blitz (1993) instead of the more recent one of Reid et al. (2009) because it still the best sampled in terms of Heliocentric and Galactocentric distances and Galactocentric azimuth.

4.2. Excitation temperatures

The excitation temperature (T_ex) is derived from the ratio R_ij of the integrated line intensities of the transitions i → i − 1 and j → j − 1, ${\mathit{R}}_{\mathit{ij}}\mathrm{\equiv }\frac{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}\mathit{,}\mathit{i}\mathrm{\to }\mathit{i}\mathrm{-}\mathrm{1}}^{}\mathrm{d}\mathit{V}}{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}\mathit{,j}\mathrm{\to }\mathit{j}\mathrm{-}\mathrm{1}}^{}\mathrm{d}\mathit{V}}\mathrm{=}\frac{\mathit{I}\mathrm{(}\mathit{i}\mathrm{\to }\mathit{i}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{I}\mathrm{(}\mathit{j}\mathrm{\to }\mathit{j}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(1)where I(i → i − 1) stands for the integrated intensity of the molecular transition i → i − 1 and the main beam temperature is ${\mathit{T}}_{\mathrm{MB}}\mathrm{=}\mathit{\eta }\mathrm{\left[}\mathit{J}\mathrm{\right(}{\mathit{T}}_{\mathrm{ex}}\mathrm{)}\mathrm{-}\mathit{J}\mathrm{(}{\mathit{T}}_{\mathrm{BG}}\mathrm{\left)}\mathrm{\right]}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}\mathit{,}$(2)with $\mathit{J}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathrm{B}}}\frac{\mathrm{1}}{\mathrm{(}{\mathrm{e}}^{\mathit{h\nu }\mathit{/}\mathrm{\left(}{\mathit{k}}_{\mathrm{B}}\mathit{T}\mathrm{\right)}}\mathrm{-}\mathrm{1}\mathrm{)}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(3)In Eq. (2), η is the source filling factor, T_BG is the background temperature and τ is the optical depth. To solve Eq. (1) for T_ex we neglected the background contribution J(T_BG), and assumed that the emission is optically thin, and that it traces the same volume of gas for both transitions. From the latter assumption it follows that the filling factor η cancels out in the ratio, if the angular resolution is similar, and that the lines have similar profiles. For subsample S1 R_ij, and thus T_ex, is obtained from C¹⁷O(3 − 2) and C¹⁷O(1 − 0), while for subsample S2 it is obtained from C¹⁷O(3 − 2) and C¹⁸O(2 − 1), assuming an isotopic ratio [¹⁸O]/[¹⁷O] = 4 (see Sect. 4.3) and correcting R_ij for the optical depth of C¹⁸O(2 − 1) (see Sect. 4.3). The values derived for T_ex range from ~ 5 K to ~ 70 K and are listed in Tables A.10 and A.11, with their uncertainties for individual sources. Since we do not have maps, we cannot smooth molecular data to a common resolution to derive the temperature, removing the filling factor from the ratio. This is not a problem for subsample S1, because C¹⁷O(3 − 2), observed with APEX, and C¹⁷O(1 − 0), observed with the 30-m, have similar angular resolutions (19′′ − 21′′), but it may be an issue for subsample S2 (19′′ − 28′′). As an example, if one assumes a source size of 70″ and tries to account for the different angular resolutions using the correction factor $\mathrm{(}{\mathit{\theta }}_{\mathrm{beam}}^{\mathrm{2}}\mathrm{+}{\mathit{\theta }}_{\mathrm{source}}^{\mathrm{2}}\mathrm{)}\mathit{/}{\mathit{\theta }}_{\mathrm{source}}^{\mathrm{2}}$, the ratio R_ij varies by ~ 10%. This implies a decrease in T_ex of ≲10% for a T_ex ≲ 20 K and up to 25% for T_ex ~ 70 K. For smaller sources the difference increases and reaches maximum for point sources. To this source of uncertainty, one has to add that the assumed value of the isotopic ratio [¹⁸O]/[¹⁷O] = 4 may not be appropriate for specific sources (cf. Fig. 6b). The highest excitation temperatures are found in S2, some of them possibly caused by the uncertainties discussed above. Equation (1) is implemented in JAGS² (Just Another Gibbs Sampler) to estimate T_ex and its uncertainty, directly from the measured quantities. In the procedure, only temperatures up to ~ 100 K were considered to be physically plausible. This is also motivated by the transitions used, for which the method is increasingly insensitive above ~ 30 K.

For S3 we used the average values of T_ex from S1 and S2, listed in Table 2, for the appropriate group of sources to derive the ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$ column density, and the mass from the submm continuum.

4.3. Isotopic abundance variations in the Galaxy

To derive more accurate optical depths, in Sect. 4.4, we start from a first estimate of the relative abundance of the different isotopes of C and O as a function of the Galactocentric distance D_GC in the range 2 − 8 kpc. The isotopic gradients can then be used as priors for the procedure to derive τ (see Bolstad 2007, for a description of the Bayesian approach to statistics). This can be achieved by comparing the column densities of different carbon monoxide isotopologues, for which we use the (1 − 0) transition of C¹⁸O, C¹⁷O and ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}$. If we assume that the transitions are optically thin, the integrated flux of a transition is proportional to the column density, and we can approximate the relative abundance of C or O isotopes with the ratio of the integrated fluxes of the lines, if we are using the same rotational transition (i.e. with similar excitation conditions, and likely tracing the same volume of gas). Using the same rotational transition thus implies that we do not need to know the actual value of T_ex, which is the same for both isotopologues. The ratio of the fluxes was corrected for the small differences in the transition frequency (see Linke et al. 1977). Panel (a) in Fig. 2 shows that we find no clear trend with D_GC for [¹⁸O]/[¹⁷O], and our values of the ratio show a large scatter. A constant value of this ratio in space and time is to be expected if ¹⁸O is a secondary nucleosynthesis product (those that can be produced only in the presence of pre-existing seed nuclei, generated by previous stellar generations), like ¹⁷O, and they both come from the same primary element, ¹⁶O (Wilson & Rood 1994). On the other hand, Wouterloot et al. (2008) find that the [¹⁸O]/[¹⁷O] ratio tends to increase in the outer Galaxy, out to D_GC ~ 16 kpc.

We have fit the formula $\mathit{N}{\mathrm{(}}^{}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}\mathrm{)}\mathit{/}\mathit{N}{\mathrm{(}}^{}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}\mathrm{)}\mathrm{=}\mathit{\alpha }{\mathit{D}}_{\mathrm{GC}}\mathrm{+}\mathit{\beta }$ to the data (see Fig. 2), selecting sources with low optical depth (<0.3, i.e. a correction of ~ 15% in column density; see Eq. (5)), derived from the detection equation alone. The fit gives a small negative slope, $\mathrm{~}\mathrm{-}\mathrm{0.}{\mathrm{1}}_{-0.2}^{\mathrm{+}\mathrm{0.1}}$. The slope is still consistent with zero and it is so small that the resulting values of [¹⁸O]/[¹⁷O] would change by only ~ 1 (25% of the assumed value of 4) over the whole range of D_GC of the S1 sample. This variation is of the order of intrinsic scatter of the fit ($\mathrm{~}\mathrm{1.}{\mathrm{0}}_{-0.2}^{\mathrm{+}\mathrm{0.2}}$). We note that a model with only one free parameter (the intrinsic scatter, and [¹⁸O]/[¹⁷O] = 4 across the inner Galaxy) is to be favoured with respect to a model with three free parameters (the slope, the intercept and the intrinsic scatter): a Bayesian model comparison shows that the odds ratio (see e.g. Gregory 2005) is about 15 in favour of the first model. This thanks to the fact that more complex models are automatically penalised in the Bayesian approach, and are to be preferred only if the data justify the added complexity (Occam’s Razor). Therefore, we assume an average constant value for [¹⁸O]/[¹⁷O] = 4, independently of D_GC, with an intrinsic scatter as given by the fit; this value of the ratio is consistent with the measurements of Wouterloot et al. (2008) and Wilson & Rood (1994) for the same range of D_GC.

We also investigate the variation of the [¹²C]/[¹³C] ratio as a function of D_GC, derived from both I [¹²C¹⁸O(1 − 0)]/ I [¹³C¹⁸O(1 − 0)] and I [¹²C¹⁷O(1 − 0)] /I [¹³C¹⁸O(1 − 0)]. The former gives a direct estimate of the [¹²C]/[¹³C] ratio, but it may be affected by optical depth issues, while the latter involves only optically thin transitions, but the ratio [¹⁸O]/[¹⁷O], despite not being dependent on the Galactocentric distance, may vary up to ~ 50% in different sources, thus increasing the scatter in the relation. The measured ratios and the results of fitting a straight line to the data points (using as priors the results reported by Milam et al. 2005, see below) are shown in Figs. 2b and c. From these panels, one can see that [¹²C]/[¹³C] increases with D_GC, in agreement with previous work (Langer & Penzias 1990; Wilson & Rood 1994), and expected for a primary/secondary nucleosynthesis product ratio (see Wilson & Rood 1994); in fact, primary elements are produced independently of the metallicity of the stars, while the production of secondary elements increases for higher metallicities. We find the relations ${\mathrm{[}}^{\mathrm{12}}\mathrm{C}\mathrm{]}\mathit{/}{\mathrm{[}}^{\mathrm{13}}\mathrm{C}\mathrm{]}\mathrm{=}\mathrm{6.}{\mathrm{2}}_{-2.1}^{\mathrm{+}\mathrm{1.1}}{\mathit{D}}_{\mathrm{GC}}\mathrm{+}\mathrm{9.}{\mathrm{0}}_{-6.2}^{\mathrm{+}\mathrm{9.9}}$ and $\mathrm{\left(}{\mathrm{[}}^{\mathrm{12}}\mathrm{C}\mathrm{\right]}\mathit{/}{\mathrm{[}}^{\mathrm{13}}\mathrm{C}\mathrm{\left]}\mathrm{\right)}\mathit{/}\mathrm{\left(}{\mathrm{[}}^{\mathrm{18}}\mathrm{O}\mathrm{\right]}\mathit{/}{\mathrm{[}}^{\mathrm{17}}\mathrm{O}\mathrm{\left]}\mathrm{\right)}\mathrm{=}\mathrm{1.}{\mathrm{6}}_{-0.3}^{\mathrm{+}\mathrm{0.8}}{\mathit{D}}_{\mathrm{GC}}\mathrm{+}\mathrm{1.}{\mathrm{1}}_{-1.4}^{\mathrm{+}\mathrm{1.9}}$, respectively from I [¹²C¹⁸O(1 − 0)] /I [¹³C¹⁸O(1 − 0)] and I [¹²C¹⁷O(1 − 0)] /I [¹³C¹⁸O(1 − 0)]. We note that the uncertainties reported in the relations are derived marginalising (see e.g. Jaynes & Bretthorst 2003; Bolstad 2007; Gregory 2005) the probability distributions shown in the right panels of Figs. 2b and c, but the values of the parameters of the fit α and β (y = αx + β) are strongly correlated and the yellow-shaded area in the left panels of the figure, representing the fit uncertainty, takes this into account. The above fit results give a [¹²C]/[¹³C] ratio for the solar neighbourhood of ~ 62 and ~ 60 respectively, the latter derived with a fixed [¹⁸O]/[¹⁷O] of 4. In this case as well source-to-source variations are found, resulting in an intrinsic scatter of $\mathrm{~}{\mathrm{8}}_{-4}^{\mathrm{+}\mathrm{3}}$ (~ 15% of the [¹²C]/[¹³C] at the position of the Sun), obtained from the fit procedure. The estimated scatter in isotopic ratios for sources at a given D_GC can be caused by a multiplicity of processes at work in the specific source (such as chemical fractionation and selective photodissociation) differently affecting the specific molecules, from our simple estimate of τ, or can be intrinsic if the metallicity and the star formation history are different or other processes have an important role in modifying the isotopic ratios (e.g. non-efficient mixing, radial mixing, cloud mergers). Milam et al. (2005) use three different molecular species, CN, CO, and H₂CO, to investigate the [¹²C]/[¹³C] ratio across the Galaxy, finding that all species give consistent results. The authors use this result to show that photodissociation does not have systematic and strong effects on the isotopic ratios. In addition, they compare the data with the predictions from a simple model for chemical fractionation, showing that they are at odds, thus concluding that also this process does not play a fundamental role. In the small panel in Fig. 2b we also show the average values of [¹²C]/[¹³C] reported by Milam et al. (2005) as black open circles, with the error bars showing the range in the measured [¹²C]/[¹³C] from the different species (or the given uncertainty where only one measurement is available), along with our points (grey open triangles). The comparison illustrates that the points are fully consistent with our measurements and fit. The fact that we find a smaller intercept is mainly due to the inclusion of Sgr B2 (D_GC = 0.1 pc) and WB189 (D_GC = 16.4 pc) in the work of Milam et al. (2005) and also to the unweighted fit that they use; the points with a high [¹²C]/[¹³C] and a large uncertainty (cf. their Fig. 2) could slightly increase the intercept.

4.4. Optical depths and column densities

Assuming Local Thermodynamic Equilibrium (LTE), we can compute the column density of the molecules using the excitation temperature, according to $\mathit{N}\mathrm{=}\mathit{C}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathit{f}\mathrm{\left(}{\mathit{T}}_{\mathrm{ex}}\mathrm{\right)}\mathrm{\int }{\mathit{T}}_{\mathrm{MB}}\mathrm{d}\mathit{V,}$(4)where τ is the optical depth of the transition, and f(T_ex) is a term grouping together all the constants and the terms depending on T_ex (see e.g. Kramer & Winnewisser 1991). For a precise determination of the column density we have to estimate the optical depth τ of the lines. The correction factor C(τ) for Eq. (4) can be expressed in the form (for τ ≲ 2Goldsmith & Langer 1999): $\mathit{C}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\frac{\mathit{\tau }}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}}\hspace{0.17em}\mathrm{·}\hspace{0.17em}$(5)Typically C¹⁸O and C¹⁷O have column densities so low that their emission is optically thin, meaning that the correction for τ is small and therefore the molecular column density is proportional to the integrated flux of the lines. However, it is worthwhile to check if the emission in environments having high column densities (as we expect for the clumps in the TOP100 sample) may still be assumed to be optically thin. The optical depth of a transition can be estimated by means of the detection equation or from the ratio of the integrated flux of the same transition, coming from different isotopologues, if the relative abundance is known (e.g. Hofner et al. 2000), ${\mathit{R}}_{\mathit{ij}}\mathrm{=}\frac{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}\mathit{,}\mathit{i}}^{}\mathrm{d}\mathit{V}}{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}\mathit{,j}}^{}\mathrm{d}\mathit{V}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathit{i}}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\varphi }{\mathit{\tau }}_{\mathit{i}}}}\mathit{,}$(6)where ϕ is the relative abundance of the two isotopologues (X_j/X_i). To derive the optical depth we used different methods, depending on the data available.

Subsample S1: for these sources we have several CO isotopologues observed in the same transition, namely ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}$, ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$ and ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$. We derive and make a fit to the probability distribution of the three line ratios by means of Monte Carlo Markov chains. The results of the fits to the probabilities are used to describe the distribution out of which the measured ratios are extracted. The τ of the lines of the different isotopologues are interconnected by means of

$\begin{array}{ccc}& & \\ & & \\ & & {\mathit{\tau }}_{\mathrm{13}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathrm{=}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}}\mathit{A}\mathit{/}\mathit{B}\mathrm{\equiv }{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}}\mathit{/}\mathit{C,}\end{array}$

where A and B are the relative abundances [¹⁸O]/[¹⁷O] and [¹²C]/[¹³C], respectively. The relations

$\begin{array}{ccc}& & \frac{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}\mathrm{\equiv }{\mathit{R}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathit{/}\mathit{A}}}\mathit{,}\\ & & \frac{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}\mathrm{\equiv }{\mathit{R}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathit{/}\mathit{B}}}\mathit{,}\end{array}$ $\begin{array}{ccc}& & \frac{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}{\mathrm{\int }{\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{}\mathrm{d}\mathit{V}}\mathrm{\equiv }{\mathit{R}}_{\mathrm{3}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathit{/}\mathit{A}}}{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathit{/}\mathit{B}}}\mathit{,}\\ & & {\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{}\mathrm{=}\mathit{\eta }\mathrm{\left[}\mathit{J}\mathrm{\right(}{\mathit{T}}_{\mathrm{ex}}\mathrm{)}\mathrm{-}\mathit{J}\mathrm{(}{\mathit{T}}_{\mathrm{BG}}\mathrm{\left)}\mathrm{\right]}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}}\mathrm{)}\mathit{,}\mathrm{and}\\ & & {\mathit{T}}_{\mathrm{MB}{\mathit{,}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}}^{}\mathrm{=}\mathit{\eta }\mathrm{\left[}\mathit{J}\mathrm{\right(}{\mathit{T}}_{\mathrm{ex}}\mathrm{)}\mathrm{-}\mathit{J}\mathrm{(}{\mathit{T}}_{\mathrm{BG}}\mathrm{\left)}\mathrm{\right]}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}\mathit{/}\mathit{A}}\mathrm{)}\end{array}$

were simultaneously solved with Monte Carlo Markov chains in JAGS, thus combining the information from the detection equation and the line ratios for a more accurate determination of τ and of the isotopic ratios. Hofner et al. (2000) have maps in C¹⁷O(2 − 1) for a sample of similar sources at a comparable angular resolution, and find them to be extended; therefore, for simplicity we assume η = 1 for these calculations. The priors for A and B are Gaussian curves, centred on the value expected from the fit of the isotopic ratios as a function of D_GC, with a dispersion equal to the intrinsic scatter of the fit (see Sect. 4.3); the prior for τ comes from the fit of the hyperfine structure of the C¹⁷O(1−0) emission. This procedure returns ${\mathit{\tau }}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}}$, ${\mathit{\tau }}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}}$, ${\mathit{\tau }}_{{}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}}$, A, B and C, and their respective uncertainties for each source in the subsample S1, refining our estimate of the isotopic ratios and of the optical depth. We use these values of the [¹⁸O]/[¹⁷O] and [¹²C]/[¹³C] isotopic ratios in Sect. 5.3 to derive a more accurate f_D, and in Sect. 5.2 to derive again the gradient of relative abundance as a function of D_GC (see Fig. 6). The results of this procedure for the refined estimate of the [¹²C]/[¹³C] and [¹⁸O]/[¹⁷O] ratios are discussed in Sect. 5.2 and summarised in Table A.12. We note that a significant fraction of the sources (~ 55%) shows an optical depth >0.3 for ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$. On the other hand, only one source has ${\mathit{\tau }}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}}\mathit{>}\mathrm{0.3}$, while τ ≲ 0.2 for the other sources; similar low values for the ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$ opacity were obtained by Hofner et al. (2000).

Subsample S2: for this subsample we have observations of ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{16}}\mathrm{O}$ and ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}$. Thus, we derive τ from the ratio of these two transitions. To estimate the uncertainty in the optical depth, we calculate the probability distribution of the line ratio, considering the calibration uncertainty and the noise of the spectrum. This can be translated into a distribution of τ for ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}$, using the appropriate value for [¹²C]/[¹³C] and [¹⁸O]/[¹⁶O] at the D_GC of the source. The results are listed in Table A.8.

Subsample S3: finally, for the last subsample, the only measure for τ we have comes from the fit of the hyperfine structure of ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}$. However, the satellites are too close to be spectrally resolved, meaning τ is highly uncertain. Thus, we assume that the emission is optically thin.

The column densities of the CO isotopologues are listed in Tables A.13–A.15.

4.5. Column density of molecular hydrogen

The column density of molecular hydrogen N(H₂) was calculated from the ATLASGAL peak flux, according to Schuller et al. (2009), $\mathit{N}\mathrm{\left(}{\mathrm{H}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\mathit{\gamma }\frac{{\mathit{F}}_{\mathrm{p}}}{\mathrm{2.8}\hspace{0.17em}{\mathit{m}}_{\mathrm{p}}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{B}}\hspace{0.17em}{\mathit{\kappa }}_{\mathrm{870}}\hspace{0.17em}{\mathit{B}}_{\mathrm{870}}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}}\mathit{,}$(15)where F_p is the peak flux at 870 μm, the factor 2.8 accounts for the presence of helium (~ 10% in number, e.g. Allen 1973), κ₈₇₀ = 1.8 cm² g^-1 is the dust opacity at 870 μm (derived from the numbers given in Ossenkopf & Henning 1994), m_p is the proton mass, Ω_B is the beam solid angle for a beam size of ≈19.2″ (i.e. 9.817 × 10^-9 sr), B₈₇₀(T_d) is the emission at 870 μm of a blackbody with a dust temperature T_d and γ is the gas-to-dust ratio, assumed to be 100. To derive the column density of molecular hydrogen we assume that gas and dust are coupled, and that the molecules are in LTE, thus T_d = T_K = T_ex (for a comparison of independently derived T_K and T_d for a similar source sample Giannetti et al. 2013). For the subsample S3, we do not have a direct estimate of the temperature. Thus, we assign the typical temperature of the group to which the object belongs, with an uncertainty equal to the dispersion of that group (see Sect. 5 and Table 2).

4.6. Masses

The clump mass is derived from the integrated 870 μm flux as given in the latest ATLASGAL catalogue (for details, see Csengeri et al. 2014), through $\mathit{M}\mathrm{=}\frac{\mathit{\gamma }{\mathit{S}}_{\mathrm{870}}{\mathit{D}}^{\mathrm{2}}}{{\mathit{\kappa }}_{\mathrm{870}}{\mathit{B}}_{\mathrm{870}}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathrm{\right)}}\mathit{,}$(16)where S₈₇₀ is the 870 μm integrated flux, and D is the distance. The objects are massive, ranging from ~ 10²M_⊙ to ~ 3 × 10⁴M_⊙: we are dealing with extreme sources in the ATLASGAL survey, many of which are known sites of high-mass star formation. For subsample S3 the same assumptions for the temperature are made as for the determination of column densities. The masses are listed in Tables A.16 to A.18.

Kauffmann & Pillai (2010) derived a relation between mass and radius, to separate clumps with the potential of forming massive stars from those without. A similar threshold for high-mass star formation, with a constant surface density Σ = 0.05 g cm^-2, was derived by Urquhart et al. (2013a) for massive star forming clumps selected cross-correlating the ATLASGAL catalogues with the catalogue of the methanol multi-beam survey and the CORNISH survey. Figure 3 shows that all of the clumps in this sample are found above this empirical relations, implying that these clumps are potentially forming massive stars. In the figure, we used the R_eff definition from Rosolowsky et al. (2010) and the sizes of the two-dimensional Gaussian fitted to the emission in Csengeri et al. (2014). With this definition, R_eff is approximately equal to the FWHM size of the clump.

5. Discussion

5.1. Linewidths and temperatures

Figure 4 shows the kernel density estimates (e.g. Wand & Jones 1995) applied to the linewidths of ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}$, for the four groups of sources. We used the KernSmooth package and an Epanechnikov kernel. The bandwidth was derived with the direct plug-in method (Sheather & Jones 1991). In practice, the procedure convolves the discrete data points with a kernel function, obtaining a smooth probability density function. The typical linewidth of the ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$ lines for IRB is much larger than that of the other groups. The less evolved RMS and D8 sources show similar values of ΔV, slightly larger than those of D24 and smaller than those typical of the IRB. This may indicate that (proto-)stars embedded in the molecular gas and dust deliver energy to the surrounding material, making it more turbulent than before their formation, which is reflected in the width of molecular lines. Alternatively, the IR-bright clumps may have been more turbulent from the start, and move much faster through the IR-dark phase, where we find very few clumps with ΔV ≳ 4 km s^-1.

The kernel density estimates for T_ex, shown in Fig. 5, are obtained in the same way as those for the linewidths. The distributions show that the four groups of objects have a different characteristic temperature, listed in Table 2. This indicates an increasing temperature as evolution proceeds. A similar result for gas and/or dust temperature is found in e.g. Rygl et al. (2010), Wienen et al. (2012) and Giannetti et al. (2013).

5.2. Refined estimate of the isotopic ratios

Figure 6 shows the refined estimate of the isotopic ratios for the sources in S1; we also included different velocity components observed in the spectra (indicated with crosses), with the appropriate estimate of the Galactocentric distance obtained with the rotation curve of Brand & Blitz (1993). First of all, both the isotopic ratios [¹⁸O]/[¹⁷O] and [¹²C]/[¹³C] show a large scatter for sources at similar D_GC, confirming the existence of source-to-source variations in the relative abundance of these isotopes at the same D_GC.

Here we derive the [¹⁸O]/[¹⁷O] ratio in massive clumps and find it to be consistent with the current estimate of ~ 4 for 2 kpc ≲ D_GC ≲ 8 kpc, albeit with a very large scatter that cannot be due to the uncertainty associated with the inferred values of this ratio. The [¹⁸O]/[¹⁷O] ratio is found to range between ~ 2 and ~ 6, as can be seen in Fig. 6b. A few sources have [¹⁸O]/[¹⁷O] ~ 5 − 6, especially objects in groups D8 and D24, values similar to that found in pre-solar grains ([¹⁸O]/[¹⁷O] ~ 5.5, reported in e.g. Prantzos et al. 1996). Prantzos et al. (1996) discuss how ¹⁸O poses a problem for models of chemical evolution of the Galaxy. Neither [¹⁸O]/[¹⁶O] nor [¹⁸O]/[¹⁷O] can be reproduced by simple models describing the chemical evolution of the Milky Way. The former because it appears to be larger now in the solar vicinity compared to the time of the Sun’s formation, whereas it is predicted to decrease with time; the latter because it appears to be substantially larger in the pre-solar cloud than in the present-day ISM, while it is predicted to remain constant (e.g. Prantzos et al. 1996). The most probable solution to this discrepancy is the pollution of the pre-solar cloud by a previous generation of massive stars undergoing Type II Supernovae (SNe) events, as discussed in detail by Young et al. (2011) (for [¹⁸O]/[¹⁶O] a role could also have been played by the position of the Sun at the time of its formation; see Wielen & Wilson 1997).

Young et al. (2011) show that values as high as those found in pre-solar grains can be reached if intermediate-mass Type II SNe events (from B stars) of a previous stellar generation pollute the surrounding medium, arguing that more massive progenitors are not able to produce similar values of the ratio [¹⁸O]/[¹⁷O]. It could be that the sources with a high [¹⁸O]/[¹⁷O] ratio had an evolution similar to that of the solar system in its earliest formation stages. Because several of these sources are near the inner molecular ring, many massive stars were formed in their vicinity, making this pollution scenario a reasonable one.

The procedure outlined in Sect. 4.4 also yields values of B (=[¹²C]/[¹³C], see Eq. (8)) for those sources without a ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$ detection (panel (c) of Fig. 6): we used a flux for the undetected ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}$ transition equal to 1σ with σ derived from the spectrum as an input for the Monte Carlo calculations. Some of these sources have quite high values of [¹²C]/[¹³C] given their D_GC, even if the uncertainties are large. Such large values of [¹²C]/[¹³C] may be caused by simple dilution effects if the area of ${}^{\mathrm{13}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$ emission is substantially smaller than those of ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$ and ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$, or if not, it may be caused for instance by selective depletion of the heavier molecule. Excluding these non-detections and using the values for the carbon isotopic ratio from Sect. 4.4 we find a gradient of ${\mathrm{[}}^{\mathrm{12}}\mathrm{C}\mathrm{]}\mathit{/}{\mathrm{[}}^{\mathrm{13}}\mathrm{C}\mathrm{]}\mathrm{=}\mathrm{6.}{\mathrm{1}}_{-1.8}^{\mathrm{+}\mathrm{1.1}}{\mathit{D}}_{\mathrm{GC}}\mathrm{+}\mathrm{14.}{\mathrm{3}}_{-7.2}^{\mathrm{+}\mathrm{7.7}}$, with an intrinsic scatter of $\mathrm{10.}{\mathrm{1}}_{-2.5}^{\mathrm{+}\mathrm{2.0}}$; this relation gives a most probable value for [¹²C]/[¹³C] of ~ 66 ± 12 for the solar vicinity. The relation derived is consistent with that derived from simpler estimates, but each point has reduced uncertainties (cf. Figs. 2b and 6), thanks to the better estimate of τ and the combination of all the line ratios. We note that the slope of the relation does not vary much, but the intercept has increased slightly, probably because of the correction for optical depth effects. This new estimate highlights even more the large scatter of the points.

No clear systematic difference in the isotopic ratios is found between the groups. This indicates that the isotopic ratios set by the complex interplay of the processes mentioned in Sect. 4.3 do not change appreciably on a large scale over the short time spanned by the first phases of high-mass star formation.

5.3. Depletion of CO

Figure 7 shows N(H₂) as a function of $\mathit{N}{\mathrm{(}}^{}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}\mathrm{)}$ and $\mathit{N}{\mathrm{(}}^{}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}\mathrm{)}$ for the TOP100 sample. The dashed line indicates the canonical abundance of the CO isotopologue in the solar neighbourhood (~ 1.7 × 10^-7 for C¹⁸O and ~ 4.2 × 10^-8 for C¹⁷O). This is a first, simple approach for investigating CO depletion: several sources seem to have abundances lower than the canonical one, suggesting that CO may indeed be frozen onto the surface of the dust grains. We observe that there appears to be a difference between groups D8 and D24 and groups IRB and RMS, with objects from the former two groups reaching markedly lower abundances.

Panels (a) and (b) of Fig. 8 show the CO depletion factor as a function of the excitation temperature. Here, the depletion factor, which is defined as the ratio of the expected (i.e. canonical) and observed abundance, is determined taking into account the CO abundance gradient with D_GC, deriving the expected ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{18}}\mathrm{O}$ abundance according to Miettinen et al. (2011, and references therein), ${\mathit{X}}_{\mathrm{12}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}}^{\mathrm{E}}\mathrm{\left(}{\mathit{D}}_{\mathrm{GC}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{X}}_{\mathrm{12}{\mathrm{C}}^{\mathrm{16}}\mathrm{O}}^{\mathrm{E}}\mathrm{\left(}{\mathit{D}}_{\mathrm{GC}}\mathrm{\right)}}{{\mathrm{[}}^{\mathrm{16}}\mathrm{O}\mathrm{]}\mathit{/}{\mathrm{[}}^{\mathrm{18}}\mathrm{O}\mathrm{]}\mathrm{\left(}{\mathit{D}}_{\mathrm{GC}}\mathrm{\right)}}\mathit{,}$(17)where ${\mathit{X}}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{16}}\mathrm{O}}^{\mathrm{E}}\mathrm{\left(}{\mathit{D}}_{\mathrm{GC}}\mathrm{\right)}$ and [¹⁶O]/[¹⁸O] (D_GC) are the ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{16}}\mathrm{O}$ and [¹⁶O]/[¹⁸O] Galactic gradients, respectively. In detail, the former is ${\mathit{X}}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{16}}\mathrm{O}}^{\mathrm{E}}\mathrm{\left(}{\mathit{D}}_{\mathrm{GC}}\mathrm{\right)}\mathrm{=}{\mathit{X}}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{16}}\mathrm{O}}\mathrm{(}\mathrm{8.5}\hspace{0.17em}\mathrm{kpc}\mathrm{)}{\mathrm{e}}^{\mathrm{1.105}\mathrm{-}\mathrm{0.13}{\mathit{D}}_{\mathrm{GC}}\mathrm{\left[}\mathrm{k}\mathrm{pc}\mathrm{\right]}}$, with ${\mathit{X}}_{{}^{\mathrm{12}}\mathrm{C}^{\mathrm{16}}\mathrm{O}}\mathrm{(}\mathrm{8.5}\hspace{0.17em}\mathrm{kpc}\mathrm{)}\mathrm{=}\mathrm{9.5}\mathrm{\times }{\mathrm{10}}^{-5}$, and the latter is [¹⁶O]/[¹⁸O] (D_GC) = 58.8D_GC [kpc] + 37.1. The expected ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$ was calculated with the [¹⁸O]/[¹⁷O] derived as described above for the subsamples S1 and S2 (see Sect. 4.4). Colder sources show larger CO depletion factors, suggesting that in less evolved sources CO in the gas phase is less abundant than in more evolved ones. Figure 9 shows that f_D (from ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}$, in this figure) also correlates with the peak column and volume densities (as also found by Liu et al. 2013, with a lower angular resolution), increasing for higher densities. The peak volume density was calculated according to n(H₂) = N(H₂)/(2R_eff). The average f_D decreases from the D24 to the IRB objects as visible in panels (c) and (d) of Fig. 8; the kernel density estimates for the depletion factors were derived as for the linewidths and T_ex (see Sect. 4.2). Mean and median values for f_D for each group are listed in Table 3. The vertical offset between the four groups in Fig. 9 may be caused by the depletion and evaporation timescales, that depend on density and temperature, which is different, on average, for each group. Denser objects are more prone to high-levels of molecular depletion, because the depletion timescale (i.e. the time after which depletion becomes important) decreases with increasing density, while for higher temperatures, molecules tend to evaporate more rapidly, as the evaporation timescale decreases for increasing temperatures. We note that CO abundances lower than the canonical one are observed towards warm protostellar envelopes (e.g. Fuente et al. 2012, and references therein); the reason for this is as yet unclear, and the authors suggest that conversion of CO in more complex molecules may happen on the grain surface and/or that the shocks and UV radiation illuminating the walls of the cavities produced by the protostellar outflow may be important factors in explaining the measured CO abundance.

That f_D decreases with age of the source is supported by Fig. 10, where we find a weak anti-correlation between the depletion factor and the L/M ratio, which we use as an indication of the time evolution of the source (e.g. Molinari et al. 2008; López-Sepulcre et al. 2011). We obtain monochromatic luminosities by cross-matching our sources with the MSX (Egan et al. 2003) and WISE (Wright et al. 2010) point-source catalogues (see Csengeri et al., in prep., for details); from these, we extrapolated the bolometric luminosities following Davies et al. (2011), based on the results of Mottram et al. (2011). In Fig. 10, using also the sources with upper and lower limits (for a total of 79 sources) on the IR flux, we get a Pearson correlation coefficient of −0.40, while if we exclude them, 44 objects remain, and we obtain a correlation coefficient of −0.34. This means that there is a probability of ~ 1.1% of getting similar results from an uncorrelated sample. The large scatter can be due to several reasons. The various distances to the targets imply that we probe different spatial scales, where the depletion factor may vary significantly, influencing the average value we derive. For sources in D8 and D24 the ratio between the depletion and the free-fall timescales may differ in different sources, and some of the less evolved ones may not have lived enough for the depletion to become important (see below), or have rather high temperatures (≳20 K) making the evaporation of molecules from the grains much faster (cf. Fig. 8). In addition, the luminosity estimate is uncertain due to the extrapolation from monochromatic luminosity to bolometric luminosity.

The observed trends could be caused by the assumption of T_ex = T_K = T_d. In order to check if the trends are real, we performed two different tests. In the case of Fig. 8 we also used different dust temperatures for the groups: T_d = 10 K for group D24, T_d = 18 K for group D8, and T_d = 25 K for groups IRB and RMS (cf. e.g, Giannetti et al. 2013; Sánchez-Monge et al. 2013). We then obtain an average f_D of ~ 6 for D24 and ~ 2 − 3 for the IRB. Calculating the average f_D in bins of T_ex shows that the depletion factor decreases only slightly or stays constant for T_ex ≲ 25 K (f_D ~ 7 − 5), while it decreases more rapidly for T_ex ≳ 25 K (down to ~ 2). On the other hand, to test the dependence of f_D on the density (Fig. 9), we selected the sources in a smaller range of T_ex (14 − 24 K, to remove extreme temperatures and have still ~ 20 sources) and recomputed the densities and depletion factors assuming a constant T_d = 20 K. Panel (c) of Fig. 9 shows that the trends are still visible. The symbol sizes are proportional to the distance of the sources, showing that the trends are not due to distance effects.

In recent works (Wienen et al. 2012; and in prep.) the ammonia (1, 1), (2, 2) and (3, 3) inversion transitions were observed towards a subsample of sources in the TOP100 sample. The observations were carried out with Effelsberg and Parkes telescopes, thus sampling a larger scale than our observations, though with a typical filling factor of 0.1 (Wienen et al. 2012). The temperatures derived from ammonia are typically higher than those calculated in this work for sources in groups D8 and D24, especially for the latter, with typical values ~ 20 K. We tried to reproduce the observed line intensities with RATRAN assuming a temperature of 20 K: the results are briefly discussed in Sect. 5.3.1.

For the sources in groups D8 and D24 we can derive the timescale for CO depletion τ_dep, using the expression given in Bergin & Tafalla (2007) and the peak volume density of molecular hydrogen. We find that τ_dep is in the range 10³ − 10⁵ yr, of the same order of magnitude as the free-fall timescale, which is a rough measure of the starless/IR-quiet clumps lifetime. It is interesting to note that three out of four sources with f_D< 2 indeed have the largest τ_dep/τ_ff of the sample (≳1.5). Tackenberg et al. (2012) and Csengeri et al. (2014) find that lifetimes of massive starless/IR-quiet clumps are of the order of 10⁴ − 10⁵ yr; we can compare this number with the timescale for depletion as a function of radius in clumps with an n(H₂) ∝ r^-1.5 density profile (approximately the mean value for the clumps in the ATLASGAL survey). In practice, assuming a mass for the clump unequivocally sets the volume density of molecular hydrogen at all radii. Using again the equation of Bergin & Tafalla (2007), one can assign a characteristic timescale for CO depletion at each r, using the above-cited n(H₂). We can thus derive an estimate of the size R_dep of the central depletion hole as the radius at which τ_dep matches the typical lifetime of massive starless clumps. As illustrated in Fig. 11, this procedure yields values of R_dep in the interval ~ 0.02 − 0.1 pc for a clump with a mass within 1 pc of ~ 550 M_⊙ and an age of 10⁴ yr and 10⁵ yr, respectively. The radius increases for larger masses. For a distance of 4 kpc a radius of 0.1 pc would be slightly in excess of 5″. Figure 11 also shows the quantity f_M, ${\mathit{f}}_{\mathit{M}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{clump}}}{{\mathit{M}}_{\mathrm{clump}}\mathrm{-}\mathit{M}\mathrm{(}\mathit{r}\mathrm{=}{\mathit{R}}_{\mathrm{dep}}\mathrm{)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{R}}_{\mathrm{dep}}\mathit{/}{\mathit{R}}_{\mathrm{eff}}{\mathrm{)}}^{\mathrm{1.5}}}\mathit{,}$(18)which is the ratio between the total mass and the mass outside the depletion radius. If we assume that after τ_dep(n(H₂)) all CO within R_dep is locked onto grains, this is also an alternative way of measuring depletion.

We note that assuming τ_dep = τ_ff, one can derive an expression for the critical density n_H,crit above which all CO is depleted, that depends only on the temperature T and on the mean grain cross section per H atom ⟨ σ_H ⟩, $\begin{array}{ccc}{\mathit{n}}_{\mathrm{H}\mathit{,}\mathrm{crit}}& \mathrm{=}& {\left(\frac{\mathrm{3}\mathit{\pi }}{\mathrm{32}\mathit{G\mu }{\mathit{m}}_{\mathrm{H}}}\right)}^{-1}\left(\frac{\mathrm{1}}{\mathrm{⟨}{\mathit{\sigma }}_{\mathrm{H}}\mathrm{⟩}\mathit{S}{\mathit{v}}_{\mathrm{th}\mathit{,}\mathrm{CO}}}\right)\\ & & \end{array}$(19)where the mean cross section per H atom is ⟨ σ_H ⟩ = σ_gn_g/n_H ≈ 10^-21 cm², v_th,CO is the thermal velocity for CO and S is the sticking probability. If coagulation proceeds maintaining constant the density of the grains material (and a spherical shape) ⟨ σ_H ⟩ is proportional to the inverse of the grain size a (see e.g. Flower et al. 2005), and the critical density (∝a²) would increase in the central regions where the grains are larger. The value for n_H2,crit (which is equal to 0.5n_H,crit in a fully molecular environment) derived from this simple estimate is in good agreement with studies of low-mass starless cores (see e.g. Bacmann et al. 2002).

5.3.1. RATRAN modelling

The previous analysis assumes that molecules are in LTE. In order to have more solid results for CO abundances we used RATRAN³ (Hogerheijde & van der Tak 2000) to build one-dimensional models of the clumps for typical parameters of the sample. To reproduce the observed CO line intensities and ratios, we build two grids of models, one set with a central heating source, the other with a constant temperature. For the former we vary the luminosity of the central source L (in 3 steps: 10²L_⊙, 5 × 10³L_⊙, 10⁵L_⊙, which is then translated in a power-law temperature profile; see Rowan-Robinson 1980; Wolfire & Cassinelli 1986) and the clump mass (from ~ 200 M_⊙ to ~ 45 000 M_⊙ in 6 steps), while for the latter we vary the temperature (from 5 K to 15 K in 4 steps) and the mass (as in the other grid). For the models with a central heating source we use grains with thin ice mantles, while for the isothermal clumps we use grains with a thick ice coating, as suggested by their derived f_D. The thick ice mantle on the dust grains increases their opacities, thus reducing the column density of dust (and thus of molecular hydrogen, for a fixed gas-to-dust ratio) needed to obtain a given flux. The clump outer radius, the linewidth and the index p for the power law describing the density profile were fixed at 1 pc, 4 km s^-1 and −1.5, respectively. We also assumed a typical distance of 4 kpc. The models were convolved with a Gaussian as large as the beam size for each wavelength, to simulate our observations.

Figure 12 shows the results for centrally-heated models. The data are shown as black triangles, with marker sizes proportional to the logarithm of the MSX flux in the E band (21 μm). Panels (a)–(c) refer to subsample S1, where we compare C¹⁷O(3 − 2), C¹⁷O(1 − 0) and the submm flux; on the other hand, panels (d)–(f) show the results for subsample S2, using C¹⁷O(3 − 2) and C¹⁸O(2 − 1). In the figure we draw grids with three different depletion factors (coded in three different colours), as shown above the panels. The dotted lines indicate clumps with the same mass, while dashed lines indicate models with the same luminosity of the central object. Figure 13 is the same as Fig. 12, but for isothermal models. Differently from Fig. 12, dashed lines indicate models with the same temperature, rather than luminosity of the central object. Again, three values of f_D are considered, but they are higher in this case (1−4 vs. 4−16).

The figures show that the sources with line ratios ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}{\mathit{/}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}\gtrsim \mathrm{1}$ or ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}{\mathit{/}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{18}}\mathrm{O}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}\gtrsim \mathrm{0.3}$ and large submm fluxes (corresponding also to sources with high MSX fluxes) are much better reproduced by a centrally heated clump, with only moderate depletion. On the other hand, sources with very low line ratios are consistent with cold, isothermal clumps, but need a larger depletion factor to reproduce the line fluxes. The fact that for subsample S2 the points are much more scattered around the model grid may be the result of having assumed the same [¹⁸O]/[¹⁷O] for all sources, while, as observations have different angular resolutions, the geometry of the region and/or a different size of the source than the one assumed for the grid (1 pc) may have a larger impact on the measured ratios.

For sources in D8 and in D24 we tried to reproduce the observed line ratios using temperatures of the order of those suggested by ammonia. With the simple spherical models considered here, it is not possible to do so for the ratios I [¹²C¹⁷O(3 − 2)] /I [¹²C¹⁷O(1 − 0)] ≲ 1 with a constant CO depletion. Even limiting the volume density of the clump to 1 − 3 × 10⁴ cm^-3 does not yield an acceptable combination of line intensities, ratios and peak fluxes in the submm regime. A different approach is to use a drop profile for the CO abundance, i.e. to assume that the abundance is canonical when the density is below a given value, while all CO is locked onto grains for densities above this value. In this case, we constructed a grid of large clumps with a radius ~ 2 pc, to try to simulate the contribution from the external layers. We varied the mass in 10 steps (in the range 550 − 2900 M_⊙ within 1 pc), the critical density defining the size of the central depletion hole (in the range 10⁴ − 10⁵ cm^-3) and the temperature of the depleted layers (in the range 8 − 15 K). The temperature of the external layer was assumed to be constant and equal to 20 K. Figure 14 shows the results of these models, as a yellow-shaded area, compared with the peak fluxes and line intensities of S1. The spread in line intensities is due to the different critical densities above which we assume that CO is completely depleted, while the range in submm flux is mainly caused by the different masses of the clump; the temperature of the central regions has a minor impact for the low values considered. Despite the rather high temperature in the low-density external parts of the clumps, the observed molecular-line ratios can be reproduced with this kind of model.

In conclusion, both using a constant CO abundance throughout the clump and the temperature suggested by the CO line ratio or a drop-profile model we find indications that depletion is important for cold, massive clumps. On the one hand, we have a simple order-of-magnitude estimate of f_D, that can be even 10−20 for cold sources, and on the other hand we can derive the rough size of the depletion hole, to be compared with that derived from the typical lifetimes of massive starless clumps (see Sect. 5.3). Both methods suggest depletion zones with radii ≲0.1 pc, for a typical clump.

5.3.2. RATRAN modelling of individual sources

A more detailed study is carried out for some individual sources in the sample, to properly take into account the source size, radial density distribution, linewidth, expected abundance and distance. The sources were selected from subsample S1 not to be extremely elongated from the 870 μm images. We selected more sources from among groups D8 and D24 objects to confirm and constrain better their large depletion factors. We build a grid of models with RATRAN, with 13 equally spaced steps in mass, and 11 in depletion factor. Depending on whether the clump is centrally heated or isothermal, 15 equal logarithmic steps in luminosity or 15 equal steps in temperature. The radius of the cloud was fixed at the value of R_eff as given in the ATLASGAL catalogue. For sources in the first two groups we use models with a central heating source and thin ice mantles, while for groups D8 and D24 we use models with a constant temperature and thick ice mantles for the dust.

We compare the observed values of line and peak 870 μm continuum fluxes with those predicted by the model, assigning a probability to each set of free parameters of the model (M, L or T, f_D). For the centrally heated models this probability is calculated according to $\begin{array}{ccc}\mathit{P}\mathrm{\left(}\mathit{M,L,}{\mathit{f}}_{\mathit{D}}\mathrm{\right|}\mathit{D,}\mathrm{model}\mathrm{)}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\psi }}\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}\mathit{M,L,}{\mathit{f}}_{\mathit{D}}\mathit{,}\mathrm{model}\mathrm{)}\\ & & \end{array}$(20)where P(D | M,L,f_D,model) is the likelihood, evaluated following $\mathit{P}\mathrm{\left(}\mathit{D}\mathrm{\right|}\mathit{M,L,}{\mathit{f}}_{\mathit{D}}\mathit{,}\mathrm{model}\mathrm{)}\mathrm{\propto }\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}{\mathrm{e}}^{\frac{\mathrm{-}\mathrm{(}{\mathrm{F}}_{\mathrm{i}}\mathrm{-}{\mathrm{F}}_{\mathrm{mod}\mathit{,}\mathrm{i}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{F}\mathit{,}\mathrm{i}}^{\mathrm{2}}}}\mathit{.}$(21)In Eq. (21) F_i are the observed integrated line and 870 μm peak fluxes, F_mod,i are the same fluxes, predicted by the model and convolved with the beam of the observations and σ_F,i are the flux uncertainties. The normalisation constant ψ is chosen so that $\mathrm{\int }\mathit{P}\mathrm{\left(}\mathit{M,L,}{\mathit{f}}_{\mathit{D}}\mathrm{\right|}\mathit{D,}\mathrm{model}\mathrm{)}\mathrm{d}\mathit{M}\mathrm{d}\mathit{L}\mathrm{d}{\mathit{f}}_{\mathit{D}}\mathrm{=}\mathrm{1.}$(22)Finally, P(M,L,f_D | model) is the prior. The priors for M and f_D are constants over the parameter range shown in Fig. C.1. The luminosity derived in this way may not be well constrained, since it is derived only from the ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}$ line ratio in the central regions probed by our single-pointing observations (i.e. representing the average T_ex along the line-of-sight), under the assumption of spherical symmetry. Therefore, using luminosities reported in literature, we set loosely informative Gaussian priors on L for the bright sources. Fazal et al. (2008), van der Tak et al. (2013) and Gaume et al. (1993) give a luminosity (scaled to the distance used here) of ~ 8 × 10³L_⊙, ~ 10⁵L_⊙ and ~ 6 × 10⁴L_⊙ for AGAL19.882-00.534, AGAL034.258+00.154 and AGAL049.489-00.389, respectively.

For isothermal models, the equations for assigning probability are identical to those listed above, with temperature instead of luminosity. The prior for the temperature is taken to be constant.

Figure 15 shows two examples of the probability distributions derived for the modelled sources; the complete figure is available in the electronic version (Fig. C.1).

The probability distribution as a function of mass and luminosity for the centrally heated sources, and mass and temperature for the isothermal objects is computed by integrating the probability cube along the depletion axis. Tables 4 and 5 show the parameters of the sources derived from RATRAN models and the width of their probability distribution, calculated integrating the probability cube along the other two axes. All regions include a large quantity of gas and dust, and the temperatures for groups D8 and D24 are very low.

We note that groups IRB and RMS and groups D8 and D24 indeed show different depletion factors, in the ranges between 1−3 and 3−15, respectively. These depletion factors are averaged along the line-of-sight and in the beam. Because groups D8 and D24 contain sources that are in an earlier stage of evolution than IRB and RMS objects, this result confirms that, as evolution proceeds, the molecules are evaporated from the dust grains in the gas phase, and thus the CO depletion decreases in the environments around high-mass stars. Larger values of f_D are found by Fontani et al. (2012), for a separate sample of massive clumps. In that work the column density of molecular hydrogen is calculated with the expression of Beuther et al. (2005), yielding column densities ~ 2.7 times larger than the expression used in the present paper, because a different dust opacity is used. Hernandez et al. (2011) study in detail depletion in a single IR-dark cloud, showing that a large mass of gas is affected by CO depletion. The large depletion factors found in the present paper, together with those derived by Fontani et al. (2012) show that the freeze-out phenomenon not only occurs in high-mass clumps, but it also involves large masses of gas (at least tens of M_⊙), as suggested by Hernandez et al. (2011).

It is reasonable to expect that the abundance profile is not constant within the clump, and that the depletion factors in the densest and coldest regions are larger than those derived; for example, in Zhang et al. (2009) it is claimed that depletion factors may reach values ≳1000 in the central regions of hight-mass cores. Our estimate may be, in fact, a lower limit for apparently starless objects, an average within an APEX beam along the line-of-sight, since we do not take into account such variable abundance profiles. A simple test performed with RATRAN shows that if we assume a depletion increasing with density, we find that the depletion factors in the inner regions may be much larger than those found with a constant profile, and, because the outer and less dense layers dominate the emission in this situation, one can use a larger gas and dust temperature to reproduce the observed line ratio, due to non-LTE excitation (similarly to the drop-profile models described in Sect. 5.3.1 and below), thus decreasing the mass by a factor of ~ 3 for “cold” sources. However, without a map of molecular emission we do not have enough constraints for these more sophisticated models. Mapping of the clumps in molecular lines and at higher angular-resolution may unveil the actual mass of gas suffering from depletion and constrain the abundance gradient, as well as give a more accurate determination of temperature and mass.

On the other hand, Zinchenko et al. (2009) and Miettinen et al. (2011) find CO undepleted, i.e. with the canonical abundance, towards massive clumps. In both works the sources show signs of active star formation, thus being similar to our more evolved sources. For these clumps with a central heating source the abundance may rise near to the (proto-)star, due to the energy injected in the medium.

We try to model the seven sources of groups D8 and D24 in Table 5 also using models with a drop profile for the abundance of CO. The model used in Sect. 5.3.1 shows that in these cases the temperature of the layers where all CO is frozen out on dust grains is not very important in determining the continuum, therefore we fixed it to 20 K for simplicity. The outer non-depleted layers also have T_d = 20 K (except for one case, as indicated in Table 6). The radius of the clump was fixed to twice the value of R_eff, in order to account for lower-density layers. We varied the mass in 20 equal steps, and the critical density of H₂n_crit, above which CO is completely depleted, in 10 equal logarithmic steps. The results are summarised in Table 6 and shown in Fig. C.2 (available in the electronic version), of which we give an example in Fig. 16.

Sources with low ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}{\mathit{/}}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$ are usually much better reproduced than those with ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{2}\mathrm{)}\mathit{/}$${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}\mathrm{~}\mathrm{1}$, where the models predict an excess of ${}^{\mathrm{12}}\mathrm{C}^{\mathrm{17}}\mathrm{O}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{0}\mathrm{)}$. This more detailed analysis shows that the critical density of molecular hydrogen above which all CO is depleted is very similar to that estimated through Eq. (19) and that the typical depletion radii derived from the source-specific drop profile models are usually larger than those derived from the depletion timescale as a function of radius (see Sect. 5.3); for two sources the discrepancy is even a factor of 3 − 4 (or more, for lifetimes shorter than 10⁵ yr). Given the simplicity of the models and of the estimate obtained using the expression for τ_dep with an age of 10⁵ yr (~ τ_ff; see Sect. 5.3) for the starless clumps, this discrepancy is not surprising. However, it is interesting to note that the larger radii are found for two of the most distant sources (AGAL008.684-00.367 and AGAL010.444-00.017), whereas the nearest source (AGAL034.411+00.234, 1.6 kpc) has R_dep(model) ~ R_dep(τ_dep). This might be the result of a clumpy medium within the beam, with the superposition of smaller regions where depletion is important simulating a larger R_dep. In any case, the results obtained with these models confirm that a large mass of gas appears to be affected by CO depletion, on average ~ 20% of the total mass (from ~ 5% up to ~ 50% in one case). Our estimate of R_dep is comparable to that of the heavily depleted cores observed in Zhang et al. (2009).

Another interesting thing to note is that the mass of gas within R_eff is up to a factor of 5 smaller than in the case of a constant abundance profile and T_d = T_ex (cf. Tables 5 and 6) due to the increased temperature and model clump radius (2R_eff, see above). This difference in mass may be the dominant term in the discrepancy observed in the ratio α_vir = M_vir/M between the sources in IRB and RMS, and D8 and D24 (see Sect. 5.4). In addition, it must be kept in mind that when using CO isotopologues to investigate the properties of massive clumps, the emission may come mainly from the outer layers even though the line is optically thin, whereas the inner, dense layers contribution to the molecular emission is minor because of depletion.

5.4. Stability of the clumps

A simple analysis of the gravitational stability can be performed by comparing the virial mass with the mass calculated from the submm emission by means of Eq. (16), assuming T_d = T_K = T_ex and that the clump is isothermal. MacLaren et al. (1988) give a simple expression to evaluate the virial mass in the case of a spherical clump of radius r, with a given density profile: ${\mathit{M}}_{\mathrm{vir}}\mathrm{\left[}{\mathit{M}}_{\mathrm{\odot }}\mathrm{\right]}\mathrm{=}\mathit{k}\hspace{0.17em}\mathit{r}\mathrm{\left[}\mathrm{pc}\mathrm{\right]}\mathrm{\Delta }{\mathit{V}}^{\mathrm{2}}\mathrm{[}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}\mathrm{]}\mathit{.}$(23)In our case we assumed an average radial density profile n(H₂) ∝ r^-1.5 (see Sect. 5.3.1), implying k = 170, r = R_eff (approximately the FWHM of the clump), and we use the C¹⁷O(3 − 2) line to derive ΔV, so that we use the same line for all subsamples. Figure 17 shows this comparison, with the dashed and dotted lines indicating M_vir = M and M_vir = 2M, respectively (see below): in the upper panel we show M vs. M_vir and in the lower one the virial parameter α_vir = M_vir/M as a function of M. Objects in groups IRB and RMS are mostly distributed around the virial equilibrium line, while the most massive ones in groups D8 and D24 appear to be unstable, i.e. with M_vir ≪ M. The most massive sources in groups D8 and D24 have virial masses ~ 10 times lower than the mass derived from dust emission.

Here we try to understand if these massive sources do really have virial masses much lower than their actual mass, or if it is due to an observational effect. There is the possibility that the temperature we derive is underestimated and that the mass is consequently overestimated. The temperature of apparently starless clumps was previously measured by e.g. Rygl et al. (2010), Giannetti et al. (2013) and Sánchez-Monge et al. (2013) making use of other temperature probes, typically NH₃. The characteristic T_K was found to be ~ 10 − 15 K, slightly larger than what we measure for group D24. However it is similar to that of group D8 and it is difficult to explain such a large difference in mass as being a consequence of this difference in temperature. The ammonia temperatures measured by Wienen et al. (2012) and Wienen et al. (in prep.) suggest temperatures of ~ 20 K for D8 and D24. Even allowing such temperature accounts for a change by only a factor ≲ 3. On the other hand, the temperature we measure in IR-bright clumps may be relatively high and not representative of the whole clump, thus leading one to underestimate M. Again it is difficult to explain a difference of a factor of up to 20 in mass, because T_K from ammonia observations is typically around 20 − 40 K for such sources. Depletion of CO may play a role in this: if the degree of depletion is very high in the centre, we might not be able to probe the physical conditions of the gas there, where the density is very high. In this case we would be tracing only the gas in a more external layer (see also Sect. 5.3.1). If the source is starless, we do not expect the temperature in this external layer to be significantly lower than in the centre; on the contrary, in low-mass sources starless objects seem to have a small negative gradient (e.g. Bergin & Tafalla 2007), being colder in the centre, as fewer photons can penetrate the dense layers. On the other hand, we do not measure a large average depletion for star forming sources (see Sect. 5.3). It is possible to have a larger effect on the mass if the density in the external layers is low, and the molecule is not in LTE. A simple test of the impact of this on the mass of the clump is described in Sect. 5.3.2, and it appears to be a major source of uncertainty. The presence of a thick ice coating on the dust grains changes the dust opacity; according to Ossenkopf & Henning (1994) it may increase by ~ 30% at these frequencies, with respect to a model with thin ice mantles. Using a model with thick mantles for “cold” clumps in groups D8 and D24 would decrease their masses, but not enough to significantly increase the derived M_vir/M ratios. Another source of uncertainty for α_vir is the clump size. The ratio between the major and minor axis of the source is typically between 1 and 2. Therefore, this may cause an uncertainty on the virial mass of a factor of 2 for some clumps. The assumed gas-to-dust ratio of 100 may not be applicable to all sources in the sample, but we do not know enough about its variation to estimate its influence on the derived masses. The uncertainties described above do not seem to be able to account for the largest differences between the mass derived from dust emission and the virial mass.

A clump mass measured from the submm emission larger than the virial mass was found also by Hofner et al. (2000) and Fontani et al. (2002). In the latter work, the authors use CH₃CCH as a temperature probe, which is a symmetric-top molecule. This molecule gives a reliable temperature estimate, thus reducing the uncertainties on the mass. Giannetti et al. (2013) also have temperature determinations both from ammonia and SED fitting, and find that the most massive clumps in their sample have M>M_vir. In a recent work Kauffmann et al. (2013) study in detail the virial parameter in molecular clouds, for a large collection of objects, from entire clouds to cores. These authors find that in high-mass sources the virial mass can be much lower than the mass derived from observations of the dust. The authors consider α_vir ~ 2 as a limit for the gas motions alone to prevent collapse, to take into account a wide range of shapes and density gradients. However, our definition of M_vir already takes into account a density gradient (n ∝ r^-1.5), implying values of α_vir about 20% lower than theirs. They also discuss the observational uncertainties on α_vir and conclude that the very low values of the virial parameter are likely to be real. Kauffmann et al. (2013) fit a straight line to the log (M) − log (α_vir) diagram for high-mass clumps, finding that the slope is ~−0.5, and that α_{vir, min} ~ 0.2 from the fit, for high-mass sources. An analogous fit yields very similar results in both parameters for the clumps in the present paper, excluding the sources in groups IRB and RMS, to have a comparable sample. As shown in Kauffmann et al. (2013), despite the low values of the virial parameter, clumps with M ≫ M_vir are not likely to be undergoing gravitational collapse. López-Sepulcre et al. (2010), on the other hand, find masses consistent with the virial mass. Thus, discordant claims on the virial stability of massive clumps exist in the literature.

If these objects are really unstable, the magnetic field may play a significant role in opposing the collapse. For the most massive clumps, fields of the order of a mG are needed to halt the collapse. In particular, using the expression from Bertoldi & McKee (1992) for cold, magnetised clouds, the average critical magnetic field strength ranges between 0.3 − 3 mG. Magnetic fields of this order of magnitude have been measured in regions of high-mass star formation (e.g. Crutcher 2005; Girart et al. 2009). If we accept that the derived mass may be overestimated by up to a factor of 5, a smaller magnetic field strength is enough to stabilise the clumps, and the magnetic and turbulent energy density are then approximately of the same order of magnitude even for the most unstable sources.

6. Summary and conclusions

We studied several CO isotopologues in 870 μm - bright clumps of the ATLASGAL survey, to investigate the depletion of carbon monoxide, and the [¹²C]/[¹³C] and [¹⁸O]/[¹⁷O] isotopic ratios in regions of (potential) massive star formation. This “TOP100” sample consists of 102 clumps.

The sample is selected to include the brightest submillimetre emission sources in different evolutionary stages, separated in four groups (from IR-bright to 24 μm-dark sources; Motte et al. 2007; Nguyen Luong et al. 2011, see Sect. 2).

From the ratio of different rotational transitions of the observed CO isotopologues we estimate an excitation temperature, with which we derive the molecular column density corrected for optical depth; in order to have a consistent estimate of the optical depth and a refined estimate of the isotopic ratios, we combined all available information from line intensities, ratios and hyperfine structure with a Bayesian approach, allowing the relative isotopic abundances to vary (see Sect. 4.4). Comparing the CO isotopologue column densities with those of H₂, derived from the dust emission at 870 μm, we find that a significant fraction of the sources suffer from depletion, especially the ones with a low T_ex.

The main result of this work is that we find that, just like for low-mass cores, depletion of CO is relevant also in massive sources during their early life, and varies with evolution. Groups D8 (the brightest 8 μm-dark objects) and D24 (the brightest 24 μm-dark sources), typically show larger depletion factors and lower temperatures than the more evolved groups IRB and RMS (the brightest objects of the whole survey and the brightest among the remaining sources classified as MYSOs in the RMS survey, respectively; the clumps in both classes have mid-IR emission). The depletion factors f_D in the less evolved sources may be as large as ~ 20 (see Fig. 8 and Table 5; this corresponds to f_D ~ 55 using the dust opacities adopted in Fontani et al. 2012). These estimates are likely to be lower limits at least for starless sources, since they are averaged along the line-of-sight and derived with a constant abundance profile. On the other hand, the more evolved sources show a typical depletion ~ 3 − 10 times lower (see Sect. 5.3). The larger depletion found in groups D8 and D24 was confirmed with one-dimensional models made with RATRAN, to take into account possible non-LTE effects. Therefore, massive objects seem to follow an evolution of CO depletion similar to that of low-mass objects, where carbon monoxide is frozen onto grains, before the feedback from star formation evaporates the molecules back into the gas phase. The column and volume densities are also found to correlate with the depletion factor: denser sources have a larger f_D, on average. However, among different groups the typical depletion at a given density decreases with increasing temperature, i.e. for more evolved sources. The dependence of the depletion factor on the temperature and density could be explained by the effect of these quantities on the evaporation and freeze-out timescales.

Another way to reproduce the observations is to use models with a central drop in CO abundance. We used models where all CO is locked onto grains for densities above a critical value, whereas below this density, the abundance is canonical. In this case as well the observations can be qualitatively reproduced, again suggesting that CO depletion is important in the dense layers of massive clumps. The estimated critical densities of molecular hydrogen are similar to those found in low-mass starless cores, around a few × 10⁴ cm^-3. Both these models and the comparison of typical lifetimes of starless/IR-quiet clumps with the timescale for depletion as a function of radius suggest that the radius of the central depletion hole is roughly in the range ~ 0.02 − 0.1 pc. However, when trying to model individual sources, the radius of this region is found to lie between ≈0.2 − 0.8 pc, even a factor of ~ 3 − 4 higher with respect to that estimated from the depletion timescale (see Table 6), confirming that a large mass of gas appears to be affected by CO depletion. This discrepancy is not surprising, considering that these are simple order-of-magnitude estimates. It is interesting to note that R_dep is typically 0.2 − 0.3 pc for sources with distances ≲ 6 kpc, and is 0.7 − 0.8 pc for two out of three of the more distant sources. This may simply indicate that, because the medium is clumpy, more than one condensation affected by depletion exists within the beam.

A simple test for the stability of clumps is performed comparing the total mass M derived from dust continuum emission and the virial mass. The clumps are found to be near virial equilibrium for groups IRB and RMS. On the other hand, several sources of groups D8 and D24 seem to be gravitationally unstable, in the sense that M_vir ≪ M. We consider some sources of uncertainty and conclude that this is likely to be real, at least for the sources where M_vir/M is the smallest. The presence of a central region affected by strong CO depletion may be the most important factor of uncertainty in the virial parameter α_vir. This is because the clumps used for such models have masses within R_eff 2−5 times lower than those derived with models with T_ex = T_d and a constant abundance profile.

We also investigated the [¹⁸O]/[¹⁷O] and [¹²C]/[¹³C] isotopic ratios in the inner Galaxy. We find no significant gradient for [¹⁸O]/[¹⁷O] as a function of D_GC for 2 kpc ≲ D_GC ≲ 8 kpc, and the ratios are consistent with ~ 4 with an intrinsic scatter of ~ 1. We find that a few sources with D_GC ~ 4 kpc have [¹⁸O]/[¹⁷O] values of ~ 5.5, similar to the values measured for the pre-solar cloud. [¹²C]/[¹³C] is found to increase with D_GC, as predicted by the models of chemical evolution of the Galaxy (see e.g. Prantzos et al. 1996), with [¹²C]/[¹³C] ~ 66 ± 12 in the solar neighbourhood; however, the intrinsic scatter of the relation is as large as ~ 7 − 13. Milam et al. (2005) show that this large scatter is not likely due to processes such as chemical fraction of selective photodissociation, therefore leaving intrinsic differences (e.g. metallicity, star formation history) between sources, or other processes unaccounted for (such as cloud mergers, non-efficient- or radial gas mixing) to explain it.

Online material

Appendix A: Tables

Appendix B: Spectra

Appendix C: RATRAN figures

---

Acknowledgments

This work was partly funded by the Marco Polo program (Università di Bologna), making it possible for A.G. to spend three months at the Max Planck Institute für Radioastronomie in Bonn. A.G. thanks the Max Planck Institute für Radioastronomie and F.W. for their hospitality and support. T.Cs.’s contribution was funded by the ERC Advanced investigator grant GLOSTAR (247078). This research made use of data products from the Midcourse Space Experiment. Processing of the data was funded by the Ballistic Missile Defense Organization with additional support from NASA Office of Space Science. This research has also made use of the NASA ADS, SIMBAD, CDS (Strasbourg) and NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

All Tables

All Figures