© ESO 2018

1 Introduction

The all-sky survey of the Planck satellite (Tauber et al. 2010), assisted by the IRAS 100 μm data, made it possible to identify and to catalogue cold interstellar clouds on a Galactic scale. The angular resolution of ~5′ is sufficient to identify compact sources that can be directly associated with the earliest stages of star formation. Analysis of the Planck and IRAS data led to the creation of the Planck Catalogue of Cold Clumps (PGCC; Planck Collaboration XXVIII 2016), which lists the properties of over 13 000 sources. At distances of up to ~ 500 pc, the Planck observations are able to resolve sub-parsec structures. For the more distant sources, up to 7 kpc, the observed signal may correspond to the detection of entire clouds. Because the PGCC detection method is based on the temperature of the sources, which is relatively cold compared to their environment (Montier et al. 2010), the low temperatures (typically T≲14 K) are a proxy for high column densities and also exclude the possibility of very strong internal heating. Even when they are spatially unresolved, the Planck clumps are thus likely to contain higher density substructure, including gravitationally bound pre-stellar and proto-stellar cores.

Some PGCC objects have been studied in detail using the Herschel follow-up observations made in the Herschel key programme Galactic Cold Cores (GCC). In the GCC, 116 fields were selected from an early version of the Planck clump catalogue (Planck Collaboration XXIII 2011), and the fields were mapped with Herschel PACS and SPIRE instruments (Pilbratt et al. 2010; Poglitsch et al. 2010; Griffin et al. 2010) between 100 and 500 μm (Juvela et al. 2010; Planck Collaboration XXII 2011; Juvela et al. 2012). The size of each Herschel map is ~40′. This means that the fields typically contain a few Planck clumps, but also a significant area of the surrounding cloud and even diffuse medium. This has enabled, for example, the study of the filamentary structure of the clouds around the clumps (Rivera-Ingraham et al. 2016, 2017; Malinen et al. 2016).

To study the variations in the submillimetre dust opacity, Juvela et al. (2015b) compared the 250 μm optical depths to the near-infrared (NIR) extinctions. The typical value of τ(250 μm)∕τ(J) in the fields was found to be 1.6 × 10⁻³, which is more than twice the values in diffuse medium (Planck Collaboration XI 2014), but consistent with other Planck studies of molecular clouds (Planck Collaboration Int. XIV 2014) (see also e.g. Martin et al. 2012; Roy et al. 2013; Lombardi et al. 2014). The submillimetre opacity was also found to be correlated with column density (and anticorrelated with temperature), with individual cores reaching values still higher by a factor of two.

Juvela et al. (2015a) studied the value and variations of the opacity spectral index β in GCC fields. Based on the IRAS and Planck data, the average value was found to be β = 1.84, somewhat higher than the typical Planck values for molecular clouds at larger scales (e.g. Planck Collaboration XXV 2011; Planck Collaboration Int. XIV 2014; Planck Collaboration XI 2014). The analysis of Planck data down to 217 GHz suggested a flattening of the dust emission spectrum (cf. Planck Collaboration Int. XIV 2014). There is strong evidence of spectral index variations within the individual fields, and the cold cores were associated with submillimetre spectral index values of β ~ 2.0 or even higher. Thus, the results are qualitatively consistent with many earlier observations of the T-β anticorrelation (Dupac et al. 2003; Désert et al. 2008; Rodón et al. 2010; Veneziani et al. 2010; Paradis et al. 2010; Etxaluze et al. 2011; Planck Collaboration XXII 2011; Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XIV 2014; Planck Collaboration XVII 2011; Planck Collaboration XIX 2011). In such studies, we recall that the colour temperature and the apparent β values are affected by both noise and temperature variations (Shetty et al. 2009b,a; Juvela & Ysard 2012a,b; Malinen et al. 2011; Ysard et al. 2012; Juvela et al. 2013; Pagani et al. 2015). In particular, to estimate the effect of noise on the apparent T-β anticorrelation, we need very precise knowledge of all error sources.

The above studies were limited to wavelengths below 500 μm or did include longer wavelengths but at a much lower resolution. The dust spectrum of cores can be further constrained by ground-based millimetre wavelength observations. However, for spectral index values, the results have shown a very wide range of β estimates, from ~1 to values far above 2.0 (Shirley et al. 2005; Friesen et al. 2005; Schnee et al. 2010; Sadavoy et al. 2013, 2016; Ward-Thompson et al. 2016; Chen et al. 2016). The large scatter is partly explained by the different wavelengths used, the spatial filtering resulting from sky noise reduction (when ground-based and space-borne observations are combined), and the varying nature of the sources (e.g. low values in regions with large temperature variations). Because of the additional problems in correlating the changes of the apparent dust spectra with the intrinsic dust properties, it is safe to say that our understanding of dust evolution across the star formation process is still very incomplete. This has implications for the estimates of the mass and the detailed structure of star-forming clouds. Accurate estimates of dust properties could also be related to laboratory measurements, to determine how dust properties vary according to the environment (Boudet et al. 2005; Coupeaud et al. 2011; Demyk et al. 2017).

In late 2015, we started the legacy survey SCOPE (SCUBA-2 Continuum Observations of Pre-protostellar Evolution) to this end, to map about 1000 PGCCs in 850 μm continuum emission with the SCUBA-2 instrument at the James Clerk Maxwell Telescope (JCMT; Liu et al. 2018, Eden et al. in prep.). Thousands of dense cores have been identified by the SCOPE survey, and most of them are either starless cores or proto-stellarcores with very young (Class 0/I) objects (see Liu et al. 2016; Tatematsu et al. 2017; Kim et al. 2017). These SCOPE sources are prime candidates for probing how pre-stellar/proto-stellar cores form and evolve, and for studying the very early stages of star formation across a wide variety of galactic environments. In addition, the 850 μm data are important in providing high-resolution information of the dust emission at the scales of individual cores and at a wavelength that according to previous studies lies in the regime between the steeper far-infrared part and the flatter millimetre part of the dust emission spectrum (Paradis et al. 2012).

In this paper we examine dust properties in a set of some 90 fields associated with PGCC sources. The fields were mapped at 850 μm with the SCUBA-2 instrument as part of the SCOPE pilot programme. The present sample contains the targets of the pilot programme for which Herschel photometric observationswere available. The goal of the paper is to examine dust SED variations as a function of the spatial scale and as a function of the wavelength.

The structure of the paper is the following. The observations and data reduction are described in Sect. 2 and the analysis methods are explained in Sect. 3. The main results are presented in Sect. 4. We discuss the results in Sect. 5 before listing the main conclusions in Sect. 6.

2 Observational data

2.1 SCUBA-2 observations

The target selection for SCUBA-2 observations was based on the PGCC catalogue, but also employed Herschel follow-up observations to pinpoint the exact locations of the column density maxima. Thus, all SCUBA-2 maps are not centred on PGCC coordinates, and one of the fields, G150.4+3.9A2, is not intersected by any PGCC clumps. The fields are listed in Table H.1.

The SCUBA-2 (Holland et al. 2013) maps at 850 μm were observed at JCMT between April and October 2015 (projects M15AI05 and M15BI061, PI: Tie Liu). The spatial resolution of the data is 14′′. All observations employed the CV Daisy mode (Bintley et al. 2014), resulting in maps with an approximate diameter of 12′ and a roughly constant integration time within the innermost 3′. The observations of all fields were of equal length with a total observing time of some 20 minutes per field and a median exposure time of some 80 seconds per map position. The observing conditions typically corresponded to τ(225 GHz) = 0.10 − 0.15. Maps with a 4.0′′ pixel size were created using the makemap task of the Starlink SMURF package (Jenness et al. 2011; Chapin et al. 2013). The final noise values are very similar in all fields. The approximate rms noise values are 5.6 mJy beam⁻¹ for the entire maps and 3.4 mJy beam⁻¹ for the central area within a 5′ radius.

In the map making, we used filtering scales of θ_F = 200′′ and θ_F = 500′′. SCUBA-2 data are typically reduced using external masking (or an auto-masking method) where source masks are derived from the SCUBA-2 data (Sadavoy et al. 2013; Moore et al. 2015; Pattle et al. 2015). However, we relied mainly on masks derived from Herschel maps, which provide a more complete picture of the cloud structure at lower column densities. The Herschel 500 μm maps were high-pass filtered with a Gaussian filter with an FWHM equal to 300′′, and the remaining emission was initially thresholded at 50 MJy sr⁻¹, which corresponds approximately to the 1σ noise of the SCUBA-2 maps (for typical values of dust temperature and spectral index). The resulting masks are thus significantly more extended than in auto-masking, where the threshold would be several times higher. When necessary, the mask was further adjusted so that it covered 10–30% of the map area, and each distinct masked area was at least 40′′ in size (larger than the Herschel beam). Alternative masks were created the traditional way, by thresholding preliminary 850 μm signal-to-noise ratio maps at a level of five. The masks were expanded by 15′′ because they were also used when Herschel data (with larger beam sizes) were processed through the SCUBA-2 pipeline, in order to match the spatial filtering that affects SCUBA-2 maps (Pattle et al. 2015). To check how the results depend on the mask size, we used a second set of Herschel-based source masks that were 50% smaller. In the following, these large and small Herschel masks are called LM and SM, respectively. The masks created using SCUBA-2 data are referred to as M850. The masks are shown in Fig. C.1.

In the following, we assume for SCUBA-2 a 10% relative uncertainty. This covers the uncertainty of the calibration and the uncertainty of the contamination by CO(3-2) line emission. Although the CO contribution in the 850 μm measurements can sometimes reach some tens of percent (Drabek et al. 2012), it is usually below 10% (e.g. Moore et al. 2015; Mairs et al. 2016) and can be expected to be lower for cold clumps. The uncertainties caused by CO emission are discussed in Appendix B.

2.2 Herschel observations

The Herschel SPIRE data at 250 μm, 350 μm, and 500 μm were taken from the Herschel Science Archive (HSA)¹. We used the level 2.5 maps produced by the standard data reduction pipelines and calibrated for extended emission (the so-called photometer extended map product). The quality of the SPIRE maps in HSA is good, and we find no need for manual data reduction. Out of the 96 fields, 43 were observed as part of the GCC project (Juvela et al. 2010). For consistency, we used HSA pipeline data for these fields as well.

