FAR-INFRARED DUST TEMPERATURES AND COLUMN DENSITIES OF THE MALT90 MOLECULAR CLUMP SAMPLE

We use public HSO data from Hi-GAL (Molinari et al. 2010) observed between 2010 January and 2012 November and obtained from the Herschel Science Archive. The observations were made using the parallel, fast-scanning mode, in which five wavebands were observed simultaneously using the PACS (Poglitsch et al. 2010) and the SPIRE (Griffin et al. 2010) bolometer arrays. The data version obtained from the Herschel Science Archive corresponds to the Standard Product Generation version 9.2.0.

Columns 1–4 of Table 1 list the instrument, the representative wavelength in microns of each observed band, the angular resolution represented by the FWHM of the point-spread function (Olmi et al. 2013), and the estimated point source sensitivity (${\sigma }_{{\rm{p}}}$), respectively. The point source sensitivity, assuming Gaussian beams, is given by ${\sigma }_{{\rm{rms}}}{{\rm{\Omega }}}_{{\rm{b}}}{({{\rm{\Omega }}}_{{\rm{b}}}/2{{\rm{\Omega }}}_{{\rm{pix}}})}^{-1/2},$ where ${\sigma }_{{\rm{rms}}}$ is the rms variations in intensity units, ${{\rm{\Omega }}}_{{\rm{b}}}$ is the beam solid angle, and ${{\rm{\Omega }}}_{{\rm{pix}}}$ is the pixel solid angle.⁸ The fifth column gives the noise level of the convolved and re-gridded maps (see Sections 3.1 and 3.2) and the sixth column lists the observatory where the data were taken. Throughout this work, we will refer to the data related to a specific waveband by their representative wavelength in micrometers. The position uncertainty of the Hi-GAL maps is ∼3''.

Table 1.  Observational Characteristics of the Data

Instrument	Band	FWHMa	${\sigma }_{{\rm{p}}}$	σ	Telescope
	(μm)	('')	(mJy)	(MJy sr−1)
PACS	70	9.2	18	1.1	HSO
PACS	160	12.0	24	1.5	HSO
SPIRE	250	17.0	12	0.7	HSO
SPIRE	350	24.0	12	0.7	HSO
SPIRE	500	35.0	12	0.7	HSO
LABOCA	870	19.2	25	1.5	APEX

Notes.The noise displayed for the SPIRE data is taken in nominal mode. Noise associated with data taken in bright mode is ∼4 times higher.

^aThe beam FWHM associated with each waveband of the PACS and SPIRE data was adopted from Olmi et al. (2013).

Download table as:  ASCIITypeset image

The generation of maps that combine the two orthogonal scan directions was done using the Herschel Interactive Processing Environment (HIPE) versions 9.2 and 10. Cross-scan combination and destriping were performed over 42 Hi-GAL fields of approximately $2\buildrel{\circ}\over{.} 2\times 2\buildrel{\circ}\over{.} 2$ using the standard tools available in HIPE. Columns 1–4 of Table 2 give the target name, the ID of the observation, the observing mode, and the observation dates, respectively. For the SPIRE maps, we applied the extended source calibration procedure (Section 5.2.5 from the SPIRE Handbook⁹ ) since most of MALT90 sources correspond to dense clumps that are comparable to or larger than the largest SPIRE beam size. The saturation limit of the nominal mode for SPIRE (Section 4.1.1 from the SPIRE Handbook) is approximately 200 Jy beam⁻¹. To prevent saturation, fields with longitudes $| l| \leqslant 5$° were observed with SPIRE using the bright observing mode instead of the nominal observing mode.

In addition to Hi-GAL data, we used data from three observations made using the SPIRE bright mode by the HOBYS key project (Herschel Imaging Survey of OB YSOs, Motte et al. 2010). Table 2 lists these observations' IDs. They were directed toward the NGC 6334 ridge and the central part of M17, areas which are heavily saturated in the Hi-GAL data.

Data at 870 μm were taken between 2007 and 2010 using the bolometer LABOCA (Siringo et al. 2009) installed on the APEX telescope located in Chajnantor valley, Chile, as part of the ATLASGAL key project (Schuller et al. 2009). Calibrated and reduced fits images were obtained from the data public releases made by Contreras et al. (2013) and Urquhart et al. (2014a). Table 1 displays the angular resolution, the point source sensitivity calculated as in Section 2.1 using ${\sigma }_{{\rm{rms}}}=60$ mJy beam⁻¹ and ${{\rm{\Omega }}}_{{\rm{b}}}/{{\rm{\Omega }}}_{{\rm{pix}}}=11.6$ (Contreras et al. 2013), and the typical noise of the convolved and re-gridded ATLASGAL maps. In addition to this noise, we assume a 10% uncertainty in the absolute calibration.

To a first approximation, the intensity assigned to each pixel is given by the average of the bolometer readings that covers that pixel position. The spatial sampling of the maps, on the other hand, includes ∼3 pixels per beamwidth. Observed astronomical signals vary spatially on angular scales ≳1 beamsize. Therefore, in the large fraction of the map area that is away from very strong sources, we expect that the differences between adjacent pixels are dominated by instrumental noise. In order to estimate this noise, we use the high-pass filter defined by Rank et al. (1999) to determine the distribution of pixel-to-pixel variations and filter out astronomical emission. The width of this distribution determines the typical noise through the relation $2.36\sigma =\mathrm{FWHM}.$ The advantage of this method is that it gives us an extra and relatively simple way to estimate the noise of the final maps. The noise estimation is similar to that obtained from jackknife maps, produced by taking the difference between maps generated by the two halves of the bolometer array (see Nguyen et al. 2010, for an analogous procedure).

