Far-reaching dust distribution in galaxy discs

Abstract

In most studies of dust in galaxies, dust is only detected from its emission to approximately the optical radius of the galaxy. By combining the signal of 110 spiral galaxies observed as part of the Herschel Reference Survey, we are able to improve our sensitivity by an order of magnitude over that for a single object. Here we report the direct detection of dust from its emission that extends out to at least twice the optical radius. We find that the distribution of dust is consistent with an exponential at all radii with a gradient of ∼−1.7 dex |$R_{25}^{-1}$|⁠. Our dust temperature declines linearly from ∼25 K in the centre to 15 K at R₂₅ from where it remains constant out to ∼2.0 R₂₅. The surface density of dust declines with radius at a similar rate to the surface density of stars but more slowly than the surface density of the star-formation rate. Studies based on dust extinction and reddening of high-redshift quasars have concluded that there are substantial amounts of dust in intergalactic space. By combining our results with the number counts and angular correlation function from the SDSS, we show that with Milky Way-type dust we can explain the reddening of the quasars by the dust within galactic discs alone. Given the uncertainties in the properties of any intergalactic dust, we cannot rule out its existence, but our results show that statistical investigations of the dust in galactic haloes that use the reddening of high-redshift objects must take account of the dust in galactic discs.

galaxies: ISM, galaxies: spiral, submillimetre: ISM

INTRODUCTION

Recent studies of the extinction of quasars (Ménard et al. 2010a) and the reddening of galaxies (Peek, Ménard & Corrales 2015) imply that approximately 50 per cent of the interstellar dust in the Universe lies outside galactic discs. However, we still know remarkably little about how far out the dust in galactic discs extends.

Spiral galaxies have huge discs of atomic hydrogen that extend much further than the standard optical radius (Bigiel & Blitz 2012, R₂₅, the radius where the optical brightness corresponds to a B-band brightness of 25 mag arcsec⁻²), but the atomic hydrogen is unprocessed material that may have fallen on the galaxies from outside. Studies into the extent of molecular gas using the carbon monoxide (CO) line, even when combining the results from a very large number of galaxies, have only managed to trace the molecular gas out to 1.1 R₂₅ (Schruba et al. 2011).

There have been many attempts to detect dust around nearby galaxies by using both its absorption and reflection properties in the optical, and its emission properties in the far-infrared. Observations of edge-on galaxies have been used to look for dust above the plane of the galaxy. Hodges-Kluck & Bregman (2014) found that UV emission extends between 5 and 20 kpc from the mid-plane of the galaxies in their sample, and was consistent with the emission being light from the disc scattered by dust above the disc. They found that the halo dust has an exponential form with scaleheights ∼2.5–5 kpc, much larger than typical scaleheights of stellar thick discs (0.3–1.5 kpc; Hodges-Kluck & Bregman 2014). Bocchio et al. (2016) combined Herschel and Spitzer data for the edge-on galaxy NGC 891, finding extended emission above the disc with a scaleheight of 1.44 kpc.

The extent of dust discs has been investigated by looking for the effects of dust extinction of background sources (e.g. quasars, background galaxies) through the outskirts of a nearby galaxy. Holwerda et al. (2009) have shown in a pair of occulting galaxies that dust extinction can be reliably detected out to 1.5 R₂₅. Attempts to detect dust in the outer disc from its emission were undertaken with the early infrared space telescopes, IRAS and ISO. However, they were limited by sensitivity, wavelength coverage and resolution. Nelson, Zaritsky & Cutri (1998) obtained a 2σ detection at 1.5 R₂₅ using IRAS by a stacking analysis, and Alton et al. (1998) used ISO observations of a sample of eight nearby galaxies to show that the dust has a larger exponential scalelength than the optical emission. Recently Spitzer and Herschel data have been used to create radial profiles of the dust for the nearby galaxies in the SINGS/KINGFISH sample (Muñoz-Mateos et al. 2009; Hunt et al. 2015). The KINGFISH radial profiles presented in Hunt et al. (2015) typically extend between 1.0 and 1.5 R₂₅ at 250 μm, across a variety of morphological types.

In this paper, we aim to go beyond the typical detection radius of dust for an individual galaxy, by combining the signal from a large sample of galaxies. Since the dust in the outskirts of galaxies is likely to be very cool (≤15 K) with low surface brightness, the Herschel Space Observatory (Pilbratt et al. 2010) is the only telescope with enough sensitivity to have the potential to directly detect this emission. We have selected galaxies from the Herschel Reference Survey (HRS; Boselli et al. 2010), the largest targeted survey of nearby galaxies with Herschel. These galaxies all have either spiral or irregular morphology and lie face-on, giving us the best possible physical resolution. In Section 2 we present the data used for this analysis. Section 3 describes our method for creating the averaged radial profiles of the dust and other galaxy components. Section 4 presents our results, and the conclusions are presented in Section 5.

THE SAMPLE

The HRS is a volume-limited sample of 322 galaxies with distances between 15 and 25 Mpc (Boselli et al. 2010). To ensure a reasonable detection rate with Herschel, the HRS also has a 2MASS K-band magnitude limit of K_Stot ≤ 12 for late-type galaxies. Observations with the SPIRE instrument (Griffin et al. 2010) at 250, 350 and 500 μm and with the PACS instrument (Poglitsch et al. 2010) at 100 and 160 μm exist for all the galaxies. 83 objects fall within the Virgo cluster and as such have been observed by the Herschel Virgo Cluster Survey (HeViCS; Davies et al. 2010, 2012). For all late-type galaxies the HRS observations were designed to cover an area of at least 1.5 times the optical diameter with a minimum size of 4 arcmin × 4 arcmin. In practice, the SPIRE images generally cover a much larger area around each galaxy, allowing us to trace the dust to larger radii. Each observation of a late-type galaxy has an instrumental 1σ sensitivity of 5.64, 5.65, and 6.60 mJy beam⁻¹ (0.57, 0.32, and 0.18 MJy sr⁻¹) for the 250, 350 and 500 μm bands, respectively. The instrumental noise is comparable to the confusion noise of 5.8, 6.3, and 6.8 mJy beam⁻¹ (Nguyen et al. 2010) for the 250, 350 and 500 μm bands, respectively. Full details about the survey can be found in Boselli et al. (2010).