The resolutions of the SPIRE observations are 18.4′′, 25.2′′, and 36.7′′ for the 250 μm, 350 μm, and 500 μm bands, respectively². The maps were convolved to 40′′ and fitted with modified blackbody (MBB) curves with β = 1.8. These initial SED estimates were used to colour correct the SPIRE maps to obtain monochromatic values at the nominal wavelengths. The SCUBA-2 data were later colour corrected using the same SEDs. In the temperature range of T = 10 − 20 K, the corrections are 0–2% for both SPIRE and SCUBA-2. The effect on the SED shapes and especially on the β estimates are even smaller because the corrections are correlated between the bands (see Appendix A). We adopted a relative uncertainty of 4% for the SPIRE bands, with a ρ = 0.5 correlation between the bands (Bendo et al. 2013). The assumed correlation between the photometric values decreases the formal T and β uncertainty in the fits to SPIRE data (Galametz et al. 2012). However, when SPIRE is combined with SCUBA-2, the situation may no longer be true. If the 250–500μm data have errors preferentially in the same direction, the effect on the SED shape is large when SPIRE is combined with the (presumably) uncorrelated850 μm data point.

The main analysis was made without PACS data, which were only compared to the SED fits of longer wavelength measurements. The reason is that PACS observations cover only a small fraction of the SCUBA-2 fields and including shorter wavelengths can bias estimates of the spectral index, especially in the presence of embedded radiation sources or emission from very small grains (e.g. Shetty et al. 2009b; Malinen et al. 2011; Juvela & Ysard 2012b). Of the 96 fields, 37 are fully or partially covered by PACS observations (for observation numbers, see Table 1). We used the UNIMAP versions of the 160 μm maps provided in the Herschel science archive (Piazzo et al. 2015). We assumed a relative uncertainty of 10% for PACS to cover the 5% calibration uncertainty (estimated for point sources) (Balog et al. 2014) and the map-making uncertainty.

2.3 Planck and IRAS data

We used IRAS and Planck observations to examine the dust emission at large ~ 10′ scales. We used the IRIS version of the 100 μm IRAS data (Miville-Deschênes & Lagache 2005) and the 857 GHz, 545 GHz, 353 GHz, and 217 GHz Planck maps taken from the Planck Legacy Archive³ that correspond to the 2015 data release (Planck Collaboration I 2016) where the cosmic microwave background (CMB) emission has been subtracted. We subtracted the estimated levels of the cosmic infrared background (CIB) from Planck maps. Because the targets have high column densities, the CIB correction ( ~1 MJy sr⁻¹ or less in all bands) is relatively unimportant (Planck Collaboration VIII 2016). The amount of emission from stochastically heated very small grains (VSG) at 100 μm is unknown. However, it could amount to some tens of percent in cold clouds (Li & Draine 2001; Compiègne et al. 2011). We did not correct for the VSG emission. Thus, the combined SED of IRAS and Planck data is expected to correspond to a higher colour temperature (and slightly altered and probably lower spectral index values) than the SED of the large grains alone (Juvela et al. 2015a).

Planck 353 GHz and 217 GHz bands are contaminated by CO J = 2 − 1 and J = 3 − 2 line emission, which is a significant source of uncertainty in the estimates of the long-wavelength dust emission spectrum. The effect is more significant at 217 GHz because J = 2 − 1 is usually the stronger of the two CO lines, while dust emission at 217 GHz is only about one fifth of the emission at 353 GHz (assuming β = 1.7). We corrected these two bands in the same way as in Juvela et al. (2015a). The correction made use of the Type 3 CO maps provided by the Planck Consortium and was based on the assumption of line ratios T_A (2–1)/T_A(1–0) = 0.5 and T_A (3–2)/T_A(1–0) = 0.3. The justification and uncertainty of these values is discussed in Juvela et al. (2015a).

The 217 GHz band could have additional contribution from free-free emission in the case of actively star-forming clumps. Therefore, similar to Juvela et al. (2015a), we subtracted a model for the low-frequency foregrounds⁴, although this is in practice much smaller than the CO correction.

We assumed absolute calibration uncertainties of 10% for IRAS 100 μm, 5% for the 857 GHz and 545 GHz Planck channels, and 2% for the 353 GHz and 217 GHz channels. Because of the common calibration, we furthermore assumed a correlation of ρ = 0.5 between the first two Planck bands and again between the last two bands. In the 353 GHz and 217 GHz bands, we included an additional uncertainty (added in squares) that corresponds to 30% of the applied CO correction.

2.4 Cloud distances

A distance estimate has been published for 86 of the 96 sources in the PGCC catalogue (Planck Collaboration XXVIII 2016). These are mostly based on 3D extinction mapping, the comparison between the observed stellar reddening, and the predictions derived from a model of the Galactic stellar distribution (Marshall et al. 2006, 2009). There is also considerable overlap between the current source sample and the fields observed as part of the Herschel key programme Galactic Cold Cores (GCC). The distances of the GCC fields were discussed in detail by Montillaud et al. (2015). With only one exception, these estimates agreed with the original PGCC values. The distance estimates are listed in Table H.2, which also notes some nearby molecular clouds and their distance estimates, if available in the SIMBAD database⁵ or in the cloud compilation of Dutra & Bica (2002). The distance estimates typically have uncertainties of several tens of percent and therefore are a major source of uncertainty when the masses of the sources are estimated.

We revisited the extinction-based distance estimates using the method MACHETE (Marshall, in prep.) The method is based on the observed reddening of background stars, which is compared to predictions based on the Besançon model (Robin et al. 2003) of stellar distributions and the optimised model of the dust distribution along the line of sight. Apart from the improvements in the method itself, we also used Herschel data to better separate the high-column-density regions associated with the SCUBA-2 targets. In practice, we used masks that covered an area where τ(250 μm) is above its 40% percentile value (separately within each map). The masks are typically larger than the size of the individual PGCC clumps and thus should provide better statistics for the background stars.

3 Methods

3.1 Clump extraction

Clumps were extracted from SCUBA-2 850 μm signal-to-noise ratio maps with the Fellwalker method (Berry 2015), which is also used by the SCOPE project (Liu et al. 2018). Clumps were required to contain more than 10 pixels and to rise above three times the rms value σ. The other relevant parameters of the algorithm were MINDIP = 2σ and MAXJUMP = 5 pixels⁶. The clump extraction was repeated for each version of the reduced data (e.g. the three different source masks, LM, SM, and M850). The extraction can result in spurious detections near the noisy map boundaries, but these are rejected later when the analysis is limited to structures within 5′ of the centre of each field.

3.2 Clump photometry

Clump fluxes were measured from the original Herschel maps and from SCUBA-2 850 μm and Herschel maps processed through the SCUBA-2 pipeline. To make the procedure consistent across maps of originally different spatial resolution, the photometric measurements were made on maps that were first all convolved to a common resolution of 40′′.

The measurement apertures were based on the footprints of the Fellwalker clumps, but were extended by 20′′ to reduce the flux losses resulting from the above-mentioned map convolution. The apertures are shown in Fig. C.1.

All flux measurements employed a local background subtraction. The reference annulus extended from 40′′ to 60′′ outside the original Fellwalker clump footprint. The annulus was thus non-circular and followed the clump shape. All pixels that were closer than 20′′ to the footprint of another clump were excluded. We used the median of the remaining annulus pixels as the background estimate and did not explicitly interpolate the background over the aperture.

For comparison and to create a subsample of more compact cores, alternative photometry was carried out using the same reference annuli, but 50% smaller apertures. These covered the pixels with τ(250μm) values (estimated using Herschel data, see Sect. 2.2 and Appendix C) above the median value of the original aperture. Apertures smaller than 256 arcsec² (16 pixels with a size 4′′× 4′′ each) were rejected. In the following, the original and the high-column-density apertures are called τ_all and τ_high apertures, respectively.

The statistical errors of the flux values were estimated based on the standard deviation of the pixel values in the annulus and the sizes of the annulus and the aperture. We took the fact into account that the maps were convolved to a resolution of 40′′ and scaled the standard deviation of the pixel values with the square root of the ratio between the area of the measuring aperture (or annulus) and the effective area of the 40′′ beam. Because of real surface brightness fluctuations within the reference annuli, the photometric errors tend to be overestimated. To reduce this effect, we removed the signal that is linearly correlated between the bands before estimating the noise. Each band was in turn fitted against a reference, which was the average of all the other bands. The correlated signal predicted by the least-squares fit and the reference data was then subtracted. The emission is not perfectly linearly correlated, mainly because of temperature variations. Therefore, the estimates remain an upper limit of the actual noise. In SED fits, the relative weighting of the bands also took calibration uncertainties into account (4% and 10% for SPIRE and SCUBA-2, respectively).

3.3 SED fits

We fitted SEDs using three or four bands, combining SCUBA-2 measurements and fluxes extracted from SPIRE maps processed through the SCUBA-2 pipeline. We also fitted the fluxes from the original, unfiltered SPIRE maps without SCUBA-2 data, to see how the SPIRE SEDs are affected by the filtering process itself.

The fitted model of flux density is the modified blackbody (MBB) function $$F(\nu )=F({\nu }_0)\frac{B_{\nu }(\nu ,T)}{B_{\nu }({\nu }_0,T)}{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }=\kappa (\nu )B_{\nu }(\nu ,T)M/d^2,$$(1)

where ν is the frequency, ν₀ the selected reference frequency, B_ν the Planck law, β the opacity spectral index, and d the cloud distance (Hildebrand 1983). The three fitted parameters are the flux density at the reference frequency F(ν₀), the dust emission colour temperature T, and the spectra index β. The assumption of a value for the dust opacity per unit mass κ(ν) enables estimating the clump masses M. We assumed a value $\kappa =0.1{(\nu /1\text{\hspace{0.17em}THz})}^{\beta }$ cm² g⁻¹ (Beckwith et al. 1990), as was used, for example, in papers of the Galactic Cold Cores project (Juvela et al. 2012; Montillaud et al. 2015)and of the Herschel Gould Belt Surveys (André et al. 2010).

Equation (1) describes the apparent shape of an SED in terms of its colour temperature and apparent spectral index. If the SED is assumed to be connected to real dust properties (such as the intrinsic opacity spectral index) and clump properties (such as the mass), one implicitly assumes that the emission is optically thin and dust is characterised by a single temperature and a single β value. The impact of these assumption is discussed further in Sect. 5.

We used only a subset of all the extracted clumps for the SED analysis (see Sects. 3.1 and 3.2). To avoid clumps near high-noise map boundaries, the centre of the included clumps was required to reside within 5′ of the field centre. We also rejected clumps whose 850 μm flux was below 0.1 Jy or where the MBB fit to four bands (850 μm and the filtered SPIRE data) resulted in poor fits (see below).

4 Results

We investigated differences of dust emission properties in dense, potentially prestellar clumps and the surrounding more diffuse clouds. We examined the values of the opacity spectral index β of the clumps. We sought evidence of correlated changes between the dust temperature and β and a potential wavelength dependence of β. These could be indications of dust processing during the star formation process or might be directly related to the local physical conditions such as the radiation field and the temperature. We start by examining large-scale dust emission of the selected fields (Sect. 4.1) before concentrating on the compact SCUBA-2 clumps and cores (Sect. 4.2). We compare thefits of different wavelength ranges up to the millimetre regime, examine the effect of the extension to the shorter 160 μm wavelength (Sect. 4.3), and check the correlations between young stellar objects and the clump and dust emission properties (Sect. 4.4).

