The Quest for the Missing Dust. II. Two Orders of Magnitude of Evolution in the Dust-to-gas Ratio Resolved within Local Group Galaxies

The relationship between dust and gas in the interstellar medium (ISM) is complex. There are various mechanisms by which dust is added to the ISM, and others that remove dust from it.

For instance, newly created dust grains can be added to the ISM via the deaths of stars—either by core-collapse supernovæ (Barlow et al. 2010; Matsuura et al. 2011; Gomez et al. 2012), or by asymptotic giant branch (AGB) stars (Höfner & Olofsson 2018). The dust mass of the ISM can also increase in situ, through dust grains accreting gas-phase metals in higher-density environments (Köhler et al. 2015; Zhukovska et al. 2016; Jones et al. 2017). Similarly, dust can be removed from the ISM in a number of different ways. Dust can undergo photodestruction by high-energy photons, especially in poorly shielded low-density environments; or dust can be sputtered by forward and reverse supernovæ shocks (making the status of supernovæ as net sources or sinks of dust uncertain; Jones 2004; Bocchio et al. 2014; Slavin et al. 2015).

The key metric for assessing the dust-richness of the ISM is the dust-to-gas ratio (D/G). The various processes that add or remove dust from the ISM will change the D/G, as will other processes, such as the preferential removal of the dust-richest ISM via the formation of new stars and planets (Hjorth et al. 2014; Forgan et al. 2017), and the dilution of the ISM by the infall of low-metallicity extragalactic gas (Edmunds & Eales 1998)—both of which will drive down the D/G.

Clearly, the balance between dust creation and destruction in galaxies, and therefore D/G, is sensitive to many factors. And unavoidably, the relative importance of these factors will vary greatly between different environments, both between, and within, galaxies. The balance between dust creation and destruction is especially important in light of the "dust budget crisis," arising from the fact that observed sources of dust appear insufficient to account for observed dust masses in some galaxies (Matsuura et al. 2009), particularly at high redshift (Rowlands et al. 2014), unless significant grain growth occurs in the ISM (De Vis et al. 2017b; Galliano et al. 2021).

Strong evidence for grain growth comes from the fact that the D/G is observed to increase in higher-density environments, even when controlling for metallicity. This has been seen in spectroscopic absorption measurements of the depletion of various metals from the gas phase (Jenkins 2009; Tchernyshyov et al. 2015; Roman-Duval et al. 2021), in measurements of extinction curve evolution (Gordon et al. 2003; Fitzpatrick & Massa 2005), and in resolved far-IR (FIR) observations of the dust emission within nearby galaxies (Mattsson et al. 2014; Roman-Duval et al. 2014, 2017).

It should be noted that work by Nanni et al. (2020) found that grain growth is not necessarily required to explain the chemical evolution of low-metallicity and high-z galaxies, if the initial mass function is top-heavy, and if outflow rates are high. Similarly, work by De Looze et al. (2020) indicated that low-redshift dust masses could be explained primarily by dust from stellar sources, without the need for efficient interstellar grain growth, if dust lifetimes are long (1–2 Gyr) and if the survival fraction of fresh supernova dust passing through supernova reverse shocks is high (37%–89%). However, even this framework requires that 20%–50% of the present-day dust mass of galaxies has grown in the ISM. These works highlight the necessity of constraining the role and relative importance of grain growth in galaxies' chemical evolution.

In recent years, one of the key avenues for understanding the chemical evolution of galaxies has been the relationship between D/G and metallicity ⁵  (see Figure 9 in Galliano et al. 2018; see also Rémy-Ruyer et al. 2014 and De Vis et al. 2019). Understanding how the D/G evolves with metallicity is especially important with respect to the "critical metallicity"—the knee in the D/G–Z relationship, at approximately 0.2 Z_⊙ (Rémy-Ruyer et al. 2014; De Vis et al. 2019). Above the critical metallicity, the D/G increases linearly with Z, while below it, the D/G drops sharply. A linear trend in the D/G–Z is expected in a regime where the dust-to-metals ratio remains constant, or in other words, where a constant fraction of metals is locked up in dust grains.

The critical metallicity is generally understood to reflect a threshold above which ISM grain growth becomes efficient, and starts to dominate over less-efficient stellar sources of dust, and is expected from theoretical models (Asano et al. 2013; Feldmann 2015; Zhukovska et al. 2016). Understanding the critical metallicity has important ramifications for understanding the ISM in general, particularly given the role of dust as the formation site of molecular hydrogen (Gould & Salpeter 1963), as a cooling pathway during star formation (Dopcke et al. 2011), and in shielding CO from photodissociation—thereby dictating the CO-to-H₂ conversion needed for using CO as a tracer of molecular gas (Wolfire et al. 2010; Clark & Glover 2015).

There are, however, discrepancies in our current understanding of the D/G–Z relationship. The critical metallicity break is only observed for (mostly nearby) galaxies for which dust mass measurements are derived from FIR observations, and for which gas masses are measured via H i and CO observations (Rémy-Ruyer et al. 2014; De Vis et al. 2019). In contrast, the critical metallicity break is not observed for galaxies for which D/G is derived from depletions, specifically in damped Lyα systems (DLAs). ⁶  Rather, for DLAs, the D/G–Z trend remains linear, over the entire 0.01 < Z < 1 Z_⊙ range for which measurements exist (Galliano et al. 2018; Roman-Duval et al. 2022a), with no critical metallicity break. Previously, this linear trend in D/G–Z for DLAs was found to overlap the trend for galaxies with D/G determined via FIR observations; however, recent recalibration indicates that the DLA trend is actually offset to lower D/G (Roman-Duval et al. 2022b).

Even within the Local Group (the largest members of which are shown in Figure 1, excepting the Milky Way), problems remain in the D/G–Z relationship. Specifically, the D/G estimates for the Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC) derived from FIR observations are considerably lower than D/G values derived from UV absorption line spectroscopy measurements of depletions (Roman-Duval et al. 2017). The FIR D/G values are smaller than the UV D/G values by a factor of ∼2 for the LMC, and a factor of ∼5 for the SMC. These offsets between D/G are persistent over the full range of gas densities studied within both Magellanic Clouds to date (Roman-Duval et al. 2017, 2022b).

Zoom In Zoom Out Reset image size

Are these D/G discrepancies in the Magellanic Clouds due to previous FIR measurements being affected by temperature mixing, or other systematic problems? Or are the D/G measurements from the UV depletions in error, meaning that the lower D/G values from the FIR might be tracing the beginning of the critical metallicity transition? The LMC and SMC have metallicities of 0.5 and 0.2 Z_⊙, respectively, so the SMC in particular might be in the critical metallicity regime.

There are challenges that need to be addressed when making D/G measurements with FIR data. If the physical resolution of FIR observations is not good enough, it will systematically bias dust mass estimates due to beam smearing. High-density regions get blurred into the lower-density regions that surround them, and temperature mixing can cause higher-luminosity emission from warm dust to dominate the FIR spectral energy distribution (SED) over fainter emission from greater masses of cold dust. Additionally, limited physical resolution also diminishes the range of ISM densities that can be sampled. Moreover, galaxies sufficiently nearby to overcome such resolution limits (like those in the Local Group) often suffer from large-angular-scale emission being filtered out of FIR data during the reduction process, systematically biasing the data against properly sampling dust in diffuse environments, where ISM densities will be lowest. Fortunately, by using high-resolution data and the latest reduction techniques, it is possible to tackle these issues.

The specifics of these obstacles, and how we are able to overcome them for the sample of Local Group galaxies we study in this work, are discussed fully in Section 2. In Section 3, we describe our pixel-by-pixel SED fitting, and the dust results obtained from it. In Section 4, we examine D/G in our target galaxies, and how it evolves with ISM density. In Section 5, we investigate the surprising turnover in D/G exhibited by some of the galaxies in our sample. In Section 6, we discuss how the FIR-derived D/G measurements from this (and previous) work compares to UV-derived D/H measurements. In Section 7, we provide specific details of the various maps being publicly released alongside this paper.

Note that in this paper, we specifically consider D/G in terms of the dust-to-hydrogen ratio, D/H. We do this for several reasons. First, hydrogen mass (or density) is the actual observable quantity. Frequently, authors will use a measured hydrogen mass to infer the total gas mass, by applying a correction factor to account for the faction of the gas mass made up of elements heavier than hydrogen. However, this fails to consider the fact that the mass fraction of elements heavier than hydrogen evolves as a function of ISM metallicity (Balser 2006)—varying from 1.33 (for low-metallicity galaxies where Z → 0) to 1.45 (for high-metallicity giant ellipticals where Z = 1.5 Z_⊙). Moreover, this evolution with metallicity is not purely monotonic, and will have intrinsic scatter due to differences in nucleosynthesis from different stellar populations. Directly working with D/H allows us to sidestep the question of this correction, and the ambiguity it can cause. Additionally, later in the paper we compare our dust-to-gas measurements to values calculated from depletions, where the observed hydrogen column is always used as the denominator; so here, once again, use of D/H allows for more clarity.

Physical resolution is key for investigating D/H, and other ISM properties. Superior physical resolution translates into superior density resolution. With better resolution, compact higher-density regions of ISM will appear less "smeared" together with surrounding lower-density material, allowing a wider range of environments to be probed. Achieving the best physical and density resolution therefore requires observing galaxies that are suitably nearby, using telescopes that have sufficiently high angular resolution.

Physical resolution is especially important for observations of dust emission. Attempts to measure dust masses (and dust mass surface densities) of nearby galaxies generally suffer from the problem that multiple dust components, each with different temperatures (and other properties), will be confused together in the observed FIR SED. This will bias the dust properties inferred from that SED; specifically, this temperature mixing tends to cause dust masses to be underestimated (Galliano et al. 2011; Priestley & Whitworth 2020). The best way to overcome temperature mixing is to observe galaxies using sufficient resolution that the dust in any given resolution element is mostly homogeneous. In practice, Galliano et al. (2011) found that this bias becomes minimal only when physical resolution is better than ∼150 pc. Specifically, Galliano et al. (2011) found that increasing physical resolution leads to increased dust masses from SED fitting, with an asymptote at ≈20 pc (above this, improving the resolution leads to no further change in modeled mass). However, physical resolution better than ≈150 pc results in dust masses for which the error is <10% divergent from masses obtained at the 20 pc asymptote; see Figure 6 in Galliano et al. (2011).

These requirements for high spatial resolution mean that we are effectively limited to observations of galaxies in the Local Group, if we want to obtain D/H estimates free of temperature mixing bias.

2.1. Sample Galaxies

For our investigation into the resolved properties of D/H in nearby galaxies, our sample consists of the four very well-resolved galaxies of the Local Group—the LMC, SMC, M31 (the Andromeda galaxy), and M33 (the Triangulum galaxy). Key properties for these galaxies are provided in Table 1.

Table 1. Basic Properties of the Local Group Galaxies Studied in This Work

	M31	M33	LMC	SMC
α (J2000)	1069	2346	8089	1316
δ (J2000)	+4127	+3066	−6976	−7280
Distance (kpc)	790	840	50	62
Hubble Type	SAb	SAcd	SBm	Irr
R25	89'	32'	323'	151'
R25 (kpc)	20.5	7.5	5.0	2.5
a Pos. Angle	35°	23°	170°	45°
a Axial Ratio	2.57	1.70	1.17	1.66
Inclination	77°	56°	26°	⋯
M⋆ (M⊙)	b 1011.1	c 109.7	d 109.4	e 108.3
f Z (Z⊙)	g 1.3	h 0.5	i 0.5	i 0.2

Notes. Values taken from the NASA/IPAC Extragalactic Database (NED; 2019), except where otherwise noted. Portions of this table are reproduced from Table 1 of Paper I, provided again here for convenience.

^a From the HyperLEDA database (Makarov et al. 2014). ^b Tamm et al. (2012). ^c Corbelli et al. (2014). ^d van der Marel (2006). ^e Rubele et al. (2018). ^f All metallicities are from "direct" method determinations using auroral lines in the ISM. For M31 and M33, we quote the "characteristic metallicity," as measured at 0.4R₂₅ (Pilyugin et al. 2004; Moustakas et al. 2010). ^g Zurita & Bresolin (2012); we have corrected for the −0.3 dex offset they report due to the low-metallicity bias in which Hii regions have detectable auroral lines in M31. ^h Magrini et al. (2016). ⁱ Toribio San Cipriano et al. (2017).

Download table as:  ASCIITypeset image

These galaxies exhibit a broad range of characteristics. The SMC is a 0.2 Z_⊙ (Toribio San Cipriano et al. 2017) dwarf galaxy exhibiting the distinct ISM traits that are characteristic of low-metallicity systems, such as: the absence of the 2175 Å extinction bump along most (but not all) sightlines (Gordon et al. 2003; Murray et al. 2019); being extremely irregular, having likely been disturbed by a recent collision with the LMC (van der Marel & Cioni 2001; Murray et al. 2019; Choi et al. 2022); and being highly elongated along our line of sight (Scowcroft et al. 2016). The LMC has a metallicity 0.5 Z_⊙ (Toribio San Cipriano et al. 2017), and is on the spiral/dwarf-irregular transition; it is experiencing intense star formation, especially on its eastern edge, where gas from the SMC appears to be infalling (Bekki & Chiba 2007; Fukui et al. 2017; Tsuge et al. 2019). M33 has a similar 0.5 Z_⊙ metallicity to the LMC (Koning 2015; Magrini et al. 2016), but has almost twice the stellar mass (van der Marel 2006; Corbelli et al. 2014), and a much more orderly spiral structure; it is an interacting companion of M31 (Bekki 2008; Putman et al. 2009), and has the highest star-forming efficiency in the Local Group (Gardan et al. 2007). M31 is a high-mass, L*, spiral/ring galaxy (Gordon et al. 2006) that reaches supersolar metallicity (Zurita & Bresolin 2012), and is currently passing through the green valley (Mutch et al. 2011).

This wide range of properties gives excellent scope to examine what traits may affect D/H. The sample is illustrated in Figure 1, where the variation in the relative abundance of dust and gas can be clearly seen from visual inspection alone. All four galaxies have sufficiently low radial velocities (−300 < v < 300 km s⁻¹) that redshift/violetshift is entirely negligible, allowing us to safely treat all emission as being rest-frame.

2.2. New Herschel Data

The natural choice of telescope for observing dust in Local Group galaxies is the Herschel Space Observatory (Pilbratt et al. 2010). Herschel provides exquisite sensitivity (able to reach the extragalactic confusion limit in most bands), over a 100–500 μm wavelength range that samples the vast majority of dust emission, with a 36'' limiting resolution at 500 μm that corresponds to a physical scale of <150 pc in M33 (the most distant of the Local Group's highly extended galaxies), and to <9 pc in the LMC.

However, FIR observations of Local Group targets come with their unique complications. Because the galaxies of the Local Group are so extended, they are vulnerable to the well-known problem of emission on large angular scales being filtered out of data produced by scanning detector arrays across the sky, as Herschel does (Meixner et al. 2013; Roussel 2013; Smith et al. 2017, 2021). The diffuse emission lost because of this effect will naturally tend to correspond to the diffuse, lower-density dust in a galaxy. This will systematically bias any analysis concerning ISM density.

All-sky FIR surveys, such as those by the Infrared Astronomical Satellite (IRAS; Neugebauer et al. 1984) and Planck (Planck Collaboration et al. 2011b), are free from the filtering of large-scale emission that affects Herschel. However, the resolution of IRAS and Planck is a factor of ∼10 worse than that of Herschel, severely limiting their ability to resolve the densest portions of the ISM and overcome the temperature mixing problem. Plus, these facilities provide sparser coverage of the dust SED than Herschel. These factors severely limit their use in studying the ISM in nearby galaxies.

Clearly, answering the open questions about D/H will require reliable Herschel observations of Local Group galaxies. Therefore, in the first paper of this series, Clark et al. (2021, hereafter Paper I), we produced new versions of the Herschel maps for the Local Group galaxies M31 (Andromeda), M33 (Triangulum), the LMC, and the SMC. These new maps were produced by combining Herschel observations, in Fourier space, with data from lower-resolution FIR telescopes that did not suffer from filtering, to restore the missing large-scale emission. These new maps also incorporate significant calibration corrections (up to 30% in some cases).

The new maps from Paper I cover two Herschel Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) bands, at 100 and 160 μm, and all three Herschel Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010) bands, at 250, 350, and 500 μm. As described in Section 6 of Paper I, we also produced Galactic-foreground-subtracted versions of the new maps; it is these that we use throughout this paper.

