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

Christopher J. R. Clark, Julia C. Roman-Duval, Karl D. Gordon, Caroline

Bot, Matthew W. L. Smith

► To cite this version:

Christopher J. R. Clark, Julia C. Roman-Duval, Karl D. Gordon, Caroline Bot, Matthew W. L. Smith. 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 Astrophysical Journal, 2023, 946, 10.3847/1538-4357/acbb66. insu-04056556

HAL Id: insu-04056556

https://insu.hal.science/insu-04056556

Submitted on 6 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

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

Christopher J. R. Clark<fnref n="1" />, Julia C. Roman-Duval<fnref n="1" />, Karl D. Gordon<fnref n="1" />, Caroline Bot<fnref n="2" />, Matthew W. L. Smith<fnref n="3" />, and Lea M. Z. Hagen<fnref n="4" />

Received 2022 May 27; revised 2022 December 3; accepted 2023 February 9; published 2023 March 27

## Abstract

We explore evolution in the dust-to-gas ratio with density within four well-resolved Local Group galaxies—the LMC, SMC, M31, and M33. We do this using new Herschel maps, which restore extended emission that was missed by previous Herschel reductions. Combining this sensitivity to diffuse dust emission with excellent physical resolution allows us to probe the dust-to-gas ratio across 2.5 orders of magnitude in interstellar medium (ISM) surface density. We find a significant increase in the dust-to-gas ratio with density, with the dust-to-gas ratio varying within each galaxy by up to a factor 22.4, as density changes. We explore several possible reasons for this, and our favored explanation is that it is being driven by dust grain growth in denser regions of the ISM. We find that the evolution of the dust-to-gas ratio with ISM surface density is very similar between M31 and M33, despite their large differences in mass, metallicity, and star formation rate; conversely, we find M33 and the LMC to have very different dust-to-gas evolution profiles, despite their close similarity in those properties. Our dust-to-gas ratios address previous disagreement between UV- and far-IR-based dust-to-gas estimates for the Magellanic Clouds, removing the disagreement for the LMC, and considerably reducing it for the SMC—with our new dust-to-gas measurements being factors of 2.4 and 2.0 greater than the previous far-IR estimates, respectively. We also observe that the dust-to-gas ratio appears to fall at the highest densities for the LMC, M31, and M33; this is unlikely to be an actual physical phenomenon, and we posit that it may be due to a combined effect of dark gas, and changing dust mass opacity.

*Unified Astronomy Thesaurus concepts: Dwarf galaxies (416); Local Group (929); Interstellar dust (836); Far infrared astronomy (529); Submillimeter astronomy (1647)*

## 1. Introduction

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; Tcherny-shyov 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

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218-2463, USA; cclark@stsci.edu</footnote>

<footnote n="2">²Observatoire Astronomique de Strasbourg, Université de Strasbourg, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="3">³School of Physics and Astronomy, Cardiff University, Queen's Buildings, The Parade, Cardiff, Wales, CF24 3AA, UK</footnote>

<footnote n="4">⁴Unaffiliated.</footnote>
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<fnref n="5" /> (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).<fnref n="6" /> 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).

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 fraction 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

<footnote n="5">⁵ This relationship is plotted later, in Figure 12, to display results presented in subsequent sections of this work.</footnote>

<footnote n="6">⁶ DLAs are neutral gas absorption systems with $N_{H I} > 2 \times 10^{20} \text{ cm}^{-2}$, found over a wide range of redshifts, and observed via quasar absorption spectroscopy (Wolfre et al. 2005; De Cia et al. 2016).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Three-color images of the galaxies in our target sample, showing hydrogen gas in red (as traced by 21 cm and CO emission), Herschel-SPIRE 350 µm cooler dust emission in green, and Herschel-PACS 100 µm warmer dust emission in blue. The Herschel images are foreground-subtracted maps; see Section 2. The LMC, the best-resolved galaxy of the sample, is shown at the top. The lower row shows the SMC (left), M31 (center), and M33 (right). Each channel of each image uses a logarithmic color scale that displays the map's structure over its full value range. Simple visual inspection of the changing colors within these images make clear that there is significant evolution in the dust and gas properties in each of our sample galaxies.</figcaption></figure>

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.

## 2. Sample and Data

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.

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

<table><caption>**Table 1**

Basic Properties of the Local Group Galaxies Studied in This Work</caption><thead><tr><th></th><th>M31</th><th>M33</th><th>LMC</th><th>SMC</th></tr></thead><tbody><tr><td>α (J2000)</td><td>10°69</td><td>23°46</td><td>80°89</td><td>13°16</td></tr><tr><td>δ (J2000)</td><td>+41°27</td><td>+30°66</td><td>-69°76</td><td>-72°80</td></tr><tr><td>Distance (kpc)</td><td>790</td><td>840</td><td>50</td><td>62</td></tr><tr><td>Hubble Type</td><td>SAb</td><td>SAcD</td><td>SBm</td><td>Irr</td></tr><tr><td>R<sub>25</sub></td><td>89'</td><td>32'</td><td>323'</td><td>151'</td></tr><tr><td>R<sub>25</sub> (kpc)</td><td>20.5</td><td>7.5</td><td>5.0</td><td>2.5</td></tr><tr><td><sup>a</sup>Pos. Angle</td><td>35°</td><td>23°</td><td>170°</td><td>45°</td></tr><tr><td><sup>a</sup>Axial Ratio</td><td>2.57</td><td>1.70</td><td>1.17</td><td>1.66</td></tr><tr><td>Inclination</td><td>77°</td><td>56°</td><td>26°</td><td>...</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td><sup>b</sup>10<sup>11.1</sup></td><td><sup>c</sup>10<sup>9.7</sup></td><td><sup>d</sup>10<sup>9.4</sup></td><td><sup>e</sup>10<sup>8.3</sup></td></tr><tr><td><sup>f</sup> Z (Z<sub>☉</sub>)</td><td><sup>g</sup>1.3</td><td><sup>h</sup>0.5</td><td><sup>i</sup>0.5</td><td><sup>j</sup>0.2</td></tr></tbody></table>

**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.

<sup>a</sup> From the HyperLEDA database (Makarov et al. 2014).

<sup>b</sup> Tamm et al. (2012).

<sup>c</sup> Corbelli et al. (2014).

<sup>d</sup> van der Marel (2006).

<sup>e</sup> Rubele et al. (2018).

<sup>f</sup> 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.4$R_{25}$ (Pilyugin et al. 2004; Moustakas et al. 2010).

<sup>g</sup> 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.

<sup>h</sup> Magrini et al. (2016).

<sup>i</sup> Toribio San Cipriano et al. (2017).

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 Photo-
detector 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.

## 2.3. Molecular and Atomic-gas Data

For the gas component of our target galaxies' ISM, we
created maps of the hydrogen surface density, Σ<sub>H</sub>, by adding
together maps of the atomic and molecular hydrogen, produced
from H I 21 cm hyperfine line and <sup>12</sup>C<sup>16</sup>O(1–0) rotational line
(hereafter CO) observations, respectively. In this work, we use
the same Σ<sub>H</sub> maps as described in Section 7.1 of Paper I. Here
we briefly recap the input data, and how the Σ<sub>H</sub> maps were
produced.

The data we used to create the maps of Σ<sub>H</sub> are as follows: for
the LMC, the H I data of Kim et al. (2003) and CO data of
Wong et al. (2011); for the SMC, the H I data of Stanimirovic
et al. (1999) and the CO data of Mizuno et al. (2001); for M31
the H I data of Braun et al. (2009) and the CO data of Nieten
et al. (2006); and for M33, the H I data of Koch et al. (2018)
and the CO map of Gratier et al. (2010) and Druard et al.
(2014). The atomic and molecular gas emission on all scales
should be captured by these observations; the CO data is all
from single-dish observations, and the H I data is a feathered
combination of high-resolution interferometric observations
and low-resolution single-dish observations. Our handling of
limiting resolutions is discussed in Section 2.4.

<figure><img src="image_5.png" /></figure>

**Figure 2.** Upper: feathered Herschel-SPIRE 250 µm maps of the LMC (left) and SMC (right), as produced by Paper I; these images show the foreground-subtracted versions of the maps. Lower: Σ<sub>H</sub> maps for the LMC (left) and SMC (right), made by combining atomic and molecular gas surface density maps produced from H I and CO data, respectively; these maps are shown at the limiting resolutions at which we conduct our analyses (see Section 2.4; FWHM shown by white circles in the lower-right corners).

<figure><img src="image_6.png" /></figure>

**Figure 3.** Same as for Figure 2, but depicting M31 (left) and M33 (right).

We calculated the molecular gas surface density using the
standard relation Σ<sub>H2</sub> = α<sub>CO</sub> I<sub>CO(1-0)</sub>, where I<sub>CO(1-0)</sub> is the
CO(1-0) line velocity-integrated main-beam brightness temp-
erature (in K km s<sup>-1</sup>), and α<sub>CO</sub> is the CO-to-H<sub>2</sub> conversion
factor (in K<sup>-1</sup> km<sup>-1</sup> s M<sub>☉</sub> pc<sup>-2</sup>). For the spirals M31 and M33,
we used the standard Milky Way α<sub>CO</sub> value of
3.2 K<sup>-1</sup> km<sup>-1</sup> s pc<sup>-2</sup>; for the LMC and SMC dwarf galaxies,
we used higher α<sub>CO</sub> values of 6.4 and 21 K<sup>-1</sup> km<sup>-1</sup> s pc<sup>-2</sup>,
respectively (Bolatto et al. 2013).
As described in Section 7.1 of Paper I, some authors prefer an $\alpha_{CO}$ for M33 that is higher than the Milky Way value (Bigiel et al. 2010; Druard et al. 2014). However, Gratier et al. (2010) found a value within 10% of that of the Milky Way. Plus Rosolowsky et al. (2003)—who do not measure an absolute $\alpha_{CO}$, only its apparent relative variation—found that $\alpha_{CO}$ in the highest-metallicity regions of M33, at > 1.4 $Z_{⊙}$, is not systematically different from that in regions at <0.5 $Z_{⊙}$; given that $\alpha_{CO}$ at such high metallicities should be comparable to that of the Milky Way, it suggests this is not significantly different elsewhere in M33. Regardless, Appendix B.1 shows that our primary results are mostly insensitive to changes in $\alpha_{CO}$.