4.1 Dust spectrum at large scales

Over the entire fields in which the PGCC clumps reside, we analysed IRAS, Herschel, and Planck observations in the wavelength range λ = 100 − 1380μm (from 3000 GHz to 217 GHz) to characterise the dust emission at large scales. With IRIS data and the zero-point-corrected Herschel maps, we were able to directly use surface brightness values averaged within a radius of r = 8′ of the centre of the SCUBA-2 fields. Alternatively, we could subtract the background that is estimated as the mean surface brightness in a selected reference area. In that case, we added the uncertainty of the background subtraction. This is described by the covariance matrices that are estimated from the values in the reference region, from the part of the signal that is not linearly correlated between the bands. The procedure to exclude the correlated signal is the same as in Sect. 3.2.

To ensure consistent handling of all bands, the maps were convolved to a common resolution of 6′. We carried out MBB fits using different band combinations to check the consistency of the data and to examine possible wavelength dependencies.

4.1.1 SEDs of IRAS and Planck data

Figure 1a shows the colour temperature T and the spectral index β calculated from the combination of the 100 μm IRAS data and the Planck bands at 857 GHz, 545 GHz, 353 GHz, or 217 GHz (wavelengths 350 μm, 550 μm, 850 μm, and 1380 μm, respectively). There is no subtraction of the local background. The confidence regions plotted in Fig. 1 were calculated with the Markov chain Monte Carlo (MCMC) method, using flat priors with 8 K < T < 25 K and 0.5 < β < 3.

The median values of the 100–850 μm fit are T = 15.74 K and β =1.73. The addition of longer wavelengths up to 1380 μm increases the temperature to T = 15.95 K and lowers the spectral index to β = 1.69. Because long wavelengths are less sensitive to line-of-sight temperature variations, we might have expected the colour temperature to decrease rather than increase (assuming that the opacity spectral index does not change withwavelength). However, differences are not significant compared to the statistical uncertainties (~ 1 K in temperature and ~ 0.1 in the spectral index, see Fig. 1) and might easily be caused by small systematic errors. Figure 1b illustrates the effect of CO corrections on the 100–1380 μm fits. The difference in Fig. 1a might be explained if we assume that the CO corrections are underestimated by some 10%. Although there is no particular reason to suspect that the correction is underestimated (see below), the evidence for a wavelength dependence of β is tentative at most in Fig. 1a.

To further quantify the reliability of the CO corrections, we compared Planck CO estimates to ground-based data that provide a partial coverage of about half of the fields. Meng et al. (2013) and Zhang et al. (2016) made CO(1–0) observations at the Purple Mountain Observatory (PMO) radio telescope, fully or partially covering 17 of our SCUBA-2 fields. In Fig. 2b we show the correlation with the Planck CO estimates after convolving the PMO observations to the same 7′ resolution and scaling with the assumed main beam efficiency of η_MB = 0.5. The correlation is generally good, but the ground-based values are higher on average by ~30%.

The ongoing SAMPLING⁷ (Wang et al. 2018) survey (PI. Ke Wang) at the Arizona Radio Observatory 10 m SMT telescope is mapping a large number of PGCC targets in the J = 2 − 1 transitions of CO isotopomers and thus provides more direct estimates of the CO contamination of the 217 GHz band. There are 56 SMT maps near our targets. Because the maps are only ~ 6′, we used the Planck CO maps at their original resolution of 5.5′ in the comparison. We convolved the SMT data to this resolution and compared values at the centre positions of the SMT maps (Fig. 2a). The SMT values are the sum of the ¹² CO(2–1) and ¹³ CO(2–1) lines and are plotted against 0.5 times the Planck CO(1–0) estimates. This corresponds to the assumption of the (2–1)/(1–0) = 0.5 line ratio that was used in the CO correction of the Planck 217 GHz data. The agreement is very good and shows that the performed CO corrections were at the appropriate level. The convolution of the SMT data was only possible to calculate using the extent of the observed map, which may cause some overestimation of the SMT values (assuming that the CO emissiondecreases with distance from the SCUBA-2 clump). When we assume consistent calibration of the SMT and PMO observations, the comparison of the two frames in Fig. 2 shows that the CO line ratio (2−1)/(1−0) tends to be below 0.5. This would agree with the interpretation that most of the objects are very cold and embedded in more extended and CO-bright envelopes. This also suggests low CO contamination for the SCUBA-2 850 μm data.

4.1.2 SEDs of Herschel and Planck data

Because the fields were selected based on the availability of Herschel data, we repeated the analysis using only wavelengths λ ≥ 250 μm. Figure 3 shows the results from the fits to Herschel (250 μm, 350 μm, and 500 μm) and Planck data. The maps were convolved to the resolution of 6′, and the surface brightness was averaged within a radius of r = 8′. We excluded three fields that are only partially covered by Herschel. The relative weighting of the bands follows the uncertainties discussed in Sects. 2.2 and 2.3, including a 30% uncertainty for the CO correction of the two lowest frequency Planck channels.

Results of Fig. 3a were calculated without background subtraction. Because the intensity zero-points of the Herschel maps are based on Planck data, the two data sets are not completely independent. The median values of the SPIRE fits are T = 14.70 K and β = 1.80. The results are practically identical when the Planck data are included, T =14.75 K and β = 1.79. Compared to Fig. 1, the temperatures are lower by 1 K and the spectral index values are higher by more than 0.05 units. Thesechanges are not larger than the estimated uncertainties given in Figs. 1 and 3.

Figure 3a includes examples of confidence regions. Regardless of the correctness of the magnitude of the error estimates, the plot shows the direction along which the T and β values are partly degenerate. Figure 3b shows the results when estimates of the local background were subtracted from the flux values. The background region was defined by taking a 24′ × 24′ area centred on the SCUBA-2 map and rejecting pixels outside of or closer than 3′ to Herschel map borders. The reference area was defined by the remaining pixels where the 857 GHz surface brightness was below its 10% percentile. We subtracted the mean surface brightness value of the reference area at each wavelength. The scatter of the (T, β) estimates increased, as expected, because of the smaller residual emission above the background. For the median, the β value is now lower for the 250–500 μm fits and higher for the 250–1380 μm fits. This would be contrary to the expectation of β decreasing at millimetre wavelengths.

All PGCC clumps with centre coordinates within a radius of 5′ of the SCUBA-2 map centres are marked in the maps of Appendix C. They are also listed in Table H.3, including the PGCC values of the spectral index and the colour temperature (with β =2 and with free β). The PGCC values result from fits to the cold emission component, after the subtraction of the warm component that is defined by the IRAS 100 μm values and the average SED of the region. Thus, the PGCC temperatures are significantly lower than in the analysis above, and the spectral index values are correspondingly higher. For the clumps of Table H.3 for which T and β estimates are available, the quartile points (25%, 50%, and 75% percentiles) are 10.7 K, 11.6 K, and 12.3 K for the colour temperature and 2.01, 2.18, and 2.37 for the spectral index. The β values for the PGCC clumps are thus higher than the values of the entire fields (Fig. 1). This is not yet direct evidence of a change in the intrinsic dust properties. By averaging emission over large areas, we also have averaged SEDs with different temperatures, and this can lead to lower apparent β values. Therefore, it is necessary to also investigate the spectra at higher spatial resolution.

Fig. 1 Results of MBB fits to IRAS and Planck data averaged within 8′ of the centre of the SCUBA-2 fields. Frame a shows the fits completed with 100–550 μm, 100–850 μm, and 100–1380 μm bands. Each marker corresponds to one field, and the dotted lines connect estimates of the same field. The median values calculated over all the fields are plotted with large symbols. The median values and the 1σ dispersion estimated from the interquartile ranges are quoted in the frame. Frame b shows the fits to the 100–1380 μm data without CO corrections (ΔCO = 0), with the default CO correction (ΔCO = 1), and with twice the default correction (ΔCO = 2). The median values are indicated by the three large triangles. The contours show the 67% (solid line) and 95% (dashed line) confidence regions for a 100–850 μm (frame a) and 100–1380 μm fit with ΔCO = 1 (frame b) calculated using the median error covariance matrix of the fields.

Fig. 2 Comparison of CO estimates from Planck Type 3 CO maps and ground-based observations. Upper frame compares the sum of 12CO(2–1) and 13 CO(2–1) values of 56 SMT maps to Planck CO estimates for the assumed line ratio of (2−1)/(1−0) = 0.5. Lower frame shows comparison with CO(1–0) observations from the PMO telescope, with data at 7′ resolution sampled at 3′ steps. Each colour and symbol combination corresponds to one of 17 fields. Solid lines show the one-to-one relations, and dotted lines correspond to slopes different by ±30%.

Fig. 3 Results of MBB fits to Herschel and Planck data averaged within 8′ of the centre of the SCUBA-2 fields. The results are shown for absolute surface brightness measurements (frame a) and after subtraction of the local background emission (frame b). The symbols correspond to the band combinations listed in the legend in frame b. The large symbols denote the median values over the fields. In frame a, the contours show the 67% (solid lines) and 95% (dashed lines) confidence regions for 250–500 μm (red) and 250–1380 μm (blue) fits, using error covariance matrices that are the median over all fields. For the SPIRE-only fits, the adopted error estimates clearly overestimate the true uncertainty.

4.2 Dust spectrum at small scales

We investigated the small-scale dust emission using the combination of Herschel 250–500 μm and SCUBA-2 850 μm data and clumps extracted with the Fellwalker method (see Sect. 3.1). We used four clump samples that correspond to (1) θ_F = 500′′ and LM source masks, (2) θ_F = 200′′, SM masks, and τ_high part of the clumps, (3) θ_F = 200′′ and M850 masks, and (4) data reduction with θ_F = 200′′ and without external source masks. The four alternatives were used to examine the robustness of the results with respect to the details of the data reduction and clump selection. After the rejection of clumps farther away than 5′ from the field centres or with S(850μm) values below 0.1 Jy, there remain about 200 clumps, the number depending on the data version. Many SED fits are not reliable, and we carried out a further selection based on the χ² values of the fits.

Figure 4a–c shows the β distributions for the reduction with θ_F = 500′′ and LM masks, using the full clump apertures (τ_all). The clumps were divided into subsamples where the χ² values of the fits were within the best P(χ²) = 85%, 25%, 10%, or 5%. Thus, P(χ²) = 5% corresponds to the strictest selection criteria (least clumps) with χ² values below the 5% percentile point. The first two frames show the results for fits to Herschel 250–500 μm bands, before and after filtering by the SCUBA-2 pipeline. The filtering could lead to different results if the filtering procedure introduces additional uncertainty or bias, or if the true SEDs are affected by the removal of large-scale emission.

