© ESO, 2016

1. Introduction

Tightly connected, dust and gas are the key ingredients of star formation, which is governed by their complex, mutual interplay as in a cyclic dance. Stars form in cold, dense molecular clouds and dust works as catalyst in transforming atomic hydrogen into molecular hydrogen (e.g., Wolfire et al. 1995). Gas is again expelled from stars during their lifetime in the form of winds and finally at the supernova (SN) stage. Dust is believed to be mainly produced in the envelopes of asymptotic giant branch (AGB) stars and at the end of the life of massive stars during the explosive SN phase (see Dunne et al. 2011, for a short summary). Supernovae shocks, on the other hand, destroy dust grains, which can form again in the interstellar medium (ISM) by an accretion process. Dust absorbs the ultraviolet (UV) emission of young stars, allowing gas to cool and condense to form new stars.

The dust and gas content of galaxies are linked to each other through metallicity, as shown in the local Universe (Leroy et al. 2011; Draine et al. 2007): the gas-to-dust mass ratio in the ISM increases as a function of metallicity. Dust has therefore been often adopted as a proxy to gas in absence of time-expensive, sub-millimeter spectroscopic observations (e.g., Magdis et al. 2012; Santini et al. 2014; Scoville et al. 2014; Genzel et al. 2015) at high redshift also, assuming that local relations hold. The energy absorbed by dust at short wavelengths is re-emitted in the infrared (IR) and sub-millimeter (sub-mm) regimes, where the thermal emission of grains dominates the spectral energy distribution (SED) of galaxies (~8−1000 μm). Dust emission is frequently used to trace the ongoing rate of star formation (SFR; e.g., Nordon et al. 2010; 2012; Elbaz et al. 2011; Kennicutt 1998a).

It is thus no wonder that the interest in the study of galaxy dust and gas properties has grown in the past two decades and that an always increasing number of studies based on Herschel data are dedicated to the dust emission of galaxies at all redshifts. The Herschel satellite (Pilbratt et al. 2010) provides far-infrared (FIR) photometric observations of local and distant galaxies thanks to its two onboard cameras, the Photodetector Array Camera and Spectrometer (PACS, 70−160 μm, Poglitsch et al. 2010) and the Spectral and Photometric Imaging REceiver (SPIRE, 250−500 μm, Griffin et al. 2010). The availability of a good SED coverage from 24 to 500 μm, obtained by combining Herschel and Spitzer observations, allows us to apply different types of galaxy emission models with the aim of deriving the main properties of dust, such as its mass, temperature, and, to some degree, composition. Different families of models, with different levels of complexity, can be broadly recognized: a) full radiation transfer treatments (e.g., Silva et al. 1998; Efstathiou et al. 2000; Piovan et al. 2006; Siebenmorgen & Krügel 2007) requiring geometric assumptions, very detailed observed SEDs, and large amounts of computational power; b) evolutionary and mixed stellar population synthesis, including FIR emission (e.g., da Cunha et al. 2008); c) physically motivated FIR dust emission models (e.g., Draine & Li 2007; Galliano et al. 2011); d) simple or multiple modified blackbody (MBB) SED fitting; e) template families following locally calibrated relations (e.g., the L–T relation, Chary & Elbaz 2001; Dale & Helou 2002); and f) semiempirical template fitting (e.g., Polletta et al. 2007; Wuyts et al. 2011a; Berta et al. 2013b).

Focusing on dust emission, of particular relevance and very popular are types c) and d), i.e., SED fitting with physically motivated models, comprising several dust components or MBB fitting under the assumption that dust emissivity can be described as a simple frequency power law and dust is in thermal equilibrium.

Dunne et al. (2011) applied MBB SED fitting to z< 0.5 H-ATLAS (Eales et al. 2010) star-forming galaxies selected at 250 μm and studied the evolution of the dust mass function over the last 5 billion years, showing that the dust/stellar mass ratio was three to four times larger at z = 0.4−0.5 than today. At these redshifts, PACS and SPIRE data cover the FIR peak of the SED and part of its Rayleigh-Jeans tail (RJ), thus leading to reliable M_dust estimates.

Cortese et al. (2012) combine Herschel photometry and multiwavelength data to build detailed SEDs of nearby galaxies from the Herschel Reference Survey (HRS; Boselli et al. 2010) and fit them with MBB and Draine & Li (2007, DL07) models. In combination with radio data, these authors study the dust-to-Hi and dust-to-stellar mass ratios of their targets, finding that the former increases as a function of stellar mass, while the latter tends to decrease as a function of stellar mass. Ciesla et al. (2014) apply DL07 modeling to HRS SEDs and have built a new set of reference, local SED templates.

Dale et al. (2012) and Bianchi (2013) apply DL07 and MBB fitting to the FIR-submm SEDs of KINGFISH (Kennicutt et al. 2011) local galaxies. The former authors show evidence of sub-mm (500 μm) flux excess in dwarf galaxies with respect to the expectation based on the extrapolation from Spitzer data. They also find a factor ~2 difference between DL07- and MBB-based dust masses, while Bianchi (2013) shows that this difference is avoided if using consistent assumptions for the two approaches.

Rémy-Ruyer et al. (2013) combine KINGFISH and other samples of nearby galaxies observed by Herschel, including direct measurements of molecular gas mass, to study the gas/dust vs. metallicity relation over a 2 dex metallicity range. Dust masses are in this case estimated with Galliano et al. (2011) models. Hunt et al. (2015) come back to the KINGFISH sample and present the surface brightness profiles of these nearby galaxies. These authors again study the spatially resolved SEDs and resolved properties of dust with DL07 and MBB modeling.

Eales et al. (2012) foresee the possibility of applying a single conversion factor to derive gas masses from sub-mm fluxes, based on a Milky Way (MW) observation and recalibrated on HRS nearby galaxies with CO and HI observations. They point out that this method suffers from the limitation that their sample is not necessarily representative of all galaxy populations and metallicity dependencies play an important role both for dust- and CO-based masses.

Successfully applied to local galaxies, the SED fitting techniques and scaling relations mentioned above have also been adopted by several authors to study the IR emission of distant galaxies. Broadly speaking, two main approaches can be identified: using a single sub-mm continuum observation to scale it to dust or gas mass using known relations (e.g., Eales et al. 2012; Scoville et al. 2014); or characterizing the galaxy SED including the FIR peak (e.g. Dale et al. 2012; Galametz et al. 2013). Nevertheless, when dealing with high-z galaxies, the use of these techniques is not as straightforward as in the nearby Universe. It is necessary not only to keep in mind that the properties of dust and gas might be evolving as a function of time, but also to face the limitations of the available data.

In the simplest case of all, Eales et al. (2012) discuss the limitations of SED coverage at high redshift, and show that small dust temperature gradients (see also Hunt et al. 2015) within a galaxy can lead to incorrect measurements when using a MBB fit to the global SED. Galametz et al. (2014) point out the effects of limited angular resolution on the derivation of model parameters.

Scoville et al. (2014) apply a similar concept to the RJ side of the SED of local and high-z galaxies. They suggest that gas masses can be measured using single-band sub-mm observations, adopting dust emissivities calibrated on Planck observations of the MW. They apply this method to early ALMA observations of four small samples of galaxies at four different cosmic epochs. Using simple MBB simulations, Genzel et al. (2015) discuss the systematic effects induced by this method and show that an approach involving two sub-mm bands and a proper MBB fit would be preferable.

Magdis et al. (2011) and Magnelli et al. (2012a) study the SEDs of star-forming and sub-mm galaxies that benefit from Spitzer, Herschel, and sub-mm photometry. Modeling of MBB and DL07 is applied to their SEDs with the aim of studying the dust properties of these objects and comparing these to CO observations. In this case, Magdis et al. (2011) combine dust and CO observations to obtain a measurement of the α_CO CO-to-gas conversion factor. Their method is then further refined by Magdis et al. (2012) on a larger sample of sources via stacking on FIR and sub-mm maps.

Stacking is also employed by Santini et al. (2014), who bin 24 μm-selected galaxies in a grid of z-M^∗-SFR and study the evolution of dust/stellar mass ratio and star formation efficiency (SFE) as a function of redshift (see also Genzel et al. 2015). These authors find no evidence that the M_dust evolves with redshift at a given M^∗ or SFR. They find that the SFE at z ~ 2 was a factor ~5 higher than at z ~ 0. Béthermin et al. (2015) also include sub-mm data in their stacking analysis and thus improve on the application of DL07 models to the average SEDs of distant galaxies.

Despite the wide use of SED fitting to derive dust and gas properties of distant galaxies, the underlying assumptions and implications of instrumental limitations are rarely highlighted nor discussed. In this work, we aim to show how real-world limitations affect dust and gas mass determinations based on DL07 and MBB SED fitting. To achieve our goal, we use the deepest PACS data available to date, namely the GOODS fields observations belonging to the PACS Evolutionary Probe (PEP, Lutz et al. 2011) and the GOODS-Herschel (Elbaz et al. 2011) surveys combined with SPIRE photometry from the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). The analysis consists of measurements carried out on real sources (individual galaxies and stacked photometry) and Monte Carlo simulations.

The paper is organized as follows. Sections 2 and 3 present the methods adopted to derive gas and dust masses, including the details of DL07 and MBB SED fitting. In Sect. 4 the different collections of data used in the analysis are described. Monte Carlo simulations and their results are presented in Sects. 5 and 6 in two flavors: analyzing M_dust uncertainties and systematics in the z-M^∗-SFR space and as a function of rest-frame spectral coverage. Real-world limitations and prospects for ALMA observations are discussed in Sect. 7. Finally, Sect. 8 summarizes our findings.

Throughout this manuscript, we adopt a Λ-CDM cosmology with (h₀, Ω_m, Ω_Λ) = (0.71, 0.27, 0.73); a Chabrier (2003) initial mass function (IMF); and the Pettini & Pagel (2004) metallicity scale. Moreover, as described in the next sections, we make use of the parameterization of the M^∗-Z relation by Genzel et al. (2015) and of the definition of reference (main) sequence of star formation (MS) by Whitaker et al. (2014). By the term gas mass, we mean the molecular gas component of galaxies unless otherwise specified.

2. Deriving gas mass

Molecular gas is generally quantified using measurements of the intensity of rotational transitions of the CO molecule and adopting a CO-to-H₂ conversion factor (see Bolatto et al. 2013, for a review). Local samples include several hundred targets (e.g., Saintonge et al. 2011), but CO detections of star-forming galaxies at intermediate- or high-redshift are still limited to modest statistics and mostly to luminous galaxies (Carilli & Walter 2013, and references therein), although they increasingly reach so-called main-sequence objects (Tacconi et al. 2013). Genzel et al. (2015) collected CO line emission measurements for 484 star-forming galaxies at z = 0−3, including more than 200 z ≤ 0.05 galaxies from the COLDGASS (Saintonge et al. 2011) survey and more than 80 at z = 0.5−1.5, z = 2.0−2.5 from the PHIBSS and PHIBSS2 (Tacconi et al. 2013, and in prep.) surveys.

Alternative, faster ways to gain gas masses for distant galaxies involve the use of locally calibrated relations linking gas mass to dust mass through gas/dust ratio scaling as a function of metallicity or to SFR through the scaling of depletion times. In this work, we extensively employ the former approach, based on the gas/dust ratio, and compare the results to CO- and to t_depl-based results.

2.1. Scaling of depletion times

Based on current CO observations of local and z ~ 1 objects (e.g., Saintonge et al. 2011, 2012; Tacconi et al. 2013), and integrating the relation found by Schmidt (1959) and Kennicutt (1998b) between volume (surface) density of star formation and of molecular gas mass, the relation between the total M_gas and SFR of MS galaxies can be described as a simple scaling with a depletion timescale mildly dependent on redshift ${\mathit{M}}_{\mathrm{gas}}\mathit{/}\mathit{SFR}\mathrm{=}{\mathit{\tau }}_{\mathrm{dep}}\mathrm{,}$(1)referred to as the integrated Kennicutt-Schmidt relation (KS) in what follows. Based on CO observations of z ~ 1.2 and 2.2 galaxies, compared to local galaxies, Tacconi et al. (2013) derived a dependence of τ_dep on redshift of ${\mathit{\tau }}_{\mathrm{dep}}\mathrm{=}\mathrm{1.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-1}\mathrm{[}\mathrm{Gyr}\mathrm{]}\mathrm{.}$(2)Based on their collection of ~500 CO measurements and on dust mass derivations, Genzel et al. (2015) derived a three-parameter expression of τ_dep, depending on stellar mass, M^∗, SFR, and redshift. We use the latest update of this relation (R. Genzel, priv. comm.): $\begin{array}{ccc}\log \mathrm{(}{\mathit{\tau }}_{\mathrm{dep}}\mathrm{)}& \mathrm{=}& \mathrm{0.13}\mathrm{-}\mathrm{0.37}\mathrm{\times }\log \mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}^{\mathrm{)}}\\ & & \mathrm{-}\mathrm{0.43}\mathrm{\times }\log \mathrm{(}\mathrm{s}\mathit{SFR}\mathit{/}{\mathrm{s}\mathit{SFR}}_{\mathrm{MS}}\mathrm{)}\\ & & \end{array}$(3)adopting the Whitaker et al. (2014) definition of MS and the Milky Way value of the $\mathit{\alpha }{}^{\mathrm{\right(}}\mathrm{CO}^{\mathrm{\left)}}$ conversion factor between CO luminosity and molecular gas mass.