The 1σ point source sensitivities derived from the high-pass filter method described above are typically 18 and 24 mJy for the two PACS bands at 70 and 160 μm, and 12 mJy for the three SPIRE bands at 250, 350, and 500 μm. These derived sensitivities are in good agreement with the ones expected for the Hi-GAL survey (Molinari et al. 2010) and in reasonable agreement with the sensitivities expected for the parallel mode,¹⁰ with the possible exception of the 160 μm band where we estimate about half of the expected noise. The noise value derived at 250 μm is comparable with the noise component derived by Martin et al. (2010) also from Hi-GAL data, indicating that our estimation effectively filters most of the sky emission variations, including the cirrus noise. Finally, and as expected, we find that the noise in fields observed in the SPIRE bright mode is ∼4 times larger compared to that in fields observed in nominal mode. For subsequent analyses, we consider an additional independent calibration uncertainty of 10% whenever we compare data among different bands, as for example, in the SED fitting. This 10% represents a conservative approximation of the combined calibration uncertainty of the SPIRE photometers (5.5%) and the beam solid angle (4%, see Section 5.2.13 of the SPIRE Handbook).

Multi-wavelength studies of extended astronomical objects, such as star-forming clumps and IRDCs, often combine data taken with different angular resolutions. Therefore, to make an adequate comparison of the observed intensities, it is necessary to transform the images to a common angular resolution. We accomplish this by convolving the images to the lowest available resolution, given by the 500 μm SPIRE instrument, using the convolution kernels of Aniano et al. (2011) in the case of HSO data. The ATLASGAL data were convolved by a two-dimensional (2D) Gaussian with FWHM equal to $\sqrt{35\buildrel{\prime\prime}\over{.} {0}^{2}-19\buildrel{\prime\prime}\over{.} {2}^{2}}\approx 29\buildrel{\prime\prime}\over{.} 3,$ under the assumption that the point-spread functions of the ATLASGAL and the 500 μm data are Gaussians. In addition, to compare the intensity of the HSO images with that of the APEX telescope, we need to remove from the HSO data the low spatial frequency emission that has been filtered from the ATLASGAL images. The ATLASGAL spatial filtering is performed during the data reduction, and is a by-product of the atmospheric subtraction method which removes correlated signal between the bolometers (Siringo et al. 2009). As a consequence, any uniform astronomical signal covering spatial scales larger than 25 is lost (Schuller et al. 2009).

We filter the HSO data in a similar way by subtracting a background image from each field and at each band. We assume that this background is a smooth additive component that arises from diffuse emission either behind or in front of the clump. In addition to filtering the HSO data in order to combine it with ATLASGAL, the background subtraction serves two more purposes: it separates the Galactic cirrus emission from the molecular clouds (e.g., Battersby et al. 2011), and it corrects for the unknown zero level of the HSO photometric scale. Our background model consists of a lower-envelope of the original data under two constrains: its value at each pixel has to be less than in the image, within a 2σ tolerance, and it has to vary by less than 10% over 25, which corresponds to the ATLASGAL filter angular scale.

We construct a background image for each Hi-GAL field following a slight modification of the CUPID-findback algorithm¹¹ of the Starlink suite (Berry et al. 2013). The iterative algorithm used to construct the background starts with the original image. Then, we calculate a smoothed image by setting to zero (in the Fourier transform plane) the spatial frequencies corresponding to flux variations on angular scales $\lt 2\buildrel{\,\prime}\over{.} 5.$ For each pixel in this smoothed image with a value larger than the corresponding pixel in the original image plus $2\sigma ,$ the pixel value from the smoothed image is replaced by the one in the original image, where σ is the uncertainty of the map. The remaining pixels in the smoothed image are kept unchanged. The resultant map is the first iteration of the algorithm. This first iteration replaces the starting image and the cycle repeats, generating further iterations, until the change between two consecutive iterations is less than 5% in all pixels.

Figure 1 shows an example of this process, which converges to a smooth lower-envelope of the original image. The solid black line shows a cut along $l=355\buildrel{\circ}\over{.} 8$ of the intensity measured at 250 μm. Negative intensity values away from the Galactic plane are a consequence of the arbitrary zero-level of the HSO photometry scale. Dashed lines show different iterations of the algorithm and the final adopted background is marked in red. The error bar at the center of the plot measures 25, that is, the shortest angular scale filtered by the background. Note that Figure 1 shows a cut across latitude at a fixed longitude, but the algorithm works on the 2D image, not assuming any particular preferred direction.

Zoom In Zoom Out Reset image size

We interpret the observed intensities as arising from a single-temperature gray-body dust emission model. The monochromatic intensity at a frequency ν is given by

$$\begin{eqnarray}&&{I}_{\nu }({T}_{{\rm{d}}},{N}_{{\rm{g}}})={B}_{\nu }({T}_{{\rm{d}}})(1-{e}^{-{\tau }_{\nu }}),\end{eqnarray} \tag{ 1 }$$

where ${B}_{\nu }({T}_{{\rm{d}}})$ is the Planck function at a dust temperature T_d and

$$\begin{eqnarray}&&{\tau }_{\nu }={N}_{{\rm{dust}}}{\kappa }_{\nu },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=\;{\rm{GDR}}\times {N}_{{\rm{g}}}{\kappa }_{\nu },\end{eqnarray} \tag{ 3 }$$

where ${\tau }_{\nu }$ is the dust optical depth, ${N}_{{\rm{dust}}}$ is the dust column density, and ${\kappa }_{\nu }$ is the dust absorption coefficient. The relation between the dust and gas (N_g) column densities is determined by the gas-to-dust mass ratio (GDR), which we assume is equal to 100. We also define the particle column density by ${N}_{{\rm{p}}}:= {N}_{{\rm{g}}}/(\mu {m}_{{\rm{H}}}),$ where $\mu =2.3.$ The number column density of molecular hydrogen (${N}_{{{\rm{H}}}_{2}}$) is obtained in the same way but using $\mu =2.8$ (Kauffmann et al. 2008), under the assumption that all the hydrogen is in molecular form. We assume throughout this work that N_g is measured in g cm⁻² and ${N}_{{{\rm{H}}}_{2}}$ and N_p in cm⁻². To compare the dust emission model to the data, we weight the intensity given by Equation (1) by the spectral response function of the specific waveband, in order to avoid post-fitting color corrections (see for example, Smith et al. 2012).

We exclude the 70 μm intensity from the single T_d fitting since this emission cannot be adequately reproduced by Equation (1) (see Section 3.4). This problem has been noted by several authors (e.g., Elia et al. 2010; Battersby et al. 2011; Smith et al. 2012; Russeil et al. 2013), who have provided at least three possible reasons:

