Dissecting the origin of the submillimetre emission in nearby galaxies with Herschel and LABOCA

Abstract

1 INTRODUCTION

Studying dust emission helps us to probe the star formation activity obscured by dust grains and thus to understand how galaxies evolve through time. Many studies using telescopes such as the Infrared Astronomical Satellite (IRAS) or the Spitzer Space Telescope (Spitzer) have been carried out to better understand the different grain populations contributing to the near-to-far-infrared (NIR, FIR) emission of nearby galaxies and investigate their hot and warm dust components. Recent observations of the FIR and submillimetre (submm) using the Herschel Space Observatory (Herschel) allow us to probe, with increased resolution compared to Spitzer, the cold dust emission, study its properties (temperature, grain opacity) within galaxies rather than on integrated scales and analyse how they vary with the interstellar medium (ISM) physical conditions (see for instance the analysis on nearby resolved galaxies of Boquien et al. 2011; Aniano et al. 2012; Bendo et al. 2012a; Foyle et al. 2012; Mentuch Cooper et al. 2012; Smith et al. 2012; Draine et al. 2013; Galametz et al. 2013a).

Submm observations of nearby galaxies with Herschel and from the ground have suggested that the cold dust properties could vary with the ISM physical conditions, for instance with metallicity. Deriving the dust masses of low-metallicity galaxies using standard dust properties for instance often leads to unphysical dust-to-gas mass ratios (D/G), compared to those expected from their lack of metals (Galametz et al. 2010; Meixner et al. 2010; Galliano et al. 2011). An excess at submm wavelengths is also often reported in those objects (Galliano et al. 2003; Dumke, Krause & Wielebinski 2004; Galliano et al. 2005; Bendo et al. 2006; Marleau et al. 2006; Galametz et al. 2009; Bot et al. 2010b; O'Halloran et al. 2010; Rémy-Ruyer et al. 2013). The submm ‘excess’ is usually defined as the excess emission above that predicted from FIR data (excluding the submm measurement) using a single λ^−β emissivity law. (but we will also compare our observations with predictions from a more realistic dust model in the resolved analysis, Section 4). The presence of submm excess in a substantial number of galaxies challenges the models standardly used (more particularly how the submm regime is modelled) and raises the issue of how cold grain properties vary with environment (metallicity, temperature, star formation activity, etc.). The origin of the observed excess emission is still an open question. The various hypotheses investigated so far can be classified in two categories.

(1) On the one hand, the excess could result from an additional dust component not accounted for in the current models. Some studies suggest that this excess could be linked with a reservoir of very cold dust (<15 K) distributed in dense clumps (Galliano et al. 2005; Galametz et al. 2011) or magnetic dipole emission from magnetic grains (Draine & Hensley 2012). Others propose that the excess characterizes the so-called anomalous microwave emission generally associated with ‘spinning dust’ emission (Draine & Lazarian 1998) with a peak frequency which varies with environment or grain properties (Bot et al. 2010a; Murphy et al. 2010; Ysard & Verstraete 2010; Peel et al. 2011; Planck Collaboration et al. 2011b), although Draine & Hensley (2012) have argued that rotational emi-ssion should not significantly contribute to the 870 μm emission.

(2) On the other hand, the excess could be linked with a variation of dust properties with environment, in particular with temperature or density. Analysis of spatially resolved dust properties in the Galaxy suggest, for instance, that the emissivity of dust grains seems to increase towards the coldest regions of molecular clouds (Juvela et al. 2011), supporting the hypothesis of grain coagulation processes in those environments (Paradis, Bernard & Mény 2009). Moreover, the use of amorphous carbon (AC) instead of standard graphite is sometimes preferred to model carbon dust in the spectral energy distribution (SED) fitting of thermal dust emission in dwarf galaxies (Galametz et al. 2010; Meixner et al. 2010; O'Halloran et al. 2010; Galliano et al. 2011). Indeed, AC grains have a lower emissivity index, resulting in a flatter submm spectrum, and thus require less dust to account for the same emission. This usually leads to lower dust masses and seems in those cases to reconcile the dust-to-gas mass ratios with those expected from the low metallicity of these objects. Many more investigations are necessary to link this excess with the ISM conditions and enable us to disentangle its different possible origins.

Even if mostly observed in low-metallicity objects, a submm (850 or 870 μm) excess has also been detected on global scales in solar-metallicity objects (Galametz et al. 2011). Submm excess emission is usually not (or barely) detected at 500 μm in spiral galaxies (Bendo et al. 2010; Dale et al. 2012; Smith et al. 2012; Kirkpatrick et al. 2013), suggesting that observations beyond Herschel are necessary to properly probe the presence of excess in these galaxies. We now investigate the potential excess emission beyond 500 μm in 11 nearby galaxies chosen from the KINGFISH programme (Key Insights on Nearby Galaxies: A Far-Infrared Survey with Herschel; Kennicutt et al. 2011), adding 870 μm measurements taken with LABOCA (on APEX) to the NIR-to-FIR data set. These 11 objects were selected because they are large enough to use the LABOCA mapping mode and their diffuse emission is expected to be bright enough to be detected by the instrument.

Our Spitzer+Herschel+LABOCA data set is ideal to resolve the potential 870 μm excess and understand its origin. This paper complements the work of Galametz et al. (2012) (hereafter [G12]) that studies the spatially resolved cold dust properties (temperature, emissivity index, mass) for the same sample. The paper is structured as follows. In Section 2, we describe the data sets we use. In Section 3, we study how the 870 μm emission compares with various model extrapolations from the Herschel wavebands on global scales. In Section 4, we simulate 870 μm maps using various resolved SED modelling techniques and compare them with our LABOCA observations in order to spatially analyse the differences between these maps. We summarize the analysis in Section 5.

2 OBSERVATIONS AND DATA REDUCTION

This work is based on a combination of Spitzer, Herschel and LABOCA observations available for a sample of 11 nearby galaxies. A description of the sample is available in the first part of this study ([G12]) and general properties are summarized in Table 1. We provide some details of the observations and data reduction processes in this section.

2.1 Spitzer data

Warm dust can contribute a non-negligible amount of the 70 μm emission in galaxies. In order to constrain the mid-infrared (MIR) SED, observations from the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) at 24 and 70 μm are included in the data set. The 11 galaxies have been observed with the MIPS as part of the Spitzer Infrared Nearby Galaxies Survey (SINGS; Kennicutt et al. 2003). The full-width maxima (FWHM) of the MIPS 24 and 70 μm PSFs are 6 arcsec and 18 arcsec, respectively. Some galaxies within the sample are also part of the Local Volume Legacy (LVL) survey. All the data were processed with the latest LVL pipeline¹ (Dale et al. 2009), whether as part of that sample or to make them consistent with it.

2.2 Herschel data

The 11 galaxies of the sample have been observed with Photodetector Array Camera and Spectrometer (PACS) and Spectral and Photometric Imaging Receiver (SPIRE) on-board Herschel as part of the KINGFISH programme. Those images allow us to sample the peak of the thermal dust emission and to probe the cold dust phases. PACS observed at 70, 100 and 160 μm with FWHMs of 5.2 arcsec, 7.7 arcsec and 12 arcsec, respectively (Poglitsch et al. 2010). The data are processed to Level 1 using the Herschel Interactive Processing Environment (hipe) software. Scanamorphos was then used to produce the final PACS and SPIRE maps (Roussel 2013). The scanamorphos algorithm was preferred over the standard madmap mapmaking technique because it can subtract the brightness drifts caused by low-frequency noise using the redundancy built in the observations. PACS maps are calibrated in Jy pixel⁻¹ and have a final pixel size of 1.4 arcsec, 1.7 arcsec and 2.85 arcsec at 70, 100 and 160 μm, respectively. The PACS calibration uncertainties are ∼10 per cent (see the PACS Observer's Manual²). SPIRE produced maps at 250, 350 and 500 μm, with FWHMs of 18 arcsec, 25 arcsec and 36 arcsec, respectively (Griffin et al. 2010)³. SPIRE maps are calibrated in Jy beam⁻¹. The final pixel size is 6 arcsec, 10 arcsec and 14 arcsec at 250, 350 and 500 μm, respectively, with calibration uncertainties of ∼7 per cent for the three wave bands (see the SPIRE Observer's Manual⁴). We estimate the SPIRE beam sizes using the flux calibration paper of Griffin et al. (2013). We apply their factors to convert from the point source pipeline to the extended source pipeline. Characterizing each band by a solid angle appropriate for a spectrum of the form I_ν ∝ ν⁻¹, we derive beam sizes of 469.1, 827.2 and 1779.6 arcsec² for 250, 350 and 500 μm, respectively. Fluxes are then converted in MJy sr⁻¹. More details on the data reduction of Herschel maps are available in the overview paper of the KINGFISH programme (Kennicutt et al. 2011).

2.3 LABOCA 870 |$\boldsymbol\mu$|m maps

In order to search for potential excess emission beyond 500 μm, we complement the NIR-to-FIR coverage of the dust emission with LABOCA observations of our galaxies at 870 μm. LABOCA is a multichannel bolometer array for continuum observations located on the Atacama Pathfinder EXperiment (APEX) telescope in the Atacama desert, Chile (Siringo et al. 2009). The FWHM of the point spread function (PSF) at 870 μm is ∼19.2 arcsec (Weiß et al. 2009), thus close to that of SPIRE 250 μm.

Data are taken from various projects (C-082.F-0005B, E-080.B-3057A, E-082.B-0679A, E-083.B-0813A, E-086.B-0927A, M-079.F-0129, M-081.F-0012, M-083.F-0052 and M-085.F-0058). We refer to Albrecht et al. (in preparation) for details on the iterative data reduction procedure used to produce the LABOCA maps. A summary is provided here. The mapping is performed in spiral mode (raster pattern). Correction for atmospheric attenuation was enabled through determinations of the zenith opacity via skydips (∼ every two hours). The flux calibration (achieved through observations of Mars, Uranus, Neptune and secondary calibrators) is estimated to be accurate within 9 per cent rms. The time-ordered data streams are reduced using the boa software⁵, following steps of calibration, opacity and temperature drifts correction, flat-fielding, conversion to Jy, flagging of dead or bad channels, removal of the correlated noise (on the global array or induced in groups sharing the same electronics), flagging of stationary points or data taken above a given scanning velocity and acceleration threshold, ³He temperature drift correction, median baseline removal and despiking. A final intensity map is obtained by co-adding the individual subscans (pixel size of 4 arcsec × 4 arcsec) along with a weight map. The steps of correlated noise or median baseline removal can lead to a subtraction of real flux in the maps. This can be solved by using an iterative approach. From the first map produced using the technique previously described, a source model is created by isolating the pixels superior to a given signal-to-noise ratio (areas of apparent source emission defined with polygons) and the data points corresponding to the model are excluded from the baseline determination steps while re-running the reduction pipeline. The new map produced is used to build a new model and serves as an input for the next iteration. The process is repeated until the intensity map converges. rms and signal-to-noise maps are finally produced using a jackknife technique.

The Herschel/350 μm maps are shown in Fig. 1 with the signal-to-noise of the LABOCA 870 μm maps overlaid as contours. Before display, the LABOCA maps are convolved to the resolution of the SPIRE 350 μm maps (i.e. 25 arcsec in Fig. 1) using the convolution kernels⁶ developed by G. Aniano (see Aniano et al. 2011, for a description of the technique). We restrict the black and white contours to the apertures used by Dale et al. (2012) to perform the global Herschel photometry of the galaxies. Therefore, the 870 μm contours of NGC 1317, companion of NGC 1316 and also detected with LABOCA are not shown in Fig. 1.

Figure 1.

Herschel/SPIRE 350 μm maps of the sample in MJy sr⁻¹ (colour scales are provided for each object in linear scale), with the signal-to-noise of the LABOCA 870 μm maps overlaid as contours. The signal-to-noise maps have been convolved to a resolution of 25 arcsec (that of SPIRE 350 μm) before display. Contours are, for NGC 0337: 3, 6, 9, 12 and 17σ where σ = 0.17 MJy sr⁻¹, for NGC 0628: 3, 6, 9 and 12σ where σ = 0.19 MJy sr⁻¹, for NGC 1097: 3, 6, 9, 20 and 30σ where σ = 0.13 MJy sr⁻¹, for NGC 1291: 3 and 6σ where σ = 0.13 MJy sr⁻¹, for NGC 1316: 3 and 6σ where σ = 0.20 MJy sr⁻¹. NGC 1316 is the southern object in this map, NGC 1317 the northern object. Contours are, for NGC 1512: 3, 6, 9, 12, 17 and 22σ where σ = 0.09 MJy sr⁻¹, for NGC 3351: 3, 5, 12 and 21σ where σ = 0.16 MJy sr⁻¹, for NGC 3621: 3, 6, 9, 15 and 21σ where σ = 0.16 MJy sr⁻¹, for NGC 3627: 3, 6, 9, 12, 17 and 27σ where σ = 0.22 MJy sr⁻¹, for NGC 4826: 3, 9, 15, 27 and 40σ where σ = 0.22 MJy sr⁻¹, for NGC 7793: 3, 6, 9, 12 and 17σ where σ = 0.19 MJy sr⁻¹.

Open in new tabDownload slide

Figure 1

– continued

Open in new tabDownload slide

3 THE 870 μm EMISSION ON GLOBAL SCALES

3.1 Global flux densities

Along with a detailed description of the LABOCA observations and data reduction, Albrecht et al. (in preparation) provide the global photometry of the sample. Their elliptical apertures are defined to cover the extended emission as observed in the Spitzer/MIPS 160 μm maps.

We compare the integrated LABOCA flux densities with 870 μm estimates extrapolated from observations at shorter IR/submm wavelengths. The global Spitzer and Herschel flux densities are taken from Dale et al. (2007) and Dale et al. (2012), respectively. Their photometric elliptical apertures (provided in [G12]) are chosen to encompass the total emission at each Spitzer and Herschel wavelength. Uncertainties on Spitzer and Herschel fluxes were computed as a combination in quadrature of the calibration uncertainty and the measurement uncertainty. The SPIRE beam sizes have been updated since the publication of Dale et al. (2012) and [G12]. The SPIRE global flux densities are thus rescaled to the new beam size estimates (see Section 2.2 for values), leading to a 7–9 per cent correction compared to [G12]. Discrepancies can be observed between MIPS 70 μm and Herschel/PACS 70 μm fluxes (Aniano et al. 2012) because PACS is less sensitive than MIPS to diffuse emission and because the instruments possess different filter profiles. We thus decide to use both observations in our modelling, to be conservative. Colour corrections will be applied during the SED model fitting. The contribution of C+ to the 160 μm emission is considered to be minor (C ii/L_TIR ∼ 0.15 per cent from Brauher, Dale & Helou 2008 or Malhotra et al. 2001 and L_TIR/νF_ν(160) ∼ 2.5 from Dale et al. 2009 so C ii/νF_ν(160) ∼ 0.4 per cent). The SPIRE flux densities are also not corrected for line contribution but potential line features (such as [N ii]205 μm or CO) are considered to be minor compared to continuum emission.

Comparison with Planck data. We can compare the global Herschel 350 and 500 μm flux densities and the LABOCA 870 μm flux densities with the Planck global flux densities at 350 μm, 550 μm, 850 μm (FWHMs of 4.3 arcmin, 4.7 arcmin, 4.8 arcmin, respectively). The Planck flux densities are provided by the Planck Catalogue of Compact Sources (PCCS Planck Collaboration et al. 2013) and obtained through the IPAC Infrared Science Archive⁷. Most of our sources are resolved by Planck. We choose the flux densities estimated using the APERFLUX technique (fitting of a variable circular aperture centred at the position of the source). The Planck global flux densities at 350 μm, 550 μm, 850 μm and 1.3 mm are overlaid on the global SEDs of the sample (see Section 3.2 for the description of the modelling) in Fig. 2. We observe a very good agreement between the global flux densities in the [350–850] μm range. The mean value of the Planck350/SPIRE350 μm flux density ratios is 0.96 ± 0.19 and the mean value of the Planck850/LABOCA870 μm flux density ratios is 0.87 ± 0.28. The mean value of the Planck550/SPIRE500 μm flux density ratios is 0.64 ± 0.15, but becomes 0.94 ± 0.23 when Planck 550 μm flux densities are converted to 500 μm flux densities using a (550/500)⁴ factor between the two (i.e. assuming dust with a ν² opacity and a Rayleigh–Jeans behaviour). Using the GAUFLUX values from the PCCS (fitting of a variable Gaussian model to the source), ratios are equal to 1.20 ± 0.21 (when converted as above) at 500 μm and 1.07 ± 0.30 at 870 μm, which is also in very good agreement as well. In particular for the galaxy NGC 1512, the GAUFLUX estimates are more consistent with the Herschel and LABOCA data than the APERFLUX estimates. In some galaxies, the dispersion in the Planck-to-SPIRE or Planck-to-LABOCA ratios are relatively high compared to the calibration uncertainties, i.e. 1 per cent at 850 μm and 10 per cent at 350 and 550 μm (Planck Collaboration et al. 2013), 7 per cent for SPIRE (see Section 2.2 for references) and 9 per cent for LABOCA. The dispersions in the Planck/SPIRE ratios is however consistent with what has been reported in other studies (see the work of Herranz et al. 2013, for instance). Some contamination of the Planck fluxes could originate from the large beam sizes of the telescope, in particular for NGC 0337 which is smaller than Planck's beam. This could explain the large differences we observe for this particular galaxy. Differences between the spectral responses of the Planck, LABOCA and SPIRE instruments are not taken into account in this comparison either, which could also explain the dispersions.

Figure 2.

Global SEDs obtained with a two-temperature fitting technique. The model using an emissivity index of the cold component β_c fixed to 2.0 is shown by a dotted blue line, with β_c fixed to 1.5 by a long dashed green line and with a free β_c by a black solid line. The corresponding warm components (with β_w fixed to 2.0 in the three cases) are overlaid in grey. MIPS, PACS and SPIRE data are overlaid with black rectangles. Observed LABOCA flux densities (not included in the fit) are overlaid in red and have been corrected for CO and radio continuum (see Section 3.1). We also indicate, when available, Planck global flux densities at 350 μm, 550 μm, 850 μm and 1.3 mm (grey crosses) but do not include them in the fit. The height of the grey crosses represents the uncertainties of the Planck fluxes.

Open in new tabDownload slide

3.2 Global SED modelling technique

While Klaas et al. (2001) or Dunne & Eales (2001) showed that the SEDs of nearby galaxies could be modelled with either one MBB with a variable β or multiple MBB with β = 2.0, others showed that the bulk of the FIR emission from dust grains could be fitted with dust heated by a single radiation field (Draine et al. 2007). Recently, other studies have shown that dust emitting in the 100-to-500 μm range may be heated by both the star-forming regions and the older stellar populations (Boquien et al. 2011; Bendo, Galliano & Madden 2012b), which questions the single MBB fitting of that full wavelength range. Isothermal or two-temperature fits are not fully able to model the range of grain temperatures observed in each of our resolved ISM elements. We remind the reader that what we call ‘emissivity index’ in the following analysis is thus, in fact, the ‘effective emissivity index’ (i.e. the power-law index resulting from the different grain populations) rather than the ‘intrinsic emissivity index’ (i.e. the power-law index of the MBB modelling the emission of a dust grain at a given temperature). The following analysis thus primarily investigates submm slope trends with environment rather than directly probing the intrinsic emissivity properties of cold grains.

We present the global SEDs obtained using the three different hypotheses: β_c fixed to 2.0 (dotted blue line), fixed to 1.5 (dashed green line) or free (plain black line) in Fig. 2. We overlay the warm dust component for each model in grey, PACS and SPIRE data with black rectangles and the LABOCA 870 μm flux densities corrected for free–free and synchrotron and CO line contribution (taken from Albrecht et al. in preparation) with red rectangles. As mentioned before, we also overlay (when available) the Planck global flux densities at 350 μm, 550 μm, 850 μm as well as the Planck 1.3 mm flux density but did not use them in the fit.

3.3 870 μm predictions and comparison with observations

We deduce global flux estimates at 870 μm from the different 2MBB models and estimate the uncertainties using a Monte Carlo technique: we generate 100 sets of modified 24-to-500 μm constraints, with fluxes varying randomly within the individual error bars (following a normal distribution around their nominal value). The absolute calibration is the dominating source of uncertainty and is highly correlated across the SPIRE bands (SPIRE's observers’ manual). SPIRE measurements are thus randomly modified but in a similar correlated way. We then apply our fitting technique to the 100 modified data sets. The global 870 μm estimates we extrapolate are the medians of these distributions; we take the standard deviation of each distribution as the uncertainty.

Table 2 summarizes for each galaxy the global 870 μm estimate and 1σ uncertainty obtained using our 3 different MBB fitting techniques (β_c = 2.0, 1.5 and β_c free). The ratios between the observed 870 μm flux densities of Albrecht et al. (in preparation), corrected for free–free and synchrotron emission and line contribution (called ‘870 obs corrected’ in the table), and our extrapolations from the different models (β_c = 2.0 or 1.5 and β_c free) are also provided. The update of the SPIRE beam sizes discussed in Section 3.1 does not significantly change the β values derived when this parameter is free in the models compared to the results of [G12] (variations by less than 5 per cent on average).

For 7 galaxies (NGC 1097, NGC 1291, NGC 1316, NGC 3351, NGC 3621, NGC 3627 and NGC 4826), the observed 870 μm flux densities match the values estimated using β_c = 2.0, within the uncertainties. The global 870 μm emission in those objects could thus be explained by thermal emission modelled using dust grains with standard properties. However, we note that for NGC 1316 and NGC 1097, using β_c = 1.5 does provide a better prediction of the observed 870 μm (observed-to-modelled ratio closer to unity) than the ‘β_c = 2.0’ case.

The four remaining galaxies (NGC 0337, NGC 0628, NGC 1512 and NGC 7793) show an 870 μm excess when β_c = 2.0 is used. For NGC 0628 and NGC 1512, using a flatter emissivity index β_c = 1.5 can reconcile the observed and modelled 870 μm global emission within the uncertainties, as shown by their observed-to-modelled ratio close to unity. For the two remaining galaxies NGC 0337 and NGC 7793, the submm slope constrained from Spitzer+Herschel data is already flat (β_c < 1.7 when free to vary) and we observe an excess of the 870 μm emission compared to the extrapolations, even when an emissivity β_c = 1.5 is used (∼50 per cent above the model predictions on average), even if this excess is statistically weak (inferior to a 2σ level). A cold dust component with an even lower emissivity index could be necessary in these galaxies to match the observations. Including the 870 μm data in the fitting procedure leads to β_c values of 1.33 ± 0.15 and 1.19 ± 0.14 for NGC 0337 and NGC 7793, respectively, thus at the lower end of values typically used in the literature. This questions the origin of the excess detected at 870 μm and the use of a simple isothermal component to fit the cold dust population up to 870 μm in those two objects.

As mentioned before, using β_c = 1.5 reconciles the observed and modelled 870 μm global emission for NGC 1512, within the 870 μm uncertainty. Nevertheless, when β_c is allowed to vary, the model fit to Spitzer and Herschel data favours very low β_c values (β_c = 1.16). Including the 870 μm data in the fitting procedure, we find β_c = 1.22 ± 0.15 for NGC 1512. The 870 μm observation thus confirms the shallow slope. The 870 μm observation reveals an extended structure in the south-east of NGC 1512 that is not observed in the SPIRE maps (Fig. 1) and accounts for 30–40 per cent of the global flux. Albrecht et al. (in preparation) could not firmly determine the nature of this extended emission that contaminates the map and therefore the global 870 μm flux estimates. They note that the global flux density is only a factor of 1.3 higher than the Planck value reported for this object in the Early Release Compact Source Catalogue (Planck Collaboration et al. 2011a, GAUFLUX technique). This is coherent with the contamination we estimate.

3.4 Dependence with global properties

We now analyse how the excess (or deficit) of the observed 870 μm emission compared to our extrapolations varies with the global characteristics of our objects. In Fig. 3, we plot the ratios between the observed 870 μm flux densities and the 870 μm estimates extrapolated from the β = 2 (top panel) and β = 1.5 (bottom panel) models as a function of the average metallicity expressed as oxygen abundances (left) and as a function of the total stellar mass in M_⊙ (middle). Metallicities, listed in Table 1, are taken from Kennicutt et al. (2011) and obtained using the calibration of Pilyugin & Thuan (2005). We use the stellar masses derived by Skibba et al. (2011). The observed 870 μm flux densities are corrected for free–free and synchrotron emission and CO(3–2) line contribution. In this figure, the ratio of the observed-to-modelled flux can be used as a ‘proxy’ for submm slope variations to a certain extent, since it characterizes how the observed slope varies compared to a theoretical slope of β = 2 or β = 1.5.

Figure 3.

Ratio between the observed 870 μm flux densities and the 870 μm estimates predicted from the β = 2 (top) and β = 1.5 (bottom) model as a function of, from left to right, metallicity (expressed as oxygen abundances), stellar mass in (M_⊙) and total infrared luminosity (in L_⊙). The observed 870 μm flux densities are corrected for free–free and synchrotron emission and line contribution (see Section 3.1). Labels indicate the galaxies ‘NGC’ numbers. The horizontal line indicates unity.

Open in new tabDownload slide

The difference we derive between the observed and the modelled 870 μm flux densities varies from one object to the other but does not seem to show a particular trend with metallicity, at least within the average oxygen abundance range probed by our sample. We however find significantly higher observed fluxes than those extrapolated from our models in NGC 0337 and NGC 7793, two of the three galaxies with the lowest metallicities. This could potentially be a sign of a flattening of the submm slope (or decrease of the effective emissivity index β) or a presence of ‘submm excess’ with metallicity. We observe a weak trend with the global stellar mass of the galaxies as well as with galaxy morphology (not plotted in this figure). The galaxy with the highest metallicity (and the highest stellar mass) NGC 1316 seems to be an outlier to this general trend but the uncertainty on the ratio is significant for this object.

Fig. 3 (right) shows the ratios plotted as a function of the total infrared luminosities of the galaxies, a quantity often used to derive star formation rates (Kennicutt 1998). The total infrared luminosities of the KINGFISH sample are taken from Galametz et al. (2013b) and reported in Table 1. No clear dependence is observed with this parameter.

4 THE 870 μm EMISSION ON LOCAL SCALES

In order to understand the range of emissivity index values we observe on global scales, we now try to characterize variations of the submm slope at 870 μm within our galaxies.

4.1 Extrapolations using 2MBB techniques

We convolve our Spitzer and Herschel images to the lowest resolution of the data set, namely that of SPIRE 500 μm, using the convolution kernels generated by Aniano et al. (2011). Maps are then projected to a common sample grid with a pixel size of 14 arcsec (the standard pixel size of the SPIRE 500 μm images). This size corresponds to ISM elements ranging from ∼265 pc for our closest galaxy NGC 7793 to ∼1.4 kpc for NGC 1316. We then use our two-temperature fitting technique to model the SED in each pixel. We only perform the SED fitting for ISM elements above the 2σ detection at 870 μm (according to the signal -to-noise maps produced by Albrecht et al. in preparation) but the LABOCA data are not included in the fitting process. We refer to [G12] for a detailed discussion of the various cold dust temperature and emissivity index maps derived using MBB fitting. The update of the SPIRE beam sizes mentioned in Section 3 does not affect our conclusions on the distribution of these cold dust properties.

As previously mentioned in Section 3.2, we choose to fix β_c in the following analysis. This will (i) prevent biases (T-β degeneracy, temperature mixing) from affecting our maps and (ii) enable us to study the behaviour of the submm slope on resolved scales compared to models for which dust properties are perfectly constrained. Beyond the 870 μm excess issue, we thus also probe potential variations of the submm slope in our objects.

We produce 870 μm maps of our 11 galaxies for two different cases, β_c fixed to 2.0 and 1.5. We also convolve the LABOCA images to the SPIRE 500 μm resolution using the convolution kernels developed by Aniano et al. (2011) (see Section 2.2 for web reference) and regrid the maps to our common pixel grid (pixel size: 14 arcsec). For this study, the 870 μm maps have not been corrected for CO or radio contamination. We discuss their potential contribution to the observed 870 μm fluxes in Section 4.4. Fig. 4 gives an example of the observed 870 μm map of NGC 3627 (top line) along with the 870 μm maps extrapolated from the models with β_c = 2.0 (second line, first panel) or 1.5 (second line, middle panel). The maps of the full sample derived using the 2MBB techniques are given in appendix in Fig. A1. We use the same colour scale for the observed and modelled 870 μm maps to allow a direct comparison.

Figure 4.

Top panel: flux density at 870 μm map observed with LABOCA (in MJy sr⁻¹) for the galaxy NGC 3627 not corrected for CO and radio continuum. Second line panels: modelled 870 μm maps (in MJy sr⁻¹) derived, from left to right, using our two-MBB procedure with β_c fixed to 2, fixed to 1.5 or using the [DL07] formalism. Third line panels: absolute difference between the observed and the modelled 870 μm map in MJy sr⁻¹ (defined as observed flux – modelled flux) for the different SED models (see Section 4.3). Bottom line panels: Relative difference defined as (observed flux − modelled flux)/(modelled flux). We provide the colour scales for each line on the right-hand side. Zeroes are indicated with a white marker on the colour bars. The signal-to-noise of the LABOCA 870 μm map is overlaid as contours on the absolute and relative difference maps (same levels than in Fig. 1).

Open in new tabDownload slide

4.2 Extrapolations using the Draine & Li (2007) model

We want to compare the 870 μm estimates obtained from the 2MBB resolved techniques with those predicted by a more physical (as far as dust composition and ISM physics are concerned) SED model. We thus also build pixel-by-pixel 870 μm maps using the resolved SED modelling of the KINGFISH sample provided in Aniano et al. (2012) and Aniano et al. (in preparation) where Spitzer + Herschel data up to 500 μm are modelled using the Draine & Li (2007) dust models (hereafter [DL07]) on a local basis. Maps are derived at the resolution of SPIRE 500 μm, for a 14 arcsec × 14 arcsec pixel grid. We refer to the two previously mentioned studies for details on the SED modelling technique. We will remind the reader of a few principles in this section. In addition to the KINGFISH sample, this SED fitting technique has also been used on M31, modelled out to 1.5 optical radius in Draine et al. (2013). They show that deviations from the [DL07] model appear in the centre (galactocentric distances lower than 6 kpc), suggesting steeper submm opacity spectral index but that the DL07 model can fit the emission out to 500 μm elsewhere, within the uncertainties.

In the [DL07] dust models, the dust size distribution and its composition are considered to be uniform in each resolution element for which the model is run. Sources of the infrared emission are old stars, Polycyclic Aromatic Hydrocarbons, graphite and silicate grains. The stellar emission (>3 μm) is mostly constrained by the two first Spitzer/IRAC bands and modelled by scaling a blackbody function, using a photospheric temperature of 5000 K (Bendo et al. 2006). In each pixel, the dust is assumed to be heated by a distribution of starlight intensities, but with most of the dust heated by a single starlight intensity that is interpreted as the starlight intensity in the diffuse ISM. The spectral shape of the radiation field is assumed to be that of the Galactic diffuse ISM estimated by Mathis, Mezger & Panagia (1983). The modelling procedure uses the bandpasses of the instruments to apply colour-corrections and convergence to a preferred model is reached using a χ² minimization technique.

We extrapolate 870 μm maps of the sample from the [DL07] resolved SED models of Aniano et al. (in preparation). Fig. 4 shows an example of the modelled 870 μm map for NGC 3627 (second line, third panel). We provide the 870 μm maps of the full sample derived with this SED modelling technique in appendix in Fig. A1 (bottom panels). Here, again, the same colour scale was applied between observed and modelled 870 μm maps to allow a direct comparison.

The modelled 870 μm maps obtained using the [DL07] formalism are quantitatively similar to those obtained using a 2MBB model with β_c = 1.5. In the wavelength range covered by SPIRE, the SED extrapolated from the [DL07] models is in fact similar to a single temperature MBB with β_c = 2.0 but the model can reach colder temperatures due to the mixture of temperatures it includes. The [DL07] dust models also incorporate small modifications to the amorphous silicate opacity, especially for λ > 250 μm, in order to better match the average high Galactic latitude dust emission spectrum measured by Cosmic Background Explorer-Far Infrared Absolute Spectrophotometer (COBE-FIRAS; Wright et al. 1991; Reach et al. 1995; Finkbeiner, Davis & Schlegel 1999). Thus, [DL07] models already include a ‘Galactic submm excess’ (modifications inferior to 12 per cent for the 250 μm< λ < 1100 μm range). Emission in excess of the [DL07] models thus means an excess compared to what could be explained by Galactic-like dust only.

4.3 Absolute and relative difference maps at 870 μm

When β_c is fixed to a standard value of 2, we systematically observe regions with a positive difference between the observed and the modelled 870 μm emission in our maps. This positive difference can be localized in the centre of the galaxy, like in NGC 1097, or extended across the object (NGC 0337 or NGC 7793 for instance). The 870 μm maps extrapolated from the β_c = 1.5 model are brighter than those with β_c = 2.0 due to the flatter submm slope induced by this SED model that leads to higher extrapolations at 870 μm. As a direct consequence, the absolute difference derived using β_c = 1.5 is systematically smaller than that derived using β_c = 2.0.

We now qualitatively and quantitatively describe the trends observed for each individual galaxy. Part of the excess we detect can be linked to contamination by free–free, synchrotron, or CO emission. We quantify this contribution in Section 4.4.

NGC 0337. The distributions of the observed and modelled 870 μm emission are very similar. However, quantitatively, the observed 870 μm is significantly higher than our model predictions, regardless of the model we use. The absolute excess above the model prediction at 870 μm is radially decreasing, following the distribution of the star formation activity of the galaxy as well as that of the dust mass surface density mapped in [G12]. The relative excess seems to be non-homogeneous across the galaxy. The relative excess in regions detected with LABOCA at a 5σ level is ∼70 per cent on average when β_c = 1.5, with a standard deviation of 25 per cent. This is slightly higher than what is found on global scales for this object (but we remind the reader that we are excluding the lowest surface brightness regions here). Unfortunately, because this is one of most distant galaxies of the sample, our poor resolution does not allow us a very detailed study of the 870 μm residual distribution.

NGC 0628. As one can see from Fig. 1, part of the emission in the south-east side of the galaxy is not detected in the LABOCA map. The modelled 870 μm maps are thus brighter than the emission actually observed with LABOCA, whatever the model, and the absolute difference thus shows a deficit of emission in that south-east region. For the rest of the galaxy, the distributions of the observed and modelled maps are similar, with the 870 μm emission (modelled and observed) following the distribution of the emission traced by SPIRE wavebands. We observe a strong 870 μm excess above the β_c = 2.0 model across the galaxy. The difference between the observed and modelled 870 μm emission partly but not fully disappears when we use a lower β_c = 1.5 or the [DL07] dust model. This excess emission is not homogeneously distributed across the disc of the galaxy and does not follow the galaxy structure. Indeed, the absolute and relative difference maps of NGC 0628 do not show a statistically significant excess emission in the centre but peaks in the low-surface brightness regions of the galaxy where cold dust at 870 μm emission is detected by LABOCA at a reasonable 3σ level.

NGC 1097. The observed and modelled 870 μm maps are similar, with a very strong 870 μm emission in the centre of the galaxy (the galaxy possesses a moderate AGN) and along the bar, and two fainter spiral arms. A positive difference between the observed and modelled map is seen across the galaxy when using the β_c = 2.0 MBB model. The 870 μm emission seems to be better explained when a β_c = 1.5 MBB model or the [DL07] formalism are used. An excess above the model prediction remains nonetheless in the central region when β_c = 1.5 (we obtain an average of ∼54^±35 per cent in regions with a 20σ detection) and in the few inter-arm resolved elements that are modelled.

NGC 1291. The three modelled 870 μm maps follow the distribution of the SPIRE maps, namely a bright centre elongated in the north–south direction and a diffuse cold dust ring (prominent longward of 160 μm) with localized bright spots in the ring in the north-west and south-east arcs (see Hinz et al. 2012, for a detailed analysis of this cold dust ring with Herschel). We analyse the excess in resolved elements detected at 870 μm to a 2σ level, which, unfortunately, restricts the study to the centre of the galaxy where a strong excess above the β_c = 1.5 model prediction is detected (∼59 per cent on average in the regions we modelled, with a standard deviation of 30 per cent). The major part of the cold dust ring is not detected with LABOCA (see Fig. 1), which prevents us to conclude about any trend in the ring of NGC 1291. It also implies that the global 870 μm flux density provided for this galaxy in Albrecht et al. (in preparation) should be considered as a lower limit. A higher global flux density would reconcile the 870 μm data with the integrated SED models derived in Section 3 for this object.

NGC 1316. The modelled and observed maps of the galaxy follow the distribution of the warm or cold dust emission traced by Herschel longward of 70 μm, with two peaks located on each side (north-west and south-east) of the star-forming centre (indicated by the MIPS 24 μm central peak). We observe a positive difference between the observed and the modelled 870 μm emission throughout the galaxy when the β_c = 2.0 MBB model is used. The ‘excess’ follows the distribution of star-forming regions (namely a central peak co-spatial with the 24 μm peak) rather than that of the cold dust emission or the dust mass surface density as mapped in [G12]. When using the [DL07] model prediction, the very central relative excess is about ∼70 per cent above the modelled 870 μm emission on average. We note that this galaxy possesses a moderate active galactic nucleus (AGN). We discuss potential non-dust contamination (synchrotron emission in the centre in that case) in Section 4.4.1.

NGC 1512. The LABOCA emission follows the same distribution but is systematically brighter than our model predictions (∼40 per cent brighter in the northern part of the galaxy when using the β = 1.5 or the [DL07] models). We detect a significant excess emission above our model predictions in the south-east part of NGC 1512, due to a faint extended structure detected at 870 μm (Fig. 1). Faint blobs are detected in the SPIRE maps in that part of the galaxy but it is difficult to determine if the SPIRE emission belongs to the galaxy due to confusion in the SPIRE maps. NGC 1512 is in interaction with NGC 1510 located in the south-west of the map in Fig 1. We refer to the work of Koribalski & López-Sánchez (2009) where the H i emission of the pair is analysed. One might argue that the diffuse 870 μm emission we observe to the south could be a continuation of the NW-SE diffuse 500 μm tidal arm observed between the two galaxies and result from very cold dust in compressed gas. Unfortunately, Albrecht et al. (in preparation) could not firmly conclude on whether or not the emission is real or an artefact in the 870 μm map itself.

NGC 3351. The modelled 870 μm maps show a central peak in the emission surrounded by a diffuse disc emission but only part of the emission in the disc is detected with the LABOCA instrument (see Fig. 1). For this galaxy as well, the non-detection of disc structure with LABOCA indicates that the global 870 μm flux density provided for this galaxy in Albrecht et al. (in preparation) should be considered as a lower limit. An emission in excess is systematically detected in the centre of the galaxy, whatever the model, but the relative 870 μm difference is small (<30 per cent for regions with a 3σ LABOCA detection). An 870 μm excess above the β_c = 1.5 model is also observed and radially increasing to a factor of 2 in the north-east region of the galaxy, i.e. in the only region of the disc detected with LABOCA at a 3σ level.

NGC 3621. Our modelled maps qualitatively match the observed distribution. They show a radial decrease of the 870 μm emission and two peaks in the centre of the galaxy, the southern being the strongest. However, the observed 870 μm emission is brighter than our β_c = 2.0 model prediction, with a distribution of the absolute difference map following the galaxy structure and the dust mass surface density mapped in [G12]. This trend remains when the β_c = 1.5 model or the [DL07] model are used. The excess is, however, statistically weak in these two cases, with an 870 μm relative excess of 19^±12 per cent on average in regions with a 5σ LABOCA detection. Peaks of the relative excess emission (barely visible for the β_c = 1.5 or [DL07] models) are distributed in low-surface brightness regions rather than in the centre of the galaxy.

NGC 3627. The 870 μm map follows the distribution of the SPIRE observations, namely a strong emission in the centre and peaks at the bar ends (bar in the NW-SE direction) where molecular hydrogen (Wilson et al. 2012) and large dust reservoirs ([G12]) are located. The star formation primarily occurs in these two knots where the SF efficiency is enhanced (Watanabe et al. 2011). For all models, the absolute 870 μm excess follows the spiral structure of the galaxy (and the dust mass surface density), with excess maxima corresponding to 870 μm peaks and spiral arms still distinguishable. The relative 870 μm excess above the β_c = 1.5 or [DL07] models shows moderate peaks in the centre and each end of the bar. We obtain an average relative excess of 30 per cent for these regions (detected by LABOCA at a 20σ level), with a standard deviation of 19 per cent.

NGC 4826. The map of the 870 μm absolute difference between observation and predictions, like that of the 870 μm emission itself, shows a smooth radial decrease of the excess in this object. The relative difference is non-homogeneous and statistically weak on average across the whole object in the β_c = 1.5 or [DL07] cases. This is consistent with what is found at global scale for NGC 4826.

NGC 7793. The modelled maps derived for the galaxy follow the distribution of the cold dust emission traced by SPIRE and the 870 μm observation. NGC 7793 is a flocculent galaxy with peaks of star formation spread out in the disc (but only partly resolved at the resolution we are working at in this study) and most of the disc is detected with LABOCA at a level >3σ. Like for the galaxy NGC 0337 (a galaxy that shows a similar flattening of the submm slope at global scale), we observe a strong 870 μm absolute excess throughout the galaxy. It however does not follow the spiral structure in this object and peaks are observed within the disc rather than in the centre (the relative excess is lower than 30 per cent in the nucleus). The trend is confirmed by the relative difference maps in which the 870 μm difference between the observations and the models seems to radially increase (by more than a factor of 2) towards the south of the galaxy.

For half of the sample, the global fit favours steeper slopes (β_c = 2.0) compared to the resolved study. This could primarily be linked with the fact that we are restricting our resolved study to pixels with good detections in the LABOCA data, excluding the faint outskirts where the 870 μm emission is weak. These global versus resolved discrepancies could also partly be due to effects of non-linearity in the SED models we are using and reported for instance in Galliano et al. (2011) or [G12]. Comparing estimates obtained on global and local scales, they show that the dust mass obtained from integrated fluxes could be much lower than that obtained on a resolved basis and attribute this effect to a possible dilution of the cold dust regions in hotter regions. Integrated models could consequently be unable to account for the cold dust populations, biasing the SEDs towards warmer dust (or steeper submm slope) in MBB models or higher minimum heating intensities for the [DL07] models. The global versus resolved discrepancies we observe re-inforce the importance of systematically confronting integrated results (obtained for instance for high-redshift sources that are lacking spatial resolution) and resolved properties derived in the nearby Universe.

4.4 Non-dust contribution to the 870 μm emission

4.4.1 Free–free and synchrotron contribution

As mentioned in Section 3, free–free and synchrotron emission can be a source of contamination of our 870 μm observations but does not significantly (<2 per cent) contribute, at global scale at least, to the 870 μm continuum emission in most of our objects. In the galaxy centres, contribution to the 870 μm by synchrotron emission could increase in case the galaxy hosts an AGN. Two galaxies, NGC 1097 and NGC 1316 possess moderate AGNs. Quantifying the contribution of synchrotron to the 1.3 mm emission in a large sample of galaxies (they neglect the free–free emission), Albrecht, Krügel & Chini (2007) showed that the mean contribution of AGN is on average of ∼8 per cent at this wavelength while less than 2 per cent in normal galaxies. We thus expect a more significant non-dust contribution from the synchrotron emission in the centre of NGC 1097 and NGC 1316 than for the rest of the sample. Radio contamination could also moderately affect the central resolved elements of other galaxies whose nucleus is classified as ‘AGN’ in the Moustakas et al. (2010) classification. Finally, the properties of the dust populations in the centre of NGC 1097 could also be affected by the presence of a prominent starbursting ring (Hummel, van der Hulst & Keel 1987).

4.4.2 Local CO contribution

The ¹²CO(3–2) line emission can also contribute to the 870 μm emission in some galaxies of the sample. The 870 μm maps we use to derive the difference maps are not corrected for the CO(3–2) contribution. Albrecht et al. (in preparation) provide details of the CO observations available for our sample from which they derive the global estimates of the CO line contribution at 870 μm we quote in this section. At present, there are no CO observations available for the galaxies NGC 1512 and NGC 3621.

From single-point measurements. Direct CO(3–2) single-point measurements (∼22 arcsec beam) were taken towards the centre of NGC 1097 (Petitpas & Wilson 2003) and NGC 4826 (Mao et al. 2010). For NGC 1097, the contribution of the 3–2 line to the measured 870 μm emission in this beam is |$S_{\rm CO(3{\rm -}2)}$| ∼ 101 mJy. This represents 20 per cent of the total LABOCA flux contained in the 56 arcsec central region. Added to the synchrotron contribution to the 870 μm (8 per cent), this could mostly explain the central excess (40 per cent on average) we observe. For NGC 4826, a contributing |$S_{\rm CO(3{\rm -}2)}$| flux density of 122 mJy is estimated, representing 14 per cent of the total flux, here again a non-negligible contribution to the central 870 μm excess emission. A CO(2–1) single-point measurement (22 arcsec beam) was taken towards the centre of NGC 1316 (Horellou et al. 2001). Using a brightness temperature ratio R_{3-2, 2-1} = 0.36 (mean value of the ratio derived in Wilson et al. 2012), Albrecht et al. (in preparation) estimated a contribution of S_CO(3-2) ∼ 12 mJy, which represents, like in NGC 1097, up to 20 per cent of the central 56 arcsec central region. Here, again, added to the radio contribution to the 870 μm, we can explain half of the emission in excess in the centre. Finally, only CO(1–0) single-point measurements (43 arcsec beam) are available towards the centre of NGC 1291 (Tacconi et al. 1991) and NGC 7793. Using a brightness temperature ratio R_{3-2, 1-0} = 0.18 (mean value of the ratio derived in Wilson et al. 2012), Albrecht et al. (in preparation) estimate the |$S_{\rm CO(3{\rm -}2)}$| flux density to reach 28 mJy at most in NGC 1291, minor compared to the excess emission we observe, and 240 mJy in NGC 7793. This represents in both cases ∼6.5 per cent of the 870 μm emission, not sufficient to explain the excess we observe. We note that in NGC 7793, the single-point observation is, anyhow, not sufficient to probe CO contribution in the low-surface brightness regions where the excess emission is mostly detected.

From CO maps. CO(2–1) maps of NGC 0337, NGC 0628, NGC 3351 and NGC 3627 were obtained as part of the HERA CO-Line Extragalactic Survey (HERACLES; Leroy et al. 2009). We derive maps of the CO(3–2) line contribution to the 870 μm maps, using a brightness temperature ratio R_{3-2, 2-1} = 0.36 to convert the CO(2–1) to CO(3–2) emission. Fig. 5 (left-hand column) shows the CO(3–2)-to-870 μm flux density ratios for these four galaxies. MIPS 24 μm contours are overlaid to indicate the distribution of the star-forming regions across the galaxies. We also show the absolute difference maps (Fig. 5, middle column) and the relative difference maps (Fig. 5, right-hand column) now corrected for CO(3–2) emission. The modelled 870 μm maps used in this analysis are those obtained using the [DL07] model. For NGC 0337 and NGC 0628, the CO(3–2) line contamination is estimated to represent less than 3 per cent of the 870 μm emission on average across the galaxy, with a peak at 12 per cent in the centre for NGC 0628. This contribution is not sufficient to explain the 870 μm excesses we observe above our model extrapolations for these two objects. For NGC 3351, an 8 per cent CO(3–2) line contamination to the 870 μm emission is derived in the centre, which can fully explain the weak central excess we detect for this object. Finally, an 11 per cent line contamination is derived on average for NGC 3627, but with a peak of the CO contribution to the 870 μm in the centre of the galaxy. After correction, a weak relative excess (<0.5 so 50 per cent above the modelled flux) remains in the centre and on each side of the bar of the galaxy. Directly using the CO(3–2) observation of NGC 3627 taken from the James Clerk Maxwell Telescope (JCMT) Nearby Galaxy Legacy Survey (NGLS; Wilson et al. 2012) lowers the relative excess to 40 per cent in the centre and at each end of the bar.

Figure 5.

Left column: CO(3–2) contribution to the 870 μm emission (in per cent, see Section 4.4.2). Middle column: excess emission on absolute scales (in MJy sr⁻¹) corrected for CO(3–2). The modelled 870 μm map is obtained using the [DL07] formalism. Right column: corresponding relative excess corrected for CO(3–2). MIPS 24 μm contours are overlaid: for NGC 0337: 0.4, 1.5, 4 and 7 MJy sr⁻¹); for NGC 0628: 0.25, 0.7, 1.5 and 2.1 MJy sr⁻¹; for NGC 3351: 0.1, 0.3, 1.1 and 10 MJy sr⁻¹; for NGC 3627: 0.6, 2.2, 6 and 15 MJy sr⁻¹.

Open in new tabDownload slide

4.5 Analysis

4.5.1 Average radial dependence

Resolved analysis of the IR emission of nearby galaxies have showed that dust properties are particularly expected to vary from the nuclei to the outskirts of galaxies (see for instance the radial profiles of various dust properties like the TIR-to-UV ratio, the dust-to-gas mass ratio or the Polycyclic Aromatic Hydrocarbons fraction in Muñoz-Mateos et al. (2009), of dust abundance gradients in Mattsson & Andersen (2012) or the bulge versus disc comparison of Engelbracht et al. (2010), among others). In order to study the radial dependence of our results, we plot in Fig. 6 the absolute (top left) and relative (bottom left) difference between the observed and modelled 870 μm emission (using the [DL07] formalism) as a function of galactocentric radius. For each object, we normalize the radius to the approximate major radius of the galaxy. These major radii are taken from Kennicutt et al. (2011) and approximately match the r₂₅ radii of de Vaucouleurs et al. (1991). We exclude the galaxies NGC 1291 and NGC 1512 from this analysis because of the low statistics we have in the first case and the unknown nature of the 870 μm southern structure in the second case. For the remaining objects, we restrict the analysis to resolved elements with a signal-to-noise superior to 5 in the 870 μm band. We bin our resolved elements in 10 annuli from the centre to r₂₅ (radii are not corrected for inclination) and take the medians in these annuli. The optical radius of NGC 0337 is only 1.5 arcmin, which corresponds to six resolved elements of 14 arcsec. Elements are thus binned in 5 annuli rather than 10 for this object. We calculate the median when three resolved elements at least fulfil our criteria to work on better statistics. Uncertainties are the standard deviations in each annulus.

Figure 6.

Top: absolute difference between the observed 870 μm flux densities (not corrected for CO and radio continuum) and the 870 μm estimates extrapolated from the [DL07] dust model as a function of galactocentric radius (left, see Section 4.5.1) and 24 μm surface brightness in units of W kpc⁻² (right, see Section 4.5.2). The different galaxies of the sample are colour coded. The horizontal dashed line indicates when the observation matches the model. Bottom: same with relative differences.

Open in new tabDownload slide

The absolute differences decrease with radius in most cases, with a systematic 870 μm excess (positive difference) in the centre (typically regions located within 0.3 r₂₅). Some galaxies like NGC 0337 show a smooth decrease with radius. Others have a sharp 870 μm excess profile in the very centre such as for NGC 1097, NGC 3627 or NGC 4826, essentially linked with the contribution of CO and radio to the 870 μm emission, as discussed in the previous section. The galaxy NGC 1316 hosts a low-luminosity X-ray AGN (Kim & Fabbiano 2003) but the possible radio contamination of the 870 μm emission does not lead to a similar broken profile. We observe a sharp decrease of the absolute difference in NGC 3351 as well but this could be due to the non-detection of the 870 μm emission previously mentioned, leading to an underestimate of the observed 870 μm flux density across most of the disc of the galaxy. The two non-barred objects NGC 0628 and NGC 7793 have very flat profiles, i.e. almost no dependence of the absolute excess compared to model predictions on radial distance. The results of [G12] also suggest differences in the submm spectra between barred and non-barred galaxies (when β_c is allowed to vary). We discuss this further in Section 4.5.3.

The radial profiles of most of the galaxies are confined to the ±50 per cent region, which suggests that the [DL07] model is, on average per radial bin, sufficient to explain the 870 μm emission we observe, within uncertainties. We observe a noticeable relative difference in NGC 0337 for which the 870 μm difference is systematically above the ±50 per cent threshold. Even if the relative 870 μm excess of NGC 0337 and NGC 7793 seems to increase with radius across the disc of the objects before dropping in the outskirts, the trend is, within the uncertainties, consistent with a constant relative (and positive) difference. Some of the decreases obtained in the outskirt regions could be linked with the low statistics (and low 870 μm surface brightnesses) towards r₂₅ radii.

4.5.2 Dependence on the star formation activity

Several studies have tried to correlate the submm excess with the IR emission. Studying sub-regions along the radius of NGC 4631, Bendo et al. (2006) showed, for instance, that their excess 850 μm emission significantly varies in the galaxy and that low excess values were predominantly found in faint IR regions, even if the trend is statistically weak. In Fig. 6, we plot the absolute (top right) and relative (bottom right) difference between the observed and the modelled 870 μm emission (using the [DL07] formalism) as a function of the 24 μm surface brightness. The 24 μm binning differs from the radial binning, the 24 μm emission being indeed sensitive to the warm dust populations and a good tracer for active star formation (see Calzetti et al. 2007, among others). We use the same detection criteria and require the same number of resolved elements to calculate the median in each bin of surface brightness than those applied for the radius profiles. Uncertainties are the standard deviations in each annulus.

For all galaxies, the absolute differences increase with the 24 μm surface brightness, with an 870 μm emission consistent with the [DL07] model prediction in low 24 μm surface brightnesses and an excess at high surface brightnesses. More particularly, we observe a systematic 870 μm excess emission when νL_ν(24 μm) > 5 × 10³³ W kpc⁻², with a significant 870 μm absolute excess emission in the bright star-forming regions (νL_ν(24 μm) > 5 × 10³⁴ W kpc⁻²) of NGC 1097, NGC 3351, NGC 3627 and NGC 4826. This is consistent with the results of Section 4.5.1, the bright star-forming regions being located in the galaxy centres (or close in the case of NGC 3627) and here again the central excess is partly (if not mostly) linked with contribution of CO and radio to the 870 μm emission. The relative difference profiles show a variety of behaviours within the sample, with large uncertainties, especially in regions with low 24 μm surface brightness. The 870 μm relative differences for the galaxies NGC 0337 and NGC 7793 show rather flat profiles, so no dependence with the star formation activity either. This absence of clear dependence makes it difficult to understand what could be driving the excess emission we detect at global scales. The galaxies NGC 0628, NGC 3351 and NGC 3621 show an increase towards star-forming regions, even if the radial averages have large uncertainties (partly due to non-detections in the disc for NGC 0628 and NGC 3351). We discuss these galaxies in more detail in the following section.

4.5.3 Decrease of the effective emissivity index?

In Fig. 6, we observe from the profiles of resolved elements with a 870 μm detection at 5σ averaged in annuli that the relative excess does not show a particular dependence with radius. However, we can see from Fig. A1 that the relative differences between observed and modelled 870 emission preferentially peak in the disc or the outskirts for the galaxies NGC 0628 and NGC 7793 (in NGC 3621 as well to a lesser extent) when probed in low-surface brightness regions. Whatever causes the excess emission in these objects seems to have a predominant effect in these regions rather than in the centre. Modelling the thermal dust emission in a strip of the Large Magellanic Cloud, Galliano et al. (2011) also derived a submm excess map of the region (excess at 500 μm in their study) and found that the excess is inversely correlated with the dust mass surface density, thus increasing towards low-surface brightnesses as well. In [G12], we showed that when the temperature and the emissivity are allowed to vary in the resolved 2MBB fitting process, the ‘effective’ emissivity index seems to radially decrease towards low surface brightness regions. The high 870 μm excess emission observed in the north-east part of the disc in NGC 3351 (unfortunately the only region of the disc detected with LABOCA) also corresponds to a region where low values of the emissivity index β_c were derived in [G12]. The two independent results (870 μm excess detected in the disc or outskirts when β_c is fixed versus radial flattening of the [250–500μm] slope when β_c is a free parameter) could both be explained by a radial flattening of the submm slope in those particular objects.

NGC 0628, NGC 3621 and NGC 7793 are non-barred galaxies, and three of the four late-type objects of the sample (SAc, SAd and SAd morphology, respectively). The fourth late-type galaxy is NGC 0337 in which we also detect an excess predominant in the outskirts, even if the compactness of the source does not allow us to properly analyse the excess distribution. The radial decrease of the effective emissivity index was also observed in the Scd galaxy M33 (Tabatabaei et al. 2014) or the SA(s)b galaxy Andromeda (Smith et al. 2012) but not in the SAB(s)c galaxy M83 (Foyle et al. 2012) or the other barred spiral galaxies of our KINGFISH/LABOCA sample ([G12]). The submm radial flattening thus seems to be dependent on the galaxy morphology and favoured in late-type non-barred spirals.

At our working resolution unfortunately, many physical parameters possess a radial dependence. As shown in Moustakas et al. (2010), NGC 0628, NGC 3621 and NGC 7793 possess, for instance, strong radial metallicity gradients and metallicity is, as discussed in Section 3.4, a driver of variations in the properties of cold dust (Sandstrom et al. 2012; Rémy-Ruyer et al. 2014). We also know that the dust temperature slowly decreases with radius or that the molecular gas distribution varies and could rapidly become patchy towards the outskirts of spirals. The low resolution of our data prevents us from studying where the excess is produced on sub-kpc scales. The lack of coverage (gap between 500 and 870 μm and no constraint longward of 870 μm) also does not allow us to know if the excess can purely be explained by a smooth flattening (so by variations in the dust opacity or by a smooth dust temperature range for instance) or if the 870 μm is produced by a separate population (cold dust, spinning dust, magnetic grains) that could translate into a break in the submm slope.

5 CONCLUSIONS

We combine Spitzer, Herschel and LABOCA data to model with MBBs the global IR-to-submm emission of 11 nearby galaxies. We investigate how the potential excess emission at submm wavelengths varies with various model assumptions (MBB with β_c fixed or free) and show the following.

• For seven galaxies of the sample (NGC 1097, NGC 1291, NGC 1316, NGC 3351, NGC 3621, NGC 3627 and NGC 4826), the global 870 μm emission can be explained, within the uncertainties, by thermal emission using an emissivity index β_c = 2.0.

• For the four other objects (NGC 0337, NGC 0628, NGC 1512 and NGC 7793), the integrated 870 μm emission is in excess compared to that expected from cold dust with an emissivity index β_c = 2.0. Using an emissivity index β_c of 1.5 reproduces the global excess, we observe in NGC 0628, and also lead to better predictions of the 870 μm emission in the previously mentioned two galaxies NGC 1097 and NGC 1316. For NGC 1512, the observed 870 μm flux and the 870 μm modelled using β_c = 1.5 are consistent within the error bars but the model fitted to Spitzer+Herschel data favours very low β_c values when allowed to vary (β_c = 1.16). Lower values (β_c < 1.3) would also be required to fit the 870 μm data in NGC 0337 and NGC 7793, which questions the use of an isothermal component to fit the cold dust population up to 870 μm in those objects.

We apply the same methodology on local scales and produce 870 μm maps using 2MBBs models with fixed β_c or the [DL07] dust models in order to compare them to the observed emission.

• We observe a systematic 870 μm excess when the emissivity index of β_c = 2.0 is used to model the ISM elements we selected (above a given signal-to-noise criterion). Using β_c = 1.5 can, in many cases, reconcile the observed and the modelled 870 μm emission. The (absolute and relative) difference maps derived using the [DL07] formalism are similar to those obtained using β_c = 1.5.

• Maps of the absolute difference between observed and modelled 870 μm emission show, in many cases, a decrease of this quantity with radius and an increase with star formation (traced using the 24 μm surface brightness). In NGC 1097, NGC 1316, NGC 3627 and NGC 4826, this excess can be partly or fully explained by contributions from CO(3–2) emission and/or free–free and synchrotron radiation.

• The relative excess profiles have larger uncertainties with various behaviours from one galaxy to another and no clear dependence in the case of NGC 0337 and NGC 7793 with either parameter (radius or star formation). This raises the issue of what could be the drivers of the major excesses we detect in these two objects. This study should be extended to closer (and thus more resolved) objects to increase our statistics.

• In the non-barred spirals NGC 0628, NGC 3621 and NGC 7793 however, the relative excess maps show peaks in the disc or towards low surface brightnesses. A radial flattening of the submm slope (so a decrease of the ‘effective’ emissivity index with radius) as suggested by the studies of [G12], could explain the distribution of the 870 μm excess we observe.

We would like to thank the referee for a very constructive report that helped to improve the clarity/robustness of the results. The research of CDW is supported by grants from the Natural Sciences and Engineering Research Council of Canada. PACS has been developed by MPE (Germany); UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy) and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy);IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK) and NASA (USA).

REFERENCES

APPENDIX A: 870 |$\boldsymbol \mu {\bf m}$| MAPS PREDICTED FROM OUR DIFFERENT SED MODELS AND 870 |$\boldsymbol \mu {\bf m}$| ABSOLUTE AND RELATIVE DIFFERENCE MAPS

Figure A1.

For each galaxy: top panel: flux density at 870 μm map observed with LABOCA (in MJy sr⁻¹) not corrected for CO and radio continuum. Second line panels: modelled 870 μm maps (in MJy sr⁻¹) derived, from left to right, using our two-MBB procedure with β_c fixed to 2, fixed to 1.5 or using the [DL07] formalism. Third line panels: absolute difference between the observed and the modelled 870 μm map in MJy sr⁻¹ (defined as observed flux – modelled flux) for the different SED models. Bottom line panels: relative excess defined as (observed flux – modelled flux)/(modelled flux). We provide the colour scales for each line on the right-hand side. Zeroes are indicated with a white marker on the colour bars. The signal-to-noise of the LABOCA 870 μm map is overlaid as contours on the absolute and relative difference maps (same levels than in Fig. 1).

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide

Figure A1

– continued

Open in new tabDownload slide