To illustrate the feathered FIR data, Herschel 250 μm maps for each galaxy are shown in the upper panels of Figure 2 for the LMC and SMC, and Figure 3 for M31, and M33.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.3. Molecular and Atomic-gas Data

2.4. Convolution to Limiting Resolution

To constrain the dust properties of our target galaxies, while also taking full advantage of the resolution provided by our data, we carried out pixel-by-pixel SED fitting of the FIR maps. For this, we applied a modified blackbody (MBB) model—specifically, a broken-emissivity modified blackbody (BEMBB) model (Gordon et al. 2014).

An MBB is essentially the simplest model that can be used to reliably fit FIR dust emission. It takes the following form, where the surface brightness, S(λ), of dust emission at a given wavelength, λ, can be expressed by:

$$\begin{eqnarray}&&S(\lambda )=\kappa (\lambda )\ B(\lambda ,{T}_{d})\ {{\rm{\Sigma }}}_{d}\end{eqnarray} \tag{ 1 }$$

where B is the Planck function evaluated at wavelength λ and dust temperature T_d , Σ_d is the dust mass surface density, and κ(λ) is the dust mass absorption coefficient at wavelength λ; κ(λ) varies with wavelength according to an emissivity law:

$$\begin{eqnarray}&&\kappa (\lambda )=\displaystyle \frac{\kappa ({\lambda }_{\mathrm{ref}})}{{\lambda }_{\mathrm{ref}}^{-\beta }}{\lambda }^{-\beta }\end{eqnarray} \tag{ 2 }$$

where κ(λ_ref) is the value of κ(λ) at a reference wavelength λ_ref, and β is the dust emissivity spectral index.

The MBB model is less directly physically motivated than dust grain models based on the radiative transfer and optical properties of different potential dust species (e.g., Draine et al. 2014; Jones et al. 2017; De Looze et al. 2019). However, these various dust grain models are not currently able to incorporate submillimeter excess, wherein more emission is observed in the Rayleigh–Jeans regime than would be otherwise expected from models, with the excess increasing toward longer wavelengths. The phenomena of submillimeter excess is observed primarily in dwarf galaxies, such as the Magellanic Clouds (Galliano et al. 2003; Bot et al. 2010; Rémy-Ruyer et al. 2013; Gordon et al. 2014), but is also seen in lower-density regions of more massive galaxies, including M33 (Paradis et al. 2012; Relaño et al. 2018).

It is therefore clearly preferable that we use an SED model that is able to incorporate submillimeter excess emission. ⁸  The BEMBB model achieves this by having the value of β change at some break wavelength λ_break; a shallower β at longer wavelengths captures the excess emission. As such, the BEMBB emissivity law takes the form:

$$\begin{eqnarray}&&\kappa (\lambda )=\displaystyle \frac{\kappa ({\lambda }_{\mathrm{ref}})}{{\lambda }_{\mathrm{ref}}^{-\beta }}E(\lambda )\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}E(\lambda )=\left\{\begin{array}{ll}{\lambda }^{-{\beta }_{1}} & \mathrm{for}\ \lambda \lt {\lambda }_{\mathrm{break}}\\ {\lambda }^{{\beta }_{2}-\beta 1}\,{\lambda }^{-{\beta }_{1}} & \mathrm{for}\ \lambda \geqslant {\lambda }_{\mathrm{break}}\end{array}\right.\end{eqnarray} \tag{ 4 }$$

for which β takes a value of β₁ at wavelengths <λ_break, while it takes a value of β₂ at wavelengths ≥λ_break.

We perform our SED fitting using the Dust Brute Force Fitter (DustBFF), a grid-based Bayesian SED-fitting code; DustBFF, and the mathematical formalism from which it operates, are presented in Gordon et al. (2014). We employ DustBFF in the same manner as in Paper I, so we refer the reader there for a full description of our implementation; the only differences in this work are that we used a slightly different setup of the parameter grid (to account for the wider range of densities our Herschel data can probe), and that the data we fitted were of course the Herschel 100–500 μm bands (and we therefore used covariance matrices specific to these bands, given below). We adopt a value of the dust mass absorption coefficient, at a reference wavelength of 160 μm, of κ₁₆₀ = 1.24 m² kg⁻¹, following Roman-Duval et al. (2017); this allows for ease of comparison with their previous work investigating D/H in the Local Group.

The DustBFF parameter grid we use for all of our sample galaxies is given in Table 2. Note that BEMBB implementation in DustBFF does not parameterize β₂ directly. Rather, the break in β is parameterized via the 500 μm excess, e₅₀₀, which gives the relative excess in the 500 μm flux, above what would arise from a standard MBB model (with negative values indicating a 500 μm deficit); e₅₀₀ is thus defined:

$$\begin{eqnarray}&&{e}_{500}={\left(\displaystyle \frac{{\lambda }_{\mathrm{break}}}{500\,\mu {\rm{m}}}\right)}^{{\beta }_{2}-{\beta }_{1}}-1\end{eqnarray} \tag{ 5 }$$

DustBFF uses a full covariance matrix in its model evaluations. This covariance matrix, ${ \mathcal C }$, is given by:

$$\begin{eqnarray}&&{ \mathcal C }={{ \mathcal C }}_{\mathrm{calib}}+{{ \mathcal C }}_{\mathrm{instr}}\end{eqnarray} \tag{ 6 }$$

where ${{ \mathcal C }}_{\mathrm{calib}}$ is a matrix that contains the calibration uncertainty for each band, and ${{ \mathcal C }}_{\mathrm{instr}}$ is a matrix incorporating the effect of the map's instrumental noise on the probability for each model.

The calibration covariance matrix ${{ \mathcal C }}_{\mathrm{calib}}$ is calculated by multiplying the proposed fluxes for a given model by the relative calibration uncertainty matrix ${{ \mathcal U }}_{\mathrm{calib}}$, which is itself given by summing ${{ \mathcal U }}_{\mathrm{uncorr}}$ and ${{ \mathcal U }}_{\mathrm{corr}}$, which are the matrices containing the fractional uncertainties that are uncorrelated and correlated between each band.

The diagonal values of the instrumental noise covariance matrix, ${{ \mathcal C }}_{\mathrm{instr}}$, were calculated for each galaxy in each band, by taking the median pixel value of the uncertainty map of each, added in quadrature to the uncertainty on the foreground subtraction (see Section 6 of Paper I). The uncertainty maps, as presented in Paper I, propagate the uncertainty arising from the feathering process used to restore large-scale emission to the Herschel data, along with other the instrumental noise contributions found in standard Herschel uncertainty maps. The off-diagonals elements of ${{ \mathcal C }}_{\mathrm{instr}}$ are all zero. The diagonal elements of ${{ \mathcal C }}_{\mathrm{instr}}$ for each galaxy and band are given in Table 3.

The uncorrelated uncertainty of a band essentially reflects the repeatability of its photometric measurements. The correlated uncertainties, on the other hand, arise from uncertainties in the photometric calibration. For instance, if all of the bands of an instrument were calibrated using a certain stellar model, but that model is only constrained to within a certain percentage, then all of that instruments' bands will share an uncertainty of that percentage; whatever the actual underlying error is, it will be the same between bands. Not accounting for correlated uncertainties can cause severe biases in model results (Galliano et al. 2011; Kelly et al. 2012; Veneziani et al. 2013).

For the Herschel-PACS 100 and 160 μm bands, the overall calibration uncertainty is taken as 7%; of this, we treat 2% as being correlated between bands, as this is the upper end of the quoted uncertainty on the continuum model of the five late-type giant stars used as the Herschel-PACS photometric calibrator sources (Decin & Eriksson 2007; Balog et al. 2014). We therefore assume the uncorrelated component is the quadrature subtraction of the correlated uncertainty from the total uncertainty, being $\sqrt{7{ \% }^{2}\mbox{--}2{ \% }^{2}}=6.7 \%$.

For the Herschel-SPIRE 250, 350, and 500 μm bands, we use an overall calibration uncertainty of 5.5%, of which 4% is taken to be correlated between bands (Griffin et al. 2010, 2013; Bendo et al. 2013)—i.e., the Herschel-SPIRE photometric uncertainty is in fact dominated by the correlated component, highlighting the importance of correctly accounting for this effect. This gives an uncorrelated component of $\sqrt{5.5{ \% }^{2}-4{ \% }^{2}}=3.8 \%$, again applying quadrature subtraction.

While using the full covariance matrices for fitting allows us to account for the often dramatic impacts of correlated uncertainties, a shortcoming is that it necessarily requires all uncertainties to be treated according to the same likelihood function—e.g., a Gaussian (see Section 4 of Gordon et al. 2014, specifically their Equation (16)). However, the uncertainties on absolute calibrations tend not to be Gaussian. Rather, the true value of the calibration error will be confidently contained within the stated bounds; this is the case for Herschel-PACS (Decin & Eriksson 2007) and Herschel-SPIRE (Bendo et al. 2013). In other words, unlike with a Gaussian tail, there is not a 32% chance of the true value being outside the stated ± range. Ordinarily, this would result in the total combination of the uncorrelated statistical (Gaussian) uncertainties and correlated absolute (non-Gaussian) uncertainties being smaller than if they were both Gaussian and added in quadrature. However, we also wish the total uncertainties, as stated in the diagonal elements of ${{ \mathcal U }}_{\mathrm{calib}}$, to reflect the canonical 7% and 5.5% uncertainties for PACS and SPIRE. To reconcile these two requirements, we use slightly larger uncorrelated uncertainties than the standard 2% for PACS (Balog et al. 2014) and 1.5% for SPIRE (Bendo et al. 2013), such that the quadrature sum of the uncorrelated and correlated components gives the canonical 7% and 5.5% values. These larger uncorrelated uncertainties also provide some leeway to account for additional sources of uncertainty (such as that caused by uncertainty on the beam area, which is at least 1%; Bendo et al. 2013).

We therefore use the correlated and uncorrelated relative uncertainty matrices

$$\begin{eqnarray}{{ \mathcal U }}_{\mathrm{corr}}=\left[\begin{array}{ccccc}0.02 & 0.02 & 0 & 0 & 0\\ 0.02 & 0.02 & 0 & 0 & 0\\ 0 & 0 & 0.04 & 0.04 & 0.04\\ 0 & 0 & 0.04 & 0.04 & 0.04\\ 0 & 0 & 0.04 & 0.04 & 0.04\end{array}\right]\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}{{ \mathcal U }}_{\mathrm{uncorr}}=\left[\begin{array}{ccccc}0.067 & 0 & 0 & 0 & 0\\ 0 & 0.067 & 0 & 0 & 0\\ 0 & 0 & 0.038 & 0 & 0\\ 0 & 0 & 0 & 0.038 & 0\\ 0 & 0 & 0 & 0 & 0.038\end{array}\right]\end{eqnarray} \tag{ 8 }$$

to give a final ${{ \mathcal U }}_{\mathrm{calib}}$ of

$$\begin{eqnarray}{{ \mathcal U }}_{\mathrm{calib}}=\left[\begin{array}{ccccc}0.07 & 0.02 & 0 & 0 & 0\\ 0.02 & 0.07 & 0 & 0 & 0\\ 0 & 0 & 0.055 & 0.04 & 0.04\\ 0 & 0 & 0.04 & 0.055 & 0.04\\ 0 & 0 & 0.04 & 0.04 & 0.055\end{array}\right]\end{eqnarray} \tag{ 9 }$$

for which the columns and rows contain the values for the Herschel bands in ascending order of wavelength.

In total, the parameter grid contains ≈9 million models, and across our galaxies, we fit it to ∼1 million pixels. Computing the parameter posterior probability distributions for all of these pixels for all of our galaxies took approximately 1 month using a 32 × 3.2 GHz thread computer. Propagating the full posterior distribution, with likelihoods for every model, for every pixel, throughout our analyses, would have been impractical. For each pixel, we therefore drew 1000 random samples (with replacement) from the full posterior distribution, with the probability of a given model being drawn being proportional to its likelihood. This set of random samples was then propagated as our posterior for all analyses. We show an example SED for a pixel in M33, with the 1000 posterior samples illustrated, in Figure 4, with a corner plot of that pixel's posterior in Figure 7.

Zoom In Zoom Out Reset image size

When deciding what pixels of our Herschel data to fit, we imposed a signal-to-noise ratio (S/N) threshold dictated by the Σ_H data. Because the motivation for our SED fitting here is to explore variation in the D/H ratio, there was no point fitting pixels for which there was not reliable gas data. In particular, we found that the very peripheries of the Σ_H maps, particularly for M31 and M33, seem to exhibit unphysically low gas surface density measurements, dropping abruptly to ∼0 at their edges. This seems to be caused by the sensitivity tailing off at the edge of the input maps' coverage areas. And in the Magellanic Cloud Σ_H data, there is conspicuous instrumental striping in very low-density areas; systematic bias from these artifacts seems likely to dominate over astrophysical emission (even when binning). We therefore imposed a sensitivity cut of S/N > 4 using our Σ_H maps. Because our Σ_H maps are a combination of the reprojected input maps of atomic and molecular gas, we calculated the noise ourselves, by measuring the rms noise within noise-dominated regions of the maps, where there was no apparent astrophysical signal. Only pixels above the S/N > 4 threshold underwent SED fitting and subsequent analysis. Conveniently, our Σ_H maps are very closely matched in sensitivity, with the S/N > 4 threshold corresponding to Σ_H > 2.8 M_⊙ pc⁻² for the LMC, Σ_H >2.6 M_⊙ pc⁻² for M31, and Σ_H > 2.1 M_⊙ pc⁻² for the SMC and M33 (note that these thresholds are before accounting for inclination projection effects, as discussed in Section 4.1).

3.1. SED Fitting Results

The results of our DustBFF SED fitting are shown in Figures 5 and 6, which shows maps of the median (i.e., 50th percentile) value of each parameter for each pixel, in each galaxy. Note that maps of the parameter medians for each pixel tend to be noisier than maps of each pixels's maximum a posteriori (MAP) SED parameter values; however, they also tend to be more representative of the marginalized distribution for each individual parameter. Hence, we chose to display these. While the SED fitting is not the focus of this work in-and-of itself, these outputs provide a good check of DustBFF, and show some interesting features, so we do make a few comments here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our SED fitting appears to do a good job of modeling the data. For every pixel, we computed the residuals between the MAP SED model and the data, and measured the median residual across all pixels for each band, for each galaxy, and found it to be <2.4% in every case. No particular subregion or environment appears to exhibit noticeable residuals, either.

Overall, our maps of SED parameters are in decent agreement with those of previous authors who have done similar resolved SED fitting of these galaxies. For the rest of this subsection, we perform comparisons to several specific works. ⁹

