GAMA/H-ATLAS: the local dust mass function and cosmic density as a function of galaxy type – a benchmark for models of galaxy evolution Free

ABSTRACT

1 INTRODUCTION

Cosmic dust is a significant, albeit small, component of the interstellar medium (ISM) of galaxies. Despite being less than 1 per cent of the baryonic mass of a galaxy, dust is responsible for obscuring the ultraviolet and optical light from stars and active galactic nuclei and is thought to have absorbed approximately half of the starlight emitted since the big bang (Puget et al. 1996; Fixsen et al. 1998; Dole et al. 2006; Driver et al. 2016). Measuring the dust mass in galaxies is therefore important for understanding obscured star formation (Kennicutt 1998; Calzetti et al. 2007; Marchetti et al. 2016), particularly at different cosmic epochs (Madau, Pozzetti & Dickinson 1998; Hopkins 2004; Takeuchi, Buat & Burgarella 2005). The dust mass function (DMF) is one of the fundamental measurements of the dust content of galaxies, providing crucial information on the reservoir of metals that are locked up in dust grains (Issa, MacLaren & Wolfendale 1990; Edmunds 2001; Dunne, Eales & Edmunds 2003). A measure of the space density of dusty galaxies is becoming even more relevant given the widespread use of dust emission as a tracer for the gas in recent years (Eales et al. 2010, 2012; Magdis et al. 2012; Scoville et al. 2014, 2017; see also the comprehensive review of Casey, Narayanan & Cooray 2014). This is of particular interest given difficulties in observing atomic and molecular-line gas mass tracers out to higher redshifts (Tacconi et al. 2013; Catinella & Cortese 2015).

Ground-based studies including observations at 450 and 850 |$\mu$|m with the Submillimetre Common User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope led to the first measurements of the DMF over the mass range ∼10⁷ M_⊙ < M_d < few × 10⁸ M_⊙ (Dunne et al. 2000; Dunne & Eales 2001; Vlahakis, Dunne & Eales 2005), where M_d is dust mass. Unfortunately the state of the art at that time meant fewer than 200 nearby galaxies were observed with small fields of view and selected at optical or infrared (60|$\mu$|m) wavelengths. At higher redshifts, the Balloon-borne Large Aperture Submillimeter Telescope (BLAST, observing at 250–500 |$\mu$|m) enabled a DMF to be derived out to |$z$| = 1 (Eales et al. 2009) and a valiant effort to measure at even higher redshifts (⁠|$z$| = 2.5) using SCUBA surveys was attempted by Dunne et al. (2003). These studies were hampered by small number statistics and difficulties with observing from the ground.

The advent of the Herschel Space Observatory (hereafter Herschel, Pilbratt et al. 2010) and Planck Satellite revolutionized studies of dust in galaxies, as they enabled greater statistics, better sensitivity, and angular resolution in some regimes, wider wavelength coverage and the ability to observe orders of magnitude larger areas of the sky than possible before. The largest DMF of galaxies using Herschel was presented in Dunne et al. (2011) consisting of 1867 sources out to redshift |$z$| = 0.5, selected from the Science Demonstration Phase (SDP) of the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS) blind 250-|$\mu$|m fields (Eales et al. 2010, 16  deg²). Their DMF extended down to 5 × 10⁵ M_⊙ and they derived a redshift dependent dust mass density of Ω_d = ρ_d/ρ_crit = (0.7–2) × 10⁻⁶. Subsequently, Negrello et al. (2013) and Clemens et al. (2013) published the DMF of 234 local star-forming galaxies from the all sky Planck catalogue. Clark et al. (2015) then derived a local DMF from a 250-|$\mu$|m selected sample consisting of 42 sources. These DMFs ranged from 10⁶ M_⊙ < M_d < few × 10⁸ M_⊙ and 2 × 10⁵ M_⊙ < M_d < 10⁸ M_⊙, respectively. These measurements¹ were found to be consistent with the |$z$| = 0 estimate from Dunne et al. (2011), once scaled to the same dust properties, as well as those derived from optical obscuration studies using the Millennium Galaxy Catalogue (Driver et al. 2007).

Interestingly, although the dust mass density is broadly consistent across most surveys, the shape of the DMF differs between all of these different estimates. Clark et al. (2015) demonstrated using a blind survey selected at 250 |$\mu$|m, around a third of the dust mass in the local Universe is contained within galaxies that are low stellar mass, gas-rich, and have very blue optical colours. These galaxies were shown to have colder dust populations on average (⁠|$1 2\,\rm K\lt \mathit{ T}_{\rm d} \lt 16\,\rm K$|⁠, where T_d is the cold-component dust temperature) compared to other Herschel studies of nearby galaxies, e.g. the Herschel Reference Survey (Boselli et al. 2010), the Dwarf Galaxy Survey (Madden et al. 2013; Rémy-Ruyer et al. 2013, see also De Vis et al. 2017b), and higher stellar mass H-ATLAS galaxies (Smith et al. 2012a). This led to higher numbers of galaxies in the low dust mass regime than predicted from extrapolating the Dunne et al. (2011) DMF down to the equivalent mass bins (Clark et al. 2015).

In comparison, the Clemens et al. (2013) and Vlahakis et al. (2005) DMFs are in reasonable agreement and both suggest a low-mass slope that is much steeper than the Dunne et al. (2011) function. Overall, comparing between these different measures is complex due to different selection effects; furthermore they are limited due to (i) small number statistics, and/or (ii) lack of sky coverage or volume, inflating uncertainties due to cosmic variance. We also show evidence in Section 4 that fitting the same dataset over different mass ranges can have a significant effect on the resulting best-fitting parameters. Since we probe further down the low-mass end than any literature study, this could therefore have a significant impact.

Here we further the study of the DMF by deriving the ‘local’ (⁠|$z$| < 0.1) DMF for the largest sample of galaxies to date, the sample is taken from the Galaxy and Mass Assembly Catalogue (GAMA, Driver et al. 2011). The large size of this sample reduces the statistical uncertainties and the effect of cosmic variance. We also employ statistical techniques to address selection effects in our sample, which allows us to probe further down the DMF by at least an order of magnitude compared to previous works. We present the observations and sample selection in Section 2 and the method used to derive the dust masses for the GAMA sources in Section 3. The DMF is presented in Section 4 and is compared to predictions from semi-analytical models in Section 5. In Sections 6 and 7, we split the DMF by morphological type and compare with their corresponding stellar mass functions, with conclusions in Section 8. Properties of the full GAMA sample are discussed in detail in Driver et al. (2018) and the accompanying stellar mass function of the same sample is published in Wright et al. (2017, hereafter W17). Throughout this work we use a cosmology of Ω_m = 0.3, Ω_Λ = 0.7, and |$H_0 = 70\,\rm km\,s^{-1}\,Mpc^{-1}$|⁠.

2 OBSERVATIONS AND PHOTOMETRY

2.1 GAMA

The GAMA² survey is a panchromatic compilation of galaxies built upon a highly complete magnitude limited spectroscopic survey of around 286 deg² of sky (with limiting magnitude r_petro ≤ 19.8 mag as measured by the Sloan Digital Sky Survey (SDSS) DR7, Abazajian et al. 2009). Around 238 000 objects have been successfully observed with the AAOmega Spectrograph on the Anglo-Australian Telescope as part of the GAMA survey. As well as spectrographic observations, GAMA has collated broad-band photometric measurements in up to 21 filters for each source from ultraviolet (UV) to far-infrared (FIR)/submillimetre (submm) (Driver et al. 2016; Wright et al. 2017). The imaging data required to derive photometric measurements come from the compilation of many other surveys: GALEX Medium Imaging Survey (Bianchi & GALEX Team 1999); the SDSS DR7 (Abazajian et al. 2009), the VST Kilo-degree Survey (VST KiDS, de Jong et al. 2013); the VIsta Kilo-degree INfrared Galaxy survey (VIKING, de Jong et al. 2013); the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010); and the Herschel-ATLAS (Eales et al. 2010). The motivation and science case for GAMA is detailed in Driver et al. (2009). The GAMA input catalogue definition is described in Baldry et al. (2010), and the tiling algorithm in Robotham et al. (2010). The data reduction and spectroscopic analysis can be found in Hopkins et al. (2013). An overview and the survey procedures for the first data release (DR1) are presented in Driver et al. (2011). The second data release (DR2) was nearly twice the size of the first and is described in Liske et al. (2015). Information on data release 3 (DR3) can now be found in Baldry et al. (2018). There is now a vast wealth of data products available for the GAMA survey, making it an incredibly powerful database for all kinds of extragalactic astronomy and cosmology.

K-corrections for GAMA sources are available from Loveday et al. (2012) using k-correct v4_2 (Blanton & Roweis 2007). Redshifts derived using autoz are available from Baldry et al. (2014). This work consists of data from the GAMA equatorial fields, which has a redshift completeness of >98 per cent at r_petro ≤ 19.8 mag (Liske et al. 2015). GAMA distances were calculated using spectroscopic redshifts and corrected (Baldry et al. 2012) to account for bulk deviations from the Hubble flow (Tonry et al. 2000).

For this paper, we select galaxies in the three equatorial fields of the GAMA survey, which cover ∼180 deg² of sky between them. The equatorial fields G09, G12, and G15 are located on the celestial equator at roughly 9, 12, and 15 h, respectively. We use the redshift range 0.002 ≤ |$z$| ≤ 0.1, with the upper limit matching the low |$z$| bin from the earlier DMF study of Dunne et al. (2011); this redshift range contains 20 387 galaxies (with spectroscopic redshift quality set at n_Qual ≥ 3).³ These GAMA galaxies have been further split into early types (ETGs), late types (LTGs), and little blue spheroids (LBSs) based on classifications using giH-band images from SDSS (York et al. 2000), VIKING (Sutherland et al. 2015), or UKIDSS-LAS (see Kelvin et al. 2014; Moffett et al. 2016a for more details on the classification).

2.2 Herschel-ATLAS