• 1.  
  emission at this wavelength comes from a warmer component;
• 2.  
  cold and dense IRDCs are seen in absorption against the Galactic plane at 70 μm rather than emission;
• 3.  
  a large fraction of the 70 μm emission comes from very small grains, where the assumption of a single equilibrium temperature is not valid.

For each pixel and given the observed background-subtracted intensities ${I}_{\nu ,\mathrm{obs}},$ we minimize the squared difference function,

$$\begin{eqnarray}&&{\chi }^{2}({T}_{{\rm{d}}},{N}_{{\rm{g}}})=\displaystyle \sum _{\nu }\displaystyle \frac{{({I}_{\nu ,\mathrm{obs}}-{\tilde{I}}_{\nu })}^{2}}{{\sigma }_{\nu }^{2}},\end{eqnarray} \tag{ 4 }$$

where the sum is taken over the observed frequencies (i.e., five bands) and ${\tilde{I}}_{\nu }$ is the intensity spectrum predicted by the model weighted by the respective bandpass. The best-fit dust temperature, T_d, and gas column density, N_g, minimize the ${\chi }^{2}$ value. The variance ${\sigma }_{\nu }^{2}$ is equal to the sum in quadrature of the noise (taken from Table 1) plus 10% of the background-subtracted intensity. We fit the model described in Equation (1) for all the pixels with intensities larger than $2{\sigma }_{\nu }$ in all bands.

The reduced ${\chi }^{2},$ defined as ${\chi }_{{\rm{r}}}^{2}:= {\chi }_{{\rm{min}}}^{2}/(m-p)$ (Bevington & Robinson 2003), is a simple measure of the quality of the model. Here, ${\chi }_{{\rm{min}}}^{2}$ is the minimized ${\chi }^{2}$ of Equation (4), m is the number of data-points, and p is the number of fitted parameters. In our case, we fit the dust temperature and the logarithm of the gas column density, so p = 2. Under the hypothesis that the data are affected by ideal, normally distributed noise, ${\chi }_{{\rm{r}}}^{2}$ has a mean value of 1 and a variance of $2/(m-p).$

Figure 2 shows the ${\chi }_{{\rm{r}}}^{2}$ cumulative distribution function (CDF), calculated using all the pixels for which we fit the SED. The median ${\chi }_{{\rm{r}}}^{2}$ value is 1.6. This value is less than $1+\sqrt{2/3}\approx 1.8,$ which is the expected value plus 1σ under the assumption of normal errors for any particular fit. We conclude that the SED model is in most cases adequate, or equivalently, the limited amount of photometric data does not justify a more complicated model. Note that, although the distribution of ${\chi }_{{\rm{r}}}^{2}$ has a reasonable mean and median, it has a large tail: the 95% quantile is located at ${\chi }_{{\rm{r}}}^{2}\approx 9.6.$ This value represents a poor fit to the model, which can be usually attributed to a single discordant data-point. Generally, this point corresponds to the 870 μm intensity, which illustrates the difficulties of trying to match the spatial filtering of the HSO with ATLASGAL data, despite the background correction and common resolution convolution. We re-examine the fitting when the ${\chi }_{{\rm{r}}}^{2}$ value is larger than 10 and remove from the fitting at most one data point only if its removal decreases the ${\chi }_{{\rm{r}}}^{2}$ value by a factor of 10 or more.

Zoom In Zoom Out Reset image size

3.3.1. Spectral Index of the Dust Absorption Coefficient

At frequencies $\nu \lt 1$ THz, the dust absorption coefficient curve ${\kappa }_{\nu }$ is well approximated by a power-law dependence on frequency with spectral index β (Hildebrand 1983), that is,

$$\begin{eqnarray}&&{\kappa }_{\nu }={\kappa }_{0}{(\nu /{\nu }_{0})}^{\beta }.\end{eqnarray} \tag{ 5 }$$

In principle, it is possible to quantify β toward regions where the emission is optically thin and the temperature is high enough such that the Rayleigh–Jeans (R–J) approximation is valid. In this case, from Equations (1) and (5) we deduce that

$$\begin{eqnarray}&&{I}_{{\nu }_{1}}/{I}_{{\nu }_{2}}={({\nu }_{1}/{\nu }_{2})}^{\beta +2},\end{eqnarray} \tag{ 6 }$$

which is independent of temperature.

We estimate β through Equation (6) using low frequency ($\lt 600$ GHz) data taken toward warm ($\gt 30$ K) sources to ensure that the R–J and the dust absorption coefficient power-law approximations are valid. Using this value of β we will be able to justify better the selection of a dust opacity law among the different theoretical models (e.g., Ormel et al. 2011). In order to ensure that the sources used to estimate β are warm enough for Equation (6) to be valid, we select IRAS sources that are part of the 1.1 mm Bolocam Galactic Plane Survey (Rosolowsky et al. 2010; Ginsburg et al. 2013) and the ATLASGAL catalog at 870 μm. We also require that they fulfill ${S}_{60}/{S}_{100}\gt 0.5,$ where S₆₀ and S₁₀₀ are their fluxes at 60 and 100 μm, respectively. In addition, we select sources with $| l| \gt 10^\circ$ in order to avoid possible confusion that may arise in the crowded regions around the Galactic center. We find 14 IRAS sources fulfilling these requirements: 18079−1756, 18089−1837, 18114−1825, 18132−1638, 18145−1557, 18159−1550, 18162−1612, 18196−1331, 18197−1351, 18223−1243, 18228−1312, 18236−1241, 18247−1147, and 18248−1158.

The average β calculated for these sources using Equation (6) is 1.6, but with a dispersion of 0.5 among the sources. This dispersion is large, but it is compatible with a 15% uncertainty in the fluxes. The spectral index is in agreement with the absorption coefficient law of silicate-graphite grains, with $3\times {10}^{4}\;{\rm{years}}$ of coagulation, and without ice coatings according to the dust models from Ormel et al. (2011). For the rest of this work, we use this model of dust for the SED fitting. The tables compiled by Ormel et al. (2011) also sample the frequency range of interest for this work in more detail than the frequently used dust models of Ossenkopf & Henning (1994). We refrain from fitting β together with the SED for two reasons: (i) we lack the adequate data to effectively break the degeneracy between β and T_d, that is, good spectral sampling of highly sensitive data below 500 GHz; and (ii) the range of dust models explored by fitting β includes only power laws instead of using more physically motivated tabulated dust models.