METHODOLOGY

Radial profiles and averaging

To obtain the highest sensitivity to extended dust we average the emission from a large sample of HRS galaxies to create an average surface-brightness profile. We included galaxies in the sample if they met the following criteria:

• The galaxy morphology must be classified as a spiral or irregular. The classifications are the same as used in other HRS studies (Cortese et al. 2012a). Ellipticals and S0s were not included as they are known to contain much less dust and have much smaller dust discs than spirals and irregulars (Smith et al. 2012a).

• The galaxy cannot have an inclination greater than 60°. The inclinations were calculated from the major and minor diameters measured by the HRS team (Cortese et al. 2012b) using images from the Sloan Digital Sky Survey (SDSS).

• The galaxy must have been detected by Herschel with emission that is clearly extended (Ciesla et al. 2012).

• There must not be another extended Herschel source that could affect the measured radial profile (i.e. the sources must be separated by a distance of at least 6.0 R₂₅).

• The Herschel image must not have significant contamination from diffuse emission from Galactic interstellar dust.

• The image of each galaxy must extend to a distance of at least 2.5 R₂₅ from the centre of the galaxy to make it possible to remove any residual background emission that is not from the galaxy itself.

Once all the criteria had been applied, there were 117 galaxies, providing an order of magnitude increase in sensitivity over what is possible for a single object.

We obtained a surface-brightness profile for each galaxy by calculating the distance to each pixel from the centre of the galaxy and then expressing this distance as a fraction of the optical radius (R₂₅). We calculated this distance using the coordinates of the pixel, the optical radius, the inclination (estimated from the eccentricity), the distance to the galaxy (Cortese et al. 2012a) and the optical centre. A line of fixed radius would appear as an ellipse on a SPIRE map. We extended the profile out until a radial annulus reaches the edge of the image. Before calculating the average profile for the whole sample, we removed a constant offset from each galaxy's profile using the mean value of the surface brightness between 2.2 ≤R/R₂₅ < 2.7 (a couple of profiles do not extend as far as 2.7 R₂₅, but contain enough pixels for an accurate background subtraction). This choice of background subtraction means by construction we are only sensitive to dust within 2.2 R₂₅. We then placed all the pixels from all images into radial bins and calculated the mean value in each bin and its error (⁠|$\sigma / \sqrt{N}$|⁠, where N is the number of data points and σ is the standard deviation). To ensure we account for all uncertainties (e.g. background variations), we also ran a simulation by carrying out an identical analysis using the same number of positions and the same radial bins as in our real analysis at random positions on the HeViCS map (but avoiding Virgo cluster objects). We then performed the same stacking analysis, and re-ran the process a 1000 times. By measuring the standard deviations in each bin, we get a robust estimate for the minimum error in each radial bin. We then use this uncertainty if it is larger than the error above.

The results are shown in Fig. 1. The red-dashed lines in Fig. 1 show the 3σ detection level found from the background simulations described in the previous paragraph. The median optical radius for the sample is 1.25 arcmin (R₂₅ ranges from 0.33 to 4.83 arcmin) and the median distance is 17.0 Mpc. Therefore, the median value of R₂₅ corresponds to a physical radius of 6.2 kpc.

Figure 1.

The average emission from dust plotted against the radial distance for the late-type galaxies in the Herschel Reference Survey. Before averaging, the radial distance has been normalized for each galaxy by the optical radius (R₂₅) of the galaxy. The blue data points show the average radial profile for each wavelength. The wavelength and number of galaxies that has been used in each plot are shown in the top-right corner. The green data points represent the profile of a point source processed using the same method and scaled by the median value of R₂₅ for the sample. The dashed red line shows the 3σ sensitivity limit found from simulations using the background regions of the Herschel Virgo Cluster Survey map.

Reliability

To check that our average profile is representative of the sample of galaxies and not just dominated by a small subset of objects, we performed a Monte Carlo test. We randomly split the sample into four subsets and measured a profile for each subset as described above. For each profile we measured the radius where the surface brightness falls below 0.002 mJy arcsec⁻² at 250 μm. We repeated this process 10 000 times and found no subset where the profile is significantly steeper than the one we show in Fig. 1, showing our result is statistically robust.

An important question is whether at large radii, where the surface brightness is very low, we are effectively measuring the point spread function (PSF) of the observations. To test this, we used SPIRE calibration observations of Neptune. We found the average surface-brightness profile of the PSF using the same procedure as for the galaxies, scaling it to the same physical units as the galaxies using the median R₂₅ of the sample. The Neptune profile has been normalized to the second data point of the averaged galaxy profile and is shown in green in Fig. 1. In the figure we see that in all bands the average galaxy profile is significantly higher at all radii than the expected flux from a point source. The shallow slope of the galaxy profile compared to the PSF precludes the possibility that dust detected at large R is from the Airy rings of the brighter inner parts of the galaxy. Using the median R₂₅ to scale the Neptune profile probably overestimates the effect of the PSF as the galaxies with larger angular sizes contribute more pixels to the average profile.

As a final check that we are not simply detecting at large radii the wings of the PSF, we repeated the analysis for the 45 galaxies with R₂₅ > 1.5 arcmin (the largest is 4.83 arcmin). We find a much larger difference between the average surface-brightness profile and the wings of the PSF (Fig. 2). Fig. 3 shows the morphological-type and stellar mass distributions for the entire HRS, our initial sample of 117 galaxies, and our sample of galaxies with R₂₅ > 1.5 arcmin. Our initial sample of 117 galaxies has a similar distribution to the entire HRS for both morphological type (excluding early-types) and stellar mass, but our sample with R₂₅ > 1.5 arcmin is biased towards higher stellar masses. As mentioned in Section 2, 83 of our galaxies fall within the Virgo cluster, and so the galaxies may be deficient in dust because they are deficient in H i (Cortese et al. 2010). However, we saw little difference in the results when we split the sample into classes based on measurements of the H i deficiency. To minimize the effect of the PSF, we used this sample for the rest of the analysis described in this paper. The median optical radius for this sample is 2.1 arcmin and the median distance is 17.0 Mpc. Therefore, the median value of R₂₅ corresponds to a physical radius of 10.3 kpc.