In the following, we refer mainly to the numbers for the P(χ²) = 25% samples. For the original SPIRE data (the θ_F = 500′′ case), the median value is β = 1.90. Based on the distribution of the β values, the error of the mean is smaller than 0.05. The filtering increases the width of the distributions, partly because of the larger relative uncertainty of the lower flux values and possibly via larger background fluctuations. If the densest parts of the clumps have lower colour temperatures and higher β values than their environment, we might expect the filtering to increase in the β estimates. Figure 4 shows only weak indications of such a trend. Figure 4c shows the results for fits to combinedHerschel and SCUBA-2 data, with self-consistent filtering across all bands. The median value has decreased to β =1.63. The same trend exists for the P(χ²) = 5% sample, although it is not very reliable because of the low number of clumps.

Figure 4d–f shows a second analysis that is a priori very far from the previous case, using a θ_F = 200′′ filter scale, SM masks, and fluxes integrated over the τ_high apertures. More of the large-scale emission is filtered out, and the smaller source mask and clump size further suppress the contribution of extended emission. The number of clumps is lower than in the previous case, and the β distributions are also slightly more narrow. Figure 5 shows the SEDs for the clumps in the P(χ²) = 5% sample. In six out of eight cases, the addition of the 850 μm has decreased the β estimate. Table H.4 lists thefit parameters for the P(χ²) = 25% sample.

Figure 4g–i shows a case where maps were made using M850 source masks and θ_F = 200′′ filtering scale. Compared to the previous cases, the clump selection is thus more directly based on the SCUBA-2 850 μm data. However, the result for the unfiltered SPIRE data shows that the underlying clump samples are still relatively similar. The inclusion of the 850 μm point again decreases the median β for both the P(χ²) = 5% and P(χ²) = 25% samples.

The final plots in Fig. 4j–l correspond to maps made with θ_F = 200′′ but without anyexternal masks. This is expected to favour the extraction of the most compact sources, although the total number of sources is not lower than in the previous two cases. When filtering is applied to SPIRE data, the β distribution moves towards higher values and becomes wider. We interpret this mostly as a sign of a general increase in uncertainty. In the combined fits of SPIRE and SCUBA-2 data, the spectral index values are also much higher than in the SPIRE-only fits.

A table of the fit parameters is provided in Appendix F, where the uncertainty of the 850 μm calibration is also briefly discussed.

Figure 6 uses clumps from Fig. 4a–f and compares the fits to the original SPIRE data and to the combination of SCUBA-2 and filtered SPIRE data. Figure 6a shows an example of the confidence regions estimated with the MCMC method. It indicates the direction of the error ellipses, but the size of the confidence region varies from clump to clump. The figure suggests that a large fraction of the anticorrelation between T and β can be caused by the observational uncertainties.

The dependence of β values on other parameters and directly on error estimates is examined further in Appendices E and F. The latter section also further quantifies the potential effects of the 850 μm calibration uncertainty.

4.3 SEDswith the 160 μm point

Although 37 fields are partially covered by Herschel PACS observations, the number of clumps with 160 μm data is limited. Figure 7 examines the 23 clumps in the P(χ²) = 25% sample of clumps from the maps made with θ_F = 500′′ and LM masks.

The 160 μm values are typically above the SPIRE fit (Fig. 7a), which might be expected if each source contains a wide range of temperatures. Figure 7b compares (T, β) values for Herschel 250–500 μm and 160–500 μm MBB fits. Both fits cover a similar region in the (T, β) plane, except for a few clumps where the inclusionof the 160 μm point leads to much higher temperatures and lower β values. The number of these sources is small enough so that the shift in the median parameter values remains small with Δ T = +0.4 K and Δβ = −0.1. For simplicity, both fits were carried out without considering error correlations between the different bands. The SPIRE fits are therefore not identical to those in Fig. 4a. The difference in the median value of β is 0.08 units.

Fig. 4 Distributions of β values in fits to Herschel and SCUBA-2 850 μm observations of clumps. Results are shown for original SPIRE data (first row), filtered SPIRE data (second row), and the combination of filtered SPIRE and 850 μm data (third row). The columns correspond to different versions of data reduction, as indicated above the frames. The white, blue, yellow, and red histograms correspond to P(χ2) values of 85%, 25%, 10%, and 5%, respectively (see text). The number of clumps in each sample and the quartile values (25%, 50%, and 75% percentiles) of the blue and red histograms are given in the frames.

Fig. 5 SEDs of eight sources from the P(χ2) = 5% subsampleof SM clumps in Fig. 4f. Plots include fits to the original SPIRE data (upper red curve and red numbers) and to the filtered SPIRE data plus the SCUBA-2 850 μm point. The clump numbers given after the field name are the same as in Fig. C.1. The 160 μm data points are available for one of the clumps: it is plotted in the figure, but is not part of the fit.

Fig. 6 Spectral index β vs. colour temperature T for the clump samples of Fig. 4a–c (left frame) and Fig. 4d–f (right frame). The black symbols correspond to fits to the original SPIRE data, and the blue symbols to the combination of SCUBA-2 and filtered SPIRE data. Only clumps from the P(χ2 )=25% samples are included. The median values are indicated with large crosses. Frame a includes a contour for a typical 67% confidence region in the four-band fits.

Fig. 7 Comparison of Herschel MBB fits with and without the 160 μm data point inthe case of the P(χ2) = 25% clump sample with θF = 500′′ and LM masks. The left frame shows the distribution of ratios between the 160 μm flux density and the value predicted by the 250–500 μm fit. The greyhistogram contains all clumps, and the black histogram is based on three clumps with YSO candidates. The second frame shows (T, β) values for individual clumps for the 250–500 μm (blue triangles) and 160–500 μm (red squares) fits. The circles indicate clumps with YSO candidates. The median values are quoted in the frame and are plotted withlarger symbols.

4.4 Young stellar objects

We examined the spatial distribution of young stellar objects (YSOs) using the catalogue of YSO candidates of Marton et al. (2016) that is based on the analysis of 2MASS (Skrutskie et al. 2006) and WISE (Wright et al. 2010) data.

Figure 8 shows the surface density as a function of the distance from the centre of the SCUBA-2 fields. The density of Class III candidates is practically constant, about 3–4 × 10⁻³ arcmin⁻². The density drops only within the innermost 7′, possibly because of the observational bias resulting from the increasing extinction. Except for this, Class III sources are not correlated with the clumps. The total number of Class I-II sources is higher. Their density also increases significantly towards the centre of the fields, reaching 2.8 × 10⁻²arcmin⁻² within the innermost 2′ radius. The highest concentration of Class I-II candidates is found in the field G130.1+11.0A1, where the Herschel coverage is unfortunately incomplete.

We separately calculated the YSO surface densities for all clumps and for the subsample selected for the SED analysis (see Sect. 4.2). For Class I-II candidates, the trend seen in radial profiles continues, and they are found preferentially inside the clumps (see Fig. 8). The surface density is 7.7 × 10⁻² arcmin⁻² within the P(χ²) = 85% clump sample. It is also significant that none of the Class III candidates falls inside a clump.

In the LM case, the P(χ²) = 25% sample includes 56 clumps, 5 of which had Class I-II YSO candidates. In SPIRE fits, the YSO clumps have median values of T =14.5 ± 0.5 K and β =1.62 ± 0.14, while starless clumps haveT = 11.9 ± 0.9 K and β = 2.1 ± 0.3, where the quoted uncertainties are median absolute deviations of the distributions. For the 250–850 μm fits, the corresponding values are T = 14.8 ± 1.0 K and β = 1.58 ± 0.08 vs. T = 14.0 ± 1.2 K and β = 1.59 ± 0.23. There is thus a tendency for the proto-stellar clumps to be warmer, while the evidence of lower β values is not altogether consistent. The sample of proto-stellar source candidates also comprises only six sources. In Fig. 7 there were only three clumps with YSO candidates that also had PACS observations at 160 μm. For this small sample, the clump temperatures were also higher than average. The ratio between the 160 μm measurement and the fit to the longer wavelength bands was also higher for clumps with YSO candidates. However, the effects of YSO heating appear to be small (Fig. 7a).

The observed differences in the apparent (T, β) values do not necessarily mean changes in the intrinsic dust properties. Large temperature variations, such as caused by internal heating, could lower the β values estimated from the SEDs.

Fig. 8 YSO statistics. The surface density of sources is shown as a function of the distance ( >2′) from the field centre. The mean values over all fields are shown for Class I-II (thick red line) and Class III (thick yellow line) candidates. The thin red lines correspond to Class I-II sources in individual fields. The shading corresponds to a region with less than one source for a given distance. The two overlapping stars indicate the density of Class I-II candidates in all clumps (all fields) selected by the Fellwalker algorithm (red star) and in clumps belonging to the P(χ2 ) = 85% sample (blue star).

4.5 Clump masses

Most fields have distance estimates in Table H.2, which allows calculating the physical size, mass, and volume density of the clumps. Figure 9 shows results for the LM data, the P(χ²) = 85% clump sample. The masses correspond to fluxes inside measurement apertures of area A, and the volume densities are calculated for a spherical volume with an effective radius r for which A =πr².

The mass and distance are correlated because the clump detections are affected by our sensitivity to low flux densities. In Fig. 9a the dashed line corresponds to a mass threshold that assumes a 3.4 mJy beam⁻¹ noise level and a 5σ detection over a single 850 μm beam. The dotted line corresponds to a 850 μm flux density of 0.1 Jy for both the T = 13 K and β = 1.8 spectrum. The correlation in Fig. 9a thus mainly reflects the flux density threshold combined with the geometrical effect of the distance. All sources are relatively close to the threshold, and we do not have nearby sources that would also be very massive. This is natural because low-mass objects are generally much more numerous, and on the other hand, the spatial filtering of the SCUBA-2 data and the clump detection algorithm both limit the possibility of detecting very extended sources. The mass-distance relation is also affected by the large uncertainty of the distances and the functional dependence between the distance and the mass estimates.

The temperature is not dependent on the field distance, but there is some negative correlation between the volume density and the clump temperature. However, negative errors in temperature will naturally lead to positive errors in column density, mass, and volume density. Higher temperature also directly enables the detection of lower density objects. Furthermore, the fluxes were measured from apertures that were extended 20′′ beyond the original SCUBA-2 clump detections in order to accommodate the lower resolution of the SPIRE data. The modified clump size can introduce some bias that is also dependent on the temperature. For the same sample of clumps, the correlation between the temperature and the column density remains very weak.

We have marked the clumps with Class I or II YSO candidates in the plot. The median distance of YSO sources is smaller than that of the full clump sample, but this is not significant considering the small number of YSO candidates. As noted in Sect. 4.4, there is some tendency for proto-stellar clumps to be warmer than clumps on average. In Fig. 9, the proto-stellar clumps do not stand out based on either their mass or density.