To compare our results with previous studies, which may have derived temperatures and column densities using different hypotheses, we review how different assumptions on β affect the best-fit estimation of T_d. Several studies (e.g., Shetty et al. 2009a, 2009b; Juvela & Ysard 2012) have discussed this problem in association with least-squares SED fitting in the presence of noise. They find that β and T_d are somewhat degenerate and associated with elongated (sometimes described as banana-shaped) best-fit uncertainty regions in the β–T_d plane. In this work, we stress one aspect that has not been sufficiently emphasized: there are two behaviors of the β–T_d degeneracy: one is evident when the data cover the SED peak, and the other when the data only cover the R–J part of the spectrum. In the first case the degeneracy is well described by the modified Wien displacement law

$$\begin{eqnarray}&&\displaystyle \frac{h{\nu }_{{\rm{peak}}}}{{{kT}}_{{\rm{d}}}}\approx (\beta +3),\end{eqnarray} \tag{ 7 }$$

that is, the uncertainty region of T_d and β is elongated along the curve defined by Equation (7). In Equation (7), ${\nu }_{{\rm{peak}}}$ represents the frequency where the SED takes its maximum value, which under optically thin conditions is proportional to the temperature. The proportionality constant depends on β in a complicated way, but the approximation of Equation (7) is correct within a 10% for $\beta \gt 1$ and within a 20% for all $\beta \geqslant 0.$ Note that by assuming a value of β and determining ${\nu }_{{\rm{peak}}}$ observationally we can estimate T_d using Equation (7) in a simple way. Sadavoy et al. (2013) and Shetty et al. (2009a, their 20 K case) show examples of uncertainty regions given by the iso-contours of the ${\chi }^{2}$ function which are elongated along the curve defined in Equation (7). On the other hand, if the spectral range of the data does not cover the observed peak of the SED and covers only the R–J region, the degeneracy between β and T_d is better described by the following relation,

$$\begin{eqnarray}&&\beta -\displaystyle \frac{h{\nu }_{m}}{2{{kT}}_{{\rm{d}}}}={\rm{constant}},\end{eqnarray} \tag{ 8 }$$

where ${\nu }_{m}$ is the highest observed frequency. This relation describes well the degeneracy of the high temperature curves (60 and 100 K) shown in Shetty et al. (2009a). The constant in the right-hand side of Equation (8) is approximately

$$\begin{eqnarray*}&&2+\displaystyle \frac{d\mathrm{ln}{S}_{{\nu }_{m}}}{d\mathrm{ln}\nu },\end{eqnarray*}$$

that is, 2 plus the logarithmic derivative (or spectral index) of the spectrum evaluated in the highest observed frequency. In practice, the exact value of the constants in the right-hand side of Equations (7) and (8) can be determined from the best-fit solutions. In this work, the HSO bands usually cover the peak of the SED so Equation (7) is more pertinent. Depending on the spectral sampling, we can use Equation (7) or (8) to compare temperatures between studies that assume different values of β. For example, the emission in the HSO bands from a cloud of ${T}_{{\rm{d}}}=15$ K with a $\beta =1.0$ dust absorption law is also consistent, by Equation (7) and assuming 10% uncertainty, with the emission coming from a cloud of ${T}_{{\rm{d}}}=12$ K and $\beta =2.$ In each case, the HSO bands cover the peak of the SED. We use this method to re-scale and compare best-fit temperatures obtained from the literature in Section 4.

We estimate the best-fit parameter uncertainties using the projection of the 1σ contour of the function ${\rm{\Delta }}{\chi }^{2}:= {\chi }^{2}-{\chi }_{{\rm{min}}}^{2}$ (Lampton et al. 1976). In the case of 2 fitted parameters, the 1σ uncertainty region is enclosed by the ${\rm{\Delta }}{\chi }^{2}=2.3$ contour. The parameter uncertainties for the SED fitting are given by the projections of these uncertainty regions onto the T_d and $\mathrm{log}{N}_{{\rm{g}}}$ axes. For pixels in images observed using the nominal observing mode, the projections are well described by the following equations

$$\begin{eqnarray}\delta {T}^{-} & = & {\eta }_{10}(0.3-0.4{T}_{10}+0.4{T}_{10}^{2}),\\ \delta {T}^{+} & = & {\eta }_{10}(1.1-1.3{T}_{10}+0.7{T}_{10}^{2}),\\ \delta \mathrm{log}{N}_{{\rm{g}}} & = & {\eta }_{10}(0.03-0.03\mathrm{log}{N}_{{\rm{g}}}),\end{eqnarray} \tag{ 9 }$$

where ${T}_{10}={T}_{{\rm{d}}}/(10\;{\rm{K}}),\;$ N_g is in g cm⁻², and ${\eta }_{10}$ is the flux calibration uncertainty in units of 10%. The best-fit temperature and log-column density with their 1σ uncertainties are given by ${{T}_{{\rm{d}}}}_{-\delta {T}^{-}}^{+\delta {T}^{+}}$ and $\mathrm{log}{N}_{{\rm{g}}}\pm \delta \mathrm{log}{N}_{{\rm{g}}},$ respectively. For pixels in images observed using the bright observing mode, the projections are well described by

$$\begin{eqnarray}\delta {T}^{-} & = & {\eta }_{10}(0.7-0.71{T}_{10}+0.53{T}_{10}^{2}),\\ \delta {T}^{+} & = & {\eta }_{10}(1.1-1.3{T}_{10}+0.74{T}_{10}^{2}),\\ \delta \mathrm{log}{N}_{{\rm{g}}} & = & {\eta }_{10}(0.05-0.03\mathrm{log}{N}_{{\rm{g}}}).\end{eqnarray} \tag{ 10 }$$