Gordon et al. (2014) previously performed resolved SED fitting of the Magellanic Clouds as part of the original Herschel Inventory of The Agents of Galaxy Evolution (HERITAGE; Meixner et al. 2013) key program. Although the general ranges of values and broad morphology of our parameter maps agree well, it is difficult to make more detailed comparisons; the striping/cross-hatching artifacts found in the original HERITAGE data are also strongly apparent in their parameter maps; plus, the S/N limits of the HERITAGE data, especially for the SMC, renders the modeled areas much smaller than for our data, with its restored diffuse emission and reduced noise.

We do note that like Gordon et al. (2014), we find very low values of β₁ and β₂ in the SMC, with many pixels in the galaxy's outskirts having values only slightly above 1. We do not believe that this is due to significant temperature mixing, causing a flat SED shape to be favored. In order for temperature mixing to cause β₂ to be significantly flatter than β₁, the colder dust components must have a Wein's Law peak at a wavelength > λ_break. However, > 94% of SMC pixels have λ_break > 290 μm. For dust to have a Wein's Law peak beyond this wavelength requires a temperature of <10 K. Temperatures this cold are hard to explain physically, and Bot et al. (2010) reported that in order to explain SED flattening in the SMC, the temperatures would have to be so cold as to approach the cosmic microwave background. Additionally, Gordon et al. (2014) showed that the required mass of dust this cold would exceed the available mass of ISM metals in the SMC. Similarly, if temperature mixing of dust at temperatures above 10 K was artificially flattening the SED, we would expect significant residuals between the models and observations. This is because even an SED with a shallow β will fail to closely fit the flat-peaked SED arising from significant temperature mixing, resulting in residuals (especially when averaged over large numbers of pixels, such as in our data). However, as discussed earlier in this section, there are no meaningful residuals evident in our SED fits.

The work of Utomo et al. (2019) provides for a particularly useful comparison to our own, as they carried out pixel-by-pixel SED fitting of Herschel data for all of the galaxies in our sample, using an algorithm substantially similar to DustBFF—albeit without the restoration of extended emission available in our data, and using an SED model that assumed a fixed β = 1.8 with no break at longer wavelengths. Their Σ_d and T_d parameter maps (see their Figure 1) are a good match to our own, for the areas where we both have results. The Utomo et al. (2019) maps of Σ_d , in particular, are almost a perfect match in morphology to ours. Our maps also agree with their general structure for T_d , too, although we differ on some of the specifics. For instance, we find the same areas of very high T_d in M33 around star-forming regions, and we also find the "ridge" of highest temperature in the SMC to be slightly offset northwest of the "ridge" of maximal Σ_d . On the other hand, the morphology of this SMC temperature structure differs somewhat between our respective maps, and we also find M31 to have consistently lower T_d than the other galaxies in the sample (dust-mass-weighted median of 15.9 K in M31, compared to 21.3, 22.0, and 18.4 K for the LMC, SMC, and M33, respectively; see Table 4), while Utomo et al. (2019) do not report such a difference. ¹⁰  The fact that we allow β to vary, whereas Utomo et al. (2019) kept it fixed, could account for this difference in temperature, given the strong degeneracy between those two parameters.

Table 4. Dust and Gas Parameters for Each of Our Galaxies

	M31	M33	LMC	SMC
Md (106 M⊙)	45.0	5.41	1.35	0.165
Td (K)	15.9	18.4	21.3	22.0
β1	2.12	1.81	1.73	1.49
β2	1.96	1.57	1.54	1.09
λbreak (μm)	277	283	272	283
e500	0.21	0.28	0.28	0.41
D/H	10−2.11	10−2.13	10−2.39	10−3.38

Note. The total dust mass, M_d , is the sum of the median posterior mass value of every pixel; the D/H is this value divided by the sum of the same pixels in the Σ_H map. For the other parameters, we give the dust-mass-weighted average of the pixel medians.

Download table as:  ASCIITypeset image

The colder dust temperatures we find for M31 do, however, agree well with Smith et al. (2012), Viaene et al. (2014), and Whitworth et al. (2019), ¹¹  who found ≈16 K dust temperatures over most of the disk, with a much smaller warmer region in the center, small enough that it would mostly fall within the low-Σ_H region we do not fit. Our results also agree with the finding from Smith et al. (2012) and Whitworth et al. (2019) that β decreases with radius in M31. Indeed, whereas both of those studies used an unbroken-β model, we find that both β₁ and β₂ fall with radius. Additionally, thanks to our restoration of the diffuse emission, we are able to perform our SED modeling out to a radius of over 32 kpc, compared to the 18 kpc of Smith et al. (2012), and 20 kpc of Whitworth et al. (2019), finding that the trends of decreasing T_d and β continue out to larger extreme radii. We are also able to successfully perform SED fitting in low-S/N pixels within diffuse regions between the star-forming rings of M31, which those previous studies did not model.

For M33, although we find that T_d and β₁ generally fall with radius, in line with the findings of Tabatabaei et al. (2014)—who also used a single-MBB, unbroken-β model—we find that β₂ does not fall conspicuously with radius, instead being depressed around regions of heightened star formation (the same regions with elevated dust temperature), suggestive of possible grain processing in these environments.

3.1.1. Global Dust and Gas Values

In Table 4, we provide global SED parameters for our galaxies. We report the total dust mass, M_d , where for each galaxy we take the sum of the posterior median dust mass for each pixel. For the D/H values (discussed more fully in Section 4), we divide this by the total hydrogen mass computed from the same pixels in our Σ_H maps. For the other SED parameters, we give the dust-mass-weighted average of the posterior medians for each pixel. In all cases, we only consider pixels that exceed our previously discussed Σ_H S/N criterion. ¹²

To make sure that our imposition of the S/N criterion does not significantly affect the global D/H values, we also recalculated them using all pixels in the dust and gas maps for which the Σ_H map had S/N > 1; this much weaker threshold is intended to make sure we are not excluding significant amounts of diffuse emission, while also still avoiding artifacts and other image errors in low-S/N regions introducing significant bias into the measurements. To prevent contamination from any Galactic foreground emission that had not been fully subtracted, we also only counted pixels for M31 and M33 that are within the elliptical apertures employed for photometry in Section 6.3 of Paper I. The resulting D/H were all with in 3% of those we found when using the S/N > 4 criterion, indicating a very minimal contribution from the masked pixels.

Gordon et al. (2014) also performed their SED fitting with DustBFF, using a BEMBB model. This allows for a good direct comparison of our global results to results obtained using the older, unfeathered Herschel maps of the Magellanic Clouds. The dust masses we find for the LMC and SMC are factors of 1.59 and 1.48 smaller, respectively, than the total dust masses found by Gordon et al. (2014), after correcting for our differing values of κ(λ_ref). Similarly, in comparison to the dust masses from the resolved SED fitting of Utomo et al. (2019; in which they model the illuminating interstellar radiation field, as per Draine & Li 2007), ours are smaller by factors of 3.3 for the LMC, 3.9 for the SMC, 1.46 for M31, and 3.8 for M33 (again, after scaling to our κ(λ_ref)).

The fact we tend to find find lower masses may seem surprising at first, given that the new Herschel maps from Paper I restored dust emission that was missing in older maps. However, there was more emission that needed restoring in the shorter wavelength bands (especially in the Herschel-PACS 100 and 160 μm data) than at longer wavelengths. As a result, essentially all of the dust emission from our galaxies is rendered significantly bluer in the new maps (see Figure 17 of Paper I). Bluer emission corresponds to warmer dust temperatures; and as the luminosity of a given mass of dust goes approximately $\propto {T}_{d}^{4+\beta }$, this reduces the modeled dust mass for a given FIR brightness.

To check this, we directly compared the dust-mass-weighed dust temperatures we find, to those of Gordon et al. (2014), by reprojecting their parameter maps to our pixel grid, and only comparing pixels for which both had coverage. We found that our mass-weighted-average global dust temperatures were 2.0 K warmer for the LMC, and 5.5 K warmer for the SMC. These correspond to factor of 1.11 and 1.29 differences in temperature, respectively. Assuming there were no other difference, then the simple $\propto {T}_{d}^{4+\beta }$ relation would suggest that we should expect our dust masses to be factors of 2.1 and 5.5 reduced due to this effect, assuming typical β = 1.7. We do not perform the same direct comparison with the Utomo et al. (2019) maps, due to differences between the foreground subtractions & effective apertures we use, and differences in method. For instance, by using a broken-emissivity approach (as compared to their fixed β = 1.8 method), it is possible for flux at longer wavelengths to be accounted for by the β becoming flatter, instead of by driving up the total SED normalization—and hence mass. The new Herschel maps from Paper I generally increase the total FIR flux measured in each pixel, with per-band average increases of 21%, which will counteract some of the mass reduction due to higher temperature. Plus, Paper I applied calibration corrections to the Herschel data, which will also have knock-on effects to the SED parameters (generally increasing 100–160 μm fluxes, but increasing 250–500 μm fluxes. So overall, the reduction in dust mass we find is of the order of the difference that should be expected.

On the other hand, comparing to the total dust masses reported by Chastenet et al. (2017), who used the THEMIS physical dust grain model (Jones et al. 2017), our LMC dust mass is 20% greater than theirs, although our SMC dust mass is a factor of 1.67 less, after κ(λ_ref) corrections. ¹³  Chastenet et al. (2017) noted that their masses are significantly lower than other authors have found, which they ascribe to the high fraction of more-emissive carbonaceous grains in their modeling results (which may be reasonable, given recent results from Roman-Duval et al. 2022b, who found LMC and SMC dust to be more carbon-rich than that of the Milky Way, at a given column density).

3.1.2. Submillimeter Excess

Significant submillimeter excess is apparent in the outskirts of all of our galaxies, especially the lowest-metallicity SMC. This is in line with previous work showing that greater submillimeter excess tends to be found in lower-metallicity environments (Galliano et al. 2003; Bot et al. 2010; Rémy-Ruyer et al. 2013; Gordon et al. 2014). However, several works have also found good evidence that submillimeter excess can be driven by ISM density, with more excess found in lower-density environments (Planck Collaboration et al. 2011a; Relaño et al. 2018). Our improved data put us in a good position to compare these possibilities.

Thanks to our restoration of diffuse emission, we are able to conduct SED fitting out to larger radii than previously performed. For instance, we are able to model dust emission in the H i filament to the southeast of the LMC, where it is interacting with the lower-metallicity Magellanic Stream (Nidever et al. 2008), in which we find elevated levels of submillimeter excess. Similarly, in M31, we are able to perform our SED modeling out to a radius of over 32 kpc, compared to 18 kpc in Smith et al. (2012), and 20 kpc in Whitworth et al. (2019) and Draine et al. (2014). At these previously unprobed extreme radii, we again find that submillimeter excess increases to higher levels.

Note that with our BEMBB model, submillimeter excess corresponds to β₂ < β₁. While our model grid does allow for β₂ ≥ β₁ (i.e., a submillimeter deficit), this only happens for <3.7% of pixels in the median e₅₀₀ maps across all four of our galaxies (being as few as 0.03% of pixels for the SMC). In other words, the dust SED effectively always becomes shallower toward longer wavelengths (agreeing with Planck Collaboration et al. 2014). We find that the break wavelength is in the range 260 μm < λ_break < 325 μm for 80% of pixels, with λ_break generally being lower in areas of greater ISM density, albeit with elevated λ_break often seen immediately around star-forming regions (see Figure 5).

Spearman rank correlation tests find that the strength of the submillimeter excess in M31 and M33 is more strongly correlated with deprojected radius than with Σ_H; for the LMC, we find that the correlation is stronger with density than with radius, although the difference is smaller. However we note that the LMC is much more disturbed than either of the spirals. Correlation coefficients are given in Table 5. We do not test correlation with radius for the SMC, due to its high degree of disturbance preventing us from meaningfully quantifying radii.

Table 5. Spearman Rank Correlation Coefficients between the Submillimeter Excess at 500 μm, e₅₀₀, and Various Other Parameters, for the Pixels in Our Maps

e500 versus	M31	M33	LMC	SMC
ΣH	−0.54	−0.16	−0.47	−0.39
Deprojected Radius	0.41	0.74	0.57	⋯
D/H	−0.80	−0.71	−0.57	−0.60
Σd	−0.77	−0.65	−0.67	−0.68

Note. All relations have ${{ \mathcal P }}_{\mathrm{null}}\lt {10}^{-16}$.

Download table as:  ASCIITypeset image

The fact that the correlation with radius is stronger than with ISM density for M31 and M33, but not the LMC (which lacks a significant metallicity gradient; see Section 4.1.1), potentially suggest that metallicity is playing a role, either through different dust composition, or through improved shielding. These possibilities are supported by the fact that we find even stronger correlation of e₅₀₀ with D/H and Σ_d (Table 5). We caution against over-interpretation of these final two correlations, however, as both depend upon Σ_d , which is a parameter in our SED fitting model along with e₅₀₀, and the two can be degenerate (e.g., Figure 7). We also see in M31 and M33 that there are narrow regions of elevated e₅₀₀ tracing the spiral arms (see Figure 6), despite these being higher-density regions. If shielding in dust-rich regions can drive down submillimeter excess, then it is conceivable that proximity to radiation and shocks from ongoing star formation elevates e₅₀₀ in certain parts of the spiral arms, possibly in lower-density regions carved by young stellar winds, on smaller scales than what our data can resolve in these galaxies.

Zoom In Zoom Out Reset image size

We note that the S/N of the FIR emission at larger radii tends to be very low, leading to large uncertainties on the modeled parameters for any individual pixel. However, the fact that the trend of increasing submillimeter excess is evident at all azimuths, across thousands of independent pixels, gives us reason to be confident in the trend, and that it is not an artifact of our foreground subtraction, or arising from localized cosmic microwave background fluctuations, etc.

With maps of both the dust and hydrogen content of our galaxies, we can now explore the behavior of the D/G within them. In Figure 8, we show maps of D/H for our targets.

Zoom In Zoom Out Reset image size

For the two spirals, we find that D/H is elevated in the central regions, and tends to fall with increasing radius, in line with the broad trends observed by previous authors (Draine et al. 2014; Gratier et al. 2017). This trend is particularly pronounced in M31, while in M33 elevated D/H more closely follows the spiral structure, rather than being purely radial.

For the LMC, the structure of the D/H map broadly traces the morphology of the dust mass (see Figure 8), but with some other noteworthy features. The highest D/H is often found at the very edges of supershells excavated by recent star formation. ¹⁴  These supershells are apparent especially in poor H i (Dawson et al. 2013); the low gas density within the shells means that they mostly fall below our S/N threshold in Σ_H. The high D/H at the edges of these shells might indicate that the gas is being blown away more efficiently than the dust (Draine 2011). Alternatively, our dust and gas mass tracers might behave differently here; for instance, environmental grain processing might be affecting the dust opacity (see Appendix B.6).

Another noteworthy feature in the map for the LMC is the significantly depressed D/H to the southeast. This is coincident with the portion of the disk where low-metallicity material from the Magellanic Stream is being accreted onto the LMC (Bekki & Chiba 2007; Nidever et al. 2008). This accretion is leading to the enhanced star formation of 30 Doradus and the surrounding area (Fukui et al. 2017; Tsuge et al. 2019, 2020). Although the area of star formation itself has fairly high D/H > 10^−2.3, the southeastern fringes of the LMC have the lowest D/H found in the entire galaxy. Indeed, the southernmost tip of this region features D/H as low as 10^−3.3, matching the typical D/H we find for the SMC.

The D/H map of the SMC shows it to have a conspicuously lower D/H than any of the other galaxies in the sample, as we would expect from its 0.2 Z_⊙ metallicity. Thanks to our new Herschel maps from Paper I, and the sensitive Σ_H data, we are able to trace D/H out to very large radii. Away from the main regions of star formation along the Bar and Wing, D/H falls steadily, to values of less than 10⁻⁴ in the most diffuse regions.

4.1. Evolution of D/H with Σ_H

In Figure 9, we present the central result of this paper—the evolution of D/H with Σ_H for our galaxies. As covered in Section 1, we should expect D/H to increase with greater Σ_H, due to the increased efficiency of the accretion of gas-phase metals onto dust grains in higher-density ISM. To construct this plot, we binned the values from our D/H maps according to the Σ_H of each pixel, into bins of width 0.025 dex. The plotted points in Figure 9 show the median value in each bin. The error bars incorporate the uncertainty on the median, calculated by performing 500 Monte Carlo bootstrap resamples of all of the D/H values in each bin, recomputing the median each time, then finding the standard deviation of those 500 bootstrapped medians; this was then added in quadrature to the 0.05 dex uncertainty arising from the fact that our Σ_d values were taken from a model grid with a 0.05 dex step size (see Table 2), which limits the precision of any given D/H measurement. The scatter within each bin can be quite large; across all four galaxies, the average standard deviation of the D/H values contained within a bin is 0.8 dex. Despite this, the distribution of values within each bin tends to be well behaved and Gaussian; hence, the average uncertainty on the bin medians, as measured via the bootstrap resampling, is only 0.02 dex (i.e., the uncertainty on the median for each bin is almost always dominated by the 0.05 dex contribution of the uncertainty on the grid step size).

Zoom In Zoom Out Reset image size

The x-axis surface density values in Figure 9 have had a deprojection correction applied, to account for the effect of each galaxy's inclination. ¹⁵  For an inclined disk galaxy like M31, the column of ISM sampled by a given pixel passes through a greater thickness of disk than it would for an equivalent face-on galaxy. This compromises our ability to treat observed surface density as a proxy for average volume density (as volume density is what will be dictating grain growth, etc.), and would prevent us from fairly comparing our galaxies' D/H at a given observed Σ_H. Therefore, for a given galaxy inclination i (see Table 1), we apply a deprojection correction:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}=\cos (i)\ {{\rm{\Sigma }}}_{{\rm{H}}}.\end{eqnarray} \tag{ 10 }$$

