© ESO, 2015

1. Introduction

The all-sky survey of the Planck satellite (Tauber et al. 2010) has enabled a new approach to studying the earliest stages of star formation. The sub-millimetre measurements, with high sensitivity and an angular resolution down to ~ 4.5′, have enabled detecting and classifying of a large number of cold and compact Galactic sources. These probably represent different phases in the evolution of dense interstellar clouds that leads to the formation of stars. Careful analysis of the Planck data has led to a list of more than 10 000 objects that form the Cold Clump Catalogue of Planck Objects (C3PO, see Planck Collaboration XXIII 2011). At Planck resolution, it is not possible to resolve gravitationally bound cores even in the nearest molecular clouds. The low colour temperature of most of the sources (T< 14 K) strongly suggests that the Planck clumps must have high column densities, possibly at scales not resolved by Planck, and they probably contain even dense cores. A significant fraction of the clumps may be transient structures produced by turbulent flows, however.

Within the Herschel open time key programme Galactic Cold Cores, we have carried out dust continuum emission observations of 116 fields that were selected based on Planck detections listed in C3PO. The fields, which are typically ~40′ in size, were mapped with Herschel PACS and SPIRE instruments (Pilbratt et al. 2010; Poglitsch et al. 2010; Griffin et al. 2010) at wavelengths 100–500 μm. The higher angular resolution of Herschel (Poglitsch et al. 2010; Griffin et al. 2010) enables studying the internal structure of the Planck clumps, detecting individual cores, and, in conjunction with mid-infrared (MIR) data, studying protostellar sources (Montillaud 2015). The inclusion of far-infrared (FIR) wavelengths helps to determine the physical characteristics of the regions and, in particular, to study the properties of dust emission. First results have been presented in Planck Collaboration XXIII (2011), Planck Collaboration XXII (2011) and in Juvela et al. (2010; 2011; 2012; Papers I−III, respectively). Montillaud (2015) presented the analysis of all clumps and cores found in our Herschel fields, including a comparison with the population of young stellar objects (YSOs). Further studies concentrated especially on high latitude clouds (Malinen et al. 2014; Rivera-Ingraham et al. and Ristorcelli et al., in prep.).

In this paper we concentrate on dust properties and especially on the submillimetre dust opacity. Variations of dust emission properties have been investigated with FIR and submillimetre observations of diffuse and molecular clouds, using data from IRAS, COBE, ISOPHOT, the PRONAOS balloon-borne experiment, and ground-based telescopes. It was shown that the low temperatures found in a sample of molecular clouds (Laureijs et al. 1991; Abergel et al. 1994, 1996) and the translucent Polaris Flare cirrus cloud (Bernard et al. 1999) cannot be explained by the extinction of the radiation field. An increase of the dust emissivity by a factor of 3 compared to the standard diffuse value was needed to reproduce the cold temperatures observed in the Taurus filament L1506 (Stepnik et al. 2003). In dense regions, several studies have shown an opacity increase by a factor between 2 to 4 (Cambrésy et al. 2001; Kramer et al. 2003; Bianchi et al. 2003; del Burgo et al. 2003; Kiss et al. 2006; Ridderstad et al. 2006; Lehtinen et al. 2007). More recently, similar results have been obtained with Herschel and Planck in molecular clouds and cold cores (Juvela et al. 2011; Planck Collaboration XXV 2011; Martin et al. 2012; Roy et al. 2013). In their detailed modelling of Herschel observations of the L1506 filament, Ysard et al. (2013) characterised the dust evolution toward the dense part of the filament. The dust emissivity in the outer layers of the filament was found to be consistent with standard grains from the diffuse medium, whereas the emissivity increases by a factor of ~2 above gas densities of a few times 10³ cm^-3. This change has been attributed to the formation of fluffy aggregates in the dense medium, resulting from grain coagulation, as suggested by previous studies (Cambrésy et al. 2001; Stepnik et al. 2003; Bernard et al. 1999; del Burgo et al. 2003; Kiss et al. 2006; Ridderstad et al. 2006). The average size of aggregates required to fit the FIR, submillimetre, and extinction profiles in the L1506 filament is about 0.4 μm (Ysard et al. 2013). This value is close to the smallest grain size needed to scatter light efficiently in the MIR, which produces the “coreshine” observed toward a number of dense cores, which has also been interpreted as a result of grain growth (Pagani et al. 2010; Steinacker et al. 2010).

On the theoretical side, various dust optical property calculations have predicted a significant increase in the emissivity of aggregates at long wavelengths compared to compact spherical grains (Wright 1987; Bazell & Dwek 1990; Ossenkopf 1993; Ossenkopf & Henning 1994; Stognienko et al. 1995; Köhler et al. 2011, 2012). This variation is shown to be mainly due to the increase in the porosity fraction with aggregate growth, but the shape, structure, material composition, and accretion of mantles can also all contribute.

Moreover, Malinen et al. (2011), Juvela & Ysard (2012), Ysard et al. (2012) have investigated the impact of radiative transfer on the results derived from observations under the assumption of a single average colour temperature. They showed that the mixing of different temperatures along the line of sight produces a tendency that is opposite to the observed one. They concluded that in the absence of internal heating sources, the observed emissivity increase toward dense clouds cannot be explained by radiative transfer effects. It must originate in intrinsic variations of the optical properties of the grains.

It is, however, important to note that the dust emissivity increase is not systematically observed in the interstellar medium (ISM, see Nutter et al. 2008; Juvela et al. 2009; Paradis et al. 2009). These intriguing results call for a broader investigation, making use of the large observations statistics provided by Herschel and Planck, probing different Galactic environments. The key questions are still open today: when, where, and how dust evolves between diffuse and dense regions, what the physical conditions enhancing (or preventing) the efficiency of the coagulation process are, what the time scales are, and whether the process is directly related to specific stages in the cloud or core evolution. Understanding these questions is critical since knowing the dust opacity has a direct impact on many key parameters derived from dust emission, such as the column densities, masses, and volume densities of the clouds. For this reason, it is also necessary to investigate the possible systematic effects on the emission and extinction measurements that could cause errors in the opacity estimates.

The structure of the paper is as follows: the observations are described in Sect. 2. The main results are presented in Sect. 3, including the estimates of submillimetre and near-infrared (NIR) optical depths, the correlations between these variables, and the correlations with environmental factors. The results are discussed in Sect. 4 before we list the final conclusions in Sect. 5.

2. Observations and basic data analysis

2.1. Target selection

The selection of the Herschel fields is described in Juvela et al. (2012) and an overview of all the maps is given in Montillaud (2015). We only repeat the main points here.

Planck sub-millimetre observations, together with IRAS 100 μm data, enabled the detection of over 10 000 compact sources in which the dust is significantly colder than in the surrounding regions (Planck Collaboration XXIII 2011). The detection procedure is based on this colour temperature difference and, furthermore, limits the size of the detected clumps to values below ~12′ (Montier et al. 2010). The sources are believed to be Galactic cold clumps or, at larger distance, entire clouds (Planck Collaboration XXIII 2011).

The fields for Herschel follow-up observations were selected using a binning of Planck cold clumps with respect to the Galactic longitude and latitude, the estimated dust colour temperature, and the clump mass. At the time of source selection, distance estimates existed for approximately one third of the sources in C3PO and, therefore, some sources of unknown mass were also included. The binning ensured full coverage of the clump parameter space, especially of the high Galactic latitudes and of the outer Galaxy. Galactic latitudes | b | < 1° were excluded because that area is covered by the Hi-GAL programme (Molinari et al. 2010). Similarly, regions included in other Herschel key programmes such as the Gould Belt survey (André et al. 2010) and HOBYS (Motte et al. 2010) were avoided.

A total of 116 separate fields were observed. The SPIRE maps are on average over 40′ in linear size, with an average area of ~ 1800 (arcmin)². The PACS maps are smaller, with an average area of ~ 660 (arcmin)². Most fields contain more than one Planck clump, the maps altogether covering ~350 individual Planck detections. The range of probed column densities extends from diffuse fields with N(H₂) ~ 10²¹ cm^-2 to cores with N(H₂) > 10²³ cm^-2. The fields are listed in Table 1 and the Herschel observation numbers are included in Table E.1.

2.2. Herschel data

2.2.1. Herschel data reduction

The fields were mapped with the SPIRE instrument at wavelengths 250, 350, and 500 μm and with the PACS instrument at wavelengths 100 and 160 μm. One field was observed with SPIRE alone (G206.33-25.94, part of the Witch Head Nebula, IC 2118). The Herschel observations are discussed in detail in Juvela et al. (2012) and Montillaud (2015). The SPIRE observations at 250 μm, 350 μm, and 500 μm were reduced with the Herschel Interactive Processing Environment HIPE v.10.0.0, using the official pipeline with the iterative destriper and the extended emission calibration options. The maps were produced with the naive map-making routine. The PACS data at 100 μm and 160 μm were processed with HIPE v. 10.0.0 up to level 1, and the maps were then produced with Scanamorphos v20 (Roussel 2013). In the order of increasing wavelength, the resolution of the maps is 7″, 12″, 18″, 25″, and 36″ for the five bands. The raw and pipeline-reduced data are available via the Herschel Science Archive, the user-reduced maps are available via ESA site¹.

The accuracy of the absolute calibration of the SPIRE observations is expected to be better than 7%². For PACS we assume a calibration uncertainty of 15%. This is a conservative estimate that is compatible with the differences of PACS and Spitzer MIPS measurements of extended emission³.

2.2.2. Estimating intensity zero points

To determine the zero point of the intensity scale, we compared the Herschel maps with Planck data complemented with the IRIS version of the IRAS 100 μm data (Miville-Deschênes & Lagache 2005). The Planck and IRIS measurements were interpolated to Herschel wavelengths using fitted modified blackbody curves, B_ν(T_dust)ν^β, with a fixed value of the spectral index, β = 2.0. The linear correlations between Herschel and the reference data were extrapolated to zero Planck (+IRIS) surface brightness to determine offsets for the Herschel maps. For SPIRE the uncertainties of these fits are typically ~1 MJy sr^-1 at 250 μm and at longer wavelengths smaller in absolute value. The derived intensity zero points are independent of Planck calibration and of any multiplicative errors in the comparison. For PACS the correlations are often less well defined, and the zero points were set directly based on the comparison of the average values of the Herschel maps and the corresponding interpolated Planck and IRIS data. The zero points were calculated iteratively, including colour corrections calculated using colour temperatures that were estimated from SPIRE data with a fixed spectral index of β = 2.0.

2.2.3. Calculating submillimetre optical depth