Figure 2.

The best-fitting models to the measured average surface-brightness profiles of the sample of large galaxies (45 objects). The curved red dashed line shows the 3σ limit (Section 3.1) with the blue points showing the measurements that are above the limit. As in Fig. 1, the green data points represent the profile of a point source, processed using the same method and scaled by the median value of R₂₅ for the sample of large galaxies,. The three models are described in Section 3.3 and are only fit to the blue data points. The colour key for the three models is shown in the legend. Only one model is shown at 350 and 500 μm because the other more complex models did not provide a significantly better fit to the data (see description in Section 3.3). The dashed and the solid lines show the models before and after they have been convolved with the PSF.

Figure 3.

The distributions of morphology and stellar mass for the entire HRS (blue), our initial sample of 117 late-type galaxies (green) and the sample of 45 galaxies with R₂₅ > 1.5 arcmin (red). The morphological classification is the Goldmine (Gavazzi et al. 2003) classification system, where types 3–10 represent Sa to Sdm galaxies, 11–14 are irregular/peculiar galaxies and 14–17 are BCD's. Morphological types −3 to 2 represent early-type galaxies which are not used in this study; the grey-dashed line shows the boundary between early-type and late-type galaxies. The stellar masses are taken from Cortese et al. (2012b). Note that 10 galaxies of the HRS do not have stellar mass determinations, with 7 and 1 of these in the initial and large sample, respectively, and so are not included in the plot.

Model fitting procedure

To fit analytic models to the measured average surface-brightness profiles in Fig. 2, we need to take into account the PSF, which is large at submillimetre wavelengths. We fitted three different models to the 250–500 μm emission:

• A single exponential. This has two free parameters: a scalelength and the normalization of the function.

• Double exponential (sum of two exponentials). This has four free parameters: the two normalizations and scalelength of each function.

• Broken exponential (an exponential where at a given radius the scalelength changes). This has four free parameters: the overall normalization, the scalelength of each exponential, and the radius dividing the regions fitted by the two functions.

The fitting procedure was similar to a standard profile fitting procedure with a chi-squared minimization except that we performed a full 2D-convolution by projecting our model on to an image, convolving the model image with the PSF, creating a model radial profile and then comparing it to the observed radial profile. To decide what is the best model, we compared the reduced chi-squared value of the three models. The fits were performed over the region in which there is at least a 3σ detection at 250 μm, except we did not fit the model to the central two data points (equivalent to central 30 arcsec) to ensure there is no effect from the nuclei. We found that the broken-exponential model provided a significantly better fit at 250 μm but at 350 and 500 μm there was no statistical reason to prefer the more complex models over the single exponential (Fig. 2). Virtually all iterations of the Monte Carlo test described in Section 3.2 show evidence for the shallower slope at large radii suggesting we are not being influenced by a small subset of galaxies.

Before deriving dust temperatures, we first convolved all images to the same resolution as the SPIRE 500 μm images which have the lowest angular resolution of all our images (using the method described in Aniano et al. 2011). We then calculated empirical average radial profiles in the SPIRE bands for the 45 HRS galaxies with R₂₅ > 1.5 arcmin, using the same method as above. Dust temperatures tend to increase towards the centres of galaxies (Smith et al. 2010, 2012b; Foyle et al. 2012; Mentuch Cooper et al. 2012), which means that data at shorter wavelengths than the SPIRE wavelengths are necessary for accurate temperature measurements. We determined average surface-brightness profiles at 100 and 160 μm from the images of the HRS galaxies made with the PACS camera (Poglitsch et al. 2010; Cortese et al. 2014) on Herschel. However, when estimating temperatures we only used PACS points at R < 0.5R₂₅, since background variations in the PACS images appear to dominate at larger radii. While only using three data points to find a temperature is not ideal, at R > 0.5 R₂₅ the dust is colder and so it becomes easier to estimate a dust temperature from the three SPIRE measurements alone than at smaller radii.

We estimated dust temperatures by fitting a modified blackbody to the five or three flux measurements for each radial bin, taking full account of the correlated flux errors between the different Herschel bands (Smith et al. 2012b). As we were limited to only three data points for R > 0.5R₂₅ we had to assume a constant dust-emissivity index (β). We chose to use β = 2 which is the typical value found in resolved studies of nearby galaxies (Eales et al. 2012; Smith et al. 2012b; Hughes et al. 2014; Kirkpatrick et al. 2014). The uncertainties on the dust temperature were calculated by a Monte Carlo approach; in each radial bin every data point was re-sampled using its value and uncertainty (taken correct account of the part of the calibration uncertainty that is correlated between bands) and the modified blackbody fit performed again. From the distribution of temperature values obtained from the 1000 runs of the Monte Carlo simulation, we estimated the uncertainty in the dust temperature (typically ∼2 K). Our temperature profile is shown in Fig. 4. We fitted a simple model to the data points, where the dust temperature decreases linearly with radius until it reaches a minimum. We then derived the dust surface density by combining our best-fitting model for the dust emission and our model of the dust temperature. We assume a value of 0.192 m² kg^{− 1} at 350 μm (Draine 2003) for the dust absorption coefficient (κ_ν).

Figure 4.

The relationship between dust temperature and galactocentric radius from fitting a modified blackbody to the spectral energy distribution of the dust emission (blue points). The dashed line shows a model in which the dust temperature declines linearly with radius until it reaches a minimum. The green line shows the prediction of a model of how the dust in a disc is heated by starlight (Section 4.1).

This value is chosen to match the values used in other HRS and HeViCS papers, and is based on Milky Way-type dust. For this work we have assumed that the value of κ is constant out to a radius of 2.0 R₂₅, as metallicity gradients in galaxies tend to be flat (Moustakas et al. 2010; Sánchez et al. 2014) and that we are measuring dust in the discs of galaxies. However, as very little is known about dust in galaxy outskirts it is possible that the value of κ could have a significant systematic change with radius. The absolute value of κ_ν is also an important consideration when we compare our results with optical reddening studies and so we discuss the assumptions and uncertainties of κ_ν in Section 4.4.

Fig. 5 shows the results. The two solid lines show the surface density of dust we derive when starting from the best-fitting models to the 250 and 500 μm surface-brightness profiles. The dashed lines show fits of single-exponential models to these empirical dust column-density distributions.