Hereafter, we refer to ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ in instances where comparison between galaxies makes the distinction between corrected and uncorrected Σ_H values important. The $\cos (i)$ correction factor equates to 0.22, 0.56, and 0.90 for M31, M33, and the LMC, respectively. No correction needs to be applied to the y-axis quantity of D/H, as this is a ratio of two quantities that experience the same projection, canceling out its effect.

Because the SMC's highly disturbed morphology lacks a disk, inclination simply is not an applicable concept, so we apply no correction. In addition, thanks to the SMC's extreme elongation along the line of sight (Scowcroft et al. 2016), the link between observed surface density and actual volume density is further weakened in the case of the SMC; the reader should bear in mind that the entire trend for the SMC in Figure 9 could be shifted left or right by this poorly constrained systematic.

Figure 9 shows clearly that there is very strong evolution in the D/G with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ for all of our sample galaxies. Indeed, M31 shows almost an order of magnitude of evolution in D/H between densities of 0.6–4 M_⊙ pc², while the SMC exhibits even more than that, with its D/H increasing by a factor of 20 between 3 and 150 M_⊙ pc². For all four galaxies, D/H increases steadily with density—although only up until a point for M31, M33, and the LMC. All three of these all exhibit a turnover in D/H, after which it begins to decrease with increasing Σ_H (this is explored in Section 5). Regardless, the increase in D/H with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ over most densities for all four galaxies provides evidence for significant density-dependent dust grain growth.

These results update our understanding of just how much D/H can vary within a galaxy, especially at fixed metallicity. Previous authors have found the D/G varying by a factor of 3 in the LMC (Roman-Duval et al. 2017), a factor of 5 in M31 (Draine et al. 2014), a factor of 7 in the SMC (Roman-Duval et al. 2014, 2017), and a factor of 10 in M101 (Chiang et al. 2018; Chastenet et al. 2021; albeit over a factor of 5 variation in metallicity, whereas our variation is seen even with little-or-no metallicity variation, as per Section 4.1.1). We now find that D/H spans a factor of 2.5 in M33, a factor of 3.5 in the LMC, a factor of 9.0 in M31, and a factor of 22.4 in the SMC.

4.1.1. Effect of Metallicity Gradients on Evolution in D/H

If there were significant systematic variations in the ISM metallicity of our target galaxies, then an increase in D/H might not indicate grain growth—instead, if regions of greater density also had higher metallicity, then D/H would be correspondingly greater even if the fraction of metals in dust grains was staying constant (in other words, even if there was no grain growth). However, this is not the case for our galaxies.

Neither the ISM (Pagel et al. 1978; Toribio San Cipriano et al. 2017) nor younger stellar clusters (Grocholski et al. 2006; Cioni 2009) in the LMC exhibit a significant radial metallicity gradient (<0.2 dex, less than the intrinsic scatter). Although metallicity does appear depressed in its southeast periphery, where the LMC is interacting with low-metallicity gas of the Magellanic Stream (Nidever et al. 2008; Tsuge et al. 2019, 2020; Roman-Duval et al. 2021), this is true for both high- and low-density gas in this region, so should not systematically influence the D/H evolution profile of the LMC. For the SMC, there is no significant metallicity gradient observed in the ISM (Pagel et al. 1978; Cioni 2009; Toribio San Cipriano et al. 2017).

There is a definite radial metallicity gradient in the ISM of M31; however, it is only −0.56 dex ${R}_{25}^{-1}$ (Zurita & Bresolin 2012); ¹⁶  so clearly this cannot account for the 0.9 dex evolution in D/H in M31. Moreover, metallicity does not scale monotonically with Σ_H in M31, thanks to the low ISM density of the inner ∼3 kpc—further reducing the potential for the radial metallicity gradient to give rise to increased D/H in regions of greater Σ_H.

The ISM of M33 has an even smaller radial metallicity gradient, at −0.22 dex ${R}_{25}^{-1}$ (Magrini et al. 2016); similar to M31, this is unable to account for the 0.4 dex evolution in D/H we see within this galaxy, especially because the gradient is comparable to the intrinsic scatter in metallicity at any given radius (Magrini et al. 2016).

We therefore do not believe that metallicity effects are a primary driver of the evolution in D/H observed in Figure 9.

4.1.2. Effect of Dust Destruction at Lower Densities on Evolution in D/H

While there are compelling reasons to expect D/H to increase at higher densities due to more efficient interstellar grain growth, we should also expect dust grain destruction to become more efficient at lower densities. Therefore, a trend of D/H increasing with density could potentially be driven by either effect.

The primary forces of dust destruction in the ISM are expected to be sputtering due to supernovæ shocks, and photodestruction by high-energy photons (Jones 2004; Bocchio et al. 2014; Slavin et al. 2015).

Generally, theoretical models expect the efficiency of grain destruction due to supernovæ shocks to only be relatively weakly dependent upon density. Typically, destruction efficiency is modeled to fall with increasing density, according to a power-law slope of around −0.1 (Jones et al. 1994; Dwek et al. 2007; Temim et al. 2015). In contrast, all four of our galaxies show evolution in D/H with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ with power-law slopes ranging from 0.4 for the LMC, up to 1.0 for M31 (up until the D/H turnover; see Section 5). This implies that either dust destruction due to supernovæ is not driving the vast majority of D/H evolution with density we find, or that the density-dependence of supernovæ dust-destruction efficiency has been underestimated by theoretical models.

The rate of interstellar dust destruction should also be expected to closely correlate with the density of young stars in a given environment. This is because young stars will produce high-energy photons capable of photodestruction of grains, and because young stars indicate where dust-destroying core-collapse supernovæ should be occurring at the greatest rate.

UV emission traces young stars that have formed within the past ∼100 Myr (Kennicutt 1998; Calzetti et al. 2005; Buat et al. 2011). Therefore UV observations should allow us to trace the relative rate of dust destruction we should expect in different locations. The only high-quality high-resolution wide-area UV data available for all four of our sample galaxies comes from the UltraViolet/Optical Telescope (UVOT; Roming et al. 2000, 2004) on the Neil Gehrels Swift Observatory (Gehrels et al. 2004). For the Magellanic Clouds, the Swift-UVOT data we use is the data presented in Hagen et al. (2017); for M31 and M33, we use reductions created according to the same method as in Hagen et al. (2017), with that data to be fully be presented in M. Decleir et al. (2023, in preparation; for data access see Section 7). The data for M31 and M33 cover the full stellar disks of both galaxies; the SMC data covers the whole bar and part of the wing, while the LMC data covers the bar, 30 Doradus, and the surrounding regions.

To trace the young stars, we use the Swift-UVOT W2 band maps. As it is the UVOT band with the shortest effective wavelength, at 192 nm, W2 should best trace the higher-energy photons that will lead to dust destruction. We converted the Swift-UVOT W2 count-rate data to SI units, as per the AB magnitude zero-points given in Breeveld et al. (2011), and thence to UV luminosity surface density, in L_⊙ pc⁻². We then reprojected these maps to the same pixel grid as our D/H maps, and applied deprojection corrections.

In Figure 10, we plot D/H against UV luminosity density for each of our galaxies. To our surprise, we found that regions of great UV luminosity density tend to have larger D/H. This trend is significant for all four galaxies; in every case, Spearman rank correlation tests found ${{ \mathcal P }}_{\mathrm{null}}\lt {10}^{-10}$, with correlation coefficients ranging from 0.37–0.50. On large scales, we might expect this sort of correlation, as the denser ISM where grain growth can occur is also the material that should be fueling ongoing star formation. But in data with resolution as high as ours, as good as 14 pc in the LMC, it is surprising to find that regions with recent star formation, and a corresponding higher density of supernovæ shocks and hard radiation, still manage to have higher D/H.

Zoom In Zoom Out Reset image size

The origin of this surprising finding warrants further investigation in future work. One conceivable explanation is that we are tracing the creation of fresh dust by core-collapse supernovæ, which will be preferentially found in regions of recent star formation. If so, this has major implications for the role of supernovæ as net creators, versus destroyers, of dust (Nozawa et al. 2006; Priestley et al. 2021). Dust creation by ABG stars, seems unlikely to be the cause. This is because D/H is visibly well correlated with the ISM structure in our galaxies; AGB stars, however, will have a distribution following that of each galaxy's evolved stellar population, which has a very different distribution in these galaxies.

We also note that Figure B4 in Appendix B.4 likewise suggests that greater ionized gas surface density, as traced by Hα observations of the LMC and SMC, does not appear to be associated with reduced D/H, and any ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$. However, as this data is not available for M31 and M33, and with poorer resolution than provided by Swift-UVOT, we consider this a secondary line of evidence.

For all four galaxies, these results suggest that neither dust destruction by young stars (via hard radiation or shocks from core-collapse supernovæ), nor greater efficiency of dust destruction in lower-density ISM, nor the presence of ionized gas, are the primary driver of the evolution in D/H with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ observed in Figure 9. This supports the explanation that this relationship is tracing increasing dust grain growth at higher densities.

4.1.3. Variation in D/H Evolution between Galaxies

There are some striking similarities and differences between the relationships in how D/H evolves with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ for each of the galaxies in Figure 9. The first thing of note is how extremely similar the profiles of M31 and M33 are to each other. Over the 0.6–30 M_⊙ pc⁻² range in ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, they overlap nearly perfectly. ¹⁷  They even experience their D/H turnovers at the same surface density, ≈4 M_⊙ pc⁻². This would superficially suggest that the ISM in these galaxies has similar properties. This is, however, rather surprising, given that M31 has an average metallicity more than 2.5 times greater than M33, and a stellar mass more than 25 times greater.

The representative dust SED parameters for M31 and M33 do indeed differ (see Table 4), most notably in terms of β. M31 has higher β values (i.e., a steeper Rayleigh–Jeans slope), with minimal difference between β₁ and β₂. M33, on the other hand, has much lower values of β, more characteristic of lower-mass and dwarf galaxies (Rémy-Ruyer et al. 2013). Plus, M33 has a much more pronounced break to a shallower β₂ at longer wavelength; again more characteristic of dwarf galaxies (Gordon et al. 2014). This is extremely interesting, as β can potentially trace actual physical properties of the dust grains in question. In general, higher values of β are expected from crystalline grains, silicate species, or coagulated grains; while lower values of β are expected from metallic grains, amorphous grains, and carbonaceous species (Tielens & Allamandola 1987; Köhler et al. 2015; Jones et al. 2017; Ysard et al. 2018). The well-known temperature–β degeneracy (Shetty et al. 2009; Smith et al. 2013) limits what we can infer from the parameters derived from any one pixel. However, because we sample very large numbers of pixels, and because the pixels are larger than the instrumental PSF (and hence statistically independent), the overall trends we find ought to be robust against the temperature–β degeneracy (Smith et al. 2012). This is because the degeneracy does not impart a bias to the results of the fits; any trend observed over very large numbers of pixels cannot therefore be caused by the degeneracy. Large-scale variation in β therefore hints at different dust compositions (either because of their elemental compositions, or the environments in which the grains formed; see Roman-Duval et al. (2022b). But despite this, the global D/H ratios for M31 and M33 are effectively identical, being 10^−2.11 versus 10^−2.13—matching their similarity in the D/H versus ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ profile.

In contrast, if any two galaxies in our sample might be expected to have ISM that behaves similarly, it would be M33 and the LMC. Their stellar masses, metallicities, and star formation rates (Harris & Zaritsky 2009; Verley et al. 2009) are very similar. Both galaxies are bound to massive spiral companions, with which they have undergone interactions (although M33 is more distant from its companion, and considerably less disturbed, with the last major encounter likely > 1 Gyr ago; Bekki 2008; McConnachie et al. 2009). However, despite their similarity, the two galaxies show strikingly different profiles in their D/H evolution in Figure 9. The profile of the LMC would appear depressed compared to that of M33 by a factor of 1.5–3. And while both galaxies have a turnover in D/H, the location of that turnover is at a deprojected surface density of 40 M_⊙ pc⁻² for the LMC, compared to just 4 M_⊙ pc⁻² for M33 (and M31). Moreover, M33 (and M31) shows very steep increase in D/H with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ before the turnover, with a power-law index of 1.1 over the 0.6–4.0 M_⊙ pc⁻² range; whereas the D/H increase for the LMC happens much more gradually, with a power-law index of only 0.39 over the 2.5–40 M_⊙ pc⁻² range.

The fact that two galaxies as superficially similar as M33 and the LMC can have such starkly different D/H properties is significant, especially in the context of the common practice of inferring a galaxy's gas mass from observations of dust emission (Eales et al. 2010; Scoville et al. 2014). This is a standard method for estimating the ISM mass of galaxies for which direct measurements of CO and/or H i are not available. This has been shown to be remarkably reliable for massive galaxies at high masses (Scoville et al. 2014, 2016), but clearly extending the technique to more intermediate masses will have to be done with care, given that even fundamental galaxy properties like metallicity and stellar mass are not accurate predictors of a galaxy's D/H ratio. Moreover, while D/H for M33 and the LMC only differ by a factor of 2, it should be remembered that the method of estimating gas mass from dust emission tends not to use dust mass, but instead simply use dust luminosity in some longer-wavelength band such as 500 or 850 μm—and the 500 μm luminosities of M33 and the LMC differ by a factor of >3.

Lastly, the SMC clearly follows a very different evolutionary profile than the other galaxies. Not only does D/H continue to increase over the entire factor of 50 in surface density we sample, but moreover the gradient of the profile continues to get steeper, even up to the very highest density we are able to trace. Indeed, above 150 M_⊙ pc⁻², the highest ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ available for the SMC, the D/H profile for the SMC appears as though it may intersect that of the LMC. This suggests that grain growth may have a particularly strong dependence on density in the SMC. We can also be sure that the dust grains themselves have different properties in the SMC. Not only do most regions in the SMC appear to lack the 2175 Å extinction bump seen in higher-metallicity systems (Gordon et al. 2003; Murray et al. 2019), but our SED fitting finds the SMC to have a much lower β than the other galaxies in our sample, along with the most significant flattening in β at longer wavelengths (see Table 4 and Figure 6), all indicating a different composition.

4.2. Modeling the Evolution in D/H

In order to aid our understanding of the D/H evolution profiles in Figure 9, we use the dust evolution model of Asano et al. (2013). This model traces how various galactic environmental parameters can affect accretion of metals onto dust grains—such as the density of the ISM, the metallicity of the gas, the temperature of the existing dust grains, and the characteristic lifetime of the molecular clouds where grain growth occurs, and the average grain size. The model balances this with the rate at which supernovæ can destroy dust grains, given the density and metallicity of the ISM.

Using this model, the timescale for accretion of metals onto dust grains, τ_acc, is given by:

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\mathrm{acc}} & = & 2\times {10}^{7}\left(\displaystyle \frac{\overline{a}}{0.1\,\mu {\rm{m}}}\right){\left(\displaystyle \frac{{n}_{{\rm{H}}}}{100\,{\mathrm{cm}}^{-3}}\right)}^{-1}\\ & & \times \,{\left(\displaystyle \frac{{T}_{d}}{50\,{\rm{K}}}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{Z}{0.02}\right)}^{-1}\end{array}\end{eqnarray} \tag{ 11 }$$

