- Split View
-
Views
-
Cite
Cite
C J R Clark, P De Vis, M Baes, S Bianchi, V Casasola, L P Cassarà, J I Davies, W Dobbels, S Lianou, I De Looze, R Evans, M Galametz, F Galliano, A P Jones, S C Madden, A V Mosenkov, S Verstocken, S Viaene, E M Xilouris, N Ysard, The first maps of κd – the dust mass absorption coefficient – in nearby galaxies, with DustPedia, Monthly Notices of the Royal Astronomical Society, Volume 489, Issue 4, November 2019, Pages 5256–5283, https://doi.org/10.1093/mnras/stz2257
- Share Icon Share
ABSTRACT
The dust mass absorption coefficient, κd is the conversion function used to infer physical dust masses from observations of dust emission. However, it is notoriously poorly constrained, and it is highly uncertain how it varies, either between or within galaxies. Here we present the results of a proof-of-concept study, using the DustPedia data for two nearby face-on spiral galaxies M 74 (NGC 628) and M 83 (NGC 5236), to create the first ever maps of κd in galaxies. We determine κd using an empirical method that exploits the fact that the dust-to-metals ratio of the interstellar medium is constrained by direct measurements of the depletion of gas-phase metals. We apply this method pixel-by-pixel within M 74 and M 83, to create maps of κd. We also demonstrate a novel method of producing metallicity maps for galaxies with irregularly sampled measurements, using the machine learning technique of Gaussian process regression. We find strong evidence for significant variation in κd. We find values of κd at 500 |$\mu$|m spanning the range 0.11–0.25 |${\rm m^{2}\, kg^{-1}}$| in M 74, and 0.15–0.80 |${\rm m^{2}\, kg^{-1}}$| in M 83. Surprisingly, we find that κd shows a distinct inverse correlation with the local density of the interstellar medium. This inverse correlation is the opposite of what is predicted by standard dust models. However, we find this relationship to be robust against a large range of changes to our method – only the adoption of unphysical or highly unusual assumptions would be able to suppress it.
1 INTRODUCTION
Interstellar dust provides an indispensable window for studying galaxies and their evolution. Dust, which primarily emits in the mid-infrared (MIR) to far-infrared (FIR) to submillimetre (submm) wavelength regime, can be observed in very large numbers of galaxies very rapidly, with the beneficial effects of negative k-correction enhancing our ability to detect dusty galaxies out to high redshift (Eales et al. 2010a; Oliver et al. 2012). This has made dust a standard proxy for studying galaxies’ star formation (Kennicutt 1998; Buat et al. 2005; Kennicutt et al. 2009), gas mass (Eales et al. 2012; Scoville et al. 2014; Lianou et al. 2016), and chemical evolution (Rowlands et al. 2014; Zhukovska 2014; De Vis et al. 2017a,b, 2019) – which are otherwise difficult and time consuming to observe directly.
However, many of the valuable insights that dust can provide rest upon one simple expectation – that we are able to use observations of dust emission to actually infer physical dust masses. Unfortunately, astronomers remain terrible at this. This is due to the fact that κd (variously called the dust mass absorption coefficient, or the dust mass opacity coefficient), the wavelength-dependent conversion factor used to calculate dust masses from FIR–submm dust spectral energy distributions (SEDs), is extremely poorly constrained.
Fig. 1 compiles a wide range of κd values that have been reported in the literature (all have been converted to a reference wavelength of 500 |$\mu$|m as per equation (2); we only plot values for which the original quoted reference wavelength was in the reliable 150–1000 |$\mu$|m range). Over 100 values are plotted, with a standard deviation of 0.8 dex, and spanning a total range of over 3.6 orders of magnitude. Worse still, there is no sign that values of κd reported in the literature are converging over time.
So, despite the excellent sensitivity and wavelength coverage provided by modern FIR–mm observatories, any dust masses inferred from observed dust emission remain enormously uncertain, stymieing our understanding of the interstellar medium (ISM) in galaxies. Moreover, this high degree of uncertainty means that, out of necessity, κd is often treated as being constant – even though it is well understood that this cannot be true in reality. Even the more complex, multiphase dust model frameworks, such as those of Jones et al. (2013, 2017), usually only incorporate two or three types of dust, each with a corresponding κd.
As such, understanding how kappa varies – both between different galaxies, and within individual galaxies – is clearly vital for the field.
In this paper, we use an empirical method for determining the value of κd – which we employ on a resolved, pixel-by-pixel basis in two nearby galaxies – to produce the first maps of how κd varies within galaxies, as a proof-of-concept study. The theory behind the dust-to-metals method we employ to find κd is described in Section 2. The galaxies and data we use in this work are described in Section 3. The application of the technique to produce maps of κd is Section 4. Our results are presented in Section 5, and are discussed in Section 6. For brevity and readability, ‘flux density’ will be termed ‘flux’ throughout the rest of the paper.
2 THEORY
Of the many methods proposed for estimating the value of κd, one of the most simple is that first proposed by James et al. (2002). The James et al. (2002) method is entirely empirical, and relies upon just one central assumption – that the dust-to-metals ratio in the ISM, ϵd, has a known value. If the ISM mass of a galaxy is known, along with the metallicity of that ISM, it is straightforward to calculate the total mass of interstellar metals in that galaxy; then, by assuming a fixed dust-to-metals ratio, it is possible to infer a galaxy’s dust mass a priori, without any reference to the dust emission. This a priori dust mass can then be compared to that galaxy’s observed dust emission, and hence κd can be calibrated. Here we use the ϵd notation for the dust-to-metals ratio, instead of |$\mathcal {D\, TM}$|. This maintains consistency with James et al. (2002) and Clark et al. (2016), and avoids any ambiguity arising from the fact that |$\mathcal {D\, TM}$| is often used to denote a dust-to-metals ratio normalized by the Milky Way value, whereas our quoted dust-to-metals ratios are always absolute values.
The vast majority of all reported values of ϵd lie in the range 0.2–0.6 (considering only values of ϵd that are not based upon some assumed value of κd: Issa, MacLaren & Wolfendale 1990; Luck & Lambert 1992; Pei 1992; Whittet 1992; Dwek 1998; Meyer, Jura & Cardelli 1998; Pei, Fall & Hauser 1999; Weingartner & Draine 2001; James et al. 2002; Kimura, Mann & Jessberger 2003; Draine et al. 2007; Jenkins 2009; Peeples et al. 2014; McKinnon, Torrey & Vogelsberger 2016; Wiseman et al. 2017; Telford et al. 2019). As such, it seems fair to conclude that ϵd is significantly better constrained than κd – making the former a useful tool for pinning down the value of the latter. And whilst some authors suggest larger values of ϵd (for instance De Cia et al. 2013, who find values in the region of 0.8), we can at least be confident that, by definition, no galaxy has a dust-to-metals ratio greater than 1 – no such helpful constraint exists for κd. Furthermore, thanks to observations of elemental depletions in the neutral ISM, ϵd can be determined far more directly than κd.
The formulation in equation (3) gives a combined κd value, that incorporates the contribution from all dust species present, for each temperature component (for n > 1). The problem becomes unconstrained if each dust component is treated as having a different κd. The potential impact of line-of-sight mixing of dust components at different temperatures is discussed in Section 4.2.
It is important to note that |$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$| measurements trace gas-phase metallicity in the ionized phase (predominantly H ii regions), whereas we are concerned with the metallicity of the ISM at large. This means that we must account for the fraction of interstellar oxygen mass in H ii regions depleted on to dust grains, δO, and hence missed by gas-phase metallicity estimators. We use a value of δO = 1.32 ± 0.09 from Mesa-Delgado et al. (2009), which is in good agreement with numerous other reported values (Peimbert & Peimbert 2010; Kudritzki et al. 2012; Bresolin et al. 2016). Whilst the oxygen depletion factor in the ISM at large is known to vary by at least 0.3 dex (Jenkins 2009), oxygen depletion in H ii regions is found to be remarkably constant, at ∼1.3 (i.e. ∼0.1 dex) across nearby galaxies (evaluated by comparing abundances in H ii regions to abundances in the atmospheres of nearby B stars; Bresolin et al. 2016 and references therein). Additionally, given that the elemental composition of oxygen-rich dust is found to exhibit minimal variation at intermediate-to-high metallicities (Mattsson et al. 2019), the assumption of a constant δO is valid modulo a constant ϵd – which is the central premise of our method.
The Amorín et al. (2016) rule is calibrated on a sample of galaxies spanning over an order of magnitude in metallicity (|$7.69 \lt 12 + {\rm log}_{10} [ \frac{\rm O}{\rm H} ] \lt 8.74$|), by using the star formation efficiency (SFE) and star formation rate (SFR) to infer the molecular gas supply present. They do this by employing the relation |$\frac{\alpha _{\rm CO}}{\alpha _{\rm CO_{\rm MW}}} = \tau _{\rm H_{2}} \frac{\rm SFR}{M_{\rm H_{2}}}$|; effectively inverting the Kennicutt–Schmidt law (Kennicutt 1998) to infer the molecular gas mass present, anchored by the known SFE of the Milky Way. Resolved studies such as Bigiel et al. (2011) and Utomo et al. (2019) find remarkably little variation in SFE within face-on local normal spirals like those studied in this work; this supports the reliability of using an SFE-calibrated method for estimating αCO in a resolved study such as ours. Additionally, the Amorín et al. (2016) prescription effectively traces the median of the commonly cited metallicity-dependent literature prescriptions (see fig. 11 of Amorín et al. 2016 and fig. 6 of Accurso et al. 2017 for comparisons of prescriptions), making it the choice most likely to not conflict with other works.
Regarding the Solar metallicity, we use the canonical value for the Solar oxygen abundance of |${}[12 + {\rm log}_{10} \frac{\rm O}{\rm H} ] _{\odot }= 8.69 \pm 0.05$| (Asplund et al. 2009), corresponding to a Solar metal mass fraction of |$f_{Z_{\odot }} = 0.0134$| (Asplund et al. 2009, uncertainty deemed to be negligible). In common with the literature at large, we assume that oxygen abundance traces total metallicity. Whilst this assumption has its limits, oxygen is the most abundant metal in the Universe, and a dominant constituent of dust (Savage & Sembach 1996; Jenkins 2009), making it a useful metallicity tracer for our purposes. Although the ratio of oxygen to carbon (the other main constituent of dust by mass) is known to vary with metallicity (Garnett et al. 1995), this systematic trend is no more prominent than the intrinsic scatter over the 0.4–1.0 Z⊙ metallicity range relevant to this work (Pettini et al. 2008; Berg et al. 2016).
Although a D2 term appears in equation (3), the MH i and |$M_{\rm H_{2}}$| terms are also both proportional to D2, which therefore ultimately cancels out. This renders the resulting values of κλ independent of distance, removing a potentially large source of uncertainty.
Throughout this work, when employing values from the literature, we take care to only use values that do not themselves rely upon any assumed value of κd.
For the value of the dust-to-metals ratio, ϵd, in equation (3), we take two approaches. For our fiducial analysis, presented in Section 5, we assume a constant value of ϵd = 0.4 ± 0.2. This is smaller than the value of 0.5 assumed in Clark et al. (2016), as more recent works (De Cia et al. 2016; McKinnon et al. 2016; Wiseman et al. 2017) suggest that for most galaxies with metallicities > 0.1 Z⊙, the dust-to-metals ratio is slightly below the Milky Way’s average value of 0.5 (James et al. 2002; Jenkins 2009).
The assumption of a constant dust-to-metals ratio is an approximation that will break down at some point. Therefore, in Section 6.2.1, we construct an alternate analysis where ϵd increases as a function of ISM surface density. This is a more physical treatment, as depletion of ISM metals on to dust grains is found to increase in regions of greater ISM column density (Jenkins 2009; Roman-Duval et al. 2019). This is in agreement with the fact that grain growth in the ISM is required to explain the dust budgets in many galaxies (Galliano, Dwek & Chanial 2008; Rowlands et al. 2014; Zhukovska 2014). As a result, dust grain growth in denser ISM (with the corresponding increase in ϵd) is a feature of dust evolution models such as The Heterogeneous dust Evolution Model for Interstellar Solids (THEMIS; Jones et al. 2013, 2017; Jones 2018). Unfortunately, the exact form of the relationship between ϵd and ISM (surface) density is very poorly constrained (the relationship we assume for our analysis is described in detail in Section 6). As such, the variable-ϵd model represents a more physical, but worse-constrained approach; whilst the fixed-ϵd model represents a less physical, but better constrained approach. For this reason, whilst the fixed-ϵd approach is our fiducial model, we none the less consider both scenarios.
3 DATA
An initial attempt by Clark et al. (2016) to detect variation in κd using the dust-to-metals method was unsuccessful; however, that study only considered the global dust properties of galaxies, and considered a sample of 22 objects, all of which were of similar masses, metallicities, and environments. A promising avenue for finding variation in κd is to look within well-resolved nearby galaxies. Many studies have found that dust properties can vary significantly – and sometimes dramatically – within galaxies (Smith et al. 2012; Roman-Duval et al. 2017; Relaño et al. 2018). It would be surprising if this variation did not extend to κd.
Creating a κd map of a galaxy using the dust-to-metals method requires resolved data for its dust emission, atomic gas, molecular gas, and metallicity; with the resolution provided by modern observations, it is possible to make many hundreds, or even thousands, of independent κd determinations within a galaxy. For this proof-of-concept demonstration we map κd within two nearby face-on spiral galaxies – M 74 (NGC 628) and M 83 (NGC 5236). We select these galaxies on account of their particularly extensive metallicity data (see Section 3.3), coupled with their resolution-matched multiphase ISM observations (see Section 3.4).
We obtained the bulk of the necessary data from the DustPedia archive.1 DustPedia (Davies et al. 2017) is a European Union project working towards a comprehensive understanding of dust in the local Universe, capitalizing on the legacy of the Herschel Space Observatory (Pilbratt et al. 2010). A centrepiece of the project is the DustPedia data base, which includes every galaxy observed by Herschel that has recessional velocity within |$3000\, {\rm km\, s^{-1}}$| (∼40 Mpc), has optical angular size in the range 1 arcmin < D25 < 1°, and has a detected stellar component.2
The continuum data we employ are described in Section 3.2, the metallicity data (and the process by which we use it to create metallicity maps) are described in Section 3.3, and the atomic and molecular gas data in Section 3.4.
3.1 Target galaxies
We selected M 74 and M 83 as the subject galaxies for this work; a summary of their basic characteristics is provided in Table 1. Both are very nearby, highly extended, and almost perfectly face-on, making them two of the most heavily studied galaxies in the sky, and ideally suited to serving as our proof-of-concept targets for mapping κd.
. | M 74 . | M 83 . |
---|---|---|
NGC No | NGC 628 | NGC 5236 |
RA (J2000) | 24.174° | 204.254° |
(01h 36m 41|${^{\rm s}_{.}}$| 8) | (13h 37m 01|${^{\rm s}_{.}}$| 0) | |
Dec. (J2000) | + 15.783° | −29.866° |
(+ 15°46′ 58|${^{\prime\prime}_{.}}$| 8) | (−29° 51′ 57|${^{\prime\prime}_{.}}$| 6) | |
Distance (Mpc) a | 10.1 | 4.9 |
Hubble Type | SAc | SBc |
(5.2) | (5.0) | |
D25 (arcmin) | 10.0 | 13.5 |
D25 (kpc) | 29.4 | 19.2 |
A25 (kpc2) | 683 | 290 |
M* (|${\rm log_{10}\, M_{\odot }}$|)b | 10.1 | 10.5 |
MH i (|${\rm log_{10}\, M_{\odot }}$|)c | 9.9 | 10.0 |
|${M_{\rm H_{2}}}$| (|${\rm log_{10}\, M_{\odot }}$|)d | 9.4 | 9.5 |
Md (|${\rm log_{10}\, M_{\odot }}$|)e | 7.5 | 7.4 |
SFR (|${\rm M_{\odot }\, yr^{-1}}$|)b | 2.4 | 6.7 |
FUV-KS (mag) | 2.9 | 3.4 |
NUV-r (mag) | 2.5 | 2.8 |
. | M 74 . | M 83 . |
---|---|---|
NGC No | NGC 628 | NGC 5236 |
RA (J2000) | 24.174° | 204.254° |
(01h 36m 41|${^{\rm s}_{.}}$| 8) | (13h 37m 01|${^{\rm s}_{.}}$| 0) | |
Dec. (J2000) | + 15.783° | −29.866° |
(+ 15°46′ 58|${^{\prime\prime}_{.}}$| 8) | (−29° 51′ 57|${^{\prime\prime}_{.}}$| 6) | |
Distance (Mpc) a | 10.1 | 4.9 |
Hubble Type | SAc | SBc |
(5.2) | (5.0) | |
D25 (arcmin) | 10.0 | 13.5 |
D25 (kpc) | 29.4 | 19.2 |
A25 (kpc2) | 683 | 290 |
M* (|${\rm log_{10}\, M_{\odot }}$|)b | 10.1 | 10.5 |
MH i (|${\rm log_{10}\, M_{\odot }}$|)c | 9.9 | 10.0 |
|${M_{\rm H_{2}}}$| (|${\rm log_{10}\, M_{\odot }}$|)d | 9.4 | 9.5 |
Md (|${\rm log_{10}\, M_{\odot }}$|)e | 7.5 | 7.4 |
SFR (|${\rm M_{\odot }\, yr^{-1}}$|)b | 2.4 | 6.7 |
FUV-KS (mag) | 2.9 | 3.4 |
NUV-r (mag) | 2.5 | 2.8 |
As a first-order estimate of the uncertainty on the distance, we use the standard deviation of the redshift-independent distances listed in the Nasa/ipac Extragalactic Data base (NED; https://ned.ipac.caltech.edu/ui/) for each galaxy. This gives uncertainties of 3.2 and 3.4 Mpc for M 74 and M 83, respectively.
Nersesian et al. (2019).
H i mass from total single-dish flux in the HI Parkes All Sky Survey (HIPASS; Meyer et al. 2004; Wong et al. 2006).
This work (see Section 3.4).
This work (using the pixel-by-pixel κd values calculated in produced 5).
. | M 74 . | M 83 . |
---|---|---|
NGC No | NGC 628 | NGC 5236 |
RA (J2000) | 24.174° | 204.254° |
(01h 36m 41|${^{\rm s}_{.}}$| 8) | (13h 37m 01|${^{\rm s}_{.}}$| 0) | |
Dec. (J2000) | + 15.783° | −29.866° |
(+ 15°46′ 58|${^{\prime\prime}_{.}}$| 8) | (−29° 51′ 57|${^{\prime\prime}_{.}}$| 6) | |
Distance (Mpc) a | 10.1 | 4.9 |
Hubble Type | SAc | SBc |
(5.2) | (5.0) | |
D25 (arcmin) | 10.0 | 13.5 |
D25 (kpc) | 29.4 | 19.2 |
A25 (kpc2) | 683 | 290 |
M* (|${\rm log_{10}\, M_{\odot }}$|)b | 10.1 | 10.5 |
MH i (|${\rm log_{10}\, M_{\odot }}$|)c | 9.9 | 10.0 |
|${M_{\rm H_{2}}}$| (|${\rm log_{10}\, M_{\odot }}$|)d | 9.4 | 9.5 |
Md (|${\rm log_{10}\, M_{\odot }}$|)e | 7.5 | 7.4 |
SFR (|${\rm M_{\odot }\, yr^{-1}}$|)b | 2.4 | 6.7 |
FUV-KS (mag) | 2.9 | 3.4 |
NUV-r (mag) | 2.5 | 2.8 |
. | M 74 . | M 83 . |
---|---|---|
NGC No | NGC 628 | NGC 5236 |
RA (J2000) | 24.174° | 204.254° |
(01h 36m 41|${^{\rm s}_{.}}$| 8) | (13h 37m 01|${^{\rm s}_{.}}$| 0) | |
Dec. (J2000) | + 15.783° | −29.866° |
(+ 15°46′ 58|${^{\prime\prime}_{.}}$| 8) | (−29° 51′ 57|${^{\prime\prime}_{.}}$| 6) | |
Distance (Mpc) a | 10.1 | 4.9 |
Hubble Type | SAc | SBc |
(5.2) | (5.0) | |
D25 (arcmin) | 10.0 | 13.5 |
D25 (kpc) | 29.4 | 19.2 |
A25 (kpc2) | 683 | 290 |
M* (|${\rm log_{10}\, M_{\odot }}$|)b | 10.1 | 10.5 |
MH i (|${\rm log_{10}\, M_{\odot }}$|)c | 9.9 | 10.0 |
|${M_{\rm H_{2}}}$| (|${\rm log_{10}\, M_{\odot }}$|)d | 9.4 | 9.5 |
Md (|${\rm log_{10}\, M_{\odot }}$|)e | 7.5 | 7.4 |
SFR (|${\rm M_{\odot }\, yr^{-1}}$|)b | 2.4 | 6.7 |
FUV-KS (mag) | 2.9 | 3.4 |
NUV-r (mag) | 2.5 | 2.8 |
As a first-order estimate of the uncertainty on the distance, we use the standard deviation of the redshift-independent distances listed in the Nasa/ipac Extragalactic Data base (NED; https://ned.ipac.caltech.edu/ui/) for each galaxy. This gives uncertainties of 3.2 and 3.4 Mpc for M 74 and M 83, respectively.
Nersesian et al. (2019).
H i mass from total single-dish flux in the HI Parkes All Sky Survey (HIPASS; Meyer et al. 2004; Wong et al. 2006).
This work (see Section 3.4).
This work (using the pixel-by-pixel κd values calculated in produced 5).
Both galaxies are classified as ‘grand design’ (Elmegreen & Elmegreen 1987) type Sc spirals, with M 83 also displaying a prominent bar (de Vaucouleurs et al. 1991). M 74 has a physical diameter of 29 kpc – similar to that of the Milky Way (Goodwin, Gribbin & Hendry 1998; Rix & Bovy 2013) – and about 50 per cent greater than that of M 83 (diameter defined according to the optical D25, being the isophotal major axis at which the optical surface brightness falls beneath 25 |${\rm mag\, arcsec^{2}}$|).
Despite being the physically smaller of the two, M 83 has a stellar mass 2.2 times greater, and an SFR 2.7 times greater (Nersesian et al. 2019). M 83 has a correspondingly higher surface brightness in dust emission, averaging 4.2 |${\rm MJy\, sr^{-1}}$| at 500 |$\mu$|m within its D25, compared to 1.6 |${\rm MJy\, sr^{-1}}$| for M 74. The nuclear region of M 83 is currently undergoing a bar-driven starburst, concentrated in the central 250 pc, accounting for ∼10 per cent of the galaxy’s total ongoing star formation (Sérsic & Pastoriza 1965; Harris et al. 2001; Fathi et al. 2008). The optical disc of M 83 has a minimal systematic metallicity gradient, with oxygen abundances varying by only about 0.1 dex from place to place; in contrast, M 74 has a pronounced metallicity gradient, with oxygen abundances in its centre about 0.3 dex greater than at its R25 (De Vis et al. 2019).
Many of the differences between M 74 and M 83 – such as in their stellar surface densities (and therefore interstellar radiation fields), star formation characteristics, metallicity profiles, ISM distributions, etc., – have the potential to affect dust properties, and thereby provide useful scope for us to contrast how κd can vary due to a range of factors.
The appearances of both galaxies, in various parts of the spectrum, are illustrated in Figs 2 and 3. The stellar masses and SFRs for the DustPedia galaxies, as presented in Nersesian et al. (2019), were estimated using the Code Investigating GALaxy Emission (cigale; Burgarella, Buat & Iglesias-Páramo 2005; Noll et al. 2009) software, incorporating the THEMIS dust model.
3.2 Continuum data
Multiwavelength imagery and photometry for the DustPedia galaxies (spanning 42 ultraviolet–millimetre bands), along with distances, morphologies, etc., are presented in Clark et al. (2018). Our analysis makes use of observations from several of the facilities included in the DustPedia archive.
In the submm, we use observations at 250, 350, and 500 |$\mu$|m from the Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) instrument onboard Herschel. In the FIR, we use observations at 160 and 70 |$\mu$|m from the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) instrument, also onboard Herschel (PACS did not perform 100 |$\mu$|m observations for M 83, so for consistency we make no use of the PACS 100 |$\mu$|m data for M 74). In the MIR, we use observations at 22 |$\mu$|m from the WISE.3 A compilation of the MIR–FIR–submm data for each galaxy is shown in the centre left panels of Figs 2 and 3.
Although not required for the creation of the κd maps, we use various additional data for reference and comparison, also drawn from the DustPedia archive. This includes UV observations from GALaxy Evolution eXplorer (GALEX; Morrissey et al. 2007); UV, optical, and NIR observations from the Sloan Digital Sky Survey (SDSS; York et al. 2000; Eisenstein et al. 2011); optical observations from the Digitized Sky Survey (DSS); plus NIR observations from the InfraRed Array Camera (IRAC; Fazio et al. 2004) and Multiband Imager for Spitzer (MIPS; Rieke et al. 2004) instruments onboard the Spitzer Space Telescope (Werner et al. 2004). A compilation of the UV–optical–NIR data for each galaxy is shown in the far-left panels of Figs 2 and 3.
3.3 Metallicity data
Galaxies sufficiently extended to have well-resolved global FIR–submm observations, atomic gas observations, and molecular gas observations, are generally too extended to have their UV–NIR nebular spectral emission – and hence metallicities – fully mapped by Integral Field Unit (IFU) spectrometry. Whilst some large-area IFU surveys of nearby galaxies have now been undertaken, these are still very much the exception rather than the rule, and even the very largest can currently only cover ∼50 per cent of the area of galaxies as extended as M 74 and M 83. (Rosales-Ortega et al. 2010; Sánchez et al. 2011; Blanc et al. 2013). As such, the few DustPedia galaxies with mostly complete IFU coverage do not have the well-resolved gas and dust data needed for this analysis.
However, extended nearby galaxies are popular targets for spectroscopic observation; most have had large numbers of individual slit and fibre spectra taken, supplementing partial IFU coverage like that described above. For DustPedia, De Vis et al. (2019) have compiled a sizeable data base of emission-line fluxes, collated from 42 literature studies plus all available archival Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) data that covers the DustPedia galaxies. The De Vis et al. (2019) spectroscopic data base contains emission-line fluxes from 10 000 spectra, with data for 492 (56 per cent) of the DustPedia galaxies. De Vis et al. (2019) also present consistent gas-phase metallicity measurements for all of these spectra, for five different strong-line relation prescriptions (all of which yield standard |$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$| metallicities). Following their tests of the internal consistency of the prescriptions considered, De Vis et al. (2019) find the Pilyugin & Grebel (2016) ‘S’ prescription most reliable; we therefore use these metallicities throughout the rest of this work. A recent study by Ho (2019) also supports the validity of the Pilyugin & Grebel (2016) prescriptions at the metallicities of our target galaxies. As an additional test, we also repeat the entire κd-mapping process using metallicity data produced using four other strong-line relations; this is presented in Appendix F.
M 74 and M 83 both have large numbers of metallicities in the De Vis et al. (2019) data base – 510 and 793 measurements, respectively, more than any other DustPedia galaxy (except UGC 09299, which lacks the resolved gas data we require). These metallicity points sample the entirety of both galaxies’ optical discs. The positions of these spectra, and the metallicities derived from them, are plotted in the upper left panels of Figs 5 and 6. Our region of interest for each galaxy4 extends approximately out to 0.55 R25 for M 74, and to 0.7 R25 for M 83. So whilst the bulk of the metallicity points lie within the region of interest of each galaxy, providing dense sampling, there are also sufficient points outside it to constrain the metallicity variations over larger scales.
In order to produce maps of κd, it was necessary to first have maps of the metallicity distributions of our target galaxies. The first step towards achieving this was modelling their radial metallicity profiles. The spectra metallicity points for M 74 and M 83, plotted as a function of their deprojected galactocentric radius, r, are shown in Fig. 4. As can be seen, there is significant scatter around the radial trends of both galaxies, far in excess of what would be expected if it were driven solely by the uncertainties on the individual metallicity points. Indeed, if one fits a naïve metallicity profile where the only variables are the gradient and the central metallicity, then the majority of data points would count as > |$5\, \sigma$| ‘outliers’ in M 83 (and most would count as > |$2\, \sigma$| outliers for M 74). This scatter represents localized variations in metallicity, which are not azimuthally symmetric – and which therefore cannot be captured by a one-dimensional model. Such variation becomes apparent when sampling the metallicity within galaxies at such high spatial resolution (Moustakas et al. 2010; Rosales-Ortega et al. 2010). For example, note the localized region of significantly depressed metallicity in the western part5 of the disc of M 83, visible in the upper left panel of Fig. 6.
We had to take this intrinsic scatter into account when modelling the radial metallicity profiles of our target galaxies; we therefore used a model with 3 parameters: the metallicity gradient mZ (in |${\rm dex\, r^{-1}_{25}}$|), the central metallicity cZ (in |$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|), and the intrinsic scatter ψ (in dex). We employed a Bayesian Monte Carlo Markov Chain (MCMC) approach to fit this model, the full details of which are given in Appendix A; the resulting parameter estimates, with uncertainties, are listed in Table 2.
. | M 74 . | M 83 . |
---|---|---|
mZ (|${\rm dex\, r^{-1}_{25}}$|) | −0.27 ± 0.04 | −0.14 ± 0.02 |
cZ (|$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|) | 8.59 ± 0.02 | 8.62 ± 0.01 |
ψ (dex) | 0.044 ± 0.01 | 0.048 ± 0.01 |
. | M 74 . | M 83 . |
---|---|---|
mZ (|${\rm dex\, r^{-1}_{25}}$|) | −0.27 ± 0.04 | −0.14 ± 0.02 |
cZ (|$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|) | 8.59 ± 0.02 | 8.62 ± 0.01 |
ψ (dex) | 0.044 ± 0.01 | 0.048 ± 0.01 |
. | M 74 . | M 83 . |
---|---|---|
mZ (|${\rm dex\, r^{-1}_{25}}$|) | −0.27 ± 0.04 | −0.14 ± 0.02 |
cZ (|$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|) | 8.59 ± 0.02 | 8.62 ± 0.01 |
ψ (dex) | 0.044 ± 0.01 | 0.048 ± 0.01 |
. | M 74 . | M 83 . |
---|---|---|
mZ (|${\rm dex\, r^{-1}_{25}}$|) | −0.27 ± 0.04 | −0.14 ± 0.02 |
cZ (|$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|) | 8.59 ± 0.02 | 8.62 ± 0.01 |
ψ (dex) | 0.044 ± 0.01 | 0.048 ± 0.01 |
It would technically be possible to create metallicity maps of our target galaxies using only these fitted radial metallicity profiles. However, using this simple one-dimensional approach (i.e. where metallicity varies only as a function of r) leads to very large uncertainties on the metallicity value of each pixel in the resulting maps, thanks to the considerable intrinsic scatter values (|$\psi = 0.044\, {\rm dex}$| for M 74, and |$\psi = 0.049\, {\rm dex}$| for M 83). In contrast, most of the individual spectra metallicity data points have uncertainties much smaller than this, with median uncertainties of 0.010 and 0.025 dex for M 74 and M 83, respectively (NB, spectra located in close proximity tend to have metallicities that are in good agreement – see the densely sampled area in Figs 5 and 6). In other words, there are many areas of these galaxies where the metallicity is known to much greater confidence than is reflected by the global radial metallicity gradient – therefore, relying upon the global one-dimensional model alone would mean ‘throwing away’ that information. As such, we opted to model the metallicity distributions of our target galaxies in two dimensions. To achieve this, we employed Gaussian process regression.
3.3.1 Gaussian process regression
Gaussian process regression (GPR) is a form of probabilistic interpolation that makes it possible to model a data set without having to assume any sort of underlying functional form for the model. GPR (and Gaussian process methodology in general) is a commonly applied tool in the field of machine learning – and in recent years GPR has seen increasing use in astronomy, to tackle problems where stochastic (and therefore impractical to model directly) processes give rise to complex features in data (for instance, capturing the effect of varying detector noise levels in time-domain data). For a full introduction to Gaussian process methodology, including GPR, see Rasmussen & Williams (2006); for an extensive list of works where Gaussian processes have been successfully applied to problems in astronomy, see section 1 of Angus et al. (2018).
Instead of trying to model the underlying function that gave rise to the observed data, GPR models the covariance between the data points. The covariance is modelled using a kernel, which describes how the values of data points are correlated with one another, as a function of their separation in the parameter space.
This covariance-modelling approach is well suited to the problem we face with mapping metallicity within our target galaxies. Spectra located very close together (e.g. within a few arcseconds) will tend to have very similar metallicities, whilst spectra with greater separations (e.g. arcminutes apart) will only be weakly correlated with one another (this is readily apparent from visual inspection of Figs 5 and 6).
For the covariance function, we used a Mátern kernel (Stein 1999). The Mátern function is a standard choice for modelling the spatial correlation of two-dimensional data (Minasny & McBratney 2005; Rasmussen & Williams 2006; Cressie & Wikle 2011) – especially physical data (Schön et al. 2018). In practice, a Mátern kernel is similar to a Gaussian kernel, but has a narrower peak (allowing it to be sensitive to variations over short distances) whilst also having thicker tails (letting it maintain sensitivity to the covariance over large distances). Like a Gaussian, the tails extend to infinity. The Mátern kernel has two hyperparameters: kernel scale and kernel smoothness (essentially how ‘sharp’ the peak of the kernel is).
Once the covariance has been modelled, it is used in combination with the observed data to trace the underlying distribution. The result is a full posterior probability distribution function (PDF) for the likely value of the underlying function at that location. The uncertainties in each input data point are fully considered by GPR. In regions where the input data points have large uncertainties, or where data points in close proximity disagree with one another, the output PDF will be less well constrained, reflecting the greater uncertainty on the underlying value at that location.
3.3.2 Metallicity maps via Gaussian process regression
We opted to apply the GPR to the residuals between the individual spectra metallicity points and the global radial metallicity profile (i.e. Fig. 4). By fitting to the residuals, the global radial metallicity profile effectively serves as the prior for the regression. The regression then traces the structure of the local deviations from the global radial metallicity profile. In regions where there are no data points, the GPR therefore tends to revert to the metallicity implied by the global radial profile.
This process is illustrated in the upper right panels of Figs 5 and 6 for M 74 and M 83, respectively. The circular points mark the positions of the individual spectra metallicities, colour coded to show the residual of each (the median absolute residual is 0.026 dex for both galaxies). The coloured background shows the Gaussian process regression to these residuals, similarly colour coded. We used GaussianProcessRegressor, the GPR implementation of the Scikit-Learn machine learning package for python (Pedregosa et al. 2011). The hyperprior for the kernel scale was flat, but limited to a range of 0.05–0.5 D25, to prevent the modelled regression being either featurelessly smooth, or unrealistically granular. The kernel smoothness hyperprior was set to 1.5, which is a standard choice due to being computationally efficient, differentiable, and often found to be effective in practice (Rasmussen & Williams 2006; Gatti 2015).
The final metallicity map for each galaxy was produced by adding the residual distribution traced by the GPR to the global radial metallicity profile, for each pixel. The resulting metallicity maps are plotted as contours in the lower right panels of Figs 5 and 6, for M 74 and M 83, respectively. Visual inspection indicates that the GPR does a good job of tracing the metallicity distribution as sampled by the spectra metallicity points (i.e. the contours consistently have the same levels as the points they pass through).
Our full procedure for calculating the uncertainty on the GPR metallicity in each pixel is presented in Appendix C. The resulting metallicity uncertainty maps are shown in the lower left panels of Figs 5 and 6.
We validated the reliability of the metallicities predicted by GPR by performing a jackknife cross-validation analysis, which is described in detail in Appendix B. This analysis found that the predicted values exhibit no significant bias, and the associated uncertainties are reliable.
There are areas in both galaxies where the data points suggest a steadily increasing residual in a certain direction; the GPR then extrapolates that this increase continues for some distance (defined by the modelled kernel scale) into regions where there are no data points. For instance, in the south-western part of M 74, the data points suggest that the metallicity gradient is steeper than for the rest of the galaxy (i.e. a trend of increasingly negative residuals) – the GPR extrapolates that this increased steepness will continue for a certain distance into an area where there are no metallicity points. A similar situation occurs in the north-west portion of M 83 (but instead with a positive residual). Naturally, extrapolations such as these are highly uncertain; but this is quantified by the uncertainty on the regression at these locations. This is illustrated in the lower left panels of Figs 5 and 6, which show the uncertainty for each pixel’s predicted metallicity.
Utilizing GPR provides a marked reduction in the uncertainty of our metallicity maps, relative to using the global radial metallicity profiles alone. If we were to use that simple global approach, every pixel in our metallicity map for M 74 would have an uncertainty at least as large as the intrinsic scatter of 0.044 dex (Table 2). In contrast, with our GPR metallicity map of M 74, 91 per cent of the pixels within the region of interest4 have uncertainties <0.044 dex; the median GPR uncertainty within this region is only 0.016 dex. Similarly, whereas the intrinsic scatter on the global radial profile of M 83 is 0.048 dex, the median error on the GPR metallicity map is only 0.037 dex within the region of interest; the GPR uncertainty is less than the global intrinsic scatter for 66 per cent of the pixels within this region.
There exist ‘direct’ electron temperature metallicity measurements for M 74, produced by the CHemical Abundances Of Spirals (CHAOS; Berg et al. 2015). Electron temperature metallicities are at reduced risk of systematic errors, compared to strong-line values like those provided by De Vis et al. (2019). However, the CHAOS data for M 74 only consists of 45 measurements. Whilst we trialled producing metallicity maps with these data, the sparse sampling meant that the uncertainty on the metallicity at any given point was extremely large. Maps of κd produced with these metallicity maps (as per the procedure described in Section 4) were so dominated by the resulting noise that they were not informative.
3.4 Atomic and molecular gas data
Atomic and molecular gas data for a sample of extended, face-on spiral galaxies in DustPedia – including those studied in this work – is presented in Casasola et al. (2017). For both of our target galaxies, we followed Casasola et al. (2017) and use H i data from The HI Nearby Galaxy Survey (THINGS; Walter et al. 2008), which conducted 21 cm observations of 34 nearby galaxies with the Very Large Array, at 6–16 arcsec resolution. We retrieved the naturally weighted moment 0 maps for M 74 and M 83 from the THINGS website.6 The H i maps for both galaxies are shown in the third panels of Figs 2 and 3.
To obtain CO observations for M 74 we again followed Casasola et al. (2017), and used data from the HERA Co Line Extragalactic Survey (HERACLES; Leroy et al. 2009), which performed CO(2–1) observations of 18 nearby galaxies using the IRAM 30 m telescope, at 13 arcsec resolution. We retrieved the moment 0 maps, as associated uncertainty maps, from IRAM’s official HERACLES data repository.7 The CO(2–1) map for M 74 is shown in the fourth panel of Fig. 2.
Although M 74 has been observed in CO(1–0) by various authors (Young et al. 1995; Regan et al. 2001), these observations are all lacking in either resolution, sensitivity, and/or coverage, in comparison to the HERACLES data. We therefore found it preferable to use the CO(2–1) data of HERACLES, despite the fact this requires applying a line ratio, r2:1 = ICO(2–1)/ICO(1–0), in order to find ICO(1–0), and hence calculate H2 mass as per equation (6).
In nearby late-type galaxies, r2:1 has an average value of ∼0.7 (Leroy et al. 2013; Casasola et al. 2015; Saintonge et al. 2017). However, it is also known that r2:1 varies significantly with galactocentric radius (Casoli et al. 1991; Sawada et al. 2001; Leroy et al. 2009). As such, accurately inferring the CO(1–0) distribution in M 74 using the HERACLES CO(2–1) map required a radially dependent r2:1. To produce this, we used the data presented in fig. 34 (lower right panel) of Leroy et al. (2009), where they compare the HERACLES ICO(2–1) maps to literature ICO(1–0) maps of the same galaxies produced by several other telescopes (with appropriate corrections applied to account for differences in spatial and velocity resolution). This yielded ≈450 directly measured r2:1 values, spanning radii from 0–0.55 R25, for nine of the HERACLES galaxies. Leroy et al. (2009) simply binned these points to trace the radial variation in r2:1; however, we chose to take a fully probabilistic approach, and use GPR to infer the underlying radial trend in r2:1. In Fig. 7, we plot all of the r2:1 points from fig. 34 (lower right panel) of Leroy et al. (2009). We applied a GPR to these data, using a Mátern covariance kernel. Because r2:1 is a ratio, we constructed the regression so that the output uncertainties are symmetric in logarithmic space; otherwise, output uncertainties symmetric in linear space would extend to unphysical values of r2:1 < 0 at larger radii. The resulting regression is shown in black in Fig. 7. It is in excellent agreement with the radial trend that Leroy et al. (2009) traced by binning the data, with r2:1 elevated to ∼1 in the galaxies’ centres, falling to 0.7–0.8 over the rest of the sampled region – but our approach has the added benefit over binning of providing well-constrained uncertainties on r2:1 values produced using the regression. The uncertainty associated with the regression is a factor of ≈1.3 over the 0 < R/R25 < 0.55 range in radius sampled by the HERACLES measurements, reflecting the intrinsic scatter present in the data points; beyond this, the uncertainty steadily increases, reaching a factor of ≈2 at R = R25. Given the uncertainty on αCO, this does not represent a large addition to the total uncertainty on the molecular gas masses we calculated.
M 83 was not observed by HERACLES. So we instead used the CO(1–0) observations presented in Lundgren et al. (2004), which were made using the Swedish–Eso Submillimetre Telescope (SEST) at a resolution of 42 arcsec, to a uniform depth of 74 mK (Tmb). The CO(1–0) map for M 83 is shown in the far right panel of Fig. 3.
We determined αCO pixel-by-pixel using our metallicity maps according to equation (7), and thereby produced H2 maps of our target galaxies. The total H2 masses contained in these maps are the H2 masses listed in Table 1.
4 APPLICATION
4.1 Data preparation
We background subtracted all continuum maps following the procedure described in Clark et al. (2018), using the background annuli they specify for our target galaxies.
All data (continuum observations, gas observations, and metallicity maps) were smoothed to the resolution of the most poorly resolved observations for each galaxy. This was done by convolving each image with an Airy disc kernel of full width at half-maximum (FWHM) given by |$\smash{\theta _{\rm kernel} = (\theta ^{2}_{\rm worst} - \theta ^{2}_{\rm data}})^{\frac{1}{2}}$|. We therefore convolve all of our M 74 data to the 36 arcsec resolution of the Herschel-SPIRE 500 |$\mu$|m observations. Likewise, we convolved all of our M 83 data to the 42 arcsec resolution of the SEST H i observations.
We reprojected all of our data to a common pixel grid for each galaxy, on an east–north gnomic tan projection. We wished to preserve angular resolution, ensuring that our data remain Nyquist sampled, to maximize our ability to identify any spatial features or trends in our final κd maps. We therefore used projections with 3 pixels per convolved FWHM. This corresponds to 12 arcsec pixels for M 74, and 14 arcsec pixels for M 83.
For each galaxy, we defined a region of interest, within which all required data are of sufficient quality to effectively map κd. We defined this as being the region within which all pixels in the smoothed and reprojected versions of the H i map, CO map, and 22–500 |$\mu$|m continuum maps, have SNR > 2 (as defined by comparison to their respective uncertainty maps). For both M 74 and M 83, the data with the limiting sensitivity are the CO observations. The borders of our regions of interest for both galaxies are shown in the far right panels of Figs 2 and 3.
4.2 SED fitting
As described in Section 2, the dust-to-metals method lets us establish dust masses a priori; then, by comparing this a priori dust mass to observed FIR–submm dust emission, we can calibrate the value of κd. This necessitates having a model that describes that FIR–submm dust emission. We wished to minimize the scope for potentially incorrect model assumptions to corrupt our resulting κd values. We therefore modelled the dust emission with the simplest model that is able to fit FIR–submm fluxes – a one-component MBB (i.e. equation 1, with n = 1). A one-component MBB model has been shown by many authors to break down in various circumstances (e.g. Jones 2013; Clark et al. 2015; Chastenet et al. 2017; Lamperti et al. accepted). However, these primarily concern either submillimetre excess in low-metallicity and/or low-density environments (which are not present in the regions of interest within our target galaxies), the emission from hotter dust components at short wavelengths (which we do not attempt to model; see below), or features only discernible in spectroscopy (which we are not employing). In ‘normal’ galaxies, a one-component MBB can be expected to fit FIR–submm fluxes successfully (Nersesian et al. 2019).
Note that, as a test, we also repeated the entire SED-fitting process described in this section with a two-component MBB model (i.e. equation 1, with n = 2, giving dust components at two temperatures). However, when comparing the χ2 values of both sets of fits, we found that adopting the two-component MBB approach adds little benefit to the quality of the fits. The median reduced χ2 values (of all posterior samples, from all pixels) for the one-component MBB fits were 0.61 for M 74 and 0.94 for M 83 – compared to 0.59 and 0.65, respectively, for the two-component fits. This indicates that the two-component MBB fits offer minimal improvement over the one-component fits (and, indeed, may be straying into the realm of overfitting). Given our desire to employ the simplest applicable model, we therefore opt to proceed with the one-component MBB approach for this work. None the less, in Appendix G, we verify that the choice of one- or two-component SED fitting does not result in considerable changes to our overall results.
By performing our SED fitting pixel-by-pixel, we are reducing the degree to which there will be contributions from multiple dust components at different temperatures. None the less, there will inevitably be some degree of line-of-sight mixing of dust populations. This risk will be greatest in the densest regions, where fainter emission from colder, but potentially more massive, dust components can be dominated by brighter emission from warmer, but less massive, components heated by star formation (Malinen et al. 2011; Juvela & Ysard 2012). If this does occur, then the resulting κd values will, in effect, factor in the mass of any cold dust component too faint to affect the SED (assuming the a priori dust masses calculated by the dust-to-metals method are accurate). In this scenario, the κd values we calculate may not be valid if applied to observations with good enough spatial resolution that line-of-sight mixing becomes negligible.
Although we use equation (1) to model SEDs, we assign an arbitrary value of κλ during the fitting process (as, of course, the SED fitting is being performed in order to allow us to find a value of κλ using the results). This means that the ‘mass’ parameter yielded by our SED fitting merely serves as a normalization term for the SED amplitude. This is not a problem, as the only output values actually required is the temperature of the dust, and its flux at the reference wavelength; these are needed in equation (3) to calculate values of κd.
We also incorporate a correlated photometric error parameter, υSPIRE, into our SED fitting. The photometric calibration uncertainty of the Herschel-SPIRE instrument contains a systematic error component that is correlated between bands (Griffin et al. 2010; Bendo et al. 2013; Griffin et al. 2013). This arises from the fact that Herschel-SPIRE was calibrated using observations of Neptune; however, the reference model of Neptune’s emission has a |$\pm 4 {{\ \rm per\ cent}}$| uncertainty. We account for this by parametrizing the correlated Herschel-SPIRE error as υSPIRE. The |$\pm 4 {{\ \rm per\ cent}}$| scale of υSPIRE accounts for the majority of the combined 5.5 per cent calibration uncertainty of Herschel-SPIRE.8 As such, for high-SNR sources (such as bright pixels within our target galaxies), where the photometric noise is minimal, the correlated calibration error can actually dominate the entire uncertainty budget. Moreover, the |$\pm 4 {{\ \rm per\ cent}}$| error on υSPIRE does not follow the Gaussian or Student’s t distribution typically assumed for photometric uncertainties – rather, it is essentially flat, with the true value of the correlated systematic error almost certainly lying somewhere within the |$\pm 4 {{\ \rm per\ cent}}$| range (Bendo et al. 2013; A. Papageorgiou, private communication; C. North, private communication). Explicitly handling υSPIRE as a nuisance parameter allows us to properly account for this with a matching prior. Gordon et al. (2014) highlight the significant differences that can be found in dust SED fitting when the correlated photometric uncertainties are considered, compared to when they are not.
The Herschel-PACS instrument also has a systematic calibration error, of |$\pm 5 {{\ \rm per\ cent}}$|, arising from uncertainty on the emission models of its calibrator sources, a set of five late-type giant stars (Balog et al. 2014). However, the error budget on the emission models is dominated by the |$\pm 3 {{\ \rm per\ cent}}$| uncertainty on the line features in the atmospheres of the calibrator stars (see table 2 of Decin & Eriksson 2007), which are different in each band, and hence not correlated. Only the uncertainty on the continuum component of the emission model, of ±1–|$2 {{\ \rm per\ cent}}$|, will be correlated between bands. Given the small scale of this correlated error component, and given that systematic error makes up a smaller fraction of the total Herschel-PACS calibration uncertainty than it does for Herschel-SPIRE, and given that the greater instrumental noise for Herschel-PACS means that calibration uncertainty makes up a small fraction of the total photometric uncertainty budget than it does for Herschel-SPIRE, we opt to not model the correlated uncertainty for Herschel-PACS as we do with υSPIRE.
Our one-component MBB SED model therefore has four variables: the dust temperature, Td; the dust ‘mass’ normalization, |$M_{\mathrm{ d}}^{\rm (norm)}$|; the emissivity slope, β; and the correlated photometric error in the Herschel-SPIRE bands, υSPIRE.
We treat photometric uncertainties as being described by a first-order (i.e. one degree of freedom) Student t distribution. The Student t distribution has more weight in the tails than a Gaussian distribution, allowing it to better account for outliers. This makes the Student t distribution a standard choice for Bayesian SED fitting (da Cunha, Charlot & Elbaz 2008; Kelly et al. 2012; Galliano 2018).
For the photometric uncertainty in each pixel, we used the values provided by the uncertainty maps, added in quadrature to the calibration uncertainty of each band: 5.6 per cent for WISE 22 |$\mu$|m,10 7 per cent for Herschel-PACS 70–160 |$\mu$|m,11 and 2.3 per cent12 for Herschel-SPIRE 250–500 |$\mu$|m8. Both of our target galaxies lie in regions with negligible contamination from Galactic cirrus. The WISE and Herschel-PACS backgrounds are dominated by instrumental noise, whilst the Herschel-SPIRE background has a significant contribution from the confused extragalactic background. Therefore, for the Herschel-SPIRE data, we also add in quadrature the contribution of confusion noise; for this we use the values given in Smith et al. (2017), of 0.282, 0.211, and 0.105 |${\rm MJy\, sr^{-1}}$| at 250, 350, and 500 |$\mu$|m, respectively, derived from the Herschel-ATLAS fields (although the instrumental noise level still dominates over this in all of our Herschel-SPIRE data).
We treat fluxes at wavelengths <100 |$\mu$|m as upper limits, as emission in this regime will include contributions from hot dust and stochastically heated small grains (Boulanger & Perault 1988; Desert, Boulanger & Puget 1990; Jones et al. 2013) that will not be accounted for by our MBB model. Therefore at these wavelengths, any proposed model flux that falls below the observed flux will be deemed as likely as the observed flux itself (i.e. no proposed model will be penalized for underpredicting the flux in these bands). Only for proposed model fluxes greater than the observed flux will the likelihood decrease according to the Student t distribution, as per usual.
We sample the posterior probability distribution of the SED model parameters in each pixel using the emcee (Foreman-Mackey et al. 2013) MCMC package for python. We perform 750 steps with 500 chains (‘walkers’); the first 500 steps from each chain were discarded as burn-in, and non-convergence was checked for using the Geweke diagnostic13 (Geweke 1992). Our priors are detailed in Appendix D.
Our SED-fitting routine incorporates colour corrections to account for the effects of the instrumental filter response functions and beam areas.14,15,16,17 An example posterior SED, along with the corresponding parameter distributions, are shown in Figs 8 and 9.
Figs 10 and 11 show maps of the median values of dust mass surface density, temperature, and β values for each pixel. We assume that the low temperatures and large β values found in the centre of M 83 are non-physical, and instead are due to non-thermal emission from the nuclear starburst affecting the SED fitting. This is limited to a beam-sized area, consisting of 9 pixels – we therefore exclude these pixels from analysis in later sections, where noted.
Unsurprisingly, the maps of dust mass surface density closely match the morphology of the dust emission (see Figs 2 and 3). The temperature map for M 74 is ‘blotchy’, with warmer dust being located around areas of particularly active star formation (compare to the regions of bright MIR emission in Fig. 2 in the northern and southern parts of the disc). The temperature map for M 83 more visibly traces the overall spiral structure; in particular, elevated temperatures are found on the exterior edges of the spiral arms. The β maps for both galaxies show correlations with the dust mass surface density; in M 74 this manifests as a broad global trend of beta decreasing with radius, whilst in M 83 beta again more obviously traces the spiral structure.
There is a well-known anticorrelation between temperature and β when performing MBB SED fits (Shetty et al. 2009; Kelly et al. 2012; Galliano, Galametz & Jones 2018). This is clearly in evidence in Fig. 9. However, as demonstrated by Smith et al. (2012), this does not introduce systematic errors into the results of such fits. And given this lack of systematic bias, the anticorrelation will not introduce spurious trends into resolved SED fits – because fits separated by more than one beam width will be independent, and will be no more likely to be biased one way than the other. Combined with the fact that we sample the full posterior in our SED fits, and propagate this into the final calculation of our κd maps (see Section 5), we do not believe that the temperature–β anticorrelation will compromise the validity of our final results.
Our SED-fitting code has been made freely available online as a python 3 package.18
5 RESULTS
We now have the atomic gas, molecular gas, metallicity, and dust emission data necessary for every pixel in order to create maps of κd for our target galaxies.
For every pixel within the region of interest for each galaxy, we produced a full posterior probability distribution for κd. We did this by drawing random samples from the posterior distributions provided by our SED and metallicity maps (which are independent of one another), and inputting them into equation (3) (with number of MBB SED components i = 1, as per Section 4.2). For all other input values (SH i, ICO, αCO, |$\alpha _{\rm CO_{\rm MW}}$|, |$y_{\rm \, CO}$|, r2:1, δO, |$f_{Z_{\odot }}$|, |${}[12 + {\rm log}_{10} \frac{\rm O}{\rm H}]_{\odot }$|, |$f_{\rm He_{p}}$|, |${}[\frac{\Delta f_{\rm He}}{\Delta f_{Z}}]$|, and ϵd) we drew random samples from the Gaussian distributions described by their adopted values and associated uncertainties (effectively assuming flat priors, so that these can be treated as posterior probabilities).
We calculated κd for a reference wavelength of 500 |$\mu$|m, as this is the longest wavelength for which we have data, and therefore the wavelength where emission is least sensitive to dust temperature; this minimizes the degree to which uncertainty in temperature is propagated to κd. Our resulting maps of κ500, produced by taking the posterior median in each pixel, are shown in Figs 12 and 13. These maps contain 585 and 1269 pixels for M 74 and M 83, respectively. Throughout the rest of this work, quoted κ500 values are pixel medians. The overall median across M 74 is κ500 = 0.15 |${\rm m^{2}\, kg^{-1}}$|, whilst the overall median across M 83 is κ500 = 0.26 |${\rm m^{2}\, kg^{-1}}$|.
The uncertainties on these κ500 values (defined by the 68.3 per cent quantile in absolute deviation away from the median along the posterior distribution) span the range 0.21–0.28 dex, with a mean uncertainty of 0.25 dex for both galaxies. Note that a large degree of this uncertainty is shared across all pixels, due to the contributions of systematics (such as the uncertainties on ϵd, |$\alpha _{\rm CO_{\rm MW}}$|, etc.), which is why the 0.25 dex average uncertainty is large relative to the scatter in κ500 values. We determined the contribution of the systematic components to the overall uncertainty via a Monte Carlo simulation, in which κ500 values were generated according to equation (3), but where only input parameters with systematic uncertainties were allowed to vary. The scatter on the output dummy values of κ500 was taken to represent the total systematic uncertainty. On average, we found that the systematic components contribute 0.20 dex to the uncertainty. Taking the quadrature difference between this and our average total uncertainty gives an average statistical uncertainty of 0.15 dex in κ500.
The values in our κ500 maps are not fully independent, as they have a pixel width of 3 pixels per FWHM; this will render adjacent pixels correlated. Therefore we also produced a version of the κ500 maps with pixels large enough to be independent (i.e. 1 pixel per FWHM). These maps contained 65 and 141 independent κ500 measurements for M 74 and M 83, respectively. When performing statistical analyses throughout the rest of this work, we used these maps in order to ensure the validity of the results. However, the use of larger pixels for these maps does involve throwing away spatial information. We therefore present the standard, Nyquist-sampled maps in Figs 12, 13, and elsewhere, in order to display all of the spatial information our data are able to resolve. Similarly, individual points plotted in Fig. 14 and elsewhere represent the pixels from the Nyquist-sampled maps, although the trend lines shown on these plots are derived from the independent-pixel data.
In order to calculate a robust estimate of the underlying range of κ500 values, we performed a non-parametric bootstrap resampling of the pixel medians. This non-parametric bootstrap approach will account for the statistical scatter, and not encompass the systematics. This gives a median underlying range for 0.11–0.25 |${\rm m^{2}\, kg^{-1}}$| for M 74 (a factor of 2.3 variation) and 0.15–0.80 |${\rm m^{2}\, kg^{-1}}$| for M 83 (a factor of 5.3 variation).
There is a strong relationship between κ500 and ΣISM (the ISM mass surface density, where |$\Sigma _{\rm ISM} = \Sigma _{\rm H\,{{\small I}}} + \Sigma _{\rm H_{2}} + \Sigma _{\rm d}$|) as shown in Fig. 14. Both galaxies exhibit this relation, but are curiously separated, with the relation for M 74 lying ∼0.3 dex beneath that of M 83. We are able to trace this behaviour over a much larger range of ΣISM for M 83 than for M 74 – the densest regions of M 83 are much denser than those of M 74, whilst the deeper CO data for M 83 allows us to probe to regions of lower density. This neatly accounts for the fact that we find a narrower range of κ500 values for M 74 than M 83 – whilst we probe a 1.7 dex range in density in the latter, we only probe 0.7 dex in the former. We estimated κ500 versus ΣISM power laws for each galaxy by performing a Theil–Sen regression (Theil 1992) to each set of posterior samples in our κ500 and ΣISM maps (specifically, the independent-pixel version of the maps, as discussed above). The resulting power-law slopes for both galaxies are in good agreement, with their indices being |$-0.35^{+0.26}_{-0.21}$| for M 74 and |$-0.36^{+0.04}_{-0.05}$| for M 83. As discussed in Section 6.3, this behaviour is in contradiction to positive correlation between κd and ISM density predicted by standard dust models. The median statistical uncertainty on pixel values of ΣISM is 0.13 dex; given the similarly small 0.15 dex average statistical uncertainty on κ500, we can be confident that the trend in Fig. 14, which spans 1.7 dex for M 83, is not merely a spurious noise induced correlation. The rank correlation coefficient of the relationship is τ = −0.36 for M 74, and τ = −0.57 for M 83 (from a Kendall tau rank correlation test; Kendall & Gibbons 1990).
In Fig.15, we see that it is the overall ISM density that is driving this trend, rather than the density of either the molecular gas, atomic gas, or dust components of the ISM alone, as all three have much weaker relationships with κ500 than is the case for the combined ΣISM. For |$\Sigma _{\rm H_{2}}$|, τM74 = −0.18 and τM83 = −0.55; for Σd, τM74 = −0.28 and τM83 = −0.42; for Σd, τM74 = 0.10 and τM83 = −0.34.
The relationship between κ500 and gas-phase metallicity is plotted in Fig. 16. Once again, whilst M 83 shows no correlation, there does appear to be a trend for M 74, with larger values of κ500 being associated with higher metallicities (|$\mathcal {P}_{\rm null} = 10^{-3.5}$| from a Kendall rank correlation test). On the one hand, metallicity is a parameter in equation (3), so once again there is a definite risk of spurious correlations arising. However, if all other parameters in equation (3) are held fixed, higher metallicity (therefore higher fZ) leads to lower values of κ500, meaning the trend for M 74 in Fig. 16 is being driven by the data in spite of this. Greater ISM metallicity will lead to increased grain growth (Dwek 1998; Zhukovska 2014; Galliano et al. 2018), and larger grains should give rise to larger values of κd (Li 2005; Köhler et al. 2015; Ysard et al. 2018).
We wished to assess whether local star formation has an effect on our calculated values of κ500. There are several mechanisms by which recent star formation can process dust grains in its vicinity (see review in Galliano et al. 2018). For instance, photodestruction by high-energy photons from massive (therefore young) stars can directly break down dust grains (Boulanger et al. 1998; Beirão et al. 2006), whilst the shocks produced by the supernovae of massive stars will sputter dust grains (Bocchio, Jones & Slavin 2014; Slavin, Dwek & Jones 2015). FUV emission should be a good proxy of these two environmental conditions; unobscured FUV emission is indicative of massive stars that are old enough to cleared their birth clouds, and hence represent the regions where supernovae will be occurring. And of course, regions with greater amounts of unobscured FUV emission demonstrably have an interstellar radiation field (ISRF) with greater amounts of high-energy photons. If the environmental effects of recent star formation were impacting κ500, this could manifest as a correlation with the total UV energy density, or with the UV energy density per dust mass (similar to the ‘heating parameter’ of Foyle et al. 2013), as the dust will be better shielded in areas with greater dust density. Therefore, in the two leftmost panels of Fig. 17, we plot κ500 against both the GALEX far-ultraviolet (FUV) luminosity surface density19 (ΣFUV), and against the FUV luminosity per dust mass surface density (ΣFUV/Σd). No trend is apparent in either plot; M 74, with its generally lower values of κ500, has a higher average value of ΣFUV/Σd, but this is to be expected given its bluer colours and lower submm surface brightness (see Table 1).
We also wished to assess whether the ISRF arising from evolved stars could be influencing κ500, given that radiation from evolved stars can be the dominant source of energy received by dust in certain environments (Boquien et al. 2011; Bendo et al. 2012; Nersesian in press). Observations in the NIR provide a good tracer of the evolved stellar population, and the ISRF it produces. Therefore, as with FUV, we plot κ500 against the WISE 3.4 |$\mu$|m luminosity surface density19 (|$\Sigma _{3.4\, \mu \mathrm{ m}}$|), and against the 3.4 |$\mu$|m luminosity per dust mass surface density (|$\Sigma _{3.4\, \mu \mathrm{ m}}/\Sigma _{\rm d}$|), shown in the two rightmost panels of Fig. 17. In M 74, it seems that the pixels with |$\frac{\Sigma _{ 3.4\, \mu \mathrm{ m}}}{\Sigma _{\rm d}} \gt 6\, \times \, 10^{-4}\, {\rm L_{\odot }\, M_{\odot }^{-1}}$| are exclusively associated with higher values of κ500. And most interestingly, there is for both galaxies a positive correlation between κ500 and |$\Sigma _{3.4\, \mu \mathrm{ m}}/\Sigma _{\rm d}$|). Whilst there is appreciable scatter, a Kendall rank correlation test gives |$\mathcal {P}_{\rm null} \gt 0.023$| for both – so it seems that this relationship, whilst broad, has probably not arisen by chance.20 Plus, the WISE 3.4 |$\mu$|m data played no part in our κ500 calculations, making it hard to see how this relation could have arisen spuriously from our methodology.
A downside to using 500 |$\mu$|m as the reference wavelength is that carbonaceous species are expected to have considerably larger κ500 values than silicate species at these longer wavelengths (due to the steeper β for silicates; Ysard et al. 2018). Whereas at shorter wavelengths, the difference in κd between carbonaceous and silicate dust is smaller. Thus the choice of the longer reference wavelength might be limiting our ability to use the κd maps to trace such compositional variation. We therefore also produced versions of our κd maps at a reference wavelength of 160 |$\mu$|m. These κ160 maps are presented in Appendix E; however, they exhibit no difference in structure to the κ500 maps.
6 DISCUSSION
6.1 Robustness of findings
Within M 74 and M 83, we find values of κ500 that vary by factors of 2.3 and 5.3, respectively. This is, to our knowledge, the first observational mapping of variation in κd within other galaxies. However, it is important to critically evaluate how much of this apparent variation could simply be an artefact of our method.
In a companion study to this work, Bianchi et al. (submitted) use the dust-to-metals method to calculate global κd values for 204 DustPedia galaxies. As that study uses integrated gas measurements, they are unable to directly constrain ISM density. However, they do find that galaxies with higher H2/H i ratios (typically associated with denser ISM) tend to have lower values of κd. This is what would be expected if the anticorrelation we find between κd and ΣISM continues on global scales, between galaxies. They also find large (a factor of several) scatter in their κd values between galaxies; in this context, the differences between the values we find for M 73 and M 83 are not conspicuous.
Our key assumption of a fixed dust-to-metals ratio, ϵd, deserves particular scrutiny. As mentioned in Section 2, the vast majority of directly measured21 values of ϵd lie in the range 0.2–0.6. Whilst this factor of 3 variation could notionally, in the worse-case scenario, be sufficient to nullify the factor 2.3 variation in κd we find in M 74, it could not nullify the factor 5.3 variation in M 83. Moreover, as we show in Section 6.2.1, in the physically most likely scenario where ϵd scales with density, the variation in κ500 actually increases. None the less, it is undoubtedly worth considering how, precisely, different kinds of systematic variations in ϵd within our target galaxies could be influencing our results.
There is evidence that ϵd is significantly reduced at low metallicities (Galliano et al. 2005; De Cia et al. 2016; Wiseman et al. 2017). However, there appears to be reduced variation in ϵd at intermediate-to-high metallicity. De Cia et al. (2016) and Wiseman et al. (2017) use depletions in damped Ly α absorbers to find only a factor of ∼2 variation in ϵd at metallicities above 0.1 Z⊙, with at most a weak dependence on metallicity in that regime. Given that our analysis is concerned only with environments at |${\gg} 0.1\, \mathrm{Z}_{\odot }$|, our results should be minimally susceptible to this scale of metallicity effect. Additionally, it should be noted that a number of studies have used visual extinction per column density of metals as a proxy for ϵd, and found it to be constant down to metallicities of 0.01 Z⊙, over a redshift range of 0.1 < z < 6.3 (Watson 2011; Zafar & Watson 2013; Sparre et al. 2014).
A number of simulations have addressed the question of how ϵd varies. McKinnon et al. (2016) trace ϵd in cosmological zoom-in simulations, finding it varies by up to a factor of ∼3.5 in the modern universe; however, they find minimal systematic variation within galaxies, except for enhanced values in galactic centres (see their figs 1, 2, and 14). Popping, Somerville & Galametz (2017) trace ϵd in semi-analytic models, and find that it can vary with metallicity by up to a factor of ∼2 at metallicities > 0.5 Z⊙ (with the degree and nature of this variation depending considerably upon the specific model).
However, if ϵd does indeed vary significantly with metallicity within our target galaxies, that will actually increase the amount of variation in κ500 in M 83. The highest metallicities are at the inner regions of the disc, where κ500 is already lowest; if increasing fZ in equation (3) also increases ϵd, then this will drive down κd still further. On the other hand, because the lowest values of κd in M 74 are found in the spiral arms, away from the centre, a correlation of ϵd with metallicity could indeed suppress some variation in κ500 – although M 74 already exhibits a much smaller range in κ500 than M 83.
Theoretical dust models can make specific predictions about how ϵd is expected to vary in different conditions. For instance, the THEMIS model traces how dust populations are expected to change in different interstellar environments, predicting that ϵd will increase monotonically with ISM density by a factor of ∼3.5, from 0.27 in the diffuse ISM (|$n_{\mathrm{ H}} = 10^{3}\, {\rm cm^{3}}$|) to 0.88 in the dense ISM (|$n_{\mathrm{ H}} = 10^{6}\, {\rm cm^{3}}$|), driven by the accretion of gas-phase metals on to grains (Jones 2018). We explore the potential effects of this in detail in Section 6.2.1, where we find that it would further increase the variation in κ500.
There are several observational studies that report variation of ϵd between and within galaxies, inferred from the fact that the gas-to-dust ratio is found to vary with metallicity (Rémy-Ruyer et al. 2014; Chiang et al. 2018; De Vis et al. 2019). However, these studies all rely upon an assumed value of κd to infer dust masses, and hence ϵd. Given that we, conversely, use an assumed ϵd to infer κd, it is not really possible to compare such results with ours in a valid way. However, we note with interest that these studies tend to find much larger ranges of ϵd than are suggested by either depletions, simulations, or theoretical dust models – up to 1 dex of scatter at a given metallicity, with up to 3 dex total range over all metallicities. One way to explain this discrepancy would be if κd is depressed at lower metallicity (which is potentially hinted at for M 74 in Fig. 16).
Beside a breakdown in our assumption of a fixed ϵd, it is possible that our method is being corrupted by the presence of ‘dark gas’ – H2 at intermediate densities that CO fails to trace (Reach, Koo & Heiles 1994; Grenier, Casandjian & Terrier 2005; Wolfire et al. 2010). The presence of dark gas would have the effect of causing us to underestimate the value of |$M_{\rm H_{2}}$| in equation (3), thereby artificially driving up κ500. The elevated areas of κ500 in our maps are indeed mainly associated with the interarm regions, where the fraction of dark gas is expected to be greatest (Langer et al. 2014; Smith et al. 2014). Estimates of the fraction of galactic gas mass that is dark range from 0 per cent from dust and gas observations in M 31 (Smith et al. 2012), to 30 per cent in theoretical models (Wolfire et al. 2010), to 42 per cent in hydrodynamical simulations of galactic discs (Smith et al. 2014), to 10–60 per cent from Planck observations of the Milky Way (Planck Collaboration XXI 2011), to 6–60 per cent from Milky Way γ-ray absorption studies (Grenier et al. 2005). Even assuming a worst-case scenario of a 60 per cent dark gas fraction for the interarm regions of our target galaxies (an extreme scenario, given that the 60 per cent represents the single largest fraction amongst the wide range of values reported within the Milky Way), dark gas could only reduce the variation in κ500 we find by a factor of 1.7.
In a similar vein, another potential confounder would be systematic variation in αCO. If αCO increases in denser ISM (independent of metallicity, which we account for), then this could counteract the variation in κd we find. However, evidence to date does not indicate that αCO varies systematically in this way (Sandstrom et al. 2013). This of course could be due to the fact that the uncertainty on αCO (and the scatter on the relations used to derive it) is large – however this uncertainty is propagated through our calculations.
In the course of determining κ500 for each pixel, values for the gas-to-dust ratio, G/D, are also generated. We find 176 < G/D < 277 for M 74, and 140 < G/D < 275 for M 83. Note that these are the ratios of total gas mass to dust mass. In the literature, quoted G/D values are often hydrogen to dust ratios (i.e. no factor of ξ is applied to account for the masses of helium and metals); our hydrogen-to-dust ratios, GH/D, are 127 < GH/D < 201 for M 74, and 100 < GH/D < 196. For high-metallicity systems such as of our target galaxies, these are normal values when compared to the literature (Sandstrom et al. 2013; Rémy-Ruyer et al. 2015; De Vis et al. 2017b). Indeed, we neatly reproduce the factor of 2–3 radial variation in GH/D in M 74 reported by Vílchez et al. (2019) and Chiang et al. (2018) over the 8.35–8.60 (|$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$|) metallicity range we sample; although their adoption of fixed κd limits the scope for detailed comparison. None the less we can say that our inferred dust masses yield sensible G/D values, following expected trends.
Overall, we are confident that our finding of an inverse correlation of κ500 is indeed robust against a wide range of changes to the initial assumptions of our method.
6.2 Alternate models
6.2.1 Variable dust-to-metals ratio
As discussed in Section 2, the assumption of a fixed ϵd is a simplification. Observed depletions in nearby portions of the Milky Way’s diffuse ISM indicate that in reality, ϵd increases with column density (Jenkins 2009; Draine et al. 2014; Roman-Duval et al. 2019). However, the form of this relation in extragalactic systems, where only integrated column density data are available, is not well constrained. None the less, we can still explore, in general terms, how such a model would affect the manner in which κd scales. Even if this approach requires more assumptions, it may be more physical than our fiducial model.
We therefore repeated our κ500 mapping, setting ϵd to vary linearly as a function of ΣISM, with ϵd = 0.75 at the point in each galaxy where ΣISM is highest, and ϵd = 0.25 at the point where ΣISM is lowest. This specific choice of relationship is effectively arbitrary, but approximates the trend reported by Chiang et al. (2018) within M 101, whilst also matching the range of ϵd values reported by De Vis et al. (2019) (although both of these sets of ϵd values were calculated with FIR–submm data, using an assumed value of κd, limiting scope for direct comparison).
The κ500 maps produced using the ϵd ∝ ΣISM model are shown in the left-hand panels of Figs 18 and 19, for M 74 and M 83, respectively. The trend of κ500 being depressed in the denser environments of the spiral arms remains. In fact, the anticorrelation between κ500 against ΣISM is even more exaggerated than was the case for our fiducial model, as can be seen in the left-hand panel of Fig. 20. The Kendall rank correlation coefficients for the ϵd ∝ ΣISM results are more strongly negative than those of the fiducial version, being τ = −0.66 for both M 74 and M 83. The range of κ500 values when using the ϵd ∝ ΣISM model increases to a factor 5 in M 74, and to a factor of 20 in M 83.
It appears that our choice of fixed ϵd in our fiducial model actually serves to reduce the variation in κ500, and that the (probably) more physical ϵd ∝ ΣISM model suggests a notably greater range of values. This increases our confidence that the variation in κ500 we see is a real effect. Whilst we could, for instance, construct a model where ϵddecreases with ISM density by a factor of > 5.3, this would be completely unphysical, and would represent an entirely contrived attempt to minimize the κd variation we find. Similarly, we could construct a model where ϵdincreases with radius – but whilst this would decrease the κd variation in M 83, it would increase it for M 74 (and would again be an unphysical contrivance).
6.2.2 ‘Toy’ model
To establish the degree to which our results might simply be an artefact of our method, we again repeated our κ500 mapping, using a ‘toy’ model. For this repeat, metallicity was fixed at the Solar value of |$12 + {\rm log}_{10} [\frac{\rm O}{\rm H}] = 8.69$|, αCO was fixed at the standard Milky Way value of |$3.2\, {\rm K^{-1}\, km^{-1}\, s\, pc^{-2}}$|, r2:1 was fixed at the local-Universe average of 0.7, Td was fixed at 20 K, β was fixed at 2, and ϵd was fixed at 0.4. Although this toy model is unphysical, it strips out as many assumptions as possible – allowing us to be confident that any trends that persist are not due to our GPR metallicity mapping, our SED fitting, our r2:1 prescription, etc.
The κ500 maps produced using the toy model are shown in the right-hand panels of Figs 18 and 19, for M 74 and M 83, respectively. The corresponding plot of κ500 against ΣISM is shown in the right-hand panel of Fig. 20, where it can be seen that the scatter is markedly increased for both galaxies. For M 83, the trend is none the less still present, with a Kendall rank correlation test giving |$\mathcal {P}_{\rm null} \lt 10^{-5}$|; the lowest values of κ500 are still visibly associated with the largest values of ΣISM, and vice-a-versa. For M 74, the correlation of κ500 with ΣISM is lost; however the far smaller dynamic range in ISM density for this galaxy made it more susceptible to the trend being removed by the toy model’s increase in scatter. The fact the trend with ISM density persists for M 83 despite the use of the toy model is extremely informative. It implies that the basic negative correlation is being driven by the interplay between the 21 cm data, CO data, and 500 |$\mu$|m data – not by the specifics of our method.
6.2.3 Other alternate models
To provide further methodological checks, we produced additional alternate κ500 maps. In Appendix F, we present κ500 maps generated using metallicities calculated via different strong-line prescriptions than the one employed for our fiducial κ500 maps. In Appendix G, we present κ500 maps generated fitting a two-component MBB model to the FIR–submm fluxes, as opposed to the one-component MBB model used for our fiducial κ500 maps. In all cases the resulting κ500 maps display the same general morphology as our fiducial ones, with lower values of κ500 associated with denser regions.
6.3 Implications of findings
Our finding that κ500 shows a strong negative correlation with ISM density is in direct contradiction to standard models of dust emission, which predict that the densest regions of the ISM should exhibit the highest values of κd (Ossenkopf & Henning 1994; Li & Lunine 2003; Jones 2018). This expectation arises from the fact that dust grains in the densest parts of the ISM are predicted to be larger, due to the coagulation of grains and the growth of (icy) mantles on their surfaces, and that larger grains should be more emissive per unit mass (Köhler et al. 2012; Jones et al. 2013; Ysard et al. 2018). The apparent incompatibility of our results with these predictions presents one of two possibilities.
The first possibility is that our method has some fundamental flaw that has systematically affected the results. We have made an effort to construct our method so that it only relies upon standard, widely used assumptions. If one (or more) of these assumptions breaks down systematically, in a manner such that the bias is a function of ISM density, and the bias is a factor of > 5, then this could give rise to the results we see. We have tried to inoculate our findings against even this scenario (for instance, by trying the toy model where all possible variables were kept fixed). However if, for example, dark gas represents 75 per cent of the total gas mass in interarm space (artificially suppressing our assumed |$M_{\rm H_{2}}$|), this could negate our results for M 74. If dark gas represents 75 per cent of the total gas mass in interarm space if and if H ii region oxygen depletion were a factor of 2 lower in interarm space versus other metals (compromising its use as a metallicity tracer), then our results for M 83 could be negated – however scenarios this extreme are unlikely, being unprecedented in the literature, and would have significant implications for extragalactic studies in general.
The second possibility is that κ500 truly does decrease in denser ISM. This would help explain some observational results. For instance, an excess in submm emission has been found in lower density areas within galaxies (Relaño et al. 2018), and within galaxies dominated by diffuse regions (Lamperti et al. accepted; De Looze et al. in preparation); if κ500 is indeed elevated in low-density regions, it could give rise to this effect.
It is hard to explain decreasing κ500 in denser ISM in the context of current dust physics. However it is possible to construct scenarios where it is not entirely unreasonable. Ysard et al. (2018) present a detailed exploration of how changes in various physical parameters of dust should affect κd. For instance, spherical grains are predicted to have lower κ500 than oblate or prolate grains, by up to a factor of ∼1.5 (see their fig. 5); and hydrogenated amorphous carbon grains are expected to have much lower κ500 than amorphous silicate or unhydrogenated amorphous carbon grains (by up to an order magnitude). Whilst we do not suggest that this (or any other) specific physical scenario is the cause of our observed trend, it demonstrates that it is at least possible to envisage an evolution in dust properties that does not entail an uninterrupted monotonic increase in κd with ISM density.
On that theme, we also note that our data only provides physical resolution of 590 pc pix−1 in M 74 and 330 pc pix−1 in M 83. As such, we can do no better than distinguish between arm and interarm pixels. This will have ‘smeared out’ the properties of the denser clouds within the spiral arms. When studies discuss grain growth in the dense ISM, and the associated increases in κd, the dense medium in question is typically described as having at least 1500–10000 |$n_{\rm H}\, {\rm cm^{-3}}$|, compared to 20–50 |$n_{\rm H}\, {\rm cm^{-3}}$| in the diffuse ISM (Ferrière 2001; Köhler et al. 2015; Jones 2018) – a difference in density of at least a factor of 30. However, our data only traces a dynamic range in density of a factor of 5 in M 74, and a factor of 50 in M 83 (discounting pixels within 1 beam of the nuclear starburst, where our κ500 values become unreliable, as per Section 5). So whilst we are probing a wide range of ISM conditions, we are unable to perform a ‘clean’ sampling of the densest grain-growth environments. Likewise, we only performed our analysis for pixels with sufficient SNR for all data – thereby excluding regions of particularly low ISM density, especially at the outskirts of the target galaxies. As such, it seems likely that, in practice, we are effective probing intermediate density environments.
Some grain models do indeed predict that κ500 should drop at intermediate densities, before increasing again at the highest densities. For example, the Köhler et al. (2015) description of the THEMIS model finds that the grain-mixture average κ500 should fall by a factor of 2.3 (relative to the diffuse ISM) for grains undergoing accretion at intermediate densities (|$1500\, n_{\rm H}\, {\rm cm^{-3}}$|) – with κ500 falling by a factor of 26 for amorphous carbon grains in particular. Then at even higher densities, as grains start to aggregate, κ500 will increase again, becoming even higher in the densest regions where icy mantles can form. Again, we do not argue that these specific effects are what are responsible for the relationship we find (as we lack the density resolution, and volume density information, necessary to test this). However, THEMIS does demonstrate that it is possible to construct a physical dust framework where κ500 falls as ΣISM increases, over some intermediate transition regime.
It is also worth considering why our distribution of κ500 values for M 74 is offset from that of M 83, by about 0.3 dex. The most obvious difference in the properties of the two galaxies is the greater ISM surface density of M 83; but given the apparent anticorrelation of κ500 with ΣISM, this seems unlikely to be the driver of the in κ500. M 83 has almost 3 times the SFR of M 73 (Nersesian et al. 2019), despite being physically more compact (see Table 1), giving it an average SFR surface density that is > 6 times greater. Despite this, M 74 has bluer colours, and the relative scale lengths of the dust and stars in M 74 and M 83, as reported in Casasola et al. (2017), differ considerably – in M 74, the dust and gas have very different scale lengths (2.35 arcmin versus 1.04 arcmin), whereas in M 83, the dust and gas scale lengths are effectively identical (1.66 arcmin versus 1.68 arcmin). So there is clearly a difference in the relative geometries of the stars and ISM in these galaxies. When comparing resolved observations of spiral galaxies, it is well established that there can be appreciable differences in ISM properties, even at a given surface density (Usero et al. 2015; Gallagher et al. 2018; Sun et al. 2018). Therefore, it is not necessarily surprising that κd may also have different values in different galaxies, at a given surface density.
X-ray observations of M 74 and M 83 indicate that their interstellar media contain diffuse hot gas components (Owen & Warwick 2009) that span much of their discs. Such gas could process the dust in a galaxy, sputtering the grains, and (in standard models) therefore decreasing the grains’ κd (Galliano et al. 2018). That said, we find decreased κd in the denser ISM, where grains should be more shielded from X-ray gas. None the less, it is possible that the trends we find may not be applicable to galaxies with less prominent X-ray gas content.
7 CONCLUSION
Using a homogenous data set assembled as part of the DustPedia project (Davies et al. 2017), we have produced the first maps of the dust mass absorption coefficient, κd, within two nearby galaxies: M 74 (NGC 628) and M 83 (NGC 5236).
Our method for finding κd is empirical, and avoids making any assumptions about the composition or radiative properties of the dust. Instead, our approach exploits the fact that the ISM dust-to-metals ratio seems to exhibit minimal variation at high metallicity. With this one assumption, we can use gas and metallicity data to determine dust masses a priori; by comparing these masses to observed dust emission, we are able to calibrate values for κd. Given that the value of the dust-to-metals ratio is much less uncertain than the value of κd, we are able to leverage the one to explore the other.
As a proof-of-concept demonstration, we have applied this method on a resolved, pixel-by-pixel basis to M 74 and M 83, two nearby face-on spiral galaxies, that have well-suited atomic gas, molecular gas, dust emission, and ISM metallicity data available. We have produced gas-phase metallicity maps for these galaxies, using the many hundreds of available spectra measurements, via a novel application of Gaussian process regression, with which we infer the underlying metallicity distribution.
We find strong evidence for significant variation in κ500 within both galaxies – by a factor of 2.3 within M 74 (0.11–0.25 |${\rm m^{2}\, kg^{-1}}$|), and by a factor of 5.3 within M 83 (0.15–0.80 |${\rm m^{2}\, kg^{-1}}$|).
We examine whether κd shows variation with other measured and derived properties of the target galaxies. We find that κd exhibits a distinct negative correlation with the surface density of the ISM, following a power-law slope of index |$-0.36^{+0.26}_{-0.21}$| (although the power laws for the two galaxies are offset by 0.3 dex). This trend appears to be dictated by the total ISM surface density, as opposed to the surface density of either its atomic, molecular, or dust components. This trend is the opposite of what is predicted by most dust models. However, the relationship is robust against a wide range of changes to our method – only the adoption of unphysical or highly unusual assumptions would be able to suppress it. We discuss possible ways of reconciling this finding with the current understanding of dust physics – such as the possibility that our combination of resolution and sensitivity means that we biased towards probing regimes of intermediate density where the broader expected correlation between density and κd may not hold true.
We also find tentative indications of correlation of κd with other properties, such as metallicity, NIR radiation field intensity, and dust emissivity slope β. However, the evidence for these is less conclusive (and some of these parameters were inputs to our κd calculations), so we are more cautious about the significance of these relationships.
This study lays the groundwork for a wide range of future work. An expanded study of resolved κd is possible with the DustPedia data set, but at present the availability of well-resolved metallicity data would limit it to a sample of only 10–20 galaxies. But in future, large IFU surveys of highly extended nearby galaxies, especially the SDSS-V Local Volume Mapper (Kollmeier et al. 2017) will dramatically improve this situation. Simultaneously, data now exists to apply the dust-to-metals method to large, statistical samples of galaxies on a global basis; in particular, the Jcmt dust and gas In Nearby Galaxies Legacy Exploration (JINGLE; Saintonge et al. 2018), which is assembling consistent high-quality CO, H i, dust, and IFU data for almost 200 galaxies, would be well suited to this task.
Most importantly, many of the questions raised could be tackled by conducting a similar analysis at improved spatial resolution. For this reason, we have begun work on applying this method as part of an analysis of several Local Group galaxies – including the Large and Small Magellanic Clouds, where we enjoy particularly exquisite resolution. Most significantly, better resolution will allow us to cleanly probe a larger range of density, from the densest grain-grown regions, down to the most diffuse ISM. We will thereby test if the surprising anticorrelation between κ500 and ΣISM holds true. Another benefit to expanding our analysis to the Magellanic Clouds is that they are the subjects of ongoing work to perform the first extragalactic depletion analyses (Jenkins & Wallerstein 2017; Roman-Duval et al. 2019). Exploiting that data will allow us to use in situ measurements of the dust-to-metal ratio, removing the single largest source of uncertainty we presently face, and allowing us to produce the most reliable empirical κd determinations available with current data.
ACKNOWLEDGEMENTS
The DustPedia project22 (Davies et al. 2017) has received funding from the European Union’s Seventh Framework Programme (FP7) for research, technological development, and demonstration, under grant agreement 606824 (PI Jon Davies).
The authors thank the anonymous referee whose comments have materially improved the quality of this work.
CJRC acknowledges financial support from the National Aeronautics and Space Administration (NASA) Astrophysics Data Analysis Program (ADAP) grant 80NSSC18K0944. CJRC thanks Andreas Lundgren and Tommy Wiklind for providing reduced SEST CO data for M 83 (Lundgren et al. 2004). CJRC also thanks Philip Wiseman, Bruce Draine, Julia Roman-Duval, Karl Gordon, Rosie Beeston, and Phil Cigan for helpful discussions and input.
This research made use of astropy,23 a community-developed core python package for Astronomy (astropy Collaboration 2013, 2018). This research made use of astroquery,24 an astropy-affiliated python package for accessing remotely hosted astronomical data (Ginsburg et al. 2019). This research made use of reproject,25 an astropy-affiliated python package for image reprojection. This research has made use of numpy26 (van der Walt, Colbert & Varoquaux 2011), scipy27 (Jones et al. 2001), and matplotlib28 (Hunter 2007). This research made use of aplpy,29 an open-source plotting package for python (Robitaille & Bressert 2012). This research made use of the pandas30 data structures package for python (McKinney 2010). This research made use of the scikit-image31 image processing package for python and the scikit-learn32 machine learning package for python (Pedregosa et al. 2011). This research made use of emcee,33 the MCMC hammer for python (Foreman-Mackey et al. 2013). This research made use of the pymc334 MCMC package for python. This research made use of the corner35 scatterplot matrix plotting package for python (Foreman-Mackey 2016). This research made use of ipython, an enhanced interactive Python (Pérez & Granger 2007). This research made use of python code for working in the luminance-chroma-hue colour space, written by Endolith,36 kindly made available free and open-source under the BSD License,37 and copyright 2014 Endlolith.
This research has made use of topcat38 (Taylor 2005), an interactive graphical viewer and editor for tabular data, which was initially developed under the UK Starlink project, and has since been supported by the Particle Physics and Astronomy Research Council (PPARC), the VOTech project, the AstroGrid project, the Astronomical Infrastructure for Data Access (AIDA) project, the Science and Technology Facilities Council (STFC), the German Astrophysical Virtual Observatory (GAVO) project, the European Space Agency (ESA), and the Gaia European Network for Improved data User Services (GENIUS) project. This research made use of ds9, a tool for data visualization supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEASARC) with support from the James Webb Space Telescope (JWST) Mission office at the Space Telescope Science Institute for 3D visualization.
This research made use of montage,39 which is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the NASA Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology.
This research made use of the VizieR catalogue access tool40 (Ochsenbein, Bauer & Marcout 2000), operated at CDS, Strasbourg, France. This research has made use of the Nasa/ipac Extragalactic Data base41 (NED), operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.
This research made use of data from the Swedish-Eso Submillimetre Telescope, which was operated jointly by the European Southern Observatory (ESO) and the Swedish National Facility for Radio Astronomy, Chalmers University of Technology.
Much of the model fitting performed in this work benefitted from the invaluable guidance provided in Hogg, Bovy & Lang (2010).
Footnotes
As defined according to detection by the Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010), at its all-sky sensitivity, in 3.4 |$\mu$|m (its most sensitive band).
Whilst 24 |$\mu$|m Spitzer data does exist for these galaxies, the background is better behaved in the WISE data, due to the superior mosaicking permitted by the larger field of view.
Centred at approximately: α = 204.20°, δ = −29.87°.
SPIRE Instrument & Calibration Wiki: https://herschel.esac.esa.int/twiki/bin/view/Public/SpireCalibrationWeb
Standardized to allow modes and widths other than zero, as per the scipy (Jones et al. 2001) implementation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html
WISE All-Sky Release Explanatory Supplement (Cutri et al. 2012): https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec4_4h.html
PACS Instrument & Calibration Wiki: https://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb
2.3 per cent being the non-correlated component of the Herschel-SPIRE calibration uncertainty, separate from υSPIRE.
Comparing the means of the last 90–100 per cent quantile of the combined chains to the 50–60 per cent quantile.
WISE colour corrections from Wright et al. (2010).
Spitzer-MIPS colour corrections from the MIPS Instrument Handbook, version 3 (Colbert 2011): https://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/51/#_Toc288032329
Herschel-PACS colour corrections from the PACS Handbook, version 4.0.1 (Exter et al. 2019): https://www.cosmos.esa.int/documents/12133/996891/PACS+Explanatory + Supplement
Herschel-SPIRE colour corrections from the SPIRE Handbook, version 3.1 (Valtchanov et al. 2017): https://herschel.esac.esa.int/Docs/SPIRE/spire_handbook.pdf
Maps were reprojected to the same pixel grid as the κ500 maps, then background subtracted in the same manner as the continuum maps in Section 4.1. We manually masked pixels containing obvious foreground Milky Way stars.
Spearman and Pearson rank correlation tests similarly both give |$\mathcal {P}_{\rm null} \lt 0.025$|, with correlation coefficients > 0.2.
By ‘direct’, we refer to those measurements where ϵd is determined from observing the mass fraction of metals depleted from the gas phase.
https://docs.pymc.io/; Salvatier, Wiecki & Fonnesbeck 2016
Standardized as per the scipy (Jones et al. 2001) gamma distribution implementation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html
REFERENCES
APPENDIX A: RADIAL METALLICITY PROFILE FITTING
We determined the posterior probability of our variables of interest – mZ, cZ, and ψ – in a Bayesian manner, sampling the posterior PDF using the pymc3 (Salvatier et al. 2016) MCMC package for python.
To inform the priors, we first performed a simple, preliminary least-squares fit, with only the gradient and central metallicity as free parameters. The priors on all three parameters then took the form of normal distributions. For cZ, the mean of the prior was set to the central metallicity found by the preliminary least-squares fit, and the standard deviation on the prior was set to the standard deviation of all the input metallicity values. For mZ, the mean of the prior was set to the gradient found by the preliminary least-squares fit, and standard deviation of the prior was set to the absolute value of the gradient found by the preliminary least-squares fit. For ψ, both the mean and standard deviation of the prior were set to the root-mean-square of the residuals between the input metallicity values and the preliminary least-squares fit.
APPENDIX B: UNCERTAINTIES ON GPR METALLICITY MAPPING
To determine the uncertainty of the GPR metallicity maps, we repeated the regression procedure 1000 times. For each iteration, we draw a random sample from the posterior PDF of our radial metallicity profile model, and used that sample to calculate the residual on each data point; we then applied the GPR to these residuals in the same manner as described above. For each iteration, the GPR produced a full posterior PDF for the predicted metallicity in each pixel (and by definition, Gaussian process regression yields Gaussian posterior PDFs).
Having repeated this process for the 1000 iterations, we had 1000 posterior PDFs for each pixel; these are then combined to give each pixel’s final metallicity PDF. To quantify the uncertainty in each pixel, we take the 63.3 per cent quantile around the posterior median; these are the uncertainty values plotted in the lower left panels of Figs 5 and 6. As can be seen, the uncertainty on the regression is low (<0.05 dex) for pixels that have plenty of spectra metallicities; whilst for pixels more distant from any spectra, making the predicted values more dependent upon extrapolation, the uncertainty is much larger (> 0.25 dex). Indeed, for pixels with few or no spectra metallicities in the immediate vicinity, relying upon the metallicity predicted by a one-dimensional globally fitted gradient could provide a false sense of confidence – especially in M 83, where the metallicity data are concentrated in a central band. We therefore argue that in these areas, the larger uncertainties predicted by our GPR approach are likely to be more realistic.
APPENDIX C: VALIDATION OF GPR METALLICITY MAPPING
To verify that our GPR metallicity mapping technique is reliable, and not generating spurious features in the final metallicity maps, we used a Monte Carlo jackknife cross-validation analysis. For this, we performed 500 repeats of the GPR metallicity mapping; for each repeat, half of the spectra metallicity points were selected at random to be excluded from the fitting, to serve as a control sample for later reference. The GPR was then computed using the remaining half of the points (but otherwise following the modelling process as laid out above). By comparing the metallicity values of the masked spectra to the metallicities predicted by the GPR method at their positions, we can evaluate the accuracy of the generated metallicity maps.
If the metallicities predicted via GPR suffer from no systematic offset, then the mean χ should be |$0 \pm n^{-\frac{1}{2}}$| (where n is the number of control spectra). Similarly, if the uncertainties on the GPR metallicities are Gaussian and accurate, then 68.3 per cent of the values of χ should lie in the range −1 < χ < 1.
The distribution of jackknife χ values we find for both galaxies are shown in Fig. C1. The distributions are symmetric, near-Gaussian, and centred close to zero. The mean jackknife χ values are 0.0079 ± 0.0029 and −0.0057 ± 0.0023 for M 74 and M 83, respectively. These offsets are >2σ, suggesting that there tends to be a small systematic offset (positive for M 74, and negative for M 83) between the metallicity predicted by the GPR, and the actual metallicity of the spectra. But whilst technically significant, these systematic offsets are none the less vanishingly small in terms of actual metallicity – the mean jackknife deviation in |$12 + {\rm log}_{10} [ \frac{\rm O}{\rm H}]$| units is 0.00056 for M 74, and −0.00037 for M 83. We are satisfied that systematic effects at this scale are minute enough to have no appreciable impact on any of our results.
For M 74, 80.5 per cent of the jackknife χ values lie in the −1 < χ < 1 range; for M 83 the fraction is 85.4 per cent. These are both somewhat larger than the expectation of 68.3 per cent, which suggests that our GPR maps are actually somewhat more precise than suggested by their uncertainties. In other words, it appears that the GPR uncertainties are overestimated by factors of approximately 1.18 and 1.25 (for M 74 and M 83, respectively) – a small enough difference that we judge it unnecessary to attempt a post-hoc fine tuning of the output uncertainties.
Additionally, see the discussion in Section 6.1 of the effect of the metallicity maps upon our resulting maps of κd.
APPENDIX D: DUST SED PRIORS
Our SED-fitting procedure, detailed in Section 4.2, has six free parameters: dust temperature, Td; dust ‘mass’ normalization, |$M_{\mathrm{ d}}^{\rm (norm)}$|; emissivity slope, β; and correlated photometric error in the Herschel-SPIRE bands, υSPIRE. The prior probability distributions for all these free parameters are shown in Fig. D1.
D1 Temperature prior
For Td, the modal value of 20 K corresponds to the approximate average of the cold dust temperatures seen in nearby galaxies, in both global (Galametz et al. 2012; Clemens et al. 2013; Ciesla et al. 2014) and resolved (Smith et al. 2012; Gordon et al. 2014; Tabatabaei et al. 2014) analyses. Across the 13–30 K temperature range, |$\mathcal {P}(T_{\mathrm{ c}}) \gt 0.8\mathcal {P}(\check{T}_{\mathrm{ c}})$|; this corresponds to range spanned by the lower cold dust temperatures seen in blue dust- and gas-rich galaxies (Clark et al. 2015; Dunne et al. 2018), to the higher cold dust temperatures seen in dust-poor dwarf galaxies (Rémy-Ruyer et al. 2013; Izotov et al. 2014). In other words, temperatures across this ‘standard’ temperature range are only slightly less favoured than |$\check{T}_{\mathrm{ c}}$|. Outside this range, there is an increasing penalty – especially towards lower temperatures, where |$\mathcal {P}(T_{\mathrm{ c}}\lt 5) = 0$|, to rule out unphysically cold dust.
D2 Mass prior
As described in Section 4.2, our SED-fitting procedure uses an arbitrary placeholder value of κd (because the whole purpose of the SED fitting is to find values of κd); as a result, the ‘mass’ variable being fitted simply serves as a normalization parameter.
Equation (D3) is a purely empirical relation, derived from the SED fitting of Herschel Reference Survey (Boselli et al. 2010) galaxies, as performed in Clark et al. (2015) and Clark et al. (2016).
The prior on M(norm) shown in Fig. D1 is for an example pixel from our M 83 data (processed as per Section 4.1), centred at α = 204.2841°, δ = −29.8559°. The brightest band for this pixel is 160 |$\mu$|m, where the flux is 1.18 Jy. Given a distance to M 83 of 4.9 Mpc, that corresponds to a priors centred at |$\check{M}^{\rm (norm)} = 5.45\, {\rm log_{10}\, M_{\odot }}$|, as per equation (D3).
The mass normalization prior is designed to be very weak. This is because the strong M ∝ T4 + β dependence of mass on temperature (for a given luminosity) means that the fitted value of the mass normalization term is often driven primarily by the fitted temperature.
D3 β prior
The prior on β takes the form of a standardized gamma distribution, identical to equations (D1) and (D2), except with β replacing Td, and |$\check{\beta }$| replacing |$\check{T}$|. The parameters for our β prior take values of α = 2.75, l = 0, and |$\check{\beta } = 1.75$|.
For nearby galaxies and the Milky Way, β is typically found to lie in the range 1.5–2.0, with resolved analyses finding values spanning 1.0–2.75 (Smith et al. 2012; Kirkpatrick et al. 2013; Planck Collaboration XI 2014). We therefore construct our prior such that it peaks at |$\check{\beta } = 1.75$|, with |$\mathcal {P}(\beta) \gt 0.8\mathcal {P}(\check{\beta })$| across the 1.0–2.75 range. To exclude dubiously physical low β values, |$\mathcal {P}(\beta \lt 0) = 0$|.
D4 υSPIRE prior
As discussed in Section 4.2, the calibration uncertainties on Herschel-SPIRE photometry has a correlated systematic error component, which we term υSPIRE, arising from uncertainty on the emission model of Neptune, the instrument’s primary calibrator. υSPIRE has a value of |$\pm 4{{\ \rm per\ cent}}$|; the true value of the systemic error is believed to be equally likely to lie anywhere in that range, with minimal likelihood (|${\sim}5{{\ \rm per\ cent}}$|) of the value lying outside it (Bendo et al. 2013; A. Papageorgiou, private communication; C. North, private communication).
We therefore use a prior for υSPIRE that takes the form of a boxcar function convolved with a Gaussian distribution. The boxcar function has a value of 1 over the range −0.04 to 0.04, with a value of 0 beyond this. The Gaussian with which it was smoothed has a standard deviation of 0.005. In the resulting prior, as shown in the lower right panel of Fig. D1, 95 per cent of the probability density is contained within the −0.04 to 0.04 range.
APPENDIX E: κd MAPS AT 160 |$\mu$|m
As discussed in Section 5, we calculated our κd maps at a reference wavelength of 500 |$\mu$|m. This is the longest wavelength at which we have data, making it less sensitive to uncertainties in temperature derived from the SED fitting. However, many authors opt to present κd at 160 |$\mu$|m, as this is the wavelength regime at which the κd of carbonaceous and silicate dust is most comparable.
For completeness, we therefore also produced κ160 maps, which are shown in Fig. E1. The maps are noisier than those computed at 500 |$\mu$|m, but the overall morphology of κ160 in both galaxies is none the less the same as that of κ500.
Using the same independent-pixel non-parametric bootstrap approach as in Section 5, we find a median underlying range of κ160 values of 0.74–2.4 |${\rm m^{2}\, kg^{-1}}$| for M 74 (a factor of 3.2 variation), and 2.1–12 |${\rm m^{2}\, kg^{-1}}$| for M 83 (a factor of 5.7 variation).
APPENDIX F: κ500 MAPS USING DIFFERENT STRONG-LINE METALLICITY PRESCRIPTIONS
As described in Section 3.3, our metallicity maps were produced using metallicities calculated using the ‘S’ strong-line prescription of Pilyugin & Grebel (2016). To ensure that our specific choice of metallicity prescription was not driving our results, we also repeated our κ500 mapping using metallicity maps produced using four other strong-line prescriptions; the O3N2 prescription of Pettini & Pagel (2004), the N2 prescription of Pettini & Pagel (2004), the prescription of Tremonti et al. (2004), and the IZI prescription of Blanc et al. (2015). As with our fiducial Pilyugin & Grebel (2016) ‘S’ prescription values, these metallicities are all taken from the standardized data base produced by De Vis et al. (2019). The resulting κ500 maps for all four prescriptions for both galaxies are presented in Figs F1 and F2. These κ500 maps all display the same general morphology as the fiducial maps in Figs 12 and 13 – with lower values of κ500 associated with regions of denser ISM. The exceptions to this are the maps produced using the Tremonti et al. (2004) prescription, which causes a negative radial gradient in κ500 to dominate over the density-anticorrelated variations; but none the less, at a given radius, areas of lowest κ500 are associated with the same areas of denser ISM as seen in the other maps.
APPENDIX G: κd MAPS FROM TWO-COMPONENT MBB SEDS
where subscripts c and w denote the cold and warm dust components, respectively. There are therefore six free parameters for the two-component MBB modelling: Tc, |$M_{\mathrm{ c}}^{\rm (norm)}$|, Tw, |$M_{w}^{\rm (norm)}$|, β, and υSPIRE. Having performed this SED fitting, computing the corresponding values of κ500 simply requires setting n = 2 in equation (3) and providing T and Sλ for both MBB components.
Expanding the method to incorporate two dust components includes the tacit assumption that both dust components have the same dust-to-metals ratio. This is perhaps unlikely, as warmer dust will generally be associated with recent star formation and more intense ISRFs, where shocks and high-energy photons might destroy grains and return their metals to the gas phase. However, for dust SEDs with two distinct components at different temperatures, the total dust mass is invariably dominated by the colder component (da Cunha et al. 2008; Kirkpatrick et al. 2014; Clark et al. 2015); therefore the resulting value of κd will primarily reflect the κd of the dominant component, insulating this approach against differences in ϵd.
The resulting maps are shown in Fig. G1. The map for M 74 shows some increase in κ500 in the centre relative to the one-component MBB approach, whilst the map for M 83 is practically identical.
The median for M 74 is κ500 = 0.20 |${\rm m^{2}\, kg^{-1}}$|, and the median for M 83 is κ500 = 0.25 |${\rm m^{2}\, kg^{-1}}$|. The ranges of values (estimated via same the non-parametric independent-pixel bootstrap method as used in Section 5) are 0.13–0.28 |${\rm m^{2}\, kg^{-1}}$| for M 74 (a factor of 2.2 variation), and 0.12–0.72 |${\rm m^{2}\, kg^{-1}}$| for M 83 (a factor of 6.0 variation). The differences between these values and their counterparts for our fiducial maps are all much less than the average 0.15 dex statistical uncertainty on each pixel’s κ500 value (and well within the 0.2 dex systematic uncertainty).