Figure 5.

The surface density of dust plotted against radius, obtained by combining our fitted dust emission models with the temperature model shown in Fig. 4. The solid lines show our best estimates of the dust surface density for our best-fitting dust emission models at 250 and 500 μm, with the shaded regions showing the 1σ uncertainties. The dashed cyan and yellow lines are the best-fitting exponential models to the empirical dust surface-density distributions. The grey vertical dashed line shows the radius below which the data points are not fit to ensure there is no effect from galactic nuclei. The green data point is the one overlapping data point from the quasar reddening study of Ménard et al. (2010a), and the green dotted line is their estimate of the relationship between dust surface density and radius, which varies with radius as r^−0.8 (green line). If Ménard et al. (2010a) does represent the true distribution of dust our disagreement in surface densities may be the result of the range of radii that contribute to the Ménard et al. (2010a) data point or our background subtraction.

Radial variation of other galaxy phases

This work requires comparing our dust profile to the other constituents which make up galaxies. For the molecular gas and total gas we used the results available in the literature (Schruba et al. 2011; Bigiel & Blitz 2012), derived from the HERACLES (Leroy et al. 2009) and THINGS (Walter et al. 2008) samples. This could introduce a potential bias due to the different selection criteria used for HERACLES and THINGS relative to the HRS, but equivalent data are not available for the HRS sample.

To derive the radial distributions of stellar surface density and star-formation rate (SFR) surface density, we used two alternative data sets. Spitzer and GALEX data exist for many of the galaxies (see Table 1), and we used this data to determine empirical average radial distributions of stellar surface density and SFR surface density to compare with the radial distributions of the other galaxy components. However, since 24 μm Spitzer data does not exist for the whole sample and a few of the 3.6 μm images are smaller than ideal, we also estimated profiles using WISE data, which does exist for the entire sample. In the case of the WISE data, we first convolved all the images to the same resolution as the 250 μm image to make the comparison with the dust surface density as fair as possible. We did not do this when using the Spitzer data as the higher resolution profile provides an independent check of the model fitting applied to the Herschel data. Table 1 summarizes the data sets. We obtained very similar results whichever data set we started with.

Table 1.

Exponential scalefactors.

Component	Data	Number	Matched	Gradient
	source	of objects	processing	(dex |$R_{25}^{-1}$|⁠)
Dust	(250 μm and T)	45	Y	−1.72 ± 0.01
Dust	(500 μm and T)	45	Y	−1.75 ± 0.02
Stellar mass	WISE, 3.4 μm	45	Y	−1.78 ± 0.07
Stellar mass	S4G, 3.6 μm	44	N	−1.90 ± 0.03
SFR	WISE 22 μm and NUV	44	Y	−2.69 ± 0.05
SFR	MIPS 24 μm and NUV	28	N	−2.58 ± 0.05
Total Gasa	H i and H2	33	–	−0.71
H2a	–	33	–	−2.17
Metallicitya	–	306	–	−0.16 (with σ = 0.12)

Notes. the gradients from exponential fits to data of different components of a galaxy. The ‘Matched Processing’ column lists whether the data have been convolved to the 250 μm resolution (or already have similar resolution) and have had the same procedure applied as for the dust profiles, or if the model has been fit to high-resolution data. Components marked ^a are gradients taken from the literature (Schruba et al. 2011; Bigiel & Blitz 2012; Sánchez et al. 2014) for different samples of galaxies.

To find the radial variation of the surface density of the stars, we assumed the 3.6 μm emission is proportional to the column density of the stars. We used Spitzer 3.6 μm observations which are available for 44 of the 45 in our sample of large galaxies and were obtained as part of the S⁴G survey (Sheth et al. 2010). Our second data set was the WISE 3.4 μm data which exists for all 45 galaxies in the sample. In the case of this data set, we first convolved the images to the same resolution as the 250 μm images before estimating the average radial distribution of stellar surface density.

To find the radial variation of the SFR surface density, we first created a SFR map for each galaxy by combining an obscured star-formation tracer (either 22 or 24 μm emission) and an un-obscured star formation tracer, the ultraviolet emission, using published recipes (Leroy et al. 2008). Ideally we would have used the far-UV (UV) data from GALEX, but many of the observations are too shallow for this analysis. Instead, we used GALEX near-UV (NUV) images collected for the HRS (Cortese et al. 2012b) that exist for all but one of our sample, and assumed a conversion factor of 0.657 to convert NUV fluxes to FUV fluxes (the conversion factor was calculated from the global UV fluxes for the sample of large galaxies). To estimate the obscured component, we used the Spitzer 24 μm data (Bendo, Galliano & Madden 2012), which exists for about half the sample, and WISE 22 μm data, which exists for the whole sample. In the case of the WISE data, we first convolved the WISE and NUV images to the same resolution as the 250 μm images before producing the SFR maps for each galaxy. We also tested the effect of applying a small correction, using the 3.6 μm emission, to account for the stellar contribution to both the infrared and ultraviolet data (Leroy et al. 2008). This correction was found to have a negligible effect on our conclusions. To summarize, we were able to create high-resolution SFR maps (resolution of ∼5 arcsec) for 28 of our 45 galaxies and maps at similar resolution to the SPIRE 250 μm images for 44 out of 45 galaxies.

Before producing the average distributions of either stellar surface density or SFR surface density, we masked any bright objects on the images, which is especially important due to contamination from stars in these wavebands, and also removed a linear background from the images. In creating the distributions we used very similar methods to the one used to create the radial distributions of submillimetre emission. The smaller sizes of the Spitzer 3.6 μm and 24 μm images meant that our default background annulus was reduced to 2.0–2.5 R₂₅. However, for about half the galaxies the 3.6 μm images were still too small for the background to be determined from this annulus. For these galaxies, we instead estimated the background on each image from an annulus with 1.0 <R/R₂₅ < 1.5 and then added a corrective offset to the profile calculated in the same region from the average profile for the other galaxies.

RESULTS AND DISCUSSION

Dust radial profile