We note that while clumps with embedded proto-stars may be typically among the warmer cores in the sample, they obey a similar temperature/density relation as the starless clumps. This is in contrast to the behaviour seen on the scale of individual star-forming cores (with sizes ~10⁻² pc): combined SCUBA-2 and Herschel observations of starless cores in the Taurus molecular cloud (Ward-Thompson et al. 2016) and in the Cepheus Flare (Pattle et al. 2017) have shown that cores heated by nearby or embedded protostars are significantly warmer than cores of the same density that are not affected by proto-stars. The lack of a systematic elevation in temperature in Planck clumps with embedded proto-stars relative to starless Planck clumps of similar densities suggests that the heating effect of a protostar on its environment is limited to a relatively small volume surrounding the proto-star itself, and so that embedded proto-stars do not significantly raise the average dust temperature of their parent clumps.

Fig. 9 Clump mass vs. field distance (frame a) and average volume density vs. temperature (frame b). The plots correspond to the LM data version and the P(χ2 ) = 25% clump sample. The fluxes and temperatures are from the 250–850 μm fits. The data are further divided into three samples based on temperature (frame a) or mass (frame b). The clumps with YSO Class I or II candidates are marked with circles. Frame a shows the thresholds corresponding to a 17 mJy per 850 μm beam (dashed line) and to a 850 μm flux density of 0.1 Jy f(dotted line; see text).

5 Discussion

We have investigated submillimetre dust emission in 96 fields observed with the SCUBA-2 instrument at 850 μm and with at least partial Herschel coverage. In the following, we discuss the results obtained at large scales, based on IRAS, Planck, and Herschel data, and at small scales, based on the combination of Herschel and SCUBA-2 observations.

5.1 Dust SEDs of PGCC host clouds

PGCC sources represent the coldest part of the interstellar medium (ISM), and the environments examined in Sect. 4.1.2 should therefore consist mostly of parts of molecular clouds. This is confirmed by the column densities. Based on IRAS and Planck data at 100–550 μm and the dust opacity adopted in Sect. 3.3, the average values are N (H ₂) = 3.5 × 10²¹ cm⁻² within the r = 8′ annuli, but only 1.1 × 10²¹ cm⁻² immediatelyoutside this region. Typical SEDs measured with IRAS and Planck data correspond to a colour temperature of T =15.7 K and a dust opacity spectral index of β = 1.7, with some variations depending on the actual bands used (see Fig. 1). Most colour temperature estimates are within 1 K of the quoted median values.

When the IRAS 100 μm point is replaced with 250–550 μm Herschel measurements, the median colour temperature is about 14.7 K and the spectral index is close to β =1.8 (see Fig. 3a). Bootstrapping gives purely statistical uncertainties of 0.13 K in temperature and 0.015 in β (sample median values). The difference between the fits could be affected by very small grain emission in the 100 μm band or by the fact that wavelengths below the emission peak are more sensitive to warm dust and lead to lower β estimates when the temperature varies. When the analysis was repeated with background-subtracted flux measurements (Fig. 3b), the spectral index estimate remained similar, which could support the latter explanation.However, we also repeated the fits using only Planck data. In spite of the absence of short wavelengths λ < 350μm, these resulted in lower spectral index values β ~ 1.6 (see Appendix D). This result is, of course, affected by larger uncertainties. We note that in Malinen et al. (2014), the use of different wavelength ranges (from 100–500 μm to 160–1380 μm) did not result in significant changes in the (T, β) estimates of the proto-stellar and starless cores of the cloud LDN 1642 either.

The colour temperatures of the selected fields are lower and the spectral indices higher than in the general ISM. For example, Planck Collaboration Int. XIV (2014) reported in the 100–850 μm range values β = 1.54 for regions dominated by atomic gas and β = 1.66 for mostly molecular clouds (see also Planck Collaboration X 2016). That our SEDs are somewhat steeper (with higher β) reflects the fact that the selected fields represent the densest parts of molecular clouds. Juvela et al. (2015a) found T = 16.1 K and β = 1.84 for a sample of PGCC field, also using the combination of IRAS and Planck bands at 100–850 μm. The values are comparable to those derived for the present sample, although the median spectral index in Fig. 1 is slightly lower, β = 1.73. The values also depend on the relative weights given for the fitted bands, which are not identical between these two studies. In the present study, the weight of the 857 GHz and 545 GHz relative to both IRAS and lower frequency Planck bands was higher than in Juvela et al. (2015a).

The results of Fig. 3a and b could be affected by temperature variations within the apertures. By removing the warmer diffuse component, background subtraction should make the remaining emission more isothermal, which could in turn increase the observed β values. However, no such effect was observed.

5.2 SEDs of SCUBA-2 clumps

Dust emission of dense clumps and cores was examined with Herschel and SCUBA-2 data. Depending on the P(χ²) value, the average clump area of the samples is 2.1–2.5 arcmin² or 0.19–0.3 pc². The spectral index distributions were examined in Fig. 4 for alternative reductions, clump selections, and wavelengthranges. These comparisons help to establish the full uncertainty of the results that may not be captured by the a formal photometric and calibration errors alone. In spite of possible differences in the reliability (partially reflected in the width of the distributions), the main results are relatively robust. Most estimates of the submillimetre spectral index are in the range β =1.6 − 1.9. The 850 μm data point has a clear effect on the spectral index estimates. While the fits to the SPIRE data gave median β values slightly above 1.8, the joint fits of SPIRE and SCUBA-2 850 μm decreased these to β ~ 1.6. However, as demonstrated by Fig. 4l, there can be error sources that are not yet all fully understood.

The SPIRE fits can be compared with the results of Juvela et al. (2015a). The spectral index values there were higher at a scale of 10′. The median estimates were T = 14.9 K and β = 2.04 compared to the values T = 14.7 K and β = 1.8 in Fig. 3a. Juvela et al. (2015a) did not report statistics on dense clumps, but presented T and β maps, which sometimes showed clear spatial β variations with maximum values rising above β = 2.0. This was true mostly for well-resolved nearby clumps. The median distance of the fields in Juvela et al. (2015a) was 620 pc. In this paper the median distance of the fields is only 230 pc, but, as shown by Fig. 9, most of the detected clumps are at larger distances. The nearby fields tend to be at higher Galactic latitudes, have lower column densities, and therefore have fewer SCUBA-2 detections. Based on Juvela et al. (2015a), we can expect significant dispersion in the β histograms. The β estimates are in the range 1.0 < β < 2.5, and we found no real evidence for extreme values of either β ≤ 1 or β ≥ 3.

Data at wavelengths below 250 μm were mostly excluded from the analysis because of the incomplete coverage of the SCUBA-2 fields and because of the potential effects of temperature mixing. In Sect. 4.3 the addition of the 160 μm point changed the median values by ΔT = +0.4 K and Δ β = −0.1. The exact values again depend significantly on the relative weighting of the different bands (see e.g. Shetty et al. 2009b).

The β values derived in this work can be compared with recent results on starless and proto-stellar clumps. Sadavoy et al. (2013) analysed Herschel and SCUBA-2 observations of Perseus (clump B1) and found typical spectral index values of β ≈ 2. Towards proto-stellar sources, the values were lower, but still mostly above β = 1.6. The values are thus similar to this study, although we did not observe a significant difference between clumps with and without YSO candidates that could be interpreted as signs of systematic dust evolution (Fig. 7). Chen et al. (2016) investigated the Perseus cloud further and found notable differences between the various cloud regions. The pixel-by-pixel histograms covered the range β ~ 1.0 − 2.7. The lower β values were found to be correlated with local peaks of colour temperature, but also with CO outflows and thus directly with YSOs. This still raises the question to what extent the changes in the observed β values are the result of dust evolution or of the more complex temperature structure of the sources. Bracco et al. (2017) mapped spectral indices of the Taurus B213 filament, including both pre-stellar and proto-stellar clumps. The observations consisted of 1.15 mm and 2.0 mm ground-based observations that were analysed together with temperature constraints obtained from Herschel. A pre-stellar core was found to have a high spectral index of β ~ 2.4 that also remained radially constant. Because the values resulted from a 3D inversion, they are already corrected for the line-of-sight temperature mixing and should thus be somewhat higher than values dirctly derived from the SEDs. In the case of proto-stellar sources, the inversion of the radial profiles suggested a decrease from β ~1.5 − 2.0 in the outerparts to values 0.5 lower in the warmer centre. The lowest values at a core centre were β ~1.0, but the values derived from the basic SED analysis would be expected to remain significantly higher. Li et al. (2017) specifically examined a sample of Class 0 YSOs with interferometric observations from submillimetre to centimetre wavelengths. Most sources exhibited values β < 1.7 (down to β ~ 0.6). The authors noted that in addition to grain growth, the result may be affected by a strong temperature mixing when the high-resolution observations probe the emission from the circumstellar disk. The comparison with the Herschel and SCUBA-2 results is also affected by the difference of the wavelength ranges (if submillimetre SEDs are generally steeper) and the interferometric observations filtering out much of the extended emission.

In summary, our results of T ~ 14 K and β ~1.6 for the joint SPIRE and SCUBA-2 fits are very similar to the previous studies. However, in our sample, the difference between starless sources and proto-stellar sources (clumps with YSO candidates) is small. Our YSO statistics may also be incomplete and affected by the sensitivity and resolution of the WISE survey relative to the data used in the other studies.

In our study, data reductions with different Herschel-based source masks and different filtering scales in Fig. 4a–i gave relatively consistent results. Larger filter scales and larger source masks resulted in the recovery of more 850 μm extended emission. The maps generally agreed well with the SPIRE observations, which are sensitive to much lower column densities. The use of θ_F = 500′′ and large LM masks leads to more than twice as many clump detections as the combination of θ_F = 200′′ and SM. On the other hand, large masks can sometimes increase random fluctuations in the maps. For example, in Fig. C.1, the main clump of G167.215.3A1 is associated with a negative feature that is located inside the LM mask. The feature disappears when a smaller filter scale (200′′ instead of 500′′) or a tighter source mask (SM instead of LM) is used. While a large filter scale and liberal source masks work well in statistical studies, one may wish to be more conservative when studying individual objects.

5.3 Wavelength dependence of the spectral index

One of the main goals of this paper was to investigate the potential flattening of dust SEDs towards millimetre wavelengths. In this respect, the most interesting comparisons are the fits with and without the Planck 217 GHz band (1380 μm), and at smaller scales, the fits with and without the SCUBA-2 850 μm band.

According to the IRAS and Planck analysis in Fig. 1, the 100–1380 μm fit leads to a marginally lower value of the spectral index (β = 1.69) than the fits at λ ≤ 550 μm (β = 1.74) or λ ≤ 850 μm (β = 1.73). The corrections for line emission in the 217 GHz and 353 GHz bands was identified as the main source of uncertainty. Based on Fig. 1, such small differences in β could be explained by a systematic 10% error in the CO corrections. The comparison with the sum of the ground-based ¹² CO(2-1) and ¹³ CO(2-1) measurements showed that the Planck CO estimates, which were used to correct the continuum data, are correct to within ~30% and do not show any significant bias (Fig. 2a). However, the possibility of systematic errors at 10% level cannot be excluded. The calibration accuracy of ground-based observations is also not much better than this. Considering the uncertainty of the line contamination, which may include contributions of other lines as well (Planck Collaboration X 2016), and because of other potential systematic errors, the change of β as a function of wavelength was not significant.