Here, star formation rates include both IR and UV contributions. The SFR_IR,UV are computed with the Kennicutt (1998a) calibrations, scaled to the Chabrier (2003) IMF. Infrared luminosities between 8 and 1000 μm, L(IR), are derived by integrating the results of FIR SED fitting (see Sect. 3.3). An alternative L(IR) estimate was obtained by fitting SEDs with the Berta et al. (2013b) templates library, leading to equivalent results.

2.2. The dust method

As described by Magdis et al. (2011, 2012) and Eales et al. (2012), and as widely found in the literature, it is possible to derive gas masses of (distant) galaxies by relying on measurements of dust mass, M_dust, to be transformed into M_gas, once the gas-to-dust ratio (GDR or δ_GDR) of the galaxy is known, i.e., ${\mathit{M}}_{\mathrm{gas}}\mathrm{=}{\mathit{M}}_{\mathrm{dust}}\mathrm{\times }{\mathit{\delta }}_{\mathrm{GDR}}\mathit{.}$(4)A number of assumptions can be adopted to estimate the GDR, from simply adopting the MW value (which turns out to work well for local Virgo cluster galaxies; Eales et al. 2012), to adopting observed relations, such as the Leroy et al. (2011) local dependence of GDR on metallicity. Here we use a revisited form of the Leroy et al. (2011) relation, consistently recalibrated to the Pettini & Pagel (2004, PP04) metallicity scale by Magdis et al. (2012) as follows: $\log \mathrm{\left(}{\mathit{\delta }}_{\mathrm{GDR}}\mathrm{\right)}\mathrm{=}\left(\mathrm{10.54}\mathrm{\pm }\mathrm{1.0}\right)\mathrm{-}\mathrm{(}\mathrm{0.99}\mathrm{\pm }{\mathrm{0.12}}^{\mathrm{)}}\mathrm{\times }\left(\mathrm{12}\mathrm{+}\log \mathrm{(}\mathrm{O}\mathit{/}{\mathrm{H}}^{\mathrm{)}}\right)\mathit{,}$(5)with a scatter of 0.15 dex. In absence of direct spectroscopic measurements, galaxy metallicity can be derived from the known stellar mass-metallicity relation and its redshift dependence. Recently, Genzel et al. (2015) combined the M^∗-Z parameterizations of Erb et al. (2006), Maiolino et al. (2008), Zahid et al. (2014), and Wuyts et al. (2014) into a single relation $\mathrm{12}\mathrm{+}\log {\mathrm{(}\mathrm{O}\mathit{/}{\mathrm{H}}^{\mathrm{)}}}_{\mathrm{PP}\mathrm{04}}\mathrm{=}\mathit{a}\mathrm{-}\mathrm{0.087}\mathrm{\times }{\left(\log {\mathit{M}}^{\mathrm{\ast }}\mathrm{-}\mathit{b}\right)}^{\mathrm{2}}\mathit{,}$(6)where a = 8.74 and b = 10.4 + 4.46 × log ⁽1 + z⁾−1.78 × [log ⁽1 + z⁾]². In what follows, we adopt this relation, which provides metallicities in the PP04 scale. Alternative parameterizations or the so-called M^∗-Z-SFR fundamental relation (Mannucci et al. 2010, 2011) give similar results, for z ≤ 2 massive star-forming galaxies with Z ~ Z_⊙, which dominate the samples discussed below, modulo the adopted metallicity scale.

We stress that it is very important to use a self consistent set of relations, i.e., the δ_GDR vs. Z scaling and the M^∗ vs. Z relation must be calibrated to the same metallicity scale to produce meaningful results.

3. Deriving dust mass through FIR SED fitting

Dust mass can be estimated by fitting the FIR SED of galaxies in several ways, e.g., by means of a MBB function or with dust models such as those prepared by Draine & Li (2007, DL07). In this section we provide more details about these two approaches.

3.1. Modified blackbody

The FIR SED of galaxies can be simplistically reproduced with a single-temperature modified blackbody, assuming that emission comes from a single-temperature δ distribution of dust. In this case, the rest-frame SED is modeled as (e.g., Blain et al. 2002) ${\mathit{L}}_{\mathit{\nu }}\mathrm{=}\mathrm{\Omega }{\mathit{ϵ}}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathit{,}$(7)where Ω is the solid angle of emission, ϵ_ν is the emissivity coefficient, and B_ν(T_dust) is the Planck function ${\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{h}}{{\mathit{c}}^{\mathrm{2}}}\frac{{\mathit{\nu }}^{\mathrm{3}}}{\exp \mathrm{(}{\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{dust}}}}^{\mathrm{)}}\mathrm{-}\mathrm{1}}$(8)in units of [erg s^-1 Hz^-1 m^-2 sr^-1].

For a uniform medium (e.g., a dust cloud in our case) of optical depth, τ_ν, radiative transfer implies that radiation is reduced by a factor exp⁽−τ_ν⁾ (e.g., Rybicki & Lightman 1979). Therefore, we can also write ϵ_ν = 1−exp⁽−τ_ν⁾ (e.g., Benford et al. 1999; Omont et al. 2001). We allow ν₀ to be the frequency where τ = τ₀ = 1. The optical depth at the generic frequency ν can be approximated as a power law ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathrm{(}\mathit{\nu }\mathit{/}{\mathit{\nu }}_{\mathrm{0}}\mathrm{)}}^{\mathit{\beta }}\mathrm{.}$(9)Therefore, ${\mathit{ϵ}}_{\mathit{\nu }}\mathrm{=}\mathrm{1}\mathrm{-}\exp \mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathit{\nu }}\mathrm{)}\mathrm{=}\mathrm{1}\mathrm{-}\exp \mathrm{[}\mathrm{-}{\mathit{\tau }}_{\mathrm{0}}\mathrm{(}\mathit{\nu }\mathit{/}{\mathit{\nu }}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}\mathrm{]}\mathrm{.}$(10)On the other hand, the frequency-dependent optical depth can be written as (e.g., Beelen et al. 2006) ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}{\mathrm{\int }}_{\mathit{s}}\mathit{\rho }\hspace{0.17em}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}\mathrm{,}$(11)where κ_ν is the mass absorption coefficient of dust at rest frequency ν and ρ is the total mass density. If the emission is optically thin (i.e., τ ≪ 1), we can expand ϵ_ν in its Taylor series expression¹ and thus obtain ${\mathit{ϵ}}_{\mathit{\nu }}\mathrm{~}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{M}}_{\mathrm{dust}}\mathrm{.}$(12)Therefore, the emergent luminosity from a given dust mass is given by ${\mathit{L}}_{\mathit{\nu }}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{M}}_{\mathrm{dust}}{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}\mathrm{,}$(13)where we have assumed that emission is isotropic over a spherical surface and κ_ν = κ₀⁽ν/ν₀^)β. For a galaxy FIR dust emission, β is in the range 1.5−2.0 (e.g., Magnelli et al. 2012a, and references therein). Here M_dust is in units of [kg], κ_ν is given in [m² kg^-1], and therefore L_ν turns out to have units of [erg s^-1 Hz^-1].

Making the dependence of frequency explicit, Eq. (13) can be rewritten as ${\mathit{L}}_{\mathit{\nu }}\mathrm{\propto }\frac{{\mathit{\nu }}^{\mathrm{3}\mathrm{+}\mathit{\beta }}}{\exp \mathrm{(}{\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{dust}}}}^{\mathrm{)}}\mathrm{-}\mathrm{1}}\mathrm{,}$(14)which is the usual simplified form of the modified blackbody function.

Inverting Eq. (13), it is possible to obtain an estimate of dust mass, once the rest-frame luminosity of a galaxy at a given FIR wavelength is known, a MBB fit to the data is possible, and a value of κ_ν (or κ₀ and ν₀) for dust is adopted, as ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}\frac{{\mathit{L}}_{\mathit{\nu }}}{\mathrm{4}\mathit{\pi }{\mathit{\kappa }}_{\mathit{\nu }}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{dust}}\mathrm{\right)}}\mathrm{·}$(15)This holds when working at rest-frame frequencies and needs to be transformed into an observed frame in the usual way. Finally, combined with the knowledge of the gas-to-dust ratio, δ_GDR, Eq. (15) leads to an estimate of gas mass.

3.2. The use of κ_ν

To compute dust masses through a modified blackbody fit, one needs to assume a value of κ_ν at a given wavelength, possibly within the rest-frame wavelength range covered by actual photometric data, or assume a value of κ₀ at the reference frequency ν₀, where τ_ν = τ₀ = 1.

One viable approach is to adopt the set of κ_ν computed by Li & Draine (2001, see their Table 6), and either use the value appropriate for the given rest-frame wavelength from their Table (or interpolated) or apply Eq. (9).

It is common (e.g., Magdis et al. 2012; Magnelli et al. 2012b) to follow the second approach, using the values of κ_ν tabulated by Li & Draine (2001) and applying κ_ν = κ₀⁽ν/ν₀^)β. In this procedure, it is important to properly apply k correction to all terms in Eq. (15) and use the correct values of κ₀, ν₀, and β. In fact, Eq. (9) assumes that the slope β used while fitting the observed SEDs (see Eq. (14)) is the same β adopted in the computation of κ_ν.

Often, β is let free to vary while fitting, if enough datapoints are available, or if there are not enough available, a fixed value (e.g., β = 1.5) is adopted. On the other hand, the value of β adopted by Draine & Li (2007) is ~2 (see Cortese et al. 2012) and the values of κ_ν in Table 6 of Li & Draine (2001; in its electronic version with the old 2001 computation and for the MW dust mixture with q_PAH = 0.47) turn out to show β_LD01,a = 2.07 between λ = 100 and 600 μm, and β_LD01,b = 1.68 at λ> 600μm.

The direct consequence of the mismatch between the adopted value of β and the β ~ 2 describing the actual Li & Draine (2001)κ_ν is that one should apply an additional correction factor when applying Eq. (9) (i.e., Eq. (6) in Magdis’ paper and Eq. (4) in Magnelli’s paper) or there will be a dependence of M_dust on the value of ν₀ used. This dependence on frequency has a power of roughly −0.5 when adopting β = 1.5.

Bianchi (2013) points out this effect and also shows that using the κ_ν values tabulated by Draine (2003), i.e., those actually used in the Draine & Li (2007) models, is more appropriate and leads to consistent results when comparing MBB-based and DL07-based dust masses. Also, Bianchi points out that adopting a different value of β ≠ 2.08 leads to a dependence of model normalizations and, hence, of M_dust, on the reference wavelength adopted when computing κ_ν starting from κ₀. Finally, Bianchi (2013) shows that after accounting for the needed corrections on κ_ν, a mild dependence on blackbody shapes is still left because that the best-fit dust temperature also depends on β (see also Bianchi et al. 1999).

In practice, these subtle differences in the treatment of κ_ν and β have led to MBB-based dust masses changing by a factor 3−5 for the same SED and the same basic reference dust opacity (see also Sect. 7.2).

3.3. Draine & Li (2007) models

A more sophisticated and physically motivated approach to FIR SED fitting and derivation of dust masses is to adopt the Draine & Li (2007, DL07) dust models, which are an upgrade of those originally developed by Draine & Li (2001), Li & Draine (2001), and Weingartner & Draine (2001).

In short, these models describe interstellar dust as a mixture of carbonaceous and amorphous silicate grains, whose size distributions are chosen to mimic the observed extinction law in the MW, Large Magellanic Cloud (LMC), or Small Magellanic Cloud (SMC) bar region. As described by Draine & Li (2007), carbonaceous grains have the properties of polycyclic aromatic hydrocarbons (PAH) molecules and ions when the effective radius is a< 5.0 [nm], the properties of graphite spheres when a ≫ 10 [nm], and optical properties intermediate between those of PAH particles and graphite particles for 5 <a< 10[nm]. The ionization fraction x_ion(a) of the PAH particles is assumed to be the average for the diffuse ISM. The properties of grains are parameterized by the PAH index q_PAH, which is defined as the fraction of dust mass in the form of PAH molecules.

The dust distribution is divided in two components: the diffuse interstellar medium (ISM), usually constituting the majority of the dust budget and dust enclosed in photodissociation regions (PDRs). The former is heated by a radiation field of constant intensity U_min. The latter, representing a fraction γ of the total amount of dust, is exposed to starlight with intensities ranging from U_min to U_max. Although PDRs usually provide a small fraction of the total dust mass, they can contribute to the majority of dust emission in mid-IR SEDs. The mass dM of dust exposed to starlight intensities in the range [U,U + dU] is given by the following power-law distribution: $\mathrm{d}{\mathit{M}}_{\mathrm{dust}}\mathrm{=}\mathrm{const}\mathrm{\times }{\mathit{U}}^{-2}\hspace{0.17em}\mathrm{d}\mathit{U}\mathrm{for}{\mathit{U}}_{\min }\mathit{<}\mathit{U}\mathit{<}{\mathit{U}}_{\max }\mathrm{.}$(16)The diffuse ISM is modeled by setting U_max = U_min.