where $\overline{a}$ is the average dust grain size in microns, n_H is the volumetric number density of hydrogen in cm⁻³, and Z is the metallicity (in units of the metal mass fraction). With τ_acc, one can then compute the fraction of the metal mass that is in dust grains (see Zhukovska et al. 2008), using:

$$\begin{eqnarray}&&{f}_{d}=\displaystyle \frac{{e}^{{t}_{\mathrm{growth}}/{\tau }_{\mathrm{acc}}}{f}_{{d}_{0}}}{1-{f}_{{d}_{0}}+({f}_{{d}_{0}}{e}^{{t}_{\mathrm{growth}}/{\tau }_{\mathrm{acc}}})}\end{eqnarray} \tag{ 12 }$$

where t_growth is the average duration of episodes of grain growth (if grain growth happens predominantly in molecular clouds, then t_growth corresponds to the average molecular cloud lifetime), and ${f}_{{d}_{0}}$ is the minimum fraction of metals that can be locked up in dust grains. Given this, the D/H for a given set of parameters is simply given by:

$$\begin{eqnarray}&&\displaystyle \frac{D}{H}={f}_{d}\,Z\end{eqnarray} \tag{ 13 }$$

which naturally gives the result that D/H will peak at a value equal to the metallicity, when all metals are locked up in dust grains at f_d = 1. This grain growth trend in the Asano et al. (2013) model is qualitatively similar to that found in the Nanni et al. (2020) model, also.

Because we are observing surface densities, we are unable to directly measure n_H. We therefore need to incorporate a conversion factor H_Σ⇒n , in units of ${\mathrm{cm}}^{-3}\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$, such that:

$$\begin{eqnarray}&&{n}_{{\rm{H}}}={H}_{{\rm{\Sigma }}\Rightarrow n}\times {{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}.\end{eqnarray} \tag{ 14 }$$

Clearly, there are many potentially tunable parameters involved. Fortunately, we are not necessarily interested in ascertaining the "true" values of these various parameters. Rather, we are concerned with identifying the general shape of the evolutionary trend we'd expect, and what relative differences in the parameters could potentially lead to the differences in D/H evolutionary profiles we find.

In order to find the best fit of this model framework to each of our galaxies, we use a model grid consisting of reasonable values of each of the parameters. The grid parameters are given in Table A1; a full discussion of motivation of the parameter ranges is given in Appendix A.

We assume an average grain radius of $\overline{a}=0.1\,\mu {\rm{m}}$ (Inoue 2011; Nozawa et al. 2011; Asano et al. 2013). For each bin, we took the median of the median dust temperature from the SED fitting of every pixel in the bin, and used this for that bin's T_d value.

For each bin, we also assume a known metallicity. Given the lack of systematic variation in metallicity within the LMC and SMC, as discussed above, we use the fixed values of 0.5 and 0.2 Z_⊙, respectively, for all density bins in these galaxies. For M31 and M33, we also use fixed metallicities for all density bins, using the representative global metallicities, of 1.3 and 0.5 Z_⊙, respectively, even though these galaxies do have metallicity gradients. This is because each bin contains pixels representing a complex distribution of metallicities. For instance, in M31, there are low-density pixels with ${{\rm{\Sigma }}}_{{\rm{H}}}\approx 0.1\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$ located in the galaxy's gas-poor center where Z > 1.5 Z_⊙, and in the galaxy's gas-poor outskirts where Z < 0.75 Z_⊙. The mixture of metallicities at each density varies considerably, and taking the average, for instance, of the metallicities within each bin leads to pathological model behavior. So instead, we use the fixed metallicities, which is in keeping with the understanding that this modeling is intended to be representative, as opposed to accurately capturing the specific physical conditions in these systems.

Using the given ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, Z, T_d , and $\overline{a}$ values for each bin, we performed a χ²-minimizing grid search to find the best-fit parameters for each galaxy. Because the model will ultimately plateau at D/H = Z above a certain density threshold, where all metals are found in dust grains, the model cannot incorporate the turnover in D/H observed for M31, M33, and the LMC; we therefore only fit the model to the points before the turnover in these cases (below 4 M_⊙ pc⁻² for M31 and M33, and below 40 M_⊙ pc⁻² for the LMC). A fuller discussion of this divergence will be conducted in Section 5; however, we present the modeling here first, so that in the following sections we can then explore how to potentially reconcile the data with the models. The resulting best-fit parameter values are given in Table A2, with the models shown on a plot of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ versus D/H in Figure 11.

Zoom In Zoom Out Reset image size

Starting with M31 and M33, we can see in Figure 11 that—up until the turnover—the Asano et al. (2013) model almost perfectly traces the observed evolution in D/H with ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$. For M31, the model matches the observed profile across almost an order of magnitude in both D/H in ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, from 0.6–4 M_⊙ pc⁻², and successfully replicates the increasing steepness of the D/H evolution, followed by a leveling-out. Our data do not probe down to surface densities quite as low as this for M33, but otherwise, the agreement between model and data for ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}\,\lt \,4\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ is similarly excellent.

For the LMC, the agreement between model and data in the 2.2–40 M_⊙ pc⁻² range, before the D/H turnover, is not quite as close a match, but nonetheless does broadly trace the increase in D/H, followed by a leveling-out.

For the SMC, the model does a similarly mixed job of fitting the data—and the lack of a D/H turnover or plateau means that our model can be fit to the full range of densities. The D/H of the SMC over the 8–50 M_⊙ pc⁻² range is somewhat elevated over what is expected from the best-fitting model. But otherwise, the model fits the rest of the profile for the SMC reasonably well, capturing the gentle slope at lower densities, and the steepening at higher densities.

So it would appear that the Asano et al. (2013) model does a good job of capturing the broad strokes of evolution in the D/H. Our data reflect the expectation in the Asano et al. (2013) framework that grain growth gets increasingly efficient as density increases; but only up until where all metals are locked into dust grains, where it plateaus. For M31, M33, and the LMC in particular, we find it to be very reassuring that the levels of the model plateaus are all close to the levels of the peak D/H for each galaxy. Because the D/H plateau indicates when 100% of ISM metals are locked up in dust grains, the D/H at which it occurs is not free to vary in the model, but rather happens at a fixed level of D/H = Z. For instance, in the case of the LMC, the plateau happens at D/H = 0.5 Z_⊙ =0.5 × 0.014 = 7 × 10⁻³, which is very close to the actual highest value of D/H we measure in the LMC, of 5.5 × 10⁻³. The agreement is similarly reasonable for the other galaxies. Given that the values D/H we measure do not "know" the metallicity of their galaxy's ISM, and therefore the D/H at which the model plateau will occur, this agreement suggests that our measured D/H values are sensible.

Where the data and model diverge radically, however, is in the apparent turnover in D/H that three of our galaxies exhibit. That is the focus of the following section.

The most surprising feature in Figure 9 is undoubtedly the fact that, for M31, M33, and the LMC, D/H does not plateau as expected from modeling, but rather turns over and starts to decrease above a certain value of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$.

This is extremely surprising. As ISM density increases, the efficiency of dust grain growth should also increase. The greater the density of the ISM, the more frequently a dust grain will encounter gas-phase metals, and therefore the more frequently that dust grain will increase its mass through accretion of those metals, driving up D/H. Moreover, higher-density ISM will provide superior shielding from the forces of grain destruction, such as shocks and high-energy radiation. Therefore, D/H should continue to increase with density, until grain growth starts to saturate as fewer and fewer metals are left in the gas phase to accrete, ultimately reaching f_d = 1—hence the plateau predicted by the Asano et al. (2013) model.

And there are further reasons to doubt whether this turnover could actually be present for our target galaxies. For instance, at ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ > 150 M_⊙ pc⁻², it appears that the D/H profile for the SMC will intersect the profile of the LMC. And it would be extremely surprising if the SMC were to have a D/H that exceeds that of the LMC at a given ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, given the fact their metallicities differ by a factor of 2.5. Indeed, the peak D/H in the LMC is ≈2.5 times greater than that in the SMC.

In short, it is hard to conceive of any physical mechanism by which D/H could actually fall as density increases. So we examined the question of what could be causing this to appear to be the case in our data. We considered six possible reasons: (1) variations in α_CO; (2) noise-induced anticorrelation; (3) physical resolution effects; (4) dust destruction by supernovae and high-energy radiation in high-density environments due to star formation; (5) the presence of dark gas; and (6) varying dust mass opacity. Our full exploration of these possible explanations is presented in Appendix B. For readers not wishing to explore this investigation in full, the results are summarized as follows:

We are confident that we can rule out overestimation of α_CO, or elevated dust destruction due to environmental effects, as causes of the turnover. It also appears that neither physical resolution limitations nor noise-induced anticorrelation could be causing the turnover, either.

It seems likely that the presence of dark gas could be contributing to the appearance of the D/H bump and turnover for M31 and M33, and driving up D/H at the peak values in the LMC and SMC. However, at the same time, we are confident that dark gas could not be causing significant bias in the overall D/H evolution profile for the LMC, nor causing the turnover. This suggests that dark gas may not be the only effect acting in M31 and M33, either.

The possibility that the dust mass absorption coefficient, κ, has a decreasing value at higher densities provides a viable solution to the apparent turnover, instead changing the high-density portion of the D/H evolution profile into a plateau, as we would expect based on dust evolution modeling. This agrees with the empirical result reported in Clark et al. (2019), but conflicts with predictions from physical dust grain models (in which κ is expected to increase with rising density).

Overall, we consider a combination of dark gas and varying κ to be the most likely explanation, but with the expectation that the contributions of each, and of other possible factors, is likely varying between environments.

A major motivation for this work is the conspicuous disagreement of D/H measurements for the LMC and SMC derived from UV absorption line spectroscopy of elemental depletions, when compared to D/H measurements derived from FIR and radio observations. As mentioned in Section 1, the FIR D/H estimates previously reported for the Magellanic Clouds are much smaller than those determined from depletions. This disagreement is clear in Figure 12, with the previous FIR D/H measurements of Roman-Duval et al. (2017) much lower than the UV D/H measurements found by the Hubble program Metal Evolution, Transport, and Abundance in the Lmc (METAL; Roman-Duval et al. 2021, 2022b). ¹⁸  The Roman-Duval et al. (2017) D/H measurements were performed using FIR data from Planck and IRAS. The disagreement between the D/H determined using the two methods is striking—being a factor of ∼2 for the LMC, and a factor of ∼5 for the SMC.

Zoom In Zoom Out Reset image size

The lower, FIR-derived value suggests the SMC has a significantly depressed D/H relative to its metallicity, and would indicate that the SMC is in the midst of the critical metallicity transition; the higher UV-derived value, however, suggests that the SMC has a D/H reasonably in line with the trend seen at higher metallicities, assuming a roughly constant dust-to-metals ratio. The situation is similar for the LMC; the lower FIR value would hint that the LMC is just starting to experience depressed D/H; the higher UV value would instead place the LMC on the constant dust-to-metals relation.

If the UV-derived D/H values are correct, that would mean that there are severe shortfalls in the dust masses determined for local galaxies, with the vast majority of the dust mass in the SMC being unaccounted for by previous studies such as Chastenet et al. (2017) and Roman-Duval et al. (2017). On the other hand, it is possible that the UV-derived values are overestimates.

While there should be less scope for such significant systematic error in the UV measurements, it is conceivable. For instance, Roman-Duval et al. (2022b) calculated the fractions of hard-to-observe carbon and oxygen (major dust constituents by mass) depleted into the dust phase in the Magellanic Clouds, by using the apparently invariant relation between the depletion of these elements, and other more easily observed elements, with increasing H column (Jenkins 2009). It is possible that this relation becomes unreliable in the Magellanic Clouds, thereby affecting the fraction of metals in the dust phase at given densities.

We also note that directly comparing the Σ_H measurements derived from UV, to those derived from the 21 cm and CO radio data, can be troublesome. The target stars for the UV sightlines will each be located at a different depth in that galaxy's disk. On average, a given sightline will sample half of the thickness of the disk; we could therefore try applying a constant factor of 2 correction to the Σ_H of each UV sightline, to attempt to correct for this. However, the true (unknown) correction will vary wildly between sightlines. This would place us in the undesirable situation of applying a "correction," where the uncertainty on that correction is much larger than the correction itself. Moreover, each UV sightline has the same radius as the target star, and so will therefore only be sampling a pencil beam of ISM tens of millions of kilometers across, in contrast to the tens of parsecs resolution of the 21 cm and CO observations we use (the effects of this difference are analyzed in depth in Section 7 of Roman-Duval et al. 2021). We therefore opt to only apply our standard deprojection correction to the surface densities derived from the UV data, and advise that the reader remains aware that the radio- and UV-derived Σ_H values cannot be compared in a wholly "apples-to-apples" manner.

While the use of FIR all-sky survey data by Roman-Duval et al. (2017) ensures that no diffuse emission was missed (in contrast to the old HERITAGE reductions of the Magellanic Cloud Herschel data), their use of Planck and IRAS data may have limited the accuracy of their results. For instance, they used IRAS data from the Iras Sky Survey Atlas (ISSA; Wheelock et al. 1994), which suffers from a very nonlinear detector response, that varies as a function of both the surface brightness and the angular scale of the emission observed (leading the ISSA explanatory supplement to suggest 100 μm photometric uncertainty of up to 60%; Wheelock et al. 1994), whereas our new feathered Herschel maps are pegged to the absolute calibration of the COsmic Background Explorer (COBE; Boggess et al. 1992). Also, the peak of the dust temperature distribution lies in the > 0.5 dex gap in wavelength coverage between the 100 μm IRAS band and the 350 μm Planck band, likely limiting the accuracy of the Roman-Duval et al. (2017) SED fitting, whereas our data samples this regime with the 160 and 250 μm bands. Plus, unlike us, Roman-Duval et al. (2017) did not allow for a broken β in their SED fitting, which may have compromised their results in areas with significant submillimeter excess (see Figure 6).