Equations (9) and (10) were derived by fitting the upper and lower limits of T_d and $\mathrm{log}{N}_{{\rm{g}}}$ projections of the uncertainty region. These approximations for the uncertainty are valid for T_d between 7 and 40 K, for $\mathrm{log}{N}_{{\rm{g}}}$ between −3.4 and 1.1 (equivalent to $\mathrm{log}{N}_{{{\rm{H}}}_{2}}$ between 19.9 and 24.4), and for values of ${\eta }_{10}$ between 1 and 2, which correspond to 10% and 20% calibration errors, respectively. Figure 3 shows an example of the prediction of Equations (9) compared to the ${\rm{\Delta }}{\chi }^{2}$ contours. Within their range of validity and for data taken in the nominal mode, the intervals $[{T}_{{\rm{d}}}-\delta {T}^{-},{T}_{{\rm{d}}}+\delta {T}^{+}]$ and $[\mathrm{log}{N}_{{\rm{g}}}-\delta \mathrm{log}{N}_{{\rm{g}}},\mathrm{log}{N}_{{\rm{g}}}+\delta \mathrm{log}{N}_{{\rm{g}}}]$ correspond to the projections of the uncertainty ellipse onto the T_d and $\mathrm{log}{N}_{{\rm{g}}}$ axes within 0.2 K and 0.02 dex, respectively.

Zoom In Zoom Out Reset image size

Equations (9) also indicate that, while the confidence interval (CI) of $\mathrm{log}{N}_{{\rm{p}}}$ is symmetric and roughly constant, the temperature uncertainties grow rapidly above 25 K and they are skewed toward the higher value (Hoq et al. 2013). This is produced mainly because of the absence of data at wavelengths shorter than 160 μm. As explained in Section 3.3, we do not use the 70 μm data in our model. Including the 70 μm data increases the median of the ${\chi }_{{\rm{r}}}^{2}$ distribution to ∼3.5. The temperature and log-column density uncertainties are also correlated along the approximate direction where the product ${T}_{{\rm{d}}}$ × ${N}_{{\rm{g}}}$ is constant, indicated in Figure 3. The better the R–J approximation is for the SED, the better will be the alignment of the major axis of the ellipse with the line ${T}_{{\rm{d}}}\times {N}_{{\rm{g}}}=\mathrm{constant}.$

The HSO detectors in SPIRE and PACS have saturation limits that depend on the observing mode. The saturation intensities for the nominal observing mode are 220 and 1125 Jy beam⁻¹ for the 70 and 160 μm PACS bands¹² , respectively, and 200 Jy beam⁻¹ for SPIRE. Saturation is most problematic in the 250 μm SPIRE band.

There are 46 MALT90 sources whose Hi-GAL data are affected by saturation. Of these, six are covered by HOBYS observations made by using the bright mode (Section 2.2) that gives reliable 250 μm intensities. For the remaining 40 sources, we replace the saturated pixels with the saturation limits given above, and we fit the SED taking these values as lower bounds.

MALT90 and other previous studies (e.g., Molinari et al. 2008; López-Sepulcre et al. 2010; Sanhueza et al. 2012; Giannetti et al. 2013; Sánchez-Monge et al. 2013; Csengeri et al. 2014) use mid-IR observations as a probe of star formation activity. However, deeply embedded, early star formation activity could be undetected at mid-IR yet be conspicuous at far-IR. Quantitatively, we expect that the 24–70 μm flux density ratio of HMYSOs to vary between ${10}^{-6}$ and 1 for a wide range of molecular core masses (60–240 M_⊙) and central star masses over 1 M_⊙ (Zhang et al. 2014). Therefore, despite MIPSGAL having ∼50 times better point source sensitivity at 24 μm than that of Hi-GAL at 70 μm, it is possible to detect embedded protostars in 70 μm images that would otherwise appear dark at 24 μm and would be classified as quiescent. The 70 μm data thus allow us to further refine the MALT90 classification since a truly quiescent clump should lack 70 μm compact sources (e.g., Beuther et al. 2013; Sanhueza et al. 2013; Tan et al. 2013).

We examined the Hi-GAL 70 μm images of the 616 quiescent sources and found 91 (15%) that show compact emission at 70 μm within 38''—or one Mopra telescope beamsize—of the nominal MALT90 source position. Hereafter, we consider these sources as part of the protostellar sub-sample. We also found 83 sources that appear in absorption at 70 μm against the diffuse Galactic emission. We refer to these clumps as far-IR dark clumps (far-IRDCs). The remaining 442 quiescent sources are either associated with diffuse emission not useful for tracing embedded star formation, or they are confused with the 70 μm diffuse emission from the Galactic plane.

The dust temperature is fixed by the balance between heating and radiative cooling of the grain population. If the density in a molecular cloud is greater than $5\times {10}^{4}$ cm⁻³ (Goldsmith 2001; Galli et al. 2002), we expect the dust temperature to be coupled to the gas temperature. We test this hypothesis by comparing the average dust temperatures with the gas kinetic temperatures determined from ammonia observations. Figure 6 shows the ammonia temperature derived by Dunham et al. (2011, 23 matching sources), Urquhart et al. (2011, 10 matching sources), Wienen et al. (2012, 106 matching sources), and Giannetti et al. (2013, 19 sources) versus the dust temperature, separated by evolutionary stage. Dunham et al. (2011) and Urquhart et al. (2011) performed the NH₃ observations using the Green Bank Telescope at 33'' angular resolution. Wienen et al. (2012) used the Effelsberg Radiotelescope at 40'' angular resolution. These are comparable to the resolution of our dust temperature maps (35''). On the other hand, Giannetti et al. (2013) used NH₃ data obtained from ATCA with an angular resolution of $\sim 20^{\prime\prime} .$ In this last case, we compare their ammonia temperatures with ${T}_{{\rm{d}},{\rm{P}}}$ instead of $\bar{{T}_{{\rm{d}}}}$.

Zoom In Zoom Out Reset image size