Summing the diffuse and PDR components, we obtain $\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{=}\left(\mathrm{1}\mathrm{-}\mathit{\gamma }\right)\mathit{\delta }\mathrm{(}\mathit{U}\mathrm{-}{\mathit{U}}_{\min }\mathrm{)}\mathrm{+}\mathit{\gamma }{\mathit{M}}_{\mathrm{dust}}\frac{\mathit{\alpha }\mathrm{-}\mathrm{1}}{{\mathit{U}}_{\min }^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathrm{-}{\mathit{U}}_{\max }^{\mathrm{1}\mathrm{-}\mathit{\alpha }}}{\mathit{U}}^{\mathrm{-}\mathit{\alpha }}\mathrm{,}$(17)for α ≠ 1; here δ is the delta function representing the diffuse interstellar radiation field of intensity U = U_min = U_max.

Using the numerical methods described in Draine & Li (2001), Li & Draine (2001), and Draine & Li (2007), temperature distribution functions are computed for all particles small enough for quantized heating to be important; large grains are treated as having steady-state temperatures determined by starlight heating and radiative cooling equilibrium. With temperature distributions and dust absorption cross sections, the time-averaged IR emission for a given grain size and type is computed, and finally the power radiated per unit frequency and unit mass for the given dust mix exposed to the radiation field U, is computed summing over all grain types and sizes. Integrating between U_min and U_max, this gives the specific power per unit mass p_ν(q_PAH,U_min,U_max,α) (see Eq. (9) in Draine et al. 2007).

Draine & Li (2007) thus propose to fit the emission of galaxies through a linear combination of emission from diffuse ISM dust and PDRs emission. The emission spectrum, expressed as emissivity per hydrogen nucleon, ${\mathit{j}}_{\mathit{\nu }}\mathrm{=}{}^{\mathrm{\right(}}\frac{{\mathit{M}}_{\mathrm{dust}}}{{\mathit{M}}_{\mathrm{H}}}^{\mathrm{\left)}}{\mathit{m}}_{\mathrm{H}}\frac{{\mathit{p}}_{\mathit{\nu }}}{\mathrm{4}\mathit{\pi }}$, can then be simply written as $\begin{array}{ccc}{\mathit{j}}_{\mathit{\nu }}\left({\mathit{q}}_{\mathrm{PAH}}\mathit{,}{\mathit{U}}_{\min }\mathit{,}{\mathit{U}}_{\max }\mathit{,\alpha ,\gamma }\right)& \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\gamma }\mathrm{)}{\mathit{j}}_{\mathit{\nu }}\left[{\mathit{U}}_{\min }\mathit{,}{\mathit{U}}_{\min }\right]\\ & & \mathrm{+}\mathit{\gamma }{\mathit{j}}_{\mathit{\nu }}\left[{\mathit{U}}_{\min }\mathit{,}{\mathit{U}}_{\max }\mathit{,\alpha }\right]\mathrm{,}\end{array}$(18)in units of [erg s^-1 Hz^-1 sr^-1 H^-1]; j_ν is the quantity contained in the model files available online².

As shown by Draine et al. (2007), the total dust luminosity is given by ${\mathit{L}}_{\mathrm{dust}}\mathrm{=}\mathrm{⟨}\mathit{U}\mathrm{⟩}{\mathit{P}}_{\mathrm{0}}{\mathit{M}}_{\mathrm{dust}}\mathit{,}$(19)where the dust-weighted mean starlight intensity, or mean radiation field, scale factor is $\mathrm{⟨}\mathit{U}\mathrm{⟩}\mathrm{=}\left[\left(\mathrm{1}\mathrm{-}\mathit{\gamma }\right){\mathit{U}}_{\min }\mathrm{+}\frac{\mathit{\gamma }\ln \mathrm{(}{\mathit{U}}_{\max }\mathit{/}{\mathit{U}}_{\min }\mathrm{)}}{{\mathit{U}}_{\min }^{-1}\mathrm{-}{\mathit{U}}_{\max }^{-1}}\right]\mathrm{,}$(20)if α = 2 as in our assumptions and P₀ is the power absorbed per unit dust mass in a radiation field U = 1.

In principle, the model includes six free parameters: q_PAH; U_min; U_max; α; γ; M_dust. Dust mass is basically the model’s normalization. In practice, studying local galaxies in the Spitzer Nearby Galaxy Survey (SINGS) and fitting their FIR SEDs, Draine et al. (2007) demonstrate that some of the parameters can be limited to a restricted range of values or even fixed to a single value. The overall fit is not very sensitive to the dust model adopted (MW, LMC, SMC), and can be limited to MW models alone, implying only seven values of q_PAH. Also, the actual values of α and U_max do not significantly influence the performance of models, and fixing them to α = 2 and U_max = 10⁶ is a reasonable choice to successfully reproduce SINGS galaxies. Finally, since small values of U_min correspond to dust temperatures below 15 K that cannot be constrained by FIR photometry alone, in the absence of rest-frame sub-mm data, they suggest limiting U_min to the range 0.7 ≤ U_min ≤ 25. A direct consequence is a possible underestimate of dust mass, if a significant amount is stored in a very cold component, but Draine et al. (2007) conclude that omitting sub-mm data increases the scatter on M_dust by 50%, but does not induce any significant systematic offset. Most observational papers in the literature (e.g., Magdis et al. 2012; Magnelli et al. 2012a; Santini et al. 2014) adopt this optimized choice of parameters, relying on a setup that has only been thoroughly verified for local galaxies near solar metallicity (Draine et al. 2007).

Figure 1 shows how varying each parameter (q_pah, U_min, γ) independently, and fixing the others, induces modifications on the shape of the modeled SED. We also indicate also a set of IR rest-frame bands, which are available in surveys of nearby or distant galaxies.

The choice of the actual PAH abundance, parameterized by q_PAH and limited to MW models, mostly influences the short-wavelength SED (1−20 μm), and is not making a difference at λ> 20μm. Therefore FIR surveys of galaxies are not very sensitive to this parameter and cannot put tight constraints on its value.

The parameter U_min regulates the range of starlight intensities that are heating dust. Since the amount of dust exposed to starlight intensities in between U and U + dU is modeled as a power law ∝U^{− α}, the smaller U_min, the larger is the amount of dust subject to low-energy radiation. As a consequence, the dust component dominating the SED is colder for smaller values of U_min and the FIR peak shifts to longer wavelengths; at the same time, at a given wavelength and for a given total dust mass, the model intensity decreases. Thus, we expect degeneracies and possibly some amount of systematics between M_dust and U_min.

The fraction γ of dust locked in PDRs influences the repartition of the emitted energy between the U_max = U_min component (ISM) and the rest. When fitting the SEDs of local SINGS galaxies, Draine et al. (2007, see their Table 4) found that the SEDs of vast majority of galaxies were successfully reproduced with γ ≤ 0.15 with only a couple of cases extending to γ ~ 0.3−0.4.

Once the FIR SED is fitted with these models, its dust mass can be computed as ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}\left(\frac{{\mathit{M}}_{\mathrm{dust}}}{{\mathit{M}}_{\mathrm{H}}}\right){\mathit{m}}_{\mathrm{H}}\frac{{\mathit{L}}_{\mathit{\nu }}}{\mathrm{4}\mathit{\pi }{\mathit{j}}_{\mathit{\nu }}}\mathrm{,}$(21)where the right-hand term is evaluated at a given rest-frame frequency ν of choice, which is at best covered by observations; L_ν is the rest-frame luminosity at frequency ν; m_H is the mass of the hydrogen nucleon; and M_dust/M_H is the dust-hydrogen ratio characteristic of the adopted dust model (and tabulated together with q_PAH). Making sure to use the right units for all terms, we use m_H = 1.67262178 × 10^-27[kg], L_ν in units of [erg s^-1 Hz^-1], and therefore dust mass turns out to have units of [kg].

4. Available data

The PEP survey (Lutz et al. 2011) has covered the most popular blank fields with the PACS (Poglitsch et al. 2010) instrument aboard Herschel (Pilbratt et al. 2010) at 100 and 160 μm. In addition, the same fields have been observed by the HerMES survey (Oliver et al. 2012) with SPIRE (Griffin et al. 2010) at 250, 350, 500 μm. Furthermore, the GOODS-Herschel (Elbaz et al. 2011) survey has provided deeper PACS coverage of the GOODS fields, reaching the confusion limit at 100 and 160 μm. Finally, the PEP team has combined all the available PACS data in the GOODS-N and GOODS-S fields, including science verification (SV) and ECDFS observations, producing the deepest FIR maps obtained so far (Magnelli et al. 2013).

Here we make use of PEP first data release (DR1³, Lutz et al. 2011), the HerMES DR2 and DR3⁴ (Oliver et al. 2012; Roseboom et al. 2010, 2012), and combined PEP + GOODS-H data (included in PEP DR1, Magnelli et al. 2013).

4.1. Individual sources

The deepest PACS fields, GOODS-N and GOODS-S benefit from an extensive coverage at all wavelengths. We use the PACS 70, 100, 160 μm data presented by Magnelli et al. (2013); the SPIRE HerMES data (Oliver et al. 2012; Roseboom et al. 2010); the multiwavelength catalogs built by Berta et al. (2011) and by the MUSIC team (Grazian et al. 2006); the collection of spectroscopic redshift by Barger et al. (2008), Balestra et al. (2010), and Berta et al. (2011); and the photometric redshifts by Berta et al. (2011) and Wuyts et al. (2011a). We refer to the dedicated publications for more details on each dataset and to Berta et al. (2011, 2013b) for an overview.

We select objects with at least a 3σ detection in the PACS 160 μm band, which for the redshift range of interest turns out to be the best single-band proxy of star formation rate (Elbaz et al. 2011; Nordon et al. 2013).

To this generic sample of ~1400, 160 μm sources, we add the list of 61 sub-mm galaxies (SMGs) compiled by Magnelli et al. (2012a) and distributed in the main PEP fields. For these sources, a rich multiwavelength dataset is available, ranging from the optical to the sub-mm as compiled by Magnelli et al. (2012a).

4.2. Stacked photometry

Magnelli et al. (2014) studied the dust temperature, T_dust, in star-forming galaxies as a function of their position in the z-M^∗-SFR space.

These authors used the M^∗, SFR, z estimates by Wuyts et al. (2011a,b) in the GOODS-N, GOODS-S, and COSMOS fields. Stellar masses are based on optical-NIR SED fitting carried out adopting Bruzual & Charlot (2003) templates. Star formation rates are based on a ladder of indicators (SED fitting, mid-IR photometry, FIR photometry), calibrated on Herschel observations of ~7000 galaxies in PEP fields. Redshifts are the combination of a collection of spectroscopic measurements and photometric estimates (see above).

Magnelli et al. (2014) produced a grid in the M^∗-SFR-z space, binned a K_s-selected sample accordingly, and finally produced an average SED (70, 100, 160, 250, 350, 500 μm) for each bin by means of stacking on Herschel maps.

We make use of the stacked photometry by Magnelli et al. (2014), plus similar Spitzer/MIPS 24 μm data, taking care to limit our analysis to those regions of the M^∗-SFR-z space, where good stacked SEDs are available (see Figs. 4 and 5 in Magnelli et al. 2014). Table 1 summarizes the z, M^∗, sSFR ranges of interest.

4.3. CO-detected sources

A direct measurement of M_gas comes from the intensity of spectral lines of molecular gas tracers, such as CO, modulo the conversion factor from the given molecular species to total gas. Here we collect the FIR and sub-mm photometry of galaxies that were observed in CO-line spectroscopy to compare CO-based and dust-based M_gas estimates directly.

The following samples of CO-detected objects are taken into account:

• PHIBSS z ~ 1 galaxies (Tacconi et al. 2013);

• other star-forming galaxies by Magnelli et al. (2012b) and Daddi et al. (2010);

• lensed galaxies (Saintonge et al. 2013); and

• sub-mm galaxies (SMGs) by Bothwell et al. (2013).

Appendix A describes the publicly available data for each case. In summary, the number of CO-detected objects that have enough FIR photometry to allow a DL07 SED fitting and that have a M^∗ estimate is: 23 sources in EGS from the PHIBSS catalog, 5 BzK galaxies by Daddi et al. (2010), six sources by Magnelli et al. (2013), ten lensed galaxies by Saintonge et al. (2013), and seven SMGs by Bothwell et al. (2013); this is a total of 51 objects.

Stellar masses of these CO-detected sources come from the respective publications dealing with each sample. They were all computed with Chabrier (2003) IMF and BC03 models or are consistent with these assumptions. For SMGs, the adopted stellar masses are similar to the scale of Hainline et al. (2011). The adopted cosmological parameters are all in line with those assumed in this work.

5. Accuracy of M_dust and T_dust, as allowed by Herschel surveys