To add together a galaxy’s maps of atomic and molecular hydrogen, to create the combined $\Sigma_{H}$ map, we first convolved the higher-resolution map to the resolution of the lower-resolution map, assuming a Gaussian point-spread function (PSF), as per the beam size information in the header of each map.

Our $\Sigma_{H}$ maps for each galaxy are shown in the lower panels of Figure 2 for the LMC and SMC, and Figure 3 for M31, and M33.

## 2.4. Convolution to Limiting Resolution

For all analyses, the Herschel and $\Sigma_{H}$ maps were convolved to a common resolution, to match whatever the worst resolution was among all of the data for each galaxy. We created conversion kernels for this with the Python package `photutils`. For the input Herschel PSFs, we used the azimuthally averaged beams from Aniano et al. (2011).<fnref n="7" /> For the input $\Sigma_{H}$ PSFs, we assumed Gaussian beams, as per Section 2.3.

For M31 and M33, the limiting resolution was the 36" of the 500 µm maps; for the LMC and SMC, the limiting resolutions were dictated by the gas observations (see Section 2.3 above), being 1' for the LMC, and 2.6' for the SMC. We also rebin all maps to use a pixel width equal to that galaxy’s limiting resolution, so that every pixel is statistically independent.

These limiting angular resolutions correspond to physical resolutions of 14 pc for the LMC, 47 pc for the SMC, 137 pc for M31, and 147 pc for M33.

## 3. SED Fitting

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(\lambda)$, of dust emission at a given wavelength, $\lambda$, can be expressed by:

$$S(\lambda) = \kappa(\lambda) B(\lambda, T_d) \Sigma_d \quad (1)$$

where $B$ is the Planck function evaluated at wavelength $\lambda$ and dust temperature $T_d$, $\Sigma_d$ is the dust mass surface density, and $\kappa(\lambda)$ is the dust mass absorption coefficient at wavelength $\lambda$;

$\kappa(\lambda)$ varies with wavelength according to an emissivity law:

$$\kappa(\lambda) = \frac{\kappa(\lambda_{\text{ref}})}{\lambda_{\text{ref}}^{-\beta}} \lambda^{-\beta} \quad (2)$$

where $\kappa(\lambda_{\text{ref}})$ is the value of $\kappa(\lambda)$ at a reference wavelength $\lambda_{\text{ref}}$, and $\beta$ 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.<fnref n="8" /> The BEMBB model achieves this by having the value of $\beta$ change at some break wavelength $\lambda_{\text{break}}$; a shallower $\beta$ at longer wavelengths captures the excess emission. As such, the BEMBB emissivity law takes the form:

$$\kappa(\lambda) = \frac{\kappa(\lambda_{\text{ref}})}{\lambda_{\text{ref}}^{-\beta}} E(\lambda) \quad (3)$$

where

$$E(\lambda) = \begin{cases} \lambda^{-\beta_1} & \text{for } \lambda < \lambda_{\text{break}} \\ \lambda^{\beta_2-\beta_1} \lambda^{-\beta_1} & \text{for } \lambda \ge \lambda_{\text{break}} \end{cases} \quad (4)$$