For both the LMC and SMC, the Roman-Duval et al. (2017) relationships between D/H and ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ trace similar slopes to our own, over the range we have in common. However, their trends are offset to much lower D/H (despite the fact they use the same κ as we do). This is presumably due to the improvements in our FIR data and SED fitting.

As a result of the consistently larger D/H values we find, there is now much less conflict between our FIR estimates of D/H, and the UV estimates of D/H from Roman-Duval et al. (2021, 2022b). For the LMC, the conflict has essentially been resolved entirely. Indeed, as can be seen in Figures 12 and 13, our D/Hs are now slightly greater than those determined from the UV. To evaluate the scale of the remaining difference, we calculated what offset factor would have to be applied to our D/H profile to minimize the difference between it the UV measurements. For each UV measurement of D/H, we calculated the D/H implied by our profile at that same Σ_H (by interpolating between the points in our profile that bracket the Σ_H of the UV point), and found the difference between them. Having repeated this for every point, we then used a χ²-minimizing routine to find what offset factor, applied to our D/H profile, would minimize the overall difference between the two data sets.

Zoom In Zoom Out Reset image size

For the LMC, we find that reducing our D/H values according to a factor of 0.87 would lead to the closest agreement between them and the UV measurements. Both the original and offset D/H profiles are shown in Figure 13. Moreover, not only does the absolute D/H level now match very well, but the slope of our D/H evolution profile is also an excellent match to that of the UV data. Given the calibration uncertainties on both data sets, this 13% disagreement between the UV and our FIR measurements of D/H is small enough that we feel confident in saying that the discrepancy can be deemed resolved in the case of the LMC.

Exceedingly frustratingly, however, the maximum surface density to which the Hubble data of Roman-Duval et al. (2021, 2022b) could probe D/H corresponds exactly to the location of the D/H turnover. Their UV data does not probe above this density, and hence cannot help us ascertain the nature of the turnover. The Hubble data does not probe to higher densities than this because extinction increases with Σ_H, and above this density it becomes impossible to adequately detect the OB stars used as background sources for this kind of absorption line spectroscopy. Extending UV absorption line spectroscopic depletion measures to higher densities would require a large next-generation UV telescope, such as the Large Ultraviolet Optical Infrared concept (The LUVOIR Team2019).

The fact that the discrepancy in UV versus FIR measures of D/H has been resolved for the LMC makes it all the more interesting, however, that a significant discrepancy persists in the case of the SMC. As can be seen in Figure 14, our new data only partially close the gap between FIR and UV estimates of D/H in the SMC. Previously, the difference was a factor of 5, as found by the Roman-Duval et al. (2017) low-resolution FIR analysis; with our new data, this has been reduced to a factor of 3 (with this offset calculated in the same manner as for the LMC, above).

Zoom In Zoom Out Reset image size

It is not immediately clear why our analysis would resolve the FIR "missing dust" problem for the LMC, but not the SMC. Both use the new Herschel data of Paper I, reduced and feathered in exactly the same way. The processing and analysis also proceed identically throughout. The physical resolution for the SMC is 47 pc, versus 14 pc for the LMC; however, as shown in Appendix B.3, even a factor of 10 degradation in physical resolution causes no systematic shift in the D/H evolution profile. We would therefore not expect a much smaller change in resolution to lead to any significant bias—let alone a bias as large as this.

As previously discussed, it is known that the dust properties in the SMC are fundamentally different from those in higher-metallicity galaxies in several ways. Most sightlines explored in the SMC lack the 2175 Å extinction bump seen in higher-metallicity systems (Gordon et al. 2003; Murray et al. 2019); the SMC has stronger submillimeter excess emission (Bot et al. 2010; Planck Collaboration et al. 2011c; Gordon et al. 2014); and the chemical composition of the dust in the SMC is different, too, with iron and carbon being more dominant than at higher metallicities (Roman-Duval et al. 2022b).

Such large difference in physical dust properties should inevitably lead to some change in κ. A simple explanation for the persistent D/H offset for the SMC is that it has a dust mass absorption coefficient that is a factor of 2.5 smaller than that of the LMC—and therefore also a factor of 2.5 smaller than that of the Milky Way, as the κ₁₆₀ = 1.24 m² kg⁻¹ value we use was itself calibrated using the emission and depletions of the Galactic cirrus Roman-Duval et al. (2017). Unlike the SMC, dust in both the LMC and Milky way exhibits the 2175 Å extinction bump (Gordon et al. 2003), which meshes with the prospect of the SMC dust being different from Milky Way and LMC in other ways, too. However, dust in the SMC is more carbon-rich than the more silicate-rich dust of the Milky Way (and lesser extent of the LMC; Roman-Duval et al. 2022b), and carbon dust should be more emissive than silicate-dominated dust, not less (Ysard et al. 2018). However, changes in grain morphology (such as porosity, size, and shape), for instance, could counteract the differences expected from composition alone.

We do, however, wish to restate that the link between actual volume density (the parameter that will drive grain growth, and hence D/H), and surface density Σ_H (our observable proxy for volume density) potentially have a less direct relationship for the SMC than for the other galaxies of our sample, due to the complex elongation structure of the SMC along our line of sight (Scowcroft et al. 2016). It is possible that this biases our measurements of D/H at different densities, but the specifics would depend entirely on the nature of the three-dimensional structure of the SMC.

However, even with the remaining discrepancy between our D/H for the SMC and the UV value, we can nonetheless rule out a steep drop in D/H at the metallicity of the SMC suggested by the older FIR value from Roman-Duval et al. (2017). This therefore indicates that the SMC is not in the critical metallicity regime.

A final note regarding the D/H offsets: when comparing global D/H, as plotted in Figure 12, the discrepancies between our values and the UV values of Roman-Duval et al. (2021, 2022b) are different than those in Figures 13 and 14 (requiring an offset correction of 0.77 for the LMC, and 2.47 for the SMC). This change is due to the fact that Roman-Duval et al. (2022b) calculated their global D/H values by taking the relationship between Σ_H and D/H over the range of densities sampled by their observations, and extrapolated it to higher and lower Σ_H. They could thereby assign an extrapolated D/H to the full range of Σ_H in the LMC and SMC, and so estimate global D/H. We, on the hand, have measured D/H over the full range of Σ_H in the data (and observe the D/H turnover at densities higher than those probed by the UV measurements, for example) and, hence, the differences in our global D/H values. We consider the offset calculated from our binned D/H profile to be the best indicator of the difference between the two, as it is calculated over the range of densities for which both data sets have actual measurements.

6.1. Comparing Local D/H Residuals between FIR and UV Measurements

In Roman-Duval et al. (2021), the authors note an interesting trend in the residuals between the iron depletion (a proxy for the fraction of metals locked up in dust, which directly correlates with D/H at fixed metallicity) determined from each individual UV spectra, and the depletion that would be expected for that spectra, given that spectra's Σ_H. In general, Roman-Duval et al. (2021) of course found that sightlines that sampled higher Σ_H tended to have greater D/H (as inferred from the amount of iron depletion). However, they also found that in the southeastern portions of the LMC, sightlines would usually have lower D/H than would be expected for their Σ_H; conversely, in the northwestern portions of the SMC, sightlines often had higher D/H than would be expected from their Σ_H. Roman-Duval et al. (2021) suggested that star formation triggered by the gas accretion from the Magellanic Stream onto the southeast of the LMC could be causing an increased amount of grain processing and destruction through radiation, shocks, etc.

We were curious whether the same trend was visible in our data. Using our D/H versus Σ_H relation for the LMC (as in Figure 9), we calculated the D/H we would expect for each pixel in the LMC, given its Σ_H. We then found the residual between the measured D/H in each pixel, and the predicted D/H. We plot a map of these residuals in Figure 15.

Zoom In Zoom Out Reset image size

Like Roman-Duval et al. (2021), we find that D/H is significantly depressed in the southeast of the LMC, relative to what we would expect based on Σ_H alone. Figure 15 shows that this depression is most conspicuous along the edge of the H i tail associated with the gas being accreted from the Magellanic Stream. Notably, the most depressed D/H is found farthest along this tail, to the southern edge of the LMC disk. This is farther from the regions of enhanced star formation being trigged by the infall, centered around 30 Doradus. This suggests that dust processing and destruction due to the effects of star formation is not the dominant cause of the lowered D/H in the southeast of the LMC; otherwise, we would expect D/H to be lower in areas closer to 30 Doradus, and the other areas of heightened star formation. However, grain processing due to gas collision and the resulting shocks, arising from the infall, could be the cause. The absence of a metallicity gradient in the young stars of the LMC argues against a metallicity effect in general (Roman-Duval et al. 2021); however, metallicity will likely be depressed at larger radii, along the tail, where the infalling gas dominates (Nidever et al. 2008; Tsuge et al. 2019, 2020).

The largest positive residuals in Figure 15 seem to be located along the edges of various large low-density features in the LMC (compare to Figure 2), which correspond to known supershells (Meaburn 1980), carved by recent star formation. It is not immediately obvious why this might be the case. It is conceivable that the winds from young OB stars and recent supernovae are able to blow away the less-dense ISM more easily than the denser ISM where more grain growth has occurred, preferentially leaving behind material with a greater dust content. Similarly, elevated D/H at the edge of star-forming regions, due to the evacuation of dust from the centers via wind-driven dust grain drift, is expected based on the work of Draine (2011).

Figure 15 also shows the positions of the individual UV sightlines probed by Roman-Duval et al. (2021), and the residual they found for that sightline. We found the residual at each of those positions in our data by taking the mean of the values in our residual map in a 3 × 3 pixel square aperture centered on the coordinates of each UV sightline. ¹⁹  The Roman-Duval et al. (2021) D/H residuals correlate with those we find at the same positions. This is shown explicitly in Figure 16, which plots our D/H residuals against theirs. While the correlation is not strong, it would appear to be significant. According to a Kendall's Tau rank correlation test (Kendall & Gibbons 1990), the probability of the null hypothesis, of no correlation, is ${ \mathcal P }(\mathrm{null})=0.05$.

Zoom In Zoom Out Reset image size

It is not surprising that there is not an especially tight correlation between the D/H residuals we find, and those of Roman-Duval et al. (2021). Their Hubble spectra sample a pencil beam only as wide as the star being measured, compared to the 15 pc resolution (and therefore 45 × 45 pc aperture) of our data. Plus, each UV spectrum's pencil beam will only sample the ISM on the near side of the star being used—and different stars in their sample will be located at different depths into the LMC disk, as viewed from Earth. Our data, on the other hand, sample the thickness of the entire LMC disk. The fact that our respective residuals do indeed seem to correlate—and that there is most certainly large-scale structure in these residuals, as clearly apparent in Figure 15—demonstrates that there are significant large-scale variations in the ISM properties of the LMC being driven by its ongoing interaction with the SMC through the Magellanic Stream.

Alongside this work, we are releasing the feathered Herschel maps presented in Paper I. These are provided as Flexible Image Transport System (FITS; Wells et al. 1981; Hanisch et al. 2001) files, with one FITS file for each band, for each galaxy. These FITS files contain four extensions. Extension 1 (IMAGE) provides the standard feathered map. Extension 2 (UNC) provides the uncertainty map. Extension 3 (MASK) provides a binary mask map indicating the portion of the data where reliable, fully feathered high-resolution coverage is available. Extension 4 provides the foreground-subtracted version of the feathered map (FGND_SUB), the header of which also describes the uncertainty on that subtraction. All maps are in units of MJy sr⁻¹ (except for the binary masks). Full details of the creation, testing, and properties of these maps can be found in Paper I.

We also provide data products from the analysis performed in this work. For the SED fitting, we provide maps of the median value of each parameter in each pixel (i.e., the maps shown in Figures 5 and 6), and maps of the uncertainties on those medians (being the 68.3% quantile around the median). This data is provided in FITS format; each of these FITS files contain two extensions. Extension 1 (median) provides the map of pixel parameter median values. Extension 2 (uncert) provides the map of uncertainties on those medians. Additionally, we provide the full posterior probability distribution for all SED parameters, consisting of 1000 posterior samples, for all pixels, in the form of an FITS file containing a four-dimensional hypercube, with axes corresponding to R.A., decl., parameters (in order: Σ_d , T_d , β₁, β₂, λ_break, and e₅₀₀), and samples.

We also provide our maps of Σ_H, and of D/H. Note that none of the provided maps have had deprojection corrections applied

Lastly, we provide the Swift-UVOT maps used in Section 4.1.2. This data is provided for Swift-UVOT bands W1, W2, and M2. For each band, we provide a FITS file containing three extensions. Extension 1 (SURF_BRI) provides the map of the surface brightness in MJy sr⁻¹ (converted using the Swift-UVOT zero-points given in Breeveld et al. 2011). Extension 2 (RATE) provides the map of the count rate (in photons per second). Extension 3 (EXP) provides the map of the exposure time (in seconds). The maps for the LMC and SMC are those presented in Hagen et al. (2017). The maps for M31 and M33 were reduced following the same process as those in Hagen et al. (2017), and will be fully presented in M. Decleir et al. (2023, in preparation), but are provided here for the purposes of reproducibility.

This full data set is available on Zenodo at doi:10.5281/zenodo.7392275. The feathered Herschel maps, as presented in Paper I, can also be accessed at the NASA/IPAC Infrared Science Archive at doi:10.26131/IRSA545 (Clark 2021).

In this paper, we have explored the relationship between dust and gas in the Local Group galaxies M31, M33, the LMC, and the SMC, a sample which spans a wide range in mass, metallicity, and other properties. Our investigation has taken advantage of new Herschel maps of these galaxies, presented in Paper I of this series. Previous Herschel data for these galaxies suffered from severe filtering of extended emission (due to the Herschel data reduction process; Meixner et al. 2013; Roussel 2013; Smith et al. 2017, 2021), systematically biasing that data's ability to detect diffuse dust, as well as compromising the scope for accurate foreground subtraction, along with other adverse effects.

The new data from Paper I combined the Herschel maps, in Fourier space, with data from Planck, IRAS, and COBE. Those other telescopes, while having resolution a factor of >10 worse than Herschel, did not filter out the diffuse emission. By merging the data, the Paper I maps preserve the exquisite angular resolution of Herschel, while also recovering the previously missed diffuse dust emission on large scales. By allowing us to probe the widest possible range of physical scales—and hence densities—in the ISM of our target galaxies, the Paper I data have allowed us to investigate how dust evolves in the ISM.

Previous work has found strong evidence that dust grains grow in denser regions of the ISM (Fitzpatrick & Massa 2005; Jenkins 2009; Roman-Duval et al. 2014, 2017, 2021; Tchernyshyov et al. 2015), but the specifics of how much grain growth happens in different environments, and especially at different metallicities, is very much an open question. In particular, galaxies with greater metallicities typically have higher D/H; however, below a certain metallicity, there is evidence that D/H drops sharply. This "critical metallicity" suggests that dust grain growth in the ISM only becomes efficient once metallicity reaches a certain level (Asano et al. 2013; Feldmann 2015; Zhukovska et al. 2016). However, the location and nature of the critical metallicity transition remains unclear.

We used the new Herschel data from Paper I to perform resolved fitting of the FIR dust SEDs of our galaxies. Comparing the resulting maps of dust mass surface density, to maps of the hydrogen surface density Σ_H (derived from 21 cm and CO observations), has allowed us to examine the dust-to-gas ratio (D/G; which we quantified using the dust-to-hydrogen ratio, D/H) over an unparalleled 2.5 orders of magnitude in surface density for our sample of Local Group galaxies, and across the factor ∼6 variation in metallicity they represent. We have therefore been able to study the evolution of the dusty ISM of galaxies to a far greater level of detail than has been previously possible.