Magnelli et al. (2014) provide a tool to derive T_dust of a galaxy once its redshifts, M^∗ and SFR are known, as well as to obtain its expected SED based on the Dale & Helou (2002) template library. Genzel et al. (2015) re-analyze the data by Magnelli et al. (2014) and produce new scalings linking T_dust and M_gas to M^∗, sSFR, and z.

In this section, we would like to study how the observational limitations of Herschel photometry affect the derivation of M_gas and T_dust. To this aim, we adopt the signal-to-noise ratios (S/Ns) of two case studies: the deepest field of the PEP survey, i.e., GOODS-S; and the widest, but ~8 times shallower, field of PEP, i.e., COSMOS. In what follows, we produce a list of artificial sources, characterized by T_dust, M_gas, and SEDs given by the relations by Magnelli et al. (2014) and Genzel et al. (2015), and fit them with MBB and DL07 models. We limit the analysis to the parameter ranges listed in Table 1. Our simulation is structured as follows:

• 1)

  We adopt the z-M^∗-sSFR grid by Magnelli et al. (2014) and we limit the analysis to the range of parameters over which the relation by Magnelli et al. (2014) holds (see Sect. 4.2).

• 2)

  The recipe by Magnelli et al. (2014) produces a value of T_dust given the position in the z-M^∗-sSFR grid along with the typical FIR SEDs of a galaxy in that position. These SEDs are based on (and limited to) the Dale & Helou (2002) templates library.

• 3)

  These SEDs are convolved with photometric filters: MIPS 24 μm; PACS 70, 100, and 160 μm; and SPIRE 250, 350, and 500 μm.

• 4a)

  In the case of the DL07 simulation, the convolved SEDs are fit with the DL07 models, adopting a ≤1% photometric uncertainty in all bands. In this way, a so-called input catalog is produced, defining the DL07 parameters and a reference M_dust to be associated with each z-M^∗-SFR bin.

• 4b)

  In the case of MBB simulation, Genzel et al. (2015) provide scaling relations to define the input value of T_dust and M_gas (to be transformed into M_dust adopting a gas/dust mass ratio) at any position in the M^∗-sSFR-z space. These relations are calibrated on Magnelli et al. (2014) data, which hold for β = 1.5. Thus, using the relation by Magnelli et al. (2014) leads to similar results. In case a value β = 2.0 was adopted, temperatures need to be increased by 4 K (Magnelli et al. 2014) and dust masses corrected as discussed in Sect. 3.2.

• 5)

  Real noise levels are associated with each band (see Magnelli et al. 2013; Berta et al. 2013b). We use the PEP/GOODS-H/HerMES noise values for individual detections of two different fields: GOODS-S (PEP plus GOODS-H and HerMES depths; see Magnelli et al. 2013; Oliver et al. 2012) and COSMOS (PEP and HerMES depths; see Lutz et al. 2011; Oliver et al. 2012). Table 2 lists the adopted depths. Two independent simulations are run: one for each set of depths. Only bands with S/N ≥ 3 are taken into account.

• 6)

  We fit the noisified catalog with DL07 or MBB models. In the DL07 case, photometric points at λ ≥ 8μm (rest frame) are considered; in the MBB case, only bands at λ> 50μm (rest frame) are used.

• 7)

  Solutions are found both through χ² minimization and through 1000 Monte Carlo (MC) realizations for each entry in the synthetic catalog and are obtained by letting the photometry vary within the “observed” noise. In the two cases, uncertainties are computed based on Δχ² or as the dispersion of all MC realizations, respectively. The two approaches lead to comparable results. The analysis and figures presented here are based on the MC approach.

5.1. Results of DL07: Relative uncertainties on M_dust

We first verify down to what precision our procedure is able to determine dust masses as a function of position in the redshift, M^∗, SFR space, and at the depths reached by Herschel for individual detections in the GOODS-S and COSMOS fields.

Figure 2 shows the M^∗-SFR plane in different redshift slices, color coding each bin on the basis of its average σ(M_dust) /M_dust value. Bins indicated in red have an average M_dust/σ(M_dust) smaller than 3. It was not possible to run SED fits for bins missing at the low SFR side because too few photometric bands are available. Figure 3 collapses this diagram along the M^∗ axis for the GOODS-S case, thus allowing for a more straightforward view as a function of redshift and SFR.

As expected, the relative uncertainty on M_dust becomes worse, as redshift increases (due to Malmquist bias) and SFR decreases. Redshift plays a double role by dimming fluxes and thus raising the luminosity threshold; and by pushing the rest-frame SED coverage to shorter wavelengths (see also Sect. 6). At the sensitivity of PEP and HerMES, dust masses are retrieved with a S/N ≥ 3 for galaxies on the MS of star formation down to M^∗ ~ 10¹⁰[M_⊙] up to z ~ 1. It is possible to obtain an estimate of M_dust at comparable stellar masses up to z ~ 2 for objects lying increasingly above the MS.

5.2. Results of DL07: Systematics on M_dust

It is now possible to study how well M_dust is recovered by comparing input and output dust masses. Figure 4 presents the mere comparison of input and output M_dust for all bins in the considered z-M^∗-SFR space. Scatter increases and the incidence of outliers becomes larger as the depth of the available bands become less balanced, i.e., in the GOODS-S simulation, where MIPS and PACS are much deeper than SPIRE.

The color coding in Fig. 4 is based on redshift, stellar mass, specific SFR, and on M_dust relative uncertainty. Critical outliers lie at intermediate to high redshift and are characterized by low M^∗. They turn out to be sources with few photometric points available and thus with limited wavelength coverage. At the adopted depths, it is not possible to derive M_dust with an accuracy better than 3σ for them, therefore, they are indicated in red in Fig. 4.

Ignoring these critical cases, milder systematic offsets seem still to occur for those bins with larger relative uncertainties on M_dust, which tend to lie preferentially below the 1:1 line. However, for these sources with M_dust measured at >3σ, such systematic offsets are well below a factor 2.

Figure 5 puts systematic offsets in the context of the z-M^∗-SFR space in the case of GOODS-S. With PACS and SPIRE depths strongly unbalanced, there is a trend to underestimate M_dust at intermediate sSFR. This effect is more critical at z ≃ 0.5−1.5, where SPIRE bands are most important to constrain the peak of the FIR SED. At lower redshift, deep PACS data provide a good constraint already, and at higher redshift the peak and M_dust are likely equally poorly constrained both in the “IN” and “OUT” cases.

Working on a small sample of high-redshift galaxies, Magdis et al. (2012) showed that the presence of photometric data at rest-frame wavelengths larger that 100−160 μm should avoid large systematics in the estimate of dust masses. Similar results were obtained locally by Ciesla et al. (2014). Imposing the condition to have at least one 3σ detection at λ_rest ≥ 160μm slightly mitigates these systematics detected in our simulations, but does not solve the problem.

5.2.1. The choice of templates and the role of β

The results presented above in Sects. 5.1 and 5.2 were obtained by fitting, with DL07 templates, the synthetic photometry computed by convolving the Magnelli et al. (2014) SEDs with photometric filters. It is worth recalling that Magnelli et al. (2014) assigned to each Dale & Helou (2002, DH02) template a value of dust temperature, T_dust, based on an MBB fit to the template itself. This T_dust was obtained fixing β = 1.5. Then they used DH02 templates to fit the Herschel stacked photometry of each z-M^∗-SFR bin, thus linking each position in this space to T_dust. In short, we are fitting the DH02-based synthetic photometry of each z-M^∗-SFR bin with DL07 templates.

The DL07 models were computed by adopting a fixed slope of the dust IR emissivity, β = 2.08. On the other hand, in DH02 models, β varies as a function of the intensity of the radiation field, U, following the relation β = 2.5−0.4log U at λ> 100 μm. This intensity varies in the range 0.3 <U< 1.e5, thus one obtains β = 2.5, 2.1, 1.3 for U = 1, 100, 1000.

We now test whether the systematic trends seen in Fig. 5 and described in Sect. 5.2 could be, at least partially, driven by the different adopted models in input and output. To this aim, the DL07 simulation (Sect. 5) is now repeated with a new photometry based on DL07 models themselves. This is carried out by convolving DL07 models with photometric filters at step number 4a) in the previous scheme, and adopting this new photometry during the SED fitting process instead of the photometry computed at step number 3).

Figure 6 briefly shows the results for GOODS-S (our worst case, see Sect. 5.2), i.e., reproduces Fig. 4 with the new setup. As expected, the majority of systematics are no longer found with the exception of the low S/N region of the parameters space, where dust masses tend to be overestimated. Their relative uncertainty on M_dust exceeds 33% (red points). Similar results are obtained at the COSMOS depth.

The lesson to be learned is that the value of these kind of simulations and their ability to describe systematics is limited to the ideal case that the SED shape adopted in the simulation, and to fit the real sources is the same as that of real-life galaxies. In our case, this is exemplified by the value of β. Small differences can lead to misleading results or, more worrying, to physical overinterpretations of numerical effects. Vice versa, the absence of systematics in these kind of simulations might not always be indicative of absence of systematics on real sources if the real SEDs differ from those adopted to model them.

5.3. Results of MBB: Constraints on T_dust

Following the scheme introduced for the DL07 simulations, we now analyze the results of MBB-based SED fitting of synthetic catalogs. The procedure followed is similar to the DL07 case and is described in Sect. 5. We fix the value of β to ease the comparison of our simulations to recent literature works and later to the analysis of individual real sources for which few photometric points are available.

We now deal with the case at COSMOS depth, which represents the shallowest and worst case scenario in this analysis. Figure 7 compares input and output dust temperatures. Data points in the four panels are color coded as a function of redshift, M^∗, sSFR, and T_dust relative error. Dust temperature is retrieved within few percents in most cases. Only few catastrophic failures are recorded: they are limited to bins with very poor S/N and, thus, with large uncertainties on T_dust. Only ~4% of the cases turn out to have T_dust overestimated by more than 10%. These statistics become even better at the GOODS-S depth.

When T_dust is not correct, M_dust can also be overestimated, up to a factor of ~5. Nevertheless, in these cases, the photometric S/N is not good enough to guarantee a relative error on M_dust smaller than 3 σ.

5.3.1. Biases in MBB T_dust determination