The Herschel maps were converted into estimates of dust optical depth at 250 μm. The surface brightness maps were convolved to a common resolution of 40″, and colour temperatures were calculated by fitting the spectral energy distributions (SEDs) with modified blackbody curves with a constant opacity spectral index of β = 2.0. The 250 μm optical depth was obtained from $\begin{array}{ccc}\mathit{\tau }\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}\frac{{\mathit{I}}_{\mathit{\nu }}\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}}{{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\mathit{,}& & \end{array}$(1)using the fitted 250 μm intensity I_ν(250 μm) and the colour temperature T. The calculations were made with 250–500 μm data and 160–500 μm data. The fits were weighted according to 15% and 7% error estimates for PACS and SPIRE surface brightness measurements, respectively (see Sect. 2.2.1).

The assumed value of β = 2.0 may be appropriate for dense clumps, although at lower column densities the average value is lower, β ~ 1.8 (Boulanger et al. 1996; Planck Collaboration XXV 2011), and the value of β may further depend on the Galactic location and the wavelength range (e.g. Planck Collaboration Int. XIV 2014). If the true value of β were 1.8 instead of 2.0, our colour temperature estimates would be higher by ~1 K and the τ(250 μm) values lower by ~30%. Furthermore, if the values of β were correlated with column density, the slope of τ(250 μm) vs. τ_J would be similarly affected. We return to these effects in Sects. 3.7 and 4.

To estimate the statistical uncertainty of τ(250 μm) values, we used Markov chain Monte Carlo (MCMC) runs. The prior distribution of temperature values is flat but limited between 5.0 K and 35 K. In addition to the relative errors quoted above, we included the uncertainty of the intensity zero points. These are typically much smaller than the assumed relative errors, but may be important at low column densities, especially at 160 μm. The zero-point errors are systematic but are included simply as another component of statistical noise. Their effect is thus reflected in the error estimates of individual pixels. The error distribution of τ(250 μm) is nearly Gaussian, and we used the standard deviation of the MCMC τ(250 μm) samples as the error estimates. These estimates were calculated separately for each pixel of the τ(250 μm) maps.

Because of line-of-sight temperature variations, the derived τ(250 μm) estimates probably systematically underestimate the true values (Shetty et al. 2009; Malinen et al. 2011). We cannot directly determine the magnitude of these errors but, with some assumptions, we can use radiative transfer modelling to estimate the magnitude of the bias. The simulations, described in detail in Appendix C, were used to derive bias maps that are taken into account when the data were correlated with τ_J values.

2.3. Near-infrared data

We used the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006) to derive estimates of dust column density that are independent of dust emission. We used the method NICER (Lombardi & Alves 2001) and the standard extinction curve (Cardelli et al. 1989) to convert the reddening of the background stars to estimates of J-band optical depth, τ_J. Because the calculations involve only NIR bands, the results are expected to be insensitive to the value of the ratio of total to selective extinction, R_V (Cardelli et al. 1989). The shape of the NIR extinction curve is believed to be relatively stable, even at high extinctions (e.g. Draine 2003a; Indebetouw et al. 2005; Lombardi et al. 2006; Román-Zúñiga et al. 2007; Ascenso et al. 2013; Wang & Jiang 2014). Some variations are observed with Galactic location and/or density, but generally only at a level of 5% of the NIR power-law index (e.g. Stead & Hoare 2009; Fritz et al. 2011). This question is discussed in more detail in Sect. 4.2. The τ(J) values are derived using both the J−H and H−K colours but, with the extinction curve used, we have the correspondence of E_J−K = 0.65 τ(J). Flags in the 2MASS catalogue were used to avoid galaxies (ext_key not null or gal_contam not zero) and sources with uncertain photometry (ph_qual worse than C).

Five of our fields are fully covered in the VISTA Hemisphere Survey, VHS (McMahon et al. 2013), which has more than ten times the sensitivity of 2MASS (in H band VHS has a 5σ detection threshold of 19.0, compared to a 2MASS point source catalogue completeness limit of ~16 mag). One of these fields is too distant to obtain a reliable extinction map, but the data for the four other fields were analysed and the results compared with those obtained with 2MASS data. The fields are G4.18+35.79 (LDN 134), G21.26+12.11, G24.40+4.68, and G358.96+36.75. For the fields G21.26+12.11 and G24.40+4.68, only J- and Ks-band data exist. The data are available in VISTA Survey Archive⁴, and the VISTA Data Flow System pipeline processing and science archive are described in Irwin et al. (2004) and Hambly et al. (2008).

Extinction maps are produced by averaging extinction estimates of individual stars with a Gaussian weighting function with FWHM = 180″. We also tested a higher resolution of FWHM = 120″. For distant sources, the extinction of the target clouds cannot be reliably reproduced because of the poor resolution and the increasing number of foreground stars. This is the main factor that limits the number of fields where the ratio τ(250 μm) /τ(J) can be reliably estimated. The extinction measurements can be significantly biased even in nearby fields if these contain steep column density variations. No special steps were taken to eliminate the contamination by foreground stars (see, e.g., Schneider et al. 2011), apart from the sigma clipping procedure that is part of the NICER method and was performed at 3σ level. The reliability of the extinction maps and the bias caused by sampling problems and the presence of foreground stars was examined with simulations (see Sect. B). The results of these simulations are used to derive maps of the expected uncertainty and the bias of the τ(J) values for each field.

2.4. Correlations between sub-millimetre and NIR opacity

The ratio k of sub-millimetre opacity τ(250 μm) and the NIR opacity τ_J was estimated for all 116 fields. The τ(250 μm) maps were convolved to the 3′ resolution of the τ_J maps. The τ(250 μm) and the τ(J) data were read at 90″ steps (half-beam sampling), excluding the map borders where the result of the convolution to 3′ resolution is poorly defined. For local background subtraction, only areas where the signal was more than 2σ above the average value of the reference area were used (see Sect. 2.2.2). Here σ is the standard deviation of the values in the reference region. This is a conservative limit because part of the fluctuations is caused by real surface brightness variations and not by noise alone.

For τ_J the error estimates were provided by the NICER routine. For τ(250 μm) these were obtained from MCMC calculations (see Sect. 2.2.3). The comparison between the different cases (for example, regarding the use of 160 μm data, background subtraction, or gradient corrections) provides information on the uncertainty caused by some sources of systematic errors.

The τ(250 μm) vs. τ(J) points of individual fields and samples of fields were fitted with a linear model to derive the ratio τ(250 μm) /τ(J). These total least-squares fits take into account the uncertainties in both variables. The fits were made using either all data points or only data below or above a given τ(J) limit. To reduce the bias caused by these cuts, the data were divided with the help of a preliminary linear fit to all data points (see Sect. 3.1 for details). The limiting value of τ(J) thus corresponds to a position on this line, and the cut was performed using a line that is perpendicular in a coordinate system where the average uncertainties of the two variables are equal.

The τ(250 μm) /τ(J) ratios were also calculated for alternative versions of the τ(250 μm) data, using local background subtraction or using ancillary data in an attempt to correct for possible large-scale errors in the surface brightness data. These alternative data are discussed in Appendix A.

3. Results

3.1. Apparent τ(250μm)/τ(J) values

We calculated τ(J) and τ(250 μm) maps of the 116 fields as described in Sects. 2.2 and 2.3. The correlations between τ(J) and τ(250 μm) were calculated at a resolution of 180″. In addition to the full range of column densities, the relationships were examined separately below and above the limit of τ(J) = 0.6 (see Sect. 2.4). This corresponds to visual extinctions A_V ~ 2.3 mag and A_V ~ 2.0 mag for the R_V values of 3.1 and 5.0, respectively (Cardelli et al. 1989). Instead of a higher limit, we selected the relatively low number of τ(J) = 0.6 to maximise the number of fields where a linear fit could also be made above the τ(J) threshold. The τ(250 μm) values were derived from Herschel data with either 250−500 μm or 160–500 μm (see Appendix A for analysis with additional alternative data sets).

In a given field, the number of points either below or above the τ(J) limit is often insufficient to determine any reliable value for the slope k = Δτ(250 μm)/Δτ(J). In a few fields no reliable value of k can be determined at all, mainly because of the low quality of the τ(J) data. This especially affects the most distant fields because of the contamination by foreground stars and because the structures are too small to be resolved with the 3′ beam. The formal errors of the k parameter were used to exclude the clearly unreliable fits. The criterion δk/k< 0.1 leaves in the default case 106 fits to all data in a field, 103 fits below τ(J) = 0.6, and 38 fits above τ(J) = 0.6. These fits appear relatively reliable also based on visual inspection.

Figure 1 shows an example of the recovered dependence between τ(J) and τ(250 μm) values, including linear fits to the three τ(J) ranges. In this example, the slope appears to become steeper as τ(J) increases. This might be an indication of an increase in the dust submillimetre opacity, which in turn might be attributed to grain growth (e.g., Ossenkopf & Henning 1994; Stepnik et al. 2003; Ormel et al. 2011; Ysard et al. 2013). However, before drawing any such conclusions, we must consider the systematic effects that affect the two parameters. Figure 2 shows a summary of all the Δτ(250 μm)/Δτ(J) values where τ(250 μm) values are based on Herschel 250−500 μm data. Before any bias corrections (see below), the values are seen to cluster around ~ 2.0 × 10^-3, with some tendency for higher values in the higher τ(J) range.

Figure 3 summarises the statistics of the dust opacity measurements as histograms, including all fits where the formal error of the slope of the least-squares fit τ(250 μm) vs. τ(J) is below 10%.

We need to observe a sufficient number of background stars for each resolution element, even at high column densities. This means that τ(J) estimates and the comparison with submillimetre emission can be made only at a low resolution (2–3′). Averaged over such large areas, the statistical uncertainty of Herschel data is very small. In Sect. 3.2 we show that the bias is probably also dominated by the errors in τ(J). In the following, we rely mainly on the Herschel data set that consists of observations 250−500 μm (the “default” data set). There are three reasons. First, in theory, the inclusion of the 160 μm data reduces statistical uncertainty of the colour temperature estimates but increases systematic errors caused by line-of-sight temperature variations (Shetty et al. 2009; Malinen et al. 2011). Second, because of the smaller (and, for parallel mode, different) area covered by the PACS observations, the use of the 160 μm band significantly reduces the area where correlations with τ(J) can be calculated. Third, 160 μm data may be affected by additional systematic effects related to the relative calibration of the two instruments, uncertainties in the zero-point determination (interpolation between IRAS and Planck channels and the contribution of stochastically heated grains in the IRAS 100 μm band) and to imperfections in the map making that could be increased by the smaller size of the PACS maps (see Sect. A.2). We are particularly interested in the coldest regions where 250–500 μm data provide adequate constraints on the dust temperature.

3.2. Bias in τ(J) values

Bias in τ(J) values is very likely a significant problem, especially for distant fields in which all high column density structures are not resolved and the results begin to be affected by foreground stars. Both effects decrease the τ(J) estimates, especially towards column density peaks. We estimated the extent of the problem with simulations using the stellar statistics in low-extinction areas near each field. The contamination by foreground stars was evaluated with the help of the Besançon model of the Galactic stellar distribution (Robin et al. 2003). We used the Herschel column density maps as a model of the column density structure, simulated the distribution of foreground and background stars, analysed the simulated observations with NICER routine, and compared the results with the known input τ(J) map. The procedure is described in detail in Appendix B. We obtained for each field a map of the expected systematic relative error in τ(J) that gives a multiplicative correction factor for the τ(J). The simulations do not consider the effect of cloud structures at scales below 18″ but the procedure probably provides a reasonable estimate of the magnitude of the effect.

We repeated the analysis of the previous section and replaced the original τ(J) maps with bias-corrected estimates. Figure 4 compares the Δτ(250 μm)/Δτ(J) distributions for the default case with and without bias correction. The statistics include all fits for which the formal error of the slope is below 10%. This corresponds to Fig. 3, but the number of points is different. Because the bias corrections depend on the cloud distance, fields without distance estimates had to be dropped. However, in the τ(J) > 0.6 interval the number of fields fulfilling the δk/k< 0.1 criterion has doubled.

3.3. Bias in τ(250μm) values

The systematic errors in τ(250 μm) values were estimated with radiative transfer modelling. The line-of-sight temperature variations are expected to be the main source of error that, in standard analysis, leads to overestimation of the mass-averaged dust temperature and, subsequently, to underestimation of τ(250 μm) (e.g. Ysard et al. 2012).

We constructed for each field a three-dimensional radiative transfer model that covered a projected area of 30′× 30′ with a 10″ pixel size. The modelling assumed spatially constant dust properties, the dust model (Draine 2003b) corresponding to R_V = 5.5 (see Appendix C for details). The density distribution and external heating were adjusted until the model exactly reproduced the observed 350 μm surface brightness and, for the area above median column density, the average 250 μm/500 μm ratio. The model-predicted surface brightness maps were analysed as in the case of the actual observations, to produce maps of τ(250 μm). To estimate the bias, these values were compared to the actual τ(250 μm) values of the model to derive multiplicative correction factors.

The results depend on the assumed cloud structure in the line-of-sight direction (Juvela et al. 2013). In our models, the line-of-sight density distribution only has one peak. This enhances temperature contrasts and increases our bias estimates. On the other hand, the densest observed cores are probably even more compact, and in their case we may be underestimating the bias. If the clouds contain embedded sources, the actual bias may again locally be very different and often lower than predicted by our models. Although the bias estimation is more difficult than for τ(J), the models should again provide a reasonable estimate of the magnitude of the effect. The relative systematic errors are smaller in τ(250 μm) than in τ(J) so that their effect on the final result is less strong.

The k = Δτ(250 μm)/Δτ(J) values were re-calculated including bias corrections in both variables. The resulting histograms are included in Fig. 4. The bias correction makes the distributions significantly narrower. Figure 4 also shows that the corrections are much stronger for τ(J) than for τ(250 μm). As a result, the average value of Δτ(250 μm)/Δτ(J) now decreases with increasing τ(J). The number of fields fulfilling the δk/k< 0.1 criterion has doubled to 76 fields (using SPIRE bands). Compared to the original data, the median values of k × 10⁴ have decreased from 20.2, 21.1, and 23.4 to 15.3, 16.0, and 12.2, the numbers corresponding to the full τ(J) range, data below τ(J) = 0.6, and data above τ(J) = 0.6, respectively. The strong change in the k values suggests that the uncertainty of k is probably often several tens of per cent, especially in the τ(J) > 0.6 interval. Therefore, Fig. 4 does not exclude a systematic increase of k as the function of τ(J), if that becomes visible only at high column densities.

Figure 5 displays the slopes k = Δτ(250 μm)/Δτ(J) as function of distance and Galactic location. The left frames show the relations without and the right frames with the bias corrections applied to τ(J) and τ(250 μm). Only fits with δk/k< 0.1 are included. The original data showed some trends, including an increase in k as a function of distance and galactocentric distance. The first is visible especially in the high τ(J) interval, but was expected because τ(J) values of distant sources can be severely underestimated. After bias corrections, the scatter of the k values is significantly reduced. This suggests that the corrections are of correct magnitude. The distance dependence has changed so that in the corrected data there is a slight decrease in k values as a function of distance. This could point to some over-correction of the τ(J) estimates, although the bias correction should not only depend on distance, but even mainly on the cloud structure. However, the decrease of k can be an indication of selection effects or direct resolution effects. For example, higher k values might be found in individual dense clumps that are only resolved at short distances.

There is little difference between the k values found in the three τ(J) intervals. In the next section we examine in more detail the global τ(J) dependence of the τ(250 μm) /τ(J) ratios, especially regarding the highest observed column densities.

3.4. Global relation τ(250μm) vs. τ(J)

To further test the hypothesis that τ(250 μm) /τ(J) ratios may change systematically as a function of column density, we carried out non-linear fits τ(250 μm) = A + B × τ(J) + C × τ(J)². The fits were first performed using the combined data of all fields. To reduce the mismatch in the zero levels of individual fields, we subtracted from each τ(J) and τ(250 μm) map the local background using the off regions listed in Table 1. The off regions are not completely void of emission but provide a common reference point for the quantities. Thus, the relation is expected to develop via the origin for each field separately, the parameter A being close to zero for the combined data as well. The sign of the fitted parameter C indicates the possible increase or decrease of k = Δτ(250 μm)/Δτ(J) as the function of column density.

Figure 6 shows the results obtained with the bias-corrected data. We included data from all fields in which individual linear fits had δk/k< 0.2, thus relaxing the previous constraint of δk/k< 0.1. The second-order polynomial was fitted to all data and separately to data points with τ(J) > 1.0. In the previous section a threshold value of τ(J) = 0.6 was used. However, some 70% of all points are below τ(J) = 0.6 and, when included, they dominate the fits that systematically underestimate the data above τ(J) ~ 5. With the combined data set, there are enough high column density data points so that the lower limit can be moved upwards. The use of the τ(J) = 1.0 threshold enables an adequate fit to all data with higher τ(J) values. The τ(J) calculations employ a different off region for each field. These may contain different amounts of extinction, which leads to small relative shifts along the τ(J) axis. Based on dust emission, the extinction in the off regions is typically τ(J) = 0.2−0.4. The uncertainty of the relative zero points contributes to the scatter in Fig. 6, but the effect is weaker than the total dispersion and the non-linearity seen at high extinctions.

To prevent the τ(J) cut itself from biasing the fits, the data were selected using lines perpendicular to a linear least-squares line fitted to all data (cf. Sect. 2.4). Thus, the quoted τ(J) limits correspond to a point on the fitted line, and the cut itself is perpendicular to the fitted line. All fits take into account the uncertainties in both variables, which are assumed to be uncorrelated. The error distributions of the parameters A–C were calculated with an MCMC.

The sign of the parameter C depends on the range of τ(J) values but is less dependent on the field selection, for example regarding the δk/k limit that was used to select the fields. Most fields with high τ(J) values (and thus with a wide dynamical range) have low values of δk/k. When all points are included, the value of the parameter C is negative, but the fit is very poor at high τ(J). When pixels τ(J) < 1.0 are excluded, C becomes positive (see Fig. 6a), which points to an increase in the submillimetre opacity above τ(J) ~ 1. Systematic additive errors in either parameter might also explain the different behaviour at very low τ(J). When the fit is made using all data τ(J) > 0.6 (not shown), the parameter C is marginally positive, but beyond τ(J) = 10 the fitted line is below all the data points. The second frame of Fig. 6 shows the fits when pixels with colour temperatures above 14 K are excluded. The values of C are now higher and positive even when data τ(J) < 1.0 are not excluded. The best fit to the high τ(J) end of the relation is still obtained by excluding data with τ(J) < 1.0, this results in the relation $\mathit{\tau }\mathrm{(}\mathrm{250}\mathit{\mu }\mathrm{m}\mathrm{)}\mathrm{=}\mathrm{0.73}\mathrm{\times }{\mathrm{10}}^{-3}\mathrm{+}\mathrm{1.25}\mathrm{\times }{\mathrm{10}}^{-3}\hspace{0.17em}\mathit{\tau }\mathrm{\left(}\mathit{J}\mathrm{\right)}\mathrm{+}\mathrm{0.11}\mathrm{\times }{\mathrm{10}}^{-3}\hspace{0.17em}\mathit{\tau }\mathrm{(}\mathit{J}{\mathrm{)}}^{\mathrm{2}}\mathit{.}$(2)The formal error estimates of the parameters A–C are of the order of 5%, probably lower than the systematic uncertainties. All data beyond τ(J) ~ 5 are affected by large bias corrections and, consequently, the value of C also depends on the accuracy of these corrections. Thus, Fig. 6 strongly suggests but does not yet provide a final proof of the variations of the ratio τ(250 μm) /τ(J). The positive offset A = 0.7 × 10^-3 results from the facts that at low column densities the relation is linear, the curvature increases only beyond τ(J) ~ 5, and the lowest data points τ(J) < 1.0 are not part of the fit.

Figure 7 shows the error distributions of the parameters A–C. The fit was made to data τ(J) > 1.0 for all fields with a linear fit accuracy δk/k< 0.2. In most fields, the formal error estimates of δτ(J) and δτ(250 μm) are smaller than the actual scatter of points. Therefore we used the residuals of the linear fits before the MCMC calculation to determine a scaling factor, typically 2.0–3.0, that makes the error estimates in each field consistent with the actual scatter. Even after this increase of uncertainties, MCMC gives a 100% probability for a positive value of C. In reality, the result is not that strong because the uncertainty may be dominated by systematic errors. The sign of C was already seen to change depending on the range of τ(J) values fitted. The result also depends on a relatively small number of fields with data above τ(J) > 5. Therefore, we must consider the τ(250 μm) vs. τ(J) relation in individual fields in more detail.

3.5. Correlations in selected fields

Figure 6 showed hints of an increase of the τ(250 μm) /τ(J) values as the column density increases, but the global statistics may be confused by the mix of different fields. Furthermore, the sign of the C parameter is determined by the highest τ(J) points that originate in a small number of individual fields. Each field may be affected by different systematic effects related to the surface brightness zero points, distance uncertainty (via bias correction), and differences in the local radiation field. Several diffuse fields are even entirely below τ(J) ~ 1.0. Therefore, we also need to examine the fields individually.

Three criteria were used to select a subset of fields. We required that (1) the uncertainty of the fitted parameter C is below 0.3 × 10^-4; (2) there are at least ten data points (selected from the maps at 90″ steps) with τ(J) above 0.6; and (3) the bias corrections change the slope of the linear fit of τ(250 μm) vs. τ(J) by less than 30%. The first two criteria ensure that there are enough data points at large τ(J) with a small scatter to gain some insight about the column density dependence of the Δτ(250 μm)/Δτ(J). The third criterion excludes distant fields for which the uncertainty of the bias correction of τ(J) renders the results uncertain, even for the apparently well-defined relation between τ(J) and τ(250 μm). The selection leaves 23 fields with distances mostly in the range 100−500 pc. There are two exceptions, G216.76-2.58 and G111.41-2.95, for which the estimated distances are 2.4 kpc and 3.0 kpc. We kept the two fields in the sample even though the results are known to be unreliable because of the large distance.

To avoid underestimating the fit errors, we scaled the error estimates of τ(J) and τ(250 μm) up to correspond to the actual scatter of points (see Sect. 3.4). All observations, sampled at 90″ steps, are fitted with a linear model τ(250 μm) = b + k × τ(J) using total least-squares. We continued to use as the default data set one with τ(250 μm) derived from SPIRE bands alone, including bias corrections in τ(250 μm) and τ(J). However, for comparison, we also examined results obtained with four Herschel bands (160−500 μm; including the bias corrections) and, finally, with three Herschel bands but without any bias corrections.

The non-linear fits were made using MCMC (with 2 × 10⁵ samples per field) and bootstrap sampling (2000 realisations per field). Both fits used total least-squares⁵ and the error estimates of the individual points. The parameter values and the uncertainties derived with these two methods should be similar, except for rare cases in which the result depends on a small number of influential points, which are always present in the MCMC calculation, but not in all bootstrap samples.

Figures 8 show the results of these fits. Each frame shows the values of the linear slopes for the three cases discussed above. The parameters of the non-linear fits are shown together with the error estimates calculated with the bootstrap method. In addition to the fit to the default data set (solid red curve), the dashed magenta lines show the effect of the distance uncertainty. Using the distance uncertainties δd listed in Table 1, we also calculated the bias corrections for τ(J) for distances d−δd and d + δd. Thus, the upper dashed line corresponds to distance d−δd and a smaller bias correction.

Figure 9 shows the linear slopes k and the values of parameter C for this sample of fields. The fits were performed to all data, without a τ(J) threshold. If the data below τ(J) = 0.6 were removed, the median slope τ(250 μm) /τ(J) = 1.6 × 10^-3 did not change appreciably (by less than 0.1 × 10^-3). The median value of C increased to 2.0 × 10^-4, which in that case is still lower than the scatter. The differences between the fields are larger than the estimated formal uncertainties (including the statistical errors of τ(J) given by NICER and τ(250 μm) derived from the uncertainty of the surface brightness measurements). The error bars only reflect the statistical errors of the fits and do not include the uncertainty of the bias correction, for example. Figure 8 demonstrates that the uncertainty of the distances can be a significant source of error. In Fig. 9, the shaded areas show the difference between the values obtained with distances d−δd and d + δd, as listed in Table 1. These were estimated directly by repeating the analysis using these two distance values. A smaller distance corresponds to smaller bias correction in τ(J) and, thus, to a steeper slope k and typically a higher value of C. In some cases, the value of C obtained with the default distance d is outside the shaded region, showing that the effect is not always this simple. The distance uncertainty is not yet enough to explain all the scatter in k and especially in C. A change in the distance estimate results at first approximation in a nearly linear scaling of τ(J) values (see Fig. 8, comparison of the dashed magenta lines). In reality, the situation may be more complex. In particular, if a field contains cloud structures at different distances, this might result in large errors in both k and C. Figure 9 also shows that in spite of the small formal errors of the least-squares fits, we cannot constrain the opacity values in the last two fields with distances exceeding 2 kpc.

In the sample of Fig. 9, the median value of C is close to zero with a number of fields with negative values. The positive values of C in Fig. 6 are due to a small number of fields, and the increase of τ(250 μm) /τ(J) values was only visible above τ(J) ~ 5. There are only 20 fields with any data points above τ(J) = 5. Only six fields have ten or more data points above this limit: G6.03+36.73, G70.10-1.69, G82.65-2.00 G92.04+3.93 G107.20+5.52, and G202.02+2.85. Of these, only G6.03+36.73 is included in the sample of Fig. 9. All the others were excluded because the bias correction changed the slope k by more than 30%. Thus, a clear steepening of the relation τ(250 μm vs. τ(J) is seen exclusively in fields with the highest column densities, which for the same reason also have the largest uncertainty regarding the bias corrections.

3.6. Maps of τ(250μm) /τ(J) ratio

We also examined the ratios τ(250 μm) /τ(J) in the form of maps. This is useful if k changes in small regions that have little effect when all data of a field are fitted. Unlike in Fig. 8, where the offset between τ(250 μm) and τ(J) is a free parameter, the appearance of the ratio maps depends on the consistency of the τ(250 μm) and τ(J) zero points. Because we do not have an absolute zero point for τ(J), we used the reference areas listed in Table 1 and subtracted from τ(250 μm) and τ(J) the average value found in the reference area. This limits the region where a reliable ratio can be calculated, excluding regions of low column density. The details of the calculations are given in Appendix E, where we also show the figures of selected fields. As an example, Fig. 10 shows the field G4.18+35.79 (LDN 134), where the ratio τ(250 μm) /τ(J) is strongly correlated with column density. In Fig. 10, the increase of k remains clear even in maps of (τ(250 μm) ± δτ(250 μm)/(τ(J) ± δτ(J)). The error estimates δτ(250 μm) include the statistical errors due to Herschel photometry and uncertainty of the surface brightness zero point (see Sects. 2.2.3 and 2.2.2). The parameter δτ(J) corresponds to the uncertainty of the τ(J) zero point (see Appendix E). For τ(J) the formal error estimates calculated with NICER are below 10%.

For high-opacity sources like G4.18+35.79 the bias correction of τ(J) is very important. If no background stars are visible through some part of the core, the values of k naturally remain more uncertain (see, however, Sect. 3.7). Conversely, the zero-point uncertainty only becomes important in diffuse regions but might even reverse the correlation with column density. The ratio maps are also affected by the assumption of a constant value of β and of potential errors in the bias corrections. However, we argue in Sect. 3.7 that these are mainly multiplicative errors (that do not affect the morphology of the maps) and/or tend to decrease the variations seen in the ratio maps. Therefore, we are confident that the increase of submillimetre opacity that is seen in some of the maps is real.

Based on the maps, the submillimetre opacity is correlated with column density in the fields G4.18+35.79, G6.03+36.73, G111.41-2.95, G161.55-9.30, G151.45+3.95, and G300.86-9.00 (see Appendix E). In G4.18+35.79 and G6.03+36.73 the values rise to close to 4 × 10^-3. In some fields the background subtraction reduces the available map area to such an extent that no conclusions can be drawn.

So far, all NIR extinction maps were calculated at 180″ resolution. Depending on the number of background stars, extinction map could be derived at a higher resolution and possibly with smaller bias. This especially applies to the four fields for which VISTA observations are available. We recalculated the extinction maps at 120″ resolution, repeating Monte Carlo simulations to estimate the bias of τ(J). The smaller beam increases the noise per resolution element, but does not yet cause holes in the extinction maps. The analysis was repeated for the four fields with VISTA data with resolutions of 180″ and 120″. The results are summarised in Fig. 11. The resolution has no strong systematic effect on the parameters. The largest differences appear when parameter C is estimated excluding low column density points. However, even in that case the difference between the results with a resolution of 120″ and 180″ is smaller than the effect of excluding low τ(J) data from the fits. When VISTA data are available, the results are close to those obtained with 2MASS. Because we used the same Herschel data and bias corrections derived in the same way, the results are not independent. However, because the uncertainty of τ(J) is expected to be one of the most significant sources of error, this gives us some confidence that the observed differences between the fields are real.

The τ(250 μm) /τ(J) ratio in the field G4.18+35.79 was shown in Fig. 10. Figure 12 shows the corresponding figure obtained with VISTA NIR data and a spatial resolution of 120″. The highest value of τ(250 μm) /τ(J) has increased from 3.6 × 10^-3 in Fig. 10 to 6.7 × 10^-3. This is mainly attributed to the increased spatial resolution, although also at the 180″ resolution the peak value is ~25% higher than in Fig. 10. Even in VISTA data there are only ~10 stars within the 2′ × 2′ area centred on the τ(250 μm maximum, and therefore the peak value of τ(250 μm) /τ(J) is subject to some uncertainty.

3.7. Potential systematic errors

Because of the significance of the bias corrections (Sects. 3.2 and 3.3), we tried to characterise the effects that systematic errors in these corrections could have on the τ(250 μm) /τ(J) ratios. Furthermore, the assumption of a constant dust emissivity spectral index may be incorrect. Below we examine the possible systematic effects caused by these factors.

The bias correction made to the τ(250 μm) values is in itself small (see Fig. 4) and, consequently, the errors made in that correction are expected to be small. The correction was derived from radiative transfer models with dust properties corresponding to R_V = 5.5 (Draine 2003b). If the submillimetre dust emissivity is in fact higher, the dust column density will be overestimated and models will also overestimate the cloud opacity at visual and NIR wavelengths, thus exhibiting stronger temperature variations than the real clouds. Because the τ(250 μm) bias is related to the line-of-sight temperature variations, we could in this case systematically overestimate the bias in τ(250 μm). To check the magnitude of the effect, we repeated the modelling using a dust model with higher long-wavelength emissivity. We used a dust model from Ossenkopf & Henning (1994) with coagulated grains with thin ice mantles accreted in 10⁵ yr at a density of 10⁶ cm^-3. Compared to NIR (the wavelengths contributing to most of the heating deep inside a cloud), this dust model has an emissivity higher by ~50% at SPIRE wavelengths. The results are shown in Fig. 13, the open circles corresponding to this alternative modelling. Because of the smaller estimated bias, these should be below the k values of our previous analysis (red symbols). The median value of k is 1.58 × 10^-3 and thus practically unchanged. The strongest change is seen for G6.03+36.73 (LDN 183), where the estimate has been reduced by more than 20%. This is a source with very high column density and thus, in its central parts, a very uncertain estimate of τ(250 μm). However, the uncertainty of τ(250 μm) bias must also be considered in connection with the uncertainty of the τ(J) correction, which could compensate for some of the change (see below).

The τ(J) bias corrections are more significant than the τ(250 μm) bias corrections. The estimation of the τ(J) bias is, in principle, more reliable because it only depends on the assumption of the τ(J) structure of the clouds. In our calculations (see Sect. 3.2), τ(J) was derived from Herschel observations, dividing τ(250 μm) by the constant factor of k = 1.5 × 10^-3 to obtain a template map of τ(J). There are two possible sources of error. First, if the targets contain much structure below the 18″ resolution of the τ(250 μm) maps, we will underestimate the bias and our k estimates will be too high. We cannot directly estimate this error, but it is expected to be a small fraction of the total bias estimate. This is because 18″ ≪ 3′ and, at least for the closest fields, Herschel already resolves most of the cloud structure. The second potential source of error is again connected with dust emissivity. If the local ratio between τ(250 μm) and τ(J) is higher than 1.5 × 10^-3 (strongly increased submillimetre opacity), we have overestimated the cloud opacity at NIR wavelengths in the modelling, the bias correction of τ(J) is too large, and we underestimate the true value of k. Thus, a change in the value of k that is used to estimate the τ(J) bias will change the recovered value of k in the same direction. To examine this potential problem more quantitatively, we repeated the bias estimation using k values of 1.0 × 10^-3 and 3.0 × 10^-3. The resulting values of k and C parameters are shown in Fig. 13 as triangles. The initial assumption of k = 1.0 × 10^-3 leads to a recovered median value of k = 1.26 × 10^-3. The initial assumption of k = 3.0 × 10^-3 leads to a recovered median value of k = 1.92 × 10^-3. In both cases the input and output values are inconsistent, unlike in our previous analysis, where an assumption of 1.5 × 10^-3 led to a recovered value of 1.6 × 10^-3. Furthermore, an error in the assumed value of k leads to a systematic error in the recovered value that is about half of the original error, even lower if the true value of k was initially overestimated.

The previous test shows that the estimates of k will be biased towards the selected value of 1.5 × 10^-3. Our previously recovered median value of 1.6 × 10^-3 is thus not significantly affected (bias lower than 0.05), but the effect can be stronger for individual fields. For example, in G4.18+35.79 the estimate was 1.8× 10^-3, but the true value is probably higher by ~10%. The calculations could be iterated, field by field, to carry out the bias correction self-consistently with the final k estimate. However, the errors are typically below 10%, and rough estimates of their magnitude can be seen in Fig. 13.

There is a specific consequence of the way the τ(250 μm) and τ(J) bias corrections are implemented. If a cloud included regions of such a high opacity that no 2MASS stars were visible through the cloud, the ratio of τ(250 μm) and τ(J) would normally be overestimated. The resulting apparent increase of submillimetre opacity could thus be an artefact resulting from errors in extinction values. However, in our analysis we also correct τ(J) in this case based on the assumed opacity derived using the τ(250 μm) input map and the extinction calculated with simulated 2MASS stars. If the gap in the distribution of background stars is increased, the ratio τ(250 μm) /τ(J) does not continue to increase, but instead tends towards the assumed ratio, 1.5 × 10^-3. Higher values should thus not be the result of gaps in extinction data.

To investigate the potential effects of a spatially varying spectral index, we repeated the analysis using an ad hoc β(T) law to introduce β variations in all of our maps. We took the temperatures calculated with β = 2.0 and fixed new β values pixel by pixel using a functional dependence β = 2.0 × (T/ 15.0)^-0.24. We then repeated the full analysis, starting with the zero point and colour corrections and continuing with the calculation of colour temperatures and submillimetre opacity. We did not solve the (T, β) values (which are very susceptible to noise effects), but simply assumed that β could vary in a systematic way so that the values are higher when the dust temperature is lower. The parameters of the β(T) formula were selected so that β changes from ~ 1.8 in warm regions with T ~ 23 K to ~ 2.2 in the coldest spots T ~ 10 K (cf. Dupac et al. 2003; Désert et al. 2008; Planck Collaboration XXIII 2011; Paradis et al. 2010; Veneziani et al. 2010; Juvela et al. 2011). The average β is still close to the original β = 2.0, and we mainly examined the effects of correlated changes of β rather than the effects of absolute β values that can be estimated more directly. The crosses in Fig. 13 show the slopes k = Δτ(250 μm)/Δτ(J) and the parameters C obtained with the spatially varying β. The values are typically within ~10% of the values obtained with β = 2.0. The strongest changes of k are seen in the two closest fields, G4.18+35.79 and G6.03+3673, which have well-resolved clumps with very low temperatures. The increase of the slope values is ~20%. Note that this is mostly consistent with the general dependence between τ(250 μm) and β and not necessarily an effect of the spatial variation of β. In this respect, it is very interesting to note that the effect on the parameter C is weak. In other words, if the variations of β are as assumed above, this will be reflected in the slope τ(250 μm) /τ(J), but without a noticeable non-linearity in the τ(250 μm) vs. τ(J) relation.

All the above suggests an uncertainty of 10–20% in k and ~0.1 units in C. In particular, if the increase of dust opacity is associated with values β> 2.0, our highest k estimates could still systematically underestimate the true values of k because of the lower assumed value of β and because the τ(J) correction biases the k values towards 1.5× 10^-3.

4. Discussion

We have examined the submillimetre opacity by correlating the 250 μm optical depth τ(250 μm) with the NIR optical depth, τ(J), assuming the latter to be an independent tracer of the total dust column density. Because the comparison was made at a resolution of 2′, it is not sensitive to dense cores (A_V ≫ 10 mag), at which both the τ(250 μm) and τ(J) estimates would become very uncertain. Nevertheless, corrections for systematic bias in τ(250 μm) and especially in τ(J) are important.

The sample consists of the heterogeneous set of 116 Galactic fields that were mapped with Herschel as part of the Galactic Cold Cores project. The main objectives were to estimate the typical ratio of τ(250 μm) /τ(J) and to search for variations of this quantity, between the fields and as function of column density. Such variations were then related to differences in the properties of interstellar dust grains. For the present sample, high column densities also imply low dust temperatures and thus conditions where submillimetre dust opacity is expected to be enhanced by grain aggregation. The limited resolution means that we did not probe the full range of opacity variations, if these are partly limited inside compact cores.

4.1. Main results and their reliability

By restricting the analysis to ~20 fields for which the results appeared most reliable, we derived a median value of k = τ(250 μm) /τ(J) = 1.6 × 10^-3 (see Sect. 3.5 and Fig. 9). In Fig. 9, 50% of the most nearby fields had positive values of C, the multiplier of the second-order term, indicating some degree of positive correlation between column density and submillimetre emissivity measured by τ(250 μm). In the maps of the ratio k, the same tendency was very clear in only six cases. The low percentage is partly caused by the noise in τ(J), whose magnitude is strongly correlated with the cloud distances. For two of the best examples, G4.18+35.79 and G6.03+36.73, the ratio τ(250 μm) /τ(J) increases to ~ 4 × 10^-3, almost a factor of three higher than the median value. The peak values are uncertain because of the large bias corrections, but because of the low spatial resolution used, the strongest effect can be even greater.

No dependence on either Galactocentric distance or on Galactic height was observed (Fig. 5). The reliability of the estimates decreases with distance, but nevertheless, our sample extends over more than ~4 kpc in Galactocentric distance. No trends are seen; if they were present at 10% level, they should still be visible over the scatter of individual data points. The result might be affected by systematic errors in the bias correction. However, these probably depend either on the distance or on the morphology of the field and at first approximation are not expected to be different in the inner and in the outer Galaxy.

In Fig. 9, the only clear trend is the decrease of k as a function of distance (also visible as larger scatter around the Galactocentric distance of 8.5 kpc). As mentioned in Sects. 3.2 and 3.3, there are two possible explanations. First, the bias correction of τ(J) might be overestimated so that as the distance increases, the error increases and k becomes underestimated. This is probably not the main reason because the bias correction depends as much on column density values and column density gradients as on distance. The trend probably is a combination of selection and resolution effects. At 1 kpc the 180″ resolution corresponds to almost 1 pc in linear scale. Therefore, we measure the mean cloud properties for the most distant fields. In the nearby fields individual clumps are resolved and the slopes k reflect more the contrast between diffuse regions and compact, core-sized objects. Thus, the trend would be compatible with the hypothesis that τ(250 μm) /τ(J) increases in the densest and coldest regions of interstellar clouds.

The global non-linear fits of Fig. 6 also indicated an increase of the ratio k = τ(250 μm) /τ(J) as the function of column density. The plots show clear deviations from a linear dependence beyond τ(J) ~ 4. For τ(J) ~ 10, k is already twice as high at low column densities. The results are dependent on a few fields with the highest column densities for which the bias corrections reach about ten per cent. However, if the distance dependence of Fig. 5 means that the bias correction of τ(J) is probably overestimated and not underestimated, the true values of k might be even higher.

Figure 9 concentrated on selected fields for which linear and non-linear fits were more reliable than on average. Within this subset, no clear distance dependence was visible in k values, but the scatter is larger than the estimated uncertainties. The linear slopes are sensitive to the highest column densities within a field. The parameter C of the non-linear fits is expected to be even more sensitive to these data points, which probe the variations of τ(250 μm) /τ(J) inside the fields. For the whole sample, the median value of C was very close to zero. Nevertheless, it may be significant that when the two fields at kiloparsec distances are excluded, the value of C appears to decrease somewhat systematically with the distance of the field. There is some preference for the most nearby fields to have positive values; here, the densest parts of the clouds are resolved.

Two of the first fields with high k values are G6.03+36.73 and G4.18+35.79, better known as LDN 183 and LDN 134. At the Herschel scale, the clouds have a simple morphology, each consisting of a single column density maximum that is partially resolved by the 180″ beam. This is illustrated by Figs. 10 and E.1. The ratio k = τ(250 μm) /τ(J) closely follows the morphology of the column density distribution, making these the best examples of clumps with increased submillimetre opacity. The highest values are ~ 4 × 10^-3, almost three times the average value over all fields.

In our results, some of the main sources of uncertainty are the corrections made for the expected bias in τ(J) and, to lesser extent, in τ(250 μm). According to Fig. 9, the net effect of all bias corrections is a ~20% decrease in the value of k. This number applies to nearby fields but is more dramatic if all fields are taken into account. Figure 5 shows that for fields at ~1 kpc distance the correction is almost a factor of two. The corrections have been remarkably successful in decreasing the scatter of τ(250 μm) /τ(J) values, which strongly suggests that they are approximately of the correct magnitude.

In the statistical sense, estimating the τ(J) bias is straightforward, using a model of the Galactic stellar distribution and the higher resolution Herschel data as a template for the column density structure of the field (see Appendix B and Sect. 3.7). The correction is large, but on average, the correction itself probably does not suffer from major systematic errors. Errors could arise either from incorrect distance estimates or from the use of an incorrect model of the NIR opacity distribution. The distances are uncertain, and through the τ(J) bias, their effect on the main parameters is shown in Fig. 9 (the grey bands) and in Fig. 8. For a sample of fields, this is mainly a statistical and not a systematic error.

The model of τ(J) distribution in each field was derived from Herschel data at 18″ resolution. If there is still significant structure below this scale, our correction of τ(J) values will be too low. The difference of the 2−3′ scale and the 18″ scale is so large, however, that most of the effects of column density variations are already included. Thus, after the distance, the other main error in the τ(J) bias correction arises from scaling the Herschel estimates of τ(250 μm) into a template of the J-band opacity. We showed in Sect. 3.7 that this amounts to an uncertainty of ~10% in the final listed values of k = Δτ(250 μm)/Δτ(J). For the sample in Fig. 9, the value of k assumed during the bias correction is consistent with the recovered value (to be precise, the assumed value of k = 1.5 × 10^-3, the recovered value of k = 1.6 × 10^-3). For individual fields, the values are biased towards the initial assumed value. Figure 13 showed that if bias correction assumed a value of k = 1.0 × 10^-3, the recovered median value was still higher than k = 1.2 × 10^-3. This shows that k cannot be significantly overestimated because of an erroneous bias correction in τ(J). Thus, the result that the average ratio k = Δτ(250 μm)/Δτ(J) is clearly higher than the normal value found in diffuse medium remains robust.

Because of the bias in the τ(J) correction, our estimates are probably conservative for fields where values k> 1.5 × 10^-3 were obtained. The dependence on the assumed value of k decreases as the true value of k increases. This is because higher k means a lower NIR opacity and lower overall bias in τ(J). Nevertheless, for fields like G4.18+35.79, we might systematically underestimate k by ~10% because of this error in the τ(J) bias correction. At the highest column densities Herschel emission data themselves will underestimate the cloud opacity, which leads to the corrections discussed in Sect. 3.3. For the optical depth ratio, the effect of τ(250 μm) bias is weaker than that of τ(J). Nevertheless, the errors made in the corrections of τ(250 μm) and τ(J) (for example, those associated with a change of submillimetre opacity) would partly cancel each other out.

In the analysis we assumed that apart from the problems associated with the sampling provided by the background stars, NIR reddening is an independent and reliable measure of column density. Unlike in the optical range, the NIR extinction curve is often assumed to be constant for a wide range of column densities (Cardelli et al. 1989; Martin & Whittet 1990; Roy et al. 2013). Nevertheless, some cloud-to-cloud variations are observed, the ratio E_J−H/E_H−K ranging from values lower than 1.5 to higher than 2.0 (Racca et al. 2002; Draine 2003a). Clear changes take place at high optical depths, above A_V ~ 20 mag, but only in the form of the flattening of the MIR extinction curve, at wavelengths above 3 μm (Indebetouw et al. 2005; Cambrésy et al. 2011; Ascenso et al. 2013). Recently, Whittet et al. (2013) studied cloud LDN 183 (which is also included in our sample). Comparison with the 9.7 μm silicate absorption feature suggested that the NIR colour excess might not be a perfect tracer of the total dust column. The ratio between E_J−K and the 9.7 μm feature was observed to increase as the column density exceeded A_V ~ 20 mag. The observation might partly arise because the 9.7 μm feature is dampened by the formation of ice mantles (see also Chiar et al. 2007). However, if we assume that NIR extinction indeed overestimates the total dust column in this object, the ratio of τ(250 μm) relative to column density would increase by ~20% (see Whittet et al. 2013, Fig. 12). This is still a weaker effect than the increase by more than a factor of two that was observed in τ(250 μm) /τ(J). We recall that at high column densities, above A_V ~ 20 mag, the values of τ(250 μm) are also uncertain and can be underestimated by a significant fraction (Pagani et al. 2004 discussed a similar limitation at the nearby wavelength of 200 μm).

Finally, we recall that τ(250 μm) values were derived using a fixed value of the spectral index, β = 2.0. By assuming β = 1.8 instead, the τ(250 μm) and k values would decrease by up to 30%. Conversely, if the value of β increased towards dense and cold regions (Dupac et al. 2003; Désert et al. 2008; Planck Collaboration XXIII 2011; Paradis et al. 2010; Veneziani et al. 2010; Juvela et al. 2011), we underestimate the dust opacity changes if we use a constant value of β. In this sense (and regarding possible bias in τ(J) corrections), Figs. 6–8 give conservative estimates of the possible increase of submillimetre dust opacity. Clearly, the assumed values of β must also be taken into account when comparing our results with other studies (see Sect. 4.2). On the other hand, our tests indicate that spatial variations of β probably do not have a strong additional effect on k (see Sect. 3.7).

4.2. Comparison with other studies

The submillimetre dust opacity has previously been studied in relation to both dust extinction (Terebey et al. 2009; Flagey et al. 2009; Juvela et al. 2011; Martin et al. 2012; Roy et al. 2013; Malinen et al. 2013, 2014) and to HI (Boulanger et al. 1996; Lagache et al. 1999; Planck Collaboration XXIV 2011; Planck Collaboration XXV 2011; Martin et al. 2012). Boulanger et al. (1996) compared COBE observations of FIR dust emission with the hydrogen 21 cm observations. At high latitudes, where HI is a good tracer of the full gas column density, the derived value was τ/N_H = 1.0 × 10^-25 (λ/ 250 μm) cm² H^-1. Similar values were obtained in the first studies using Planck data (Planck Collaboration XXIV 2011).

In the most recent Planck papers, somewhat lower values were reported, corresponding to τ(250 μm) /N_H ~ 0.55 × 10^-25 cm² H^-1 (Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014). The change is associated with the revised calibration of the highest frequency channels and, correspondingly, a lower value of β. In Planck Collaboration Int. XVII (2014) the spectral index above 353 GHz was found to be β ~ 1.65, lower than the value of 1.8 assumed in Planck Collaboration XXIV (2011) or the value of 2.0 used in Boulanger et al. (1996). The higher dust opacity value of Boulanger et al. (1996) is thus largely explained by the higher value of β.

The Planck result for the diffuse medium can be compared to our result most directly by converting their N_H to τ(J). We can use the conversion factor N(H₂) /A_V = 9.4 × 10²⁰ cm^-2 mag^-1 derived by Bohlin et al. (1978) at low extinction, E(B−V) < 0.5 (see also Nozawa & Fukugita 2013). Some studies (Rachford et al. 2009; Planck Collaboration XI 2014; Liszt 2014) have found N_H/E(B−V) values that are 10–30% higher than in Bohlin et al. (1978), these results apply partly to even more diffuse lines of sight (E(B−V) ≲ 0.1). On the other hand, Gudennavar et al. (2012) examined a sample of lines of sight with E(B−V) extending to values higher than one. The result, N(H) /E(B−V) = (6.1 ± 0.2) × 10²¹ H cm^-2 mag^-1, was close to that of Bohlin et al. (1978). With the Bohlin et al. (1978) relation and R_V = 3.1 extinction curve (Cardelli et al. 1989), the Planck result τ(250 μm) /N_H ~ 0.55 × 10^-25 cm² H^-1Planck Collaboration XI (2014) corresponds to τ(250 μm) /τ(J) = 0.41 × 10^-3. Our median value of τ(250 μm) /τ(J) = 1.6 × 10^-3 is thus 3.9 times higher, and our highest local values, 4.5 × 10^-3 in G4.18+35.79, are higher by one order of magnitude. In Planck Collaboration XI (2014) the estimated value of N(H /E(B−V)) was ~20% higher than in Bohlin et al. (1978). With this value, our median value of τ(250 μm) /τ(J) would be ~ 3.3 times higher than the Planck estimate.

We can alternatively convert our result into τ(250 μm) /N_H, but in dense regions the shape of the extinction curve (dependence on R_V) and the ratio between visual extinction and total column density are more uncertain. Using the R_V = 3.1 and N(H₂) /A_V from Bohlin et al. (1978), our median value corresponds to τ(250 μm) /N_H = 2.16 × 10^-25 cm² H^-1 (again, of course, 3.9 times the Planck value). However, R_V is expected to be higher than 3.1 in dense clouds. Using R_V = 5.5 instead of R_V = 3.1 in converting τ(J) into visual extinction, the values of τ(250 μm) /N_H would decrease by ~ 15% (change in A_V/E(J−K)). However, the scaling we used above (corresponding to the R_V=3.1 extinction curve and the ratio of N(H₂) /A_V taken from Bohlin et al. 1978) corresponds to N(H) /E(J−K) = 11.0 × 10²¹ cm^-2 mag^-1, which is very close to the value of 11.5 × 10²¹ cm^-2 mag^-1 that Martin et al. (2012) derived in Vela molecular cloud using 2MASS NIR data up to E(J−K) ~ 0.55 mag (E(B−V) ~ 1.1 mag).

When R_V is modified, the NIR extinction curve remains practically unchanged and the differences take place between the optical and NIR wavelengths. Compared to shorter (optical) and longer (MIR) wavelengths, the extinction curve is considered relatively constant in the NIR regime (e.g. Indebetouw et al. 2005; Lombardi et al. 2006; Román-Zúñiga et al. 2007; Ascenso et al. 2013; Wang & Jiang 2014) and the power-law slope of the NIR extinction curve typically only varies at a level of ~5% (Stead & Hoare 2009; Fritz et al. 2011). In very dense cores the formation of ice mantles and the grain growth could have an additional impact. Román-Zúñiga et al. (2007) examined the cloud Barnard 59 up to A_V = 59 mag, but found no significant changes in the NIR extinction law (compatible with R_V = 5.5). Similarly, Lombardi et al. (2006) found no changes in the Pipe nebula in the reddening NIR law up to E(H−K) ~ 1.5 mag (A_V ~ 9 mag). Therefore, it is not likely that our results are significantly biased because of the assumed NIR extinction curve. At the resolution of our extinction maps, practically all our data points are lower than A_V ~ 10 mag and are thus not affected by extreme optical depths. However, Whittet et al. (2013) observed a decrease in L183 in the ratio of E(J−K) and the 9.7 μm silicate absorption feature. Above E(J−K) ~ 1.0 (A_J ~ 1.7 mag, A_V ~ 5.0 mag for R_V = 5.0), the ratio deviated from diffuse medium value by ~ 20%. Stronger deviations were only seen beyond E(J−K) ~ 3 (A_J ~ 5 mag, A_V ~ 15 mag for R_V = 5.0). This is above the range probed by our measurements. It is not clear that changes in this ratio would only be caused by the NIR extinction curve. However, if, in some sources, the extinction curve does flatten above A(J) ~ 2, this might contribute to the observed increase of τ(250 μm) /τ(J) (possibly at the 20% level).

The effect of the bias corrections on τ(250 μm) /τ(J) is typically ~20% or weaker, and their uncertainty is smaller. The largest uncertainties in τ(250 μm) /N_H are caused by β and possibly by R_V. With β = 1.8 and R_V = 5.5, the τ(250 μm) /τ(J) values would be 40% lower than with β = 2.0 and R_V = 3.1. At this lower limit, our median value of τ(250 μm) /τ(J) would not be 3.9 times, but ~2.3 times the value found in the diffuse high-latitude sky. The 40% uncertainty may be a realistic 1σ lower limit for the linear fits that include all pixels in the maps. However, it is a conservative estimate for the clumps where the average value of β is expected to be clearly higher than 1.8. In fact, preliminary results indicate that the average value of β in our fields (including the more diffuse regions) is close to β = 1.9, and that with β = 2.0 we underestimate the τ(250 μm) values of many clumps (Juvela, in prep.).

Increased FIR and submillimetre opacity has been reported by many authors (Kramer et al. 2003; Lehtinen et al. 2004; del Burgo & Laureijs 2005; Ridderstad & Juvela 2010; Bernard et al. 2010; Suutarinen et al. 2013, etc.). One famous example is the Taurus filament L1506, for which the models of Stepnik et al. (2003) suggested an increase by more than a factor of three. With the recent detailed modelling of the dust properties and the structure of this filament, Ysard et al. (2013) estimated the increase of the 250 μm opacity to be ~2. Martin et al. (2012) studied the relation in the Vela cloud, comparing BLAST and IRAS data with the reddening of 2MASS stars. The properties of the examined areas corresponded to the average properties of our fields, with column densities extending to 10²² cm^-2 and with typical dust temperatures of ~ 15 K. They found a very similar range of dust opacities, τ(250 μm) /N_H = (2−4) × 10^-25 cm² H^-1 (assuming β = 1.8). In the Orion A cloud, a comparison of Herschel and 2MASS data led to the detection of a dependence of N^0.28 on the 250 μm opacity (Roy et al. 2013). The range of column densities in Orion A was similar to our fields, and the derived dust opacities were mainly in the range of τ(250 μm) /N_H = (1−3) × 10^-25 cm² H^-1. These estimates were derived with β = 1.8 and would become ~ 20% higher (depending on the details of the fitting) if a value of β = 2.0 were used. More recently, Lombardi et al. (2014) derived ratios A(K) /τ(850 μm) of 2640 mag and 3460 mag for the Orion A and B molecular clouds, respectively. The 850 μm optical depth was derived from Herschel observations rescaled using a comparison with Planck measurements, and the NIR extinction was calculated with the method NICEST (Lombardi 2009). The modified blackbody fits used β values that were estimated in Planck Collaboration XI (2014) at a resolution of 30′. With the reported average value of β = 1.8 and adopting a ratio 0.40 between A(K) and A(J), the results correspond to τ(250 μm) /τ(J) values of 1.56 × 10^-3 and 1.19 × 10^-3 for Orion A and B, respectively. These fits were made for data τ(850 μm) < 2 × 10^-4 (τ(250 μm) ≲ 1.8 × 10^-3). The optical depth range is similar to many of our fields, and the result for Orion A is close to our median value for the fits concerning entire fields. By adopting β = 1.8, our median value would fall between the Orion A and Orion B estimates of Lombardi et al. (2014).

4.3. Implications for dust evolution

As shown by simulations, an observed increase of the dust FIR/submm opacity towards dense regions cannot be due to radiative transfer effects, but must originate in intrinsic variations of dust properties (Malinen et al. 2011; Juvela & Ysard 2012; Ysard et al. 2012). Theoretical studies have shown that such an increase, coupled with a decrease in dust temperature, can be explained by the formation of large aggregate particles (Ossenkopf & Henning 1994; Stognienko et al. 1995; Ormel et al. 2011; Köhler et al. 2011, 2012). Köhler et al. (2012) showed that the coagulation of just four big grains of the diffuse ISM type, coated by smaller carbon grains, already leads to an increase in the opacity at 250 μm of a factor 2.6. These authors also showed that these aggregates can form within a typical cloud lifetime of 10 million years (Walmsley 1991). Consequently, we interpret the observed increase in τ(250 μm) in our cold core sample as grain growth in dense molecular regions.

In our sample, the clouds LDN 183 and LDN 134 (G6.03+36.73 and G4.18+35.79) were the most convincing examples of increased submillimetre dust opacity. LDN 183 has been studied thoroughly in both continuum and line emission (e.g. Juvela et al. 2002; Pagani et al. 2005, 2007). A slight increase of 200 μm opacity was already reported based on ISOPHOT observations (Juvela et al. 2002), but ISOPHOT data λ ≤ 200 μm and, to some extent even Herschel observations, are not sufficient to fully probe the inner parts of the cloud where the high visual extinction is approaching 100 mag and the dust temperature drops well below 10 K (Pagani et al. 2004, 2015). LDN 183 was the first object where enhanced MIR light scattering, the so-called coreshine phenomenon, was detected in Spitzer data (Steinacker et al. 2010; Pagani et al. 2010). The effect was also seen in Juvela et al. (2012), where WISE MIR observations were analysed and both LDN 183 and LDN 134 were found to be sources of strong MIR emission. If the signal is interpreted as scattering of the interstellar radiation field, it seems to imply the presence of very large, micrometre-sized dust particles (Steinacker et al. 2010, 2014; Lefèvre et al. 2014). Thus, our detection of increased submillimetre opacity in these clouds agrees with the evidence provided by MIR wavelengths.

5. Conclusions

We have examined dust optical depths by comparing measurements of submillimetre dust emission and the reddening of the light of background stars in the NIR. The goal was to measure the value of dust submillimetre opacity and to search for variations that might be correlated with the physical state and the environment of the cloud. The study led to the following conclusions:

• For a subsample of 23 fields with well-defined correlationbetween the two variables, we obtained a median ratio ofτ(250 μm) /τ(J) = (1.6 ± 0.2) × 10^-3. This is more than three times the value that was derived from Planck data for the diffuse medium at high Galactic latitudes. Assuming β = 1.8 instead of β = 2.0, the value decreases by 30%, but is still more than twice the diffuse value.

• The conversion to τ(250 μm) /N_H involves more assumptions. Using the R_V = 3.1 extinction curve and the N(H₂) /A_V ratio of Bohlin et al. (1978), our median estimate corresponds to τ(250 μm) /N_H = 2.16 × 10^-25 cm² H^-1.

• The fit to all data above τ(J) = 1 gives a relation τ(250 μm) ~ 1.25 × 10^-3τ(J) + 0.11 × 10^-3τ(J)². The positive second-order coefficient C = 0.11 × 10^-3 is determined by a small number of fields that, because of the high column density, are subject to large uncertainty in the bias corrections.

• For the same sample, the scatter in the coefficients of the second-order terms C is ~ 2 × 10^-4 and the median value is consistent with zero. Spatial variations of the ratio τ(250 μm) /τ(J) are only seen in a few fields.

• From the maps of τ(250 μm) /τ(J), we identified six fields where the ratio appears to increase further at the location of the main column density peaks. The highest values in the fields G4.18+35.79 and G6.03+36.73 are ~ 4 × 10^-3 at the resolution of 180″. Thus, although the densest clumps are associated with the largest uncertainties, we consider this increase of submillimetre opacity to be real.

Online material

Appendix A: Alternative data analysis

To investigate the robustness of the results to details of the data reduction, we calculated a set of alternative τ(250 μm) maps. In addition to the analysis of three or four Herschel bands (see Sect. 2.2), we considered the use of local background subtraction and the possibility of making a correction for residual mapping artefacts with the help of other all-sky surveys. The resulting τ(250 μm) /τ(J) ratios were compared to the values found in Sect. 3.1. The comparison was carried out without bias corrections, comparing the results with the 160−500 μm fits of Sect. 2.4.

Appendix A.1: Local background subtraction

Our default analysis is based on Herschel data for which the absolute zero points were derived from a comparison with Planck and IRAS maps. As an alternative, we used Herschel surface brightness maps from which the diffuse background was subtracted using the reference regions listed in Table 1. Thus, the average surface brightness of the reference region was subtracted from each Herschel surface brightness map separately, before calculating the colour temperatures and the dust optical depths.

The local background subtraction might be a more reliable way to ensure a consistent zero level for the compared quantities. However, it also means that colour temperature and column density can only be estimated in the part of the map in which the surface brightness values are significantly higher than those of the reference area. We masked the area in which the signal is lower than twice the estimated statistical uncertainty of the surface brightness in each band. The final mask is a combination of these masks and the original mask that eliminated the map boundaries for which the convolution to the resolution of the τ(J) data is only poorly defined.

Appendix A.2: Checks for mapping artefacts

Although Herschel data are usually of very good quality, there can still be some small artefacts that affect some parts of the maps. Errors might result from data reduction or from instrumental effects such as striping or general gain changes (Xu et al. 2014; Paladini et al. 2013). If processing includes high-pass filtering, the large-scale surface brightness gradients may be affected and the contrast between faint and bright regions may be decreased. Our maps often contain significant emission up to the map boundary. Without a flat border region with very low emission, it is difficult to estimate whether the baseline assumed for the scans is correct. Such effects could be more important and more difficult to detect for small maps. Thus, this could mostly affect PACS maps, for which the signal-to-noise ratio is also typically lower than in the SPIRE data (Juvela et al. 2010).

These effects were investigated with the help of independent FIR and submillimetre data. At 100 μm, 350 μm, and 500 μm we can compare Herschel data almost directly with the corresponding IRIS and Planck bands. The Planck 545 GHz data were corrected to 500 μm using a modified black body with the Herschel colour temperature map and β equal to 2.0. Keeping the other parameter constant, an error of 2 K in temperature or an error of 0.2 in β would both correspond to only a ~2% error in the extrapolated value. At 160 μm and 250 μm we used values interpolated from IRIS 100 μm, AKARI 140 μm (Murakami et al. 2007), and Planck 350 μm channels. Here Δβ ~ 0.2 translates into a change of less than 1% at 160 μm. For a direct extrapolation from the 350 μm to 160 μm, an error of Δβ = 0.1 would result in an error in the interpolated value that is still lower than 8%. Typically interpolation errors should thus be below the statistical errors. The same considerations apply to the zero-point corrections of Sect 2.2.2, with the difference that they are not affected by any multiplicative errors.

The reference data were compared with the original Herschel maps at 6′ resolution to derive an additive correction that leaves the median value of the maps unchanged and only affects scales larger than ~ 6′. We also checked similar multiplicative corrections, assuming that the zero points of the different surveys are compatible, and even calculating corrections where linear fits between Herschel and reference data were first used to estimate the differences in zero-point and gain calibration. The last alternative would avoid the assumptions of consistent calibration and zero points between the data sets. In most cases, there are no significant differences between the three choices. A typical map has no clear artefacts, and the proposed correction will probably decrease the data quality. However, when the local artefacts are clear (e.g., excessive surface brightness some corner of a Herschel map), the corrected map should give a better description of the true surface brightness. Thus, we do not believe that the corrected maps represent a clear improvement, but the difference between the corrected and uncorrected data should give some idea of the potential effect that artefacts of that magnitude could have.

Appendix A.3: Comparison of the data sets

We compared the τ(250 μm) /τ(J) values (without bias corrections) among six data sets: (1) three bands at 250-500 μm (our default data set); (2) four bands at 160–500 μm; (3) three bands with local background subtraction; (4) three bands with the corrections of Sect. A.2; and (5) four bands with the corrections of Sect. A.2. Table A.1 shows the results, comparing the dispersion of the obtained τ(250 μm) /τ(J) values between fields and the change in the values compared to the default case where three bands and the absolute zero points were used.

The second column of the table lists the number of fields where the formal error of the slope τ(250 μm) /τ(J) is lower than 10%. This number is smaller when PACS data are included, 82–83 fields compared to the 102−105 fields with SPIRE data alone. The numbers do not include the Witch Head Nebula, for which we have no PACS data and which therefore was excluded from this comparison. The background subtraction and the Sect. A.2 corrections both decrease the mean value, but not significantly. Note that on physical grounds one could have expected the values to increase with background subtraction (if dense material has higher τ(250 μm) /τ(J)) and to decrease with the inclusion of 160 μm channel (if the inclusion of shorter wavelengths increases estimated colour temperatures). The last column shows the difference relative to our default case. In this column we only list the common set of 81 fields where the error estimates are lower than 10% for all the five cases. On average, the changes in the τ(250 μm) /τ(J) values are less than 1.0 units (lower than 5%). The dispersion between different SPIRE analyses is lower than 2.0 units and somewhat higher when analyses of four and three bands are compared.

Appendix B: Simulated NIR extinction maps

In addition to photometric errors, the reliability of NIR optical depth estimates is mainly affected by the sampling provided by the background stars and the possible contamination by foreground stars and galaxies.

We used 2MASS catalogue flags to eliminate most of the obvious galaxies. In addition to requiring a photometric quality corresponding to ph_qual in classes A–C, we excluded all point sources that were extended (flag ext_key is set) or were flagged with gal_contam. These only remove part of the galaxies. The increased dispersion of intrinsic colours caused by galaxies is taken into account in the error estimates provided by the method NICER. Because our simulations use actual 2MASS data near the target fields, this extragalactic contamination is also automatically present in the simulations described below.

The simulations are based on dust 250 μm optical depth maps derived from Herschel observations. Using only SPIRE data, we derive column densities at 25″ resolution as a combination of $\mathit{\tau }\mathrm{=}\mathit{\tau }\mathrm{\left(}\mathrm{500}\mathrm{\right)}\mathrm{+}\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathrm{350}\mathrm{)}\mathrm{-}\mathit{\tau }\mathrm{(}\mathrm{350}\mathrm{\to }\mathrm{500}\mathrm{\left)}\mathrm{\right]}\mathit{,}$(B.1)where τ(500) is calculated using 250, 350, and 500 μm maps at the lowest common resolution, τ(350) is calculated similarly from 250 μm and 350 μm maps, and τ(350 → 500) is the latter convolved to the resolution of the 500 μm observations (see Palmeirim et al. 2013). Thus, the expression in square brackets describes structures that are seen at the resolution of 350 μm data (25″), but not at the resolution of 500 μm data (36″). The calculations assume a fixed value β = 2.0. The differences to the τ(250 μm) maps used in Sect. 3 are small and, furthermore, we only consider differences between these input maps and the values recovered by NICER. On the other hand, we wish to retain the highest resolution possible (18″ instead of 36″) because the bias in NICER estimates is probably linked to the amount of small-scale structure.

The Besancon model (Robin et al. 2003) was used to create a simulated catalogue of stars over a 0.5° × 0.5° area centred on each target field. The catalogue includes stellar distances and H-band magnitudes, and together with the distance estimates listed in Table 1, this was converted into the probability that a star of given magnitude resides between the cloud and the observer. In the simulation, the corresponding fraction of stars was assumed to be located in front of the cloud.

To simulate NIR observations, we used the same 2MASS data that were used to calculate the actual τ(J) maps of the fields. The stars in the reference area (see Table E.1) were used to determine an empirical probability distribution of H-band magnitudes and the dependence between the J, H, and Ks magnitudes and their uncertainties, as given in the 2MASS catalogue. This reference area may be affected by small amounts of absorption by diffuse dust, but gives a good approximation of the brightness distribution of stars that are unextincted by the main cloud.

We simulated a uniform distribution of stars over the Herschel field, generating the magnitudes from the empirical H-band magnitude distribution and matching the average stellar density of the reference region. The J and Ks magnitudes of each star were generated using the J−H and H−Ks colours of a random star selected from the reference region. This ensures that the distribution of intrinsic colours is realistic and that the simulations reproduce proper correlations (also in errors) between the bands.

Based on the Besancon model, a fraction of stars was assumed to reside in front of the cloud and to be unaffected by extinction. For the remaining stars, the line-of-sight optical depth in J band was calculated using the input column density map and a fixed ratio of τ(250 μm) /τ(J) = 1.5 × 10^-3. The optical depths in H and Ks bands then follow from the Cardelli et al. (1989) extinction curve. The magnitudes were adjusted according to the line-of-sight optical depths. The magnitudes and their uncertainties in the reference area were used to calculate the typical photometric uncertainties as a function of magnitude, and they were used in the NICER calculation. Because the intrinsic colours of the stars were generated based on observed stars, no additional scatter needed to be added for intrinsic colours. However, because the extinction makes the stars fainter, the typical photometric errors increase as well. This was taken into account by adding normal distributed noise to the magnitudes, which corresponds to the difference in the typical uncertainties between the original and the extincted magnitudes.

The simulated stellar catalogues were fed to the NICER algorithm to derive extinction maps with the same parameters as in Sect. 2.3. For each target field, one hundred realisations of the τ(J) maps were calculated to obtain maps for the standard deviation and the bias of the estimated τ(J) values. Figure B.1 shows one example.

Appendix C: Simulated Herschel observations

The estimation of the τ(250 μm) bias is more uncertain than the estimation of τ(J) bias. Because of line-of-sight temperature variations, the colour temperature overestimates the mass-averaged dust temperature, which translates into too low estimates of τ(250 μm). The magnitude of the effect cannot be estimated precisely because the line-of-sight temperature distribution is unknown. Order of magnitude estimates can be obtained with radiative transfer modelling by making assumptions of the radiation field, dust properties, and the cloud structure.

We carried out radiative transfer calculations individually for all the 116 fields. We assumed constant dust properties corresponding to Milky Way dust with a selective extinction of R_V = 5.5 (Draine 2003b). The initial radiation field was assumed to correspond to the Mathis et al. (1983) model of solar neighbourhood. The spectrum of the illuminating radiation has an impact on the temperature contrasts. We have no way to independently determine the shape of the spectrum of the radiation field. However, this typically remains a second-order effect, and the main effect, the level of the radiation field intensity, is part of the modelling. One exception are the possible stars that heat clouds from the inside and may locally have a strong effect. These are considered later in the analysis, but are not part of the simulations.

For each field, we built a model that attempts to reproduce the 250–500 μm observations of that field. From the background-subtracted surface brightness maps we selected the central 30′ × 30′ area. The model cloud had the same angular dimensions and was discretised onto a 181³ cell grid. Each cell of the model therefore corresponds to an angular size of 10″, but the linear scale depends on the distance of the cloud. In the line-of-sight direction we assumed a Gaussian density distribution with a FWHM equal to 25% of the field size. The linear size again depends on the cloud distance because in more distant fields we are also probably concerned with larger structures. In the line-of-sight direction the density peak is always in the central plane of the model cube. This increases mutual shadowing of dense regions, which in reality can reside at different distances and increases the temperature contrasts in the model. On the other hand, for a field at 200 pc distance, the selected line-of-sight FWHM extent is ~0.4 pc, which is larger than the size of typical cores. Therefore, it is difficult to say whether the selection of this particular line-of-sight density profile leads to an over- or underestimation of the final τ(250 μm) bias. These are usually second-order effects (Juvela et al. 2013) except for very dense clumps that can remain practically invisible in τ(J) maps as well.

We carried out radiative transfer calculations to produce synthetic surface brightness maps at Herschel wavelengths that were then convolved to the resolution of the observations. The ratio of observed and modelled 350 μm maps was used to adjust the column densities, applying the same multiplicative factor to all cells along the same line-of-sight. The intensity of the external radiation field was scaled based on the ratios between the observed and modelled 250 μm and 500 μm surface brightness. The aim is to also reproduce the average shape of the spectra. The full procedure was iterated until the model matched the observed 350 μm map at ~1% accuracy and the average ratio 250 μm/500 μm were correct within the same tolerance.

The final model takes into account the range of column density, the morphology, and the radiation field intensity of a field. It is not necessarily a perfect match to all surface brightness maps, but is a good facsimile of the pixel-by-pixel column density structure. We analysed the synthetic surface brightness maps as in Sect. 2.2.3. The comparison of the obtained τ(250 μm) estimates and the true values known for the model cloud gives a 30′ × 30′ map of the τ(250 μm) bias (resolution 40″). The remaining border areas usually have a low column density and therefore low bias. However, to extend bias estimates over the whole map area, we assigned values calculated using the average bias vs. column density relation estimated from central 30′ × 30′ area to the remaining pixels.

Figure C.1 shows one example, the surface brightness maps produced by the model and the bias in the τ(250 μm) values estimated from these synthetic observations.

Appendix D: Additional checks of correlations between τ(250 μm) and τ(J)

In addition to the factors examined in Sect. 2.4, we checked the importance of two additional factors, the technical implementation of the least-squares fits, and the importance of internally heated regions.

The total least-squares fits of Sect. 2.4 used the formal error estimates of τ(J) and τ(250 μm). The former were obtained from NICER routine, the latter were estimated with MCMC calculations starting with the assumption of 7% (SPIRE) or 15% (PACS) relative errors in the surface brightness data. If the correlation is poor, the estimated slope becomes sensitive to the error estimates. For example, if the true errors of τ(250 μm) were much larger, for example because of some artefacts in map making, the use of too low error estimates would increase the slope estimates. We checked this by repeating the analysis using twice the original τ(250 μm) error estimates. In an extreme case, we can ignore the error estimates altogether and perform unweighted least-squares fits. Based on the error estimates used, the true relative uncertainty is significantly larger in τ(J) than in τ(250 μm). Therefore, the unweighted least-squares fit probably underestimates the true slope.

The third test concerns the internally heated regions in which because of strong temperature variations combined with compact, high column density clumps, both optical depth estimates are particularly uncertain. Furthermore, the estimates of τ(250 μm) bias are probably incorrect in the same areas. This is caused by two factors. First, without the internal heating source,

the models are unable to produce sufficient surface brightness values and result in very high column densities and thus high estimates of the bias (Juvela et al. 2013). Second, in the real clump the internal heating may help to decrease the actual bias if the clump centre remains warm instead of being too cold to be registered in Herschel bands (Malinen et al. 2011). We repeated the analysis of Sect. 2.4 by masking warm regions. We first masked all pixels for which the dust colour temperature was higher than 20 K. The mask was then extended to cover areas in which after convolution to 180″ resolution, the influence of the T> 20 K region was more than 10% of the convolved value.

Figure D.1 compares the result with the bias-corrected results already shown in Fig. 4. It shows that the shape of the distribution is not sensitive to the fitting procedure, nor is it significantly affected by the warm regions.

Appendix E: Maps of τ(250 μm) /τ(J) ratio

In the least-squares fits of Sect. 3.4, the interesting parameters, k and C were independent of additive errors in the correlated quantities. The highest τ(J) points were particularly important, both visually and regarding the fitted parameters. As mentioned in Sect. 3.5, we also calculated maps of the τ(250 μm) /τ(J) ratio. Because the column densities are typically low over most of the mapped area, the visual appearance of the maps is dominated by regions with low τ(250 μm) and τ(J) values where, by definition, the results become sensitive to any zero-point mismatch.

The maps were calculated by first correcting the τ(250 μm) and τ(J) maps for the bias that was estimated with modelling (see Sect. B and Sect. C). The maps of τ(250 μm) were then convolved to the resolution of the τ(J) map. A 2′ wide region near the Herschel map borders was masked because the convolved values would be affected by data outside our map coverage. Before calculating the ratio τ(250 μm) /τ(J), we subtracted from both quantities the values in the reference regions that are listed in Table 1. A typical diameter of these reference areas is ~ 6′. For τ(250 μm) the statistical error of the reference value is very small. The true uncertainty of the remaining τ(250 μm) is completely dominated by systematic errors. Because the main purpose of the ratio maps is to determine variations in the τ(250 μm) /τ(J) ratio, we are not very concerned with multiplicative errors.

We expect the statistical errors to be more significant in τ(J) and, because the variable is in the denominator, we need to mask areas with a low SN of τ(J). NICER has provided error maps for τ(J), but here we estimated the uncertainty using the following procedure: We took the data plotted in Fig. 8, selected 20% of the lowest τ(J) points, and subtracted from them the prediction of the non-linear fit (red line in Fig. 8). The uncertainty of the reference value, Δτ(J), was calculated as the standard deviation of the residuals, scaled by the ratio of (90″)² and the area of the reference region. This should be a very conservative estimate because it assumes that in Fig. 8 the scatter would be due to τ(J) errors alone.

The results are shown in Fig. E.1. The first frames show the τ(250 μm) maps, the contours indicating the region with a SN of each parameter higher than one. The second frames show the calculated maps of τ(250 μm) /τ(J). The regions where either parameter falls below SN = 0.5 were masked. The remaining frames show the extreme cases corresponding to τ(J) ± Δτ(J) and τ(250 μm) ± δτ(250 μm), where δτ(250 μm) is the estimated error map (see Sect. 2.2.3).

---

Acknowledgments

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. M.J., J.M.o., V.M.P., and J.M.a. acknowledge the support of the Academy of Finland Grant No. 250741.

References

All Tables

All Figures