Our key findings are as follows:

• 1.  
  The dust-to-gas ratio, D/H, shows very significant evolution with gas surface density, Σ_H, for all four of the galaxies in our sample. We find that D/H increases with density by a factor of 22.4 in the SMC, a factor of 9.0 in M31, a factor of 3.5 in the LMC, and a factor of 2.5 in M33. This is considerably more than found by any previous study of dust-to-gas variation within galaxies.
• 2.  
  Because M31 and M33 have shallow metallicity gradients ²⁰  and because the LMC and SMC have little-to-no metallicity gradient, we can be confident that this evolution in D/H is not simply due to a metallicity effect.
• 3.  
  We examine whether greater efficiency of dust destruction at lower densities could be driving the evolution in D/H. We consider: dust destruction by recently formed stars (from hard radiation, and from the corresponding core-collapse supernovæ) as traced by UV emission; dust destruction due to ionized gas as traced by Hα emission, and models of improved efficiency of dust destruction in lower-density ISM. It does not appear that any of these can account for the observed evolution of D/H with ISM density. On the contrary, we are surprised that greater UV luminosity density, and ionized gas surface density, are correlated with higher D/H.
• 4.  
  In light of the above, our favored explanation for the strong evolution in D/H with Σ_H is that it is being driven by increasingly efficient dust grain growth at higher ISM densities.
• 5.  
  The D/H versus Σ_H evolution profiles of M31 and M33 agree extremely well. The peak D/H for M31 of 0.01 is 20% higher than the peak for M33, but otherwise they follow each other very closely, with integrated D/H differing by only 5%. This is somewhat surprising, as M31 has 2.6 times higher metallicity, and 25 times more stellar mass, than M33.
• 6.  
  Conversely, the large differences between the D/H evolution profiles of M33 and the LMC are very surprising, given these galaxies' close similarity in mass, metallicity, and star formation rate. Nonetheless, the peak D/H of M33 is 45% greater than that of the LMC, and occurs at a surface density a factor of 10 smaller, with the integrated D/H of M33 being 82% greater than the LMC's. This has implications for the common technique (often applied at high redshift) of estimating a galaxy's gas mass from its dust emission; while an observer would likely feel confident in applying the same conversion factor to two galaxies as apparently similar as the LMC and M33, this would not in fact be reliable.
• 7.  
  The D/H evolution profiles of M31, and M33, and the LMC share a confusing trait, whereby after steadily increasing with Σ_H, they then turn over and decrease at higher densities. There is no physical reason to expect this; dust evolution modeling predicts that D/H should plateau at higher densities, not turn over.
• 8.  
  After extensive investigation, we rule out overestimation of α_CO, elevated dust destruction due to star formation, physical resolution effects, and noise-induced anticorrelation, as the cause of the D/H turnover. We find that dark gas (atomic and/or molecular) could cause a turnover to appear for M31 and M33 (and artificially steepen the relationship between D/H and density in the SMC). However, while dark gas is probably elevating D/H in the LMC at intermediate densities, it could not be the driver of the turnover observed for the LMC. We find a fall in the dust mass absorption coefficient (κ) with density to be a plausible explanation for the turnover in the LMC; if this is the case, we would expect this effect to also contribute to the turnover in M33 and M33.
• 9.  
  Previous FIR-based estimates of D/H in the LMC and SMC disagreed with those derived from UV absorption line spectroscopy measurements of elemental depletions. This disagreement was a factor of 2 for the LMC, and a factor of 5 for the SMC, implying the existence of previously missed dust in these galaxies. Our new D/H estimates resolve this tension for the LMC, but only reduce it to a factor of 2.5 disagreement for the SMC. Given the otherwise close agreement in the D/H evolution profiles between the FIR and UV results, we propose this suggests that the dust mass absorption coefficient, κ, is a factor of 3 lower in the SMC than for the sample's other galaxies, causing its D/H measurements (and dust mass) to be underestimated.

This research was improved by the insightful comments of the anonymous referee, whose input improved this work, and spurred additional lines of investigation that proved fruitful.

C.J.R.C. and J.R.-D. acknowledge financial support from the National Aeronautics and Space Administration (NASA) Astrophysics Data Analysis Program (ADAP) grant 80NSSC18K0944.

This research benefited throughout from the constructive, thoughtful, and friendly input (and research environment), provided by the ISM*@ST group, ²¹  whose help made this a better paper; with particular thanks to C. Murray, J. Wu, and C. Zucker.

This research made use of Astropy, ²²  a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018). This research made use of reproject, ²³  an Astropy-affiliated Python package for image reprojection. This research made use of Photutils, ²⁴  an Astropy-affiliated package for detection and photometry of astronomical sources (Bradley et al. 2020). This research made use of NumPy ²⁵  (van der Walt et al. 2011; Harris et al. 2020), SciPy ²⁶  (Jones et al. 2001; Virtanen et al. 2020), and Matplotlib ²⁷  (Hunter 2007). This research made use of the pandas ²⁸  data structures package for Python (McKinney 2010). This research made use of corner, ²⁹  a python package for the display of multidimensional samples (Foreman-Mackey 2016). This research made use of iPython, an enhanced interactive Python (Pérez & Granger 2007).

This research made use of sequential color-vision-deficiency-friendly color maps from cmocean ³⁰  (Thryng et al. 2016) and CMasher ³¹  (van der Velden 2020)

This research made use of TOPCAT ³²  (Taylor 2005), an interactive graphical viewer and editor for tabular data.

This research made use of the SIMBAD database ³³  (Wenger et al. 2000) and the VizieR catalog access tool ³⁴  (Ochsenbein et al. 2000), both operated at CDS, Strasbourg, France. This research has made use of the Nasa/ipac Extragalactic Database ³⁵  (NED), operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. This research has made use of the NASA/IPAC InfraRed Science Archive ³⁶  (IRSA), which is funded by NASA and operated by the California Institute of Technology. This research made use of the HyperLEDA database ³⁷  (Makarov et al. 2014).

Facilities: Herschel - European Space Agency's Herschel space observatory, Planck - , IRAS - , COBE - , Swift - , Parkes - , Arecibo - , VLA - , ATCA - , Effelsberg. -

In Section 4.2, we describe how we model the evolution of D/H with Σ_d , using the model framework of Asano et al. (2013). In Table A1, we give the parameter grid we use to fit the model to our data.

Table A1. D/H Evolution Model Grid Parameter Ranges and Step Sizes, for Our Implementation of the Asano et al. (2013) Model

Parameter	Minimum	Maximum	Step
HΣ⇒n (${\mathrm{cm}}^{-3}\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$)	0.05	5.0	0.02
${f}_{{d}_{0}}$	0.0	0.30	0.002
tgrowth (Myr)	1	100	0.05 dex

Note. The grid for t_growth is spaced logarithmically, with step size therefore given in dex.

Download table as:  ASCIITypeset image

For the hydrogen surface-to-volume density conversion, H_Σ⇒n , we use a grid ranging from 0.05–5.0 ${\mathrm{cm}}^{-3}\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$. For context, assuming the disk of the LMC has a thickness of 100 pc (Elmegreen et al. 2001), and taking from our Σ_H maps a mean hydrogen surface density for the LMC disk of approximately 7 M_⊙ pc⁻², this implies a mean volume density of 2.8 cm⁻³, and hence ${H}_{{\rm{\Sigma }}\Rightarrow n}=0.4\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$. Our 0.05–5.0 ${\mathrm{cm}}^{-3}\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$ grid, spanning over an order of magnitude, is therefore roughly centered on this value in log space. Our best-fit parameters for H_Σ⇒n tend to be a factor of 8–12 greater than this; we postulate that this may reflect the fact that dust grain growth along a given sightline is driven not by the average volume density being sampled, but rather the areas of greater density in particular. H_Σ⇒n will be very dependent upon what fraction of the ISM along a given sightline is found in the denser molecular phase, versus the more vacuous atomic phase. This will be especially true for a galaxy like the SMC, being highly elongated along our line of sight, such that a given sightline could have a very high observed surface density, but still have little to none of that material in environments of greater volume density. In practical terms, increasing H_Σ⇒n has the effect of decreasing the density at which the D/H plateau is reached.

For the minimum fraction of the metal mass locked up in dust grains, ${f}_{{d}_{0}}$, we use a grid ranging from 0.0–0.30. Given that the average value of f_d is approximately 0.5 in high-metallicity spiral galaxies (James et al. 2002; Jenkins 2009; Chiang et al. 2018; Telford et al. 2019), and given that elemental depletions (and therefore f_d ) vary with density by an order of magnitude (Jenkins 2009; Roman-Duval et al. 2021), it is unlikely that the very lowest value of f_d in a galaxy would exceed 0.30. In practical terms, ${f}_{{d}_{0}}$ sets the asymptote in D/H that is reached at low densities. This also has the effect of dictating the steepness of the relationship between Σ_H and D/H. This is because the maximum possible D/H is set by metallicity (being when all metals are in dust), a lower ${f}_{{d}_{0}}$ requires a steeper evolution of D/H with Σ_H (the range of densities from which D/H evolves from its minimum to maximum value is not adjustable by the free parameter). This appears to be the main reason why the LMC requires a larger best-fit value of ${f}_{{d}_{0}}$ than the other galaxies, as it has a shallower evolutionary trend between Σ_H and D/H.

For the average duration of episodes of grain growth, t_growth, we use a grid ranging from 1–100 Myr. If grain growth happens predominantly in molecular clouds, then this corresponds to the average molecular cloud lifespan. The average molecular cloud lifespan is thought to be around 10 Myr, with estimates varying by a factor of a few (Kruijssen et al. 2015; Meidt et al. 2015; Chevance et al. 2020). The fact our data show D/H increasing steadily over a wide range of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ suggests that growth is not just happening in the very densest regions (i.e., not only in molecular clouds), and/or that dust destruction takes place to different degrees over a wide range of densities. Observationally, an increased rate of dust destruction in less-dense environments would manifest in this model framework as a reduction in t_growth. So in practice, t_growth encompasses the typical duration of episodes of grain growth, weighted by how efficiently destruction occurs when grains are not undergoing growth. Ultimately, our best-fit values of t_growth are within a factor of a few of typical estimates of molecular cloud lifetimes. In practical terms, t_growth is degenerate with H_Σ⇒n , with greater values of t_growth decreasing the density at which the D/H plateau is reached. Essentially, this is due to the fact that more grain growth can occur when episodes of grain growth last longer, and/or when the volume density of the ISM is greater.

Despite t_growth being degenerate with H_Σ⇒n , we opt to keep both parameters in the model. The primary reason for this is that both could reasonably be expected to vary a great deal between between galaxies, and we do not want to force an entirely unphysical value of either parameter to be adopted. This degeneracy also means that there are combinations of these parameters, other than those given in Table A2, which give fits that are effectively just as good.

Table A2. Best-fit Parameters for the D/H Evolution Model, for Each of Our Galaxies

Parameter	M31	M33	LMC	SMC
HΣ⇒n (${\mathrm{cm}}^{-3}\,{M}_{\odot }^{-1}\,{\mathrm{pc}}^{2}$)	4.71	4.09	3.37	4.71
${f}_{{d}_{0}}$	0.062	0.084	0.292	0.05
tgrowth (Myr)	45	18	2	6

Note. These are the parameters used for the models plotted in Figure 11.

Download table as:  ASCIITypeset image

The model grid contains just over 1.2 million models. The best-fit parameters we found for each galaxy are provided in Table A2. We stress again that we do not suggest that these are the best estimates for the "true" values of these parameters. But rather, that they are intended to be reasonable values that yield models that fit the data as well as possible, and indicate the sort of D/H evolutionary profile we should expect.

Here we present our in-depth investigation into the possible causes of the apparent turnover in the D/H at higher values of Σ_H for the LMC, M31, and M33, as presented in Section 5. Specifically, we explore six possible explanations for the apparent turnover: (1) differences in α_CO; (2) noise-induced anticorrelation; (3) physical resolution effects; (4) dust destruction by star formation; (5) the presence of dark gas; and (6) varying dust mass opacity.

B.1. Variations in α_CO?

If we were to have overestimated α_CO for our galaxies, then we would be overestimating the amount of molecular gas present, which would mean that we would get erroneously low D/H values—specifically at the higher densities where molecular gas represents a larger faction of the total dust budget.

Very suggestively, the surface densities at which molecular gas starts to dominate over atomic gas in our ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ data are 2.1, 3.6, and 35 M_⊙ pc⁻² for M31, M33, and the LMC, respectively—all very close to the densities at which their turnovers happen, being ≈4 M_⊙ pc⁻² for M31, M33, and ≈40 M_⊙ pc⁻² for the LMC. Meanwhile, for the SMC, which does not have a turnover, there is likewise no surface density at which the measured molecular gas surface density exceeds that of the atomic gas. ³⁸

We examined to what degree α_CO would need to be reduced, in order to get rid of the turnover. To do this, we repeated our D/H analysis, instead using values of α_CO that were modified by factors of 0.5, 0.2, and 0.1, relative to the fiducial values we use for each galaxy (namely 3.2 K⁻¹ km⁻¹ s pc⁻² for M31 and M33, 6.4 K⁻¹ km⁻¹ s pc⁻² for the LMC, and 21 K⁻¹ km⁻¹ s pc⁻² for the SMC; see Section 7.1 of Paper I for specifics on our choice of α_CO values). The resulting D/H evolutionary profiles are shown in Figure B1.

Zoom In Zoom Out Reset image size

The plots in Figure B1 show that reducing α_CO by over a factor 2–5 is able to mostly remove the turnover in the case of the LMC. But even a factor of 10 reduction in α_CO cannot remove the turnover for M31 or M33; indeed in the case of M33, the turnover becomes more distinct. This may seem counterintuitive, but the reasons for it are twofold.

First, reducing α_CO will increase the D/H measured for a given pixel, while decreasing the ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$. This will have the effect of moving values upward and leftward on Figure B1, therefore tending to cause measurements to move along the post-turnover portion of a profile, as opposed to changing the shape of a profile.

Second, for the factor 0.2 and 0.1 changes to the fiducial α_CO, we are effectively "removing" the vast majority of the molecular gas mass from these galaxies. This has the effect of making it very hard to accurately distinguish the relative gas surface densities of pixels. As an example, consider two pixels, for which the molecular gas content has been calculated using the fiducial α_CO: an atomic-gas-dominated pixel with Σ_{H I} =1.6 and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}=0.4;$ and a molecular-gas-dominated pixel with Σ_{H I} = 1.3 and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}=3.0$. Here, the molecular-dominated pixel has a considerably larger total mass, while both pixels have similar atomic-gas content, due to any higher-density gas tending to enter the molecular phase. Now consider the same pair of pixels when α_CO is multiplied by 0.1. Now, the originally atomic-dominated pixel will have Σ_{H I} = 1.6 and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}=0.04$, while the formerly molecular-dominated pixel will have Σ_{H I} = 1.3 and ${{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}=0.3$, meaning that the two pixels would now be measured as having near-identical gas content. In short, because differences in atomic-gas content between pixels become relatively smaller at higher densities, significantly reducing α_CO makes high-density pixels appear to have similar gas content to one another. As a result, whereas higher dust densities were previously very strongly associated with pixels that have higher gas densities, that relationship is now more mixed. The worse this mixing gets, the more the trend in D/H versus ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ will tend toward a gradient of −1, in a classic instance of noise-induced anticorrelation, ³⁹  which is indeed what we see happening for M33 in Figure B1.

Given that significantly reducing α_CO gives rise to this sort of spurious effect in our data, it would seem that reducing α_CO by such larger factors is likely unphysical, which should not be surprising. For instance, reducing the α_CO of the LMC by a factor of 5 would make it only 1.28 K⁻¹ km⁻¹ s pc⁻²—less than half the standard Milky Way value of 3.2 K⁻¹ km⁻¹ s pc⁻², and much less than would be expected for a galaxy with less than half solar metallicity (Bigiel et al. 2010; Gratier et al. 2010; Druard et al. 2014). Additionally, reducing α_CO does not even manage to get rid of the D/H turnover, although it is possible that it may minimize the turnover partially in the case of the LMC, even for less-aggressive reduced values of α_CO.

As such, we are confident that overestimation of α_CO is not the cause of the D/H turnover observed for our galaxies.

B.2. Noise-induced Anticorrelation?

As seen in Appendix B.1, noise-induced anticorrelation drove negative gradients in D/H evolution when using very low values of α_CO. So it is clearly worth checking whether noise-induced anticorrelation could also be giving rise to the turnover in the first place. To test this, we constructed simulated versions of our data, to explore how noise with different behaviors can affect the trends we observe.

To start with, we created a noiseless toy model in which D/H evolves from 0.001 at ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}=1\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, to 0.01 at ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}=10\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, at which D/H then plateaus. This toy model is plotted in Figure B2, and is intended to be the simplest possible version of the increase-then-plateau trend in dust evolution we might expect (Asano et al. 2013).