When the same analysis was based on the combination of Herschel and Planck data, the 217 GHz band caused no significant change in β (Fig. 3a). Furthermore, in the background-subtracted measurements, the β value was higher when the 217 GHz band was included. This shows that if β does change with wavelength, the effect is at most at the level of the observational uncertainties.

The analysis of Herschel and SCUBA-2 data gave clear indications that the 850 μm point would favour lowerβ values (Figs. 4 and 5). The difference with and without the 850 μm point resulted in a change of about Δβ ~0.2, but the result may still depend on the details of the data reduction, as indicated by Fig. 4j–l. On the other hand, a relative calibration error of 10% between SPIRE and SCUBA-2 (or similar CO correction of the 850 μm band) would correspond to a change of only Δβ = ±0.03 in our fits.

5.4 Relation between T and β

The dependence between dust temperature and spectral index is important because it is related to dust evolution during the star formation process and is a potential tracer of the past history and current physical state of clumps. The simultaneous determination of T and β is difficult because observational noise produces anticorrelation between the parameters. Moreover, temperature variations within the beam make the observed colour temperature a biased estimate of the physical grain temperature, and spectral index estimates will be similarly biased.

It is therefore not surprising that the observed relations in Figs. 1, 3, and 6 all exhibit some negative correlation between T and β. In the following, we examine whether these can be caused by noise alone. It is more difficult to make the final connection between the apparent and the real dust properties. This would require detailed modelling and is beyond the scope of this paper.

We first investigate the SEDs of IRAS and Planck data shown in Fig. 1a. Fits of a function β(T) = A × T^B (see Désert et al. 2008; Juvela et al. 2015a) give for B values of −0.53 and −0.69 for the 100–550 μm, and 100–1380 μm data, respectively. We examined whether the data are consistent with a constant β value, the negative correlation being caused by noise alone. We fitted the data of each field with an MBB function B_ν (T)×ν^1.7 and rescaled the flux values onto this SED curve. In other words, at this point, the simulated photometry exactly matches the assumed β(T) law. We then added photometric noise to the flux density values (see below) and determined the apparent T and β values with new fits of the simulated SEDs. In the above we assumed B_ν(T) × ν^1.7 spectra, but because the goal is to simulate the relative difference between the true and the observed SEDs, the procedure is not sensitive to the selected value of β. The results remain essentially identical when a value of β = 1.4 or β = 2.0 is used instead of β = 1.7

Because a direct comparison between observations and simulations would require precise knowledge of the uncertainties, we included the noise level as an additional parameter k (noise). Thus we used the same error covariance matrices as in Sect. 4.1.1, but scaled all the matrix elements with the square of k (noise). We examined the effect of k (noise) on the β(T) relation, and on the other hand, the scatter around the fitted β = A × T^B curve. Comparison with the observed scatter should provide the most direct empirical upper limit on the noise effects. It is an upper limit because the real emission cannot be expected to follow a single β(T) relation, and thus only part of the observed scatter is caused by noise.

For the 100–550 μm fits (Fig. 10a, c), the simulations match the observed scatter over a range of k (noise) values, where the value of the B parameter of the β(T) relation is also similar to the observed value. The data are thus compatible with a constant spectral index. The situation is different for the fits to 100–1380 μm data, which have a smaller uncertainty on the β parameter (Fig. 10b, d). The observed scatter is matched near k(noise) = 1.0, at which point B is –0.35 instead of the value B = −0.69. The result would thus suggest that the intrinsic β(T) relation is indeed a decreasing function of temperature. Half of the apparent decrease would be caused by noise, but only assuming that the intrinsic β(T) relation itself has no scatter. The result remains qualitatively the same, although with different k (noise) values, if the relative weighting of the bands is modified (such as doubling the error estimates of the IRAS band or of the two longest wavelength Planck bands). Figure 10 contains some inconsistency because the 100–550 μm fits suggest a lower value of k (noise). However, if the relative error estimates of the 353 GHz and 217 GHz bands are increased, it is possible to fit both 100–550 μm and 100–1380 μm ranges even with the same value of k (noise).

In Fig. 11 we show the same test applied to the combination of Herschel and Planck data (cf. Fig. 1a). In the simulations, the noise was generated using the same error covariance matrices as in Sec. 4.1.2, but scaling the elements with k (noise) squared. For this shorter wavelength baseline of 250–850 μm, the observed scatter is reproduced around k(noise) ≈ 0.45, at which point the predicted B is only slightly above its observed value. For the longer wavelength baseline 250–1380 μm, the scatter and observed B can both be reproduced with k(noise) = 0.8, and the data are consistent with a flat intrinsic β(T) relation. This does not completely exclude a decreasing β(T) function if the clumps do have some real scatter relative to the average β(T) relation.

For the clumps, we can only use the shorter wavelength baseline provided by Herschel and SCUBA-2 data. The error estimates that were used to derive the confidence regions in Fig. 6 may have been too large because the 67% contour already covers almost all of the observations. Nevertheless, the figure shows that the observed distribution of points is consistent with the orientation of the error ellipse, and the observed distribution of (T, β) points might in principle be caused by noise.

Figures 10 and 11 only weakly support the T −β anticorrelation, although by attributing all the observed (T, β) scatter to observational errors, we underestimate the probability of a decreasing β(T) relation. The results at least exclude the possibility of an increasing (apparent) β(T) relation.

If the intrinsic β of dust grains is constant, the apparent values (derived from the SEDs) could still appear to be positively correlated. This is the case for samples of externally heated clumps of different masses (see e.g. Juvela & Ysard 2012b, Fig. 6). If the clumps have temperaturevariations caused by proto-stellar sources, they tend to cause a negative correlation between the observed (apparent) T and β values (Shetty et al. 2009b; Malinen et al. 2011; Juvela & Ysard 2012b). When the noise effects are excluded, this would still not constitute proof of similar behaviour for the actual dust opacity spectral index of the dust grains if we were to conclude that the observation exhibits a negative correlation between the apparent T and β values. Figure 10b gave the strongest indications of an apparent T −β anticorrelation, but because the 100 μm data were included, the result may be affected by temperature variations within the aperture or more directly by the emission of very small grains.

The effect of noise on the β(T) relation and the joint estimation of the noise level and the β(T) relation are examined further in Appendix G.

Fig. 10 Comparison of (T, β) relations in observations and in simulations with a constant value of β = 1.7. The observations consist of the same IRAS and Planck data as in Fig. 1 and represent the average emission of the fields. Left and right frames show results for 100–550 μm and 100–1380 μm fits, respectively. In the upper frames, black symbols show the real observations, and green crosses show one realisation of simulations with k (noise) equal to 0.4 (left frame) and 1.0 (right frame). The black lines are fits β = A × TB to the observations (solid line) and to the green crosses (dashed line). The lower frames show the change in the rms value of the residuals of the β(T) fit (solid blue lines, left axis) and the value of B (dashed red lines and right axis) as a function of k (noise). The horizontal lines correspond to the values in the actual observations.

Fig. 11 Same as Fig. 10, but for the combined Herschel and Planck data from Fig. 3. The fitted data are 250–850 μm (frames a and c) and 250–1380 μm (frames b and d). The green crosses correspond to one realisation with k (noise) = 0.4 (frame a) and k (noise) = 0.7(frame b).

6 Conclusions

We have used IRAS, Planck, Herschel, and SCUBA-2 observations to examine the dust emission spectra in 96 fields that were originally selected based on the presence of prominent PGCC cold clumps. SCUBA-2 data, in conjunction with Herschel SPIRE observations, were used to examine the submillimetre spectra of the densest clumps and cores. The study led to the following conclusions:

• The fields have typical colour temperatures of 14–18 K and spectral index values of β = 1.5–1.9.

• The clumps extracted from SCUBA-2 maps are characterised by median values of T ~13 K and β ~ 1.7. Different options in the data reduction lead to uncertainty at a level of δβ ~±0.2.

• The use of large source masks and large filter scales was successful in recovering more extended emission in the 850 μm SCUBA-2 maps. This also resulted in a larger number of clump detections. Because very extended masks can lead to undesirable map fluctuations in individual maps, they may be most useful in statistical studies.

• At large scales, there is no consistent proof for the dust spectrum becoming flatter at millimetre wavelengths. We are limited by the uncertainty of the corrections for CO line emission in the 353 GHz and 217 GHz Planck bands.

• At clump scales, the SCUBA-2 850 μm flux measurements tend to be above the SED fitted to the shorter wavelength SPIRE bands. The effect on the β values in the joint fits is Δβ ~ 0.2. The result is not entirely robust with respect to the decisions made in the data reduction.

• The positions of the SCUBA-2 submillimetre clumps are positively correlated with Class I and II YSO candidates. For Class III sources, the correlation is negative.

• Compared to starless clumps, protostellar clumps have a 1–2 K higher median colour temperature and marginally lower β values. The two populations mostly overlap, however. This also applies to the mass and volume density distributions.

• Most of the observed T-β anticorrelation may be caused by noise. Evidence for the anticorrelation is found mainly at large scales.

Acknowledgements

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of PR China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. This research made use of Montage, funded by the National Aeronautics and Space Administration’s Earth Science Technology Office, Computational Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. The code is maintained by the NASA/IPAC Infrared Science Archive. The data presented in this paper are partly based on the ESO-ARO programme ID 196.C-0999(A). M. Juvela and V.-M. Pelkonen acknowledge the support of the Academy of Finland Grant No. 285769. T. Liu is supported by KASI fellowship and EACOA fellowship. J. Malinen acknowledges the support of ERC-2015-STG No. 679852 RADFEEDBACK. V.-M. Pelkonen acknowledges the financial support from the European Research Council, Advanced Grant No. 320773 entitled Scattering and Absorption of Electromagnetic Waves in Particulate Media (SAEMPL).

Appendix A: Colour corrections

The maps were colour-corrected pixel by pixel, using the MBB fits B(T) × ν^1.8 to the original, unfiltered SPIRE data. Because the fits that include the 850 μm point can lead to different SEDs, we examined how the colour correction factors vary over the plausible parameter range. Figure A.1 compares the ratios of colour correction factors⁸ relative to these factors at 350 μm. The normalisation with the 350 μm values is included because we are interested only in the SED shape, not its normalisation. In the plot, the ratio is further normalised so that it is one for T = 15.0 K and β = 1.5. If the ratios significantly change between the position where colour correction was estimated and the true SED shape, it would bias the fit results.