The FIR and submm imaging data, which are necessary to derive dust masses, are provided via H-ATLAS⁴ (Eales et al. 2010), the largest extragalactic Open Time survey using Herschel. This survey spans ∼660 deg² of sky and consists of over 600 h of observations in parallel mode across five bands (100 and 160|$\mu$|m with PACS – Poglitsch et al. 2010, and 250, 350, and 500 |$\mu$|m with SPIRE – Griffin et al. 2010). H-ATLAS was specifically designed to overlap with other large area surveys such as SDSS and GAMA. The GAMA/H-ATLAS overlap covers around 145  deg² over the three equatorial GAMA fields, G09, G12, and G15. Photometry in the five bands for the H-ATLAS DR1 is provided in Valiante et al. (2016) based on sources selected initially at 250 |$\mu$|m using madx (Maddox et al. in preparation) and having S/N > 4 in any of the three SPIRE bands. Bourne et al. (2016) present optical counterparts to the H-ATLAS sources, identified from the GAMA catalogue using a likelihood ratio technique (Smith et al. 2011). In this paper, we use the aperture-matched photometry from Herschel based on the GAMA r-band aperture definitions using the lambdar package (Wright et al. 2016), this method is described briefly in Section 2.3.

Given the requirement for H-ATLAS and GAMA coverage, the final sample for this work consists of 15 951 galaxies, this number includes a selection on r_petro ≤ 19.8 and the fact that due to the shapes of the H-ATLAS and GAMA fields, some of the GAMA sources were not covered by Herschel.

2.3 Photometry with lambdar

The Lambda Adaptive Multi-band Deblending Algorithm in R (lambdar)⁵ is an aperture photometry package developed by Wright et al. (2016), which performs photometry based on an input catalogue of sources. Aperture-matched photometry can be implemented on any number of bands and for each band the apertures are convolved by the PSF of the instrument. lambdar also deblends sources occupying the same on-sky area, this is achieved by sharing the flux in each pixel between all overlapping apertures. The fractional splitting is done iteratively and, depending on user preference, can be based on the mean surface brightness of a source, central pixel flux, or a user-defined weighting system. Each source is considered in a postage stamp of the input image focused on the source, the size of which depends upon the size of the aperture itself. All known sources within the postage stamp are deblended, including an optional list of known contaminants specified by the user. For this paper this includes H-ATLAS detected sources from Valiante et al. (2016) which do not have a reliable optical counterpart. These are assumed to be higher redshift background sources.

The sky estimate for each source is calculated by randomly placing blank apertures with dimensions equal to the object aperture on the postage stamp, using the number of masked pixels in each blank aperture to weight its contribution to the background estimate. Furthermore, during flux iteration, if any component of a blend is assigned a negative flux then it is rejected for all subsequent iterations (and any negative measurement is set to zero). There are a very small number of sources which end up with negative fluxes at the final iteration and, for consistency, the lambdar pipeline sets these to zero also. For the purposes of this work, we note that the fluxes for 11 210 (70.3 per cent) sources are not above the 3 σ level at 250 m; however, even galaxies which fall below 3 σ do have a valid measurement and error estimate in five Herschel bands and thus provide information for deriving dust masses. We discuss potential biases and tests in later sections. For further details on the lambdar software and data release see Wright et al. (2016).

3 DERIVING GALAXY PROPERTIES WITH magphys

For each galaxy we take the dust and stellar properties from Driver et al. (2018), who used the magphys⁶ package (da Cunha, Charlot & Elbaz 2008) to fit model spectral energy distributions (SEDs) to the 21-band lambdar photometry. magphys uses libraries containing 50 000 of model SEDs covering both the UV-NIR (Bruzual & Charlot 2003) and MIR/FIR (Charlot & Fall 2000) components of a galaxy’s SED along with a χ-squared minimization technique to determine physical properties of a galaxy, including stellar mass, dust mass, and dust temperature. magphys imposes energy balance between these components, so that the power absorbed from the UV-NIR matches the power re-radiated in the MIR/FIR. In the FIR-submm regime, two major dust components are included in the libraries: a warm component ( 30–60 K) associated with stellar birth clouds; and a cold dust component (  15–25 K) associated with the diffuse ISM. A dust mass absorption coefficient of |$\kappa _{850} = 0.077\, \rm m^{2}kg^{-1}$| is assumed, with an emissivity index of β = 1.5 for the warm dust, and β = 2 for cold dust, where κ_λ ∝ λ^−β. This is consistent with the κ values derived from observations of nearby galaxies (James et al. 2002; Clark et al. 2016, see also Dunne et al. 2000) and ∼2.4 times higher than the oft-used Draine (2003) theoretical values (based on their κ scaled to 850 |$\mu$|m with β = 2).⁷ Using the latest values for κ in the diffuse ISM of the Milky Way from Planck Collaboration XXIX (2016) would give dust masses 1.6 times higher than quoted here. For each galaxy magphys uses all of the lambdar measurements to find the best-fitting combination of optical and FIR model SEDs, and outputs the physical parameters for this combined SED. We do not apply any signal-to-noise cuts, but low signal-to-noise measurements clearly do not contribute strong constraints in the fitting. So long as the estimated fluxes and uncertainties are unbiased, this makes maximum use of the information available. magphys also generates a ‘probability distribution function’ (PDF) for each parameter by summing |$\mathrm{ e}^{-\chi ^2/2}$| over all models. The PDF for each parameter is used to determine the acceptable range of the physical quantity, expressed as percentiles of the probability distribution of model values. The results from magphys for the GAMA equatorial regions are presented in Driver et al. (2016), and we use them throughout this work. For our analysis we use the median value for each parameter, because this is more robust than the estimate from the best-fitting model combination. Where uncertainties are required, we use the 16th and 84th percentiles, which correspond to a 1 σ uncertainty for a Gaussian error distribution.

Our version of magphys is slightly modified compared to the default distribution available online. We use the most up-to-date estimates of the Herschel band-pass profiles for both the PACS and SPIRE instruments. Also in our version, the model photometry for each of the Herschel pass bands is calibrated to the nominal central wavelength of each band, as described in the SPIRE Handbook⁸ (Griffin et al. 2010, 2013), rather than the effective wavelength, which is the case for other photometry. Running the code with and without these changes does not highlight any systematic error in the FIR-based magphys output; however, it does change individual measurements by up to a few per cent.

A large fraction of the GAMA sources have measurements with signal-to-noise ratio below 3 σ in the FIR bands: for the |$z$| < 0.1 sample that we use here 32 per cent have fluxes >3 σ. Given that lambdar assigns a zero flux for each blend component that returns a negative flux at any iteration, the error distribution of faint sources becomes one-sided. If we assume that the errors are Gaussian and consider sources which have a true flux much less than σ, then the bias introduced is the mean value of the positive half of a Gaussian i.e. |$\sigma /\sqrt{2{\rm \pi} } \approx 0.4\,\sigma$|⁠. Sources with more positive fluxes will have a smaller bias.

3.1 Temperatures

The normalized distribution of dust temperatures output by magphys for lambdar sources with fluxes above 3 σ in one, two or three Herschel-SPIRE bands is shown in Fig. 1 (top panel). Where we have sources with Herschel fluxes >3 σ in one or more bands, the temperature is well constrained (±∼1 K), and has a tendency to be fairly cold, ∼18 K. There is also a tendency for the galaxies with Herschel fluxes >3 σ in all three bands to be colder than those with only one or two bands; this is not unexpected given that the combination of the shape of the SED of a modified blackbody, and the more sensitive bluer SPIRE bands. The temperature histogram for these sources appears to continue to rise at temperatures below 17 K, with a peak at 16 K. This potentially suggests that a colder dust prior than the  15−25 K used in this work might be needed for a small fraction of galaxies (e.g. Smith et al. 2012a; Viaene et al. 2014; De Vis et al. 2017a). We will return to this below.

Figure 1.

Top: The normalized distribution of the cold ISM dust temperature for the low-redshift sample (⁠|$z$| ≤ 0.1). The red, blue, and green histograms show galaxies with >3 σ fluxes in one, two, or three SPIRE bands, respectively. Each histogram is normalized to a total count of one: the fraction of sources in each histogram is 32, 17, and 6 per cent, respectively. Bottom: The distribution of uncertainties on the dust mass estimates. The uncertainties are calculated as half the difference between the 84th and 16th percentiles of the PDF; if the uncertainties are Gaussian, they correspond to 1σ. The black, red, blue, and green histograms show galaxies with >3 σ flux measurements in zero, one, two, or three SPIRE bands, respectively.

Open in new tabDownload slide

For the galaxies that have fluxes below 3 σ in all of the Herschel SPIRE bands, we have poor constraints on the cold dust temperature. For these galaxies, the temperature PDF follows the underlying flat temperature prior used in the magphys code with limits from  15 to 25 K. Since the temperature estimate is the median of the PDF, this tends towards the median of the prior as the constraints become weaker.⁹ Despite this, the combination of UV and optical photometry and the FIR measurements do provide useful information on the dust masses for those galaxies with FIR fluxes <3 σ in all Herschel bands. This can be seen in Fig. 1 (bottom panel), which shows the distribution of estimated dust mass uncertainties for galaxies with >3 σ in zero, one, two, or three SPIRE bands. For the subsets in one, two, or three bands, the corresponding uncertainties in mass are 0.18, 0.14, and 0.1 dex. Galaxies with <3 σ in any SPIRE band typically have dust mass uncertainties of 0.4 dex on average.

As a further, though indirect, check that the estimated uncertainties are reasonable we look at the distribution of dust mass and stellar mass of the GAMA |$z$| ≤ 0.1 sample, as shown in Fig. 2. The sources with fluxes >3σ in at least one band are shown in green (as expected, these are the more dusty galaxies), with the entire sample shown by the grey bins. We see that the distribution shows a marked bimodality in this plane, clearly visible even for sources without fluxes >3σ in any of the FIR bands. To investigate this further, Fig. 2 highlights the morphological classifications of the galaxies, split into ETGs and LTGs (Moffett et al. 2016a).¹⁰ The ETGs have many fewer >3 σ sources than the LTGs, even for bright optical sources, and this is as expected given that ETGs contain an order of magnitude less dust than late-type galaxies of the same stellar mass (see e.g. Bregman et al. 1998; Clemens et al. 2010; Skibba et al. 2011; Rowlands et al. 2012; Smith et al. 2012b; Agius et al. 2013, 2015). If the true uncertainties in M_d were larger than 0.5 dex, the bimodal structure in Fig. 2 would be smeared out, suggesting the errors in magphys do reasonably represent the uncertainties.

Figure 2.

The distribution of dust mass and stellar mass in GAMA galaxies. The black underlying points show the whole low-redshift (⁠|$z$| ≤ 0.1) sample. The green points show galaxies with >3 σ fluxes in one or more SPIRE bands. Contours show the demarcation into ETGs (black/red contours) and LTGs (black/blue contours) – see the text for details.

Open in new tabDownload slide

3.2 Dust masses and the temperature prior