All but eight MALT90 sources with ammonia temperature estimation are classified in one of the four evolutionary stages. The Spearman correlation coefficient of the entire sample (158 sources, including these eight with Uncertain mid-IR classification) is 0.7, with a 95% CI between 0.6 and 0.8, indicating a positive correlation between both temperature estimators. The scatter of the relation is larger than the typical temperature uncertainty, and it grows with the temperature of the source. For sources below 22 K, ammonia and dust temperatures agree within ±3 K. Above 22 K, the uncertainties of both temperature estimators become larger (see Equations (9) and, for example, Walmsley & Ungerechts 1983), consistent with the observed increase of the scattering. In addition, higher temperature clumps are likely being heated from inside and therefore associated with more variations in the dust temperatures along the line of sight, making the single-temperature approximation less reliable. The slopes of the linear regressions performed in the data are 0.7 ± 0.1, 0.8 ± 0.1, 0.7 ± 0.1, and 0.9 ± 0.3 for the quiescent, protostellar, H ii region, and PDR samples, respectively. The ammonia and dust temperatures relation agrees in general with that found by Battersby et al. (2014), except that we do not find a systematically worse agreement on quiescent sources compared with the other evolutionary stages.

Figure 7 shows maps of smoothed 2D histograms of the distributions of $\bar{{T}_{{\rm{d}}}}$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ of the MALT90 clumps for each mid-IR classification. In the following analysis we focus on these two quantities and their relation with the evolutionary stage. We use ${N}_{{\rm{g}},{\rm{P}}}$ instead of $\bar{{N}_{{\rm{g}}}}$ because ${N}_{{\rm{g}},{\rm{P}}}$ is independent of the specific criterion used to define the extension of the clump and because the column density profiles are often steep ($\propto {s}^{-0.8},$ where s is the plane-of-the-sky distance to the clump center, Garay et al. 2007), making the average $\bar{{N}_{{\rm{g}}}}$ less representative of the clump column density values. On the other hand, dust temperature gradients are shallower ($\propto {r}^{-0.4},$ where r is the distance to the clump center, see van der Tak et al. 2000) and $\bar{{T}_{{\rm{d}}}}$ has less uncertainty compared with the temperature calculated toward a single point. We include in the protostellar group those sources that are associated with a compact source at 70 μm (Section 3.6). The most conspicuous differences between the evolutionary stages are evident between the quiescent/protostellar and the H ii region/PDR populations (Hoq et al. 2013). The main difference between these groups is the temperature distribution. Most of the sources in the quiescent/protostellar stage have temperatures below 19 K, while most H ii region/PDR sources have temperatures above 19 K. We also note that the quiescent, protostellar, and H ii region populations have peak column densities $\gtrsim 0.1$ g cm⁻², equivalent to $2.13\times {10}^{22}$ H₂ molecules cm⁻², while the PDR population has peak column densities of typically half of this value.

Zoom In Zoom Out Reset image size

These differences are also apparent in Figure 8, where solid lines display the marginalized CDFs of $\bar{{T}_{{\rm{d}}}}$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ for each evolutionary stage. The dashed lines show the distributions of the Uncertain group. It is clear from these plots that the median temperature increases monotonically with evolutionary stage, and that the protostellar and PDR clumps are the stages associated with the largest and smallest column densities, respectively. Figure 9 shows Tukey box plots (Feigelson & Babu 2012, Section 5.9) of the marginalized distributions of $\bar{{T}_{{\rm{d}}}}$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ separated by evolutionary stage. In these plots, the boxes indicate the interquartile range (half of the population), the thick horizontal line inside each box indicates the median, and the error bars encompass the data that are within 1.5 times the inter-quartile distance from the box limits. The remaining points, in all cases less than 4% of the sample, are plotted individually with small circles and we refer to them formally as outliers. Figure 9 shows that the $\bar{{T}_{{\rm{d}}}}$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ interquartile range shifts with evolutionary stage. This is evidence of systematic differences between the different populations, despite the large overlaps. In practice, the overlap between populations implies that it is unfeasible to construct sensible predictive criteria that could determine the evolutionary stage of a specific source based on its temperature and peak column density and it also reflects that star formation is a continuous process that cannot be precisely separated into distinct stages. Nevertheless, the fact that the proposed evolutionary stages show a monotonic increase in mean temperature demonstrates that the classification scheme has a legitimate physical basis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the following, we focus our analysis on the quiescent and protostellar populations. Figures 7–9 show that these two samples are similarly distributed and they exhibit the most important overlap. We test the statistical significance of the quiescent and protostellar differences in $\bar{{T}_{{\rm{d}}}}$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ by comparing these differences with their uncertainties. Table 4 shows the medians, means, and rms deviations of $\bar{{T}_{{\rm{d}}}},\;$ ${T}_{{\rm{d}},{\rm{P}}},\;$ $\mathrm{log}\bar{{N}_{{\rm{g}}}},$ and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ for each population. In general, these dispersions are larger than the uncertainties of the individual values, indicating that the dispersions are intrinsic to each population and not due to the fitting uncertainties. The means $\bar{{T}_{{\rm{d}}}}$ for the quiescent and protostellar populations are 16.8 and 18.6 K, respectively. The protostellar population has a mean $\bar{{T}_{{\rm{d}}}}$ larger by $+1.8$ K compared with the quiescent population. We can estimate the expected uncertainty of this mean difference using the dispersion and size of each population, which gives $\sqrt{({3.8}^{2}/464)+({4.4}^{2}/788)}\approx 0.24$ K. Therefore, the difference is more than seven times the expected uncertainty. On the other hand, the difference of the means of $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ is $+0.17$ in favor of the protostellar population. The expected uncertainty in this case is $\sqrt{({0.25}^{2}/464)+({0.34}^{2}/788)}=0.017,$ that is, ten times smaller than the observed difference. We conclude that the observed differences of the quiescent and protostellar populations are statistically significant. Furthermore, the differences in temperature and column density are orthogonal to the expected uncertainty correlation (Figure 3), giving us more confidence that we are observing a real effect in both parameters. We confirm the significance of the difference using more sophisticated statistical tests in appendix. Note that either the temperature or the column density difference between other population pairs are larger than those between the quiescent and protostellar.