According to Fig. A.1, the relative variations are of about 2% or smaller over most of the plausible parameter space. Since the differences between the fitted and the unknown true SED shape should be smaller, and if the differences are caused by noise, they are mostly along the line between the high T and low β and the low T and high β combinations. Thus, the colour correction uncertainties are likely to be at most ~2% and thus remain much smaller than the typical photometric errors, which are closer to 10%. A 2% relative error between 350 μm and 850 μm corresponds to a spectral index error smaller than Δβ = 0.05.

Fig. A.1 Colour correction factors at the 250 μm, 500 μm, and 850 μm bands relative to the same at 350 μm. The ratios are normalised to one at T = 15 K and β = 1.5. The contours denote the ratio one (solid lines) and deviations from one by ± 1% (dashed lines) and ±2% (dotted lines).

Appendix B: CO contamination in the 850 μm data

We do not have high-resolution observations of the CO J = 3−2 lines to directly correct the SCUBA-2 850 μm data for linecontamination, but we can use lower transition data to estimate the possible effect.

The ¹²CO(1-0) maps of Meng et al. (2013) and Zhang et al. (2016) cover a number of SCUBA-2 targets. The CO observations were made with the Purple Mountain Observatory telescope, which at 115 GHz has a beam size of 52′′. This is sufficiently close to 40′′, the resolution of the maps where flux densities were measured, so that a direct comparison is meaningful. The ratio between the J = 3−2 and J = 1−0 lines is not known. The line ratios were discussed in Juvela et al. (2015a), where they were found to be mostly ~0.3 or below, both based on Planck measurements of the velocity-integrated CO emission and based on examples of ground-based measurements. The low line ratios are explained by the nature of our sources, PGCC being a selection of the coldest dust emission sources on the sky. This does not mean that the line ratios could not be much higher in individual, actively star-formingclumps. In the following we assume a line ratio of one, which thus very likely overestimates the typical CO contamination. For example, Fig. B.1 suggested that the line ratios (2–1)/(1–0) would be lower than 0.5.

We scaled the CO emission to 850 μm flux density with a constant 0.70 mJy beam⁻¹ per K km s⁻¹ (Drabek et al. 2012), and similar to continuum photometry, calculated the difference between clump apertures and reference annuli. The exercise was made using the LM data version. We did not process the CO maps through the SCUBA pipeline. This should lead to an overestimation of the CO effects and not to an underestimation because we ignored the loss of CO intensity that is caused by the spatial filtering.

Figure B.1 shows the ratios between CO-corrected and original 850 μm flux densities. The mean value is close to one, and for many clumps, the CO correction would increase the 850 μm flux density estimates. This means that the CO correction was larger outside the clump, which is not expected given that the cold clumps are not strong CO emitters. Molecular depletion and the anticorrelation with local heating sources can both lead to a situation where the CO emission is stronger outside the cold clump.

Fig. B.1 Ratio of CO-corrected and original 850 μm flux values. The values are for a subset of clumps from the P(χ2) 85%, 25%, 10%, and 5% samples (white, blue, yellow, and red, respectively) with the LM data version. The mean and standard deviation values are listed in the figure.

Figure B.1 suggests that even assuming a line ratio (3–2)/(1–0) of one, the effect of the CO contamination on the 850 μm flux density would typically be about two percent. Even if the CO correction were a systematic − 5% for most clumps, the net effect on themedian β of the 250–850 μm fits would be only Δβ = +0.10.

Appendix C: Plots of Herschel and SCUBA-2 maps

Figure C.1 shows maps of 250 μm optical depth and 850 μm surface brightness for all the 96 fields.

Fig. C.1 Maps of 250 μm optical depth and 850 μm surface brightness. In τ(250μm) maps, white and black contours show the extent of the LM and M850 source masks, respectively. Dashed ellipses show the FWHM size of PGCC clumps (numbers at the PGCC centre coordinates refer to the numbering in Table H.3). The stars stand for the Class I-II (cyan stars) and Class III (yellow stars) YSO candidates (Marton et al. 2016). In the S(850 μm) maps, the clumps of the P(χ2) = 85% sample are indicated with white contours and the other clumps with grey contours.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

Fig. C.1 continued.

The optical depths τ(250 μm) are based on MBB fits of the SPIRE 250 μm, 350 μm, and 500 μm data, without background subtraction. The intensity zero-points of the input maps (obtained from HSA) were set through a comparison with Planck data. For optical depths, our SED fits use a constant value of β =1.8 and a map resolution of40′′. The SCUBA-2 850 μm maps correspond the LM version of data. (see Sect. 2.1). These are shown at the full resolution of 14′′ and are scaled to surface brightness units MJy sr⁻¹.

The contours in the τ(250 μm) maps indicate the extent of the LM and M850 source masks that were used in the SCUBA-2 pipeline. The clumps extracted with the Fellwalker method are indicated in the 850 μm maps with white contours. Only the sources with dashed black contours were used in investigations of the spectral index values.

Appendix D: SED fits of Planck data alone

In Sect. 4.1.1 the mean emission of the fields was fitted with IRAS and Planck data. The shortest Planck wavelength of 350 μm is near the maximum of the cold dust emission spectrum. Nevertheless, we can attempt SED fits also without the IRAS data point (Fig. D.1).

Fig. D.1 MBB fits to Planck data averaged within 8′ of the centre of the SCUBA-2 fields. The wavelength range is either 350–850 μm or 350–1380 μm. Each marker corresponds to one field, and the dotted lines connect estimates of the same field. The median values are plotted with large symbols, and they are quoted in the frame together with the 1σ dispersion estimated from the interquartile ranges. The plot assumes the default CO corrections.

Although the scatter is significantly larger than in Fig. 1, the median values are similar to the fits where the IRAS 100 μm point was also used. There is no evidence that the latter would be significantly affected by very small grain emission or by temperature variations inside the measurement aperture. As in Fig. 1a, the extension of the fitted wavelength range to 1380 μm slightly decreases the median β value.

Appendix E: Dependence between dust parameters and error estimates

The dependence of β estimated from SPIRE and SCUBA-2 data is compared to some other parameters in Fig. E.1. The estimates are almost independent of the (non-reduced) χ² values of the fits and are weakly correlated with the colour temperature that is estimated from 250–500 μm data using a fixed value of β = 1.8. Finally, Fig. E.1c shows some negative correlation with the source brightness. However, none of the correlations are statistically significant.

Fig. E.1 Dependence of β estimates on χ2 values (frame a), SPIRE-derived colour temperatures for β = 1.8 (frame b), and 250 μm flux densities (frame c). The colours are the same as in Fig. 4. The linear correlation coefficients r and median-β curves are shown for the full set of plotted points.

It would seem natural to analyse β distributions using the sources with the lowest β error estimates. However, the (T, β) anticorrelation leads to a narrow, curved χ² valley in the (T, β) plane. In the case of large uncertainties, it runs at low T runs parallel to the β axis and at high T parallel to temperature axis. Thus, the selection of sources with a low β uncertainty would give preference to warm sources with low β values. The effect is illustrated in Fig. E.2, where the error estimates δ(T) and δ(β) correspond to half of the interquartile range estimated from MCMC runs. As expected, the figure shows a correlation between β and the error estimate of the β uncertainty. Subsamples selected based on the χ² values of the fits are not significantly correlated with δ(T) or δ(β).

Fig. E.2 Estimates of β as the function of the T and β MCMC error estimates. The symbol colours are as in Fig. E.1, and the black lines show the median values as a functionof the x-axis variable. The linear correlation coefficients r are given in the figure.

Appendix F: Effect of map making and clump selection on β

The obtained β distributions depend on our decisions regarding the map making and the data calibration. In map making, the options were to use 200′′ or 500′′ filtering scale. Furthermore, maps could be produced using either large or small source masks based on Herschel 500 μm data (SM and LM, respectively) or masks based on SCUBA-2 850 μm data themselves (M850). For illustration of the calibration uncertainty of the 850 μm data, we compared the result for the flux conversion factor FCF = 2.41 Jy (pW)⁻¹, estimated from calibration measurements near the observing times, to the results obtained with the default calibration value of FCF = 2.34 Jy (pW)⁻¹ arcsec⁻².

The following table describes the β distributions. These include the β distributions calculated using SPIRE data alone (unfiltered and filtered versions) because the clump selection also changes from case to case.

The data show that the two cases of the 850 μm calibration differ by less than Δβ = 0.1, except for the M850 case, where the larger FCF factor results in a median spectral index value lower by 0.2 units.

Appendix G: Effect of noise on the β(T) relation

In Sect. 5.4 we used simulations to test the significance of the observed negative correlation between the colour temperature and the spectral index. The noise was kept as a free parameter because otherwise the result would critically depend on the precision of the a priori estimates of photometric errors. Here we examine the method further with purely synthetic observations.

We assumed an intrinsic dependence of $\beta (T)={(T/15\text{\hspace{0.17em}K})}^B$ and simulated observations at 250 μm, 350 μm, 500 μm, and 850 μm. The default photometric errors were 4% for the first three and 10% for the 850 μm band. We generated a sample of 100 clumps. The temperatures followed a normal distribution with a mean of 15 K and a standard deviation of 3 K. Together with the β(T) relation, this defines the source fluxes, to which we added normally distributed photometric errors. The fitting of a set of synthetic observations gives the observed values of B and the rms dispersion with respect to the fitted β(T) curve. The analysis of these synthetic observations proceeded as described in Sect. 5.4. We generated a similar sample of observations under the assumption of a flat β(T) relation. This was obtained by fitting the syntetic observations with a fixed value of β and replacing the flux density values by the values predicted by the fit. We then examined the B parameter and the rms dispersion of β vs. the fitted β(T) relation of this simulated sample as a function of the noise level. The noise level was again scaled by the parameter k(noise). When the rms noise in the simulation matched the rms noise of the (this time synthetic) observations, we checked whether the corresponding value of B in the simulation was higher than the value of B in the observations. When the difference was positive and significant, we concluded that the observations exhibited a non-constant β(T) relation (β decreasing with temperature) that could not be explained by the noise.

Figure G.1 shows the results for two values of B and k (noise). The simulations were made for 100 values between k (noise) = 0.1and k (noise) = 2.0, which thus also indicate the dispersion between different realisations. The left-hand frames correspond to cases where the observational noise is half of the default. The observed rms dispersion of β values is matched at a point where the simulated B is still significantly above the observed value. The test thus shows that the T −β anticorrelation is real in the (synthetic) observations. Comparison of frames a and c shows that the significance of the detection (of non-constant β(T)) decreases as the intrinsic relation becomes flatter.