The cold dust temperature prior is clearly going to impose some limits on the dust mass uncertainty from the fits. However, we argue that the prior temperature range from magphys used in this work is appropriate for a number of reasons. (i) A range of cold dust temperatures between  15 and 25 K is in fact a good description of the observed range of cold dust temperatures in galaxies (Dunne & Eales 2001; Skibba et al. 2011; Smith et al. 2012c; Clemens et al. 2013; Clark et al. 2015). (ii) Smith et al. (2012a, Appendix A) investigated whether a broader temperature prior should be used in magphys fitting. They found that changing the prior range suggested that only 6 per cent of their Herschel detected sources were actually colder than 15 K. They also demonstrated that adopting a wider temperature prior is not always appropriate given the non-linear increase in dust mass when the temperature falls below 15 K (where the SPIRE bands are no longer all on the Rayleigh–Jeans tail). At T < 15 K, symmetrical errors in the fitted temperature produce a very skewed PDF for the dust mass and result in a population bias to higher dust masses for a distribution of Gaussian errors in cold temperature. Furthermore, in relation to SED fitting, a very cold dust component contributes very little to the luminosity in the FIR per unit mass, so it can be included by a fitting routine with very little penalty in χ² when the photometry in the FIR and submm is of low SNR. Indeed Smith et al. (2012a) use simulated photometry to show that galaxy dust masses can be overestimated by (in excess of) 0.5 dex when widening the prior to below 15 K; they therefore strongly caution in using wider temperature priors for sources with weak submm constraints (as is the case here). (iii) Though some galaxies have been shown to require colder dust temperatures than 15 K (Viaene et al. 2014; Clark et al. 2015; De Vis et al. 2017a; Dunne et al. accepted), the fraction of our sample with >3 σ in at least one band that have dust temperatures |$\lt 16\,\rm K$| is <9 per cent.

As an example to illustrate the potential size of the effect, consider the case that 6 per cent of our galaxies had a true dust temperature of 12 K but instead we fit a temperature of 15 K due to the limited prior. We would underestimate the dust mass for this population by a factor of ∼2.6 (i.e. 0.4 dex). However, 94 per cent of galaxies have true temperatures in the range 15–25 K and since most of them do not have >3σ FIR fluxes they will have errors on the fitted temperature of order ±5 K. Widening the prior to extend to 12 K would mean that 16 per cent of sources would be erroneously returned a temperature which was below 15 K resulting in a large positive bias to their dust masses.

Appendix  A presents a more thorough investigation of the effects on the DMF that result from poorly constrained cold dust temperatures for galaxies with low signal to noise in the FIR.

4 THE DUST MASS FUNCTION

4.1 Volume estimators

We use two methods toestimate the accessible volume for each galaxy. First we derive V_max for each galaxy by estimating the maximum redshift at which that source would still be visible given the limiting magnitude of the survey. This requires taking into account both the optical brightness of each galaxy and the K-correction required as the galaxy SED is redshifted. The maximum redshift is not allowed to exceed the user-imposed redshift range of the sample (here we use |$z$| < 0.1). Using this maximum redshift and the area of the survey, an accessible comoving volume can be calculated. These maximum volumes are the same as used in W17. We refer to this method as pV_max, since it is based on the simple photometric selection of the survey.

The second method we use to estimate the V_max for each galaxy is based on a bivariate brightness distribution (BBD). This involves binning the data in terms of the two most prominent selection criteria, and aims to account for the selection effects that they introduce. Since our sample is optically selected, we choose the absolute r-band magnitude, and for the second axis we choose surface brightness in the r band (Loveday et al. 2012, 2015). We have estimated fluxes in all other bands for all galaxies, even if they are not significantly detected, so we do not directly apply any further selection criteria.

This method follows closely the format of the galaxy stellar mass function (GSMF) produced by W17 for the same sample; see also Fig. 3, which is a diagrammatic representation of the BBD method. For each 2D r-band magnitude/surface brightness bin (Fig. 3a), the volume enclosed by the median luminosity distance of the galaxies in the bin and the on-sky area of GAMA is calculated (Fig. 3b) and doubled in order to find an ‘accessible volume’ for all of the galaxies in that bin (Fig. 3c). Using twice the median value will provide an effective V_max that, at some level, corrects for the incompleteness at large distances whatever the cause of the incompleteness. Thus the BBD method has the benefit that it can correct for selection effects in two parameters at once. Using the median volume to determine the effective V_max has the advantage that it is more statistically robust than the actual maximum volume observed in a given bin. However, this estimator is only strictly valid when the underlying galaxy distribution in any given bin is randomly and evenly distributed in space, so the average V/V_max = 0.5. Given the large density fluctuations seen in the galaxy distribution, we cannot state that it is always the case, particularly for local, low-mass galaxies, which are hampered by small-number statistics and strongly affected by cosmic variance. It is more likely to be the case that |$V/V^{\prime }_{\rm max} = 0.5$|⁠, i.e. the maximum volume weighted by density. To allow for these density fluctuations, we find a median weighted by the inverse of the density correction factors δ_n/⟨δ_f⟩, defined below. Galaxies in overdense regions are given less weight in the median compared to galaxies in underdense regions, so any bias in the median volume from density fluctuations should be minimized. We note that in order to reduce noise introduced into the DMF from BBD bins with poor statistics we perform a Monte Carlo (MC) simulation whereby we perturb the quantities used for the two ‘axes’ of our BBD within their associated uncertainties and recalculate the BBD 1000 times and find the median BBD V_max associated with each bin. In essence, this smooths the BBD by the estimated errors, and reduces the uncertainty in the BBD V_max.

Figure 3.

The BBD for our sample with surface brightness and r-band magnitude as the two ‘axes’ (W17) with (a) raw counts in surface brightness/r-band magnitude bins, (b) median volume in surface brightness/r-band magnitude bins, (c) Weighted counts, i.e. volume density in the surface brightness/r-band bins. Each of the panels represents the BBD resulting from the median of 1000 Monte Carlo simulations where we perturb the r-band magnitude and surface brightness within their associated uncertainties.

Open in new tabDownload slide

A direct comparison of the maximum volumes derived from both the pV_max and BBD methods is shown in Fig. 4 with the points coloured by the average number of galaxies in the BBD bin containing that galaxy across all the MC simulations. The largest deviation from the 1:1 line is seen for galaxies that lie in bins with a small number of galaxies contributing to the median volume. These volumes are generally low, meaning they are also strongly affected by cosmic variance. The pV_max values are systematically higher by 0.8 per cent on average than those derived from the BBD method, which translates to an average offset of 1 per cent in the binned DMF values when determined by the median weighted by the error on the measurement.

Figure 4.

The maximum effective volumes for our galaxies at |$z$| < 0.1 derived using the pV_max method (x-axis), and BBD method using r-band magnitude and surface brightness as the two selection features (y-axis). The colour of the points is determined by the number of galaxies in the BBD bin that each galaxy resides in (Fig.3), as shown by the colour bar in the top left corner. We note that the number of galaxies per bin is the median resulting from 1000 Monte Carlo simulations, where we perturb the r-band magnitude and surface brightness within their associated uncertainties.

Open in new tabDownload slide

Since we compare to the GSMF from W17, who use stellar mass and surface brightness as the BBD axes, it may be argued we should use the same approach. We consider this in Appendix  B, and conclude that the Schechter parameters are consistent with the r-band and surface brightness BBDs within uncertainties. We opt to use the r-band magnitude for our second axis here as it is more in line with the optical pV_max, and does not depend on stellar properties directly.

Density fluctuations in the GAMA equatorial fields e.g. Driver et al. (2011) and Dunne et al. (2011) have a pronounced effect on the DMF, and so we apply density corrections as calculated by W17 to account for the over- or underdensities present in each of the equatorial fields (see e.g. Loveday et al. 2015). These multiplicative corrections were derived as a function of redshift by determining the local density of the survey at the redshift of the galaxy in question. This is achieved by simply finding the running density as a function of redshift, and convolving this function with a kernel of width 60 Mpc. These were compared to the fiducial density, taken from a portion of the GAMA equatorial fields with stellar masses above 10¹⁰ M_⊙ and 0.07 < |$z$| < 0.19. This subset was chosen because of its high completeness level, uniform density distribution, and low uncertainty due to cosmic variance. To correct the effective volume for galaxy n, V_eff,n, we simply multiply by a factor of δ_n/⟨δ_f⟩ to obtain |$V^{\prime }_{\rm max}$|⁠.

To remove any spuriously low |$V^{\prime }_{\rm max}$| values introduced either by the density correction factor or by uncertainties in the calculated |$V^{\prime }_{\rm max}$|⁠, we employ a clipping technique. We split the galaxies into 100 stellar mass bins and remove 5 per cent of the most spurious |$V^{\prime }_{\rm max}$| values, and up-weight the remaining galaxies accordingly giving a final sample size of 15 750. For consistency with W17, the 5 per cent clipping is performed on the total sample, i.e. before the imposition of the requirement of H-ATLAS coverage, translating to the removal of ∼200 galaxies from the sample requiring H-ATLAS coverage that we use for this work. W17 perform a one-sided clipping, since higher |$V^{\prime }_{\rm max}$| values tend to be more stable than lower ones since brighter galaxies tend to have better constraints. Galaxies with high |$V^{\prime }_{\rm max}$| values also contribute less volume density and therefore tend to add less noise to their given bin than faint galaxies. Once removed, the weights of the remaining galaxies are scaled by the fraction of removed galaxies.¹¹

4.2 The shape of the DMF

The DMF, derived for the largest sample of galaxies to date, based on the optically selected GAMA sample, is shown in Fig. 5 using the two methods described in Section 4.1 to calculate volume densities. We have extended the function well below the low dust mass limit of all previous studies; indeed we extend to dust masses ∼10⁴ M_⊙ whilst dust masses above 10^4.5 M_⊙ are well constrained. We have therefore extended the observed range of the DMF by ∼2 dex in |$\rm \mathit{ M}_{d}$| compared to e.g. Dunne et al. (2011) and significantly reduced the statistical uncertainty compared to previous measurements (with ∼70 × the sample size, see Section 4.3 for more details). The offset at the low-mass end of the DMF seen between the two methods can be attributed to the differences shown in Fig. 4, the sources with the lowest dust mass tend to be those which are nearby and faint, and so most likely to be affected by small number statistics when calculating the BBD V_max values.

Figure 5.

The pV_max (purple) and BBD (blue) DMFs for the GAMA/H-ATLAS sources at |$z$| < 0.1. The data points show the observed values corrected for over- and underdensities in the GAMA fields (see W17). The solid lines are the best-fitting (minimum χ²) single Schechter functions from our SB measurements. Error bars are derived from our PB measurements. The total number of sources in each bin is shown in the top panel.