The surface-brightness profile produced from combining the radial profiles from the 117 images (see Section 3.1) obtained with the Herschel SPIRE camera is shown in Fig. 1. Emission from the dust is detected at >3σ to R = 1.975 R₂₅ at 250 μm and at >3σ to R = 1.75 R₂₅ at 350 and 500 μm (the 3σ sensitivity limit is shown by the red-dashed line in the figure). While an increase to R ∼2.0 R₂₅ may not sound a significant improvement over studies that traced dust to R ∼1.2 R₂₅, this corresponds to over a factor of 10 increase in sensitivity.

The average surface-brightness profile for the sample of large galaxies and the best-fitting models to the profile are shown in Fig. 2. We find that at 250 μm the average surface-brightness profile is best-fitted with a broken exponential, but that at 350 and 500 μm a single exponential is sufficient (for details of the models see Section 3.3).

The continuum emission from dust depends both on the column density and temperature of the dust, and so to determine the column density, we determined the average line-of-sight temperature of the dust using the method described in Section 3.3. On the assumption that the dust is heated by starlight, both the light from the old stellar population and the hot young stars in star-formation regions, we have used our stellar surface density, SFR surface density and dust-mass profiles to predict the dust temperature. We used the stellar kinematics including radiative transfer (skirt) code (Camps & Baes 2015), which models the radiative transfer of starlight through dust using a Monte Carlo based approach. The results of the temperature prediction are shown in Fig. 4.

The best-fitting temperature (Fig. 4) declines linearly from ≃24 K in the centre to approximately 15 K at a radius of ≈R₂₅ where the temperature appears to reach a minimum. Our radiative transfer model of the dust in a galaxy predicts a similar decline in the dust temperature with radius but does not predict this plateau in the dust temperature (Fig. 4). We note that fairly high dust temperatures well outside discs have been seen before, with Bocchio et al. (2016) finding that the dust around the edge-on galaxy NGC 891 still has a temperature of ∼19 K 4 kpc above the plane.

Kashiwagi & Suto (2015) estimated that the cosmic UV background (Xu et al. 2005) would produce a minimum dust temperature of ∼10 K. However, given dust radiates with a T⁶ dependence, the UV background would have to be higher by ∼1.5⁶ to explain our minimum temperature of 15 K. Another possible heating source is from collisional heating in a hot X-ray gas. Bocchio (2014) found dust that had been pulled out of NGC 4338 was heated to 20–30 K in the hot X-ray gas in the Virgo Cluster, and Yamada & Kitayama (2005) show that dust could reach ∼30 K in the centres of clusters. However, although spiral galaxies are known to have haloes of hot gas (e.g. Anderson, Churazov & Bregman 2016) it is unclear whether this gas could heat the dust to ∼15 K.