The SEDs of real galaxies are indeed not single-temperature modified blackbodies: as seen before, the FIR emission is produced by a mixture of dust grains, of different sizes and shapes, which turn into a nontrivial distribution of temperatures. Although convenient, the MBB approximation is an over-simplification in most of the cases. More complex descriptions include using multitemperature blackbodies, DL07 modeling (see previous sections), or mixed approaches as that by Magnelli et al. (2014, described in Sect. 5.2.1. In this section we study possible biases in MBB-based T_dust due to changes in β and to redshift.

As in Magnelli et al. (2014), each template in the DH02 library is convolved with PACS/SPIRE passbands (70 to 500 μm) and the broadband photometry thus obtained is fitted with MBB models with varying dust temperature, T_dust. The dust emissivity index, β, is fixed to 1.5 and 2.0 to test different scenarios. In this way, a value of T_dust is associated with each DH02 template. This is the reference dust temperature.

Each template is then shifted to increasing redshifts, in the range z = 0−3, and convolved again with filters. The z> 0 synthetic photometry is then refit with a MBB, taking care to only use those bands with a rest-frame wavelength λ_rest> 50μm.

For each template, characterized by the value of the parameters α (Dale & Helou 2002), Figure 8 plots the values of T_dust obtained with β = 1.5 or 2.0 at z = 0 (black squares and open circles) and those derived at z> 0 with β = 1.5 (light blue grid). The inset includes a zoom on a single template, which is aimed at showing the details of the derivation of T_dust on redshift. Each redshift point was artificially shifted by a small amount along the x-axis to avoid overlapping.

For reference, the results of applying Wien’s displacements law to DH02 templates, thus deriving dust temperatures from the wavelength of the SED’s peak (Casey et al. 2012), are shown (black crosses).

Using different values of β produces a systematic change in the inferred value of T_dust, a smaller β implying a larger T_dust. At higher redshift, T_dust further increases, but the change in temperature does not follow a linear trend, showing a saw-toothed pattern instead. This reflects the fact that as redshift increases, the rest-frame wavelength sampled by each photometric band decreases, and the fitting procedure tends to be biased toward higher dust temperatures because the SED is constrained at shorter wavelengths. When one band falls shorter than 50 μm rest frame, it is discarded and T_dust shows a jump toward colder values because suddenly a short-λ band is missing. This effect happens periodically, every time one band is discarded as a result of the λ_rest> 50μm requirement.

In other words, for sources lying at different redshift, the discreteness of the sampled wavelengths and their shift to the rest frame cause a few K systematic difference in the determination of T_dust.

6. Rest-frame simulations

To understand the results presented in Sect. 5, we ran an additional simulation, which is limited to rest-frame SEDs, aimed at studying the effects of losing one or more photometric bands on the short- or long-wavelength side of the SED. At the same time, we also study the effects of having different amounts of noise in the data by randomly assigning S/N values to each band, independently. This means that, for each entry in the synthetic catalog, the distribution of relative photometric uncertainties among bands is neither flat nor based on the ratio of noise levels in selected Herschel surveys, but is a random combination.

The new simulation is structured as follows:

1. First of all, a library of models (DL07 or MBB) is generated based on a grid of input parameters (q_PAH, U_min, and γ in one case, and T_dust in the other). We limit the analysis for DL07 models to the setup suggested by Draine et al. (2007). In the MBB case, the parameter β is fixed to a value of 1.5 and T_dust spans the range 10−50 K. All models are renormalized to a fixed total dust mass of 10⁸[M_⊙] and then convolved with mid- and FIR photometric bands. We use ten bands for the DL07 models: Spitzer IRAC 8 μm, IRS 16 μm, MIPS 24 μm, Herschel PACS 70, 100, 160 μm, SPIRE 250, 350, 500 μm, and SCUBA 850 μm. We also performed an additional run, limited to six bands between 8 and 250 μm. Only eight bands are used for MBB models: a square box filter centered at 40 μm and with a half width of 10 μm, the 70, 100, 160, 250, 350, 500 μm PACS bands, and the SCUBA 850 μm passband. All entries in the catalog are at redshift z = 0. This is the so-called input catalog.

2. The input catalog is degraded in two ways:

   • An increasing number of bands is removed, progressively one by one, until only two are left. This procedure is repeated twice: first photometric bands are removed from the long-wavelength side of SEDs⁵; then the number of bands is reset to eight (i.e., the full SED is re-installed) and bands are removed one by one from the short-wavelength side of SEDs⁶.

   • The uncertainty on the photometry is increased. Relative uncertainties are increased randomly and independently for each band, spanning the range between 5% and 50% in relative errors.

In practice, the degraded catalog contains all entries of the input catalog, each one modified eight times getting rid of bands on the short-λ side and eight times getting rid of bands on the long-λ side. Each of the N × 16 entries is modified ten times with random S/N values; each time, a random level of noise is assigned to each band, independently from the other bands.

By removing bands from the long-wavelength side of the SED, we are simulating the case of deep surveys, such as the deepest fields in PEP/HerMES (Lutz et al. 2011; Oliver et al. 2012), benefiting from deep PACS photometry for the majority of sources, and progressively missing SPIRE detections because of depth, confusion, and blending effects. Vice versa, when removing bands from the short-wavelength side, we are simulating the case of shallow surveys, such as H-ATLAS (Eales et al. 2010): SPIRE quickly reaches the relatively high confusion limit at these wavelengths, but a fast PACS observation cannot fully take advantage of the smaller beam and deeper confusion limit, and remains limited to bright and mostly lower redshift sources.

6.1. Reading key

The results of the rest-frame simulations are presented with two main flavors of diagrams.

First the relative uncertainty on parameters, as obtained with Monte Carlo runs, is shown. In these diagrams, each displayed dot belongs to a full MC run per one artificial object, i.e., the value and error associated with each dot are computed as the average and standard deviation of 1000 evaluations of the same entry of the synthetic catalog. Color coding is based on the average relative photometric uncertainty of all available bands, if not otherwise specified. Darker colors refer to smaller average photometric uncertainties. The x-axis shows the maximum or minimum wavelength covered by the data, λ_max (or λ_min), depending on whether photometric bands were removed at the long- or short-wavelength side. In correspondence of the wavelength of each photometric band, a column of dots is plotted. These dots represent all those entries in the artificial catalog with that value λ_max (or λ_min), therefore, the given object in the given column benefits from available photometry in all bands shortward (longward) of the λ_max (λ_min) it belongs to. This is exemplified by the left arrows (right arrows) at the bottom of each diagram.

The second type of diagram (see also Appendix B.1) probes for systematic effects. This diagram presents the comparison of output derived quantities to their known input values as a function of λ_max,min. Each dot again belongs to a given object in the synthetic catalogs, i.e., the dot represents the average of 1000 MC evaluations.

6.2. Results of DL07: Relative uncertainties

Figure 9 presents the relative uncertainty on the free DL07 parameters and M_dust as a function of λ_max,min for cases obtained removing long-wavelength and short-wavelength bands. We limit the diagrams describing relative uncertainties to specific values of the input parameters to avoid the piling up of too many cases.

Trivially, objects with larger photometric uncertainties have larger relative uncertainties on the derived quantities. However, since the photometric S/N of each band is independent from the others and the color coding based on an average photometric relative uncertainty is computed over all bands, in Fig. 9 there is not a smooth gradient of colors, but a mixing of cases is present.

Unless the rest-frame mid-infrared SED is sampled by the available photometry (8−24 μm in the specific case), and with good S/N, constraints on the value q_PAH and γ are weak. As long as the peak of the FIR SED is sampled, and if the average S/N is good (darker dots and lower envelope of the relative error distribution in Fig. 9), it is possible to have an estimate of U_min within a ~30% uncertainty. A broader peak (induced by larger input values of U_min) is more difficult to constrain and the consequent uncertainty on U_min is larger.

The large uncertainties on the γ and U_min parameters have strong consequences on the determination of ⟨ U ⟩. Unless the SED is robustly constrained (see below), the propagation of γ and U_min errors onto Eq. (20) significantly hinders the estimate of ⟨ U ⟩. On the other hand, ⟨ U ⟩ can also be derived from L_dust and M_dust (see Eq. (19)) under the assumption that the value of P₀ is known (e.g., Magdis et al. 2012) and that the two quantities can be measured independently enough.

The best constrained and most stable parameter turns out to be dust mass, M_dust. When several bands are available and the λ> 250μm SED is sampled, dust mass is constrained within a 30% uncertainty even with poor photometry, but as soon as sub-mm data are missing, the maximal uncertainty (poor S/N, light dots) on M_dust immediately reaches 70−80%, even with the full 8−250μm wavelength range covered. When removing short-wavelength bands, the uncertainty on M_dust remains more stable than in the MBB case (see Sect. 6.4) on both the high and low S/N sides. This might be because of the choices made on the basis of SINGS results (Draine et al. 2007), effectively limiting the freedom of choosing extreme models that would cause drastic changes of normalization. In the best case of high quality photometry, the uncertainty on M_dust can be as low as 10−20% even with only three FIR bands (rest-frame 100, 160, 250 μm) available.

Upper limits are simulated with relative photometric uncertainties larger than 0.33 (see Fig. 9). If we compare the lightest blue dots in a given column to those in the subsequent column that has the given band fully removed, we find that the use of upper limits can produce an improvement in the determination of M_dust, but this improvement is nevertheless modest.

Expectations for a typical PEP galaxy at z = 1.5−2.0, and observed photometry covering the rest-frame 8−160 μm (i.e., λ_obs = 24−350 μm) can be read from the results of simulations carried out with six photometric bands and removing long-wavelength data (Fig. 9, right). The column of dots at λ = 160μm represents results for sources with rest-frame 8−160 μm photometry available. In such cases, the relative uncertainty on M_dust can be as low as ~10−20% for good S/N data (darker dots and lower envelope of the dots distribution), the exact value depending on the value of input parameters. Nevertheless, it easily grows beyond 50% for 3σ-only detections. Results only get worse by a factor ~2 for λ_rest up to 100 μm.

6.3. Results of DL07: Systematics

Possible systematic offsets are studied by comparing output and input values. We focus on dust mass estimates as the other parameters are much more poorly constrained (see Sect. 6.2). Moreover, only a short compendium of results is presented here, while a thorough description of details is given in Appendix B.1.

When removing short wavelength bands, M_dust is usually easily retrieved and there are no very significant trends related to γ_in in over- or underestimating M_dust. This also holds in varying q_PAH,in and U_min,in. When dealing with a limited number of bands, i.e., with MC runs using only six filters in total, systematics as function of U_min,in can be triggered. This is because not only bands at the short wavelength side are progressively missing, but also the long wavelength SED is sampled only up to 250 μm rest frame.

When removing long-wavelength photometric bands, systematics show trends as a function of the position in the q_PAH,in, U_min,in, and γ_in parameter space. At low U_min,in there is a tendency to underestimate M_dust if the band coverage is poor, while at high U_min,in there is a tendency to overestimate it. A larger value of γ_in produces an increased chance to overestimate M_dust. Similarly, the larger q_PAH,in, the more M_dust can be overestimated.

The three effects can add up or compensate each other. If there is a general tendency to under-estimate M_dust (e.g., because U_min,in is small), then the larger values of γ_in or q_PAH,in mitigate it. In contrast, large values of U_min,in combined to larger γ_in or q_PAH,in can induce a tendency to significantly overestimate M_dust when the SED coverage is poor at long wavelengths.

Generally speaking, if the maximum covered rest-frame wavelength is λ_max< 160−200μm systematics on M_dust can become significant, especially when the S/N of the available FIR photometry is poor.

6.4. Results of MBB: Relative uncertainties

We now analyze the case of MBB models with fixed β and how relative errors on T_dust and M_dust depend on the available photometry and on its S/N. Figure 10 shows the behavior of relative uncertainties on T_dust and M_dust, for all entries in the artificial catalog built on MBB models without distinction of input dust temperature.

The left panel shows that, as the long-wavelength part of the SED is progressively less sampled, the uncertainty on T_dust and M_dust increases, as expected. The parameter T_dust is always constrained relatively well within a 25% uncertainty. As the average photometric uncertainty increases (brighter color tones), naturally the uncertainty on the derived quantities increases. This is particularly true for M_dust. When rest-frame sub-mm data are available, with the exception of a few outliers dust masses can generally be retrieved with relative uncertainties, within 20−30% in case of poor photometry. As the maximum covered wavelength decreases, the uncertainty on M_dust explodes when the available bands do not sample the SED at λ ≥ 350 μm (restframe) anymore and the photometry is poor. This turns out to be slightly worse than DL07 modeling, which produces more stable M_dust uncertainties even with SEDs only limited to ~160 μm rest frame (see Sect. 6.2).

The right-hand panel shows instead what happens when short-wavelength bands are progressively removed. The uncertainty on T_dust is now larger than in the previous case, reflecting the fact that while the long-wavelength side of the MBB has a weak dependence on T_dust, there is a stronger dependence on the short-wavelength side. Uncertainties on M_dust are not significantly enhanced with respect to the previous case.

As the temperature increases the minimum uncertainty on T_dust tends to become larger, reflecting the fact that the peak of the MBB moves to shorter wavelengths. A similar trend is also detected for M_dust, although with smaller amplitude.

The case of a PEP z = 1.5−2.0 galaxy can be studied in Fig. 10, focusing on the left panel, at λ = 40−160μm. The observed photometry (e.g., at 100, 160, 250, 350 μm) samples the rest-frame 40−200 μm wavelength range. The simulation shows that it is possible to constrain T_dust within a <20% accuracy, but the uncertainty on M_dust could easily grow up to 80% with poor S/N. In the best case, the smallest uncertainty possible (~15%) is reached if photometric accuracies of ~15−20% are available in almost all bands (350 μm included).

6.5. Results of MBB: Systematics

Here the possible systematics on dust mass are tested for the MBB modeling along with T_dust reliability. As in the DL07 case, we only summarize the main results, while we defer to Appendix B.2 for fine details.

When the photometric accuracy is poor, M_dust can be systematically overestimated if bands are missing on the long-wavelength side of the SED; the problem is relatively milder when removing bands on the short-wavelength side.

The estimate of T_dust is not very sensitive to the lack of long wavelength datapoints. On the contrary, it can be systematically underestimated or overestimated when removing data on the blue side of the SED. This happens because the peak of the MBB emission shifts to longer wavelengths as dust temperature decreases, and the effect of poor SED sampling (due to the removed datapoints) is thus amplified. Overall, only ~10% of cases turn out to have M_dust,out> 2 × M_dust,in.

7. Discussion

In the previous sections, we built expectations on DL07 and MBB SED fitting applied to Herschel-detected sources. It is now time to derive dust and gas masses for real sources and stacked data (see Sect. 4) and to study their actual uncertainties and the properties of Herschel galaxies in this context.

7.1. The real world and its limitations

We apply the methods described in Sect. 3 to the GOODS-N and GOODS-S data, taking care to limit the MBB fit to λ_rest ≥ 50μm to avoid contamination from warm dust. We adopted a value of β = 2.08 to simplify the comparison of MBB results to DL07 dust masses, and we used the Draine (2003) revisitation of κ_ν (see also Li & Draine 2001).

The errors analysis is only focused on M_dust. In fact, simulations have already shown that via MBB fitting, T_dust can be determined within a 30% uncertainty as long as the short-wavelength side of the SED is constrained; the values of the γ, q_PAH, U_min, and DL07 parameters can hardly be constrained by the available photometry. These results are also shared in the SED fitting of real SEDs discussed here and therefore are not covered further in this section.

The left panel of Fig. 11 shows the relative uncertainty on M_dust for a MBB fit to all 160 μm-selected GOODS-N/S sources as a function of the maximum available rest-frame wavelength, λ_max, of each object. The right panel of Fig. 11 deals with DL07 models. In both cases, color coding is based on the maximum S/N of the available photometry for each object, with darker symbols indicating a higher S/N.

At larger rest-frame λ_max, the number of available bands is also larger, mostly because of redshift effects. In the case of the MBB fit, we emphasize that sources with smaller $\mathit{\sigma }{}^{\mathrm{(}}\mathit{M}_{\mathrm{dust}}{}^{\mathrm{)}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ tend to benefit from a higher quality photometry, i.e., a higher (maximum) photometric S/N at λ_obs ≥ 100μm. On the other hand, no such trend is seen in the DL07 case.

The relative uncertainty on M_dust increases as the maximum covered rest-frame wavelength decreases, thus confirming the findings of Monte Carlo simulations (see Sect. 6). In the case of GOODS fields, the use of photometric upper limits on the long-wavelength side does not provide a significant advantage in SED fitting because SPIRE data are much shallower than PACS data (up to a factor of ~10; see Table 2).

The performance turns out to be better in the DL07 case: The trend of $\mathit{\sigma }{}^{\mathrm{(}}\mathit{M}_{\mathrm{dust}}{}^{\mathrm{)}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ vs. λ_max is flatter, and the distribution is characterized by a smaller scatter at a given value of λ_max,rest. For example, at a rest-frame λ_max = 200 μm with DL07 modeling, it is always possible to constrain M_dust to a ~40% relative uncertainty, while the relative error on the MBB-based M_dust can reach values as large as 70%.

In the best case scenario of five to six bands available (covering up to 500 μm in the observed frame), the uncertainty on dust mass is typically on the order of 20% but can reach up to 30−40%, depending on the S/N of the available photometry and on the adopted modeling.

7.2. Comparison of DL07 and MBB M_dust estimates

Figure 12 reports on the comparison of M_dust estimated with MBB and DL07 SED fitting. A median systematic offset of ≤50% between the two estimates is found. Studying objects with fully sampled SEDs (up to sub-mm wavelengths), Magdis et al. (2012, 2013) report a systematic offset of a factor ~2. Magnelli et al. (2012a) report a factor ~3 discrepancy. Comparing these different results is not as straightforward as it might seem. In fact, the underlying assumptions for the MBB modeling differ significantly. We have already mentioned in Sect. 5.2.1 there is a discrepancy between the β = 1.5 emissivity adopted by Magnelli et al. (2012b) and the β used by the DH02 models, which depends on the intensity of the radiation field. Finally, the κ_ν (Li & Draine 2001; Draine 2003) adopted by most authors has a frequency dependency on the power of β = 2.08. Bianchi (2013) shows that when fitting local NGC galaxies and using β = 2.08, any offset should be washed away (see Sect. 3.2). Although we are using a consistent (β,κ_ν) set for the two approaches, this incongruence with respect to Bianchi (2013) results might come from possible mismatched temperatures due to less well-sampled SEDs of high-z objects. It is therefore always important to report on the adopted setup when referring to this delicate comparison.

Color coding the ratio of the two estimates by redshift, we see that the scatter in the distribution of points becomes very large above z ~ 1. At this redshift, one begins to lose PACS bands, on the blue side of the SED, because of k correction and sensitivity issues; at the same time, the SPIRE photometry also becomes poorer and more affected by confusion noise.

Similar results are obtained fitting the stacked photometry by Magnelli et al. (2014), which has the advantage of a more extensive SED coverage on the long-wavelength side (Fig. 12, right). In this case, it is seen that objects with longer λ_rest,max lie closer to the 1:1 locus. On the contrary, for poorly sampled SEDs, the discrepancy becomes larger (see also Sect. 6). We conclude that two concomitant effects contribute to the difference in MBB/DL07 mass ratios obtained by different authors: the underlying MBB model assumptions and the available spectral coverage.

7.3. Dust-based M_gas

Following Sect. 2.2, the DL07-based dust masses are converted into the molecular gas content of galaxies following δ_GDR-Z scaling of Magdis et al. (2012). Metallicities are computed from stellar masses with the parameterization of the M^∗-Z relation by Genzel et al. (2015). Both of these relations are calibrated to the PP04 metallicity scale.

Figure 13 compares different M_gas estimates for the datasets in hand, all inclusive of the helium contribution. CO-based masses have been rescaled to a common Milky Way conversion factor α_CO,MW = 4.36, and we have taken a metallicity correction into account as well. The adopted correction is the geometric mean of the Bolatto et al. (2013) and Genzel et al. (2012) dependencies of α_CO on metallicity (see Genzel et al. 2015).

Gas masses based on the scaling of depletion times have been obtained adopting Eq. (3) in Sect. 2.1 (Genzel et al. 2015, and priv. comm.) and the Whitaker et al. (2014) definition of star formation main sequence.

In the case of CO-detected galaxies, the CO-based and τ_dep-based M_gas estimates are consistent overall within a factor of ~2 with only few exceptions. For the low-mass lensed galaxies, Saintonge et al. (2013) adopted the metallicity-dependent values of α_CO derived using the Genzel et al. (2012) relation, nevertheless, the metallicities of most of their sources are out of the range where the Genzel et al. (2015)τ_dep scaling was calibrated (and holds). It is thus no surprise that these sources show an offset between CO- and τ_dep-based determinations of M_gas (bottom right panel of Fig. 13).

As far as CO- and dust-based estimates are compared (top-right panel of Fig. 13), for roughly one third of the sources the two M_gas estimates differ by more than 0.3 dex with dust systematically providing a higher M_gas than CO. The SEDs of these sources are missing some long-wavelength bands and, therefore, the M_dust estimate is affected by systematics, as seen in previous sections.

Dust-based M_gas of the stacked photometry by Magnelli et al. (2014) are in good agreement with τ_dep-based estimates with a relatively small scatter of ~0.23 dex, which is on the same order of the scatter in the adopted τ_dep⁽z,M^∗,sSFR⁾ scaling relation (Genzel et al. 2015). The situation for individual objects is more complex: at the low-mass end, the two estimates of M_gas are in fair agreement, but, at high masses, the dust-based estimate is significantly larger than τ_dep results. The mismatch between the two M_gas determinations becomes larger as the number of available bands decreases, mainly driven by limited wavelength coverage. This holds both for PACS-selected sources and for CO-detected objects.

For the general GOODS-N/S FIR population, λ_max(rest) ≤ 200μm and can be as low as ≤100 μm even at intermediate-low redshift because of the large noise in SPIRE bands (see Table 2). According to our simulations, in this case dust masses can be overestimated (see Sect. 6.3), thus explaining the difference in M_gas in these cases. At higher redshift, λ_max(rest) becomes even shorter if no sub-mm detections are available, and the overestimate of masses becomes more critical (see Sect. 5.2). Figure 15 limits the results on GOODS-N/S sources to λ_max(rest) ≥ 160μm.

When an object benefits from sub-mm observations (e.g., as is the case for most CO-detected galaxies and SMGs), then the maximum rest-frame wavelength available is on average longer than for other PACS-selected galaxies. If a sufficient number of bands is available, and the SED coverage is fine enough, then τ_dep- and DL07-based masses turn out to be in good agreement. Outliers with M_gas(DL07) ≫ M_gas(τ_dep) suffer from as poor SED coverage as analogous GOODS-N/S cases.

It is worthwhile recalling that the adopted τ_dep scaling was calibrated for MS galaxies of nearly solar metallicity (Genzel et al. 2015). Therefore, we also expect a contribution to the M_gas(DL07) /M_gas(τ_dep) mismatch from the adopted τ_dep scaling. Figure 14 exemplifies the possible trends of this ratio as a function of distance from the MS of star formation in terms of Δ(log ⁽sSFR⁾)_MS. The Whitaker et al. (2014) definition of MS has been adopted. The dependence of M_gas on the distance from the star-forming MS found by Genzel et al. (2015) does not play a role in this case because it is factorized out by taking the M_gas,DL07/M_gas,τdep ratio.

In the case of the Magnelli et al. (2014) stacked photometry, the data show a trend of the τ_dep-based M_gas estimates of bins above the MS to be larger than those based on dust. Limiting to Δ(log ⁽sSFR⁾)_MS ≤ ± 0.5 dex (vertical dashed lines), we restrict the analysis to the MS proper locus. The scatter of the M_gas(DL07) /M_gas(τ_dep) decreases to 0.19 dex for the stacked photometry. As a result of stacking, the SED is covered up to λ_max(rest) ~ 160μm and up to z = 2. Our simulations (Sect. 6.3) prove that this is enough to avoid underestimates of M_dust and keep overestimates to less than a factor of 2 (details depending on the actual values of q_pah, U_min, and γ).

As for individual sources, the majority of catastrophic failures lying at > 3σ (measured above using the result for stacking) above M_gas(DL07) /M_gas(τ_dep) = 1 turn out to be on the high-z tail of the Herschel sample. The SEDs of these galaxies are poorly sampled and are only covered up to λ_max(rest) ≤ 100μm.

In summary, the analysis of the SEDs of real sources and stacked photometry by means of DL07 fitting confirms the expectations from the MC exercise discussed before. When the mid- and FIR rest-frame SEDs are covered up to at least 160 μm, the DL07 dust-based estimate of M_gas is reliable, not affected by systematics, and consistent with independent estimates based on τ_dep scalings or CO observations. If the SED wavelength coverage is poorer, i.e., limited to shorter wavelengths, dust masses can be overestimated. This might happen because the SED is redshifted and therefore the available bands sample shorter wavelengths and also because of a decreased S/N and increased confusion in the case of faint, high-z galaxies. If the maximum rest-frame wavelength available falls shorter than 100 μm, DL07-based dust masses are severely overestimated and should be not recommended anymore.

7.4. The dust-to-gas ratio

One assumes a value of the gas-to-dust mass ratio, δ_GDR to compute the molecular gas content of galaxies from their dust masses. As seen in Sect. 2.2, it is common to adopt the local scaling of δ_GDR with the metallicity derived by Leroy et al. (2011) or some a variant of it (e.g., Magdis et al. 2012). This procedure implicitly assumes that the relation holds regardless of redshift.

The sample of CO-detected sources in hand benefits from fully independent measurements of gas and dust masses. We combine CO-based M_gas with the DL07-based determinations of M_dust to derive δ_GDR. All gas masses have been renormalized to the Galactic α_CO,MW = 4.36 and a metallicity dependence of α_CO has been included as well (Genzel et al. 2015). The M^∗-Z relation described in Eq. (6) is used with stellar masses from BC03 SED fitting (Wuyts et al. 2011a) to produce metallicities on the PP04 scale.

We only consider sources with ${\mathit{\lambda }}_{\max }{}^{\mathrm{\right(}}\mathrm{rest}^{\mathrm{\left)}}\mathrm{\ge }\mathrm{200}$μm hereto minimize the uncertainties on M_dust and avoid systematics (see Sects. 5 and 6). The left panel in Fig. 16 shows the resulting trend of δ_GDR vs. $\mathrm{12}\mathrm{+}\log {}^{\mathrm{\right(}}\mathrm{O}\mathit{/}\mathrm{H}^{\mathrm{\left)}}$ compared to the Magdis et al. (2012) local relation and its 0.15 dex scatter. A ±1.0 dex systematic uncertainty on its zero point is also reported by Magdis et al. (2012).

The data lie close to the local relation and are consistent with it within the uncertainties and possible systematics. At the high metallicity end, the scatter is very large and there exist cases with very small δ_GDR, which are significantly below the locus occupied by local galaxies. The right-hand panel of Fig. 16 seeks possible dependencies of δ_GDR on other derived quantities. No significant trend is found as a function of M^∗, SFR, τ_dep, or redshift. The most critical outlier turns out to be a z ≃ 1.4 galaxy with a very high M_dust ~ 8 × 10⁹ [M_⊙] and an SED with 10−205 μm rest-frame coverage.

7.5. Synergies between Herschel and ALMA

The photometric determination of dust masses through SED fitting allows us, in principle, to obtain an estimate of gas masses for a very large number of objects via the δ_GDR scaling. This approach offers the advantage that photometric observations are still significantly faster than sub-mm spectroscopy, even with the last generation of FIR or sub-mm facilities, and especially for galaxies at z> 1.

Far-infrared photometry also has the advantage of sampling the dust emission of galaxies near the SED peak, providing a calorimetric measurement of their SFR (Elbaz et al. 2011; Nordon et al. 2010) for a relatively cheap time investment. Nevertheless, the analysis has shown the effects of the limited wavelength coverage: when the SED extends only up to λ_rest,max ≤ 200μm, the uncertainty on the derived dust mass can be very large and systematics might affect this measurement in a fraction of cases (see Sects. 7.1, 6).

On the sub-mm side, assuming the Rayleigh-Jeans (RJ) tail of SEDs is observed, Scoville et al. (2014) developed a strategy aimed at deriving gas (dust) masses on the basis of a one-band sub-mm continuum measurement. If confirmed reliable, this approach would be extremely competitive in terms of exposure time with respect to CO spectroscopy. Genzel et al. (2015) showed that this method can lead to incorrect (up to a factor >3) gas masses, even when applying corrections to take the RJ approximation into account. In fact, the proposed scaling of sub-mm luminosity holds for a specific dust temperature T_dust = 25 K, which is a condition that is only met in a few cases. Magnelli et al. (2014) have shown that T_dust varies across the MS of star formation, and increases as a function of Δ(log ⁽sSFR⁾)_MS. A dust temperature measurement is hence necessary to reach a correct estimate of M_dust from the sub-mm emission, even in the simple MBB, RJ approximation.

If not using the available T_dust⁽z,M^∗,sSFR⁾ scaling (e.g., Magnelli et al. 2014), Genzel et al. (2015) suggest a strategy based on continuum observations in two distinct sub-mm bands. In this context, using ALMA bands 6 and 7 (centered at 1100 and 850 μm), the best relative uncertainty of M_dust that can be reached is on the order of 60%, for galaxies at z ~ 2 with 10σ detections in the two sub-mm bands. The uncertainty on dust mass becomes worse at lower redshift. In practice, the SED coverage is too small and too far off the SED peak to allow for a robust determination of T_dust and M_dust. Adopting bands 7 and 9 (at 850 and 450 μm) lowers $\mathit{\sigma }{}^{\mathrm{(}}\mathit{M}_{\mathrm{dust}}{}^{\mathrm{)}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ to 30−40% at best and using bands 6 and 9 finally brings it to 20−30%, i.e., a better than 3σ estimate of M_dust. However, band 9 is particularly demanding in terms of atmospheric conditions. At z< 1, this approach is not necessarily faster than ALMA CO spectroscopy. We defer to Genzel et al. (2015) for further details.

It is interesting to study how FIR and sub-mm observations complement each other and how combining them improves M_dust measurements. We simulate the ALMA performance in conjunction with Herschel data in the following way. For each galaxy detected by PACS in the GOODS-S field, we assume that the best-fit DL07 model obtained in Sect. 7.1 represents the actual emission of the galaxy. This best-fit model is then convolved with a box ALMA passband, centered at ~1100 μm (i.e., in ALMA’s band 6). The artificial photometry thus obtained is added to the real Spitzer and Herschel data and the extended SED is fitted again following the usual procedure. We assume that a S/N of 5 holds in ALMA band 6.

Figure 17 compares the uncertainties on M_dust obtained with and without the artificial band 6 photometric point as a function of the maximum rest-frame wavelength covered by real data, λ_rest,max, and dust mass (see also Fig. 11). The availability of good quality (S/N = 5 in this case) photometry at λ_obs = 1100 μm reduces the uncertainty on M_dust to ≤33% for >85% of sources. For comparison, only ~20% of the sample have similarly good quality M_dust estimates when no ALMA data are available. Similar results are obtained for the Magnelli et al. (2014) stacked data, adding the ALMA band 6 filter to the simulation described in Sect. 5.

We then use the artificial ALMA band 6 flux density obtained above to derive a gas mass expectation applying the recipe by Scoville et al. (2015; see their Appendix A and also Scoville et al. 2014). As usual, DL07-based dust masses are transformed into M_gas adopting the Genzel et al. (2015)M^∗-Z-z and the Magdis et al. (2012)δ_GDR-Z relations (see Sect. 2.2). The M_gas expectation based on band 6 flux only (Scoville et al. 2015) includes the RJ correction prescribed by the authors. Figure 18 compares the two.

The DL07-based M_gas estimate tends to be systematically larger than the expectation based on the recipe by Scoville et al. (2015). For a small fraction of the sample the opposite trend holds. The overestimate becomes larger at the high-mass end. The trend seen here is the combination of the systematic higher M_dust obtained with DL07 modeling with respect to MBB (see Sect. 7.2), the fact that dust temperature is varying as a function of the position in the M^∗-SFR-z space (see Genzel et al. 2015; Magnelli et al. 2014), a slight difference in the β values underlying the two methods (β = 2.08 vs. 1.8), and possibly other systematics (e.g., related to metallicity dependencies). This is different from the problem encountered in Fig. 13 at large M_gas because now we are reasoning in relative terms (DL07 vs. Scoville methods), and the long wavelength side of the SED is constrained by the artificial ALMA photometry.

Following a similar path to Genzel et al. (2015), we now estimate the on source exposure time, t_exp, needed to reach S/N = 5 in band 6 for our Herschel-detected star-forming galaxies based on DL07 models.

To this aim, the ALMA sensitivity calculator⁷ is used with the standard Cycle 3 configuration via an array of 36 12 m antennas, requesting an angular resolution of 1 arcsec. The left-hand panels of Fig. 19 show t_exp for individual sources and the cumulative t_exp for sources below a given dust mass (solid histograms) along with the number distribution of objects (dotted histograms). The dispersion of the points in the upper panel reflects the spread in redshift of the sample. However, the sub-mm negative k correction implies that exposure times for a given mass are of the same order of magnitude at all redshifts.

Targeting all Herschel galaxies detected by the PEP/GOODS-H survey in the GOODS-S field down to M_dust = 10⁹ [M_⊙] would require a few minutes of on-source exposure, i.e., without accounting for overheads. Reaching down to M_dust = 10^8.5 [M_⊙] requires roughly 20−25 h on source with ALMA band 6 in the above mentioned configuration. Using ALMA band 7 would need a similar amount of time.

The right-hand panel of Fig. 19 includes the distribution of the expected ALMA band 6 fluxes as a function of dust mass. Instead of requesting a fixed S/N, a 0.38 mJy and 0.18 mJy depths are shown, similar to recent ALMA Cycle-3 approved surveys⁸. At these depths, ALMA detects the majority of Herschel sources above M_dust = 10^8.5 [M_⊙] with a fainter tail extending down to ~10⁸ [M_⊙], and with the added value of possible undetected PACS objects, which were not included in the current analysis.

Assuming an average dust/gas ratio of 0.01, the limits mentioned above correspond to M_gas = 10^10.5−11 [M_⊙]. Berta et al. (2013a) showed that the Schechter characteristic mass, ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{\ast }}$ of the molecular gas mass function lies between M_gas = 10^10.3−10.9 [M_⊙] at z = 0.2−2.0. A measurement of the molecular gas mass function based on dust observations, modulo the δ_GDR scaling, and on the dust mass function itself down to ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{\ast }}$ and up to z ~ 2, is within reach of ALMA in 20−25 h spent on source.

8. Conclusions

We have exploited the deepest FIR blank field maps available to date from the PEP, GOODS-H, and HerMES surveys to study the feasibility of deriving dust and gas masses via SED fitting for individually detected sources. In parallel, we built extensive Monte Carlo simulations to study the limitations of real data, and to understand how they influence the uncertainties and systematics on such dust mass determinations. We focused the analysis on two popular modeling approaches: SED fitting by means of a single-temperature modified blackbody and by means of the Draine & Li (2007) model. The main results of the analysis based on artificial sources and MC sampling are:

• FIR SED fitting recovers dust mass consistently as long as the wavelength coverage offered by the data extends at least up to 160−200 μm (rest frame). In this case, no systematics are expected on M_dust. In contrast, if the >3σ detections fall shorter than this wavelength limit, M_dust can be severely overestimated with the amplitude of the systematic effect depending on the details of the model.

• The uncertainty of M_dust also strongly depends on wavelength coverage. As a rule of thumb, it is not possible to reach a ≥ 3σ determination of M_dust if the available photometry does not extend at least up to ~200 μm in the rest frame.

• The determination of dust temperatures based on MBB fitting is rather stable, even if long-wavelength bands are missing. In addition, T_dust is always constrained to better than a 10−20% uncertainty as long as the blue side of the FIR SED is constrained. Nevertheless, a small offset on T_dust can induce a large systematic error on M_dust.

• The discreteness of the sampled wavelengths (photometric bands), combined with redshift, can cause significant residual systematics in the determination of T_dust, up to few K.

• Reliability tests based on artificial catalogs tend to minimize systematics on M_dust if the SED shape adopted for SED fitting is consistent with the one used to produce the mock catalog. In other words, these kinds of simulations assume that the adopted SED library is consistent with real-world’ SEDs. On the other hand, small differences in SED shape (e.g., DH02 vs. DL07 models) can induce strong and not easily predictable systematic effects.

The GOODS fields benefit from an extensive multiwavelength coverage from the UV to FIR and sub-mm wavelengths. The SEDs of galaxies detected by Herschel can reach up to 200 μm rest frame and beyond. The SEDs of individual sources have been fitted with DL07 and MBB models. At the same time, stacked photometry of NIR selected galaxies has also undergone the same analysis, binned in z-M^∗-SFR space. The main results of these pieces of analysis are:

• At the depth of the deepest Herschel extragalactic surveys (GOODS-S as observed by PEP, GOODS-H, and HerMES), it is possible to retrieve dust masses with a S/N ≥ 3 for galaxies on the main sequence of star formation down to M^∗ ~ 10¹⁰ [M_⊙] up to z ~ 1. At higher redshift (z ≤ 2), the same goal is achieved for objects only lying increasingly above the MS at similar stellar masses or for galaxies at the tip of the MS (i.e., with higher M^∗). At shallower depths (e.g., in the case of the COSMOS field), this reasoning shifts to even higher values of SFR. Dust temperatures (based on MBB fit) can be constrained within a 10% accuracy in most of the cases across the z-M^∗-SFR space, modulo residual systematics due to the discreteness of SED sampling (see above).

• As in the case of simulated data, Spitzer and Herschel data alone are not sufficient to produce an estimate of M_dust to better than 30% uncertainty if the maximum rest-frame wavelength covered by the data is shorter than ~160−200 μm.

• Comparing MBB- and DL07-based masses, the average offset between the two, regardless of redshift, is a factor ~1.5. At z> 1, the scatter of the MBB/DL07 mass ratio becomes very large, mainly because photometric points are progressively missed. We stress that to allow a direct and meaningful comparison, it is paramount to adopt a consistent set of parameters, taking special care in the value of β, i.e., the power of the dependence of dust emissivity on frequency.

• Dust masses estimated with DL07 modeling are more robust than those based on MBB: relative errors are more mildly dependent on the maximum covered rest-frame wavelength and less scattered.

• Dust mass estimates, based on DL07 SED modeling and on scaling of depletion times, τ_dep, are consistent with each other as long as the data guarantee a sufficient wavelength coverage. Applying the local dependence of δ_GDR on metallicity to transform M_dust into M_gas, these estimates are overall consistent with CO-based estimates for a small sample of star-forming galaxies. While comparing M_gas estimates based on different methods, it is important to adopt a consistent set of relations, calibrated to a common metallicity scale.

• Using CO-based M_gas renormalized to α_CO,MW = 4.36, M_dust obtained through DL07 modeling, and metallicities computed with the M^∗-Z-z relation (Genzel et al. 2015), the δ_GDR of z> 1 galaxies depends on metallicity in a similar manner as for local galaxies within uncertainties and systematics.

CO-based M_gas estimates, which also represent an important anchor for validation of indirect methods despite α_CO uncertainties, are still limited at z> 1 (see also Genzel et al. 2015). More CO observations of individual galaxies covering a wide range of parameters are highly desirable and are becoming more easily accessible with ALMA and NOEMA. In parallel, also the determination of gas-phase abundances and, therefore, of the redshift evolution of ISM physical properties (e.g., Kewley et al. 2013; Shapley et al. 2015), will undergo significant progress over the next few years thanks to NIR multiobject spectroscopy (e.g., with KMOS, MOSFIRE).

Finally, as M_dust estimates based on Spitzer and Herschel photometry are limited to cases with high quality SEDs available, we recalled the advantages and limitations of estimates including only sub-mm data (e.g., Genzel et al. 2015). For example, a scaling of sub-mm fluxes into M_dust can be affected by strong systematics if the characteristic dust temperature of the SED is not known. Continuum observations in two sub-mm continuum bands might eventually help, but they are time consuming and the results are still limited by large uncertainties because the wavelength range covered is relatively small.

We therefore explored a combined IR plus sub-mm approach, combining existing Herschel data to expected ALMA 850 or 1100 μm continuum fluxes. These single band observations allow one to reduce the uncertainties on M_dust down to <30% for virtually all Herschel-detected galaxies in the GOODS-S field. A direct measurement of the molecular gas mass function based on dust observations up to z ~ 2 will be soon within reach.

---

Acknowledgments

The authors would like to thank the referee, Dr. Georgios Magdis, for his useful and effective comments that improved the quality of this manuscript, as well as for his positive and fast response. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada), NAOC (China), CEA, LAM (France), IFSI, University of Padua (Italy), IAC (Spain), Stockholm Observatory (Sweden), Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK), Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK) and NASA (USA).