Open in new tabDownload slide

We estimate uncertainties on thevolume densities calculated here using three techniques. First, using a jackknife method in two different ways: (i) taking random subsamples of the data, and (ii) by splitting the sample by on-sky location. Secondly, we perform 1000 bootstrap resamplings on our volume densities to determine the sample errors. We refer to this as the simple bootstrap or SB method. Thirdly, we use the bootstrap technique but for each realization, we also perturb each dust mass by a Gaussian random deviate with σ set according to the 16–84 percentile uncertainty from magphys (hereafter the PB method). Unsurprisingly, Poisson noise estimates agree with all these techniques at the high-mass end (M_d > 10^7.5 M_⊙), but underestimate the uncertainty in the low dust mass bins (M_d < 10⁶ M_⊙). The random jackknife and SB error estimates agree very well (within 0.5 per cent), whereas the on-sky jackknife uncertainty is around 5 per cent higher. This is not unexpected since this method will include a component of uncertainty from cosmic variance within the survey volume. By disentangling the statistical uncertainty from the cosmic variance uncertainty, the larger uncertainty in the on-sky jackknife suggests an error due to cosmic variance of at least 7 per cent assuming that the difference is due only to cosmic variance. The cosmic variance estimator from Driver & Robotham (2010)¹² suggests an error of 16.5 per cent for the full survey volume. This is significantly higher than the effective cosmic variance that we measure, because we make corrections for the density variations within the survey volume. For the rest of this work, we use the simple bootstrap method without perturbation of the dust mass (SB) for the data points. For the uncertainties we use the bootstrap with additional perturbation using the magphys dust mass uncertainties (PB) since this takes into account both the variation within the sample and the uncertainty in the dust mass estimations themselves. As discussed in Section 4.4 the PB is likely to give biased estimates of the best-fitting parameters, but since it includes our mass uncertainties, it provides a better estimate of the uncertainties on the best-fitting parameters.

We fit a Schechter function to each of our bootstrap realizations, and use the median of the resulting values as the best-fitting value for each parameter. We use the standard deviation between the values to estimate uncertainty on the parameters. The parameters for both the pV_max and BBD fits are quoted in Table 1. Note that cosmic variance will introduce further uncertainty in our measurements. This will mostly be seen as an increased uncertainty on ϕ*, though both M* and α will also have slightly larger errors.

4.3 Comparing the dust mass function with previous work

We compare the SSF fit parameters derived here with single SF fits in the literature (Fig.6 left and Table 1). We also compare the confidence intervals for our derived parameters in Fig. 7 with previous work. For the first time we are able to directly measure the functional form at masses below 5 × 10⁵ M_⊙ and determine the low mass slope of the DMF, α. We see that there is a good overall match at the high-mass end with the Dunne et al. (2011) DMF, but at the faint end, the DMF is steeper than predicted from the Dunne et al. (2011) function suggesting larger numbers of cold or faint galaxies than expected. We note that the Dunne et al. (2011) sample is different to our DMF in two ways (i) it is a dust-selected (or rather 250-|$\mu$|m-selected) sample rather than optically selected and (ii) was drawn from the H-ATLAS science demonstration phase data, which is only 16 sq deg of the GAMA09 field at |$z$| < 0.1 and is known to be underdense compared to the other GAMA fields (Driver et al. 2011). Our DMF is also similar to the optically selected Vlahakis et al. (2005)¹³ SSF at the highest masses, though we find a higher space density of galaxies around the ‘knee’ of the function potentially due to the higher redshift limit probed in this study and improvement in statistics in this work. In general, the 2-d parameter comparisons in Fig. 7 show that the DMF in this work has intermediate values of α, M*, and ϕ* in comparison to the Clemens et al. (2013),¹⁴ Vlahakis et al. (2005), and Dunne et al. (2011) parameters but here we have tighter constraints due to the larger sample of sources. Differences could also arise because of the variation of best-fit parameters with the minimum mass limit of the fit since all the surveys have different mass ranges. We discuss the implications of changing the minimum mass limit of our fits in Appendix  C.

Figure 6.

Comparison of the (left) DMF and (right) dust mass densities Ω_d from this work with those from the literature. We compare with (i) the blind, local |$z$| < 0.01 galaxy sample from Clark et al. (2015), (ii) the all-sky local star-forming galaxies from the bright Planck catalogue from Clemens et al. (2013), (iii) the ground-based submm measurements of local optical galaxies from Vlahakis et al. (2005), and (iv) the 222 galaxies out to |$z$| < 0.1 from the H-ATLAS survey (Dunne et al. 2011). SF fit parameters are listed in Table1. The dust density parameter (Ω_d) measurements are scaled to the same cosmology, with diamonds representing dust-selected measurements and circles representing optically selected samples. Our work is shown as pV_max and BBD for the single SF fit – SSF to each. The solid error bars on Ω_d indicate the published uncertainty derived from the error in the fit whilst the transparent error bars indicate the total uncertainty derived by combining the published uncertainty and the cosmic variance uncertainty estimate for that sample (where known). We note that the solid error bars indicating the uncertainty from our bootstrap analysis lie within the point itself for both our BBD and pV_max values. The Dunne et al. (2011) DMF includes the correction factor of 1.42 for the density of the GAMA09 field (Driver et al. 2011) whilst our data points have been weighted by density correction factors from W17. The shaded region emphasizes the range of Ω_d derived from our observed SSF fits to the DMFs with width showing the error from cosmic variance.

Open in new tabDownload slide

Figure 7.

The confidence intervals for the pV_max single SF DMF fit parameters derived in this work (blue ellipses) showing the correlation between the fit parameters (insets) and comparison with previous values (note that ϕ* is in units of |$\rm Mpc^{-3}\ dex^{-1}$|⁠). Error bars on our fit parameters are taken from the Δχ² = 1 for each parameter (these are consistent with errors derived from the bootstrap process described in Section 4). The contours are from the 1σ, 2σ, 3σ values of Δχ² for the parameter slice centred on the best fit for the non-plotted third parameter. Green denotes Vlahakis et al. (2005), orange represents Dunne et al. (2011), and grey shows Clemens et al. (2013). We note that the error bars on the Vlahakis et al. (2005) values were derived using Poisson statistics, and so may be an underestimate of the error in the measurements.

Open in new tabDownload slide

The integrated dust mass density parameter Ω_d at |$z$| ≤ 0.1 is derived by using the incomplete gamma function to integrate down to M_d = 10⁴ M_⊙ (our lower limit on measurement of the form of the DMF). This gives (1.11 ± 0.02) × 10⁻⁶ for both the pV_max and BBD methods. For comparison, our Ω_d values calculated without imposing this limit are (1.11 ± 0.02) × 10⁻⁶ and (1.06 ± 0.01) × 10⁻⁶ for the pV_max and BBD methods, respectively, so the difference is very small. Previous measurements of Ω_d are shown in Fig. 6 (right) (all scaled to samecosmology and κ), we also recalculate the literature values using the SSF fit parameters from Table 1, this ensures that they are integrated down to our mass limit. Our measurement is consistent with Dunne et al. (2011), Vlahakis et al. (2005), Clemens et al. (2013) and with the lower range of Driver et al. (2007) but smaller than the Clark et al. (2015) values. However, the latter measurement is subject to a large uncertainty due to cosmic variance (46.6 per cent, Driver et al. 2007) in comparison to the  7–17 per cent for this work.¹⁵ Further discussion on the evolution of the dust properties over cosmic time is provided inDriver et al. (2018).

4.4 Eddington bias in the dust mass function

Here we check whether our DMF is biased due to the dust mass errors from magphys. Since the scatter due to the mass error could moves galaxies into neighbouring bins in either direction, and as the volume density is not uniform, this could have the effect of introducing an Eddington bias (Eddington 1913) into the DMF. Loveday et al. (1992) showed that this bias effectively convolves the underlying DMF with a Gaussian with width equal to the size of the scatter in the variable of interest (here dust mass) to give the observed DMF. This is valid assuming that the parameter uncertainties, and hence resulting errors, have a Gaussian distribution. Here we test whether we can correct for the Eddington bias in the DMF by deconvolving our observed DMF and attempt to extract the underlying ‘true’ DMF. We expect that any bias in the overall cosmic dust density will be small since galaxies with at least one measurement over  3σ in one of the Herschel SPIRE bands contribute around four times as much to the dust density of the Universe than those without a  3σ measurement in the FIR regime.

We fit an SSF convolved with a Gaussian, where we estimate the width of the Gaussian using two methods. First, we derive the width of the convolved function by calculating the mean dust mass error from magphys as a function of mass (the varying error method). Secondly, we take the mean value of the error in dust mass around the knee of the single SF where the convolution will have the strongest effect (where the mean error is |$0.11\,\rm dex$|⁠, the constant error method). Both produce very similar deconvolved SF fit parameters that are in agreement with the traditional SF method within a few per cent. The deconvolved fit parameters derived with constant error are listed in Table1; this produces a dust mass density of (1.08 ± 0.02) × 10⁻⁶ for both the pV_max and BBD DMFs. We find that the traditional SSF is a better fit (Δχ² ∼ 0.75) than the deconvolved constant error function, and the varying error method produces a comparable goodness of fit to the traditional SSF without deconvolution. The reason that the best fit is insensitive to the mass errors is that the mass errors are a strong function of mass: for low mass galaxies, the errors are large ( ∼0.5 dex); while for higher masses (∼M*), the errors are small (<0.1 dex). At low masses the DMF is a power law, the slope of which is unchanged when convolved by a Gaussian. At higher masses, near the exponential cut-off, the errors are small, and so the effect on the knee is negligible. We therefore conclude that there is no strong argument for choosing to use the deconvolved SSF fits instead of the original SSF fits, therefore we include the results here for completeness but continue using the original SSF fits throughout the paper.

In principle, the difference between the DMFs simple bootstrap method (the SB) and the bootstrap method where we perturbed the data by the underlying uncertainties in dust mass (the PB) provides us another method to test whether the DMF is biased and could provide a way to correct for this. Fig. 8 compares the data and resulting SF fits for the observed pV_max DMF (SB DMF) and the perturbed (including the uncertainties in the dust mass from magphys) pV_max DMF (PB DMF). We see that the two DMFs are very similar with fit properties differing by only a few per cent. The largest differences in the DMF are seen at the noisier low dust mass end suggesting the biases are indeed small, we believe this is because the uncertainties in the DMF around the knee aresmall.