The frames b and d correspond to cases where the observational noise is 50% higher than the default values. In frame d, with B = −0.2, the simulations match the observed values at exactly the same k (noise) value, which means that the data are compatible with a flat β(T) model. Only when the intrinsic relation is steeper with B = −0.6 is the detection of the T − β anticorrelation still found to be statistically significant (frame b). In these cases, we could proceed with further simulations to quantify the B value of the intrinsic β(T) relation. However, plots like Fig. G.1 already indicate which fraction of B is caused by noise.

Fig. G.1 Analysis of a set of synthetic observations that follow an intrinsic β(T) ∝TB relation. The parameters B and k (noise) of the input simulation are given in the frames and are indicated by solid vertical and horizontal lines. The horizontal dot-dashed and dashed lines correspond to the values of B and the β rms dispersion, respectively, which are estimated from one realisation of synthetic observations. The other curves show B (dot-dashedred line) and the rms dispersion (blue dashed line) from simulations as a function of the noise scaling k (noise).

The recovered k (noise) value (when the rms value of β values vs. the fitted β(T) relation matches the rms value in the observations) is a biased estimator of the photometric errors of the observations. For example, in Fig. G.1b the difference is a factor of 1.5 between the input value k (noise) = 1.5 and the fitted value k(noise) ≈1.0. The effect is smaller for low absolute values of B because the simulations explicitly start with the B = 0 assumption.

Figure G.1 suggests that the precision of SPIRE and SCUBA-2 observations should be sufficient to make a detection of the anticorrelation, assuming that the default error estimates are correct and the value of B is −0.2 or lower. However, the simulations are highly idealised. They do not take the possibility of systematic or non-Gaussian errors or the effects of temperature mixing into account. One of the assumptions is that all sources would follow the same β(T) relation. In Fig. G.2 we repeat the previous analysis, but allow a dispersion of σ(β) = 0.1 with respect to the β(T) relation in both the synthetic observations and in the subsequent simulations. This increases the noise in the plotted relations, but does not significantly affect the main features of Fig. G.1.

Fig. G.2 Same as Fig. G.1, but allowing for a dispersion σ(β) = 0.1 of the intrinsic β(T) relation.

Appendix H: Additional tables

References

All Tables

All Figures

Fig. 1 Results of MBB fits to IRAS and Planck data averaged within 8′ of the centre of the SCUBA-2 fields. Frame a shows the fits completed with 100–550 μm, 100–850 μm, and 100–1380 μm bands. Each marker corresponds to one field, and the dotted lines connect estimates of the same field. The median values calculated over all the fields are plotted with large symbols. The median values and the 1σ dispersion estimated from the interquartile ranges are quoted in the frame. Frame b shows the fits to the 100–1380 μm data without CO corrections (ΔCO = 0), with the default CO correction (ΔCO = 1), and with twice the default correction (ΔCO = 2). The median values are indicated by the three large triangles. The contours show the 67% (solid line) and 95% (dashed line) confidence regions for a 100–850 μm (frame a) and 100–1380 μm fit with ΔCO = 1 (frame b) calculated using the median error covariance matrix of the fields.

In the text

Fig. 2 Comparison of CO estimates from Planck Type 3 CO maps and ground-based observations. Upper frame compares the sum of 12CO(2–1) and 13 CO(2–1) values of 56 SMT maps to Planck CO estimates for the assumed line ratio of (2−1)/(1−0) = 0.5. Lower frame shows comparison with CO(1–0) observations from the PMO telescope, with data at 7′ resolution sampled at 3′ steps. Each colour and symbol combination corresponds to one of 17 fields. Solid lines show the one-to-one relations, and dotted lines correspond to slopes different by ±30%.

In the text

Fig. 3 Results of MBB fits to Herschel and Planck data averaged within 8′ of the centre of the SCUBA-2 fields. The results are shown for absolute surface brightness measurements (frame a) and after subtraction of the local background emission (frame b). The symbols correspond to the band combinations listed in the legend in frame b. The large symbols denote the median values over the fields. In frame a, the contours show the 67% (solid lines) and 95% (dashed lines) confidence regions for 250–500 μm (red) and 250–1380 μm (blue) fits, using error covariance matrices that are the median over all fields. For the SPIRE-only fits, the adopted error estimates clearly overestimate the true uncertainty.

In the text

Fig. 4 Distributions of β values in fits to Herschel and SCUBA-2 850 μm observations of clumps. Results are shown for original SPIRE data (first row), filtered SPIRE data (second row), and the combination of filtered SPIRE and 850 μm data (third row). The columns correspond to different versions of data reduction, as indicated above the frames. The white, blue, yellow, and red histograms correspond to P(χ2) values of 85%, 25%, 10%, and 5%, respectively (see text). The number of clumps in each sample and the quartile values (25%, 50%, and 75% percentiles) of the blue and red histograms are given in the frames.

In the text

Fig. 5 SEDs of eight sources from the P(χ2) = 5% subsampleof SM clumps in Fig. 4f. Plots include fits to the original SPIRE data (upper red curve and red numbers) and to the filtered SPIRE data plus the SCUBA-2 850 μm point. The clump numbers given after the field name are the same as in Fig. C.1. The 160 μm data points are available for one of the clumps: it is plotted in the figure, but is not part of the fit.

In the text

Fig. 6 Spectral index β vs. colour temperature T for the clump samples of Fig. 4a–c (left frame) and Fig. 4d–f (right frame). The black symbols correspond to fits to the original SPIRE data, and the blue symbols to the combination of SCUBA-2 and filtered SPIRE data. Only clumps from the P(χ2 )=25% samples are included. The median values are indicated with large crosses. Frame a includes a contour for a typical 67% confidence region in the four-band fits.

In the text

Fig. 7 Comparison of Herschel MBB fits with and without the 160 μm data point inthe case of the P(χ2) = 25% clump sample with θF = 500′′ and LM masks. The left frame shows the distribution of ratios between the 160 μm flux density and the value predicted by the 250–500 μm fit. The greyhistogram contains all clumps, and the black histogram is based on three clumps with YSO candidates. The second frame shows (T, β) values for individual clumps for the 250–500 μm (blue triangles) and 160–500 μm (red squares) fits. The circles indicate clumps with YSO candidates. The median values are quoted in the frame and are plotted withlarger symbols.

In the text

Fig. 8 YSO statistics. The surface density of sources is shown as a function of the distance ( >2′) from the field centre. The mean values over all fields are shown for Class I-II (thick red line) and Class III (thick yellow line) candidates. The thin red lines correspond to Class I-II sources in individual fields. The shading corresponds to a region with less than one source for a given distance. The two overlapping stars indicate the density of Class I-II candidates in all clumps (all fields) selected by the Fellwalker algorithm (red star) and in clumps belonging to the P(χ2 ) = 85% sample (blue star).

In the text

Fig. 9 Clump mass vs. field distance (frame a) and average volume density vs. temperature (frame b). The plots correspond to the LM data version and the P(χ2 ) = 25% clump sample. The fluxes and temperatures are from the 250–850 μm fits. The data are further divided into three samples based on temperature (frame a) or mass (frame b). The clumps with YSO Class I or II candidates are marked with circles. Frame a shows the thresholds corresponding to a 17 mJy per 850 μm beam (dashed line) and to a 850 μm flux density of 0.1 Jy f(dotted line; see text).

In the text

Fig. 10 Comparison of (T, β) relations in observations and in simulations with a constant value of β = 1.7. The observations consist of the same IRAS and Planck data as in Fig. 1 and represent the average emission of the fields. Left and right frames show results for 100–550 μm and 100–1380 μm fits, respectively. In the upper frames, black symbols show the real observations, and green crosses show one realisation of simulations with k (noise) equal to 0.4 (left frame) and 1.0 (right frame). The black lines are fits β = A × TB to the observations (solid line) and to the green crosses (dashed line). The lower frames show the change in the rms value of the residuals of the β(T) fit (solid blue lines, left axis) and the value of B (dashed red lines and right axis) as a function of k (noise). The horizontal lines correspond to the values in the actual observations.

In the text

	Fig. 11 Same as Fig. 10, but for the combined Herschel and Planck data from Fig. 3. The fitted data are 250–850 μm (frames a and c) and 250–1380 μm (frames b and d). The green crosses correspond to one realisation with k (noise) = 0.4 (frame a) and k (noise) = 0.7(frame b).
In the text

	Fig. A.1 Colour correction factors at the 250 μm, 500 μm, and 850 μm bands relative to the same at 350 μm. The ratios are normalised to one at T = 15 K and β = 1.5. The contours denote the ratio one (solid lines) and deviations from one by ± 1% (dashed lines) and ±2% (dotted lines).
In the text

	Fig. B.1 Ratio of CO-corrected and original 850 μm flux values. The values are for a subset of clumps from the P(χ2) 85%, 25%, 10%, and 5% samples (white, blue, yellow, and red, respectively) with the LM data version. The mean and standard deviation values are listed in the figure.
In the text

Fig. C.1 Maps of 250 μm optical depth and 850 μm surface brightness. In τ(250μm) maps, white and black contours show the extent of the LM and M850 source masks, respectively. Dashed ellipses show the FWHM size of PGCC clumps (numbers at the PGCC centre coordinates refer to the numbering in Table H.3). The stars stand for the Class I-II (cyan stars) and Class III (yellow stars) YSO candidates (Marton et al. 2016). In the S(850 μm) maps, the clumps of the P(χ2) = 85% sample are indicated with white contours and the other clumps with grey contours.

In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

	Fig. C.1 continued.
In the text

Fig. D.1 MBB fits to Planck data averaged within 8′ of the centre of the SCUBA-2 fields. The wavelength range is either 350–850 μm or 350–1380 μm. Each marker corresponds to one field, and the dotted lines connect estimates of the same field. The median values are plotted with large symbols, and they are quoted in the frame together with the 1σ dispersion estimated from the interquartile ranges. The plot assumes the default CO corrections.

In the text

	Fig. E.1 Dependence of β estimates on χ2 values (frame a), SPIRE-derived colour temperatures for β = 1.8 (frame b), and 250 μm flux densities (frame c). The colours are the same as in Fig. 4. The linear correlation coefficients r and median-β curves are shown for the full set of plotted points.
In the text

	Fig. E.2 Estimates of β as the function of the T and β MCMC error estimates. The symbol colours are as in Fig. E.1, and the black lines show the median values as a functionof the x-axis variable. The linear correlation coefficients r are given in the figure.
In the text

Fig. G.1 Analysis of a set of synthetic observations that follow an intrinsic β(T) ∝TB relation. The parameters B and k (noise) of the input simulation are given in the frames and are indicated by solid vertical and horizontal lines. The horizontal dot-dashed and dashed lines correspond to the values of B and the β rms dispersion, respectively, which are estimated from one realisation of synthetic observations. The other curves show B (dot-dashedred line) and the rms dispersion (blue dashed line) from simulations as a function of the noise scaling k (noise).

In the text

	Fig. G.2 Same as Fig. G.1, but allowing for a dispersion σ(β) = 0.1 of the intrinsic β(T) relation.
In the text