References

Appendix A: Ancillary data for CO samples

The PHIBSS survey (Tacconi et al. 2013) performed CO spectroscopy of galaxies at z = 1−2, and included additional data from past work, for a total of 73 objects. The survey observed 38 sources in the Extended Groth Strip (EGS) field, also targeted by PEP and HerMES. Another 18 objects belong to the “Q-fields” (Steidel et al. 2004), Tacconi et al. (2013) also included six BzK sources by Daddi et al. (2010) and six sources by Magnelli et al. (2012b), plus a few other additional objects.

FIR photometry is sought in the EGS from the PEP and HerMES data (Lutz et al. 2011; Oliver et al. 2012), using a closest-neighbor algorithm and visual inspection of multiwavelength maps (IRAC, MIPS 24 μm, PACS, and SPIRE). The sources by Magnelli and Daddi, and few other isolated objects studied by Magnelli et al. (2012b) and Saintonge et al. (2013), and their photometry can be found in their works.

Out of the 38 PHIBSS sources in EGS, we have the following detection statistics: 37, 18, 15, 14, 15, and 13 are detected at 24, 100, 160, 250, 350, and 500 μm, respectively. Five out of the six BzK sources by Daddi et al. (2010) included in PHIBSS have 24 μm to 500 μm photometry by PEP and HerMES (Lutz et al. 2011; Oliver et al. 2012). The object BzK21000 benefits from millimeter photometry (1.3, 2.2, 3.3 mm) by Dannerbauer et al. (quoted as in prep. by Magdis et al. 2011), Carilli et al. (2010), and Daddi et al. (2009) as collected by Magdis et al. (2011). CB58 and the Cosmic Eye have photometry collected by Saintonge et al. (2013). For the other sources in the PHIBSS sample (mainly the Q-fields), no FIR photometry has been retrieved from the literature yet.