Figure 8.

Top: The pV_max DMF (purple) from the SB measurements compared to the DMF derived using the bootstrap perturbed (PB) by the uncertainties in the dust mass estimates from magphys (the PB DMF, green). The data points show the volume densities in each mass bin and the solid lines are the best-fitting (χ²) SSF, to the data. Bottom: Comparison of the SSF with the DSF including the major and minor components. The data points show the volume densities in each mass bin, the major and minor components are shown in grey and purple, respectively, the overall DSF is shown in magenta. Error bars are derived from a bootstrap analysis and the data points have been corrected for over and under densities in the GAMA fields (see W17). The total number of sources in each bin is shown in the top panel.

Open in new tabDownload slide

We also perform another test to quantify the bias introduced to the DMF by the inclusion of sources with poor FIR constraints. We use the distribution of temperatures of sources with high total FIR signal to noise to define a new temperature prior. Then for each bootstrap sample we draw new temperatures from this prior and adjust the dust masses accordingly. In this way we perform another kind of perturbed bootstrap in which each realization has a temperature distribution that matches the high signal-to-noise galaxies. We find that the bias introduced to the DMF in this way is very small, and so we believe our DMF is robust. This is discussed in more detail in Appendix  A.

4.5 A double Schechter fit to the DMF

It is tempting to link the two SF components to early- and late-type galaxies, but the parameters of the minor component of the DSF do not match those of the early types (see Section 6). This suggests that the two components of the DSF do not represent physically distinct populations, and so does not provide a better representation of the data.

5 THEORETICAL PREDICTIONS FROM GALAXY FORMATION MODELS

Next we compare the SSF fit to the DMF from Section 4.4 with theoretical predictions for |$z$| = 0 from the dust models of Popping, Somerville & Galametz (2017) and McKinnon et al. (2017). Popping et al. (2017) derive DMFs from semi-analytic models (SAMs) of galaxy formation based on cosmological merger trees from Somerville, Popping & Trager (2015) and Popping, Somerville & Trager (2014) and include prescriptions for metal and dust formation based on chemical evolution models. They predict DMFs at different redshifts using dust models with dust sources from stars in stellar winds and supernovae (SNe), grain growth in the ISM, and dust destruction by SN shocks and hot halo gas (see also Dwek 1998; Morgan & Edmunds 2003; Michałowski, Watson & Hjorth 2010; Dunne et al. 2011; Asano et al. 2013; Rowlands et al. 2014; Feldmann 2015; De Vis et al. 2017b). Note that for consistency, we have scaled the Popping DMFs down by a factor of 2.39 in dust mass since their |$z$| = 0 models were calibrated on dust masses for local galaxy samples from the Herschel Reference Survey (Boselli et al. 2010; Ciesla et al. 2012; Smith et al. 2012b) and KINGFISH (Skibba et al. 2011) where Draine (2003) dust absorption coefficients are assumed. After this scaling, their DMF (based on their SAMs) is consistent with a Schechter function with |$M^* \sim 10^{7.9}\, \rm \mathrm{M}_{\odot }$|⁠. In Fig. 9 (top) we compare three of their |$z$| = 0 DMF models as defined in Table 2: the so-called fiducial, high-cond, and no-acc models. Their fiducial model assumes 20 per cent of metals from stellar winds of low-intermediate mass stars (LIMS) and SNe are condensed into dust grains, with interstellar grain growth also allowed. The high-cond assumes that almost all metals available to form dust that are ejected by stars and SNe are condensed into dust grains, with additional interstellar grain growth. The no-acc model assumes 100 per cent of all metals available to form dust that are ejected by stars and SNe are condensed into dust grains, with no grain growth in the ISM.

Figure 9.

Top: A comparison with the predicted |$z$| = 0 DMFs from Popping et al. (2017) (P17) and McKinnon et al. (2017) (McK17) with the SSF fits derived from the BBD and pV_max methods, see also Table 2. We include three models from P17: the fiducial, no-acc, and high-cond models which consist of varying dust condensation efficiencies in stellar winds, supernovae and grain growth in the ISM, respectively. The McK17 histogram is their L25n512 simulation at |$z$| = 0 (their fig. 2). Bottom: Comparing the |$z$| = 0 stellar mass functions for the GAMA sources (W17, in blue) with that derived using the SAMs of Popping et al. (2017) (in black). W17 is the SMF of the same optical sample from which our DMF is derived. The vertical line shows the boundary at which W17 fit their data with a Schechter function.

Open in new tabDownload slide

The fiducial and high-cond models overpredict the number density of galaxies in the high dust mass regime, >10^7.5 M_⊙. The no-acc model is the closest model to the observed high-mass regime of the DMF, though underestimates the volume density around M* compared to our DMF (dotted lines in Fig. 9). Both the no-acc and high-cond models are better matches at low masses (<10⁷ M_⊙), while the fiducial model underpredicts the volume density in this regime. This likely suggests that LIMS and SNe have to be more efficient than the fiducial model at producing dust in low dust mass systems, i.e. the dust condensation efficiencies in both stellar sources need to be larger than 0.3, or that the dust destruction and dust grain growth time-scales in the fiducial model need to be increased and decreased, respectively. At high masses, the fiducial and high-cond models appear to be forming too much dust. This implies that dust production and destruction are not realistically balanced in these models. This is likely due to the model introducing too much interstellar gas and metals, which allow for very high levels of grain growth in the ISM.

We note that the no-acc P17 model (without grain growth in the ISM) is likely not a valid model as it assumes 100 per cent efficiency for the available metals condensing into dust in LIMS and SNe which is unphysically high, see e.g. Morgan & Edmunds (2003) and Rowlands et al. (2014). Hereafter we no longer discuss this model even though by eye it appears to be an adequate fit to the observed DMF at masses below |$\rm 10^{7}\,\mathrm{M}_{\odot }$|⁠.

To investigate the discrepancy between the observed DMF in this work and the predicted SAM DMF from Popping et al. (2017), we first check that the stellar mass function from the SAMs is consistent with the observed GSMF for the GAMA sample in W17 (Fig. 9 bottom). The SMFs at the high-mass end are in agreement though the model SMF has a slight overdensity of galaxies in the range 10⁸ < M_s (M_⊙) < 10^9.4, where M_s is stellar mass. If this overdensity of sources were responsible for the discrepancy between the predicted and observed DMFs in the high M_d regime, those intermediate stellar mass sources would have to have dust-to-stellar mass ratios of ∼0.5 which is again unphysical. We can see this is not the case when comparing the dust-to-stellar mass ratios of the Popping et al. (2017) fiducial |$z$| = 0 model in Fig. 10 (as mentioned earlier, this is based on Herschel observations of local samples of galaxies).

Figure 10.

The dust to stellar mass ratio for galaxies in the local Universe. The data from this work are shown in the underlying grey points with mean dust masses (±standard error) in each stellar mass bin (black). We include a compilation of Herschel results for local galaxies including the stellar-mass selected HRS (Boselli et al. 2010), the dust-selected sample of Clark et al. (2015), the H i-selected sample from De Vis et al. (2017b) and the dwarf galaxy survey from Rémy-Ruyer et al. (2013). Overlaid are the local Universe relationships (⁠|$z$| < 0.12) based on stacking on Herschel maps for 80 000 galaxies from Bourne et al. (2012) in three different g − r colour bins (their fig. 16). All of these samples have been scaled to the same the κ value with parameters derived using magphys (see De Vis et al. 2017b) or modified blackbody fitting but scaled to the same κ (DGS and stacked samples). The median dust and stellar masses from Popping et al. (2017) are shown by the grey line with 16 and 84 percentile error bars (scaled by a factor of 1/2.39 in dust mass).

Open in new tabDownload slide

In Fig. 10 we plot the dust and stellar masses from the compilation of local galaxy samples collated in De Vis et al. (2017a,b)¹⁶ and compare with P17 and our sample of ∼15 000 sources. Here we can clearly see the cause for the discrepancy between the observed DMF from this work and the model: the model overpredicts the amount of dust in high stellar mass sources, well above any dust-to-stellar ratios observed locally. Although the observations show a flattening of dust mass at the highest M_s regime (where early type galaxies are dominating), this is not the case in the SAM. In general the SAM prediction assumes a constant dust-to-stars ratio of ∼0.001 across all mass ranges. The observations however suggest that there is a roughly linear relationship until M_s > 10¹⁰ M_⊙, after which the slope flattens, with M_d/M_s ∼ 0.001.

Fig. 10 also suggests that M_d/M_s increases to ∼0.025 in low-stellar-mass galaxies (in agreement with Santini et al. 2014; Clark et al. 2015; De Vis et al. 2017a). This is further supported by the stacking analysis carried out in Bourne et al. (2012) whose dust-to-stellar mass trends in different bins of optical colour are added to Fig. 10. These were derived by stacking on ∼80 000 galaxies in the Herschel maps, revealing that low stellar mass galaxies had higher dust-to-stellar mass ratios, consistent with these sources having the highest specific star formation rates. Our binned data (black points) are in agreement with local galaxy surveys and the Bourne et al. (2012) trends: we see that the slope of dust-to-stellar mass flattens at high masses, and that there exists a population of dusty low-stellar-mass sources that the SAM does not predict.

Alternative predictions for a local DMF are provided by McKinnon, Torrey & Vogelsberger (2016, hereafter McK16) and McKinnon et al. (2017, hereafter McK17). In these models, dust is tracked in a hydrodynamical cosmological simulation with limited volume. The McK16 dust model is similar to the P17 high-cond model (including interstellar grain growth and dust contributed by both low mass stellar winds and SNe) but has no thermal sputtering component. The updated model from McK17 reduces the efficiency of interstellar grain growth and includes thermal sputtering (see Table 2). The DMF from McK17 (their L25n512 simulation at |$z$| = 0) is shown in Fig. 9 (top). Their values have been scaled to the same cosmology as used here (they use the same κ and Chabrier IMF as this work). We can see that McK17 predicts fewer massively dusty galaxies than P17 and our observed DMFs. Although their DMF fails to produce enough galaxies in the highest mass bins in Fig. 9, the simulated DMF becomes more strongly affected by Poissonion statistics in this regime due to the small volume of thesimulation.

