A Dusty Locale: Evolution of Galactic Dust Populations from Milky Way to Dwarf-Mass Galaxies

Observations indicate dust populations vary between galaxies and within them, suggesting a complex life cycle and evolutionary history. Here we investigate the evolution of galactic dust populations across cosmic time using a suite of cosmological zoom-in simulations from the Feedback in Realistic Environments (FIRE) project, spanning $M_{\rm vir}=10^{9-12}M_{\odot};\,M_{*}=10^{6-11}\,M_{\odot}$. Our simulations incorporate a dust evolution model that accounts for the dominant sources of dust production, growth, and destruction and follows the evolution of specific dust species. All galactic dust populations in our suite exhibit similar evolutionary histories, with gas-dust accretion being the dominant producer of dust mass for all but the most metal-poor galaxies. Similar to previous works, we find the onset of efficient gas-dust accretion occurs above a `critical' metallicity threshold ($Z_{\rm crit}$). Due to this threshold, our simulations reproduce observed trends between galactic D/Z and metallicity and element depletion trends in the ISM. However, we find $Z_{\rm crit}$ varies between dust species due to differences in key element abundances, dust physical properties, and life cycle processes resulting in $Z_{\rm crit}\sim0.05Z_{\odot},\,0.2Z_{\odot},\,0.5Z_{\odot}$ for metallic iron, silicates, and carbonaceous dust, respectively. These variations could explain the lack of small carbonaceous grains observed in the Magellanic Clouds. We also find a delay between the onset of gas-dust accretion and when a dust population reaches equilibrium, which we call the equilibrium timescale ($\tau_{\rm eq}$). The relation between $\tau_{\rm eq}$ and the metal enrichment timescale of a galaxy, determined by its recent evolutionary history, can contribute to the scatter in the observed relation between galactic D/Z and metallicity.


INTRODUCTION
Observations of the Milky Way (MW) and the Large and Small Magellanic Clouds (LMC & SMC) reveal significant variations in their respective dust populations.The 2175 Å feature/bump in dust extinction curves, which correlates with the abundance of small carbonaceous grains 1 , shows a strong decrease between the Milky Way and LMC, and an almost complete lack of the feature for the SMC (e.g.Pei 1992;Weingartner & Draine 2001).Mid-infrared dust emis-★ email: cchoban@iu.edu 1 While it is currently unknown whether PAHs or small carbonaceous grains (i.e.amorphous graphite) are the dominant carrier for the 2175 Å feature, their life cycles could be intimately linked.One proposed formation mechanism for PAHs is through photoprocessing, where hydrogen atoms are removed from small carbonaceous dust grains (i.e.aromatization; Rau et al. 2019;Hirashita & Murga 2020).However, other possible formation mechanisms via AGB winds (Galliano et al. 2008) or directly through low-temperature chemical channels in dense molecular clouds (Parker et al. 2012) have not been ruled out.
sion produced by polycyclic aromatic hydrocarbons (PAHs) shows a dramatic decrease in the average fraction of dust mass comprised of PAHs ( PAH ;Draine & Li 2007) from the Milky Way (∼4.6%) to the LMC (∼3.3%) and SMC (∼1.0%) (Li & Draine 2001;Weingartner & Draine 2001;Chastenet et al. 2019).Gas-phase element depletions in the MW, LMC, and SMC exhibit an increasing fraction of metals locked in dust (dust-to-metals ratio; D/Z) with gas surface density but an overall decrease in galactic D/Z corresponding with their relative metallicities (Roman-Duval et al. 2022a).These observations also reveal a considerable variation in dust chemical composition; however, the lack of observed C and O depletions in the LMC and SMC results in an incomplete picture.
Looking further afield, observations of local galaxies tell a similar story.Dust emission surveys find an overall increase of galaxyintegrated D/Z with metallicity along with a large (>1 dex) scatter at any given metallicity, but the exact relation varies between studies (Rémy-Ruyer et al. 2014;De Vis et al. 2019).These studies also find a dependence of  PAH with logarithmic metallicity, with a threshold galactic metallicity of 12 + log 10 (O/H) ∼ 8.0 above which  PAH rapidly increases (Draine et al. 2007;Rémy-Ruyer et al. 2015;Aniano et al. 2020).Local and high-z dust attenuation curves in galaxies exhibit a wide range of 2175 Å feature strengths (e.g.Salim & Narayanan 2020;Shivaei et al. 2022).However, the strength/absence of this feature is not solely determined by the underlying dust population-derived extinction curve.In particular, radiative transfer effects caused by interstellar medium (ISM) clumpiness, star-dust geometry, and adopted dust albedo can modify the expected attenuation feature (Granato et al. 2000;Panuzzo et al. 2007;Seon & Draine 2016;Narayanan et al. 2018).Observations of damped Ly systems (DLAs) show increasing gas-phase element depletions with metallicity and no clear relation with redshift (Péroux & Howk 2020;Roman-Duval et al. 2022b).These observations suggest an evolving dust population in both amount and chemical composition, with a decreasing total fraction of metals locked in dust and a scarcity of carbonaceous dust in low-metallicity environments.The exact cause of this is unknown, but detailed dust evolution models that account for the main processes of the dust life cycle coupled with galaxy formation models can help elucidate this issue.
The current dust life cycle paradigm posits all dust grains begin their life in the ejecta of Type II supernovae (SNe II) and stellar winds of asymptotic giant branch (AGB) stars, where a fraction of metals within these winds coalesce into dust (Draine 1990;Todini & Ferrara 2001;Dwek 2005;Ferrarotti & Gail 2006).The grains are then primarily destroyed by SNe shocks.In isolation, stellar dust production alone cannot explain the dust content of either the Milky Way (e.g.McKee 1989;Draine 2009) or LMC (Zhukovska & Henning 2013), requiring another growth/production mechanism.Observations of gas-phase element depletions (Jenkins 2009;Roman-Duval et al. 2021) and dust emission (Roman-Duval et al. 2014, 2017;Chiang et al. 2018;Clark et al. 2023) show D/Z increases from diffuse to dense regions of the Milky Way and Local Group galaxies.These provide strong evidence that dust grows via the accretion of gasphase metals in dense environments and could be the main producer of dust mass within these galaxies.However, the exact details of the accretion process are poorly understood, and many questions have yet to be completely answered.When does gas-dust accretion become efficient, and does it differ between dust species?How long does it take for accretion to build up a sizable dust mass within a galaxy?Can accretion explain dusty galaxies at high redshift?
Analytical galactic dust evolution models were first developed to answer these questions, giving rise to the concept of a 'critical' metallicity ( crit ) threshold (Inoue 2011;Asano et al. 2013;Zhukovska 2014;Feldmann 2015;Popping et al. 2017;Triani et al. 2020).While the exact definition varies between works, they all agree there is a 'critical' metallicity above which gas-dust accretion becomes the dominant source of dust mass, rapidly increasing the expected galactic D/Z 2 .Inoue (2011) found that  crit is determined by the competition between growth via accretion and destruction via SNe shocks.However, equilibrium models developed by Feldmann (2015) suggest the dilution of dust via dust-poor gas inflows dominates over destruction via SNe shocks when determining  crit , especially for dwarf galaxies.Zhukovska et al. (2008) developed analytical models that track the evolution of specific dust species and found that for a Milky Way-mass galaxy, silicates, carbonaceous, and metallic iron dust species have different  crit .In later works, this model was used 2 These findings are not unanimous among all works.In particular, Priestley et al. (2022) suggests the contribution of dust growth via accretion may be overestimated if both high stardust creation efficiencies and increased SNe dust destruction in low-metallicity environments are assumed.
to investigate the evolution of dust within late-type dwarf galaxies experiencing episodic starbursts (Zhukovska 2014).They concluded that such galaxies have lower  crit than galaxies with constant star formation due to long quiescent periods between starbursts where dust has ample time to grow.
More recently, numerical dust evolution models integrated into galaxy formation and evolution simulations have been put to the task (e.g.Bekki 2015;McKinnon et al. 2016;Zhukovska et al. 2016;McKinnon et al. 2017;Aoyama et al. 2020;Granato et al. 2021).The dust evolution models utilized in these works differ in their methodology and included physics and the galaxy simulations employed primarily encompass cosmological volumes, relying on varying sub-resolution star formation, feedback, and ISM prescriptions.Despite these differences, they all agree with the 'critical' metallicity concept (albeit with a large range  crit ∼ 0.03 − 0.2 Z ⊙ ) above which the average galactic D/Z rapidly increases due to dust growth via accretion (Hou et al. 2019;Li et al. 2019;Graziani et al. 2020;Parente et al. 2022).A majority of these works also find that  crit and the resulting relation between D/Z and Z have little evolution over redshift3 .However, the exact cause of the large scatter in observed D/Z at any given Z is debated.Parente et al. (2022) predicts this is due to the fraction of cold gas in each galaxy, with higher cold gas fractions allowing for the faster build-up of dust via accretion.Li et al. (2019) proposes a more complex dependence on galactic metal enrichment history reflected by  and  * and evolutionary stages quantified by gas content and depletion timescales.Li et al. (2021) adds that variance in galactic mass-averaged grain sizes (and thus gas-dust accretion rates) can also drive this scatter.Regarding the differences in the evolution between silicate and carbonaceous dust species, Granato et al. (2021) and Parente et al. (2022) predict that silicate growth lags behind carbonaceous growth due to the overall lower abundance of Si compared to C, with silicate dust only dominating the dust mass past  crit ≳ 0.15Z ⊙ ,  * ≳ 10 8.5  ⊙ .
While these dust evolution models generally agree in terms of galaxy-integrated dust observations, they have major limitations that make them ill-suited tools for investigating the gas-dust accretion process.Models integrated into cosmological galaxy simulations do not resolve the multi-phase ISM, which is necessary to accurately model each process in the dust life cycle.Furthermore, most models do not track the evolution of chemically distinct dust species, instead assuming a single, chemically ambiguous dust population, affecting how gas-dust accretion is modeled.Choban et al. (2022, C22 hereafter) developed two separate dust evolution models to investigate the effects of both these assumptions.They found that only a model that tracks the evolution of chemically-distinct dust species and incorporates a physically motivated, two-phase dust growth routine can reproduce the observed scaling relation between individual element depletions and D/Z with column density and local gas density in the Milky Way.This relation results from the equilibrium between SNe dust destruction and dust growth via accretion.Therefore, such a model is crucial for investigating the evolution of galactic dust populations and testing the importance of gas-dust accretion in that evolution.However, C22 only utilized simulations of an idealized, non-cosmological Milky Way-like galaxy.Due to their idealized nature, these simulations do not capture the formation and hierarchical growth of galaxies in fully cosmological settings across cosmic time and do not include a realistic baryon cycle in the IGM/CGM/ISM.Consequently, these predictions may not hold for the entire evolutionary history of a galaxy or lower-mass galaxies.
In this work, we present a subset of cosmological zoom-in simulations of Milky Way to dwarf-halo mass galaxies from the Feedback in Realistic Environments (FIRE) project4 rerun with the integrated "Species" dust evolution model presented in C22.This model tracks the evolution of specific dust species with set chemical compositions and incorporates a physically motivated dust growth routine.With these simulations, we investigate how galactic dust populations evolve in both their amount and chemical composition, focusing on the determining factors for when gas-dust accretion dominates and how it differs between dust species.We find that gas-dust accretion is the dominant producer of dust mass for all but the most metal-poor galaxies and, in the case of the MW, dominates for the majority of the galaxy's life.We discover that the onset of rapid growth via gasdust accretion differs between dust species, arising from differences in their key element abundances, physical properties, and life cycle processes.These differences can explain the variable dust population, in both amount and composition, in the MW, LMC, and SMC.We also find a delay between the onset of rapid dust growth via accretion and when a dust population reaches equilibrium between growth and destruction processes.The relation between this delay and the metal enrichment timescale of a galaxy can contribute to the scatter in observed D/Z at any given metallicity.
This paper is organized as follows.In Section 2, we provide a brief overview of our simulation sample along with the galaxy formation and dust evolution model used.In Section 3, we present the results of our simulations, focusing on the evolution of the galactic and dust population properties for each galaxy in Section 3.1 and comparing them with local observations in Section 3.2 & 3.3.We discuss the concept of a 'critical' metallicity threshold, which marks the onset of efficient gas-dust accretion, in Section 4.1 and an 'equilibrium' timescale, which is the time over which a dust population's mass builds up via accretion to an effective equilibrium, in Section 4.2.Finally, we present our conclusions in Section 5.

METHODOLOGY
To study the evolution of galactic dust populations within Milky Way to dwarf-mass galaxies, we reran a subset of cosmological simulations from the FIRE suite presented in Hopkins et al. (2018b) and El-Badry et al. (2018).We selected 7 galaxies with a broad range of stellar (106  ⊙ <  * < 10 11  ⊙ ) and halo (10 9  ⊙ <  vir < 10 12  ⊙ ) masses at present day.We give the exact details of the  = 0 halo virial mass, virial radius, stellar mass, stellar half-mass radius, and mass resolution for each simulation in Table 1.
All simulations in this work are run with the GIZMO code base (Hopkins 2015) in the meshless finite-mass (MFM) mode with FIRE-2 (Hopkins et al. 2018b) model of star formation and stellar feedback (see Appendix A for a subset of simulations run with FIRE-3;Hopkins et al. 2023).FIRE-2, an updated version of FIRE (Hopkins et al. 2014), incorporates multiple sources of stellar feedback, specifically stellar winds (O/B and AGB), ionizing photons, radiation pressure, and supernovae (both Types Ia and II).Gas cooling is followed selfconsistently for T = 10 -10 10 K including free-free, Compton, metalline, molecular, fine-structure, and dust collisional processes while gas is also heated by cosmic rays, photo-electric, and photoionization heating by both local sources and a uniform but redshift dependent meta-galactic background (Faucher-Giguère et al. 2009), including the effect of self-shielding5 .Star formation is only allowed in cold, molecular, and locally self-gravitating regions with number densities above  H = 1000 cm −3 .
Each star particle represents a stellar population with a known mass, age, and metallicity assuming a Kroupa (2002) initial mass function (IMF) from 0.1 − 100 M ⊙ .The luminosity, mass loss rates, and SNe II rates of each star particle are calculated based on the STARBURST99 (Leitherer et al. 1999) libraries, and SNe Ia rates following Mannucci et al. (2006).Metal yields from SNe II, Ia, and AGB winds are taken from Nomoto et al. (2006), Iwamoto et al. (1999), and Izzard et al. (2004) respectively.Evolution of eleven species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) is tracked for each gas cell.Sub-resolution turbulent metal diffusion is modeled as described in Su et al. (2017) and Escala et al. (2018).FIRE-2 adopts the older Anders & Grevesse (1989) solar metal abundances with  ∼ 0.02 so whenever we mention solar abundances we are referring to the Andres & Gravesse abundances. 6 FIRE is ideally suited to investigate galactic dust evolution over cosmic time given its success in matching a wide range of observations related to galaxies, including the mass-metallicity relation and its evolution over redshift (Ma et al. 2016;Marszewski et al. in prep.) and the Kennicutt-Schmidt star formation law (Hopkins et al. 2014;Orr et al. 2018;Gurvich et al. 2020).This success is owed to the high resolution, star formation criteria, cooling to low temperatures, and multi-channel stellar feedback of FIRE, all of which result in a reasonable ISM phase structure and giant molecular cloud (GMC) mass function (Benincasa et al. 2020).These also lead to the selfconsistent development of galactic winds that eject large amounts of gas (Muratov et al. 2015;Anglés-Alcázar et al. 2017) and metals (Muratov et al. 2017;Hafen et al. 2019;Pandya et al. 2021) out of galaxies, preventing excessive star formation and leading to a plausible stellar-halo mass relation.
Our simulations utilize the integrated "Species" dust evolution model presented in C22, which we refer to the reader for full details.This model includes the dominant sources of dust production, tracking and differentiating between dust created from SNe Ia and II, AGB stars, and dust growth from gas-phase metal accretion in the ISM.It includes the dominant dust destruction mechanisms, accounting for dust destroyed by SNe shocks, thermal sputtering, and astration (dust destroyed during the formation of stars).These processes are modeled self-consistently, owing to the FIRE model's indepth treatment of the multi-phase ISM and capacity to time-resolve individual SNe events (Hopkins et al. 2018a).Notably, we restrict gas-dust accretion to cool ( ≤ 300 K) gas and destroy dust locally around individual SNe events, allowing us to track the local variability of dust in the ISM.We also follow the evolution of specific dust species (carbonaceous, silicates, and silicon carbide) and theoretical oxygen-bearing (O-reservoir) and nanoparticle metallic iron (Nano- iron) dust species with set chemical compositions.Consequently, this means each dust species has a key element7 that limits their individual accretion growth rates and the maximum formable amount of said dust species.We also incorporate sub-resolution turbulent dust diffusion, which follows the metal diffusion prescription in FIRE, and a dense molecular gas scheme.This scheme is critical to account for Coulomb enhancement of gas-dust accretion in atomic/diffuse molecular gas and the reduction in carbonaceous dust accretion due to the lock-up of gas-phase C into CO in dense molecular gas.

RESULTS
We first showcase mock HST ugr composite images of each galaxy to highlight the breadth of galaxy types and their varying dust structure contained in our simulation suite.Fig. 1 shows face-on images of all galaxies at  = 0.These images use STARBURST99 (Leitherer et al. 1999) to compute the stellar spectra for each star particle given their age and metallicity.These are then ray-traced through the ISM using the tracked D/Z for each gas cell produced by our dust model along with assumed MW, LMC, or SMC dust opacities from Pei (1992) depending on the median gas cell silicate-to-carbonaceous dust mass (Sil-to-C) ratio 8 in each galaxy (MW∼3, LMC∼5., SMC∼10).We then volume-render the observed images in each band and construct a ugr composite image as seen by a distant observer.m12i, m11d, m11v_halo2, m11v_halo0, and m11e show prominent and detailed dust structure to varying degrees, with dark patches produced by the attenuation of UV/optical light from dust tracing the dense gas in each galaxy.m11i and m10q have no distinguishable features produced by dust due to extremely low D/Z for m11i and the complete lack of both dust and gas for m10q.
We present the evolution of each galaxy in terms of galactic properties and their dust population, in both composition and amount, in Sec.3.1.We then compare our simulations to present-day observations, focusing on gas-phase element depletions in Sec.3.2 and galaxy-integrated and spatially-resolved IR dust emission in Sec.3.3.

Galactic Properties & Dust Population
In Fig. 2, we show the evolution of various properties for each galaxy, specifically total gas mass, stellar mass, star formation rate (SFR), median gas cell 9 metallicity, median D/Z, and median Sil-to-C ratio.These values are determined from star particles and gas cells within 10 kpc of the galactic center for all galaxies 10 .We provide a brief description of each galaxy's evolutionary history, noting any merger events, in Table 1.
In Fig. 3 & 4, we show a detailed breakdown of each galaxy's metal and dust population evolution, including the median metallicity, median D/Z, dust creation source contribution, and dust species composition for all gas within 10 kpc of the galactic center.We also include the breakdown for each gas phase: cold neutral gas ( < 10 3 K), warm neutral gas (10 3 K ≤  < 10 4 K), and warm/hot ionized gas ( ≥ 10 4 K).While the evolution of each galaxy varies considerably, their dust populations follow similar evolutionary trends, which we describe below.
1. SNe-Dominated: Initially, as the first massive stars die, SNe II are the dominant producer of dust mass.This 'stardust' population is dominated by carbonaceous dust, arising from the high SNe II carbonaceous dust production efficiency adopted in our model ( SNII carb = 15%).This produces an extremely low median D/Z ∼ 0.01.

AGB-Dominated:
As stellar populations age and low-mass stars transition to the AGB phase, dust production by AGB stars becomes a sizable component of the dust mass.However, this only dominates the dust population for the lowest-mass dwarf galaxies, which have little or extremely late dust growth via accretion.For our model, lowmetallicity AGB stars produce primarily carbonaceous dust and only x2-3 more dust than SNe II over the lifetime of a stellar population (see Fig. 2 in C22).For these reasons, the median D/Z and dust where  is the number abundance of the element and  is the number of atoms of the element in one formula unit of the dust species under consideration. 8We define this ratio as Sil-to-C = (  sil + iron + O−res )/ carb in our simulations. 9Whenever we mention a median value for gas cells, we specifically mean the mass-weighted median.While these values can differ due to deviations in individual gas parcel mass from the stated resolutions in Table 1, caused by SNe ejecta and stellar wind mass injection, these deviations are minor. 10We avoid using an evolving radius, such as  vir , as our cutoff when determining galactic properties in order to avoid jumps due to merger events.A radius of 10 kpc encompasses the ISM of each galaxy while not including a sizable fraction of the galactic halo, which is sufficient for our focus on interstellar dust.Face-on images of our spiral and dwarf galaxies at  = 0 (note the differences in scale).Each image is a mock Hubble Space Telescope ugr composite accounting for dust attenuation utilizing the D/Z produced by our model along with assumed MW, LMC, or SMC dust opacities depending on the median gas cell Sil-to-C ratio.m12i: A compact, Milky Way-mass spiral galaxy that experienced multiple minor mergers during its life.m11d: An LMC-mass galaxy with bursty star formation and multiple mergers throughout its life.m11v_halo2: An LMC-mass galaxy with faint spiral arms that evolved in isolation for a majority of its life.m11v_halo0: An LMC-mass galaxy that evolved in isolation for a majority of its life and experienced a recent burst in star formation.m11e: An irregular dwarf galaxy that spent most of its life as a compact dwarf, similar to the SMC in mass, until it experienced a recent major merger producing prominent stellar shells.m11i: A sub-SMC-mass galaxy that experienced multiple merger events and a bursty star formation history that has successively stripped gas from the galaxy.m10q: A low-mass dwarf galaxy that evolved in isolation and experienced a major blowout event.Figure 2. Evolution of galactic gas, stellar, and dust properties for our simulated galaxies for all gas/stars within 10 kpc of the galactic center.Specifically, total gas mass (top left), median gas metallicity (top right), total stellar mass (middle left), median D/Z ratio (middle right), star formation rate averaged in 10 Myr intervals (bottom left), and median silicate-to-carbonaceous dust mass ratio (bottom right).Our simulation suite covers a wide range of evolutionary histories, star formation rates, and metallicities.However, the resulting dust population evolution is relatively similar for all galaxies.Initially, the dust population is dominated by carbonaceous dust produced by SNe II, and later AGB stars, which results in a low D/Z ∼ 0.01 and Sil-to-C ∼ 0. Eventually, gas-dust accretion becomes efficient for first metallic iron dust, then silicates, and finally carbonaceous dust.The largest increase in D/Z occurs when silicate dust grows efficiently, as can be seen by the rapid increase in Sil-to-C due to silicates dominating the dust mass.When carbonaceous dust begins to grow efficiently, Sil-to-C correspondingly decreases and eventually settles at the typical MW value of ∼3.Note that m10q experiences a complete blowout event at ∼10 Gyr, evacuating almost all gas and dust from the galaxy.
composition exhibit little change between SNe and AGB-dominated dust regimes.

Onset of Accretion:
Once the gas within a galaxy reaches a 'critical' metallicity threshold (i.e.key element number abundance), dust growth via gas-dust accretion becomes more efficient than dust destruction by SNe shocks.The exact 'critical' metallicity varies between dust species (metallic iron, silicates, carbonaceous) due to differences in their key element abundances, physical properties, and life cycle processes, which we discuss in Sec.4.1.Metallic iron is the first dust species to grow efficiently via accretion, then silicates, and last carbonaceous dust.Each successively increases the median D/Z, up to D/Z ∼ 0.02, ∼ 0.2, and ∼ 0.3 respectively.

Build Up & Equilibrium:
Once a dust species begins to grow efficiently via accretion, its dust mass builds up over time until an m12i m12i m12i m12i Median values are given for all gas within 10 kpc of the galactic halo (solid), cold neutral gas (dashed), warm neutral gas (dotted), and warm/hot ionized gas (dash-dotted).The dust population of all galaxies is initially dominated by carbonaceous dust from SNe II stars, but this produces an extremely low D/Z.Carbonaceous dust from AGB stars eventually takes over as the dominant dust creation source as long as accretion remains inefficient due to low galactic metallicity, but D/Z changes little compared to the SNe II stardust-dominated population.If the galactic metallicity passes a critical threshold, accretion becomes efficient, becoming the dominant dust creation source and increasing the median D/Z in cool, dense gas.The onset of rapid growth via accretion for metallic iron, silicates, and carbonaceous dust species can be seen in the rapid increase of their corresponding species mass fraction.Overall, the growth of silicate dust via accretion produces the largest change in median D/Z.
equilibrium is reached, which we discuss in Sec.4.2.Once all dust species grow via accretion, the dust population of the galaxy reaches saturation, with the maximum amount of each dust species effectively forming.This produces an equilibrium dust population similar in composition to the Milky Way (D/Z ∼ 0.4 and Sil-to-C ∼ 3:1) and is relatively constant over time.
To gauge the accuracy of this predicted trend in dust population evolution, we must compare to observations of galactic dust population amount and chemical composition and observed dust population variability within galaxies.

Element Depletions & Aggregate D/Z
Element depletions provide a detailed accounting for how much of each element is locked in dust and are a strong constraint for our model.Observationally, element depletions are measured using UV spectral absorption features from sight lines to bright background sources, usually O/B type stars or quasars.Sight line column densities for ionized refractory elements and neutral hydrogen are determined from fits to their absorption profiles, and relative abundances between each element and hydrogen are determined.These relative abundances are then compared to known reference abundances and any missing elements from the gas-phase are then assumed to be m11e m11e m11e m11e Note that bursty star formation and resulting stellar feedback in the low-mass dwarf galaxies evacuates a majority of their cold gas which produces sudden increases in Z and decreases in dust for only cold gas.
locked in dust.The gas-phase depletion of element X is represented logarithmically as and linearly as where  X and  H are the gas-phase column density of element X and column density of neutral hydrogen (  (Jenkins & Wallerstein 2017) show similar trends to those in the MW but are offset correlating with the relative differences in metallicity between the galaxies 11 Element depletions have been observed for damped Ly systems (DLAs) via quasar absorption lines (Péroux & Howk 2020).However, DLAs probe a variety of systems (Prochaska & Wolfe 1997;Wolfe et al. 2005 To compare directly to these observations, we created a set of sight lines for each galaxy that originate from young star particles (formed within <10 Myr).For m12i, the sight lines terminate within the galactic disk at a distance of 0.1 − 2 kpc from the star to simulate depletion observations within the Milky Way.For the dwarf galaxies (m11v_halo2, m11d, m11e, m11v_halo0, and m11i), the sight lines terminate ∼50 kpc outside the galaxy to simulate observations of the LMC and SMC.For each sight line we then calculated  H,neutral and C, O, Mg, Si, and Fe depletions.The specifics of our sight line methodology are described in detail in Appendix B. We binned these sight lines in logarithmic  H,neutral bins and calculated the median values and 16-/84-percentiles for C, O, Mg, Si, and Fe depletions.The resulting relation between sight line element depletion and  H,neutral for each element can be seen in Fig. 5 for each galaxy.We include the 243 observed Milky Way sight line depletions from Jenkins (2009) and the carbon depletions from Parvathi et al. (2012), the 32 sight line depletions of the LMC from Roman-Duval et al. (2022a), and the 18 sight line depletions of the SMC from Jenkins & Wallerstein (2017).We modify the carbon depletions from Jenkins (2009) and Parvathi et al. ( 2012) similar to Sec. 3.2.1 in C22 to correct for differences in depletion estimates determined from strong and weak C ii transition lines and total C abundances.We also include a range of expected maximum C depletions in dense environments based on observations of 20% to 40% of C in CO in the Milky Way (e.g.Irvine et al. 1987;van Dishoeck et al. 1993;van Dishoeck & Blake 1998;Lacy et al. 1994).
We can also directly measure the depletion for each element as a function of physical gas density  H,neutral since our simulations track the total abundance of each element locked in dust for all gas cells.This is valuable for understanding sight line depletions since individual sight lines probe various gas phases and sight lines with similar  H,neutral can probe vastly different gas environments.We therefore bin the gas cells in logarithmic neutral gas density and calculate the median values and 16-/84-percentiles for C, O, Mg, Si, and Fe depletions.The resulting relation for each element depletion and neutral gas density,  H,neutral , is shown in Fig. 6 for each galaxy.We also include fits to the sight line depletion trends in the Milky Way from Jenkins (2009) assuming mean sight line density is the physical density as a lower bound, and using Zhukovska et al. (2016Zhukovska et al. ( , 2018) ) mean sight line density to physical density fit (see C22 Sec.

for details).
O, C, Mg, Si, and Fe show similar depletion trends for all galaxies, transitioning from a shallow to steep slope with increasing density ( H,neutral and  H,neutral ) similar to observations.This is due to the cycling of gas into cold, dense regions, where metallic iron, silicates, and/or carbonaceous dust grow via gas-dust accretion, and out to hot, diffuse regions where all species are destroyed via sputtering and SNe shocks similar to the results found in C22 for an idealized MW.The trends for C and Mg, in particular, flatten at high density ( H,neutral > 10 2 cm −3 ; not visible at high  H,neutral ) for all galaxies but for different reasons.C depletions flatten due to our sub-resolution dense molecular cloud prescription.This prescription assumes gas-phase C is rapidly converted to CO in dense molecular gas, halting carbon dust growth.Mg depletions flatten due to the relative abundances of Si and Mg and silicate dust being the only depletion source of Mg in our model.Si is the key element for silicate dust for all galaxies, so there will always be some leftover gas-phase Mg once the maximum amount of silicate dust has formed 12 .
Comparing predicted trends between galaxies, all depletions exhibit a staggered offset from one another.These offsets roughly correlate with each galaxy's median gas-phase metallicity, with m11i having the weakest and m12i having the strongest depletion trends, similar to what is seen in the MW, LMC, and SMC.These offsets are caused by differences in the average gas-dust accretion rate for each galaxy, which scales with gas-phase metallicity (see Eq.20 in C22).Galaxies with higher average gas-phase metallicity will have higher gas-dust accretion rates in dense environments and, therefore, stronger element depletions in these environments.C depletion trends exhibit the largest offset between galaxies and decreasing slopes at high densities.This decrease is so dramatic that only m12i, m11d, and m11v_halo2 produce noticeable increases in depletion with increasing density while m11e, m11v_halo2, and m11i have entirely flat relations.This is due to the turn-off of efficient carbonaceous dust growth via accretion.Galaxies with metallicities below a 'critical' metallicity threshold needed for efficient carbonaceous dust growth element for silicates.This was due to the idealized nature of the simulations, whose initial conditions assumed uniform gas and stellar metallicities of  =  ⊙ and a uniform stellar age distribution over 13.8 Gyr.The metal abundances arising from such a simulation are invariably different from those produced by a realistic stellar population successively built over the age of the universe.
will produce flat depletion relations set by the creation of dust in SNe and AGB winds.This can also be seen in the flat O, Mg, and Si depletion relations for m11i due to the lack of silicate dust growth via accretion.
We further aggregate the element depletions into total D/Z to determine the resulting distribution of D/Z with  H,neutral and temperature for each galaxy.We binned the gas cells in logarithmic neutral gas density and temperature bins, calculating the median D/Z and 16-/84percentiles shown in Fig. 7. Aggregating the observed Milky Way element depletions, we include an upper bound on the expected D/Z relation with  H,neutral along with an estimate for the WNM13 based on depletions from Jenkins (2009).We also include a reasonable estimate of D/Z in dense environments following the relation from Zhukovska et al. (2016Zhukovska et al. ( , 2018)).
The exhibited D/Z relations with  H,neutral mirror the individual element depletion trends for each galaxy, transitioning from a shallow to steep slope with increasing density and offset from one another in proportion to each galaxy's relative metallicity.This transition is less prominent in m12i, but it matches reasonably well with observations of the Milky Way.In regards to temperature, the resulting D/Z relation for all galaxies can be broken into three regimes.(i) Cool/Warm Gas (T < 10 4 K): As hot gas cools, gas-dust accretion eventually 'turns on' and dust begins to grow, steadily increasing D/Z with decreasing temperature.Our model assumes this 'turn on' point occurs when  < 300K, which can be seen as a moderate increase in D/Z.(ii) Ionized Gas (T = 10 4 − 10 5 K): This regime is primarily comprised of gas that has recently experienced SNe feedback, with dust residing in gas cells immediately around the SNe being destroyed.This produces a relatively flat D/Z for all galaxies, which lies at ∼40% of the maximum D/Z at low temperatures for each galaxy.This matches surprisingly well with the set dust destruction efficiency of ε ≈ 40% assumed for gas shocked to   ≥ 100 km/s in our SNe dust destruction routine.This indicates that, on average, our routine predicts all gas in this regime has been shocked to   ≥ 100 km/s once.However, the large scatter in D/Z indicates this varies considerably between gas cells, with some being shocked multiple times and others not shocked at all.(iii) Coronal Gas (T > 10 5 K): Dust destruction via thermal sputtering becomes efficient in this regime, with D/Z rapidly dropping with increasing temperature.For gas with  ≳ 10 6 K, on average, all dust is destroyed.

Extragalactic Dust Emission
Direct measurement of element depletions in galaxies beyond the Milky Way and its satellites is currently not possible, aside from DLAs.An alternative method to probe galactic dust populations is through the combination of multi-wavelength estimates of dust mass, gas mass, and metallicity.This method infers a total dust mass by fitting dust emission models to observed dust emission spectra.The assumed dust model can vary in complexity from a simple modified black body to physical dust models with specific dust grain sizes and chemical compositions, affecting the total dust mass inferred (Chastenet et al. 2021).This method was applied to a large sample of galaxies by Rémy-Ruyer et al. (2014), who studied the galaxy-integrated D/Z for 126 local galaxies, including dwarf, spiral, and irregular type galaxies, covering a large metallicity range (12 + log 10 (O/H) = 7.1 − 9.1;  ≈ 0.02 − 2.3 Z ⊙ ).(left) and temperature (right) for gas within 10 kpc for each galaxies at  = 0.All galaxies with efficient dust growth via accretion produce a sloped relation which increases with increasing  H,neutral and decreasing temperature.All galaxies with efficient silicate growth via accretion also have sizable D/Z even in low density, hot gas ( > 10 4 K or  H,neutral < 1 cm −3 ), which is exposed to SNe shocks.For gas with  > 10 5 K, D/Z quickly drops due to the onset of efficient dust destruction via thermal sputtering.

Galaxy-Integrated D/Z
We first show the resulting evolution between galaxy-integrated D/Z and metallicity for each simulated galaxy in Fig. 8 2019)14 along with their fitted relations.We define the galaxy-integrated D/Z as the median D/Z for cool, neutral ( < 1000 K) gas, following the assumption that galaxy-integrated dust emission studies primarily probe dense environments.In particular, Aniano et al. (2012Aniano et al. ( , 2020) ) directly compared observations of spatially-resolved and galaxy-integrated dust emission SEDs for individual galaxies and found that a majority of galaxy-integrated dust emission is produced by dust in dense regions exposed to diffuse radiation fields.Chiang et al. (2021Chiang et al. ( , 2023) also found that regions with good SNR (high dust emission) correlate with regions of CO detection (molecular gas), which further supports this assumption.We define the galaxy-integrated metallicity as the median 12 + log 10 (O/H) for gas with 7000 <  < 15000 and  H > 0.5 cm −3 to match the properties of nebular regions typically probed by empirical strong emission line methods used in these studies (e.g.Pilyugin & Grebel 2016).We also account for the depletion of O into dust by only considering gas-phase O instead of total (gas+dust) O abundance and include a -0.2 offset to correct for differences in reference O abundances assumed in our simulations (Anders & Grevesse 1989) and observations (Asplund et al. 2009).
We investigate other definitions of 12 + log 10 (O/H) in Appendix C and find that accounting for only gas-phase O has the largest effect, but this only produces an offset of ≲ 0.2 primarily for high values of 12 + log 10 (O/H).However, this depends entirely on our model's prescription for O depletion into dust.We also include the predicted relation between D/Z and metallicity from the equilibrium chemical and dust analytical model presented in Feldmann (2015), which tracks the evolution of a single, chemically-ambiguous dust species and accounts for stardust production, dust growth in the ISM, dust destruction by SNe shocks, and dust dilution by inflowing, pristine gas.We specifically plot the results from their model with slight et al. ( 2019) unique definition of D/Z = Σ dust /(Σ dust + Σ metals ) to account for possible depletion of metals into dust.
adjustments to their fiducial parameters.We set the maximum depletion limit for all metals to  dep = 0.6 to match the maximum D/Z predicted by our model.We modify the dust injection via AGB stars and SNe II of a stellar population to  D = 2 × 10 −3 to match our typical D/Z for stardust dominated galaxies.Finally, we set the ratio between the molecular gas depletion by star formation timescale and the dust growth via accretion timescale to  = 1.3 × 10 4 to match our model's predicted rapid increase in D/Z above a critical metallicity threshold.
At low metallicity, all of our galaxies exhibit similar D/Z relations.These galaxies are dominated by stardust production, which yields very low D/Z ∼ 10 −2 that is roughly constant with metallicity.Metallic iron dust growth becomes efficient for galaxies above 12 + log 10 (O/H) ∼ 7.6, increasing D/Z by roughly a factor of two.Above 12 + log 10 (O/H) > 7.9, the relations begin to diverge due to the onset of efficient silicate and later carbonaceous dust growth.Some galaxies (m11v_halo2 and m11v_halo0) follow the predicted equilibrium relation from Feldmann (2015) to surprising accuracy.Others (m12i, m11d, and m11e) fall below this relation to varying degrees and at different metallicities, producing a maximum scatter in D/Z of ∼0.5 dex at 12 + log 10 (O/H) ∼ 8.3.These varying relations are caused by delays in the buildup of dust mass once gas-dust accretion becomes efficient, which we discuss in Sec.4.2.At high metallicity, the relations begin to reconverge onto the equilibrium track where D/Z saturates at ∼0.4.At these metallicities, dust growth is so strong that a majority of all refractory elements have been locked into dust.Another interesting feature is the apparent 'backslides' in metallicity shown by m12i.These correspond to flyby and merger events which cause the infall of low-metallicity gas.This infalling gas decreases the metallicity of warm/hot gas, as can be seen in Fig. 3, including the nebular gas we use for our 12 + log 10 (O/H) definition.Even with our small sample of galaxies, we are able to recreate a reasonable amount of the scatter in D/Z for certain metallicities, which could explain some of the observed scatter15 .However, compared to the observed sample from Rémy-Ruyer et al. (2014) Figure 8. Resulting relation between galaxy integrated D/Z and metallicity for our galaxies over time.Each connected point is +1 Gyr apart in time, starting at ∼3 Gyr up to present day.The galaxy-integrated D/Z is set as the median D/Z for cool, neutral ( < 1000 K) gas, similar to the type of environments probed by observations of dust SEDs.The metallicity is set as the median gas-phase oxygen abundance (only O not depleted into dust) for gas with 7000 <  < 15000 and  H > 0.5 cm −3 , similar to regions probed by nebular emission lines used in galactic observations.12 + log 10 (O/H) is also offset by -0.2 to account for the differences in assumed oxygen abundances between Anders & Grevesse (1989) and Asplund et al. (2009).We compare with galaxy-integrated observations of local galaxies from Rémy-Ruyer et al. ( 2014) (orange circles) and De Vis et al. (2019) (gold triangles), along with the fits presented in both studies (dashed lines).We also include predictions from the equilibrium chemical and dust model presented in Feldmann (2015) (blue dash-dotted).All galaxies start with low D/Z and Z, being dominated by stardust production.For galaxies that surpass 12 + log 10 (O/H) ∼ 7.9, D/Z rises steeply due to the onset of efficient silicate, and later carbonaceous, dust growth via accretion.The differences in the steepness of this relation are due to the time it takes for the dust population to build up through successive cycling of gas into and out of cold, dense clouds.m12i also experiences a temporary decrease in the metallicity due to a minor merger which decreases the median metallicity of warm gas as seen in Fig. 3. appear to be shifted to higher metallicity, indicating our 'critical' metallicities may be too high.

Spatially-Resolved D/Z and D/H
We further compare our simulations at  = 0 to the spatially resolved observations of Chiang et al. (2023)  16 in Fig. 9, examining the relation between D/Z and inclination corrected galactocentric radius, neutral gas surface density (Σ gas,neutral ), molecular gas surface density (Σ H 2 ), stellar surface density (Σ star ), star formation rate surface density (Σ SFR ), and 12 + log 10 (O/H) gas 17 .We limit our comparison to the 17 galaxies in their 2 kpc resolution sample that have direct metallicity measurements and the derived D/Z values using the Bolatto et al. ( 2013)  CO prescription ( B13 CO ), which Chiang et al. ( 2021) argued yields the most reasonable D/Z.To match the observational resolution, we bin each simulation in 2 kpc face-on square pixels and calculate D/Z = Σ dust /Σ metals for each pixel.We then bin these 16 Note that the specific  B13 CO data products shown here are not explicitly presented in Chiang et al. (2023), but can be found in Chiang (2021). 17Our spatially-resolved O abundances follow the same definition as defined for galaxy-integrated O abundances.pixels across each property and calculate the median D/Z values and 16-/84 percentiles for each 18 .
All galaxies produce a weakly slopped relation between D/Z and all galactic properties.Specifically, increasing D/Z with decreasing galactocentric radius and increasing density and metallicity.All relations are roughly offset downwards by each galaxy's relative metallicity, besides the relation between D/Z and 12 + log 10 (O/H), which suggest all galaxies follow the same trend across their respective metallicity ranges.Compared to observations, m12i falls near the middle of the observed range for all galactic properties.The dwarf galaxies fall either at the bottom or below the range of observations, but the observational sample includes only one dwarf galaxy (M33/NGC598), so this is not unexpected.While our simulations disagree with the observed trends of decreasing D/Z with increasing Σ gas,neutral and Σ H 2 for many galaxies, this decreasing trend is likely due to observational biases at low column densities which can increase the expected D/Z.Specifically, CO-dark H 2 gas is not ac-18 While the observations from Chiang et al. (2023) primarily probe H 2dominated regions (employ a minimum  CO detection threshold), we do not include a molecular gas mass fraction (  H 2 ) cutoff for simulated pixels.We find that including such a cutoff has little effect on the resulting relations for m12i, primarily truncating the relations at lower surface densities and 12 + log 10 (O/H), but excludes the majority of pixels for our dwarf galaxies., star formation rate surface density (bottom middle), and gas-phase oxygen abundance (bottom right) in 2 kpc face-on bins for our galaxies at  = 0. We compare with dust emission-based observations of spatially-resolved D/Z of local galaxies from Chiang et al. (2023), which is an extension of the technique used in Chiang et al. (2021), with 2 kpc resolution and  B13 CO conversion factor.All of our simulated galaxies produce a weakly slopped D/Z relation, which decreases with galactocentric radius and increases with density and metallicity.m12i falls near the middle of, and the dwarf galaxies fall at the bottom or below the observed range for all galactic properties, but this is not unexpected since the observed sample includes only one dwarf galaxy.Our results disagree with the observed relations for Σ gas,neutral and Σ H 2 , which suggest D/Z decreases in denser environments.However, these observations disagree with high-resolution observations (Clark et al. 2023), so this may be due to resolution effects.counted for, which will decrease the expected Σ gas,neutral /Σ H 2 and thus increase D/Z.Furthermore, SNR pixel cuts for dust emission will preferentially exclude low D/Z pixels for low Σ gas,neutral /Σ H 2 .We also highlight that high-resolution dust emission observations from Clark et al. (2023) exhibit the opposite trend, with D/Z increasing with Σ gas,neutral .
We similarly compare our simulations to high-resolution observations of dust-to-neutral-gas mass ratio (D/H neutral = Σ dust /Σ Hi+H 2 ) for Local Group galaxies from Clark et al. (2023) in Fig. 10.We utilize the same face-on projection technique outlined above but with 144 pc pixel sizes to match the lowest observational resolution in their sample.We also only consider pixels with a neutral gas mass fraction  neutral > 0.5 to avoid Hii-dominated regions since Clark et al. (2023) does not include Hii gas in their gas mass accounting.Overall, our simulations predict D/H neutral increases with surface density for galaxies in which silicate and/or carbonaceous gas-dust accretion is efficient and match the typical D/H neutral observed at high Σ gas,neutral for all observed galaxies19 .However, the steepness of this relation is underpredicted, with D/H neutral only increasing by a factor of ≲2 in our simulations compared to the >3 factor observed for all galaxies.This suggests that our model's SNe dust destruction prescription is too weak, leaving too much dust in diffuse gas.However, we cannot entirely rule out our model predictions since the observations have an average standard deviation of 0.8 dex across all galaxies.We also note that when compared to the expected D/H neutral trends derived from element depletions (again using inferred C and O depletions), Resulting D/H neutral versus Σ gas,neutral for 144 pc resolution pixels from face-on projections of each galaxy at  = 0.Only pixels with neutral gas mass fraction >50% are included since pixels dominated by ionized gas produce artificially high D/H neutral at low Σ gas,neutral .We compare to deprojected 'face-on' observations of M33, M31, LMC, and SMC from Clark et al. (2023), which vary in resolution from 14-144 pc.Note the error bars presented for these observations are the uncertainty of the median and not the scatter.The average standard deviation for D/H neutral across all galaxies is 0.8 dex, but well-behaved and Gaussian.We caution that these observations may underpredict D/H neutral , when compared to element depletion observations.Overall, our simulations predict D/H neutral increases with surface density and match the typical values observed at high Σ gas,neutral .However, they underpredict the steepness of this relation, suggesting our SNe dust destruction prescription may be too weak.
these observations underpredict D/H neutral for the SMC but agree with the LMC.Our simulations also indicate that efficient silicate dust growth may have occurred relatively recently in the SMC, given its extremely steep relation that overlaps with m11v_halo0 at high Σ gas,neutral and m11i at low Σ gas,neutral .We caution that our approximation of observables begins to break down at these high resolutions.Specifically, individual phases of the ISM are resolved, resulting in pixels that primarily probe diffuse, ionized gas.These pixels have misleadingly high D/H neutral at low neutral surface densities since a majority of dust resides in ionized gas.This produces an upswing in D/H neutral at the lower end of our predicted relation and is even more pronounced when we do not include the  neutral > 0.5 cutoff.Whether observations would detect these pixels due to their possibly weak IR emissions is unknown and would require the creation of mock dust SEDs using radiative transfer codes.Due to this complexity, we save a more detailed comparison to spatially-resolved dust emission studies for future work.

DISCUSSION
For all but the most metal-poor galaxies in our simulation suite, gasdust accretion is the dominant producer of dust mass for a majority of each galaxy's life.Therefore, understanding when gas-dust accretion becomes efficient within a galaxy and when that galaxy's dust population transitions from a low to high D/Z ratio are critical for investigating galactic dust evolution.However, these two transitions have typically been lumped together under the moniker of 'critical' metallicity (e.g.Zhukovska 2014;Feldmann 2015).Below we break down these two transitions in our simulations, redefining the meaning of 'critical' metallicity, and investigate their implications for observations.

Efficient Gas-Dust Accretion
The transition to efficient gas-dust accretion for a dust species occurs when, on average, more dust mass is created by accretion in dense regions than is destroyed by destruction processes (SNe, thermal sputtering, astration) or removed from the galaxy by outflows.While this implies that numerous, interwoven processes on both a local and galactic scale determine this transition, in practice, the galaxyaveraged metallicity is found to be the primary determinator.In particular, a 'critical' metallicity threshold has been proposed above which the predicted galactic D/Z rapidly increases due to accretion (e.g.Zhukovska 2014;Feldmann 2015;Triani et al. 2020).For our purposes, we define the 'critical' metallicity as the point at which the species mass fraction of a given dust species begins to increase as seen in Fig. 3 & 4. Using this definition, all of our simulated galaxies show the same 'critical' metallicity threshold, above which the mass of a given dust species begins to increase.They all also indicate that each dust species has its own 'critical' metallicity, with  crit ∼ 0.05Z ⊙ for metallic iron,  crit ∼ 0.2Z ⊙ for silicates, and  crit ∼ 0.5Z ⊙ for carbonaceous dust.Three factors cause these differences in 'critical' metallicity: (1) Key Elements: For a dust species with a set chemical composition, its accretion rate scales with the number abundance of that species' key element.For metallic iron this is Fe, for silicates this is Si 20 , and for carbonaceous this is C. Assuming solar abundances 20 Depending on the relative abundances of Mg and Si, Mg can be the key element for silicates.For our simulations, we find Si is the key element for nearly the entire history of each galaxy.(Anders & Grevesse 1989or Asplund et al. 2009), C is 10 times more abundant than either Si or Fe.Atomic C is also lighter, so it will have thermal velocities 1.5 and 2.2 times greater than Si and Fe respectively.If we assume dust grows purely from hard-sphere type encounters, and all dust species are the same in all other respects, then carbonaceous dust should grow ∼15 times faster than silicates and ∼22 times faster than metallic iron dust, which is similar to the predictions of Granato et al. (2021).
(2) Physical Properties: Differences in physical properties between dust species (i.e. grain sizes, geometry, grain charging, sticking efficiency) can alter their effective accretion rates.In particular, our model assumes different average grain sizes (⟨⟩ 3 ; see Eq. 3 in C22) for each dust species due to differences in grain size distributions and effects of Coulomb enhancement 21 in atomic and diffuse molecular gas.Overall this produces a relative difference in accretion timescales of 1:150:900 in atomic/diffuse molecular gas for metallic iron, silicates, and carbonaceous dust, respectively, and 1:10:10 in dense molecular environments.
(3) Life Cycles: Differences in dust life cycle processes (creation, growth, destruction) between species affect the net change in dust mass during one cycle of a gas parcel into and out of dense regions.For carbonaceous dust, our model accounts for the rapid formation of gas-phase CO in dense molecular environments, which halts dust growth.This decreases the amount of time carbonaceous dust has to grow compared to other dust species.For metallic iron, we assume a fraction of the dust population is locked inside silicate dust as inclusions (  incl = 0.7), protecting them from destruction by SNe and thermal sputtering and reducing the amount of dust grain surface area available for accretion.
We caution that the included physics in our dust evolution model are in no way complete and the predicted 'critical' metallicities could change with the incorporation of new physics.However, the predictions of varying 'critical' metallicities between dust species have some support from observations.As shown in Fig. 5, the LMC and SMC show large depletions of Si, Mg, and Fe in dense environments, suggesting the efficient growth of silicate and a theoretical iron-bearing dust species which agree with our model.While C depletions are currently unobservable in these galaxies, dust extinction and emission observations show decreasing amounts of small carbonaceous grains and PAHs (Weingartner & Draine 2001;Chastenet et al. 2019), suggesting carbonaceous dust growth may not be efficient in such galaxies.A sizable carbonaceous dust population dominated by large grains could exist and be indiscernible by current observational techniques.However, this seems unlikely as it would require carbonaceous grains to be far less prone to shattering into smaller grains compared to other dust species.

Dust Buildup and Equilibrium
Once gas-dust accretion becomes efficient for a given dust species, the dust mass will build through successive cycles of gas into and out of dense environments until it reaches an equilibrium between dust growth and dust destruction via SNe.This can be seen in Fig. 2, 3, & 4 as sharp increases in species mass fraction that eventually plateau (a similar but weaker trend can also be seen in D/Z).We label the time between the onset of gas-dust accretion and the equilibrium plateau as the 'equilibrium' timescale ( equil ) and find this timescale varies 21 Coulomb enhancement arises from a grain size dependent electrostatic enhancement factor which accounts for the change in interaction cross section between ionized gas-phase metals and charged dust grains.
between dust species and between galaxies, with  iron equil ∼ 0.75−4 Gyr,  sil equil ∼ 1.0 − 1.5 Gyr, and  carb equil ≳ 3 Gyr for metallic iron, silicates, and carbonaceous dust respectively.The variation in  equil between dust species is due to the differences in dust life cycle processes outlined in Sec.4.1, while the range in  equil for a given dust species is due to differences in gas cycling and SNe destruction timescales between galaxies.We can confirm this with predictions of  equil from analytical models (Eq.20 and Eq.21 in C22).In Fig. 11, we show the predicted evolution of the fraction of Si locked in silicate dust (  Si ) from this model.We assume an initial degree of condensation of  o,Si = 0.01, cold cloud lifetime of  cloud ∼ 10 Myr estimated from Milky Way-mass galaxy simulations in Benincasa et al. (2020), and a constant accretion growth timescale of  grow ∼ 40 Myr which is the median value for gas in which accretion can occur in our model (T<300 K) for all galaxies in our sample at the onset of efficient silicate dust growth.We vary the SNe dust destruction timescale ( d ; Eq. 18 in McKee 1989) and the fraction of the ISM with gas  < 300 K ( cloud ), which determines the average time it takes to cycle all ISM material from the cold cloud phase through the diffuse/warm ISM phases and back into cold clouds  cycle =  cloud 1− cloud  cloud , to match typical values for dwarf, Milky Way, and massive/starburst galaxies.By varying these two parameters, we predict a spread in  equil similar to those exhibited in our simulations.However, these predictions are an upper bound on  equil , especially for galaxies with high star formation rates, since  grow inversely scales with metallicity and will thus decrease over time.We, therefore, also consider a simple case of decreasing  grow /increasing metallicity.If we assume the metal mass of a galaxy only changes due to metal injection from SNe and metal depletion via astration (metal locked up in stars) and a constant galactic gas mass, with no inflows or outflows, then the galactic metallicity will have the general form  ∼  −  exp −/ , where , , and  are constants that depend on the SNe metal mass yields, star formation rate, and galactic gas mass.We show the results for the expected  equil with this decreasing  grow in Fig. 11 choosing values of  and  such that the growth timescales at  = 0 and  = ∞ are  grow,o = 40 Myr and  grow,∞ = 10 Myr respectively, and  = 5 Gyr.Overall, this inclusion only slightly decreases the variation in  equil but still overlaps with our simulation predictions.We caution that our simulations' resulting equilibrium timescales depend on our model's prescriptions for accretion and SNe dust destruction.In particular, our accretion routine assumes a constant size distribution, which will underpredict accretion timescales, and the theoretical predictions for SNe dust destruction efficiency vary substantially (e.g.Hu et al. 2019;Kirchschlager et al. 2019).
The equilibrium timescale has important implications for the expected relation between galactic D/Z and metallicity when it is compared to the galactic metal enrichment timescale,  metal , which we define as the time it takes for the galactic metallicity to double.For  equil ≲  metal , the dust population within a galaxy is able to reach a relative equilibrium D/Z for a given metallicity before it changes appreciably.This is similar to the model predictions presented in Feldmann (2015), which assumes equilibrium galactic metallicity and D/Z.When  equil ≳  metal , the dust population is still in the process of building up its mass and will 'lag' behind the predicted equilibrium D/Z at a given metallicity.However, since  equil inversely scales with metallicity, as the galactic metallicity increases,  equil will eventually fall below  metal , causing these two trends to converge.This effect can be seen in the evolution between galactic D/Z and Z for each galaxy shown in Fig. 8, with the emergence of two evolutionary states we call 'equilibrium' and 'lagging'.The 'equilibrium' state is exhibited by galaxies that have long quiescent periods, evolving in isolation for the majority of their life and/or experiencing relatively constant and low star formation rates (i.e.m11v_halo2, m11v_halo0, and m11i) resulting in large  metal .These galaxies all follow the same relation between D/Z and metallicity, similar to the equilibrium predictions from Feldmann (2015).The 'lagging' state is exhibited by galaxies that have bursty episodes, experiencing major and/or multiple minor merger events and rapid changes in star formation rates (i.e.m12i, m11d, and m11e) resulting in instances of small  metal .These galaxies exhibit relations between D/Z and metallicity that occasionally fall below the 'equilibrium' trend, which roughly corresponds with recent bursty periods.These results suggest that recent galactic history plays a critical role in the expected galactic D/Z for galaxies in which gas-dust accretion is efficient and could explain a portion of the scatter in observed D/Z at a fixed metallicity.However, our sample of galaxies is small, and the one galaxy that exhibits the largest deviation from equilibrium predictions is the only Milky Way-mass galaxy in our sample.More simulations are therefore needed to make a more statistically robust conclusion.

CONCLUSIONS
In this work, we investigate the evolution of galactic dust populations, in both composition and amount, across cosmic time utilizing a suite of 7 cosmological zoom-in simulations run with the FIRE-2 model (Hopkins et al. 2018b) for stellar feedback and ISM physics and the "Species" dust evolution model (C22).This dust evolution model accounts for dust creation in stellar outflows, growth from gas-phase accretion, destruction from SNe shocks, thermal sputtering, and astration, and turbulent dust and metal diffusion in gas.It tracks the evolution of specific dust species (silicates, carbon, silicon carbide), treating each uniquely depending on their chemical composition, along with theoretical nano-particle metallic iron (Nano-iron) dust species and an oxygen-based (O-reservoir) dust species.It also incorporates a physically motivated dust growth routine which accounts for Coulomb enhancement and CO formation in dense molecular environments.
The 7 galaxies we selected cover a broad range of stellar (10 6  ⊙ <  * < 10 11  ⊙ ) and halo (10 9  ⊙ <  vir < 10 12  ⊙ ) masses at present day and showcase a variety of growth histories.We summarize our findings of dust population evolution below: (i) Despite the variety of galactic evolutionary histories probed in our sample, all galactic dust populations follow similar evolutionary trends.Initially, they are dominated by dust production from first SNe II and later (∼ 2 Gyr) AGB stars, resulting in a low D/Z ∼ 0.01.Above a 'critical' metallicity threshold, gas-dust accretion becomes efficient and eventually increases the galactic D/Z > 0.1 (Fig. 2, 3,  & 4).This suggests that gas-dust accretion is the main producer of dust mass for all but the most metal-poor galaxies and, in the case of Milky Way-mass galaxies, for a majority of the galaxy's life.
(ii) Differences between key element abundances, physical properties, and life cycle processes between dust species result in varying 'critical' metallicities of  crit ∼ 0.05Z ⊙ ,0.2Z ⊙ , and 0.5Z ⊙ for metallic iron, silicates, and carbonaceous dust respectively.These differences reproduce observed depletions of Si, Mg, and Fe in the MW, LMC, and SMC (Fig. 5 & 6), and suggest that silicate and a theoretical iron-bearing dust species grow efficiently by accretion in the LMC and SMC.They also suggest that C depletion decreases rapidly with galactic metallicity, which could explain the reduced amount of small carbonaceous grains observed in the LMC and SMC.
(iii) After the onset of efficient gas-dust accretion, there is a characteristic equilibrium timescale over which dust mass builds up over time until an equilibrium between dust growth and SNe dust destruction is reached (Fig. 11).This equilibrium timescale is ≳ 1 Gyr in our sample and varies between galaxies.For galaxies with quiescent histories (long metal enrichment timescales), their dust populations are almost always in relative equilibrium and produce similar galaxyintegrated D/Z-Z trends.For galaxies with bursty episodes due to mergers/rapid changes in SF (short metal enrichment timescales), their dust populations can 'lag' behind the expected equilibrium values for a time, producing lower D/Z.This 'lagging' state can explain part of the large scatter in the observed relation between galaxyintegrated D/Z and metallicity (Fig. 8).
(iv) Extragalactic observations of spatially-resolved D/Z in spiral galaxies are roughly consistent with our Milky Way-mass galaxy (Fig. 9).When compared to high-resolution observations of D/H neutral in Local Group galaxies, our model overpredicts the amount of dust in diffuse neutral gas (Fig. 10), but this may be due to our approximations of observables.A more direct comparison with mock dust emission SEDs is needed for a stronger test of our model.B2.Face-on gas density projection of m11v_halo0 (top left), m11d (top right), m11e (middle left), m11v_halo2 (middle right), and m11i (bottom) with all young stars (gold stars) overplotted.Three sight lines are created for each young star particle to simulate the view of an outside observer at different orientations.One is face-on with the average angular momentum vector of young stars in the galaxy, one is orthogonal to the angular momentum vector, and one is halfway between the two.Note that the number of young stars varies by up to 1 dex between the different galaxies however, they all roughly trace the densest parts of each galaxy.(right) Relation between oxygen abundances for definition (1) and ( 4), median gas-phase O abundance (accounting for O depleted into dust) for all gas with 7000 <  < 15000 K and  H > 0.5 cm −3 .The data points shown are the same as those in Fig. 8.Note the disagreement between definitions (1) and ( 2) at low 12 + log 10 (O/H) is due to our choice of only considering gas within 10 kpc of the galactic center.At early times, when 12 + log 10 (O/H) is low, 10 kpc encompasses a large amount of ionized halo gas for our dwarf galaxies.If we only consider gas within 0.1 vir this disagreement disappears entirely.

Figure 1 .
Figure1.Face-on images of our spiral and dwarf galaxies at  = 0 (note the differences in scale).Each image is a mock Hubble Space Telescope ugr composite accounting for dust attenuation utilizing the D/Z produced by our model along with assumed MW, LMC, or SMC dust opacities depending on the median gas cell Sil-to-C ratio.m12i: A compact, Milky Way-mass spiral galaxy that experienced multiple minor mergers during its life.m11d: An LMC-mass galaxy with bursty star formation and multiple mergers throughout its life.m11v_halo2: An LMC-mass galaxy with faint spiral arms that evolved in isolation for a majority of its life.m11v_halo0: An LMC-mass galaxy that evolved in isolation for a majority of its life and experienced a recent burst in star formation.m11e: An irregular dwarf galaxy that spent most of its life as a compact dwarf, similar to the SMC in mass, until it experienced a recent major merger producing prominent stellar shells.m11i: A sub-SMC-mass galaxy that experienced multiple merger events and a bursty star formation history that has successively stripped gas from the galaxy.m10q: A low-mass dwarf galaxy that evolved in isolation and experienced a major blowout event.

Figure 3 .
Figure3.Time evolution of median metallicity (top), D/Z (second from top), dust creation source mass fraction (second from bottom), and dust species mass fraction (bottom) for m12i (left), m11v_halo0 (middle), and m11d (right).Median values are given for all gas within 10 kpc of the galactic halo (solid), cold neutral gas (dashed), warm neutral gas (dotted), and warm/hot ionized gas (dash-dotted).The dust population of all galaxies is initially dominated by carbonaceous dust from SNe II stars, but this produces an extremely low D/Z.Carbonaceous dust from AGB stars eventually takes over as the dominant dust creation source as long as accretion remains inefficient due to low galactic metallicity, but D/Z changes little compared to the SNe II stardust-dominated population.If the galactic metallicity passes a critical threshold, accretion becomes efficient, becoming the dominant dust creation source and increasing the median D/Z in cool, dense gas.The onset of rapid growth via accretion for metallic iron, silicates, and carbonaceous dust species can be seen in the rapid increase of their corresponding species mass fraction.Overall, the growth of silicate dust via accretion produces the largest change in median D/Z.

Figure 4 .
Figure 4. Same as Fig.3for m11e (left), m11v_halo2 (middle left), and m11i (middle right) and m10q (right).Note that bursty star formation and resulting stellar feedback in the low-mass dwarf galaxies evacuates a majority of their cold gas which produces sudden increases in Z and decreases in dust for only cold gas.
; Faucher-Giguère & Kereš 2011; Rhodin et al. 2019; Stern et al. 2021) and their reference abundances can be difficult to determine due to dust depletions (Roman-Duval et al. 2022b).For these reasons, we forgo comparing to DLA depletion observations in this work.(Roman-Duval et al. 2022a).However, direct observations of C and O depletions are limited only to the Milky Way due to either saturated or extremely weak absorption lines, so C and O depletion trends for the SMC and LMC are usually inferred from the MW relation between Fe and C or O depletions respectively (Roman-Duval et al. 2022a,b).
Figure10.Resulting D/H neutral versus Σ gas,neutral for 144 pc resolution pixels from face-on projections of each galaxy at  = 0.Only pixels with neutral gas mass fraction >50% are included since pixels dominated by ionized gas produce artificially high D/H neutral at low Σ gas,neutral .We compare to deprojected 'face-on' observations of M33, M31, LMC, and SMC fromClark et al. (2023), which vary in resolution from 14-144 pc.Note the error bars presented for these observations are the uncertainty of the median and not the scatter.The average standard deviation for D/H neutral across all galaxies is 0.8 dex, but well-behaved and Gaussian.We caution that these observations may underpredict D/H neutral , when compared to element depletion observations.Overall, our simulations predict D/H neutral increases with surface density and match the typical values observed at high Σ gas,neutral .However, they underpredict the steepness of this relation, suggesting our SNe dust destruction prescription may be too weak.

Figure B1 .
Figure B1.Face-on (top) and edge-on (bottom) gas density projection of m12i with a subset of our sight lines (white lines) to young stars overplotted.The endpoints (purple crosses) of the sight lines are randomly selected within the disk.The start of the sight lines (gold stars) are up to 10 young star particles within 0.1 − 2kpc of the endpoints.Our sight lines probe stars both within and above the disk plane, voids, and dense spiral arms.

Figure C1 .
Figure C1.Resulting oxygen abundances from different definitions for our suite of simulations.(left) Relation between oxygen abundances for definition (1), median O abundance for all gas within 10 kpc, and definition (2), median O abundance for all gas with 7000 <  < 15000 K and  H > 0.5 cm −3 .(middle) Relation between oxygen abundances for definition (1) and (3), median gas-phase O abundance (accounting for O depleted into dust) for all gas within 0.1 vir .(right) Relation between oxygen abundances for definition (1) and (4), median gas-phase O abundance (accounting for O depleted into dust) for all gas with 7000 <  < 15000 K and  H > 0.5 cm −3 .The data points shown are the same as those in Fig.8.Note the disagreement between definitions (1) and (2) at low 12 + log 10 (O/H) is due to our choice of only considering gas within 10 kpc of the galactic center.At early times, when 12 + log 10 (O/H) is low, 10 kpc encompasses a large amount of ionized halo gas for our dwarf galaxies.If we only consider gas within 0.1 vir this disagreement disappears entirely.

Figure C2 .
Figure C2.Resulting relation between median O abundance and median metallicity for all gas within 10 kpc for our suite of simulations.The data points shown are the same as those in Fig. 8.Our simulations follow the O abundances expected from Anders & Grevesse (1989) which are offset by ∼ 0.2 dex from Asplund et al. (2009) abundances usually assumed by observations.
H,neutral =  H i + 2 H 2 ) respectively and ( X / H ) ref is the assumed reference abundance.Due to the high spatial and spectral resolution needed for such observations, they are mainly limited to Milky Way, LMC, and SMC 11 .Observations of Milky Way depletions compiled by Jenkins (2009) show C, O, Mg, Si, and Fe depletion increases with increasing neutral H column density.Mg, Si, and Fe depletion observations in the LMC (Roman-Duval et al. 2021) and SMC Zhukovska et al. (2016Zhukovska et al. ( , 2018) ),neutral for gas within the 10 kpc of each galaxy at  = 0.The 16-/84-percentile ranges are represented for each by shaded regions.We compare with observed elemental depletion relations in the Milky Way from Jenkins (2009) assuming mean sight line density is the physical density (black-dashed), which can be treated as a lower limit, and usingZhukovska et al. (2016Zhukovska et al. ( , 2018) )mean sight line density to physical density fit (black-solid).For O, Mg, Si, and Fe we include estimates for the WNM depletions (diamond) along with an interpolation to the Jenkins' relation (black-dotted).Si, O, and Fe which transition from a shallow to steep slope for increasing  H,neutral .Mg and C depletions show an additional transition to a flat slope at the highest densities.For Mg, this is due to the saturation of silicate dust, which has run out of available Si to grow further.For C, this is due to the rapid formation of CO, which takes up any remaining gas-phase C. m11i only exhibits a slopped trend for Fe since only metallic iron grows efficiently.
Clark et al. (2023)Chiang et al. ( , 2023) )  H,neutral , but each relation is offset roughly corresponding with each galaxy's median metallicity.The C relations exhibit a weaker transition and a larger offset, with the lowest metallicity galaxies (m11e, m11v_halo0, and m11i) producing a nearly flat relation.This is due to the 'turn-off' of carbonaceous growth via accretion since these galaxies lie below the carbonaceous 'critical' metallicity threshold.m11iexhibitsflatdepletionrelationsforallbutFesince only metallic iron dust grows efficiently in this galaxy.Overall, the predicted trends are a good match with observations, but the lack of C or O depletion measurements outside the Milky Way paints an incomplete picture.DeVis et al. (2019)further expanded upon this sample, adding ∼ 500 local DustPedia galaxies, which cover a more limited metallicity range (12 + log 10 (O/H) = 7.9 − 8.7;  ≈ 0.15 − 0.9 Z ⊙ ) and are mainly limited to spiral galaxies.Chiang (2021);Chiang et al. (2021Chiang et al. ( , 2023) )extended this method to ∼2 kpc spatially-resolved studies of individual galaxies, compiling spatially resolved relations between D/Z and local environments for 46 nearby galaxies.Their sample is limited to primarily spiral galaxies with  star > 10 9  ⊙ and 12 + log 10 (O/H) ≳ 8.4 ( ≳ 0.5 Z ⊙ ).Recently,Clark et al. (2023)utilized this method with improved Herschel data to produce high resolution (14 -144 pc) maps of dust-to-neutral-gas (D/H neutral ) ratios for Local Group galaxies (LMC, SMC, M33, and M31).
Predicted evolution of the fraction of Si locked in silicate dust through successive cycles of gas into dense clouds, where the dust grows via accretion, and out into diffuse clouds, where it is destroyed by SNe, using Eq.20 and Eq.21 in C22.We assume an initial degree of condensation  Si = 0.01, dust growth timescale  grow ∼ 40 Myr, and cold cloud lifetime  cloud ∼ 10 Myr (solid).We vary the SNe dust destruction timescale ( d ) and the fraction of the ISM in cold clouds ( cloud ), selecting values typical for (red) a MW-mass galaxy, (blue) a dwarf-mass galaxy, and (green) a massive/starburst galaxy.While the resulting equilibrium fractions are similar, the time to equilibrium varies dramatically  equil ∼ 0.4 − 3 Gyr, matching the spread in  equil seen in our simulations.However, these predictions are an upper bound on  equil since  grow will decrease over time as the metallicity of a galaxy increases.We, therefore, also consider the case where  grow decreases/metallicty increases (dashed).We assume the galactic metallicity has the general form  ∼  −  exp −/ , choosing  and  such that  grow,o = 40 Myr initially and  grow,∞ = 10 Myr at  = ∞, and  = 5 Gyr.This modification shortens  equil for each typical galaxy, but the variability remains large with  equil ∼ 0.3 − 2 Gyr.