Magnelli et al. (2012b) included in their analysis six PEP sources, two HDF sources, and the well-studied object GN20. They actually included nine sources in their observations, but only six are used by Tacconi et al. (2013). For the six PEP objects, 24-to-500 μm photometry has been retrieved from PEP and HerMES catalogs (Lutz et al. 2011; Oliver et al. 2012). The two HDF169 and HDF242 objects also have 24-to-500 μm photometry Magnelli et al. (2012b). In addition to the 24-to-500 μm imaging, GN20 also has 850, 1100, 2200, 3300, 6600 μm photometry collected by Magdis et al. (2011) and obtained by Pope et al. (2006), Perera et al. (2008), Dannerbauer et al. (2009), Daddi et al. (2009), Carilli et al. (2011). For this source, a 1160 μm photometry is also available, obtained by Penner et al. (2011) from the combination of 1100 μm Atzec and 1200 μm MAMBO maps.

Saintonge et al. (2013) have observed 17 lensed galaxies, out of which ten have CO detections (and derived physical quantities). The available photometry includes 100-to-500 μm PACS and SPIRE data, and 1200 μm IRAM/MAMBO measurements for five of them.

Bothwell et al. (2013) observed 40 SMGs in CO and detected 32 of them.

They also included some additional data taken from the literature. Their sample is spread over several different sky areas: three sources in the Subaru deep field, or UDS; three sources in smaller fields; five sources in the Lockman Hole East (LH-East); 16 sources in the HDF; three sources in SSA-13; seven sources in ELAIS-N2 (seven sources); and three sources in SSA-22. The Subaru Deep field (SXDF or UDS) is part of HerMES; the LH-East is part of PEP and HerMES; the HDF is part of PEP and HerMES (GOODS-N); ELAIS-N2 is part of HerMES, but has not been released yet (at the time of writing, in DR2 and DR3). In the UDS field, only 250, 350, 500 μm bands from HerMES DR2, plus the 850 μm by Bothwell et al. (2013) are available.

