Over the last decade the Atacama Large Millimeter Array ALMA has revealed massive extremely dusty starforming galaxies at zgtrsim5 and the James Webb Space Telescope JWST is primed to uncover even more information about them. These extreme observations both need dust evolution theory to provide context and are excellent benchmarks to test this theory. Here we investigate the evolution of galactic dust populations at cosmic dawn using a suite of cosmological zoomin simulations of moderately massive highredshift Mgtrsim109 Modot zgtrsim5 galaxies from the Feedback in Realistic Environments FIRE project the highest resolution of such simulations to date. 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 allowing it to replicate a wide range of presentday observations. We find similar to other theoretical works that dust growth via gasdust accretion is the dominant producer of dust mass for these galaxies. However our fiducial model produces Mrm dust that fall gtrsim1 dex below observations at any given M which we attribute to reduced accretion efficiencies caused by a combination of low metallicities and extremely bursty star formation in these galaxies. Modest enhancements i.e. within observationaltheoretical uncertainties to accretion and SNe II dust creation raise Mrm dust by lesssim1 dex but this still falls below observations which assume Trm dustsim25 K. One possibility is that inferred dust masses for zgtrsim4 galaxies are overestimated and recent observational works that find Trm dustsim50 K along with metallicity constraints tentatively support this.
['Astrophysics - Astrophysics of Galaxies']
A Dusty Dawn Galactic Dust Buildup at zgtrsim5
{'A Dusty Dawn: Galactic Dust Buildup at 𝑧 ≳ 5': 'Caleb R. Choban 1 ★ , Samir Salim 1 , Dušan Kereš 2 , 3 , Christopher C. Hayward 4 , and Karin M. Sand-  \nstrom 2  \n- 1 Department of Astronomy, Indiana University, Bloomington, IN 47405, USA\n- 2 Department of Astronomy & Astrophysics, University of California at San Diego, La Jolla, CA 92093, USA\n- 3 Department of Physics, University of California at San Diego, La Jolla, CA 92093, USA\n- 4 Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010 USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'Over the last decade, the Atacama Large Millimeter Array (ALMA) has revealed massive, extremely dusty star-forming galaxies at 𝑧 ≳ 5, and the James Webb Space Telescope (JWST) is primed to uncover even more information about them. These extreme observations both need dust evolution theory to provide context and are excellent benchmarks to test this theory. Here, we investigate the evolution of galactic dust populations at cosmic dawn using a suite of cosmological zoom-in simulations of moderately massive, high-redshift ( 𝑀 ∗ ≳ 10 9 𝑀 ⊙ ; 𝑧 ≳ 5) galaxies from the Feedback in Realistic Environments (FIRE) project, the highest resolution of such simulations to date. 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, allowing it to replicate a wide range of present-day observations. We find, similar to other theoretical works, that dust growth via gas-dust accretion is the dominant producer of dust mass for these galaxies. However, our fiducial model produces 𝑀 dust that fall ≳ 1 dex below observations at any given 𝑀 ∗ , which we attribute to reduced accretion efficiencies caused by a combination of low metallicities and extremely bursty star formation in these galaxies. Modest enhancements (i.e., within observational/theoretical uncertainties) to accretion and SNe II dust creation raise 𝑀 dust by ≲ 1 dex, but this still falls below observations which assume 𝑇 dust ∼ 25 K. One possibility is that inferred dust masses for 𝑧 ≳ 4 galaxies are overestimated, and recent observational works that find 𝑇 dust ∼ 50 K along with metallicity constraints tentatively support this.  \nKey words: methods: numerical - dust, extinction - galaxies: evolution - galaxies: ISM', '1 INTRODUCTION': 'Until recently, our understanding of dust and its abundance in the early universe was \'shrouded\' in mystery. The advent of the Atacama Large Millimeter Array (ALMA) in the last decade has parted the proverbial curtain on rest-frame far-infrared (FIR) light at this epoch, allowing for unprecedented detection of dust continuum sources at 𝑧 > 4. In the succeeding period, numerous massive, extremely dusty star-forming galaxies (DSFGs) have been discovered at 𝑧 ≳ 5, starting with individual detections (Vieira et al. 2013a; Hezaveh et al. 2013; Watson et al. 2015; Spilker et al. 2016; Laporte et al. 2017; Strandet et al. 2017; Miller et al. 2018; Tamura et al. 2019; Reuter et al. 2020; Witstok et al. 2023b) and recently expanding to ∼ 100 galaxies with the ALMA Large Program to INvestigate [C /i.pc/i.pc] at Early Times (ALPINE; Le Fèvre et al. 2020; Béthermin et al. 2020; Faisst et al. 2020) 𝑧 ∼ 5 survey and the Reionization Era Bright Emission Line Survey (REBELS; Bouwens et al. 2022) 𝑧 ∼ 7 survey. However,  \nthese samples are primarily UV-selected 1 and are, therefore, biased towards a subset of UV-bright galaxies. Furthermore, the information these surveys provide is primarily restricted to estimates of galactic dust and stellar mass. These limitations make determining the exact process responsible for such dusty galaxies difficult. In particular, both observations (Rémy-Ruyer et al. 2014; De Vis et al. 2019) and simulations, which we discuss later, find the dust content of local galaxies has a strong correlation with galactic metallicity. Therefore, any predictions of dust evolution in these high-z galaxies are largely predicated on their assumed chemical evolution (e.g. Palla et al. 2020, 2024).  \nFortuitously, the recently launched James Webb Space Telescope (JWST) probes the rest-frame optical at 𝑧 > 4, opening another observational window for high-z DSFGs. Notably, JWST can resolve prominent nebular emission lines, providing measurements of galactic gas-phase metallicity. This has already led to estimates of the  \n1 ALMA has also observed galaxies detected in the submillimetre by the South Pole Telescope (Vieira et al. 2013b; Hezaveh et al. 2013; Spilker et al. 2016; Miller et al. 2018; Reuter et al. 2020), finding they are generally the most massive DSFGs (median SFR ∼ 10 3 𝑀 ⊙ / yr and 𝑀 dust ∼ 10 9 𝑀 ⊙ ). However, these have no accompanying stellar mass estimates.  \ngas-phase mass-metallicity relation, which suggests high-z galaxies are more metal-poor than their present-day stellar mass counterparts (Nakajima et al. 2023; Curti et al. 2024; Chemerynska et al. 2024) 2 . Furthermore, JWST can resolve attenuation curves, revealing information on dust population composition, such as the presence of the 2175 Å bump produced by small carbonaceous grains (Witstok et al. 2023a). JWST can also identify extremely red galaxies that are heavily dust-obscured at rest-frame UV wavelengths (Akins et al. 2023), avoiding the UV-selection bias inherent in many previous observations.  \nGiven current and forthcoming DSFGs observations, dust evolution theory needs to be taken to task to provide context for observations, further our understanding of the dust life cycle, and make subsequent predictions that can be followed up with JWST. Indeed, not long after the first detections of high-z DSFGs, analytical models showed that the creation of dust by supernovae (SNe) alone 3 cannot produce the dust masses observed due to the subsequent destruction of dust by supernovae shocks (e.g. Michałowski 2015; Leśniewska &Michałowski 2019). In the local universe, there is ample evidence that preexisting dust grains can grow from the accretion of gas-phase metals (gas-dust accretion). Notably, observations of the Milky Way (MW) and nearby galaxies find the fraction of metals locked in dust (dust-to-metals ratio; D/Z) increases with local gas surface density (Jenkins 2009; Roman-Duval et al. 2021, 2014, 2017; Chiang et al. 2018; Clark et al. 2023). Given this, gas-dust accretion is generally believed to be responsible for the large dust masses seen in high-z DSFGs, but our understanding of this process is limited due to the inherent difficulties of experimental study.  \nIn recent years, dust evolution models integrated into semianalytical and cosmological simulations have been utilized to further our understanding of gas-dust accretion and all other processes in the dust life cycle (e.g. Bekki 2015; McKinnon et al. 2016; Li et al. 2019; Granato et al. 2021; Choban et al. 2024). Despite variations in methodologies and included physics these models agree in broad strokes, finding that dust growth via accretion is responsible for the bulk of dust content of the MW and local galaxies (e.g. McKee 1989; Draine 2009). However, details in the dust life cycle important for the predictions of high-z DSFGs vary considerably between works. This is due to a focus on present-day observations, which are the result of dust life cycle processes integrated over ∼ 13 Gyrs of evolutionary history. In particular, simulations predict that accretion becomes efficient (more dust is grown than is destroyed) within a galaxy once that galaxy\'s mass-weighted median metallicity passes a certain threshold (critical metallicity threshold; 𝑍 crit ) (Inoue 2011; Asano et al. 2013; Zhukovska 2014; Feldmann 2015; Popping et al. 2017; Hou et al. 2019; Li et al. 2019; Graziani et al. 2020; Triani et al. 2020; Parente et al. 2022; Choban et al. 2024), but the predicted thresholds vary considerably ( 𝑍 crit ∼ 0 . 03 -0 . 5 Z ⊙ ). Works that track the evolution of chemically distinct dust species further suggest a separate 𝑍 crit for each dust species (Granato et al. 2021; Parente et al. 2022; Choban et al. 2024). There is also a general consensus that accretion steadily increases the galactic D/Z, and dust mass, up to an equilibrium value once 𝑍 crit is reached. However, the exact timescale of this buildup varies between galaxies, depending on the ISM phase structure, star formation rate (SFR), and the initial amount of dust (e.g. Zhukovska  \net al. 2008; Zhukovska 2014), and could be ≳ 1 Gyr for a MW-mass galaxy (Choban et al. 2024). Furthermore, the assumed efficiency of dust creation by SNe II varies considerably between models since it is found to have little effect on present-day observations due to the dominance of accretion.  \nFewer works have focused on dust evolution in high-z DSFGs, utilizing semi-analytical models (Popping et al. 2017; Vĳayan et al. 2019; Triani et al. 2020; Dayal et al. 2022; Mauerhofer & Dayal 2023) and cosmological simulations (Graziani et al. 2020; Esmerian & Gnedin 2022, 2024; Lewis et al. 2023; Di Cesare et al. 2023; Lower et al. 2023, 2024), and their predictions vary considerably. While many authors agree that accretion is the dominant producer of dust mass at this epoch, others predict accretion is inefficient and that SNe II dominate dust production (Triani et al. 2020; Dayal et al. 2022). Furthermore, some models struggle to reproduce the high dust masses observed for high-z DSFGs (Vĳayan et al. 2019; Triani et al. 2020; Dayal et al. 2022), while others reproduce them with relative ease (Popping et al. 2017; Di Cesare et al. 2023; Lewis et al. 2023; Lower et al. 2023). A major limitation of the above-mentioned simulations is their inability to resolve the multi-phase ISM and self-consistently model dust life cycle processes. In particular, they utilize \'tuned\' sub-resolution prescriptions for star formation, affecting the predicted galactic star formation history (e.g. Iyer et al. 2020). They also use sub-resolution gas-dust accretion routines and do not track the evolution of chemically distinct dust species, affecting the predicted dust population evolution and spatial variability within galaxies (Choban et al. 2022, 2024).  \nThere has also been a growing trend of post-processing high-z simulations originally run without an explicit dust evolution model by using radiative transfer codes and assumed dust populations to make observational predictions for ALMA and JWST (Cochrane et al. 2019; Ma et al. 2019; Liang et al. 2018, 2019; Vogelsberger et al. 2020; Shen et al. 2020; Parsotan et al. 2021; Shen et al. 2022; Vĳayan et al. 2022; Cochrane et al. 2022, 2023, 2024). The dust populations assumed are typically derived from MW observations and are relatively simplistic (e.g. D/Z = 0 . 4 for gas with 𝑇 ≲ 10 6 K). However, high-z dust populations can be drastically different from the MW, affecting predicted observables (Lower et al. 2024), and so accurate dust population predictions from high-z simulations are needed.  \nIn this work, we present a subset of cosmological zoom-in simulations of moderately massive, high-redshift ( 𝑀 ∗ ≳ 10 9 𝑀 ⊙ ; 𝑧 ≳ 5) galaxies from the Feedback in Realistic Environments (FIRE) project 4 , originally simulated in Ma et al. (2018, 2019), rerun with the integrated \'Species\' dust evolution model presented in Choban et al. (2022, C22 hereafter). This model is able to reproduce numerous galaxy-integrated and spatially-resolved observations of dust in the local universe due to its ability to track the evolution of specific dust species with set chemical compositions and incorporation of a physically motivated dust growth routine (Choban et al. 2024, C24 hereafter). Notably, this is the first application of such a dust evolution model in high-redshift simulations that resolve the multi-phase ISM and giant molecular clouds. We find that accretion is the dominant producer of dust in these massive, high-z galaxies, but their extreme burstiness and low metallicity reduce the efficiency of dust mass buildup via accretion. Compared to observations, our model produces systematically lower dust masses ( ≳ 1 dex) than inferred dust masses of 𝑧 ≳ 5 DSFGs. Even with reasonable modifications (i.e. within observational and theoretical uncertainties) to accretion timescales  \nTable 1. Parameters describing properties of simulated galaxies, all run with 7100 𝑀 ⊙ mass resolution for gas and star particles. The "Fiducial" group of simulations utilized our "Species" dust evolution model as presented in C22. The "Enhanced Stardust & Accretion" increases the SNe II dust creation efficiencies and decreases the gas-dust accretion timescale by a factor of 4. (1) Name of simulation. (2) The redshift the simulation is run to. (3) Virial mass of dark matter halo at 𝑧 final . (4) Virial radius of dark matter halo at 𝑧 final . (5) Stellar mass within virial radius at 𝑧 final . (6) The stellar half mass radius within virial radius at 𝑧 final . (7) Gas mass within virial radius at 𝑧 final . (8) Metal mass within virial radius at 𝑧 final . (9) Dust mass within virial radius at 𝑧 final . (10) Star formation rate averaged over the last 10 Myr within the virial radius.  \nand SNe II dust creation routines, we can reproduce only a subset of observed DSFGs. However, many observed dust masses lie above metal-budget upper limits from our simulations, which match the observed high-z mass-metallicity relation (Marszewski et al. 2024). Therefore, we argue that the dust masses of some DSFGs are likely overestimated due to uncertainties in assumed dust temperature. We also find that these galaxies are dominated by silicate dust and have a large spread in D/Z across ISM phases, both critical aspects that should be considered when post-processing high-z simulations.  \nThis 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 stellar, metal, and dust population properties for each galaxy in Section 3.1 and comparing them with high-z observations in Section 3.2. We discuss observational uncertainties and biases for high-z DSFGs in Section 4.1, uncertainties in our dust model in Section 4.2, and compare our findings with other high-z simulations in Sec. 4.3. Finally, we present our conclusions in Section 5.', '2 METHODOLOGY': "To study the evolution of dust in massive galaxies at high redshift, we reran a subset of cosmological simulations from the HiZ FIRE-2 suite presented in Ma et al. (2018, 2019), selecting galaxies from the suite with stellar masses 𝑀 ∗ ≳ 10 9 𝑀 ⊙ at 𝑧 = 5, 7, or 9 which matches the lowest stellar masses with measured dust masses observed at 𝑧 ≳ 5. The exact details for each simulation's final redshift ( 𝑧 final ), resulting galactic properties at said redshift, and mass resolution are provided in Table 1.", '2.1 Galaxy and Dust Evolution Models': "All simulations in this work are run with the GIZMO code base (Hopkins 2015) in the meshless finite-mass (MFM) mode with FIRE-2  \n(Hopkins et al. 2018b) model of star formation and stellar feedback. FIRE-2, an updated version of FIRE (Hopkins et al. 2014), incorporates multiple sources of stellar feedback, including stellar winds (O/B and AGB stars), ionizing photons, radiation pressure, and supernovae (both Types Ia and II). Gas cooling is followed for 𝑇 = 10 -10 10 Kincluding free-free, Compton, metal-line, 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 an uniform but redshift dependent meta-galactic background (Faucher-Giguère et al. 2009), including the effect of self-shielding. Star formation is only allowed in cold, molecular, and locally self-gravitating regions with 𝑛 H ≥ 1000 cm -3 .  \nEach 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 𝑀 ⊙ . 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 any future mention of solar abundances refers to the Andres & Gravesse abundances.  \nFIRE is ideally suited to investigate galactic dust evolution at high redshift given its success in matching a wide range of observations related to galaxies and their evolution, including the mass-metallicity relation up to 𝑧 ≳ 10 (Ma et al. 2016; Feldmann et al. 2023; Bassini et al. 2024; Marszewski et al. 2024) 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  \nFigure 1. Images of our galaxies at 𝑧 final (see Table 1) run with our 'Enhanced' dust evolution model (note the differences in scale). Each image is a mock rest-frame Hubble Space Telescope ugr composite accounting for dust attenuation utilizing the D/Z produced by our model along with assumed MW, LMC, or SMCdust opacities depending on the median gas cell Sil-to-C ratio. These images highlight the complex and varying morphologies of galaxies in our suite. The 𝑧 final ∼ 5 galaxies have varying shapes, ranging from compact to extended, due to their bursty star formation. The 𝑧 final > 5 galaxies are extremely complex, with extended structures due to their chaotic environments. Regardless of their morphology or 𝑧 final , all galaxies exhibit prominent dusty structures in the inner halo.  \n<!-- image -->  \nand giant molecular cloud (GMC) mass function (Benincasa et al. 2020). These also lead to the self-consistent development of galactic winds that eject large amounts of gas (Muratov et al. 2015; AnglésAlcá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 mass-halo mass relation.  \nOur simulations utilize the integrated 'Species' dust evolution model presented in C22, which we refer to the reader for full details. This model includes the present-day 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,  \nand astration (dust destroyed during the formation of stars). These processes are modeled self-consistently, owing to the FIRE model's in-depth treatment of the multi-phase ISM and time-resolved individual SNe events (Hopkins et al. 2018a). Notably, we restrict gas-dust accretion to cool ( 𝑇 ≤ 300 K) 5 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  \ntheoretical oxygen-bearing (O-reservoir) and nanoparticle metallic iron (Nano-iron) dust species with set chemical compositions. Consequently, each dust species has a key element 6 that limits 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. The above-listed details enable this model to match a wide range of dust observations in the local universe. In particular, the relation between galactic D/Z and metallicity along with its scatter, gas-phase element depletion trends, and varying dust population chemical compositions seen in the MW and Large and Small Magellanic clouds (C24).  \nWe also reran all of our simulations with an 'enhanced' version of our dust evolution model to investigate how uncertainties in the dust life cycle could affect our results. In particular, the 'enhanced' model increases the SNe II dust creation efficiency to 20% for all dust species and decreases the gas-dust accretion timescale by a factor of 4. We label these two model versions as Fiducial and Enhanced in Table 1 and use the \\_enh suffix to distinguish individual simulations run with the 'enhanced' version. We discuss the reasoning for the changes used in the 'enhanced' model in Sec 4.2 and showcase other tested changes to the dust model for one galaxy in Appendix A.  \nWe highlight that there are some differences between the final redshifts and resulting galactic properties presented here and those from the original simulation runs listed in Ma et al. (2019). These differences are due to changes in FIRE-2 and how final redshifts were originally chosen. Specifically, the simulations in Ma et al. (2019) were run with a version of FIRE-2 that included erroneous heating from the cosmic ray background at high redshift. This extra heating suppressed star formation in the intergalactic medium, restricting it to the densest regions in the galactic halo (Su et al. 2018; GarrisonKimmel et al. 2019). We find that the removal of this extra heating has little effect on the star formation histories of each galaxy, but star formation is generally more dispersed, occurring in less dense regions as well, but still above the threshold density for star formation. The main consequence of this dispersed star formation is a significant slowdown of the simulation, primarily for the earliest forming massive galaxies. This then affects the 'final' redshift of the simulations, which were originally chosen to be the point at which the simulations became too computationally expensive to continue. The end result is a higher 'final' redshift for simulations originally run to 𝑧 final = 7 and 9, and slight variations between reruns of the same simulation.", '2.2 Model Limitations': "Theoretical models cannot incorporate every important physical process due to resolution and computational constraints, and our models are no different. Below, we list the main physics limitations of our galaxy and dust evolution models in regards to missing/excluded physics.  \nAGN: Our galaxy formation model does not include the formation and feedback of active galactic nuclei (AGN) nor does our dust evolution model include the possible nucleation of dust grains in  \nAGNoutflows(Elvis et al. 2002). However, AGN are unlikely to form in the types of galaxies we simulate (Harikane et al. 2023; Maiolino et al. 2023; Matthee et al. 2024) and the current understanding of AGN dust nucleation is highly uncertain (Schneider & Maiolino 2024, see Sec. 6).  \nCold ISM: While our simulations resolve the multi-phase ISM, an advantage compared to previous theoretical works discussed in Sec. 4.3, the predicted structure of the cold, dense ISM, especially the dense molecular phase (Keating et al. 2020), is physics and resolution-dependent. This then affects our dust evolution model's predictions for the gas-dust accretion process, which occurs in the cold neutral and molecular media. To test our model's sensitivity to variations in ISM evolution and structure, we compare the predicted results of the 'enhanced' dust evolution model for one galaxy run with the FIRE-2 and FIRE-3 (Hopkins et al. 2023) stellar feedback and ISM physics models in Appendix A. Notably, FIRE-3 includes improved modeling of cold gas from the CNM to molecular clouds and removes the density threshold criteria for star formation. Overall we find the predicted 𝑀 dust and median D/Z are slightly lower with FIRE-3 due to a reduction in cold gas mass.  \nDust Physics: We note that all cooling and heating processes and radiative transfer modeled in our simulations are not coupled with our dust evolution model and instead follow the default assumptions in FIRE-2. Specifically, dust heating and cooling and radiative transfer assume a constant D/Z ratio, and metal-line cooling assumes no metals are locked in dust. This choice was made to avoid possible changes to galaxy evolution and resulting 𝑧 final galactic properties, which could affect our dust evolution predictions. In future work, we will investigate what effects the full integration of our dust evolution model with FIRE physics has on predicted galaxy evolution.  \nGrain Size: The primary limitation of our dust evolution model is the exclusion of evolving grain sizes. Since most dust processes depend on a dust population's effective grain size, our model's predictions will be weakest in environments where the dust population is expected to deviate significantly from an MRN (Mathis et al. 1977) grain size distribution.", '3 RESULTS': "We first showcase mock rest-frame HST ugr composite images of each galaxy to highlight the breadth of galaxy morphologies and variable dust structures contained in our simulation suite. Fig. 1 shows images of all galaxies at 𝑧 final with our 'enhanced' dust evolution model. 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 7 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. The 𝑧 final ∼ 5 galaxies vary in morphology, from compact systems with branches of infalling gas ( z5m11d and z5m12b ) to extended complexes with  \nFigure 2. Evolution of total stellar mass ( top left ), star formation averaged over 10 Myr intervals ( top right ), gas mass ( bottom left ), and median metallicity in cool ( 𝑇 < 1000 K) gas ( bottom right ) within 0.2 𝑅 vir for our simulated galaxies. For clarity, we only include one simulation for each galaxy run with the 'enhanced' dust evolution model since reruns of the same galaxy have near identical evolution. We also note that the short drops in stellar mass are due to major merger events during which the galactic center jumps between the two merging galaxies. The resulting metallicities for our galaxies are quite low compared to present-day stellar-mass analogs. These low metallicities delay the onset of efficient accretion, which is determined by a critical metallicity threshold, to later points in the galaxy's life and reduce the rate of dust growth via accretion, both critical for the buildup of dust mass. Furthermore, these galaxies are extremely burst, experiencing multiple rapid rises in their star formation which evacuate large amounts of gas, metals, and dust from the inner part of the halo. This disruption of cold gas on a galactic scale delays the build-up of dust mass in the galaxy since accretion only occurs in such environments.  \n<!-- image -->  \nprominent substructures ( z5m12a and z5m12d ). These variations largely arise from bursts in star formation that evacuate large amounts of gas from the galactic center, shutting down star formation until gas recollapses back onto the center, triggering the next burst. The 𝑧 final ∼ 7 and 9 galaxies are uniform in their heterogeneity, having complex and extended structures due to the chaotic interactions of multiple infalling systems that have yet to fully settle into one central galaxy. Regardless of these differences, all galaxies exhibit dusty structures at their centers.", '3.1 Evolution of Galactic Properties': "In Fig. 2, we present the evolution of various galactic properties, specifically the stellar mass, star formation rate averaged over 10 Myr intervals, gas mass, and median metallicity of cool ( 𝑇 < 1000 K) gas for each galaxy. These values are determined from star particles  \nand gas cells within 0 . 2 𝑅 vir of the galactic center for all galaxies 9 . For clarity, we only show simulations run with the 'enhanced' dust evolution model. While there are stochastic variations between simulations of the same galaxy, the overall evolution is very similar as can be seen in Appendix A. We highlight two galactic properties that are critical for the buildup of galactic dust mass.  \n(1) Metallicity: The galactic metallicity determines both the rate of dust buildup via accretion of gas-phase metals onto preexisting dust grains (accretion timescale inversely scales with the local metallicity) and when dust growth via accretion becomes more efficient than dust destruction via SNe (above a given 𝑍 crit ). Therefore, galaxies that reach higher metallicities earlier will experience a faster build-up of their dust mass and have more time to do so. As seen in Fig. 2, the metallicity of our galaxies is low compared to present-day stellar  \n9 We choose 0 . 2 𝑅 vir as the outer boundary for our galaxies due to the tendency of high-redshift galaxies to have more expansive stellar populations relative to their virial radii as compared to galaxies at lower redshift, where 0 . 1 𝑅 vir is commonly used. This is also consistent with other works that analyze the same simulation suite (Sun et al. 2023a,b; Marszewski et al. 2024) and we find only a small amount of the stellar mass within the virial halo is excluded by this cutoff, except for during mergers.  \nFigure 3. Evolution of the total dust mass ( top ) and median D/Z in cold neutral ( 𝑇 < 1000 K) gas ( bottom ) within 0 . 2 𝑅 vir for our simulated galaxies produced by the fiducial ( left ) and 'enhanced' ( right ) dust evolution models. Galaxy labels are the same as Fig. 2. We also highlight the 'canonical' MW D/Z = 0 . 4 ( grey line ). While both models initially predict low dust masses determined by the SNe II dust creation efficiency, the fiducial model produces ∼ 1 dex lower dust masses. Eventually, the galactic dust masses begin to increase due to the onset of efficient dust growth via accretion. However, for the fiducial model, this buildup occurs later and is slower on average due to longer accretion timescales resulting in higher critical metallicity thresholds. Ultimately, the fiducial model produces ≲ 1 dex lower dust masses when compared to the 'enhanced' model.  \n<!-- image -->  \nmass analogs 10 , with cool, dense gas reaching 𝑍 ∼ 0 . 4Z ⊙ for the most metal-rich galaxy in our sample. Examining the set of 𝑧 final ∼ 5 and 7 galaxies, z5m12b and z7m12c stand out since they reach higher metallicities at earlier times compared to their 𝑧 final counterparts, suggesting these galaxies are more likely to be dust-rich.  \n(2) Burstiness: Bursts in star formation can lead to galaxy-wide blow outs of cold gas where gas-dust accretion occurs ( 𝑇 < 300 K for our model), hindering dust growth and ultimately slowing the buildup of dust. As evident from their star formation histories in Fig. 2, all galaxies in our suite are extremely bursty, experiencing multiple, rapid increases in their SFR. The resulting feedback from these bursts evacuates large amounts of gas, metals (Muratov et al. 2015, 2017), and dust from the inner halo.  \nThe compounding effects of low metallicity and burstiness can be seen in the evolution of the galactic dust mass and dust population composition. In Fig. 3, we compare the resulting evolution of the galactic dust mass and median D/Z produced by the fiducial and 'enhanced' dust evolution model for all galaxies in our suite. In Fig. 4, we highlight the dustiest galaxies in the 𝑧 final ∼ 5 and 7 subsets, z5m12b and z7m12c , comparing a detailed breakdown of each galaxy's metal and dust population evolution produced by the two dust evolution models. We include the median metallicity, median D/Z, dust creation source contribution, and dust species composition. We also include the breakdown for each gas phase: cold neutral  \ngas ( 𝑇 < 10 3 K), warm neutral gas (10 3 K ≤ 𝑇 < 10 4 K), and warm/hot ionized gas ( 𝑇 ≥ 10 4 K). Again, only gas within 0 . 2 𝑅 vir of the galactic center is considered. We provide a similar comparison for all other galaxies in our simulation suite in Appendix B.  \nFor the fiducial model, galaxies are initially extremely dust-poor (D/Z < 0 . 01; 𝑀 dust ∼ 10 3 𝑀 ⊙ ) and entirely dominated by carbonaceous dust due to an assumed low SNe II dust creation efficiencies for all species besides carbonaceous which are taken from Zhukovska et al. (2008). Accretion eventually takes over as the dominant dust producer, growing metallic iron and later silicates, as seen by the increase in their respective species mass fractions. This buildup is relatively quick ( ∼ 0 . 2 Gyr), but the equilibrium values achieved are quite low with a maximum 𝑀 dust ∼ 3 × 10 6 𝑀 ⊙ and median D/Z ∼ 0 . 08 in cold neutral gas for our dustiest galaxy at 𝑧 final . Overall, the fiducial model struggles to build up a sizable dust mass for all galaxies in our simulation suite due to low SNe II dust creation efficiencies coupled with possibly overestimated accretion timescales, as we discuss in Sec. 4.2. In contrast, the 'enhanced' model produces larger initial dust masses (D/Z ∼ 0 . 05; 𝑀 dust ∼ 5 × 10 4 𝑀 ⊙ ) dominated by silicate dust, with appreciable amounts of carbonaceous and metallic iron dust, due to the increased SNe II dust creation efficiencies for all dust species. Similar to the fiducial model, accretion eventually takes over as the dominant producer of dust mass. However, the onset of efficient dust growth occurs at earlier times (i.e. lower 𝑍 crit ), and the build-up is faster due to the model's 4 × decrease in accretion timescales for all dust species. Metallic iron is the first dust species to grow, and then silicates, and finally carbonaceous (in the case of z5m12b ) as can be seen in the increase of the species  \nce Mass  \nSour  \nFigure 4. 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 z5m12b and z7m12c run with the fiducial ( left/middle right ) and 'enhanced' ( middle left/right ) dust evolution models. Median values are given for all gas within 0 . 2 𝑅 vir of the galactic halo ( solid ), cold neutral gas ( dashed ), warm neutral gas ( dotted ), and warm/hot ionized gas ( dash-dotted ). Note the different x-axis scales for z5m12b and z7m12c . We also highlight the spikes in Z are due to blowout events where almost the entire gas mass is ejected outside of 0 . 2 𝑅 vir , as can be seen in Fig. 2, with the majority of gas left being hot, metal-rich SNe ejecta. For both models, the galactic dust populations are initially dominated by SNe II created dust but differ in dust amount and chemical composition. In particular, the fiducial model predicts almost entirely carbonaceous dust with D/Z < 0 . 01, while the 'enhanced' model predicts predominantly silicate dust, with some carbonaceous and metallic iron, and D/Z ∼ 0 . 05. Eventually, the galactic metallicity reaches 𝑍 crit for first metallic iron and later silicates, causing accretion to be the dominant producer of dust mass and increasing the median D/Z. The rise in D/Z once accretion becomes dominant is slightly faster for the 'enhanced' model, and reaches a higher median D/Z in cold, neutral gas, D/Z ∼ 0 . 3 compared to D/Z ∼ 0 . 06 for the fiducial model. Both models also predict a large 1 dex spread in D/Z across gas phases despite only a factor of ≲ 2 difference in metallicity.  \n<!-- image -->  \n0.8  \nmass fraction or metallic iron, O Reservoir 11 , and carbonaceous, respectively. Ultimately this leads to ≲ 1 dex higher dust mass for all galaxies with a maximum 𝑀 dust ∼ 10 7 𝑀 ⊙ and median D/Z ∼ 0 . 3 in the cool neutral gas of our dustiest galaxy at 𝑧 final . Another notable result shared by both models is a large spread in D/Z, up to 1 dex, across gas phases. This is particularly relevant for works that post-process simulations with an assumed dust population which we discuss in Sec. 4.3.", '3.2 Comparisons and Predictions for 𝑧 ≳ 5 Observations': "Observations of galactic dust properties at 𝑧 ≳ 5 are limited in quantity, with only ∼ 100 direct detections of galactic rest-frame IR dust continuum. Furthermore, the dust population information these observations provide is generally limited to estimates of 𝑀 dust . Due to this limited data, the most commonly examined relation at these high redshifts is the relation between the galactic dust mass and stellar mass (dust-to-stellar mass ratio; 𝑀 dust / 𝑀 ∗ ), which is believed to measure the efficiency of dust production or the survival capability of dust grains depending on the literature (e.g. Calura et al. 2017).  \nFigure 5. Resulting relation between galactic dust mass and stellar mass within 0 . 2 𝑅 vir of the galactic center for our galaxies. We include ( top ) the exact results from our simulations and observations and ( bottom ) a schematic representation for ease of reading which shows shaded regions encompassing our simulation results and observational samples with error bars. We present results from our simulations run with the fiducial ( small colored ) and 'enhanced' ( large colored ) dust evolution model. Each connected point is -1 redshift apart starting at 𝑧 = 10, including an additional point at the simulation's final redshift, with the color denoting the redshift. We also provide two upper limits on the possible dust mass. One ( small grey ) is an upper limit for each galaxy assuming a constant D/Z = 0 . 6 for all gas within 𝑅 vir . Note the upper bound stellar masses increase slightly due to the change in radius. The other ( grey dashed ) is an upper limit determined by the total amount of metals produced by a galaxy of a given stellar mass assuming the same IMF and SNe yields as used in FIRE-2 ( 𝑀 metals ∼ 0 . 02 𝑀 ∗ ). We include observations from ALMA surveys ALPINE, from Faisst et al. (2020); Schaerer et al. (2020); Pozzi et al. (2021), ( hollow circles ) and REBELS, from Inami et al. (2022); Topping et al. (2022); Sommovigo et al. (2022a), ( hollow squares ) and group together numerous small number observations from Watson et al. (2015), Laporte et al. (2017), Strandet et al. (2017), Tamura et al. (2019), Witstok et al. (2023a), Akins et al. (2023), and one 𝑧 = 10 . 6 upper limit estimate from Fudamoto et al. (2024) ( hollow triangles ). Alternate ALPINE 𝑀 dust estimates from Sommovigo et al. (2022b) are also included ( filled circles ; ALPINE* in schematic). Both models predict a roughly constant relation with redshift with ∼ 0 . 5 dex scatter at any given stellar mass. However, the fiducial model is unable to reproduce the large dust masses observed at any stellar mass or redshift. The 'enhanced' model is able to reproduce observed dust masses for the 𝑀 ∗ ∼ 10 10 𝑀 ⊙ end of the REBELS 𝑧 ∼ 7 sample, but falls ∼ 1 dex below the ALPINE 𝑧 ∼ 5 sample. However, the 𝑀 dust from the ALPINE sample may be overpredicted, as discussed in Sec. 4.1. Moreover, given the metal content of our galaxies, only an unrealistic assumption of D/Z = 0 . 6 everywhere in the galactic halo could produce 𝑀 dust value similar to those in the ALPINE sample.  \n<!-- image -->  \nWe, therefore, first present the resulting dust-to-stellar mass ratio from our galaxies in Fig. 5, showing both the exact results from our simulations and observations and a schematic representation for ease of reading. We calculate 𝑀 dust and 𝑀 ∗ for each galaxy run with our fiducial and 'enhanced' dust evolution model starting at 𝑧 = 10 and stepping -1 in redshift up to and including 𝑧 final , with the color denoting the redshift. All values are determined using gas cells and star particles within 0 . 2 𝑅 vir of the galactic center. We further include  \ntwo extreme upper bounds on this relation. The first upper bound is determined by the available metal mass in the entire galactic halo of each galaxy. In particular, we determine a maximal dust mass by assuming a constant D/Z = 0 . 6, the maximum D/Z predicted by our model and the typical D/Z seen in the densest regions of the MW (Jenkins 2009), for all gas within 𝑅 vir of the galactic center. We also recalculate 𝑀 ∗ within 𝑅 vir for this upper bound, but this produces a relatively minor increase. However, this upper bound is dependent on  \nFigure 6. Resulting relation between galactic dust mass ( top ) and median D/Z for cool, neutral ( 𝑇 < 1000 K) gas ( bottom ) with respect to 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. Galaxy labels and redshift colors are the same as Fig. 5.  \n<!-- image -->  \nthe fraction of metals and dust ejected from the galactic halo due to the aforementioned bursty star formation. Therefore, the second upper bound is more stringent, representing all metals possibly formed in a galaxy of a given stellar mass assuming a Kroupa (2002) IMF and SNe yields used in FIRE-2 ( 𝑀 metals ∼ 0 . 02 𝑀 ∗ ).  \nComparing with observations, we include results from the ALMA surveys ALPINE (Faisst et al. 2020; Schaerer et al. 2020; Pozzi et al. 2021), and REBELS (Inami et al. 2022; Topping et al. 2022; Sommovigo et al. 2022a; Palla et al. 2024) and various small sample observations from Watson et al. (2015), Laporte et al. (2017), Strandet et al. (2017), Tamura et al. (2019), Witstok et al. (2023b), and Akins et al. (2023). Alternative estimates of 𝑀 dust for the ALPINE sample from Sommovigo et al. (2022b) are also included. These observations range in redshift from 𝑧 ∼ 4 -8 . 4 and are similarly colorcodeddependingontheir redshift. We further include an upper-bound galactic dust mass estimate for GN-z11 at 𝑧 = 10 . 6 from Fudamoto et al. (2024). We caution that the methods used for determining dust and stellar mass differ between works, leading to larger uncertainties than the error bars shown, as can be seen with the two ALPINE 𝑀 dust estimates. We discuss these caveats in detail in Sec. 4.1.  \nOverall, both the fiducial and 'enhanced' models predict a roughly power-law 𝑀 dust / 𝑀 ∗ relation ( 𝑀 dust ∝ 𝑀 𝛼 ∗ ; 𝛼 ∼ 1 . 3) for 𝑀 ∗ ≳ 10 9 𝑀 ⊙ irrespective of redshift, with, on average, a ∼ 0 . 5 dex scatter at any given stellar mass. This scatter is due to the bursty star formation of these galaxies, which periodically evacuates large amounts of gas and dust out of the central part of the galaxy. However, the 'enhanced' model predicts a relation ∼ 0 . 8 dex higher than the fiducial model and a larger scatter for 𝑀 ∗ ≲ 2 × 10 9 𝑀 ⊙ due to lower 𝑍 crit . Compared to observations, the fiducial model underpredicts 𝑀 dust  \nacross the entire observed 𝑀 ∗ range, falling ∼ 2 -0 . 5 dex below the REBELS 𝑧 ∼ 7 observations and even more so for the ALPINE 𝑧 ∼ 5 observations. The 'enhanced' model is a better match with observations, overlapping with the 𝑀 ∗ ≳ 6 × 10 9 𝑀 ⊙ subset of the REBELS, but still falls 1 dex below the ALPINE observations. This large discrepancy with the ALPINE sample is likely due to overestimated dust masses, as we discuss in Sec. 4.1, with alternative estimates from Sommovigo et al. (2022b) reducing ALPINE 𝑀 dust values by ∼ 0 . 8 dex on average. This is further bolstered by our upper bound estimates on dust mass, which still fall short of the ALPINE sample, and overlap with the small amount of 𝑀 ∗ ≳ 2 × 10 9 𝑀 ⊙ galaxies in the REBELS sample. These upper bounds are entirely unrealistic, requiring both maximal dust growth and minimal dust destruction everywhere in the galactic halo. Higher, centrally concentrated galactic metallicities in our simulations or a top-heavy IMF could alleviate this tension, but as shown by Marszewski et al. (2024), the galaxies from this suite match the observed gas-phase mass-metallicity relation up to the highest redshift where data exists, 𝑧 ∼ 10. Therefore, based on metal-budget constraints alone, the ALPINE 𝑀 dust estimates appear overestimated.  \nCurrently, no galactic metallicity measurements exist for the observed high-z DSFGs presented, which has been shown to be the main determinator of efficient dust growth. In light of JWST's ability to measure galactic metallicity at high-z, we showcase predictions for the relation between galactic dust mass and median D/Z for cool, neutral ( 𝑇 < 1000 K) gas with respect to the galactic metallicity in Fig. 6. 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. Nakajima et al. 2023; Curti et al. 2024). 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). Our simulations predict 𝑀 dust and D/Z increase with metallicity, but with differing trends. 𝑀 dust exhibits a steep rise with metallicity and a large ∼ 1 dex scatter at any given metallicity. Meanwhile, D/Z increases more slowly with metallicity with a minor amount of scatter.", '4 DISCUSSION': 'As shown in Sec. 3.2, our fiducial simulations cannot reproduce observations of extremely dusty galaxies at 𝑧 ≳ 5. Uncertainties in our dust evolution model can explain some of this discrepancy, and differing observational methodologies and biases may explain the rest. We discuss each in detail below and compare our results with other theoretical works.', '4.1.1 Sensitivity to Dust Temperature': 'Observations of dust emission SEDs at 𝑧 ≳ 5 are generally limited to one ALMA Band 6 ( ∼ 1250 𝜇 m) or 7 ( ∼ 950 𝜇 m) photometric measurement, which probes the rest-frame FIR. This typically lies at or near the Rayleigh-Jeans (RJ) tail of the dust continuum, which is well described by a modified blackbody (MBB) function (e.g. Hildebrand 1983). Assuming the FIR is optically thin, the relation for the observed flux density at frequency 𝜈 ( 𝑆 𝜈 ) is therefore  \n𝑆 𝜈 ∝ 𝜈 𝛽 𝑀 dust 𝐵 𝜈 ( 𝑇 dust ) ,  \nFigure 7. Diagram showing the relation between the observed flux density at 870 𝜇 m to inferred dust mass ratio ( 𝑆 870 𝜇 m / 𝑀 dust ) and the assumed dust temperature for set 𝛽 using the RJ approximation ( solid ) and not using the RJ approximation for 𝑧 = 4 ( dashed ) and 𝑧 = 7 ( dotted ) sources. Note the y-axis units are arbitrary. At high-z, the RJ approximation is not valid, so small changes in 𝑇 dust produce large changes in the inferred 𝑀 dust for a given 𝑆 870 𝜇 m .  \n<!-- image -->  \nwhere 𝛽 is the spectral emissivity index determined by dust population properties (grain size and composition), 𝑀 dust is the total dust mass, 𝐵 𝜈 ( 𝑇 ) is the Planck function, and 𝑇 dust is the effective dust temperature 12 . Thus, given 𝑇 dust and 𝛽 , 𝑀 dust can be determined and vice versa. Since they are difficult to derive observationally, 𝑇 dust and 𝛽 are typically set to assumed values or determined from fitting codes. All 𝑀 dust observations presented in Sec. 3.2, besides REBELS sample predictions from Sommovigo et al. (2022a) which we discuss below, utilize one of the above-mentioned methods, with assumed/fitted values in the range of 𝑇 dust ∼ 25 -60Kand 𝛽 ∼ 1 . 5 -2. These differences in assumed dust properties, specifically 𝑇 dust 13 , can have a significant impact on the inferred dust mass. In particular, the FIR SEDs probed at these redshifts are not entirely in the RJ regime (i.e. cannot approximate 𝐵 𝜈 ( 𝑇 dust ) ∝ 𝑇 dust ), and so relatively small changes to 𝑇 dust can lead to much larger changes to 𝑀 dust (Casey 2012; Cochrane et al. 2022). We highlight this sensitivity in Fig. 7, which shows the relation between the observed flux density to inferred dust mass ratio ( 𝑆 𝜈 / 𝑀 dust ) and the assumed dust temperature for sources at 𝑧 = 4 and 7 and using the RJ approximation. We set 𝜆 = 870 𝜇 m(ALMA Band 7) and assume a constant 𝛽 .  \nWe, therefore, suggest that part of the discrepancy between our simulations and observations is due to too low 𝑇 dust assumptions at these redshifts, which inflate the derived dust mass. To showcase this, we highlight differing dust mass predictions for the ALPINE sample, the sample our simulations are most in tension with as shown in Sec 3.2. Pozzi et al. (2021) made the first dust mass estimates for this sample by assuming a constant 𝑇 dust = 25 K ( 𝛽 = 1 . 8), which  \nis a common high-z assumption derived from local observations (Scoville et al. 2014, 2016, 2017). However, more recent estimates by Sommovigo et al. (2022b) predict a mean 𝑇 dust ∼ 48K (with set 𝛽 = 2 . 0). This higher 𝑇 dust corresponds to a factor of ∼ 7 reduction in the inferred 𝑀 dust compared to Pozzi et al. (2021), putting them in line with the REBELS sample and our simulations. These higher 𝑇 dust estimates are determined utilizing a novel method introduced by Sommovigo et al. (2021), which employs the [C /i.pc/i.pc] 158 𝜇 m line emission as a proxy for 𝑀 dust in order to break the degeneracy between 𝑀 dust and 𝑇 dust . However we caution that this [C /i.pc/i.pc] method relies on an assumed dust-to-gas ratio (and with it an assumed metallicity) and nontrivial [C /i.pc/i.pc]-to-neutral gas mass conversion factor ( 𝛼 C /i.pc/i.pc = 𝑀 gas / 𝐿 C /i.pc/i.pc ) 14 . Burgarella et al. (2022) also predicts higher dust temperatures and lower dust masses for the ALPINE sample ( 𝑇 dust ∼ 54 K; 𝛽 ∼ 0 . 87). Their method uses a composite/stacked IR SED, which incorporates additional MIR Herschel data fitted with a MIR power-law + MBB fit (PL-MBB; Casey 2012) to determine 𝑇 dust . However, the effective 𝑇 dust given by a PL-MBB fit is typically higher than that given by a MBB fit due to the sensitivity of the MIR, and total IR luminosity, to the small amount of hot dust near newly formed stars. In closing, the intricate details of different dust temperature definitions are beyond the scope of this work, and we point the reader to Liang et al. (2019) for an in-depth comparison of these definitions.', '4.1.2 UV-Bright Selection Bias with Detectable IR-Continuum': 'Another important aspect to consider is the inherent bias of current high-z observations towards a subset of galaxies that our suite of simulations may not include. The vast majority of observed 𝑧 ≳ 5DSFGs are UV-selected galaxies, which biases them towards less-dusty galaxies. However, these galaxies also have corresponding IR dust continuum measurements, indicating they are on the dustier side. In particular, only ∼ 20% and ∼ 40% of the original ALPINE and REBELS UV-selected samples have corresponding IR detections, respectively. Furthermore, stacking analysis utilizing non-detections for the ALPINE sample further indicates that the median dust masses for the entire sample is ≲ 0 . 5 dex lower than the direct detections (Pozzi et al. 2021). These suggest that current observations are possibly limited to a middle ground in dust content, not being the most dust-rich or the most dust-poor. However, this could also be a viewing angle effect, with dusty galaxies possibly appearing either UVbright or dark depending on their orientation (Cochrane et al. 2024). We also highlight that, due to their bursty SF, our galaxies spend relatively brief periods in the high SF regime (SFR ≳ 10 𝑀 ⊙ /yr), as can be seen in Fig. 2. This burstiness affects when these galaxies would be observable in the rest-frame UV (Sun et al. 2023a), limiting their comparability with observations. In particular, the ALPINE and REBELS samples are limited to SFR ≳ 20 𝑀 ⊙ /yr and ≳ 40 𝑀 ⊙ /yr, respectively. Ultimately, more in-depth comparisons with observable quantities will require the creation of mock SEDs using radiative transfer codes, which we leave for future work.', '4.2 Dust Model Uncertainties - Fiducial and Enhanced Model Assumptions': "While our fiducial dust evolution model reproduces some present-day galaxy-integrated and spatially resolved observations (C24), there is ∼ 10 Gyrs of evolutionary history between 𝑧 = 5 and now, which can obfuscate uncertainties in the dust life cycle. This has led to a vast array of dust evolution models which vary in their assumptions, methodology, and included physics, but all reproduce some present-day observations (e.g. Bekki 2015; McKinnon et al. 2016; Zhukovska et al. 2016; McKinnon et al. 2017; Aoyama et al. 2020; Granato et al. 2021). Therefore, it may be one or multiple uncertainties/assumptions in our model which is the cause of these low dust masses, but would not greatly affect present-day observations. In particular, two uncertainties stand out:  \n(1) Gas-dust accretion timescales are overestimated, leading to a delay in dust growth via accretion and lower D/Z. Most dust evolution models, including our own, agree that above a 'critical' galactic metallicity threshold ( 𝑍 crit ), dust growth via accretion becomes efficient and the average galactic D/Z rapidly increases (Hou et al. 2019; Li et al. 2019; Graziani et al. 2020; Parente et al. 2022). They also agree that once accretion becomes efficient, an equilibrium D/Z is eventually reached between dust growth via accretion and dust destruction via SNe shocks. However, 𝑍 crit and the equilibrium D/Z are determined by the average accretion timescale in a galaxy. If the accretion timescale is systematically overestimated, then 𝑍 crit and equilibrium D/Z will increase and decrease, respectively. The consequence of this is two-fold. First, galaxies will need to reach higher metallicities for the onset of efficient accretion, resulting in less time to build sizable dust masses by 𝑧 ≳ 5. Second, galaxies will have overall lower dust masses at a given metallicity even when accretion is the dominant producer of dust mass. While the accretion routine in our fiducial model is self-consistent and physically motivated, it is in no sense complete and makes many strong assumptions, such as neglecting dust-gas clumping factors and grain size evolution. Overall, our routine likely overestimates the gas-dust accretion timescale, as seen when comparing with local universe observations of the galactic D/Z and metallicity relation. In particular, our model predicts a rise in D/Z at ∼ 0 . 4 dex higher metallicities than observed as can be seen in Fig. 8 of C24.  \n(2) SNe II dust creation efficiency is underestimated, increasing the time it takes for accretion to build up dust mass. As previously mentioned, once accretion becomes efficient, the galactic dust mass increases until an equilibrium between dust growth and destruction is reached. This equilibrium timescale ( 𝜏 equil ; see Sec 4.2 in C24) can vary between galaxies, depending on the SNe dust destruction timescale, the fraction of the ISM mass in cold clouds where gas-dust accretion occurs, and the lifetime of those clouds. However, 𝜏 equil can be universally shortened for all galaxies if the initial amount of dust produced by stars is increased. At 𝑧 ≳ 5, SNe II are the primary producers of dust, so higher SNe dust creation efficiencies could shorten 𝜏 equil . In the case of our fiducial model, we assume extremely low SNe II dust production for all dust species besides carbonaceous dust (see Table 1 in C22) taken from Zhukovska et al. (2008). These creation efficiencies are determined by comparing in situ abundance ratios of presolar dust grains from SNe and AGBs found in meteorites, but these are limited to a handful of grains. Furthermore, both observations (e.g. Schneider & Maiolino 2024) and simulations (Kirchschlager et al. 2019, 2023, 2024) predict both extremely high and extremely low SNe II dust production.  \nTo test these uncertainties in our fiducial model, we reran simulations of z5m12a with various reasonable changes (i.e. within  \ntheoretical/observational uncertainty) to our accretion and SNe dust production routines, which are presented in Appendix A. In summary, we find that both decreasing the accretion timescale by a factor of 4 and increasing the SNe II dust creation efficiency 15 to 20% for carbonaceous, silicates, and metallic iron (set to 15%, 0.035%, and 0.1% respectively in the fiducial model) are likely needed, which is the basis for the 'enhanced' model. The decreased accretion timescale allows for efficient accretion to commence well before 𝑧 = 7 and results in a ∼ 1 dex higher equilibrium D/Z. The additional increase in the SNe creation efficiency increases D/Z by ∼ 1 dex at early times and decreases the time it takes to reach an equilibrium D/Z. As shown in Fig. 3, this results in an average ≲ 1 dex increase in the galactic dust mass for our galaxies at 𝑧 ∼ 5 -7. In the end, however, we are mainly limited by the metal budget of our galaxies. All of our galaxies have median 𝑍 ≲ 0 . 4Z ⊙ at 𝑧 final , and so even assuming the most extreme case of D/Z = 0 . 6 everywhere in the halo, we cannot reproduce a large number of 𝑀 dust observations as shown in Fig. 5. A top-heavy IMF could alleviate this issue, but given the HiZ FIRE-2 simulation suite matches the observed stellar mass-metallicity relation up to 𝑧 ∼ 10 (Ma et al. 2016; Marszewski et al. 2024), it seems unlikely that our galaxy evolution model is to blame. This gives us further credence that many high-z 𝑀 dust measurements are overpredicted, as discussed in Sec. 4.1.  \nWe highlight that the above changes would only marginally affect the local universe results presented in C24. In particular, the increased SNe II dust creation efficiencies would result in slightly higher D/Z for low-metallicity galaxies, which are dominated by SNe and AGB dust production, and shorten the predicted 𝜏 equil by < 50%. The reported 𝑍 crit 's for each dust species would also decrease by roughly the same factor as their accretion timescales. This could cause issues for C24 finding that the high 𝑍 crit for carbonaceous dust explains the lack of small carbonaceous grains in low-metallicity galaxies. However, we stress that only the accretion timescale for silicate dust needs to be shorter since silicates represent the majority of the maximum formable amount of dust.", '4.3 Previous Theoretical Works': 'Here, we compare our simulation predictions to the theoretical field at large and provide insights for assumed high-z dust populations used in simulation post-processing to produce observable quantities.', '4.3.1 Comparisons to Other High-z Dust Evolution Models': "In recent years, an increasing number of works utilizing dust evolution models integrated into semi-analytical models and galaxy simulations have focused on explaining local universe dust observations. However, relatively few of these have been turned towards 𝑧 ≳ 5 observations. The current high-z galaxy simulation landscape includes Graziani et al. (2020) and Di Cesare et al. (2023) cosmological simulations run with DUSTYGADGET, Esmerian & Gnedin (2022, 2024) post-processed simulations from the Cosmic Reionization on Computers project (Gnedin 2014), Lewis et al. (2023) DUSTiER cosmological simulation run with RAMSES-CUDATON (Ocvirk et al. 2016), and Lower et al. (2023, 2024) Cosmic Sands cosmological zoom-in simulations run with SIMBA (Davé et al. 2019). Despite the  \ndiffering galaxy formation models used, all of these works predict that gas-dust accretion is the dominant producer of dust mass for galaxies in the observed stellar mass range ( 𝑀 ∗ ≳ 10 9 𝑀 ⊙ ), similar to our findings. However, they all produce a similar dust-to-stellar mass relation, with relatively little scatter within each study, that is ≳ 0 . 5 dex higher than our 'enhanced' model predictions, and tentatively agree 16 with both low stellar mass ALPINE and most REBELS dust masses.  \nThe cause of this discrepancy with our results is likely due to the lower resolution of previous works ( ≳ 2 dex lower than our own). These simulations do not resolve the multi-phase ISM, instead relying on sub-resolution prescriptions that result in the overprediction of 𝑀 dust at high-z as we describe below. Firstly, these works employ sub-grid star formation schemes, which result in systematically less bursty star formation histories (i.e. Iyer et al. 2020). Case in point, the CosmicSandssimulations, the highest resolution simulations utilized in previous works, predict an effectively monotonic increasing SFH (Lower et al. 2023, see their Fig. 3). This lack of bursty star formation means these simulations do not experience strong blowouts of their cold gas reservoir as indicated by our simulations, and so dust growth via accretion is likely more efficient on average. Furthermore, this can suppress any predicted scatter in the stellar-to-dust mass relation. Secondly, these works employ a one-phase gas-dust accretion routine. This means accretion is not restricted to cold gas, and numerous 'sub-grid' physical processes are grouped into an overall accretion timescale normalization factor which is resolution and implementation dependent 17 . As shown in C22, when in the accretion-dominated regime, such a routine results in a relatively constant D/Z for all but the hottest gas phase within a MW-like galaxy. The combined effect of these two assumptions is an increased D/Z in all gas phases to the point where many of these simulations reach their maximum D/Z ≳ 40% for most gas in the galactic halo (Graziani et al. 2020; Di Cesare et al. 2023; Lewis et al. 2023).  \nIn contrast, works utilizing semi-analytical models have more mixed results. Notably, the DELPHI model used by Dayal et al. (2022) and Mauerhofer & Dayal (2023) uniquely predicts that accretion is subdominant to dust creation by SNe II at these epochs. Meanwhile, the Vĳayan et al. (2019) L-GALAXIES model and Triani et al. (2020) Dusty SAGE model find accretion is the dominant producer of dust mass. Despite these differences, all of these works underpredict dust masses relative to observations similar to our results. Furthermore, some works conducted a 'maximal dust model' gedankenexperiment based on the available metal budget, similar to our 𝑀 dust upper bound, and find many observations fall on or above this maximum limit (Vĳayan et al. 2019; Dayal et al. 2022).", '4.3.2 Insights for Post-Processing Observables': 'Due to the numerous observational uncertainties at high-z, as we discuss in Sec. 4.1, and upcoming observations with JWST, there has been a recent surge in post-processing existing simulations that do not include dust evolution, including the FIRE simulations, to  \npredict observables. This process utilizes radiative transfer codes coupled with assumed dust populations to create mock-SEDs from which direct observational predictions can be made (Ma et al. 2019; Liang et al. 2019; Vogelsberger et al. 2020; Shen et al. 2022; Vijayan et al. 2022; Cochrane et al. 2022, 2023, 2024). The typical dust populations assumed are quite simplistic, with the standard being a constant D/Z for all gas within the galactic halo below a set temperature above which dust is presumed to be instantly destroyed by sublimation/sputtering (typically D/Z ∼ 0 . 4; 𝑇 ≲ 10 6 K). However, as we show in Fig. 4, high-z galaxies likely have lower D/Z on average, along with a large variation between gas phases, which can have a large impact on observables.  \nFor example, current theoretical works that attempt to constrain 𝑇 dust in the early universe vary considerably due to differing dust population assumptions. Post-processed 𝑧 ≳ 2 FIRE simulations, which utilized an assumed dust population, suggest 𝑇 dust ∼ 25 ± 7 K with little evolution over redshift (Ma et al. 2019; Liang et al. 2019). Meanwhile, post-processed SIMBA simulations that use self-consistently evolved dust populations predict higher dust temperatures on average along with a large variation between galaxies ( 𝑇 dust ∼ 30 -80 K; Lower et al. 2024). However, the SIMBA simulations have > 2 dex lower resolution and do not resolve the multi-phase ISM, so this is not a one-to-one comparison. Therefore, we advise that future post-processing works utilize realistic dust populations derived from simulations like those presented here. We also plan to investigate the sensitivity of post-processed predictions in future works.', '5 CONCLUSIONS': "In this work, we investigate the evolution and buildup of dust in moderately massive galaxies at 𝑧 ≳ 5 utilizing a suite of 8 cosmological zoom-in simulations, originally simulated in Ma et al. (2018, 2019) with the FIRE-2 model (Hopkins et al. 2018b) for stellar feedback and ISM physics, rerun with the 'Species' dust evolution model (C22). This dust evolution model accounts for dust creation in stellar outflows (SNe II/Ia and AGBs), 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-bearing (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. Due to these details, this model is able to replicate a wide range of present-day observations (C24).  \nThe 8 galaxies we selected cover a range of stellar masses (10 9 𝑀 ⊙ ≤ 𝑀 ∗ ≤ 10 10 𝑀 ⊙ ) and final redshifts 𝑧 final ∼ 10 -5 similar to dusty, star-forming galaxies observed by the ALMA ALPINE and REBELS surveys. We summarize our findings on high-z dust buildup below:  \n- (i) Accretion is the dominant producer of dust mass for these galaxies. However, our fiducial model cannot reproduce observed dust masses at any redshift or stellar mass (Fig. 5) due to two factors. First, our simulations predict these galaxies have low metallicity ( 𝑍 < 0 . 5Z ⊙ ) and extremely bursty star formation (Fig. 2). These attributes delay the onset of efficient dust growth via accretion, which is determined by a critical metallicity threshold ( 𝑍 crit ), and limit the efficiency of accretion, due to the continuous disruption of cold gas where accretion occurs. Second, our model assumes low SNe II dust creation efficiencies, resulting in extremely low initial D/Z < 0 . 01.  \nThis prolongs the buildup of dust mass once accretion becomes efficient (Appendix A).  \n- (ii) Our 'enhanced' model incorporates modest changes that are within theoretical/observational uncertainty (i.e decreasing both the accretion timescale by a factor of 4 and increasing SNe II creation efficiencies to 20%). These changes result in a ∼ 1 dex increase in both the initial D/Z and the 𝑀 dust at 𝑧 final (Fig. 3 and 4). This agrees with observed 𝑀 dust values for 𝑀 ∗ > 6 × 10 9 𝑀 ⊙ from the REBELS survey but still falls ∼ 1 dex below Pozzi et al. (2021) values from the ALPINE survey. These results suggest SNe II are efficient dust producers and highlight that our accretion routine is likely missing an important physical process.  \n(iii) Given the low predicted metallicities of our simulated galaxies, the Pozzi et al. (2021) 𝑀 dust estimates for the ALPINE survey may be overestimated (Fig. 5). Post-processing with radiative transfer codes will be needed to determine the possible cause of this, but low 𝑇 dust ∼ 25 𝐾 assumptions may be to blame (Fig. 7). Alternative 𝑀 dust estimates from Sommovigo et al. (2022b) which predict 𝑇 dust ∼ 50 𝐾 also tentatively support this.  \n- (iv) Future observations to determine the burstiness of high-z galaxies are needed to constrain the efficiency of gas-dust accretion in the early Universe. Previous high-z cosmological simulations incorporating dust evolution models rely on sub-resolution prescriptions for star formation and gas-dust accretion routines. These prescriptions produce a less bursty star formation history and do not restrict accretion to cold gas, resulting in a higher accretion efficiency.\n- (v) Works utilizing radiative transfer codes to post-process simulations typically assume a dust population with a constant D/Z, but a spatially variable D/Z should be considered to understand its effects on observables. In particular, we find a maximum D/Z ∼ 0 . 3 in cool ( 𝑇 < 1000 K) gas with a ≲ 1 dex scatter across gas-phases in these high-z galaxies (Fig. 4).\n- (vi) JWST allows for metallicity measurements of these high-z DSFGs, which is the primary determinator of dust evolution. We, therefore, provide predictions for the relation between galactic dust mass and D/Z with galactic metallicity (Fig. 6).", 'ACKNOWLEDGEMENTS': "We thank Denis Burgarella for sharing details of his work and Laura Sommovigo, Guochao Sun, and Robert Feldmann for insightful discussions. DK was supported by NSF grant AST-2108324. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute. The authors acknowledge the Indiana University Pervasive Technology Institute (Stewart et al. 2017) for providing supercomputing, database, and storage resources that have contributed to the research results reported within this paper. We ran simulations using: the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562; Frontera allocations AST21010 and AST20016, supported by the NSF and TACC; Big Red 200 at the Indiana University Pervasive Technology Institute. The data used in this work were, in part, hosted on facilities supported by the Scientific Computing Core at the Flatiron Institute, a division of the Simons Foundation. This work also made use of MATPLOTLIB(Hunter 2007), NUMPY (Harris et al. 2020), SCIPY (Virtanen et al. 2020), and NASA's Astrophysics Data System.", 'DATA AVAILABILITY STATEMENT': 'The data supporting the plots within this article are available on reasonable request to the corresponding author. A public version of the GIZMO code is available at http://www.tapir.caltech. edu/~phopkins/Site/GIZMO.html .', 'REFERENCES': '- Watson D., Christensen L., Knudsen K. K., Richard J., Gallazzi A., Michałowski M. J., 2015, Nature, 519, 327\n- Witstok J., et al., 2023a, Nature, pp 1-3\n- Witstok J., Jones G. C., Maiolino R., Smit R., Schneider R., 2023b, Monthly Notices of the Royal Astronomical Society, p. stad1470\n- Zhukovska S., 2014, Astronomy and Astrophysics, 562, A76\n- Zhukovska S., Gail H.-P., Trieloff M., 2008, Astronomy and Astrophysics, 479, 453\n- ZhukovskaS., Dobbs C., Jenkins E. B., Klessen R. S., 2016, The Astrophysical Journal, 831, 147', 'APPENDIX A: EFFECTS OF DUST MODEL VARIATIONS AND ISM EVOLUTION': "To test how uncertainties in our dust evolution model affect the resulting dust population evolution and final dust mass, we reran simulations of one galaxy in our suite ( z5m12a ) with mostly reasonable variations to our dust evolution model. In particular, we tested (1) decreasing the accretion timescale by a factor of 4, (2) increasing the SNe II dust creation efficiency of each dust species to 90%, (3) same as (1) and increasing the SNe II dust creation efficiency of each dust species to 20% (labeled the 'enhanced' model in the main text), and (4) same as (3) and increasing the temperature cutoff of gas-dust accretion ( 𝑇 cut = 1000 K). We also test the robustness of our model to variations in ISM evolution by running one simulation with both the 'enhanced' dust evolution model and the FIRE-3 version (Hopkins et al. 2023) of the FIRE code. FIRE-3 makes a variety of improvements to the stellar inputs and numerical methods, focusing in particular on updating the stellar evolution tracks used for stellar feedback and nucleosynthesis with newer, more detailed models, as well as improving the detailed thermochemistry of cold atomic and molecular gas, and adopting the newer Asplund et al. (2009) protosolar reference abundances with Z ⊙ ∼ 0 . 014. The resulting galactic properties at 𝑧 final for each of these reruns are provided in Table A1.  \nFig. A1 shows the evolution of various galactic properties for each rerun, specifically the stellar mass, star formation rate average over 10 Myr, gas mass, median metallicity of cool ( 𝑇 < 1000 K) gas, dust mass, and median D/Z of cool gas. In Fig. A2 we compare a detailed breakdown of each simulation's metal and dust population evolution. Focusing first on galactic properties, it can be seen that despite stochastic variations between each rerun with FIRE-2, there is little variation in the overall galactic evolution besides the exact timing of star formation bursts. The rerun with FIRE-3 exhibits higher star formation rates at early times, similar to what is seen in Hopkins et al. (2023), resulting in a higher stellar mass and metallicity early on, but roughly matches FIRE-2 galactic properties by 𝑧 final .  \nLooking at each galaxy's dust population evolution, we can see modest variations in the dust evolution model result in considerable evolutionary changes and dust masses at 𝑧 final . Decreasing the accretion timescale results in an earlier onset of efficient accretion (i.e. lower 𝑍 crit ) and a higher overall dust mass and median D/Z. However, when this is coupled with low SNe dust creation efficiency, as with model (1), the timescale of dust buildup via accretion can be long ( ≳ 0 . 4 Gyr), which may be too long to produce 𝑧 ∼ 7 DSFGs. Increasing SNe II creation efficiency alone, and to its most extreme case as with model (2), results in an initially high median D/Z ∼ 0 . 2, but this decreases over time as more dust is destroyed by SNe shocks then can be created by SNe II and accretion. This highlights the difficulty of producing 𝑧 ∼ 5 DSFGs with primarily SNe-created dust, even when considering near maximum SNe dust creation. Modest decreases and increases to the accretion timescale and SNe II dust creation efficiency respectively, as with model (3), result in both an  \nTable A1. Parameters describing galactic properties at simulation end for reruns of z5m12a with various changes to our dust evolution model and/or with FIRE-2 vs FIRE-3 ISM physics and feedback model. Columns are the same as Table 1, with one additional column (11) Notes for each simulation on the specific changes made to the model. Overall, the final galactic properties, besides dust mass, change little between reruns with FIRE-2. FIRE-3 produces higher stellar mass and lower metal mass due to enhanced star formation at early times and updated stellar yields, which are more in line with Asplund et al. (2009) abundances.  \nearlier increase in D/Z, due to efficient accretion, and a faster increase in D/Z, due to the initially higher D/Z produced by SNe dust creation. Further increasing the temperature cutoff for gas-dust accretion, as with model (4). has little effect on the overall evolution and only slightly increases the median D/Z.  \nComparing the predicted dust evolution from FIRE-2 and FIRE-3 with model (3), FIRE-3 predicts an initially carbonaceous-dominated dust population due to its updated stellar yields (from the synthesis of Nomoto et al. 2013; Pignatari et al. 2016; Sukhbold et al. 2016; Limongi & Chieffi 2018; Prantzos et al. 2018), which have higher C and lower Si and Fe SNe II yields. Despite producing a higher metallicity at earlier times, FIRE-3 predicts a rise in D/Z via accretion at roughly the same time as FIRE-2 due to lower Si and Fe abundances at a given metallicity (Anders & Grevesse (1989) vs Asplund et al. (2009)). FIRE-3 also predicts a lower median D/Z due to an overall decrease in cold gas compared to FIRE-2. In particular, FIRE-3 predicts ∼ 2 × less gas mass with 𝑇 < 1000 K compared to FIRE-2.  \nFigure A1. Evolution of total stellar mass ( top left ), star formation averaged over 10 Myr intervals ( top right ), gas mass ( middle left ), median metallicity in cool ( 𝑇 < 1000 K) gas ( middle right ), dust mass ( bottom left ), and median D/Z of cool gas ( bottom right ) within 0.2 𝑅 vir for reruns of z5m12a with the varying modifications to our dust evolution model along with one simulation rerun with the FIRE-3 ISM physics and feedback model. We note that despite stochastic variations between FIRE-2 reruns, the resulting galactic evolution is largely unchanged. On the other hand, FIRE-3 has a higher initial SFR, resulting in a higher metallicity at early times, but produces similar galactic properties at 𝑧 final as FIRE-2. Various minor changes to our gas-dust accretion routine and SNe II dust creation efficiencies can increase the final galactic dust mass by ∼ 1 dex. Compared to FIRE-2, FIRE-3 produces a lower galactic dust mass primarily due to lower amounts of cold, dense gas due to its updated star formation criteria.  \n<!-- image -->  \n5  \n6  \n7  \nz  \n8  \n10  \n5  \n6  \n7  \nz  \n8  \n10  \n5  \n6  \n7  \nz  \n8  \n10  \n5  \n6  \n7  \nz  \n8  \nz  \n10  \n5  \n6  \n7  \nz  \n8  \n10  \nz5m12a\\_4acc A ccr etion SNe IIA GB SNe Ia 0.8 1.0 1.2 T ime [Gyr] Silicates Carbon Ir on O R eservoir SiC 0.4 8 7 6 5 10 ar iations to our dus t e v olution ∼ 1 de x higher D/Z at 𝑧 final . Ho higher D/Z, but if this is not in D/Z and a slightl y higher  \n.  \nw  \na  \nFigur  \n<!-- image -->  \n0.6  \n0.4  \n1.2  \n1.0  \n0.8  \nime [Gyr]  \nT  \n0.6  \n0.4  \n0.0  \nv  \nwith  \nz5m12a  \nof  \nuns  \nrer  \nor  \nf  \n2  \nig.  \nF  \nas  \nSame  \nA2.  \ne  \nand  \naccretion  \nefficient  \nof  \nonset  \nlier  \near  \nan  \nin  \nresults  \nalone  \nincrease  \nter  \nas  \nf  \ny  \nginall  \nmar  \na  \nin  \nresults  \naccretion  \nt  \nas-dus  \ng  \nor  \nf  \ny  \ninitiall  \nproduce  \nefficiencies  \ncreation  \nII  \ne  \nSN  \nIncreasing  \nz5m12a  \n0  \n0  \n1  \n]  \n1  \n0  \n1  \nZ  \n[  \nZ  \n2  \n0  \n1  \nK  \n3  \n0  \n1  \n<  \nall T  \nK  \n4  \n0  \n1  \n<  \nT  \nK  \n3  \n0  \n1  \nK  \n4  \n0  \n1  \nT  \n1  \n0  \n1  \nD/Z  \n2  \n0  \n1  \n0  \n0  \n1  \n1  \n0  \n1  \nF  \nraction  \nSour  \nce Mass  \n2  \n0  \n1  \n1.0  \n0.8  \n0.6  \n0.4  \nF  \nraction  \nSpecies Mass  \n0.2  \nas.  \ng  \ndense  \ncool,  \nof  \namount  \nthe  \nin  \nreduction  \ntematic  \ns  \nsy  \nDust Buildup at  \nMNRAS  \n𝑧  \n≳  \n000  \n5  \n19  \n, 1-20 (2024)  \nFigure B1. Same as Fig. 4 for z5m11d and z5m12d .  \n<!-- image -->", 'APPENDIX B: ADDITIONAL FIGURES': "This appendix shows the detailed evolution of each galaxy's metal and dust population, comparing results from the fiducial and enhanced dust evolution model similar to Fig. 4 for the rest of the galaxies in our simulation suite. Fig. B1, B2, and B3 show the evolution for z5m11d and z5m12d , z5m12a and z7m12b , and z7m12a and z9m12a respectively.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure B2. Same as Fig. 4 for z5m12a and z7m12b .  \n<!-- image -->  \nFigure B3. Same as Fig. 4 for z7m12a and z9m12a .  \n<!-- image -->"}