Possible explanations for the difference between the predicted (P17, Mck17) and observed DMFs at large dust masses are (i) the efficiency of thermal sputtering due to hot gas in the halo has been under or overestimated in these highest stellar mass sources; (ii) the fiducial and high-cond dust models of P17 allow too much interstellar grain growth in highest stellar mass galaxies due to the assumed time-scales or efficiencies of grain growth being too high; (iii) the predicted highest stellar mass galaxies have too little (McK17) or too much (P17) gas reservoir potentially due to feedback prescriptions being too strong/not strong enough, respectively. If the gas reservoir is too high, interstellar metals can continue to accrete onto dust grains and increase the dust mass. Conversely if it is too low, then the contribution to the dust mass via grain growth will be reduced. We will address each of these possibilities inturn.

• We can test if the amount of dust destruction by thermal sputtering in hot (X-ray emitting) gas could explain the differences in the predicted and observed DMF at the high mass end as McK16 and McK17 already compared the results using dust models without and with thermal sputtering respectively. They find that including thermal sputtering only makes small changes to the shape of DMF since this affects dust in the halo rather the ISM, this is therefore not likely to be responsible for the discrepancy.

• Comparing the dust models in P17, McK16, and McK17 allow us to test the effect of changing the grain growth parameters. The time-scale for grain growth is shortest in P17 and McK16 and both those models produce too much dust in the high dust mass regime of the DMF. McK17 has a longer grain growth time-scale (⁠|$t_{\rm acc,0} =40\,\rm Myr$|⁠, Table 2) than both P17 and McK16 and this change indeed reduces the volume density of the highest dust mass sources. McK17 also compares the DMFs from the same simulation methods with different dust models and they find that a significantly reduced DMF at the high mass end can be attributed to the longer grain growth time-scales.

• Earlier we showed that the galaxies that are responsible for the highest dust mass bins in the P17 DMF have too much dust for their stellar mass (Fig. 10). To test whether they have too much dust due to the gas reservoir of the SAM massive galaxies being too high (hence leading to more interstellar grain growth) we refer to the predicted and observed gas mass function comparison in Popping et al. (2014). There they showed that these are not as discrepant as we see here with the modelled and observed DMFs and therefore are likely not responsible for the discrepancy in the DMF.

We therefore conclude that it is likely that the interstellar grain growth in these massive galaxies is simply too efficient/fast in the P17 and McK16 dust models. In this scenario, the few largest stellar mass galaxies are allowed to form too much dust in the ISM at a rate that is not observed in real galaxies. However, the growth time-scale may also be too slow in the McK17 model. All of the P17 high-cond, McK16 and McK17 dust models assume very high dust condensation efficiencies in AGB stars and Type Ia and II SNe. We propose therefore that the most realistic dust model must lie somewhere in between these and the fiducial P17 model, with stardust condensation efficiencies larger than 0.3 but lower than 0.8 and a similar dust grain growth time-scale as assumed in P17.

Neither the P17 fiducial, nor the McK16 and McK17 dust models provide reasonable matches to the low dust mass regime (∼10^7.5 M_⊙) of the DMF. McK16 and McK17 overpredicts the dust masses in the low mass regime and P17 fiducial model underpredicts the DMF suggesting again that stardust condensation efficiencies may be intermediate between the three models with a grain growth time-scale similar to P17. Only the P17 high-cond model provides an adequate match to this regime. Based on the observed DMFs, we will explore ways to improve the theoretical models in future work.

6 THE DMF BY MORPHOLOGICAL TYPE

We can test the standard prejudice that spiral galaxies are full of dust, and ellipticals have very little dust by using the DMF to quantify the difference in dust content of ETGs and LTGs. We create ETG and LTG subsets for our sample of galaxies based on classifications carried out by GAMA in Driver et al. (2012, hereafter D12) and Moffett et al. (2016a, hereafter M16a). For both studies visual classifications were based on three colour images built from H (VIKING), i, and g (SDSS) bands. Classifications were based on three pairs of classifiers in which there was an initial classifier and a classification reviewer. D12 classified the entire sample out to |$z$| ≤ 0.1 and split the sample only into ‘Elliptical’ and ‘Not Elliptical’ galaxies, which we hereafter refer to as ETGs and LTGs (later type galaxies). The classifications from M16b were carried out on the same sample as D12, but limited to |$z$| ≤ 0.06. In M16a they attempted to produce an updated set of morphological classifications using classification trees with a finer binning system than D12. However, for consistency with the D12 classifications (and because here we do not want to split the DMF into finer morphological classes), we include the M16a Ellipticals in the ETG class, and we group all remaining galaxy types apart from little blue spheroids (LBSs) into the LTG category (following M16a we keep the LBS sources separate). We note that if LBSs are included in the LTG category there is very little change in either M* or ϕ*; however, the low-mass slope is steepened by ∼6 per cent, overall the difference in the dust mass density made by including LBSs in the LTG category is only around 2 per cent.

We next use these morphologies in order to investigate the shape and dust mass density of ETGs and LTGs. We choose to limit our redshift range to |$z$| ≤ 0.06 for two reasons (i) the finer, updated classifications of D12 provided in M16a is limited to this range and beyond this range visual classifications become more uncertain (ii) with increasing redshift, the sample will suffer more from incompleteness at lower masses. This is demonstrated in Fig. 12 where we compare the dust and stellar masses of the ETGs and LTGs from D12 at |$z$| ≤ 0.06 and 0.06 < |$z$| ≤ 0.1. We see a dearth of galaxies below M_d ∼ 10^5.5 and stellar masses below M_s ∼ 10^7.5 in the higher redshift bin.

In our final |$z$| < 0.06 sample, a total of 5736 sources were classified by D12, 588 of which were ETGs, 4837 as LTGs, and 474 as LBS. In the same redshift range, M16a classified 5765 galaxies with 639 ETGs, 4599 LTGs, and 690 LBSs. There are 773 disagreements between the two sets of classifications (13 per cent of the overall sample). The resulting DMFs from the two classification methods are displayed in Fig. 11 ( using SSF fits), the red and blue points are ETG and LTG respectively, and translucent and solid represent the DMFs for the early and late populations as defined by D12 and M16a, respectively. The DMFs for D12 and M16a agree well out to |$z$| ≤ 0.06 showing that although the individual classifications are not an exact match, the shapes of the DMFs of the ETG and LTG populations appear to be consistent between the two different classification methods. The SF fit parameters and dust mass densities for ETGs and LTGs are listed in Table3. (From now on we choose to discuss only the SF fit parameters arising from the M16a classifications as these are the most recent). Unsurprisingly we see there is an order of magnitude more dust mass contained within LTGs than ETGs at |$z$| ≤ 0.06, with dust mass densities of Ω_d = (0.88 ± 0.03) × 10⁻⁶ and Ω_d = (0.060 ± 0.005) × 10⁻⁶, respectively. LTGs are mostly responsible for the dust content of the Universe. However, Fig. 11 demonstrates that the ETG DMF is not well described by a Schechter function, indeed there is a significant downturn in the volume density of ETGs below dust masses of 10⁶ M_⊙. We believe this is a real effect since we have no reason to believe that incompleteness may be biasing our measurements at this redshift range. The downturn is also in line with the shape of the GSMF for this subset as measured by M16a. Since the DMF for the ETGs clearly does not match a Schechter function, the dust mass density in Table 3 may be overestimated. We therefore also calculate the dust mass density for these galaxies by summing the contribution from each galaxy and derive a density of Ω_d = (0.060 ± 0.005) × 10⁻⁶. This is consistent with the integrated Schechter function since although that fit overpredicts the total dust mass density in low dust mass sources, it also underpredicts the dust mass density at the high-mass end. Comparing our ETG and LTG SF fits with the double component fit to the entire sample from Section 4.5, we find that the major component of the double fit matches the high mass end but slightly overshoots the volume densities derived for the LTGs at intermediate masses 10⁵ < M_d(M_⊙) < 10⁷ whereas the second component has a peak volume density at higher dust masses than the ETGs. However, this may be due to misclassification of galaxy types or simply due to the degeneracy in what the fitting routine assigns to each component at the faint end.

Figure 11.

The pV_max DMFs for the GAMA/H-ATLAS sources at 0.002 < |$z$| < 0.06. Here the opaque lines show the sample for the fitted range split into ETGs and LTGs by Moffett et al. (2016a) – M16a, and the translucent lines show the sample for the fitted range as split into ETG and LTG by Driver et al. (2012) – D12. Red denotes ETGs and blue the LTGs. The data points show the observed values and the solid lines are the best-fitting (χ²) single SF fits to the data for their respective fitted regions, beyond this we show extrapolations down to 10⁴ M_⊙ as dashed lines. Error bars are derived from a bootstrap analysis and the data points have been corrected for over and under densities in the GAMA fields (see W17).

Open in new tabDownload slide

7 COMPARISON WITH THE GALAXY STELLAR MASS FUNCTION

Scaling relations between dust and stellar mass can reveal the relation between internal galaxy properties and the dust content and whether there is a simple prescription that can tell us how much dust exists in galaxies given a unit of stellar mass (e.g. Driver et al. 2018). Cortese et al. (2012) and Smith et al. (2012b) investigated M_d−M_s scaling relations in local galaxies using the HRS, finding that larger stellar mass galaxies have lower dust-to-stellar mass ratios (see also Santini et al. 2014). This was further confirmed in the larger statistical study of Bourne et al. (2012) from H-ATLAS (Fig. 10). As we saw in Section 5, the Popping et al. (2017) SAMs produce a trend in M_d−M_s which does not agree with the scatter seen in the observations of local galaxies (due to colour, morphological type, environment, etc.) and their models produce too much dust in the highest stellar mass discs. Further modelling of dust and stellar mass scaling relations carried out in Bekki (2013) using chemodynamical simulations reproduced the M_d−M_s trend roughly (with massive disc galaxies more likely to have smaller dust-to-stellar mass ratios), but could not reproduce the M_d−M_s at stellar masses >10¹⁰ M_⊙. In Fig. 12 we bin galaxies in stellar mass, and plot the mean dust mass as a function of mean stellar mass in each bin, with error bars calculated from the standard error on the average. The low-mass end of the LTG dust and stellar mass comparison is actually fairly well-represented by a linear relationship, but diverges at higher masses. The change in slope at high M_s is largely caused by LTGs with M_d similar to ETGs at similar M_s. It is difficult to disentangle whether this effect is entirely physical due to dust-poor LTGs or due to the misclassification of ETGs as later type galaxies. It is clear from this plot that the ETGs do not follow a linear trend.

Figure 12.

The mean dust to stellar mass ratio for galaxies in our high- and low-redshift samples (Top: 0.002 ≤ |$z$| ≤ 0.06, Bottom: 0.06 ≤ |$z$| ≤ 0.1) in bins of galaxy stellar mass for ETGs (red), LTGs (blue) in comparison to the ETG and LTG subsamples. We also show in blue a straight-line fit to the LTG data to illustrate the approximate linear scaling of the LTG dust-to-stellar mass ratio at low masses.