for which $\beta$ takes a value of $\beta_1$ at wavelengths $<\lambda_{\text{break}}$, while it takes a value of $\beta_2$ at wavelengths $\ge\lambda_{\text{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 $\kappa_{160} = 1.24 \text{ m}^2 \text{ kg}^{-1}$, 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 $\beta_2$ directly. Rather, the break in $\beta$ is parameterized via the 500 µm excess, $e_{500}$, which gives the relative excess in the 500 µm flux, above what would

<footnote n="7">⁷ https://www.astro.princeton.edu/~draine/Kernels.html</footnote>

<footnote n="8">⁸ As a test, we repeated the main analyses we present later in this work, but using a single-MBB model instead of BEMBB. We found that all of the trends in D/H we discuss later in this work also appear when using single-MBB fitting, albeit with much worse quality fits in the lower-density and lower-metallicity environments.</footnote>
<table><caption>**Table 2**
SED Model Grid Parameter Ranges and Step Sizes, Used for Pixel-by-pixel SED Fitting with *DustBFF*, for Our BEMBB Model</caption><thead><tr><th>Parameter</th><th>Minimum</th><th>Maximum</th><th>Step</th></tr></thead><tbody><tr><td>Σ<sub>d</sub> (M<sub>☉</sub> pc<sup>-2</sup>)</td><td>10<sup>-6</sup></td><td>10<sup>1</sup></td><td>0.05 dex (12%)</td></tr><tr><td>T<sub>d</sub> (K)</td><td>10</td><td>70</td><td>0.04 dex (5.9%)</td></tr><tr><td>β<sub>1</sub></td><td>0</td><td>3.5</td><td>0.175</td></tr><tr><td>λ<sub>break</sub> (μm)</td><td>125</td><td>525</td><td>33.3</td></tr><tr><td>e<sub>500</sub></td><td>-0.5</td><td>2.0</td><td>0.2</td></tr></tbody></table>

**Note.** For the logarithmically spaced parameters $Σ_d$ and $T_d$, we also give the percentage difference between grid steps.

arise from a standard MBB model (with negative values indicating a 500 µm deficit); $e_{500}$ is thus defined:

$$e_{500} = \left( \frac{\lambda_{\text{break}}}{500 \ \mu\text{m}} \right)^{\beta_2 - \beta_1} - 1 \quad (5)$$

*DustBFF* uses a full covariance matrix in its model evaluations. This covariance matrix, *C*, is given by:

$$C = C_{\text{calib}} + C_{\text{instr}} \quad (6)$$

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

The calibration covariance matrix $C_{calib}$ is calculated by multiplying the proposed fluxes for a given model by the relative calibration uncertainty matrix $U_{calib}$, which is itself given by summing $U_{uncorr}$ and $U_{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, $C_{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 $C_{instr}$ are all zero. The diagonal elements of $C_{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

<table><caption>**Table 3**
Diagonal Elements of $C_{instr}$ Instrumental Noise Covariance Matrix for Our SED Fitting (Off-diagonals Are All Zero)</caption><thead><tr><th>Galaxy</th><th>100 µm</th><th>150 µm</th><th>250 µm</th><th>350 µm</th><th>500 µm</th></tr></thead><tbody><tr><td>LMC</td><td>11.98</td><td>13.85</td><td>2.64</td><td>0.91</td><td>0.43</td></tr><tr><td>SMC</td><td>10.51</td><td>9.30</td><td>1.23</td><td>0.45</td><td>0.29</td></tr><tr><td>M31</td><td>6.25</td><td>4.24</td><td>2.03</td><td>0.71</td><td>0.39</td></tr><tr><td>M33</td><td>6.06</td><td>4.01</td><td>2.11</td><td>0.75</td><td>0.44</td></tr></tbody></table>

**Note.** All values are in map units of MJy sr$^{-1}$.

therefore assume the uncorrelated component is the quadrature subtraction of the correlated uncertainty from the total uncertainty, being $\sqrt{7\%^2-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 $U_{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

$$U_{\text{corr}} = \begin{bmatrix} 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{bmatrix} \quad (7)$$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** SED for example pixel in M33 located at $\alpha = 23^{\circ}6486$, $\delta = 30^{\circ}6592$. The observed surface brightness values from our feathered Herschel data are plotted, along with the 1000 model SEDs from each set of DustBFF posterior sample parameters, and the DustBFF maximum a posteriori (MAP) fit SED.</figcaption></figure>

and

$$ U_{uncorr} = \begin{bmatrix} 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{bmatrix} \quad (8) $$

to give a final $U_{calib}$ of

$$ U_{calib} = \begin{bmatrix} 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{bmatrix} \quad (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.

When deciding what pixels of our Herschel data to fit, we imposed a signal-to-noise ratio (S/N) threshold dictated by the $\Sigma_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 $\Sigma_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 $\Sigma_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 $\Sigma_H$ maps. Because our $\Sigma_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 $\Sigma_H$ maps are very closely matched in sensitivity, with the S/N > 4 threshold corresponding to $\Sigma_H > 2.8 M_{\odot} pc^{-2}$ for the LMC, $\Sigma_H > 2.6 M_{\odot} pc^{-2}$ for M31, and $\Sigma_H > 2.1 M_{\odot} pc^{-2}$ 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' 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.

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.<fnref n="9" />

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 $\beta_1$ and $\beta_2$ 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

<footnote n="9">⁹ Several authors, such as Draine et al. (2014) and Chastenet et al. (2017), have done resolved SED fitting of galaxies in our sample using physical dust models, with different parameter sets to the MBB-based model we use. Therefore, other than $\Sigma_d$, we do not have parameters in common that can be compared directly.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Results of our pixel-by-pixel SED fitting for our target galaxies, with the median value for each parameter shown in each pixel. This figure shows Σ<sub>d</sub>, T<sub>d</sub>, and λ<sub>break</sub>; Figure 6 shows β<sub>1</sub>, β<sub>2</sub>, and ε<sub>500</sub>. SED fitting was only performed for pixels where our Σ<sub>H</sub> maps had an S/N > 4. For each parameter, the same color scale is used for all galaxies, to allow for direct comparison.</figcaption></figure>

mixing to cause β₂ to be significantly flatter than β₁, the colder
dust components must have a Wein's Law peak at a
wavelength > λ<sub>break</sub>. However, >94% of SMC pixels have
λ<sub>break</sub> > 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
<figure><img src="image_9.png" /><figcaption>**Figure 6.** As per Figure 5, except for SED parameters $\beta_1$, $\beta_2$, and $e_{500}$. Identical color scales are used for both $\beta_1$ and $\beta_2$.</figcaption></figure>

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 $\beta$ 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.
<table><caption>**Table 4**
Dust and Gas Parameters for Each of Our Galaxies</caption><thead><tr><th></th><th>M31</th><th>M33</th><th>LMC</th><th>SMC</th></tr></thead><tbody><tr><td>M<sub>d</sub> (10<sup>6</sup> M<sub>☉</sub>)</td><td>45.0</td><td>5.41</td><td>1.35</td><td>0.165</td></tr><tr><td>T<sub>d</sub> (K)</td><td>15.9</td><td>18.4</td><td>21.3</td><td>22.0</td></tr><tr><td>β<sub>1</sub></td><td>2.12</td><td>1.81</td><td>1.73</td><td>1.49</td></tr><tr><td>β<sub>2</sub></td><td>1.96</td><td>1.57</td><td>1.54</td><td>1.09</td></tr><tr><td>λ<sub>break</sub> (μm)</td><td>277</td><td>283</td><td>272</td><td>283</td></tr><tr><td>e<sub>500</sub></td><td>0.21</td><td>0.28</td><td>0.28</td><td>0.41</td></tr><tr><td>D/H</td><td>10<sup>-2.11</sup></td><td>10<sup>-2.13</sup></td><td>10<sup>-2.39</sup></td><td>10<sup>-3.38</sup></td></tr></tbody></table>

**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 $\Sigma_H$ map. For the other parameters, we give the dust-mass-weighted average of the pixel medians.

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 $\beta = 1.8$ with no break at longer wavelengths. Their $\Sigma_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 $\Sigma_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 $\Sigma_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.<fnref n="10" /> The fact that we allow $\beta$ 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.

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),<fnref n="11" /> 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-$\Sigma_H$ region we do not fit. Our results also agree with the finding from Smith et al. (2012) and Whitworth et al. (2019) that $\beta$ decreases with radius in M31. Indeed, whereas both of those studies used an unbroken-$\beta$ model, we find that both $\beta_1$ and $\beta_2$ 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 $\beta$ 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 $\beta_1$ generally fall with radius, in line with the findings of Tabatabaei et al. (2014) —who also used a single-MBB, unbroken-$\beta$ model—we find that $\beta_2$ 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 $\Sigma_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 $\Sigma_H$ S/N criterion.<fnref n="12" />

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 $\Sigma_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 $\kappa(\lambda_{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 $\kappa(\lambda_{ref})$).

The fact we tend to 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,

<footnote n="10">¹⁰ We do not fit the SED of many of the pixels in the gas-poor center of M31, where dust will be warmest, because they fell below the S/N threshold in $\Sigma_H$. However, outside this central region there is an extensive area where both we and Utomo et al. (2019) modeled the SED, where we still find ≈5 K difference.</footnote>

<footnote n="11">¹¹ Smith et al. (2012) and Whitworth et al. (2019) used an unbroken variable $\beta$ for their MBB models, while Viaene et al. (2014) adopted a fixed-$\beta$ MBB model.</footnote>

<footnote n="12">¹² For M31, the application of the $\Sigma_H$ S/N cut does mean that there are pixels in the center of the galaxy for which the dust content is not incorporated into Table 4, due to their low hydrogen column density. However as these pixels comprise <2% of the total dust mass of M31, their omission has minimal impact on the global values.</footnote>
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 $\beta = 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 $\beta = 1.8$ method), it is possible for flux at longer wavelengths to be accounted for by the $\beta$ 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 $\kappa(\lambda_{\text{ref}})$ corrections.<fnref n="13" /> 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

<caption>**Table 5**</caption>

Spearman Rank Correlation Coefficients between the Submillimeter Excess at 500 µm, $e_{500}$, and Various Other Parameters, for the Pixels in Our Maps

<table><thead><tr><th>$e_{500}$ versus</th><th>M31</th><th>M33</th><th>LMC</th><th>SMC</th></tr></thead><tbody><tr><td>$\Sigma_{\text{H}}$</td><td>-0.54</td><td>-0.16</td><td>-0.47</td><td>-0.39</td></tr><tr><td>Deprojected Radius</td><td>0.41</td><td>0.74</td><td>0.57</td><td>...</td></tr><tr><td>D/H</td><td>-0.80</td><td>-0.71</td><td>-0.57</td><td>-0.60</td></tr><tr><td>$\Sigma_d$</td><td>-0.77</td><td>-0.65</td><td>-0.67</td><td>-0.68</td></tr></tbody></table>

Note. All relations have $P_{\text{null}} < 10^{-16}$.

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 $\beta_2 < \beta_1$. While our model grid does allow for $\beta_2 \ge \beta_1$ (i.e., a submillimeter deficit), this only happens for $<3.7\%$ of pixels in the median $e_{500}$ 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 \mu m < \lambda_{\text{break}} < 325 \mu m$ for 80% of pixels, with $\lambda_{\text{break}}$ generally being lower in areas of greater ISM density, albeit with elevated $\lambda_{\text{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 $\Sigma_{\text{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.

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_{500}$ with D/H and $\Sigma_d$ (Table 5). We caution against over-interpretation of these final two correlations, however, as both depend upon $\Sigma_d$, which is a parameter in our SED fitting model along with $e_{500}$, 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_{500}$ tracing the spiral arms (see Figure 6), despite these being higher-density regions. If shielding in dust-rich regions can drive down submillimeter

<footnote n="13">¹³ The ratio of silicate to carbonaceous dust is different in every pixel modeled by Chastenet et al. (2017); as per their Section 5.2, we assume an average 2:1 ratio of carbonaceous to silicate dust, which agrees with the range of average ratios they quote for both galaxies. Using the average THEMIS emissivity slope of $\beta = 1.78$ (Nersesian et al. 2019) to convert via Equation (2), this gives an average $\kappa_{160} = 2.84 \text{ m}^2 \text{ kg}^{-1}$. We compare to their favored multi-ISRF masses.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Corner plot, showing the marginalized posterior distribution of each parameter, and their covariances, for an example pixel in M33 located at $\alpha = 23^\circ 6486, \delta = 30^\circ 6592$.</figcaption></figure>

excess, then it is conceivable that proximity to radiation and shocks from ongoing star formation elevates $e_{500}$ 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.

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.

## 4. Evolution in the Dust-to-gas Ratio

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.

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.<fnref n="14" /> 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 $\Sigma_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_\odot$ metallicity. Thanks to our new Herschel maps from Paper I, and the sensitive $\Sigma_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^{-4}$ in the most diffuse regions.

### 4.1. Evolution of D/H with $\Sigma_H$

In Figure 9, we present the central result of this paper—the evolution of D/H with $\Sigma_H$ for our galaxies. As covered in Section 1, we should expect D/H to increase with greater $\Sigma_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 $\Sigma_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 $\Sigma_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).

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.<fnref n="15" /> 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

<footnote n="14">¹⁴ As we discuss later in Section 6.1, with relation to Figure 15, D/H at the edge of these supershells is even elevated above what would be expected for their $\Sigma_H$, based on the LMC's global relationship between D/H and $\Sigma_H$.</footnote>

<footnote n="15">¹⁵ The same correction has also been applied to the x-axis $\Sigma_H$ values in Figures 11, 13, 14, B1, B2, B3, B4, and B5.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 8. Maps of the dust-to-hydrogen ratio for our target galaxies. All four are displayed using the same color scale, to allow for direct comparison.</figcaption></figure>

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 $\Sigma_{\text{H}}$. Therefore, for a given galaxy inclination $i$ (see Table 1), we apply a deprojection correction:

$$ \Sigma_{\text{H}}^{(\text{deproj})} = \cos(i) \Sigma_{\text{H}}. \qquad (10) $$

Hereafter, we refer to $\Sigma_{\text{H}}^{(\text{deproj})}$ in instances where comparison between galaxies makes the distinction between corrected and uncorrected $\Sigma_{\text{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 $\Sigma_{\text{H}}^{(\text{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_{\odot}$ pc$^2$, while the SMC exhibits even more than that, with its D/H increasing by a factor of 20 between 3 and 150 $M_{\odot}$ pc$^2$. 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 $\Sigma_{\text{H}}$ (this is explored in Section 5). Regardless, the increase in D/H with $\Sigma_{\text{H}}^{(\text{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
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Plot of dust-to-hydrogen ratio D/H against deprojected hydrogen surface density $\Sigma_{H}^{(deproj)}$, with values from individual pixels binned together into 0.025 dex wide bins of $\Sigma_{H}$, for each our of sample galaxies. Points indicate bin medians, with error bars show uncertainty on those medians. Values of $\Sigma_{H}$ have had a deprojection correction applied to account for each galaxy's inclination.</figcaption></figure>

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 \text{ dex } R_{25}^{-1}$ (Zurita & Bresolin 2012);<fnref n="16" /> so clearly this cannot account for the 0.9 dex evolution in D/H in M31. Moreover, metallicity does not scale monotonically with $\Sigma_{H}$ in M31, thanks to the low ISM density of the inner $\sim 3 \text{ kpc}$—further reducing the potential for the radial metallicity gradient to give rise to increased D/H in regions of greater $\Sigma_{H}$.

The ISM of M33 has an even smaller radial metallicity gradient, at $-0.22 \text{ 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

<footnote n="16">¹⁶ The metallicity gradient in Zurita & Bresolin (2012) was reported in terms of dex kpc⁻¹; we have converted it to dex $R_{25}^{-1}$ to ease comparison between galaxies for the reader, as per the $R_{25}$ value in Table 1.</footnote>
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 $\Sigma_{H}^{(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_{\odot} pc^{-2}$. 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 $P_{null} < 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

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Plots of D/H against deprojected Swift-UVOT W2 luminosity surface density, for each of our four sample galaxies. Binned median values are plotted with black crosses.</figcaption></figure>

corresponding higher density of supernovæ shocks and hard radiation, still manage to have higher D/H.

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\alpha$ observations of the LMC and SMC, does not appear to be associated with reduced D/H, and any $\Sigma_{H}^{(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 $\Sigma_{H}^{(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 $\Sigma_{H}^{(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_{\odot} pc^{-2}$ range in $\Sigma_{H}^{(deproj)}$, they overlap
nearly perfectly.<fnref n="17" /> 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⁻².¹¹ versus 10⁻².¹³—matching their similarity in the D/H versus Σ<sub>H</sub><sup>(deproj)</sup> 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 Σ<sub>H</sub><sup>(deproj)</sup> 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 Σ<sub>H</sub><sup>(deproj)</sup> 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 supernovae can destroy dust grains, given the density and metallicity of the ISM.

Using this model, the timescale for accretion of metals onto dust grains, τ<sub>acc</sub>, is given by:

$$ \tau_{\text{acc}} = 2 \times 10^7 \left( \frac{\bar{a}}{0.1 \, \mu\text{m}} \right) \left( \frac{n_{\text{H}}}{100 \, \text{cm}^{-3}} \right)^{-1} \times \left( \frac{T_d}{50 \, \text{K}} \right)^{-\frac{1}{2}} \left( \frac{Z}{0.02} \right)^{-1} \quad (11) $$

where ¯a is the average dust grain size in microns, n<sub>H</sub> is the volumetric number density of hydrogen in cm⁻³, and Z is the metallicity (in units of the metal mass fraction). With τ<sub>acc</sub>, one

<footnote n="17">¹⁷ We are able to probe to lower densities for M31 than M33, primarily because of M31’s greater inclination. Surface densities that would not be detectable in a less-inclined galaxy can be observed in M31, as the sampled ISM column is longer at higher inclination. Our deprojection corrections, discussed earlier in Section 4.1, allow us to compare galaxies directly despite this.</footnote>
can then compute the fraction of the metal mass that is in dust grains (see Zhukovska et al. 2008), using:

$$f_d = \frac{e^{t_{\text{growth}}/\tau_{\text{acc}}} f_{d_0}}{1 - f_{d_0} + (f_{d_0} e^{t_{\text{growth}}/\tau_{\text{acc}}})} \quad (12)$$

where $t_{\text{growth}}$ is the average duration of episodes of grain growth (if grain growth happens predominantly in molecular clouds, then $t_{\text{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:

$$\frac{D}{H} = f_d Z \quad (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_{\Sigma \to n}$, in units of cm$^{-3}$ M$_{\odot}^{-1}$ pc$^2$, such that:

$$n_H = H_{\Sigma \to n} \times \Sigma_H^{(\text{deproj})}. \quad (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 $\bar{a} = 0.1$ µ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$_{\odot}$, 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$_{\odot}$, 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 $\Sigma_H \approx 0.1 M_{\odot}^{-1}$ pc$^2$ located in the galaxy’s gas-poor center where $Z > 1.5 Z_{\odot}$, and in the galaxy’s gas-poor outskirts where $Z < 0.75 Z_{\odot}$. 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.

<figure><img src="image_14.png" /><figcaption>Figure 11. Plot of D/H against $\Sigma_H^{(\text{deproj})}$. Data points are the same as those shown in Figure 9, but now with the addition of the best-fit dust evolution model for each galaxy. The parameters that give the model for each galaxy are provided in Table A2.</figcaption></figure>

Using the given $\Sigma_H^{(\text{deproj})}$, Z, $T_d$, and $\bar{a}$ values for each bin, we performed a $\chi^2$-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_{\odot}^{-1}$ pc$^{-2}$ for M31 and M33, and below $40 M_{\odot}^{-1}$ pc$^{-2}$ 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 $\Sigma_H^{(\text{deproj})}$ versus D/H in Figure 11.

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 $\Sigma_H^{(\text{deproj})}$. For M31, the model matches the observed profile across almost an order of magnitude in both D/H in $\Sigma_H^{(\text{deproj})}$, from 0.6–4 M$_{\odot}$ pc$^{-2}$, 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 $\Sigma_H^{(\text{deproj})} < 4 M_{\odot}^{-1}$ pc$^{-2}$ is similarly excellent.

For the LMC, the agreement between model and data in the 2.2–40 M$_{\odot}$ pc$^{-2}$ 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$_{\odot}$ pc$^{-2}$ 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.

**5. The Nature of the D/H Turnover**

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 $\Sigma_{\text{H}}^{(\text{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 $\Sigma_{\text{H}}^{(\text{deproj})} > 150 M_{\odot} \text{ pc}^{-2}$, 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 $\Sigma_{\text{H}}^{(\text{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 $\alpha_{\text{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 $\alpha_{\text{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, $\kappa$, 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 $\kappa$ is expected to increase with rising density).

Overall, we consider a combination of dark gas and varying $\kappa$ 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.

**6. The Discrepancy between FIR to UV Dust-to-gas Ratios**

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).<fnref n="18" /> 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.

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

<footnote n="18">¹⁸ Not plotted on Figure 12 are even earlier estimates of D/H from Roman-Duval et al. (2014). These used the older HERITAGE maps of the Magellanic Clouds, and employed the SED fitting results presented by Gordon et al. (2014). As discussed in Section 3.1.1, the Gordon et al. (2014) dust masses are biased high, due to the large amount of missing flux in the HERITAGE Herschel-PACS 100 and 160 µm data, causing their SED fitting to output erroneously high dust temperatures.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The relationship between D/H and Z, showing observations and models from various sources. The large polygonal points indicate D/H values for the four galaxies in our sample, from this work (hexagons), from the previous low-resolution FIR measurements of Roman-Duval et al. (2017; wide diamonds), and the UV absorption line measurements of elemental depletions from the Hubble METAL program of Roman-Duval et al. (2022b; narrow diamonds). The small circular points indicate values from the sample of De Vis et al. (2019; which incorporates and standardizes measurements of galaxies from the samples of Rémy-Ruyer et al. 2014; Clark et al. 2015, 2018; Davies et al. 2017; De Vis et al. 2017a), the largest sample to date exploring this parameter space. The dotted line shows a trend of the D/H ∝ Z (i.e., a constant dust-to-metals ratio), passing through Z = Z⊙, D/H = 0.01. The thick gray lines show a selection of dust evolution models from Feldmann (2015), with a different location for each model of the “critical metallicity” at which D/H sharply increases with Z; specifically, we plot the Feldmann (2015) “equilibrium models” where the location of the critical metallicity is set by the ratio of molecular gas depletion timescale to the interstellar grain growth timescale, with lower ratios leading to higher critical metallicities (see their Figure 3).</figcaption></figure>

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 Σ<sub>H</sub> 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 Σ<sub>H</sub> 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
Σ<sub>H</sub> 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
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Plots of D/H against $\Sigma_H$ for the LMC, comparing the measurements reported at various densities from several sources. Plotted in black circles are D/H values from UV absorption line measurements of elemental depletions by Roman-Duval et al. (2022b). Plotted in pink diamonds are D/H values from IRAS–Planck FIR data by Roman-Duval et al. (2017). The values from this work are shown for each galaxy, both before and after having an offset in $\kappa$ applied, to make our FIR values match the UV values of Roman-Duval et al. (2022b). Error bars on our points are omitted for clarity, but would be the same as in Figure 9.</figcaption></figure>

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 $\beta$ in their SED fitting, which may have compromised their results in areas with significant submillimeter excess (see Figure 6).

For both LMC and SMC, the Roman-Duval et al. (2017) relationships between D/H and $\Sigma_H^{(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 $\kappa$ 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 $\Sigma_H$ (by interpolating between the points in our profile that bracket the $\Sigma_H$ of the UV point), and found the difference between them. Having repeated this for every point, we then used a $\chi^2$-minimizing routine to find what offset factor, applied to our D/H profile, would minimize the overall difference between the two data sets.

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

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Plots of D/H against $\Sigma_H$ for the SMC, comparing the measurements reported at various densities from several sources. Details are as per Figure 13, except for the SMC instead of the LMC.</figcaption></figure>

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 $\Sigma_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 Team 2019).

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).

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 $\kappa$. 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 $\kappa_{160} = 1.24 \text{ m}^2\text{kg}^{-1}$ 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 $\Sigma_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 $\Sigma_H$ and D/H over the range of densities sampled by their observations, and extrapolated it to higher and lower $\Sigma_H$. They could thereby assign an extrapolated D/H to the full range of $\Sigma_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 $\Sigma_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 $\Sigma_H$. In general, Roman-Duval et al. (2021) of course found that sightlines that sampled higher $\Sigma_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 $\Sigma_H$; conversely, in the northwestern portions of the SMC, sightlines often had higher D/H than would be expected from their $\Sigma_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 $\Sigma_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 $\Sigma_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.

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 $\Sigma_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 triggered 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
<figure><img src="image_18.png" /><figcaption>**Figure 15.** A map of the residuals between the D/H measured for each pixel in the LMC, and the D/H we would expect for each pixel (as calculated using each pixel's Σ<sub>H</sub>, and the relationship between Σ<sub>H</sub> and D/H in Figure 9). Red areas have a higher D/H than would be expected given their Σ<sub>H</sub>; blue pixels have lower D/H than would be expected. The circles indicate the UV absorption spectroscopy sightlines of Roman-Duval et al. (2021); the color with each circle indicates the residual between iron depletion (a proxy for D/H) calculated by Roman-Duval et al. (2021) for that sightline based on its UV spectra, and the depletion they predicted based on the relationship they found between depletion and Σ<sub>H</sub>.</figcaption></figure>

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.<fnref n="19" /> 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 *P*(null) = 0.05.

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.

7. Data Products

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.

<footnote n="19">¹⁹ We used a 3 × 3 pixel square aperture in order to sample multiple beams along each axis, and therefore reduce noise, while still sampling a small enough region to be comparable to the value at the specific location of the UV sightline.</footnote>
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Plot of the residuals found by Roman-Duval et al. (2021), between the iron depletion they measure for a given UV sightline (a proxy for D/H), and the depletion they would expect given that sightline's Σ<sub>H</sub>—plotted against the residual between the D/H we find for that position of that sightline in our data, as compared to the D/H we would expect given that position's Σ<sub>H</sub> (given the relationship between Σ<sub>H</sub> and D/H we find in Figure 9). The dotted line shows the 1:1 relation.</figcaption></figure>

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: Σ<sub>d</sub>, T<sub>d</sub>, β<sub>1</sub>, β<sub>2</sub>, λ<sub>break</sub>, and e<sub>500</sub>), and samples.

We also provide our maps of Σ<sub>H</sub>, 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).

8. Conclusion

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; Tcherny-shyov 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 Σ<sub>H</sub> (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, Σ<sub>H</sub>, 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<fnref n="20" /> 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 Σ<sub>H</sub> is that it is being driven by increasingly efficient dust grain growth at higher ISM densities.

5. The D/H versus Σ<sub>H</sub> 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 Σ<sub>H</sub>, 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 α<sub>CO</sub>, 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,<fnref n="21" /> whose help made this a better paper; with particular thanks to C. Murray, J. Wu, and C. Zucker.

This research made use of *Astropy*,<fnref n="22" /> a community-developed core *Python* package for Astronomy (Astropy Collaboration et al. 2013, 2018). This research made use of *reproject*,<fnref n="23" /> an *Astropy*-affiliated *Python* package for image reprojection. This research made use of *Photutils*,<fnref n="24" /> an *Astropy*-affiliated package for detection and photometry of astronomical sources (Bradley et al. 2020). This research made use of *NumPy*,<fnref n="25" /> (van der Walt et al. 2011; Harris et al. 2020), *SciPy*,<fnref n="26" /> (Jones et al. 2001; Virtanen et al. 2020), and *Matplotlib*,<fnref n="27" /> (Hunter 2007). This research made use of the *pandas*,<fnref n="28" /> data structures package for *Python* (McKinney 2010). This research made use of *corner*,<fnref n="29" /> 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*,<fnref n="30" /> (Thryng et al. 2016) and *CMasher*,<fnref n="31" /> (van der Velden 2020).

This research made use of *TOPCAT*,<fnref n="32" /> (Taylor 2005), an interactive graphical viewer and editor for tabular data.

<footnote n="20">²⁰ The metallicity of M31 and M33 falls by only a factor of 3.6 per R<sub>25</sub> for M31, and a factor of 1.65 for M33—much less than the evolution in D/H.</footnote>

<footnote n="21">²¹ https://www.ismstar.space/</footnote>

<footnote n="22">²² https://www.astropy.org/</footnote>

<footnote n="23">²³ https://reproject.readthedocs.io</footnote>

<footnote n="24">²⁴ https://photutils.readthedocs.io</footnote>

<footnote n="25">²⁵ https://numpy.org/</footnote>

<footnote n="26">²⁶ https://scipy.org/</footnote>

<footnote n="27">²⁷ https://plotlib.org/</footnote>

<footnote n="28">²⁸ https://pandas.pydata.org/</footnote>

<footnote n="29">²⁹ https://corner.readthedocs.io</footnote>

<footnote n="30">³⁰ https://plotlib.org/cmocean/</footnote>

<footnote n="31">³¹ https://cmasher.readthedocs.io</footnote>

<footnote n="32">³² http://www.star.bris.ac.uk/~mbt/topcat/</footnote>
This research made use of the SIMBAD database<fnref n="33" /> (Wenger et al. 2000) and the VizieR catalog access tool<fnref n="34" /> (Ochsenbein et al. 2000), both operated at CDS, Strasbourg, France. This research has made use of the NASA/IPAC Extragalactic Database<fnref n="35" /> (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<fnref n="36" /> (IRSA), which is funded by NASA and operated by the California Institute of Technology. This research made use of the HyperLEDA database<fnref n="37" /> (Makarov et al. 2014).

*Facilities:* Herschel, Planck, IRAS, COBE, Swift, Parkes, Arecibo, VLA, ATCA, Effelsberg.

## Appendix A

### D/H Evolution Model Parameters

In Section 4.2, we describe how we model the evolution of D/H with $\Sigma_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.

For the hydrogen surface-to-volume density conversion, $H_{\Sigma \to n}$, we use a grid ranging from $0.05-5.0 \text{ cm}^{-3} M_{\odot}^{-1} \text{ pc}^2$. For context, assuming the disk of the LMC has a thickness of 100 pc (Elmegreen et al. 2001), and taking from our $\Sigma_H$ maps a mean hydrogen surface density for the LMC disk of approximately $7 M_{\odot}^{-2}$, this implies a mean volume density of $2.8 \text{ cm}^{-3}$, and hence $H_{\Sigma \to n} = 0.4 M_{\odot}^{-1} \text{ pc}^2$. Our $0.05-5.0 \text{ cm}^{-3} M_{\odot}^{-1} \text{ 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_{\Sigma \to 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_{\Sigma \to 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_{\Sigma \to 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 $\Sigma_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}$

<caption>**Table A1**</caption>

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

<table><thead><tr><td>Parameter</td><td>Minimum</td><td>Maximum</td><td>Step</td></tr></thead><tbody><tr><td>$H_{\Sigma \to n}$ (cm<sup>-3</sup> M<sub>\odot</sub><sup>-1</sup> pc<sup>2</sup>)</td><td>0.05</td><td>5.0</td><td>0.02</td></tr><tr><td>$f_{d_0}$</td><td>0.0</td><td>0.30</td><td>0.002</td></tr><tr><td>$t_{\text{growth}}$ (Myr)</td><td>1</td><td>100</td><td>0.05 dex</td></tr></tbody></table>

**Note.** The grid for $t_{\text{growth}}$ is spaced logarithmically, with step size therefore given in dex.

<caption>**Table A2**</caption>

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

<table><thead><tr><td>Parameter</td><td>M31</td><td>M33</td><td>LMC</td><td>SMC</td></tr></thead><tbody><tr><td>$H_{\Sigma \to n}$ (cm<sup>-3</sup> M<sub>\odot</sub><sup>-1</sup> pc<sup>2</sup>)</td><td>4.71</td><td>4.09</td><td>3.37</td><td>4.71</td></tr><tr><td>$f_{d_0}$</td><td>0.062</td><td>0.084</td><td>0.292</td><td>0.05</td></tr><tr><td>$t_{\text{growth}}$ (Myr)</td><td>45</td><td>18</td><td>2</td><td>6</td></tr></tbody></table>

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

requires a steeper evolution of D/H with $\Sigma_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 $\Sigma_H$ and D/H.

For the average duration of episodes of grain growth, $t_{\text{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 $\Sigma_H^{(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_{\text{growth}}$. So in practice, $t_{\text{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_{\text{growth}}$ are within a factor of a few of typical estimates of molecular cloud lifetimes. In practical terms, $t_{\text{growth}}$ is degenerate with $H_{\Sigma \to n}$, with greater values of $t_{\text{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_{\text{growth}}$ being degenerate with $H_{\Sigma \to 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.

<footnote n="33">³³ https://simbad.u-strasbg.fr/simbad/</footnote>

<footnote n="34">³⁴ https://vizier.u-strasbg.fr/viz-bin/VizieR</footnote>

<footnote n="35">³⁵ https://ned.ipac.caltech.edu/</footnote>

<footnote n="36">³⁶ https://irsa.ipac.caltech.edu</footnote>

<footnote n="37">³⁷ http://leda.univ-lyon1.fr/</footnote>
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.

## Appendix B
### Possible Causes of the D/H Turnover

Here we present our in-depth investigation into the possible causes of the apparent turnover in the D/H at higher values of $\Sigma_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 $\alpha_{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 $\alpha_{CO}$?

If we were to have overestimated $\alpha_{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 $\Sigma_H^{(deproj)}$ data are 2.1, 3.6, and $35 M_{\odot} pc^{-2}$ for M31, M33, and the LMC, respectively—all very close to the densities at which their turnovers happen, being $\approx 4 M_{\odot} pc^{-2}$ for M31, M33, and $\approx 40 M_{\odot} pc^{-2}$ 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.<fnref n="38" />

We examined to what degree $\alpha_{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 $\alpha_{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^{-1} km^{-1} s pc^{-2}$ for M31 and M33, $6.4 K^{-1} km^{-1} s pc^{-2}$ for the LMC, and $21 K^{-1} km^{-1} s pc^{-2}$ for the SMC; see Section 7.1 of Paper I for specifics on our choice of $\alpha_{CO}$ values). The resulting D/H evolutionary profiles are shown in Figure B1.

The plots in Figure B1 show that reducing $\alpha_{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 $\alpha_{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 $\alpha_{CO}$ will increase the D/H measured for a given pixel, while decreasing the $\Sigma_H^{(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 $\alpha_{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 $\alpha_{CO}$: an atomic-gas-dominated pixel with $\Sigma_{H1} = 1.6$ and $\Sigma_{H2} = 0.4$; and a molecular-gas-dominated pixel with $\Sigma_{H1} = 1.3$ and $\Sigma_{H2} = 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 $\alpha_{CO}$ is multiplied by 0.1. Now, the originally atomic-dominated pixel will have $\Sigma_{H1} = 1.6$ and $\Sigma_{H2} = 0.04$, while the formerly molecular-dominated pixel will have $\Sigma_{H1} = 1.3$ and $\Sigma_{H2} = 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 $\alpha_{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 $\Sigma_H^{(deproj)}$ will tend toward a gradient of -1, in a classic instance of noise-induced anticorrelation,<fnref n="39" /> which is indeed what we see happening for M33 in Figure B1.

Given that significantly reducing $\alpha_{CO}$ gives rise to this sort of spurious effect in our data, it would seem that reducing $\alpha_{CO}$ by such larger factors is likely unphysical, which should not be surprising. For instance, reducing the $\alpha_{CO}$ of the LMC by a factor of 5 would make it only $1.28 K^{-1} km^{-1} s pc^{-2}$—less than half the standard Milky Way value of $3.2 K^{-1} km^{-1} s pc^{-2}$, 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 $\alpha_{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 $\alpha_{CO}$.

As such, we are confident that overestimation of $\alpha_{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 $\alpha_{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 $\Sigma_H^{(deproj)} = 1 M_{\odot} pc^{-2}$, to 0.01 at $\Sigma_H^{(deproj)} = 10 M_{\odot} 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).

<footnote n="38">³⁸ While the SMC of course has regions that are molecular-gas dominated, there is no bin of $\Sigma_H^{(deproj)}$ in which the mass contribution of molecular gas exceeds that of atomic gas. This is most likely due to a combination of our limiting resolution, and the extreme line-of-sight depth of the SMC. The contribution of a given molecular-gas-dominated region to the surface density of a given pixel can be exceeded by the contributions of very long columns of atomic gas in front of and behind it, further diluted by the limits of our spatial resolution.</footnote>

<footnote n="39">³⁹ For any plot of $a/b$ versus $b$, where there is no intrinsic correlation between $a$ and $b$, the gradient $b$ will tend toward -1 as noise in $a$ and $b$ increases.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure B1.** Plots of D/H against $\Sigma_{\text{H}}$, for which molecular gas surface densities have been recomputed using values of $\alpha_{\text{CO}}$ that have been changed by factors of 0.5 (left), 0.2 (center), and 0.1 (right), as compared to each galaxy's fiducial value of $\alpha_{\text{CO}}$ used in Figure 9.</figcaption></figure>

First, we assess the impact that generally large scatter could have when observing this trend. We created simulated data for which $10^6$ data points, evenly distributed in $\log_{10}(\Sigma_{\text{H}}^{(\text{deproj}))}$, were drawn from the toy model over a $10^{-2} < \Sigma_{\text{H}}^{(\text{deproj}))} < 10^5$ interval.<fnref n="40" /> We then applied 1 dex of Gaussian noise to the $\Sigma_{\text{H}}^{(\text{deproj}))}$ and $\Sigma_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.<fnref n="41" /> 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 $\Sigma_{\text{H}}^{(\text{deproj}))} = 10 M_{\odot} \text{ 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.<fnref n="42" /> We proceed as above, but now, instead of being a constant 1 dex, the scale of the Gaussian scatter on both $\Sigma_{\text{H}}^{(\text{deproj}))}$ and $\Sigma_d$ is a function of $\Sigma_{\text{H}}^{(\text{deproj}))}$, falling as density increases. Scatter in dex, $\sigma$, is given by $\sigma = -0.5 \log_{10}(\Sigma_{\text{H}}^{(\text{deproj}))} + 1.0$, with a floor value of $\sigma = 0.1$ dex imposed at the highest densities where $\sigma$ 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

$\Sigma_{\text{H}}^{(\text{deproj}))}$ and 0.2 dex in $\Sigma_d$ for values up to $\Sigma_{\text{H}}^{(\text{deproj}))} = 10 M_{\odot} \text{ pc}^{-2}$; above this, the noise grows sharply with $\Sigma_{\text{H}}^{(\text{deproj}))}$ according to a power law, increasing to 0.15 dex in $\Sigma_{\text{H}}^{(\text{deproj}))}$ and 1.2 dex in $\Sigma_d$ at $\Sigma_{\text{H}}^{(\text{deproj}))} = 1000 M_{\odot} \text{ 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 than 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 $\Sigma_{\text{H}}^{(\text{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

<footnote n="40">⁴⁰ This interval is larger than the $10^0 < \Sigma_{\text{H}}^{(\text{deproj}))} < 10^3$ range that we plot and discuss. The reason being that if we only drew values from within the range of interest, then some of the values would be scattered out of interval when we apply noise, biasing the binned averages at either end.</footnote>

<footnote n="41">⁴¹ Among the bins for our actual data for the sample galaxies, 70% of bins' values have standard deviations of <1 dex, and 92% have standard deviations of <1.5 dex (all bins for the LMC and SMC have standard deviations of <1.41 dex).</footnote>

<footnote n="42">⁴² Roughly speaking, typical scatter in bins falls from ≈1.5 dex at $\Sigma_{\text{H}}^{(\text{deproj}))} = 1 M_{\odot} \text{ pc}^{-2}$, down to ≈0.2 dex at $\Sigma_{\text{H}}^{(\text{deproj}))} = 50 M_{\odot} \text{ pc}^{-2}$.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure B2.** Plot of simulated trends in D/H against $\Sigma_{H}$, to examine whether noise can introduce artificial anticorrelation at high surface densities. The gray line shows the underlying truth value for our toy model. The blue squares show the medians of each 0.025 dex bin, when using artificial observed data, where each measurement has 1 dex of Gaussian scatter in both $\Sigma_{H}^{(deproj)}$ and $\Sigma_{d}^{(deproj)}$. The green right-pointing triangles show bin medians where the scatter gets smaller as $\Sigma_{H}^{(deproj)}$ increases. The pink left-pointing triangles show bin medians where the scatter is constant up until $\Sigma_{H}^{(deproj)} = 10 M_{\odot} \text{ pc}^{-2}$, then steadily increases at higher densities.</figcaption></figure>

compare fairly to both.43 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
$\Sigma_{H}^{(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.44 Other-
wise, the degraded data simply traces the same D/H evolution
profile as the original data, albeit with somewhat more scatter.

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 temper-
ature 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

<figure><img src="image_22.png" /><figcaption>**Figure B3.** Plot of D/H against $\Sigma_{H}^{(deproj)}$ in which values for the LMC and SMC have been produced using maps degraded to recreate how the Magellanic Clouds would appear in our analysis were they at the same distance as M31 and M33, thereby giving all four galaxies the same physical resolution. The original data for the LMC and SMC (as plotted in Figure 9) is indicated by the crosses, for comparison. Error bars have been omitted to assist in clarity, but are not significantly different for the degraded data than for the original, as plotted in Figure 9.</figcaption></figure>

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_{\odot} \text{ pc}^{-2}$. 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).

<footnote n="43">43 We felt no need to make any adjustments to the data for M31 and M33, given how similar their distances are already. Indeed, the nearest parts of the disk of M33 are closer to us than the farthest parts of the disk of M31, so degrading the LMC and SMC data to match the midpoint between the two should allow for an entirely fair comparison.</footnote>

<footnote n="44">44 The reason why the LMC and SMC D/H profiles still probe to greater Σ<sub>H</sub><sup>(deproj)</sup> than for M31 and M33, despite their data being degraded to the same effective distance, is because of the larger deprojection corrections applied to M31 and M33, due to their greater inclination. For instance, were the factor 0.22 correction not applied to M31, it would appear to probe up to densities of 80 M<sub>☉</sub> pc<sup>-2</sup>, closely matching the highest densities probed by the LMC and SMC in their degraded data.</footnote>
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 $\Sigma_{H}^{(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 $\Sigma_{H}^{(deproj)}$.

To test this, we use the H$\alpha$ 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$\alpha$ maps to create maps of the surface density of ionized gas, $\Sigma_{H+}^{(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 $\Sigma_{H}^{(deproj)}$ for the LMC and SMC, splitting the data into six bins of ionized gas surface density, from $\Sigma_{H+}^{(deproj)} = 10^{-1.5} M_{\odot} pc^{-2}$ to $\Sigma_{H+}^{(deproj)} = 10^{1.5} M_{\odot} pc^{-2}$, with each bin having a width of 0.5 dex. In Figure B4, the D/H evolution profile within each of these $\Sigma_{H+}^{(deproj)}$ bins is plotted, for the LMC and SMC. Note that for these plots, $\Sigma_{H}^{(deproj)}$, including the denominator of D/H, has been revised to also include the contribution of $\Sigma_{H+}^{(deproj)}$, to ensure internal consistency.45 We also plot the global D/H evolution profile for both galaxies in Figure B4, incorporating the ionized gas into $\Sigma_{H}^{(deproj)}$.

Figure B4 shows that the D/H evolution profile behaves very similarly for regions of high $\Sigma_{H+}^{(deproj)}$ as it does in regions of low $\Sigma_{H+}^{(deproj)}$. The only exception to this is for the regions of the very highest ionized gas surface density, $10^{1.0} < \Sigma_{H+} < 10^{1.5} M_{\odot} pc^{-2}$, 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 $\Sigma_{H+}^{(deproj)}$ that trace that range of $\Sigma_{H}^{(deproj)}$ (although the profile for the highest $\Sigma_{H+}^{(deproj)}$ bin does get noisy at $\Sigma_{H} > 100 M_{\odot} pc^{-2}$).

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 $\Sigma_{H}^{(deproj)}$ is not the cause of the D/H turnover.

## B.5. Dark Gas?

To measure the gas content of our galaxies, we use 21 cm and CO observations as tracers of the atomic and molecular gas components. However, if these tracers miss some fraction of the gas in certain environments, we would find incorrect D/H

<figure><img src="image_23.png" /><figcaption>**Figure B4.** Plot of D/H against $\Sigma_{H}^{(deproj)}$ for the LMC (upper) and SMC (lower), with each line showing the relationship for pixels that lie within a certain range of ionized gas surface densities, $\Sigma_{H+}^{(deproj)}$. Also plotted, in black points, is the relationship when all pixels are counted. Note that for these plots, $\Sigma_{H}^{(deproj)}$ (and therefore D/H) incorporates the contribution of $\Sigma_{H+}^{(deproj)}$ (in contrast to Figure 9, etc.).</figcaption></figure>

values, thereby introducing errors into the D/H evolution profiles.

In particular, the presence of optically thick H I (Fukui et al. 2015; Murray et al. 2018), or CO-dark H$_2$ (Reach et al. 1994; Grenier et al. 2005; Wolfire et al. 2010), could lead us to underestimate the amount of gas present in given environment, and therefore artificially inflate D/H. If this were to preferentially happen over a specific $\Sigma_{H}$ regime, then what ought to be a plateau in D/H could instead incorrectly manifest as a bump. If such a bump happened at the same surface density regime where D/H evolution transitioned from growth to plateau (see Section 4.2), then we would observe a spurious turnover in D/H.

### B.5.1. Dark Gas in the LMC?

First, we consider this possibility for the case of the LMC. The D/H for the LMC has a peak value of 0.0056, at $\Sigma_{H}^{(deproj)} = 40 M_{\odot} pc^{-2}$. Then, as density increases, D/H falls to an average of 0.0024 for $\Sigma_{H}^{(deproj)} \ge 10 M_{\odot} pc^{-2}$. If the true

<footnote n="45">45 We do note include the contribution of $\Sigma_{H+}^{(deproj)}$ when considering $\Sigma_{H}$ elsewhere in this paper, because in general we do not expect the ionized gas to be colocal with the cold dense ISM where grain growth should occur. Therefore including $\Sigma_{H+}^{(deproj)}$ would be expected to weaken our ability to trace evolution in D/H. Moreover, the D/H evolution profile does not change significantly when $\Sigma_{H+}^{(deproj)}$ is included—compare the black points in Figure B4 to the corresponding points in Figure 9.</footnote>
value of the D/H plateau is 0.0024, and the peak to D/H = 0.0056 is being caused artificially by dark gas, this would require 57% of the gas content at densities around $40 M_{\odot} pc^{-2}$ to be dark, not traced by 21 cm or CO emission.

In general, CO-dark molecular gas would be expected to be more prevalent in regions with the lowest molecular gas densities (Glover et al. 2010; Wolfire et al. 2010), while optically thick H I is most likely to occur in regions with higher atomic-gas density (Lee et al. 2015). Therefore the contribution of dark gas is likely to be greatest at intermediate densities, where these two regimes overlap—which should broadly correspond to the density where $\Sigma_{\text{H}}^{(\text{deproj})}$ transitions from being atomic-gas dominated to molecular-gas dominated.

As already discussed in Appendix B.1, this atomic-molecular transition region, at $\Sigma_{\text{H}}^{(\text{deproj})} = 35 M_{\odot} pc^{-2}$ (see also Roman-Duval et al. 2014), well matches the surface density regime where the D/H turnover starts for the LMC. At first glance, this is highly suggestive that dark gas could be contributing to the apparent turnover in LMC.

However, dark gas cannot *depress* D/H measurements, only raise them. So, given that D/H falls to $<3 \times 10^{-3}$ at the highest measured LMC densities of $\Sigma_{\text{H}}^{(\text{deproj})} = 200 M_{\odot} pc^{-2}$, this means that for dark gas to be causing the turnover, all D/H measurements above $3 \times 10^{-3}$ would have to be artificially inflated by dark gas, and in fact should be no more than $3 \times 10^{-3}$. However, the bins with D/H greater than $3 \times 10^{-3}$ span a very wide range of density, $5-200 M_{\odot} pc^{-2}$. If we have to assume D/H is in reality no more than $3 \times 10^{-3}$ over this entire range, then this would require D/H to be essentially flat over 1.6 dex in density—tantamount to saying that there is essentially no D/H evolution with density in the LMC. We rule out this possibility, especially because D/H estimates for the LMC determined via UV absorption line depletion measurements from Roman-Duval et al. (2021) agree excellently with our own for densities of $2 < \Sigma_{\text{H}}^{(\text{deproj})} < 40 M_{\odot} pc^{-2}$; they find a factor $\approx 3$ increase in D/H over this range, with a peak D/H of 0.0046, within 18% our own. This lets us be confident that D/H is not flat over this range, and therefore the D/H turnover is not predominantly a dark gas artifact. We perform a detailed comparison with the Roman-Duval et al. (2021) UV measurements of D/H in Section 6.

Even if the D/H turnover for the LMC cannot be explained by dark gas, it is nonetheless worth considering what other effect dark gas is having upon the LMC D/H evolution profile. It remains highly suggestive that D/H peaks at nearly the same $\Sigma_{\text{H}}^{(\text{deproj})}$ where atomic-to-molecular transition occurs in the LMC. Moreover, D/H shows a conspicuous “bump” at densities between 30 and $60 M_{\odot} pc^{-2}$. This bump occurs at a very similar range of surface densities as those where the simulations of Glover & Mac Low (2011) find CO-dark gas to be significant. At average visual extinction of $A_V < 1$ mag, Glover & Mac Low (2011) found that $H_2$ makes up a negligible fraction of the total H mass, while at $A_V > 3$ mag, reduced photodissociation means that effectively all gas-phase carbon is found in CO, causing the CO-to-$H_2$ conversion stabilize. Given the Glover & Mac Low (2011) conversion between column density and $A_V$, being $5.348 \times 10^{-22} Z/Z_{\odot}$ mag cm$^2$, and given that $1 M_{\odot} pc^{-2} = 1.25 \times 10^{20} cm^{-2}$, the $1 < A_V < 3$ mag regime where Glover & Mac Low (2011) found CO-dark $H_2$ to be significant corresponds to a surface density range of

$30-90 M_{\odot} pc^{-2}$—a close match to the D/H bump in the LMC. If this is indeed the origin of this bump feature, it suggests that approximately 30% of the gas in the $30-60 M_{\odot} pc^{-2}$ density range is CO-dark.

So far we have not considered dark, optically thick H I here. However, the contribution of dark H I would be expected to occur *below* the densities at which $H_2$ starts to dominate (and be traced by CO). Therefore, dark H I is not a likely explanation for the D/H turnover in the LMC, and it occurs at densities above those at which molecular gas begins to dominate.

So, while dark gas cannot be the driver of the high-density D/H turnover in the LMC, it can make it appear even more conspicuous, by elevating D/H at these intermediate densities.

### B.5.2. Dark Gas in M31 and M33?

Next we consider whether dark gas could be causing a spurious D/H turnover in the case of M31 and M33. Because of the striking similarity in how D/H evolves with $\Sigma_{\text{H}}^{(\text{deproj})}$ for both galaxies, including the turnover, we focus this analysis on M31, as it has more data, over a wider range of densities—under the expectation that explanations for our observations of M31 will also be valid for M33.

The peak D/H for M31 is 0.01, at $\Sigma_{\text{H}}^{(\text{deproj})} = 4 M_{\odot} pc^{-2}$. Then D/H falls to an average of 0.0067 for $\Sigma_{\text{H}}^{(\text{deproj})} > 10 M_{\odot} pc^{-2}$. If the true value of the D/H plateau is 0.0067, and the peak to D/H = 0.01 is being caused artificially by dark gas, this would require 33% of the gas content at densities around $4 M_{\odot} pc^{-2}$ to be dark, not traced by 21 cm or CO emission. This is fairly plausible. A dark gas fraction of 33% falls well within the range of estimates proposed by various authors (Grenier et al. 2005; Braun et al. 2009; Abdo et al. 2010; Planck Collaboration et al. 2011d; Paradis et al. 2012). Additionally, the turnover for M31 happens at the density where the gas is transitioning from atomic- to molecular-gas dominated, the regime where dark gas is most likely to have an effect.<fnref n="46" />

On the other hand, studies focused on M31 have suggested that it hosts a very minimal CO-dark molecular gas component (Smith et al. 2012; Athikkat-Eknath et al. 2021). Additionally, because dark gas is expected to be more common in lower-metallicity galaxies (Genzel et al. 2012; Madden et al. 2020), any dark gas artifact should be more prominent for the LMC than for M31—whereas in Appendix B.5.1 above, we have just established that the influence of dark gas in the LMC appears most likely limited to a range of intermediate densities, and not driving the D/H turnover at higher densities.

As such, we are unable to say with confidence whether dark gas is a major driver of the turnover in M31 and M33, but it seems a plausible explanation.

### B.5.3. Dark Gas in the SMC?

Although D/H does not turn over for the SMC, we nonetheless briefly consider what effect dark gas could be having on our measurements for this system, too.

In particular, we note that the relationship between D/H and density in the SMC gets abruptly steeper at

<footnote n="46">⁴⁶ Unfortunately, our poorer spatial resolution for M31 (and M33) means we cannot distinguish the highest-density regimes at $\Sigma_{\text{H}}^{(\text{deproj})} > 15 M_{\odot} pc^{-2}$, preventing us from performing the detailed $\Sigma_{\text{H}}^{(\text{deproj})}$-to-A$_V$ comparison to the Glover & Mac Low (2011) CO-dark $H_2$ models that we could test for the LMC.</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure B5.** Plots of D/H against $\Sigma_H$, for M31 (left), M33 (center), and the LMC (right), where we examine the effect of decreasing the dust mass absorption coefficient, $\kappa$, above a certain break density, according to a power-law slope of -0.4, as per Clark et al. (2019).</figcaption></figure>

$\Sigma_H^{(deproj)} \approx 50 M_{\odot} \text{ pc}^{-2}$. Once again considering the Glover & Mac Low (2011) simulations of CO-dark gas in different environments, we note that $\Sigma_H^{(deproj)} = 50 M_{\odot} \text{ pc}^{-2}$ corresponds to $A_V = 0.85$ mag, once again very close to the point at which they predict CO-dark $H_2$ to become significant. This also matches closely with the density at which Roman-Duval et al. (2014) found the atomic-to-molecular transition begins, at $80 M_{\odot} \text{ pc}^{-2}$.

If the steepening in D/H at $>50 M_{\odot} \text{ pc}^{-2}$ in the SMC is indeed an artifact of CO-dark gas, then it is conceivable that the true evolution of D/H at these higher densities in fact follows the shallower trend found at $<50 M_{\odot} \text{ pc}^{-2}$. If this is the case, then that suggests that >60% of the gas at $\Sigma_H^{(deproj)} > 100 M_{\odot} \text{ pc}^{-2}$ is dark. If the Glover & Mac Low (2011) simulations are correct, then we would expect the $A_V > 3$ mag regime (where CO reliably traces $H_2$) to start at $150-200 M_{\odot} \text{ pc}^{-2}$, above the densities we are able to probe with our spatial resolution. Nonetheless, the effect of dark gas could even be masking the beginning of a D/H plateau at the highest densities we can sample.

## 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, $\kappa$, varies as a function of $\Sigma_H^{(deproj)}$.

It is essentially guaranteed that $\kappa$ 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 $\kappa$ would not evolve with density to some degree.

For our SED fitting, we assumed a constant value of $\kappa$. If, however, $\kappa$ were to decrease at higher $\Sigma_H^{(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 $\kappa$ decreases in higher-density ISM? The exact nature of how $\kappa$ evolves with density is mostly

unconstrained at present. However, models consistently predict that $\kappa$ 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 $\kappa$ within nearby spiral galaxies M74 and M83, and found evidence that $\kappa$ *decreases* with increasing $\Sigma_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 $\kappa$ falls with $\Sigma_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 $\Sigma_H$ to which Clark et al. (2019) could probe in M74 and M83 was $3 M_{\odot} \text{ pc}^{-2}$ (with the highest well-sampled $\Sigma_H$ being 303 $M_{\odot} \text{ pc}^{-2}$). 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 $\Sigma_H$ range for which Clark et al. (2019) reported a fall in $\kappa$ with $\Sigma_H$ corresponds to the surface densities above which we start to see the D/H turnovers in our target galaxies.

We examined whether a $\kappa$ that falls with increasing density, as suggested by Clark et al. (2019), could explain our apparent D/H turnover. Specifically, we tried a broken $\kappa$ model, where $\kappa$ remains constant at $\kappa_{160} = 1.24 \text{ m}^2 \text{ kg}^{-1}$ (see Section 3) up until a break density, after which $\kappa$ falls with $\Sigma_H$ according to a -0.4 power law, as per Clark et al. (2019). For the location of the $\kappa$ break in each galaxy, we used the surface density at which D/H peaks: $\Sigma_H = 4 M_{\odot} \text{ pc}^{-2}$ for M31 and M33, and $\Sigma_H = 40 M_{\odot} \text{ pc}^{-2}$ for the LMC. The results of this are plotted in Figure B5.

In Figure B5, we see that introducing this broken $\kappa$ 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 $\kappa$ at higher $\Sigma_H$ could be the cause of the D/H turnover. It is
particularly interesting that the −0.4 power-law rate of decrease of $\kappa$ versus $\Sigma_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 $\kappa$ 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.<fnref n="47" />

Why the $\Sigma_H$ at which the D/H turnover happens—and hence the $\Sigma_H$ at which $\kappa$ 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 $\Sigma_H = 4 M_\odot \text{ pc}^{-2}$, 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 $\Sigma_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 $\kappa$ in the SMC, relative to the LMC. While separate from any variation in $\kappa$ with density, this nonetheless further highlights the fact that we really should not expect $\kappa$ to remain constant throughout.

## ORCID iDs

Christopher J. R. Clark † https://orcid.org/0000-0001-7959-4902

Julia C. Roman-Duval † https://orcid.org/0000-0001-6326-7069

Karl D. Gordon † https://orcid.org/0000-0001-5340-6774

Caroline Bot † https://orcid.org/0000-0001-6118-2985

Matthew W. L. Smith † https://orcid.org/0000-0002-3532-6970

Lea M. Z. Hagen † https://orcid.org/0000-0001-8918-1597

## References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 710, 133

Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218

Aoyama, S., Hirashita, H., & Nagamine, K. 2020, MNRAS, 491, 3844

Asano, R. S., Takeuchi, T. T., Hirashita, H., & Inoue, A. K. 2013, EP&S, 65, 213

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Athikkat-Eknath, G., Eales, S. A., Smith, M. W. L., et al. 2021, MNRAS, 511, 5287

Balog, Z., Müller, T., Nielbock, M., et al. 2014, ExA, 37, 129

Balser, D. S. 2006, AJ, 132, 2326

Barlow, M. J., Krause, O., Swinyard, B. M., et al. 2010, A&A, 518, L138

Bekki, K. 2008, MNRAS, 390, L24

Bekki, K., & Chiba, M. 2007, MNRAS, 381, L16

Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, MNRAS, 433, 3062

Bianchi, S., Casasola, V., Baes, M., et al. 2019, A&A, 631, A102

Bianchi, S., Casasola, V., Corbelli, E., et al. 2022, A&A, 664, A187

Bigiel, F., Bolatto, A. D., Leroy, A. K., et al. 2010, ApJ, 725, 1159

Bocchio, M., Jones, A. P., & Slavin, J. D. 2014, A&A, 570, A32

Boggess, N. W., Mather, J. C., Weiss, R., et al. 1992, ApJ, 397, 420

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bot, C., Ysard, N., Paradis, D., et al. 2010, A&A, 523, A20

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: v1.0.1, Zenodo, doi:10.5281/zenodo.596036

Braun, R., Thilker, D. A., Walterbos, R. A. M., & Corbelli, E. 2009, ApJ, 695, 937

Breeveld, A. A., Landsman, W., Holland, S. T., et al. 2011, AIP Conf. Proc. 1358, Gamma Ray Bursts 2010, ed. J. E. McEnery, J. L. Racusin, & N. Gehrels, (Melville, NY: AIP), 373

Buat, V., Giovannoli, E., Takeuchi, T. T., et al. 2011, A&A, 529, A22

Calzetti, D., Kennicutt, R. C., Jr., Bianchi, L., et al. 2005, ApJ, 633, 871

Chastenet, J., Bot, C., Gordon, K. D., et al. 2017, A&A, 601, A55

Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2021, ApJ, 912, 103

Chevance, M., Kruijssen, J. M. D., Hygate, A. P. S., et al. 2020, MNRAS, 493, 2872

Chiang, I. D., Sandstrom, K. M., Chastenet, J., et al. 2018, ApJ, 865, 117

Choi, Y., Olsen, K. A. G., Besla, G., et al. 2022, ApJ, 927, 153

Cioni, M. R. L. 2009, A&A, 506, 1137

Clark, C. J. R. 2021, Herschel Maps of Local Group Galaxies v1.0, IPAC doi:10.26131/IRSA545

Clark, C. J. R., De Vis, P., Baes, M., et al. 2019, MNRAS, 489, 5256

Clark, C. J. R., Dunne, L., Gomez, H. L., et al. 2015, MNRAS, 452, 397

Clark, C. J. R., Roman-Duval, J. C., Gordon, K. D., Bot, C., & Smith, M. W. L. 2021, ApJ, 921, 35

Clark, C. J. R., Verstocken, S., Bianchi, S., et al. 2018, A&A, 609, A37

Clark, P. C., & Glover, S. C. O. 2015, MNRAS, 452, 2057

Corbelli, E., Thilker, D., Zibetti, S., Giovanardi, C., & Salucci, P. 2014, A&A, 572, A23

Cuppen, H. M., & Herbst, E. 2007, ApJ, 668, 294

Davies, J. I., Baes, M., Bianchi, S., et al. 2017, PASP, 129, 044102

Dawson, J. R., McClure-Griffiths, N. M., Wong, T., et al. 2013, ApJ, 763, 56

De Cia, A., Ledoux, C., Mattsson, L., et al. 2016, A&A, 596, A97

De Looze, I., Barlow, M. J., Bandiera, R., et al. 2019, MNRAS, 488, 164

De Looze, I., Lamperti, I., Saintonge, A., et al. 2020, MNRAS, 496, 3668

De Vis, P., Dunne, L., Maddox, S., et al. 2017a, MNRAS, 464, 4680

De Vis, P., Gomez, H. L., Schofield, S. P., et al. 2017b, MNRAS, 471, 1743

De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5

Decin, L., & Eriksson, K. 2007, A&A, 472, 1041

Dopcke, G., Glover, S. C. O., Clark, P. C., & Klessen, R. S. 2011, ApJL, 729, L3

Draine, B. T. 2011, ApJ, 732, 100

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118

Dwek, E., Galliano, F., & Jones, A. P. 2007, ApJ, 662, 927

Eales, S. A., Smith, M. W. L., Wilson, C. D., et al. 2010, A&A, 518, L62

Edmunds, M. G., & Eales, S. A. 1998, MNRAS, 299, L29

Elmegreen, B. G., Kim, S., & Staveley-Smith, L. 2001, ApJ, 548, 749

Feldmann, R. 2015, MNRAS, 449, 3274

Fitzpatrick, E. L., & Massa, D. 2005, AJ, 130, 1127

Foreman-Mackey, D. 2016, JOSS, 1, 24

Forgan, D. H., Rowlands, K., Gomez, H. L., et al. 2017, MNRAS, 472, 2289

Fukui, Y., Torii, K., Onishi, T., et al. 2015, ApJ, 798, 6

Fukui, Y., Tsuge, K., Sano, H., et al. 2017, PASJ, 69, L5

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

Galliano, F., Hony, S., Bernard, J. P., et al. 2011, A&A, 536, A88

Galliano, F., Madden, S. C., Jones, A. P., et al. 2003, A&A, 407, 159

Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18

Gardan, E., Braine, J., Schuster, K. F., Brouillet, N., & Sievers, A. 2007, A&A, 473, 91

Gaustad, J. E., McCullough, P. R., Rosing, W., & Van Buren, D. 2001, PASP, 113, 1326

Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005

Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69

Glover, S. C. O., Federrath, C., Mac Low, M. M., & Klessen, R. S. 2010, MNRAS, 404, 2

Glover, S. C. O., & Mac Low, M. M. 2011, MNRAS, 412, 337

Gomez, H. L., Krause, O., Barlow, M. J., et al. 2012, ApJ, 760, 96

Gordon, K. D., Bailin, J., Engelbracht, C. W., et al. 2006, ApJL, 638, L87

<footnote n="47">⁴⁷ Given that the lead author of this work was also lead author of the Clark et al. (2019) study, they feel it worth mentioning that they did not enter into this current study with any intention of searching for anticorrelation between $\kappa$ and $\Sigma_H$, nor did they have any expectation that they would find evidence for it. Rather, the fact that the Clark et al. (2019) result resolves the D/H turnover found here came as a surprise.</footnote>