Are these statistical differences evident when comparing the $\bar{{N}_{{\rm{g}}}}$ and ${T}_{{\rm{d}},{\rm{P}}}$? On one hand, the difference between the mean ${T}_{{\rm{d}},{\rm{P}}}$ of the protostellar and quiescent samples is 2.6 K, while the expected uncertainty of this difference is 0.2 K. Therefore, despite the larger fitting uncertainties of ${T}_{{\rm{d}},{\rm{P}}}$ (see Table 3), we still detect a statistically significant difference between both populations when comparing only their central temperatures. On the other hand, the difference between the means of $\mathrm{log}\bar{{N}_{{\rm{g}}}}$ is 0.04 over an expected uncertainty of 0.02, that is, the difference is only two times the uncertainty. The latter is not highly significant, which is somehow expected because of the reasons explained at the beginning of this section. It is also expected from previous studies that indicate that neither the mass of the clumps (Hoq et al. 2013) nor their radii (Urquhart et al. 2014b) change conspicuously with evolutionary stage. This in turn implies that the average column density should remain approximately unchanged.

Within the quiescent sample, we identified in Section 3.6 a population of 83 clumps that appear as far-IRDCs at 70 μm. Of these, there are 77 associated with a single source along the line of sight. This sample has a mean and median $\bar{{T}_{{\rm{d}}}}$ of 14.9 and 14.7 K, respectively. The remaining quiescent population has mean and median $\bar{{T}_{{\rm{d}}}}$ equal to 17.2 and 16.4 K, respectively. Based on a Wilcoxon non-parametric test (Wall & Jenkins 2012, Section 5.4) we obtain a p-value of $4\times {10}^{-7}$ under the null hypothesis that these distributions are the same. Therefore, the temperature differences between the far-IRDCs and the rest of the quiescent clumps are significant. The column densities of the far-IRDC subsample are also larger compared to those of the rest of the quiescent sample. The far-IRDC $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ mean and median are −0.76 and −0.78, respectively, while for the remaining quiescent clumps they are −0.88 and −0.91, respectively. Again, we reject the null hypothesis (Wilcoxon p-value of $\sim {10}^{-5}$) and conclude that the far-IRDC sample is a colder and denser subsample of the quiescent population.

Finally, the Uncertain group (that is, MALT90 sources that could not be classified into any evolutionary stage) seems to be a mixture of sources in the four evolutionary classes, but associated with lower column densities (median $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}\sim 0.1$ g cm⁻²). Figure 8 shows that the $\bar{{T}_{{\rm{d}}}}$ values of the Uncertain group distribute almost exactly in between the other evolutionary stages. Neither a Wilcoxon nor a Kolmogorov–Smirnov (K–S) test can distinguish the $\bar{{T}_{{\rm{d}}}}$ distributions of the Uncertain sample from the remainder of the MALT90 sources combined with a significance better than 5%. Figure 8 also shows that the column densities of the Uncertain group are in general lower compared with those of any evolutionary stage except the PDRs. Molecular clumps with low peak column density may be more difficult to classify in the mid-IR, since they are probably unrelated to high-mass star formation. It is also possible that a significant fraction of these sources are located behind the Galactic plane cirrus emission and possibly on the far-side of the distance ambiguity, making the mid-IR classification more difficult and decreasing the observed peak column density because of beam dilution.

4.2.1. Column Density and Temperature Evolution in Previous Studies

4.2.2. Temperature and Column Density Contrasts

We analyze spatial variations of T_d and N_g by comparing their values at the peak intensity position with the average value in the clump. For each MALT90 clump, we define the temperature contrast and log-column density contrast as ${\rm{\Delta }}T=\bar{{T}_{{\rm{d}}}}-{T}_{{\rm{d}},{\rm{P}}}$ and ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}=\mathrm{log}(\bar{{N}_{{\rm{g}}}}/{N}_{{\rm{g}},{\rm{P}}}),$ respectively.

Table 5 lists the means and medians of the temperature and log-column density contrasts. Table 5 also gives 95% CIs¹⁴ for the medians of ${\rm{\Delta }}T$ and ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ per evolutionary stage, determined using the sign test (Ross 2004, Section 12.2). They were calculated using the task SIGN.test from the R statistical suite¹⁵ (version 3.1.1). The sign test is not very sensitive but it has the advantage that it is non-parametric and, in contrast to the Wilcoxon test (for example), it does not assume that the distributions have the same shape.

Table 5.  Temperature and Log-column Density Contrasts

Classification	${\rm{\Delta }}T$			${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$
	Mean	Median	95% CI	Mean	Median	95% CI
Quiescent	0.69	0.46	[0.39, 0.54]	−0.24	−0.23	[−0.24, −0.22]
Protostellar	−0.09	0.04	[−0.03, 0.10]	−0.38	−0.37	[−0.39, −0.36]
H ii region	−0.58	−0.33	[−0.43, −0.22]	−0.42	−0.42	[−0.43, −0.40]
PDR	0.79	0.54	[0.40, 0.75]	−0.36	−0.34	[−0.37, −0.33]

Note. ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}=\mathrm{log}(\bar{{N}_{{\rm{g}}}}/{N}_{{\rm{g}},{\rm{P}}})$ and ${\rm{\Delta }}T=\bar{{T}_{{\rm{d}}}}-{T}_{{\rm{d}},{\rm{P}}}$ (Section 4.2.2). For values inside the CIs we cannot reject the hypothesis that they are the median with a significance better than 5% using the sign test.

Download table as:  ASCIITypeset image

A negative ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ indicates that the clump has a centrally peaked column density profile with the absolute value of ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ being a measure of its steepness. As a reference, a critical Bonnor–Ebert sphere is characterized by ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}=-0.5515.$ Perhaps not very surprisingly, the ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ means, medians, and CIs are always negative, indicating that most of the clumps are centrally peaked. We find that the medians are significatively different between evolutionary stages, with no overlap in the CIs. The H ii region clumps are those associated with the steepest column density profiles, followed by the protostellar and the PDR clumps. Clumps in the quiescent evolutionary stage are associated with the smoothest column density profiles.