Open in new tabDownload slide

With our dataset we can test whether there is a simple scaling relation between dust and stellar mass by simply comparing the dust and stellar mass functions. Since the GSMF in W17 is fitted by a coupled DSF with shared M*, it is not equivalent to our DSF, we instead compare the SSF fit with the GSMFs from M16a who present the GSMFs for ellipticals and for later morphological types including Sab-Scd, SBab-SBcd, Sd-Irr, and S0-Sa. We also compare with Moffett et al. (2016b, hereafter M16b) who have further decomposed the sample into bulges, spheroids, and disc components. This decomposition of galaxies into bulge and disc was performed by fitting a double-Sérsic profile to those galaxies which were morphologically classified as double-component galaxies to obtain bulge-to-total luminosity ratios. From these ratios, the g − i colour and i-band absolute magnitudes for both bulge and disc were derived and used, along with the stellar mass relation of Taylor et al. (2011), to calculate bulge and disc component stellar masses. We have already seen that the DMF is dominated by the LTGs in the previous section, we now test whether most of the dust will be associated with the disc component of galaxies and whether we can scale from stellar mass to dust mass for different galaxy subsamples. In Table 4 we compare the ratio of the knee of the SF fit parameters (⁠|$M^{*}_{\rm d}/M^{*}_{\rm s}$|⁠) for the GSMF and DMF and the integrated mass densities (ρ_d/ρ_s) between different galaxy populations including LTG, ETGs, and discs.

First we compare the GSMFs to their equivalent DMFs for ETGs and LTGs. We produce a composite Schechter function from the later type GSMFs from M16a containing the same sample of galaxies as our LTG DMF. We show a version of the LTG GSMF composite Schechter function as scaled by the ratio ρ_d/ρ_s in Fig. 13. The scaled LTG GSMF fits the high-mass end of the LTG DMF; however, it diverges from the data points around 10⁷ M_⊙ where a more pronounced shoulder is seen in the GSMF than we observe in the DMF. Otherwise the composite LTG GSMF is in good agreement with our data, and therefore an estimate of the LTG DMF could be made by scaling the LTG GSMF by a factor of (8.07 ± 0.35) × 10⁻⁴.

Figure 13.

The pV_max DMFs for the GAMA/H-ATLAS sources at 0.002 < |$z$| < 0.06. The SF fits for the ETGs and LTGs are shown by the solid red and blue lines for the fitted range, beyond this we show extrapolations down to 10⁴ M_⊙ as dashed lines. We also compare the GAMA GSMFs from M16a and M16b scaled by ρ_d/ρ_s including the M16a Elliptical and late types GSMFs scaled by ρ_d,ETG/ρ_s,Elliptical and ρ_d,LTG/ρ_{s,(Sab-Scd + Irr + S0-Sa)} (purple and green, respectively) and the M16b disc GSMF scaled by ρ_d,LTG/ρ_s,disc (yellow).

Open in new tabDownload slide

M16a includes an Elliptical GSMF (equivalent to what we have defined here as ETG), and we scale this function by the ρ_d/ρ_s in order to compare with the ETG DMF. As seen in Fig. 13, the scaled Elliptical GSMF is not a good match to the ETG DMF data or Schechter function. Compared to the data, the scaled GSMF is too high for M_d > 10⁷ and M_d < 10^5.5, and too low for 10^5.5 < M_d < 10^6.5. Indeed, the Elliptical GSMF is also not well-fitted by a Schechter function, and displays the same downturn that we see in the DMF. Whilst the low-mass slope derived by M16a does show a drop-off, it is not as severe as actually observed in either the dust mass or stellar mass function data. Also, we note that the dust mass as a function of stellar mass for ETGs as seen in Fig. 12 is not consistent with a simple linear relationship. We therefore caution against using a scaling law between stars and dust that relies upon a Schechter function fit to ETGs. We also note that the ratio |$M^{*}_{\rm d}/M^{*}_{\rm s}$| for this sample is |$\left(0.94^{+0.25}_{-0.24}\right) \times 10^{-4}$|⁠, which is of order 17–25times higher than the average dust-to-stellar mass ratios of Ellipticals compared to recent Herschel studies of ETGs in the local volume (⁠|$D\lt 40\rm \,Mpc$|⁠, Smith et al. 2012c; Amblard et al. 2014). This may suggest some contamination in our E category from S0s, though we note that the Herschel studies were potentially biased to older redder sources (Rowlands et al. 2012; Clark et al. 2015) and high-density environments (Smith et al. 2012c; Agius et al. 2013).

Next we attempt to explore the relationship between stellar mass of the disc component and the dust mass of a galaxy. M16b found that the disc GSMF¹⁷ is well fitted by a Schechter function with α_disc = −1.20 ± 0.02, which is consistent with our α values both for the LTG DMF and the total DMF for the |$z$| < 0.06 sample. Based on the compatibility of the α values, there may be a simple scaling from the disc GMSF to either the total or LTG population DMFs. The scaled function is a good but imperfect fit, and we see a moderate overshoot of the high mass data points by this scaled function, reflecting the fact that the dust-to-stellar mass ratio is not constant. We therefore conclude most of the dust in galaxies is associated with their disc components, and that it is possible to obtain a reasonable representation of the DMF by scaling the disc GSMF by the ratio ρ_d/ρ_s = (10.21 ± 0.45) × 10⁻⁴. We can see that the ratios |$M^{*}_{\rm d}/M^{*}_{\rm s}$| and ρ_d/ρ_s for the disc GSMF and both the total and LTG DMFs shown in Table4 are discrepant by more than 1 σ. This provides further evidence that the scaling from the disc GSMF to the DMF of either population cannot be exactly linear. We also infer from this that the dust-to-stellar mass ratio is higher for lower mass discs. We refer the reader to De Vis et al. (2017b) whose work hints that the observed dust-to-stellar mass properties of local galaxies may require the contribution of dust sources from stars and interstellar grain growth to be different for low- and high-mass galaxies.

Given that the observed dust-to-stellar masses of galaxies are not linear across the whole stellar mass range (Figs 10 and 12), it is somewhat surprising that we can simply scale the GSMFs of LTGs and discs and obtain DMFs close to the observed. Although the binned dust masses in Fig. 12 at stellar masses greater than 10^9.5 M_⊙ depart from a linear scaling relation, on average we can assume the slopes are close to linear (especially around the knee of the SF) and therefore this simple scaling appears to work. Surprisingly, this simple scaling from stellar mass to dust mass for LTGs at z= 0 (Fig. 13) may suggest that all the different dust processes (dust condensation in stellar atmospheres and SNe, grain-growth, dust destruction) are correlated to the growth of stellar mass in galaxies. The dispersion in this scale (e.g. Figs 10 and 12) could therefore place limits on the way dust is formed and how it evolves. We will return to this in future work.

We note that the uncertainties on the data points in the GSMFs quoted in M16a and M16b are based on an on-sky jackknife analysis. As we described in Section 4.2, this estimate will include a cosmic variance uncertainty component. Since we cannot disentangle the inherent cosmic variance from their values we choose instead to use the same percentage uncertainty in the integrated stellar mass density as our percentage uncertainty in the integrated dust mass density. Since the errors in dust mass are larger than for stellar mass for our dataset, the uncertainty in the integrated stellar mass density should not be larger than for our integrated dust mass density for the same sample; as such this estimate of the error is a conservative one.

8 CONCLUSIONS

We measure the DMF for galaxies at 0.002 ≤ |$z$| ≤ 0.1 using a modified V_max method for 15 750 sources. This represents the largest study of its kind both in terms of numbers of galaxies and the volume probed at this redshift range. Dust masses are derived using the magphys SED fitting tool given photometric fluxes from lambdar. Despite the sources that have measurements below  3σ in one or more Herschel SPIRE bands, we show that the magphys derived errors of 0.4, 0.18, 0.14, and 0.1 dex for galaxies with flux  >3σ in zero, one, two, and three SPIRE bands do reasonably represent the uncertainty in dust mass.

We use a single Schechter function fit to our data to compare with previous observations and we extend the DMF down to lower dust masses than probed before, constraining the faint end slope below 10⁵ M_⊙. Our main findings are the following:

• We compare the DMF derived using the traditional pV_max method and the BBD method which allows us to incorporate selection effects, and find both ultimately produce similar results. The best-fitting single Schechter function for the |$z$| ≤ 0.1 DMF has α = −1.22 ± 0.01, |$M^* = (4.65 \pm 0.18) \times 10^7\,h^2_{70}\, \mathrm{M}_{\odot }$|⁠, |$\phi ^* = (6.26 \pm 0.28) \times 10^{-3}\,h^3_{70}\,\rm Mpc^{-3}\,\rm dex^{-1}$|⁠, and Ω_d = (1.11 ± 0.02) × 10⁻⁶. There is an additional uncertainty from cosmic variance of  7–17 per cent.

• We find that a double Schechter function is formally a better fit to the observed DMF, though including cosmic variance would make the improvement in the fit not significant. Also given the variation of errors as a function of mass and the lack of correspondence to physical subsets of galaxies, we do not have confidence that the double Schechter function represents a better description of the DMF.

• We find that there is a discrepancy between the observed and predicted DMFs derived from the semi-analytic models of Popping et al. (2017). This is largest at the high-stellar mass end, with the model predictions of M* higher by 0.5 dex compared to our observed, Schechter functions. The likely cause for the discrepancy is that the Popping model uses a relationship between dust and stellar mass which is inconsistent with properties observed in local galaxies samples such as the Herschel-ATLAS, the HRS and the DGS, and also with our sample of GAMA sources: the models produce high stellar mass galaxies with dust masses far higher than is observed. This discrepancy is alleviated somewhat when we compare with the predicted DMF from cosmological hydrodynamical simulations (McKinnon et al. 2017) who use longer grain growth time-scales. This reduces the amount of dust formed in high mass galaxies; however, McKinnon et al. (2017) underpredict the number of high dust mass galaxies compared to our observations, although the limited volume of their simulation does not allow a proper comparison. Both sets of theoretical predictions also fail to match the observed volume density of low dust mass galaxies. Our dataset thus provides a useful benchmark formodels.