We used our temperature estimates and the models that fit the surface-brightness profiles at the SPIRE wavelengths to determine the radial variation in the column density of the dust (Fig. 5). We find that the dust column density has an exponential decline with radius (gradient ∼−1.7 dex |$R_{25}^{-1}$|⁠), no matter which model of the dust emission we start with. Throughout this paper the gradients (α) are given in units of dex |$R^{-1}_{25}$| (i.e. a value of −1 corresponds to a decrease by a factor of 10 over a radial distance of R₂₅, and can be described by:

\begin{equation} \Sigma = \Sigma _{0} 10^{\alpha r}, \end{equation}

(1)

where Σ is the surface density, Σ₀ is the surface density when r = 0, and r is the radius in angular units. If we extrapolate our exponential surface density model to an infinite radius, we estimate that the dust between 1.0 R₂₅ < R ≤ ∞ constitutes only 9.3 per cent of all the total dust in the disc

Comparison with other galaxy components

Fig. 6 shows all galaxy components plotted on one diagram and Table 1 gives the gradients of the exponential models that best fit the radial distributions of stellar surface density, SFR surface density, dust surface density, molecular gas and total gas. In the case of stellar surface density and SFR surface density, we obtained very similar values for the gradient whether we started with the high-resolution images or the images convolved to the same resolution as the 250 μm images. The blue and green points in Fig. 6 show the empirical average radial distributions of stellar surface density and SFR surface density estimated from the high-resolution data sets. The blue and green dashed lines show the fits of single exponential models to the empirical radial distributions of stellar surface density and SFR surface density derived from the low-resolution data sets. In this case, the empirical distributions are not shown in Fig. 6.

Figure 6.

Average surface density versus radius for all components of a galaxy. The data points show the estimated average surface density for the HRS galaxies in each radial bin for dust (red) star-formation rate (green) and stellar mass (blue). The stellar mass and SFR data points are radial profiles measured from high-resolution data sets (≤6.0 arcsec full-width half maximum) from Spitzer and GALEX. The dashed lines of the same colour show the best-fitting exponential models for images which have been convolved to the same resolution as our 250 μm images and have had the same procedure applied to produce the average profiles. In this case we have used WISE data instead of Spitzer data as observations exist for the entire sample and cover a larger area. In the case of the total gas (black) and molecular gas (magenta) we have shown the exponential distributions derived from samples of local galaxies in the literature (Section 3.4).

As expected, the average surface density of all galaxy components falls with distance from the centre of the galaxy. As expected, given the extended distribution of atomic gas around galaxies (Leroy et al. 2008), the total column density of gas declines much more slowly than the other components. The gradients of the distribution of stellar mass and of dust are in good agreement, but the SFR falls more rapidly with radius than either of the other two.

Hunt et al. (2015) also find similar distributions of the dust and stars, finding the scalelengths of KINGFISH spirals to have a value of R₂₅₀/R₂₅ ≈ 0.37 at 250 μm versus R_3.6/R₂₅ ≈ 0.35 at 3.6 μm (1.17 dex |$R_{25}^{-1}$| and 1.24 dex |$R_{25}^{-1}$|⁠, respectively, in our gradient units). Our results show that the similarity between the distributions of dust and stellar mass extends out to R ∼ 1.8 R₂₅. Our dust gradient of −1.75 dex |$R_{25}^{-1}$| is steeper than the value of 1.17 dex |$R_{25}^{-1}$| for the KINGFISH spiral galaxies (Hunt et al. 2015), although the KINGFISH dust distributions are distributions of dust 250 μm emission rather than dust surface density. Muñoz-Mateos et al. (2009) find the median gradient of the radial distribution of the surface density of dust for the entire SINGS survey of 1.49 dex |$R_{25}^{-1}$|⁠, significantly closer to our value.

It is still uncertain whether dust is mostly formed in evolved stars (Matsuura et al. 2009; Boyer et al. 2012), supernovae (Gomez et al. 2012; Indebetouw et al. 2014) or formed by grain growth in the interstellar medium (ISM; Ossenkopf 1993; Hirashita & Kuo 2011). We note that if dust grains are formed in supernovae and then fairly quickly destroyed, we would expect the surface densities of the dust and the SFR to have the same dependence on radius, which is not found.

Dust, gas and metallicity

Metallicity gradients in galaxies tend to be small. The average metallicity gradient for a sample of 306 galaxies measured with CALIFA is −0.16 dex |$R_{25}^{-1}$| (scatter of ∼0.12 dex |$R_{25}^{-1}$|⁠; Sánchez et al. 2014), from observations of H ii regions out to 1.3 R₂₅ (there is possible evidence of a flattening at larger radii). For the galaxies in the SINGS survey the measured radial metallicity gradients are steeper, ∼−0.42 to −0.33 dex |$R_{25}^{-1}$| (depending on the calibration used; Moustakas et al. 2010). On the assumption that these metallicity measurements apply to all phases of the interstellar medium, we combine the gradient in metallicity from the CALIFA measurement and the gradient in the column density of all gas (Table 1) and find that the gradient in the column density of metals is ∼−0.9 dex |$R_{25}^{-1}$|⁠, again significantly flatter than the dust gradient of −1.7 dex |$R_{25}^{-1}$|⁠. Similar radial declines in the dust-to-metal ratio have been seen in M31 (Mattsson et al. 2014) and in SINGS galaxies (Mattsson & Andersen 2012). In closed-box chemical evolution models this radial decline can be explained by dust grain growth in the ISM (Mattsson, Andersen & Munkhammar 2012).

Are galaxies surrounded by extended haloes of dust?

Based on the angular correlation between the reddening of ∼85 000 high-redshift quasars and the positions of ∼2.4 × 10⁷ galaxies detected in the SDSS (York et al. 2000), Ménard et al. (2010a, henceforth M10) concluded that galaxies are surrounded by very extended haloes of dust. They found that the mean dust surface density varies with angular distance, θ, from the position of an SDSS galaxy as

\begin{equation} \Sigma _d(\theta ) = 2.5 \times 10^{-3} \left({\theta \over 0.1\ {\rm arcmin} }\right)^{-0.8}\ h\ {({\rm M}_{{\odot }}\ {\rm pc}^{-2})}, \end{equation}

(2)

(taken from fig. 8 of M10), with this distribution extending to a physical distance of several Mpc (h is H₀/100). Kashiwagi & Suto (2015) have suggested that the reddening measurements are actually due, not to extended haloes of dust around individual galaxies, to the dust within galactic discs and the fact that the positions of galaxies are correlated. Their evidence for this claim are the results from stacking measurements around the positions of the same SDSS galaxies, using the far-infrared image of the sky made with IRAS. They, too, detect extended dust around the SDSS galaxies, but this time in emission, and by comparing the reddening and emission measurements, they estimate the temperature of the dust is ≃20 K, which they argue is what one would expect for dust in the galactic discs rather than dust well outside galactic discs (however, see Section 4.1).

The strong practical motive for trying to distinguish between these two possibilities is the need to understand the effect of dust on experiments to measure cosmological parameters, such as experiments that use the magnitudes of Type Ia Supernovae to measure Ω_Λ and the equation of state of dark energy (Ménard, Kilbinger & Scranton 2010b). If all galaxies are surrounded by such extended haloes of dust, all measurements of the magnitudes of high-redshift objects need to be corrected for the effects of dust extinction and reddening. But, if the the conclusion of Kashiwagi & Suto (2015) is correct, and almost all dust is confined to galaxy discs, only a small percentage of lines of sight actually pass through galactic discs (we quantify this below). We now discuss what our results can contribute to resolving the controversy. We borrow the nomenclature of Kashiwagi & Suto (2015) in which the first model is referred to as the circumgalactic dust model (CGD model) and the latter model as the interstellar dust model (ISD model).

Our result that the dust column density of galaxies declines more rapidly with radius out to R ≃ 2 R₂₅ than given by equation (2) (Fig. 5) does not allow us to discriminate between the ISD and the CGD models. This is because equation (2) is only valid beyond the minimum radius for which M10 have a reddening measurement which is only slightly less than 2 R₂₅ (Fig. 5), where we are not sensitive. Therefore it is possible that the surface density of dust declines rapidly with radius out to R ≃ 2 R₂₅ but then follows the more gentle relationship given by equation (2).

We now adopt the alternative approach of combining the statistical properties of the SDSS galaxies – their number counts and clustering properties – with the average dust surface density radial profile for the HRS galaxies (Fig. 5) to test whether we can produce a relationship similar to equation (2) from the dust within galaxy discs.

M10 used SDSS galaxies with i < 21, whereas our study is ultimately based on the photometric B-band, since R₂₅ is the distance from the centre of a galaxy out to an isophote with a surface brightness in the B-band of 25 mag arcsec⁻². To use a sample of galaxies as similar as possible to that used by M10 and to connect the sample to the standard definition of optical radius we used the galaxy number counts derived from the SDSS r-band images and an average B − r colour of 1.65 (Yasuda et al. 2001). Since the mean r − i colour measured by Yasuda et al. (2001) is ≃0.4, we used the number counts with r < 21.5.

Of course, not all galaxies detected in the SDSS contain much dust, and in particular early-type galaxies contain small amounts of dust (Smith et al. 2012a). Eales et al. (2015) have recently measured the percentage of the stellar mass in the nearby Universe that is in late-type galaxies is ∼49 per cent, using a definition of ‘late-type’ that means a late-type galaxy is almost certainly dominated by a disc. Informed by this result, we multiply the SDSS number counts by 0.5 to produce an estimate of the number counts of SDSS galaxies which contain a disc, N_disc. We then make the assumption that each of these galaxies has the dust surface density radial profile shown in Fig. 5. Although these are large assumptions, they are not as important as they appear because of the larger uncertainty in the amount of evolution of dust in galaxies, which we discuss below.

On the assumption that the surface brightness of a galaxy depends on radius with the form shown in equation (1), the total flux density out to the standard optical radius is given by

\begin{equation} I_{{\rm tot}} = {2 \pi \, I_{25} \, \over {\left( \alpha \,{\rm log_e(10)} \right) ^2}} \left( \left(\alpha \,R_{25}\, {\rm log_e(10)} - 1 \right) + 10^{-\alpha \,R_{25}} \right), \end{equation}

(3)

in which I₂₅ is the surface brightness corresponding to 25 mag arcsec⁻². Because α is given in units of R₂₅ (Table 1), this equation makes it possible to calculate θ₂₅, and thus the area subtended on the sky, for a galaxy of any apparent magnitude. Pohlen & Trujillo (2006) have fit exponential models to the SDSS images of 90 face-on late-type galaxies. The mean value of α derived from the values given by them for the inner scalelength in the g-band for each galaxy (their table 3) is −1.64 ± 0.08 dex |$R_{25}^{-1}$|⁠. Given the average value of α we can use equation (3) to calculate the optical radius of a galaxy of any apparent magnitude.

A very simple thing we can then do is calculate the fraction of the sky covered by the discs. As we can only detect dust out to a radius of 2 × R₂₅, we assume that a disc only extends out to this radius. Using this assumption, the number counts for late-type galaxies, N_disc, and the average value of α from Pohlen & Trujillo (2006), we estimate that the percentage of the sky covered by discs is 0.44 ± 0.05 per cent, with the errors obtained using the error for the mean value of α given above. This low percentage confirms the importance of distinguishing between the CGD and ISD models for making the correct corrections for dust extinction and reddening when carrying out cosmological experiments, because in the ISD model a very small fraction of the sky is covered by dust.

For a galaxy in which the dust surface density varies with radius with the relationship given in equation (1) the mean column density out to twice the optical radius (2 × θ₂₅) is given by:

\begin{equation} <\Sigma _{{\rm d},{\rm galaxy}}> = {\Sigma _0 \left( 10^{2\alpha \,R_{25}} \left( 2 \alpha \,R_{25}\, {\rm log_e(10)} - 1 \right) + 1\right)\over {2 \left( \alpha \,R_{25}\,{\rm log_e(10)} \right) ^2}}. \end{equation}

(4)

We use the values for Σ₀ and α for the mean dust surface-density profile shown in Fig. 5: Σ₀ = 0.6 M_⊙ pc⁻² and α = −1.75 dex |$R_{25}^{-1}$|⁠. The expected average dust column density profile around SDSS galaxies is then given by:

\begin{equation} \Sigma _{\rm d}(\theta ) = \int _0^{21} <\Sigma _{{\rm d},{\rm galaxy}}> 4 \pi \ (\theta _{25})^2\ w(\theta ,m)\ N_{{\rm disc}}(m)\ {\rm d}m. \end{equation}

(5)

In this equation, w(θ, m) is the angular correlation function for SDSS galaxies. Connolly et al. (2002) have measured the angular correlation function for SDSS galaxies as a function of r magnitude. They found that at all magnitudes the angular correlation function is well represented by w(θ) = Aθ^−γ, with γ ≃ 0.7 and log₁₀A a linear function of r-band magnitude. In calculating w(θ, m), we assume the same value of γ and estimate A at each magnitude from a linear fit to the data shown in fig. 3 of Connolly et al. (2002).

Fig. 7 shows our estimate of Σ_d(θ) from this procedure. The red solid line in the figure shows our prediction with the mean value of α we calculated from the results in Pohlen & Trujillo (2006) and the shaded region shows the predictions if we use values for α between the ±1 σ errors on the mean. The black dashed line shows the relationship measured by M10 (equation 2). This relationship is very similar to the result from the reddening study (M10).

Figure 7.

The relationship between mean dust surface density and angular distance from the centre of an SDSS galaxy. The dashed black line shows the relationship estimated by M10 from looking at the correlation between the reddening of high-redshift quasars and the positions of SDSS galaxies. The solid red line shows our prediction for this relationship considering only the dust in galactic discs, with the shaded region showing the predictions if we use values for α in our model between the upper and lower 1σ errors in the mean value of α (see the text for details). The blue line shows our prediction if we incorporate the evolution in the dust masses of galaxies measured by Bourne et al. (2012).

This is without considering the effect of galactic evolution, as the mean redshift of the SDSS galaxies is ≃0.3. Two results, both from the Herschel-ATLAS survey (Eales et al. 2010), have shown that there is substantial evolution in the amount of dust in galaxies by this redshift. Dunne et al. (2011) have shown that the dust masses of galaxies have increased by a factor of ∼4.6 by this redshift. This result is based on the small fraction of galaxies with very high dust masses, but Bourne et al. (2012) have carried out a stacking analysis on a large sample of SDSS galaxies which is sensitive to galaxies with much lower dust masses, finding that the average dust mass of optically selected blue galaxies of the same stellar mass has increased by a factor of ≃3 by a redshift of 0.3. In Fig. 7 we have shown our estimate of Σ_d(θ) scaled by a factor of 3.0 to allow for the effect of galaxy evolution. This now exceeds the relationship found in the reddening study.

An important consideration is the absolute and relative uncertainties in the values of κ_ν used for both our analysis of dust emission and the reddening study of M10. Our value of κ_ν = 0.192 m² kg^{− 1} at 350 μm (Draine 2003) assumes Milky Way like carbonaceous and silicate dust, and is in the typical range of values used of κ_ν in local galaxies studies (see Clark et al. 2016, and references therein). Our κ_ν value is higher by a factor of ∼2 than the latest empirical value derived in this study (Clark et al. 2016); using this value would roughly double the amount of dust we measure in the discs of galaxies. As we are only measuring the dust emission in the discs of galaxies these values should be reasonable values to assume out to 2.0 R₂₅, particularly as the metallicity gradients of galaxies tends to be small. A greater unknown is what the extinction coefficient would be for dust in the IGM; M10 assume an SMC-type dust due to the lack of a 0.2 μm bump in low-ionization absorbers in the haloes of ∼L* galaxies (e.g. York et al. 2006) and the fact that extinction curves for high-redshift galaxies also do not show the 0.2 μm bump. If larger Milky Way-type dust were assumed, the halo dust surface density implied by the M10 results would approximately double. Of course, the dust in intergalactic space might be very unlike the dust within galactic discs. Sputtering is likely to be important (Yamada & Kitayama 2005) and possibly only large grains would survive being expelled from galactic discs (Croft et al. 2000). We emphasize, however, that in our model 90 per cent of the dust contributing to the average surface density shown in Fig. 7 is produced by dust within a radius less than R₂₅, where our knowledge of the properties of dust is much better.

Of course there are also large uncertainties in our model, including the division between early-type and late-types, the assumption that all late-type galaxies have the same distribution of dust, and the amount of galactic evolution. Because of these uncertainties and because of the huge uncertainty as to whether the dust in intergalactic space is like dust in galactic discs, we cannot rule out the CGD model. All we can say is that with some assumptions about the properties of dust (e.g. Milky Way-type dust grains), the results of M10 might be explained by the dust within individual galaxies rather than very extended haloes of dust around all galaxies.

Despite using results from most sensitive ever far-infrared space telescope, and then combining the results in a powerful statistical stacking analysis, we have still only managed to detect the continuum emission from dust out to a distance of twice R₂₅ from the galactic centre. Therefore the only practical way to investigate the properties of dust at greater distances from centres of galaxies is still through statistical investigations of the reddening of high-redshift objects. Over three decades, investigations of the dust associated with damped Lyα quasar absorption-line systems have produced many significant detections of this dust and have even produced detailed measurements of the properties of this dust (e.g. Fall, Pei & McMahon 1989; Khare et al. 2004; York et al. 2006; Wild, Hewett & Pettini 2007; Vladilo, Prochaska & Wolfe 2008). There have also been detections through quasar absorption-line of gas and metals in the haloes of intervening individual galaxies, with projected distances of up to 100 kpc (e.g. Steidel et al. 2002; Chen, Kennicutt & Rauch 2005), although some galaxies fail to show any absorption (Bechtold & Ellingson 1992). Investigations of the correlation of quasar absorption-line systems with galaxy redshift surveys imply that the absorption-line gas is often a long way from the centre of a galaxy (Churchill, Steidel & Kacprzak 2005 and references therein; Pérez-Ràfols et al. 2015), which suggests that much of the dust detected via reddening of quasars is also outside galactic discs. The method used by M10 is an elegant way of both measuring the properties of the dust and also its distribution around galaxies. We have shown that there could be a significant contribution to the reddening measured by M10 from dust within galactic discs.

Fortunately, there is actually a simple way of adapting the M10 technique to allow for this contamination. Even if the CGD and ISD models predict very similar relationships for how the mean reddening depends on distance from an SDSS galaxy, in the ISD case this relationship is caused by the small percentage of lines of sight that pass through galactic discs, whereas in the CGD case it is caused by all lines of sight. Therefore, a simple way to distinguish between the two hypotheses would be to measure the variance of the quasar reddening measurements, which should be much higher for the ISD model than the CGD model – a test that should be easy to carry out using the measurements described by M10.

CONCLUSIONS

In this paper we present the results of combining images of a large sample of nearby galaxies made from observations with the Herschel Space Observatory. By combining the images, we determine how the average surface density of dust depends on radius. We find the following results:

• We detect the emission of dust at >3σ threshold out to a radius of 2.0 R₂₅ at 250 μm and out to a radius of 1.75 R₂₅ at 350 and 500 μm, by combining images of 117 galaxies from the HRS. When we restrict the sample to the 45 galaxies with R₂₅ >1.5 arcmin we can trace dust at 250 μm out to a radius of 1.75 R₂₅ at > 3σ.

• We fit modified blackbody spectral-energy distributions to the SPIRE and PACS (100–500 μm) radial profiles. We find the dust temperature declines from ∼25 K in the centre of the galaxy to ∼15 K at 1.0 R₂₅ where it remains constant out to 1.75 R₂₅.

• We find the radial dust-mass distribution in a galaxy is consistent with an exponential profile out to a radius of 1.8 R₂₅. The gradient of the surface-density profile is −1.7 dex |$R^{-1}_{25}$|⁠.

• We created profiles of the stellar mass (traced by 3.4 or 3.6 μm) and SFR (UV and Infrared) using images from Spitzer, GALEX and WISE. We found the stellar gradient is in good agreement with the dust gradient with the dust gradient between 1 per cent and 10 per cent larger (in dex |$R^{-1}_{25}$|⁠) depending on the data set used. The SFR profile is significantly steeper than that of the dust with a gradient of ∼−2.6 dex |$R^{-1}_{25}$|⁠.

• Our results show a much steeper gradient in the surface density of dust compared to metals (based on typical metallicity gradients for nearby galaxies), similar to the results for individual galaxies (e.g. Mattsson & Andersen 2012). In closed-box chemical evolution models this is a signature that dust originates from grain growth in the ISM.

• We use our results for the radial distribution of dust in galaxies to see if the large extended haloes of dust around galaxies as described by M10 could have a significant contamination caused by the clustering of galaxies and the dust emission within individual galaxies. By using the number counts and angular correlation function from SDSS, we show that, assuming Milky Way-type dust, the dust within galaxy discs might explain the reddening results. Because of the large uncertainties in our analysis and the even larger uncertainty of whether the dust in intergalactic space is like dust in galactic discs, we cannot rule out the existence of intergalactic dust. We propose a simple method of testing whether there is actually intergalactic dust.

Our results from a statistical stacking analysis with images from the Herschel Space Observatory represent the most sensitive study of the extended dust emission around galaxies that will be possible until a new submillimetre space mission. It is likely to be impossible to improve on our results with ground-based telescopes because of the fluctuating emission from the atmosphere, which makes it difficult to detect low surface-brightness emission from galaxies.

We thank everyone involved with the Herschel Space Observatory.

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including: University of Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK); and Caltech, JPL, NHSC, University of 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).

HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center and the HIFI, PACS, and SPIRE consortia.

This research made use of Astropy, a community-developed core python package for Astronomy (Astropy Collaboration, 2013).

HLG acknowledges support from the European Research Council (ERC) in the form of Consolidator Grant CosmicDust (ERC-2014-CoG-647939).