The temperature contrasts are also distinct for different evolutionary stages. A positive ${\rm{\Delta }}T$ indicates that the dust temperature increases away from the clump center, that is, dust temperature at the peak column density position (${T}_{{\rm{d}},{\rm{P}}}$) being lower than the average temperature ($\bar{{T}_{{\rm{d}}}}$). On the other hand, ${\rm{\Delta }}T$ is negative for decreasing temperature profiles. ${\rm{\Delta }}T$ is positive for quiescent clumps and PDRs, it is consistent with zero for the protostellar sources (temperature at peak similar to average temperature), and positive (peak column density warmer than average temperature) for the H ii region sample.

The previous sections have presented the differences between the temperature and column densities of the MALT90 groups. These differences are qualitatively consistent with the evolutionary sequence sketched in Jackson et al. (2013) that starts with the quiescent, and proceeds through the protostellar, H ii region, and PDR evolutionary stages. As Figure 9 shows, quiescent clumps are the coldest, in agreement with the expectation that these clumps are starless and there are no embedded young, high-mass stars. The far-IRDC subsample of the quiescent population is colder and denser on average compared to the rest of the quiescent clumps, and they might represent a late pre-stellar phase just before the onset of star formation. The mean temperature and $\mathrm{log}{N}_{{\rm{g}},{\rm{P}}}$ of the far-IRDC subsample are ∼15 K and −0.78, respectively.

To establish what fraction of the quiescent clumps might evolve to form high-mass stars, we use for now criteria defined by previous authors based on distance independent information, such as the column density; a more complete analysis will be done in Y. Contreras et al. (2016, in preparation). Lada et al. (2010) and Heiderman et al. (2010) propose that the star formation rate in a molecular cloud is proportional to the mass of gas with column densities in excess of ∼120 M_⊙ pc⁻² ($\sim 2.43\times {10}^{-2}$ g cm⁻²). Since there is considerable overlap in column density between MALT90 clumps that have different levels of star formation activity, we start by assuming that this relation gives the average star formation rate over the timescale of 2 Myr adopted by Lada et al. (2010) and Heiderman et al. (2010). We find that 98% of the quiescent clumps have $\bar{{N}_{{\rm{g}}}}\gt 120$ M_⊙ pc⁻², including all of the far-IRDCs, which suggest that most of these clumps will support some level of star formation activity in the future. Urquhart et al. (2014b) propose a column density threshold of 0.05 g cm⁻² for what they denominate "effective" high-mass star formation. This same threshold was recently proposed by He et al. (2015) based on a study of 405 ATLASGAL sources. Of the quiescent sample, 78% of the clumps have an average column density above this threshold, with the percentage increasing to 92% for the far-IRDCs. Based on these criteria, we conclude that virtually all quiescent clumps will develop at least low-mass star formation activity and that a large fraction ($\gt 70\%$) will form high-mass stars. On the other hand, López-Sepulcre et al. (2010) suggest a third column density threshold based on the observed increase of molecular outflows for clumps with column densities in excess of 0.3 g cm⁻². This column density is significatively larger than the previous thresholds, and only 3% and 6% of the quiescent and far-IRDCs populations, respectively, have larger average column densities. However, half of the clump sample of López-Sepulcre et al. (2010) have diameters $\lt 35^{\prime\prime}$ (the beam size of our column density maps) and more than a third have masses $\lt 200$ M_⊙, which indicates that the 0.3 g cm⁻² threshold may be pertinent for more compact structures than the clumps considered in this work.

The temperature and temperature contrast of the quiescent clumps are qualitatively consistent with equilibrium between the interstellar radiation field and dust and gas cooling (Bergin & Tafalla 2007). We find that quiescent clumps are the coldest among the evolutionary stages, but they are typically warmer (∼17 K) than expected from thermal equilibrium between dust cooling and cosmic ray heating alone (${T}_{{\rm{d}}}\sim 10$ K). We also find that the central regions of the quiescent clumps are in general colder than their external layers (${\rm{\Delta }}T$ negative). These characteristics are consistent with quiescent clumps being heated by a combination of external radiation and cosmic-rays. The quiescent sources also have the flattest density structure, with the largest ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ among all the other evolutionary stages. This is similar to the behavior found by Beuther et al. (2002), that is, the earliest stages of high-mass star formation are characterized by flat density profiles that become steeper as they collapse and star formation ensues.

The protostellar clump sample can be distinguished from quiescent clumps based on their column density and dust temperature. Protostellar clumps have larger column densities (∼0.2 g cm⁻²) and are slightly warmer (∼19 K). The central temperatures of the protostellar clumps also increase and become comparable to the temperature in their outer regions (${\rm{\Delta }}T\cong 0$). These characteristics indicate that protostellar clumps have an internal energy source provided by the HMYSOs. According to the results presented by Hoq et al. (2013), there is no significative difference in the distribution of masses between the quiescent and protostellar population. If we assume that this is also the case for the sample presented in this work (which will be confirmed in upcoming publications, see also He et al. 2015), then the most likely reason for the larger column densities of the protostellar sample compared with the quiescent sample is gravitational contraction. Because contraction develops faster on the densest, central regions, we expect the column density profiles to become steeper at the center of the clump. This is consistent with the observed decrease of ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}$ for the protostellar clumps compared with the quiescent clumps.

The H ii region sample is associated with the most negative temperature and column density contrasts (median of ${\rm{\Delta }}T=-0.33$ K and ${\rm{\Delta }}\mathrm{log}{N}_{{\rm{g}}}=-0.42$) compared with any other population, which indicates that H ii region clumps are very concentrated and they have a strong central heating source. This picture is consistent with the presence of a young, high-mass star in the center of the clump. The slight decrease of the peak column density compared with the protostellar phase could be explained because of the expansion induced by the development of the H ii region and the fraction of gas mass that has been locked into newly formed stars.

Finally, PDR clumps have the lowest column densities and largest temperatures among the four evolutionary stages. They are also associated with having colder temperatures toward the center compared to their outer regions. PDR clumps are possibly the remnants of molecular clumps that have already been disrupted by the high-mass stars' winds, strong UV radiation field, and the expansion of H ii regions. These molecular remnants are being illuminated and heated from the outside by the newly formed stellar population, but probably are neither dense nor massive enough to be able to sustain further high-mass star formation.