In the HDF (GOODS-N) there are the 24, 100, 160, 250, 350, 500 μm data by PEP, GOODS-Herschel and HerMES, plus the 850 μm fluxes by Bothwell. In the LH-East there are the 24, 100, 160, 250, 350, 500 μm data by PEP + HerMES, plus the 850 μm fluxes by Bothwell. None of the Bothwell et al. (2013) sources is detected in the AzTEC 1.1 mm maps by Michałowski et al. (2012).

In synthesis, out of the 40 SMGs in the sample by Bothwell et al. (2013), 16 have a secure CO detection and six have a candidate detection. Out of these, 16 detected objects lie in UDF, LH, or HDF, and two candidates lie in HDF. Out of these 16+2, 13 have enough mid- and FIR photometric detections to ensure an estimate of dust mass based on DL07 SED fitting. Out of these 13, seven sources have enough data to be included in our analysis (i.e., M_gas,CO, M_dust, M_stars).

Appendix B: Detailed analysis of systematics

We dissect the fine details of possible systematic effects in the derivation of M_dust, as emerging from the MC runs with rest-frame photometry (see Sects. 6, 6.3, and 6.5).

Appendix B.1: DL07 systematics

The DL07 case is analyzed first. We only focus on dust mass estimates because the other parameters are much more poorly constrained (see Sect. 6.2).

Figure A.1 presents the comparison of M_dust,out and M_dust,in as a function of λ_max,min, i.e., for cases obtained removing long-wavelength bands (left hand diagrams) and short-wavelength bands (right). Simulations including ten photometric bands are shown. Possible dependencies on q_PAH,in, U_min,in, and γ_in are studied by fixing two parameters and splitting the analysis in bins of the third.

First of all, we focus on the top row of diagrams, showing cases at specific values of q_PAH,in = 0.47 and U_min,in = 0.70 and sampling γ_in in the allowed range. When removing long-wavelength bands, there is a tendency to underestimate M_dust, if the bands coverage and the S/N are not adequate. This tendency becomes smaller for larger values of γ_in.

This is still a special case, and it is necessary to disentangle the effect of the other two parameters q_PAH and U_min (see middle and bottom panels in Fig A.1). For example, moving to U_min,in = 2.0, the tendency is to overestimate M_dust, rather than to underestimate it. This difference is indeed mainly driven by U_min,in (see middle panels), but the overestimate becomes larger for larger γ_in. A similar effect is also shared by q_PAH (bottom panels). A larger value of γ_in produces an increased chance to overestimate M_dust. Nevertheless, if there is a general tendency to underestimate M_dust (e.g., because U_min,in is small), then this underestimate becomes milder for the larger values of γ_in. If, instead, there is the tendency to overestimate M_dust (e.g., because of the large value of U_min,in), then this tendency becomes even worse with larger γ_in.

However, here we are speaking of large values of γ_in (>0.4), which are rarely seen in real galaxies (see, e.g., Tab. 4 in Draine et al. 2007). Keeping γ_in< 0.2, differences in trends are more difficult to appreciate.

When removing short-wavelength bands, M_dust is usually easily retrieved and there are not very significant trends related to γ_in in over- and underestimating M_dust. This holds also varying q_PAH,in and U_min,in.

We now examine the middle panels in Fig. A.1, highlighting what happens when sampling U_min,in variations at q_PAH,in = 0.47 and γ_in = 0.11. When removing long-wavelength bands, there are now big differences in the in/out trends as a function of λ_max, depending on the actual value of U_min,in. At low U_min,in there is a tendency to underestimate M_dust if the band coverage is poor, while at high U_min,in there is a tendency to overestimate it.

This recalls what is happening with the MBB (see Sect. 6.5). Changing the input U_min,in, in fact, produces a change of the position of the FIR peak (see Fig. 1 and Sect. 3.3). Consequently, there is a tendency to overestimate M_dust when the combination of bands coverage and “equivalent temperature” is produces the effects described for the MBB case study (see Sect. 6.5). At lower U_min,in, the SED peaks at longer wavelengths and at some point it exits the range where overestimate of fluxes is “preferred”, causing underestimates of M_dust to dominate.

When removing short-wavelength bands the same effects are seen, although with a much smaller amplitude.

Finally, in the bottom panels of Fig. A.1, we focus on U_min,in = 0.7 and γ_in = 0.11, and follow the variation of q_PAH,in. In this case, as in the case of variable γ_in, the trend to over- or underestimate M_dust changes as a function of the value of q_PAH,in.

When removing long-wavelength bands, for the quoted choice of U_min,in and γ_in, dust masses can be underestimated for small values of q_PAH,in, but this trend becomes much milder for large values. On the other hand, remember that for U_min,in = 2.0 and same range of γ_in, the tendency is to overestimate M_dust (see middle panels of Fig. A.1), and the tendency is to have a larger overestimate when q_PAH,in becomes larger.

In summary, the larger q_PAH,in, the more M_dust can be overestimated. If there is a tendency of underestimating M_dust (e.g., driven by a small value of U_min,in), then a large q_PAH,in makes this underestimate milder. On the other hand, if U_min,in is large and thus there is already the tendency to overestimate M_dust, then a large q_PAH,in makes things worse.

Using Eq. (20) to derive the mean radiation field ⟨ U ⟩, a similar behavior of the ratio ⟨ U ⟩ _out/ ⟨ U ⟩ _in is observed. Depending on the actual values of U_min, γ and q_PAH, the chance of overstimate or underestimate ⟨ U ⟩ can be amplified or suppressed due to summing or compensating effects of the three parameters. Since ⟨ U ⟩ and M_dust is linked by Eq. (19), this explains the relative stability of the M_dust estimate on a more fundamental level.

When removing short-wavelength bands, as usual there are very small tendencies for under- or overestimates of M_dust. In this case, by varying q_PAH,in these trends do not seem to change significantly, or at least not significantly enough to be appreciated by this analysis.

Simulations limited to a maximum wavelength of 250 μm, obviously behave equivalently when removing long-wavelength bands. On the other hand, when removing short-wave bands, things are different because the three longest wavelength bands are missing. Trends as a function of γ_in are still very small and are only slightly enhanced with respect to the ten bands case, but still much smaller than when removing long-wavelength bands. Trends as a function of U_min,in also start having a larger amplitude when removing short-wavelength bands, while trends with ten bands continued to be marginal. Finally q_PAH,in trends are marginal, as in the case of ten bands.

Appendix B.2: MBB systematics

It is worthwhile to verify how reliably the MBB fit to the artificial photometry can retrieve the input values of dust temperature and mass. Figure B.1 shows the comparison of input and output quantities, as a function of λ_max,min, split in bins of input dust temperature, T_dust,in.

As the average photometric uncertainty grows (lighter blue dots), there is the chance to systematically overestimate M_dust. The overestimate worsens when removing datapoints at the long-wavelength side, while the problem is relatively milder (but still significant) when removing bands at the short-wavelength side.

When removing datapoints at the long-wavelength side, T_dust is always well recovered, but when removing datapoints at the short-wavelength side it can be systematically underestimated or overestimated. This effect reflects the fact that the peak of the MBB shifts to longer wavelengths as dust temperature decreases, combined with the poor sampling of the SED obtained by removing datapoints.

Nevertheless, the incidence of catastrophic M_dust overestimates is small: only ~10% of cases turn out to have M_dust,out> 2 × M_dust,in. Finally, there exists a 0.1% of cases with the tendency to underestimate M_dust for the lowest T_dust,in considered, when removing datapoints at the short-wavelength.

Although affecting a limited fraction of cases, it is interesting to study in detail the nature of these systematic trends to understand what causes them. Figure B.2 directly compares the ratio of input/output T_dust to the same quantity for M_dust. Color coding is now based on the number of available bands.

As long as short-λ bands are progressively removed, it is possible to overestimate T_dust at low input temperatures (T_dust,in ≤ 20 K, for example). Consequently, the shape of the MBB SED changes, but the long-λ bands still give a good constraint on the model normalization, hence, on M_dust. Something different happens at higher input dust temperatures. At T_dust,in = 25−35 K, the output T_dust is more stable but there are still variations, which start driving an overestimate of dust mass. The peak of the MBB starts moving out of the covered wavelength range and the shape of the model in the covered range of wavelengths is roughly constant (Rayleigh-Jeans regime; RJ). Nevertheless, an underestimate of T_dust implies a lower emissivity, and thus a larger normalization (i.e., a larger dust mass) is needed to reproduce the “observed” fluxes. This causes an overestimate of dust mass.

When removing long-λ bands, the trend is again to overestimate M_dust, but the cause is more subtle. In fact, in this case, T_dust,out is rather consistent with T_dust,in. This is because the short-wavelength side of the SED is always constrained by the available photometry. There is still a possibility of slightly underestimating dust temperatures, however, because by anchoring the model at the short-wavelength side, it is still possible to reproduce the short-λ colors with a lower temperature, which allows for a different normalization (i.e., a different M_dust). The smaller the number of long-wave bands, the larger is the freedom in renormalization. However, this time a small effect on T_dust translates into a big effect on M_dust. A small variation of T_dust implies a significant change in the shape of the MBB SED blueward of the peak. Now the short-wavelength side of the SED is constrained by the data because we are removing long-λ bands. Therefore a small underestimate of T_dust induces a significant overestimate of M_dust.

In the previous case (i.e., removing short-λ bands), the opposite was happening: a big variation on T_dust implied a relatively smaller variation of M_dust because we were sampling the RJ side of the SED, where shape variations due to T_dust changes are less prominent than at the short-wavelength side. For example, focusing on input temperature in the range T_dust,in = 45−50 K one can note that:

• When removing long-wavelength bands, a variation of 10% in T_dust turns out to produce a change of normalization of up to a factor of 10.

• When removing short-wavelength bands, a variation of 30% in T_dust induces a change in normalization of up to a factor of 4−5 “only”.

All Tables

All Figures