Zoom In Zoom Out Reset image size

First, we assess the impact that generally large scatter could have when observing this trend. We created simulated data for which 10⁶ data points, evenly distributed in ${\mathrm{log}}_{10}({{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})})$, were drawn from the toy model over a ${10}^{-2}\lt {{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}\lt {10}^{5}$ interval. ⁴⁰  We then applied 1 dex of Gaussian noise to the ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ and Σ_d of each randomly drawn value, recomputed D/H for these, then placed them in 0.025 dex bins, as per the real data in Figure 9, etc. A noise level of 1 dex should do a good job of probing the effects of significant scatter on the observed trend. ⁴¹  The points generated by introducing this scatter are shown as blue squares in Figure B2, and as can be seen, they still trace the underlying trend. The knee at ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}=10\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ is no longer as sharp, having been smoothed out by the scatter, but otherwise there are no adverse effects present, and certainly no sign of a spurious D/H turnover.

Next, we perform a test designed to more closely simulate the scatter present in our actual data, which tends to steadily decrease for bins at higher densities. ⁴²  We proceed as above, but now, instead of being a constant 1 dex, the scale of the Gaussian scatter on both ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ and Σ_d is a function of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, falling as density increases. Scatter in dex, σ, is given by $\sigma =-0.5{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})})+1.0$, with a floor value of σ = 0.1 dex imposed at the highest densities where σ would otherwise fall lower than this. The result of this simulation is shown with the green triangles in Figure B2. Introducing this form of scatter behavior does not appear to introduce any pathologies into the D/H evolution profile (other than the knee once again being smoothed out). This is reassuring, given that this should roughly emulate the way the scatter evolves with density in the real data.

Lastly, we test a "plausible worst-case scenario" for the evolution of scatter with density. For this, we hold the scale of the Gaussian noise at a constant low level of 0.05 dex in ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ and 0.2 dex in Σ_d for values up to ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}\,=10\,{M}_{\odot }\,{\mathrm{pc}}^{-2};$ above this, the noise grows sharply with Σ_H according to a power law, increasing to 0.15 dex in ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ and 1.2 dex in Σ_d at ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}=1000\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$. By keeping noise low before the knee, where D/H is increasing, then making noise increase considerably after the knee (i.e., in the density regime where we observe the turnover in the real data), we should be maximizing the potential for noise-induced anticorrelation, or some other artifact, to create a spurious trend in D/H at higher densities. The binned medians for this test are shown with pink triangles in Figure B2. Even with this test, we find no suspicious behavior in the resulting D/H evolution profile.

Following these tests, we are aware of no way in which noise-induced anticorrelation could give rise to an artificial D/H turnover in our data.

B.3. Physical Resolution Effects?

The physical resolution our data is able to achieve for M31 and M33 (137 and 147 pc, respectively) is up to 10 times worse that what it can achieve for the LMC and SMC (15 and 47 pc, respectively). It is therefore noteworthy that the two galaxies with the much poorer physical resolution also exhibit the sharper turnovers in D/H. So it is conceivable that limited physical resolution could give rise to spurious trends in our results. In particular, temperature mixing could cause the dust mass determined from SED fitting to be biased low. The impact of temperature mixing will be expected to increase as physical resolution worsens and as density increases—because a single pixel may then contain not only the cold, dense dust in giant molecular clouds, but also warm dust being heated by recent-formed stars, all blended together. This might cause dust mass to be underestimated progressively more as density increases, driving down D/H in a manner resembling our observed turnover.

To test this possibility, we reprocessed our Herschel and ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ data for the LMC and SMC, degrading the observations so that the Magellanic Clouds appeared as they would were they at the same distance as M31 and M33. Specifically, we degraded the data for the LMC and SMC to produce a physical resolution of 142 pc—the equivalent of 36'' angular resolution at a distance of 815 kpc. As 815 kpc is the midpoint between the distances to M31 and M33, this maximizes our ability to compare fairly to both. ⁴³  If limited physical resolution is causing spurious trends in our results, then degrading the physical resolution of the LMC and SMC by a factor of ∼10 should lead to significant changes in their D/H evolution profiles.

However, as can be seen in Figure B3, the D/H versus ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ profiles followed by the LMC and SMC are almost completely unchanged by being degraded to the same physical resolution as M31 and M33. The only significant difference is that, as would be expected, the degraded data cannot trace to surface densities as high as the full-resolution data. ⁴⁴  Otherwise, the degraded data simply traces the same D/H evolution profile as the original data, albeit with somewhat more scatter.

Zoom In Zoom Out Reset image size

We are in fact pleasantly surprised by how well the D/H evolution profiles for the LMC and SMC have been preserved. In particular, reducing physical resolution from 14 pc to 142 pc, in the case of the LMC, will lead to vastly more temperature mixing, and other merging of dust populations within each pixel (Galliano et al. 2011). The fact that dust properties can still be recovered well enough to leave the D/H evolution profile effectively unchanged is a better outcome than we might have expected. This gives us further confidence that temperature mixing is not significantly biasing the results of our SED fits, in addition to the lack of residuals described in Section 3.1.

The D/H evolution profiles for the LMC and SMC are not at all biased by degrading their data to the effective distance of M31 and M33. This suggests that the observed profiles for M31 and M33 are similarly not significantly biased versus how they would appear were they observed at a closer distance to us. We therefore conclude that the effects of physical resolution limits are probably not causing the appearance of the D/H turnover in the case of these two galaxies.

For the LMC, the degraded data is not able to probe to the highest densities, and only traces up to ≈ 60 M_⊙ pc⁻². We therefore cannot make a strong statement either way about the possibility of a spurious turnover being caused by temperature mixing at densities greater than this. However, we note that of all four galaxies we consider, the SMC should be the one vulnerable to the greatest impact from temperature mixing. The reason for this is that its extreme elongation along the line of sight significantly increases the likelihood of different dust populations, heated to different temperatures by different environmental conditions, being present along a shared column, in a single pixel. However, instead of displaying a turnover, the SMC in fact shows the D/H relation getting steeper at the highest densities. This remains true even for the degraded data in Figure B3, where we would expect the impact of any temperature mixing to be exaggerated, thanks to each pixel sampling 9× more area. This gives us some additional reason to think that temperature mixing is not seriously biasing our D/H values at the highest densities.

B.4. Dust Destruction by Environmental Effects at Higher Densities?

Various processes can destroy dust grains in the ISM. Many of these processes arise due to star formation. For instance, by the supernova shocks following the deaths of recently formed massive stars, or by direct photodestruction of grains by UV and X-ray photons, or by thermal sputtering of grains in the hot ionized gas produced by young massive stars (Bocchio et al. 2014; Slavin et al. 2015; Jones et al. 2017; Galliano et al. 2018).

It is therefore conceivable that, in areas of higher-density ISM where more star formation occurs (Kennicutt 1998), the dust destructive processes associated with recent star formation lead to D/H being depressed, manifesting as the turnover we observe. If this is the case, then for a given ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ we would expect to find lower D/H in areas with more star formation, than in areas with less star formation at the same ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$.

To test this, we use the Hα maps of Gaustad et al. (2001), which cover both the LMC and SMC at a resolution of 48'', well matched to the resolution of our dust and gas data. We use these Hα maps to create maps of the surface density of ionized gas, ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$, following the prescription of Paradis et al. (2011). These maps of ionized gas are an ideal way of tracing where dust destruction due to star formation is likely to be happening. Not only is ionized gas a proxy for star formation, but the ionized gas itself is the environment in which the resultant dust destruction will occur.

We repeated our analysis of D/H versus ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ for the LMC and SMC, splitting the data into six bins of ionized gas surface density, from ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}={10}^{-1.5}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ to ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}\,={10}^{1.5}\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$, with each bin having a width of 0.5 dex. In Figure B4, the D/H evolution profile within each of these ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$ bins is plotted, for the LMC and SMC. Note that for these plots, ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, including the denominator of D/H, has been revised to also include the contribution of ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$, to ensure internal consistency. ⁴⁵  We also plot the global D/H evolution profile for both galaxies in Figure B4, incorporating the ionized gas into ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$.

Zoom In Zoom Out Reset image size

Figure B4 shows that the D/H evolution profile behaves very similarly for regions of high ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$ as it does in regions of low ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$. The only exception to this is for the regions of the very highest ionized gas surface density, 10^1.0 < Σ_H+< 10^1.5 M_⊙ pc⁻², where in both the LMC and SMC, the D/H is conspicuously elevated above the general trend. This argues strongly against the hypothesis that the effects of nearby star formation could be depressing D/H. And most importantly of all, the D/H turnover in the LMC is still present for all of the bins of ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$ that trace that range of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ (although the profile for the highest ${{\rm{\Sigma }}}_{{\rm{H}}+}^{(\mathrm{deproj})}$ bin does get noisy at Σ_H >100 M_⊙ pc⁻²).

Moreover, Figure 10 and the corresponding discussion in Section 4.1.2 also show that more intense UV radiation fields are associated with higher D/H, not lower, in all four of our sample galaxies.

These observations strongly suggest that dust destruction due to elevated star formation in regions of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$ is not the cause of the D/H turnover.

B.5. Dark Gas?

B.6. Varying Dust Mass Opacity?

The apparent turnover in D/H can be explained if we are systematically underestimating the dust mass at higher densities. This could occur if the intrinsic luminosity of the dust changes at different densities—i.e., if the mass absorption coefficient, κ, varies as a function of ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$.

It is essentially guaranteed that κ varies with density. We know that the composition of dust varies considerably with density, with different elements depleting from the gas phase onto dust grains at different rates in different environments (Jenkins 2009; Jenkins & Wallerstein 2017; Roman-Duval et al. 2021, 2022b). Similarly, the dust grain-size distribution is believed to evolve with density (Hirashita & Kobayashi 2013; Aoyama et al. 2020), as is the structure of the grains, as they accrete layers and coagulate in different phases of the ISM (Cuppen & Herbst 2007; Jones et al. 2016, 2017; Jones 2018). In light of all of this, it is nearly inconceivable that κ would not evolve with density to some degree.

For our SED fitting, we assumed a constant value of κ. If, however, κ were to decrease at higher ${{\rm{\Sigma }}}_{{\rm{H}}}^{(\mathrm{deproj})}$, then this would lead us to underestimate D/H at those higher densities, and could cause a plateau in D/H to instead manifest as a turnover.

Is it plausible that κ decreases in higher-density ISM? The exact nature of how κ evolves with density is mostly unconstrained at present. However, models consistently predict that κ should increase in higher-density ISM, due to grains developing a "fluffy" structure as they coagulate (Ossenkopf & Henning 1994; Li & Lunine 2003; Jones 2018). On the other hand, Clark et al. (2019) used an empirical method to construct resolved maps of κ within nearby spiral galaxies M74 and M83, and found evidence that κ decreases with increasing Σ_H. Relatedly, Bianchi et al. (2019, 2022) found that higher ISM surface densities and higher molecular-to-atomic-gas ratios are associated with greater dust surface brightness per unit gas surface density (resolved within nine nearby spiral galaxies), and with greater dust luminosity per gas mass (for integrated measurements of 204 late-type galaxies), also suggesting dust becomes less emissive at higher ISM densities.

Specifically, Clark et al. (2019) found that κ falls with Σ_H according to a power-law index of −0.4. Because M73 and M83 were relatively poorly resolved in that study (590 and 330 pc resolution, respectively, compared to 14–147 pc in this work), and because those galaxies are highly star-forming (with correspondingly high gas surface densities), the lowest Σ_H to which Clark et al. (2019) could probe in M74 and M83 was 3 M_⊙ pc⁻² (with the highest well-sampled Σ_H being 303 M_⊙ pc⁻²). The much greater distances to M74 and M83 (as compared to the Local Group) also make it likely that temperature mixing will be influencing the trends observed in those galaxies (Priestley & Whitworth 2020). Nonetheless, the Σ_H range for which Clark et al. (2019) reported a fall in κ with Σ_H corresponds to the surface densities above which we start to see the D/H turnovers in our target galaxies.

We examined whether a κ that falls with increasing density, as suggested by Clark et al. (2019), could explain our apparent D/H turnover. Specifically, we tried a broken κ model, where κ remains constant at κ₁₆₀ = 1.24 m² kg⁻¹ (see Section 3) up until a break density, after which κ falls with Σ_H according to a −0.4 power law, as per Clark et al. (2019). For the location of the κ break in each galaxy, we used the surface density at which D/H peaks: Σ_H = 4 M_⊙ pc⁻² for M31 and M33, and Σ_H = 40 M_⊙ pc⁻² for the LMC. The results of this are plotted in Figure B5.

Zoom In Zoom Out Reset image size

In Figure B5, we see that introducing this broken κ model does an excellent job of resolving the D/H turnover for M31, M33, and the LMC. In each case, the D/H evolution profile at higher densities now conforms well to the plateau we would expect, based on dust evolution modeling (see Section 4.2). We therefore think it is certainly plausible that decreasing in κ at higher Σ_H could be the cause of the D/H turnover. It is particularly interesting that the −0.4 power-law rate of decrease of κ versus Σ_H, as found in Clark et al. (2019), is also the rate of decrease required to resolve the turnover in our D/H evolution profiles. That said, the fact that physical dust grain models consistently predict the opposite behavior—that κ is expected to increase with density—should count against this hypothesis. Nonetheless the fact that this explanation works so well makes us consider it very plausible. ⁴⁷

Why the Σ_H at which the D/H turnover happens—and hence the Σ_H at which κ breaks—would be an order of magnitude greater in the LMC than in M31 and M33 is unclear, as is the lack of break in the SMC. We tested changing the break density for the LMC to Σ_H = 4 M_⊙ pc⁻², to match M31 and M33, but this pushes the D/H plateau for the LMC to a level of D/H > 0.015 even higher than the peak D/H in M31, in excess of the metals available in the LMC for forming dust, so this is an unlikely scenario. As discussed in Appendix B.1, the different turnover surface densities do correspond to the point at which the ISM transitions from atomic- to molecular-gas dominated in these galaxies, which happens in LMC at a Σ_H that is an order of magnitude higher than in M31 and M33. So there are processes in the ISM of these galaxies happening at the different densities at which their turnovers happen.

Separately, in Section 6, we find evidence for a constant offset in κ in the SMC, relative to the LMC. While separate from any variation in κ with density, this nonetheless further highlights the fact that we really should not expect κ to remain constant throughout.