• Splitting our sample into early- and late-type on the basis of morphology and colour (to a redshift limit of |$z$| ≤ 0.06), results in DMFs with very different shapes. The late-type DMF is well represented by a Schechter function, whereas the ETG DMF is not. The LTG DMF has far higher space density at a given dust mass. We derive dust mass densities of Ω_d = (0.88 ± 0.03) × 10⁻⁶ and Ω_d = (0.060 ± 0.005) × 10⁻⁶ for late types and early types, respectively. In total there is ∼10 times more dust mass density in late-type galaxies compared to early-types at 0.002 < |$z$| ≤ 0.06.

• In comparing our DMF to the GSMFs from W17 and Moffett et al. (2016a,b), it is possible to scale from the LTG GSMF to the LTG DMF using a ratio of ρ_d/ρ_s = (8.07 ± 0.35) × 10⁻⁴. Similarly, we show that one can scale from the disc GSMF to the LTG DMF using the ratio ρ_d/ρ_s = (10.21 ± 0.45) × 10⁻⁴. We caution that using Schechter values derived from Schechter function fits to the ETG DMF and Elliptical GSMF may be inadvisable since neither are well-fitted by a Schechter function, although scaling from the Elliptical GSMF to the ETG DMF by multiplying by ρ_d/ρ_s = (1.00 ± 0.11) × 10⁻⁴ returns a reasonable representation of the ETG DMF around the knee.

Footnotes

ACKNOWLEDGEMENTS

Australian Astronomical Observatory (AAO) Australian Research Council (ARC, Science and Technology Facilities Coucil (STFC We acknowledge our anonymous referee for their helpful input. We acknowledge Edo Ibar and Michal Michałowski for their comments on earlier drafts of this paper and acknowledge Ryan McKinnon for sharing his DMF for this paper. RAB, HLG, SM, and LD acknowledge support from the European Research Council (ERC) in the form of Consolidator Grant CosmicDust (ERC-2014-CoG-647939, PI H L Gomez). LD, SM, and NB acknowledge support fromEuropean Research Council Advanced Investigator Grant COSMICISM, 321302. CJRC acknowledges support from the ERC in the form of the 7th framework programme (FP7) grant DustPedia (PI Jon Davies, proposal 606824). GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the United Kingdom Infrared Telescope (UKIRT) Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX Medium Imaging Survey (MIS), VST KiDS, Visible and Infrared Survey Telescope for Astronomy (VISTA) VIKING, WISE, Herschel-ATLAS, the Giant Metrewave Radio Telescope (GMRT), and the Australian Square Kilometre Array Pathfinder (ASKAP) providing UV to radio coverage. GAMA is funded by the , UK), the Australia), the , and the participating institutions. The GAMA website is http://www.gama-survey.org/. The Herschel-ATLAS is a project with Herschel, which is a European Space Agency (ESA) space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from National Aeronautics and Space Administration (NASA). The H-ATLAS website is http://www.H-ATLAS.org/. This research made use of Ned Wright’s cosmology calculator http://www.astro.ucla.edu/ wright/CosmoCalc.html (Wright 2006). Fig. 2 was made with the python seaborn package http://seaborn.pydata.org/index.html.

REFERENCES

APPENDIX A: FIR CONSTRAINTS ON DUST TEMPERATURE

In Section 3, we found that when a galaxy has poor FIR constraints, the cold dust temperature PDF from magphys simply reflects the prior temperature distribution. Since the FIR data provide the only constraints on the cold temperature, this is the correct result that magphys should output. However, taking the median of the PDF in these cases may lead to a systematic bias in the temperature, which would in turn lead to a bias the dust mass. Roughly one third of our sample has a combined FIR signal to noise of less than 1, so potentially, it is a significant effect. Here we attempt to quantify how this effect propagates into in our estimates of the DMF. In panel (a) of Fig. A1 we show the median cold dust temperature output by magphys as a function of total FIR signal to noise. As the total FIR signal to noise decreases there is a clear trend for galaxies to be assigned a cold dust temperature nearer to 20 K, the median of the magphys temperature prior. Panel (b) of Fig. A1 shows the distribution of cold dust temperatures for sources with a total FIR signal to noise of at least 7σ, where the constraints on the dust temperature are very good. Assuming the underlying temperature distribution is independent of observed signal to noise, we can use this as a new prior for the temperatures of sources which have poor FIR constraints.

Figure A1.

(a) The cold dust temperature from magphys as a function of total FIR signal to noise, in red all those sources below 7 σ, in black all those sources with 7 σ or greater total FIR flux. (b) In black the probability density function (PDF) of temperatures from magphys for those galaxies with a total FIR signal to noise greater than 7 using kernel density estimation (KDE), in blue the prior we use to describe this PDF, in red a KDE of one example of a random draw of temperatures from this prior. (c) In red we show one MC simulation of new cold dust temperatures for sources below 7 σ total FIR flux, in black the magphys cold dust temperatures for sources with 7 σ or greater total FIR flux. (d) One MC realization of the ratio of dust masses adjusted for their MC temperature as a function of their dust mass as assigned by magphys.

Open in new tabDownload slide

The analysis we perform is another form of perturbed bootstrap, where, instead of perturbing the masses based on the 16 and 84 percentiles output by magphys, we perturb their cold dust temperatures to match the new prior temperature distribution. As in the PB method we first resample from the underlying dataset to give each bootstrap realization. We keep the temperatures of the high signal-to-noise sources the same, but for all other galaxies we assign new temperatures using the prior from the high signal to noise sample. This means that the resulting temperature distribution for all sources matches the prior. This is shown for one realization in panel (c) of Fig. A1. Finally we adjust the galaxy dust masses in proportion to the change in black body luminosity at 250 |$\mu$|m given the original and re-assigned temperatures. We refer to this method as the perturbed temperature bootstrap (PTB). The ratio of the PTB and magphys temperatures as a function of dust mass from magphys for one realization can be seen in panel (d) of Fig. A1. The figure shows that the dust mass can increase or decrease as a result of the temperature MC, but generally the dust masses increase as a result of the new temperature.

Fig. A2 and Table A1 show the resulting DMF and best-fitting Schechter function parameters along with our simple bootstrap (SB) results. The differences are very small, showing that the poorly constrained temperatures do not introduce a significant bias in our DMF estimates, even though they make up a significant fraction of our sample. The derived α and M_s values both agree within uncertainties. The value of ϕ* is  2σ or ∼16 per cent higher for the PTB DMF. Overall the revised temperatures lead to a difference in the cosmic dust density of only ∼6 per cent. Thus we are confident that any bias from poor mass constraints for the lower signal to noise sources in our sample is at the level of a few per cent.

Figure A2.

The pV_max DMFs for the GAMA/H-ATLAS sources, first showing the SB DMF as seen in Section 4 in magenta, and secondly the DMF resulting from the PTB method shown in grey. The data points show the observed values and the solid lines are the best-fitting (χ²) single Schechter functions to the data for their respective fitted regions, beyond this we show extrapolations down to 10⁴ M_⊙ as dashed lines. Error bars are derived from a bootstrap analysis and the data points have been corrected for over and under densities in the GAMA fields (see W17).

Open in new tabDownload slide

APPENDIX B: CHANGING THE BIVARIATE DMF TO STELLAR MASS–SURFACE BRIGHTNESS

Wealso performed the BBD analysis described in Section 4.1 with stellar mass and surface brightness as the two ‘axes’. We show the raw counts of galaxies in M_s/μ_e,abs bins in Fig. B1(a), the median volume in Fig.B1(b), and the weighted counts (volume density) in Fig.B1(c). In comparison to Fig. 3, we can see that the shapes and trends of the BBDs with r-band/μ_e,abs or M_s/μ_e,abs as the second axis are unsurprisingly very similar. The main difference between the two is that the shape of the latter reflects a slightly stronger evolution of stellar mass compared to r-band magnitude with surface brightness.

Figure B1.

The BBD for our sample with r-band magnitude and surface brightness as the two ‘axes’ (W17). (a) Raw counts in stellar mass/surface brightness bins; (b) median volume in the stellar mass/surface brightness bins, (c) weighted counts, i.e. volume density in the stellar mass/surface brightness bins. Each of the panels represents the BBD resulting from the average of 1000 Monte Carlo simulations where we perturb the stellar mass and surface brightness within their associated uncertainties.

Open in new tabDownload slide

Fig. B2 reproduces Fig. 4 but now with stellar mass as the second ‘axis’ of the BBD. The pV_max and M_s/μ_{e, abs} are coloured by the number of galaxies in the BBD bin containing that galaxy. The largest deviations from the 1:1 line are seen when the galaxy lies in a more sparsely populated bin, generally these volumes are low and are therefore subject to large cosmic variance. The pV_max values are systematically higher by 1.2 per cent on average than those derived from the M_s/μ_e,abs BBD method, which gives an average offset of 6 per cent in the binned DMF values. This is because the largest differences in the binned DMF values come from the low dust mass end where BBD bins are more likely to be below the required 4 galaxies per bin. The scatter about the 1:1 line is large compared to Fig. 4 since a galaxy’s r-band magnitude has more reason to impact the sample selection than its stellar mass, and so ultimately we decide to show the r-band magnitude BBD in the main body of the paper.

Figure B2.

Left: The maximum effective volumes for our galaxies at |$z$| < 0.1 derived using the pV_max method (x-axis), and BBD method using stellar mass and surface brightness as the two selection features (y-axis). The colour of the points is determined by the number of galaxies in the BBD bin that each galaxy resides in (Fig.B1), as shown by the colour bar in the top left corner. Right: The DMF derived using a stellar mass and surface brightness BBD (blue) and an r-band magnitude and surface brightness BBD (green) for the GAMA/H-ATLAS sources at |$z$| ≤ 0.1, also shown are the pV_max values for comparison in magenta. The data points show the observed values and the solid lines are the best-fitting (χ²) single SF fits to the data. Error bars are derived from a bootstrap analysis and the data points have been corrected due to over and under densities in the GAMA fields (see W17). The residual between the two BBD DMFs is shown in the top panel, as points with error bars in purple.

Open in new tabDownload slide

Finally, Fig. B2 compares the resulting DMFs from using r-band magnitude or stellar mass as the second ‘axis’ of the BBD, where the surface brightness is used for the remaining axis. The SF fit parameters are compared in Table B1, and the residuals between the data points resulting from the both the M_s/μ_e,abs and r-band magnitude/μ_e,abs BBD are shown in Fig. B2. The SF fit parameters for both setsof BBD DMFs agree within uncertainties; however, the r-band magnitude BBD is closer to the pV_max fit parameters. The pV_max DMF points are offset from the stellar mass BBD by ∼6 per cent, and the r-band magnitude BBD by ∼3 per cent on average.