JWST has enabled the detection of the UV continuum of galaxies at zgt10 evidencing a population of extremely blue potentially dustfree galaxies. Interpreting the UV spectra of galaxies as they redden is complicated by the wellknown degeneracy between stellar ages dust and nebular continuum. The main goal of this paper is to develop a theoretical model for the relationship between galaxy UV slopes bursty star formation histories dust evolution and the contribution from nebular regions. We accomplish this via cosmological zoomin simulations and in specific build a layered model where we simulate the UV slopes of galaxies with increasingly complex physics. Our main results follow. i Unattenuated stellar populations with no nebular emission exhibit a diverse range of intrinsic UV slopes with values ranging from beta  3 gt 2.2 due to long delays between bursts. This is manifested by an inverse correlation between the intrinsic UV slope and sSFR for early galaxies such that higher sSFR corresponds to bluer UV slopes. ii When including dust our model galaxies demonstrate a rapid rise in dust obscuration between z  810. This increase in dust mass is due to high graingrain shattering rates and enhanced growth per unit dust mass in very small grains resulting in UVdetected galaxies at z  12 descending into ALMAdetectable galaxies by z  6. The rapid rise in dust content at z  810 leads to a systematic reddening of the UV slopes during this redshift range. iii The inclusion of nebular continuum reddens the UV slope by a median factor Delta beta  0.20.4. However when including nebular continuum our highest redshift galaxies z12 are insufficiently blue compared to observations this may imply an evolving escape fraction from HII regions with redshift.
['Astrophysics - Astrophysics of Galaxies']
The Ultraviolet Slopes of Early Universe Galaxies The Impact of Bursty Star Formation Dust and Nebular Continuum Emission
{'No Header': 'D/r.pc/a.pc/f.pc/t.pc /v.pc/e.pc/r.pc/s.pc/i.pc/o.pc/n.pc A/u.pc/g.pc/u.pc/s.pc/t.pc 27, 2024 Typeset using L A T E X twocolumn style in AASTeX631', 'The Ultraviolet Slopes of Early Universe Galaxies: The Impact of Bursty Star Formation, Dust, and Nebular Continuum Emission': "D/e.pc/s.pc/i.pc/k.pc/a.pc N/a.pc/r.pc/a.pc/y.pc/a.pc/n.pc/a.pc/n.pc, 1, 2 D/a.pc/n.pc/i.pc/e.pc/l.pc P. S/t.pc/a.pc/r.pc/k.pc, 3 S/t.pc/e.pc/v.pc/e.pc/n.pc L. F/i.pc/n.pc/k.pc/e.pc/l.pc/s.pc/t.pc/e.pc/i.pc/n.pc, 4 P/a.pc/u.pc/l.pc T/o.pc/r.pc/r.pc/e.pc/y.pc, 5 Q/i.pc L/i.pc, 6 F/e.pc/r.pc/g.pc/u.pc/s.pc C/u.pc/l.pc/l.pc/e.pc/n.pc, 7 M/i.pc/c.pc/h.pc/e.pc/a.pc/l.pc W. T/o.pc/p.pc/p.pc/i.pc/n.pc/g.pc, 3 F/e.pc/d.pc/e.pc/r.pc/i.pc/c.pc/o.pc M/a.pc/r.pc/i.pc/n.pc/a.pc/c.pc/c.pc/i.pc, 8, 9 L/a.pc/u.pc/r.pc/a.pc V. S/a.pc/l.pc/e.pc/s.pc, 10 X/u.pc/e.pc/j.pc/i.pc/a.pc/n.pc S/h.pc/e.pc/n.pc, 11 /a.pc/n.pc/d.pc M/a.pc/r.pc/k.pc V/o.pc/g.pc/e.pc/l.pc/s.pc/b.pc/e.pc/r.pc/g.pc/e.pc/r.pc 11  \n1 Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL 32611 USA  \n2 Cosmic Dawn Center at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark  \n3 Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ, 85721, USA  \n4 Department of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA  \n5 Department of Astronomy, University of Virginia, 530 McCormick Road, Charlottesville, VA 22903, USA  \n6 Max Planck Institute for Astrophysics, Garching bei Munchen, Germany  \n7 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK  \n8 Department of Physics and Astronomy 'Augusto Righi', University of Bologna, via Gobetti 93/2, 40129, Bologna, Italy  \n9 INAF, Astrophysics and Space Science Observatory Bologna, Via P. Gobetti 93/3, 40129 Bologna, Italy  \n10 Department of Physics and Astronomy, University of California, Riverside, CA, 92521, USA  \n11 Department of Physics, Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n(Received August 27, 2024; Revised; Accepted)  \nSubmitted to ApJ", 'ABSTRACT': "JWST has enabled the detection of the ultraviolet (UV) continuum of galaxies at 𝑧 > 10, evidencing a population of extremely blue, potentially dust-free galaxies. Interpreting the UV spectra of galaxies as they redden is complicated by the well-known degeneracy between stellar ages, dust reddening, and nebular continuum. The main goal of this paper is to develop a theoretical model for the relationship between galaxy UV slopes ( 𝛽 ), bursty star formation histories, dust evolution, and the contribution from nebular regions. We accomplish this via cosmological zoom-in simulations, and in specific, build a layered model where we simulate the UV slopes of galaxies with increasingly complex physics. In particular, we study the UV slopes of (i) unattenuated intrinsic stellar populations in high-redshift galaxies, (ii) reddened populations using a new on-the-fly evolving dust model, and (iii) stellar populations that include both dust and the contribution of nebular continuum. Our main results follow. Unattenuated stellar populations with no nebular emission exhibit a diverse range of intrinsic UV slopes, with values ranging from 𝛽 0 ≈ -3 → -2 . 2 due to long delays between bursts. This is manifested by an inverse correlation between the intrinsic UV slope and sSFR for early galaxies such that higher sSFR corresponds to bluer UV slopes. When including dust, our model galaxies demonstrate a rapid rise in dust obscuration between 𝑧 ≈ 8 -10. This increase in dust mass is due to high grain-grain shattering rates, and enhanced growth per unit dust mass in very small grains, resulting in UV-detected galaxies at 𝑧 ∼ 12 descending into ALMA-detectable galaxies by 𝑧 ∼ 6. The rapid rise in dust content at 𝑧 ≈ 8 -10 leads to a systematic reddening of the UV slopes during this redshift range. The inclusion of nebular continuum reddens the UV slope by a median factor Δ 𝛽 neb ≈ 0 . 2 -0 . 4. However, when including nebular continuum, our highest redshift galaxies ( 𝑧 ≈ 12) are insufficiently blue compared to observations; this may imply an evolving escape fraction from HII regions with redshift.  \nKeywords: Galaxies, Galaxy Formation, Dust formation, Interstellar Medium, JWST  \n1. INTRODUCTION  \nCorresponding author: Desika Narayanan  \ndesika.narayanan@ufl.edu  \nTherest-frame ultraviolet (UV) continuum power-law slope of galaxies has long been used as a diagnostic for the physical conditions of high-redshift galaxies (see the reviews by Dunlop 2013; Finkelstein 2016; Stark 2016; Dayal & Ferrara 2018; Robertson 2022, and references therein). Formally defined by 𝑓 𝜆 ∝ 𝜆 𝛽 (Calzetti et al. 1994), the slope 𝛽 is driven by a combination of the intrinsic ultraviolet continuum of unreddened stellar populations (Reddy et al. 2018; Calabr'o et al. 2021), potential contributions by nebular continuum (Byler et al. 2017; Cullen et al. 2017; Chisholm et al. 2022; Topping et al. 2022, 2024), and reddening by dust extinction and attenuation (e.g. Salim & Narayanan 2020, and references therein).  \nFor the vast majority of dust-enriched, star-forming galaxies, any deviation in 𝛽 from the intrinsic stellar UV continuum slope (which we refer to hereafter as 𝛽 0) is dominated by dust obscuration. Seminal work by Calzetti (1997) and Meurer et al. (1999) found that local starburst galaxies lie on a well-defined plane between the infrared excess (IRX: ratio of infrared luminosity to UV luminosity) and 𝛽 , which can be translated to a family of dust-attenuation relations (e.g. Siana et al. 2009; Popping et al. 2017; Narayanan et al. 2018; Shen et al. 2020; Liang et al. 2021). As a result, a number of studies have employed optical and near-infrared photometry to probe the rest-frame UV slopes of galaxies from 𝑧 ≈ 2 -9 in order to investigate the potential evolution of dust attenuation properties of galaxies during the first ∼ 3 -4 billion years of the Universe's history. The general sense from these surveys is that galaxies toward the Epoch of Reionization (EoR; 𝑧 > 6) are systematically bluer in their UV slopes than galaxies at 𝑧 ∼ 2, with the inference that dust attenuation becomes less impactful at earlier times (e.g. Stanway et al. 2005; Bouwens et al. 2010; Finkelstein et al. 2010; McLure et al. 2011; Bouwens et al. 2012; Dunlop et al. 2012; Finkelstein et al. 2012; Bouwens et al. 2014). At a fixed redshift, some studies additionally find bluer UV slopes with decreasing UV luminosity (Bouwens et al. 2012; Cullen et al. 2023; Tacchella et al. 2022), which is natural in a scenario where more massive galaxies /one.sup have elevated dust content (e.g. B'ethermin et al. 2015; da Cunha et al. 2015; Donevski et al. 2020). The inferred rise in dust content from rest-frame UV continuum observations from 𝑧 = 8 → 2 is complemented by depletion measurements (P'eroux & Howk 2020; P'eroux et al. 2023), extinction measurements (M'enard et al. 2010; M'enard & Fukugita 2012), SED fits to galaxies (Dunne et al. 2011; Beeston et al. 2018; Driver et al. 2018; Pozzi et al. 2020), and theoretical models (Aoyama et al. 2018; Li et al. 2018), all of which point towards a rise in the cosmic dust density from  \nearly times toward cosmic noon ( 𝑧 ∼ 2) (P'eroux & Howk 2020).  \nWhile some limited samples of UV detections of even higher redshift ( 𝑧 > 10) galaxies have been studied in the preJWST era (e.g. Ellis et al. 2013; Oesch et al. 2016; Wilkins et al. 2016; Tacchella et al. 2022), the sheer number of restframe UV detections at 𝑧 > 10 has exploded since the successful launch of JWST/NIRCam (e.g., Naidu et al. 2022; Curtis-Lake et al. 2023; Donnan et al. 2023a,b; Finkelstein et al. 2023a,b; Harikane et al. 2023; Robertson et al. 2023; Casey et al. 2024; Ciesla et al. 2024; Roberts-Borsani et al. 2024). This has enabled the investigation of the redshift evolution of the UV slope 𝛽 to increasingly early times, where the impact of dust obscuration may become relatively minimal. To wit, recent observations by Cullen et al. (2024) Topping et al. (2024), and Morales et al. (2024b) found average UV slopes of value ⟨ 𝛽 ⟩ ≈ -2 . 5 to -2 . 8 at redshifts 𝑧 > 10, which is consistent with their expectations of dust-free stellar populations within the first 500 Myr given the specifics of their population synthesis modeling.  \nAt the same time, Atacama Large Millimetre Array (ALMA) observations of galaxies at 𝑧 ∼ 6 have demonstrated significant dust reservoirs by 𝑧 = 6 (e.g. Pozzi et al. 2021; Bouwens et al. 2022; Inami et al. 2022), with dust to stellar mass ratios ∼ 10 -3 , already comparable to lower-redshift systems (Donevski et al. 2020). If the JWST UV-detected 𝑧 > 10 galaxies are related to the ALMA REBELS 𝑧 ≈ 6 galaxies (and indeed, in this paper, we demonstrate that such a connection may be credible /two.sup ), then this would imply a sudden onset of dust obscuration in the relatively short time between 𝑧 = 10 → 6. How galaxies go from inferred dust-free stellar populations to significantly dust-enriched in only ∼ 0 . 5 Gyr represents a theoretical challenge (e.g. Ferrara et al. 2023; Ziparo et al. 2023). The observed inference of the rise of dust obscuration is further complicated by the potential impact of intrinsically red UV slopes due to aged stellar populations (e.g. Topping et al. 2022; Cullen et al. 2024), as well as nebular continuum emission from free-free and free-bound processes (Byler et al. 2017; Katz et al. 2024).  \n/two.sup A very rough back-of-the-envelope calculation additionally suggests that a causal connection between the JWST UV-detected galaxies at 𝑧 > 10 and ALMA-detected galaxies at 𝑧 ∼ 6 may be reasonable. The lower bound of galaxies with dust-to-stellar ratios of 10 -3 at 𝑧 ∼ 6 have stellar masses ∼ 10 9 M ⊙ (Topping et al. 2022, see also Figure 3 in this paper as well). At the same time, some studies have inferred stellar masses ∼ 10 8 M ⊙ for 𝑧 ∼ 9 -10 UV selected galaxies (Morales et al. 2024b, acknowledging, of course, the challenges in inferring stellar masses in potentially bursty systems at such high redshifts; Tacchella et al. 2022; Whitler et al. 2023; Ciesla et al. 2024; Narayanan et al. 2024). This implies a stellar growth of ∼ 10 9 M ⊙ in ∼ 500 Myr, which is achievable with a SFR ∼ 2M ⊙ yr -1 . Main sequence star formation for these stellar masses is already SFR ∼ 1 -10M ⊙ yr -1 at these redshifts (Speagle et al. 2014; Lower et al. 2023).  \nThemaingoalofthis paper is to develop a theoretical model for the relationship between galaxy UV slopes, dust growth, bursty star formation histories, emission from nebular regions, and cosmological galaxy growth in the early Universe. Weaccomplish this via a combination of cosmological zoomin simulations (in order to generate a set of model galaxies to study), a new model for the formation, growth, and evolution of multi-species and multi-sized dust, and dust radiative transfer calculations (in order to generate mock observations). In § 2, we describe these methods in detail. In § 3, we take the reader on a tour of the physical properties of our model galaxies to help place them in a broader context. We follow this in § 4 with a study of the UV properties of intrinsic, unreddened stellar populations: this is important as aging stellar populations redden, resulting in degenerate solutions with the impact of dust and nebular continuum. We follow this in § 5 with an exploration of the onset of dust obscuration in the early Universe. In § 6, we investigate how the emission from nebular regions impacts the UV slopes of our model galaxies. In § 7 we provide discussion, and in § 8 we summarize.  \nFinally, as discussed previously in this section, there are a range of potential physical effects on 𝛽 that we will investigate in this paper. In an attempt to avoid confusion, we summarize our working definitions of the various forms of 𝛽 in Table 1. For all mock SEDs, we fit our SEDs at the center of the Calzetti et al. (1994) wavelength windows.", '2.1. Galaxy Formation Models': 'We run a suite of cosmological zoom-in simulations of massive galaxies down to 𝑧 = 6 using the /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc galaxy formation model within the moving mesh /a.pc/r.pc/e.pc/p.pc/o.pc hydrodynamic and gravity code (Springel 2010; Weinberger et al. 2020). The /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc physics model is described fully in Marinacci et al. (2019), and we refer the reader to this paper for an in-depth discussion of the model; we summarize the salient points here.  \nPrimordial cooling of the gas occurs via two-body collisional processes, free-free emission, recombination, and Compton cooling off of cosmic microwave background photons (Katz et al. 1996). Once enriched with metals, gas can additionally undergo metal line cooling with rates computed as a function of temperature and density via /c.pc/l.pc/o.pc/u.pc/d.pc/y.pc photoionization calculations (Ferland et al. 2013) as described in Vogelsberger et al. (2013). Low-temperature gas cools via a range of processes, including metal-line cooling, fine-structure emission, and molecular cooling via fits to the Hopkins et al. (2018) /c.pc/l.pc/o.pc/u.pc/d.pc/y.pc cooling tables (Marinacci et al. 2019). Selfshielding occurs in gas with densities 𝑛 > 10 -3 cm -3 via the Rahmati et al. (2013) redshift-dependent parameterization. Cosmic ray heating follows the density-dependent prescription of Guo & White (2008), and photoelectric heating follows  \nthe Wolfire et al. (2003) density, metallicity, and temperaturedependent rates.  \nStar formation occurs in gravitationally bound molecular gas (Hopkins et al. 2013) above a density threshold 𝑛 thresh = 100 cm -3 , following a volumetric Kennicutt (1998) relation / 𝑀 ∗ = 𝜖 𝑀 gas / 𝑡 ff, where / 𝑀 ∗ is the star formation rate (SFR), 𝑀 gas is the gas mass, and 𝑡 ff is the free fall time. Following Marinacci et al. (2019), we set the star formation efficiency 𝜖 = 0 . 01 (Krumholz & Tan 2007; Krumholz & Dekel 2010). The molecular fraction of the interstellar medium (ISM) is computed via the Krumholz et al. (2008) prescription linking the H2 fraction to the local gas surface density and metallicity.  \nOnce formed, stars return energy to the nearby ISM via a range of feedback processes. We compute the fraction of supernovae (SNe) from each formed star particle assuming a Chabrier (2001) initial mass function, with an additional delay time distribution for deriving the number of Type Ia supernovae events (Vogelsberger et al. 2013; Torrey et al. 2014). The details for the SNe-driven stellar mass loss rates, as well as the coupling of the energy and momentum from SNe ejecta to the ISM are provided in Marinacci et al. (2019) and Zhang et al. (2024). Similarly, feedback is included via radiation from young stars via photoionization, radiation pressure, and OB and AGB stellar winds. Taken together, these feedback processes regulate the star formation rates in our model zoomin galaxies, as well as the physical conditions in the ISM. One major consequence of the default /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc galaxy formation framework that we adopt is bursty star formation histories (SFHs), which may be reasonable for early Universe galaxies per observational constraints (Ciesla et al. 2023; Dressler et al. 2023; Endsley et al. 2023; Looser et al. 2023; Shen et al. 2023). This said, it is worth noting that the magnitude of this burstiness may be sensitive to the details of unresolved SNe feedback physics (Zhang et al. 2024).  \nIn addition to injecting feedback, aging stellar populations return mass and metals to the ISM. The enrichment of the ISM is tracked for 9 individual elements, as well as metal content as a whole. Metals that are present in the ISM are advected as passive tracers with the fluid flow. This approach allows for self-consistent co-evolution of galaxies and their metal content including, e.g., the emergence of a time and mass dependent mass-metallicity relation and internal metallicity gradients. The locally resolved free abundance of metals, in turn, can impact/modulate the accretion rates of metals onto dust grains.', '2.2. Dust Physics': 'In order to investigate the buildup of dust reservoirs in early galaxies, and their obscuring impact on stellar populations, we require a model for dust formation, growth, and destruction in galaxies. In particular, to properly model the local extinction of UV photons from individual star particles in our  \nTable 1. 𝛽 Definitions Used in this PaperTable 2. Physical properties of model zoom-in simulations considered in this paper.  \nsimulations (as well as the global attenuation), we need to know the local dust mass, grain size distribution, and chemical composition. To do this, we have implemented the newly developed dust model of Li et al. (2021) and Narayanan et al. (2023) on the fly in cosmological zoom-in simulations, which build on the algorithms developed by McKinnon et al. (2016, 2017, 2018). We refer the reader to Narayanan et al. (2023) for details of the coupling of our dust model with the /a.pc/r.pc/e.pc/p.pc/o.pc /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc galaxy formation model, and summarize the salient points here.  \nDust is produced through the condensation of metals that are ejected from evolved stars. The dust yields are taken from Schneider et al. (2014) for AGB stars, and Nozawa et al. (2010) for SNe dust production. This dust is initialized with a lognormal size distribution, though dust particles lose their memory of these initial conditions rapidly owing to coagulation and shattering processes (more on this below). Each dust particle is discretized into 16 logarithmic size bins in the range 4 × 10 -4 ≤ ( 𝑎 / 𝜇 m ) ≤ 1. The mass of dust produced by evolved stars follows the methodology of Dwek (1998), with condensation efficiencies taken from Ferrarotti & Gail (2006) for AGB stars, and Bianchi & Schneider (2007) for SNe. We assume a silicate-carbonaceous chemical breakdown of our dust particles, where the total carbon mass in a dust particle corresponds to the carbonaceous dust mass, and the remainder to silicate.  \nDust grains evolve from their initialized state by growing mass and size via metal accretion and coagulation, and being destroyed via shattering events and shocks in star-forming regions. The growth timescale is inversely proportional to  \nthe metal density and the square root of the temperature:  \n𝜏 accr = 𝜏 ref GLYPH<18> 𝑎 0 . 1 𝜇 m GLYPH<19> GLYPH<18> 1000 cm -3 × 𝑚 H × 𝑍 ⊙ 𝜌 𝑍 GLYPH<19> GLYPH<18> 10K 𝑇 𝑔 GLYPH<19> 1 2 GLYPH<18> 0 . 3 𝑆 GLYPH<19> (1)  \nwhere 𝜌 𝑍 is the metal density, 𝑇 𝑔 is the gas temperature, and 𝑆 is the sticking coefficient. We adopt ( 𝜏 ref , 𝑍 ⊙ , SI) = GLYPH<0> 0 . 224 Gyr , 7 × 10 -3 GLYPH<1> for silicates (assuming a composition of MgFeSiO4 for silicates; Weingartner & Draine 2001), and ( 𝜏 ref , 𝑍 ⊙ , C) = GLYPH<0> 0 . 175 Gyr , 2 . 4 × 10 -2 GLYPH<1> for carbonaceous grains. We additionally adopt temperature dependent sticking coefficient following Zhukovska et al. (2016) which drops at higher temperatures. This has the effect of reducing the growth rates in the warm and dense ISM that is heated by stellar feedback (Choban et al. 2024). Following Choban et al. (2022), we assume that the growth timescale is limited by the least abundant element required by the grain species. Building upon the framework that we introduced in Narayanan et al. (2023), we have further updated our dust model to also include size-dependent growth rates following Weingartner & Draine (1999): metals are preferentially depleted onto very small dust grains owing to enhanced Coulomb potentials. Functionally, we employ the ISM-phase dependent enhancement vs grain size curves from Weingartner & Draine (1999), with a log-linear interpolation between the temperature-dependent growth enhancement rates. Finally, as will be described in a forthcoming paper (Narayanan et al., in prep), our cosmological simulations do not resolve the densities comparable to dense clumps in giant clouds, where dust growth may be most efficient. As a result, we assume that gas is clumped on scales below our resolution limit. In practice, we assume that 50% of the mass of a cell has a clumping factor 𝐶 = 10, which impacts all density-dependent rates in our dust model.  \nWeconsider the impact of grain-grain collisions on the size distribution of dust grains. There are two important effects: dust shattering in high-speed collisions, which transforms large grains into smaller grains, and dust coagulation, which is the vice-versa. We follow McKinnon et al. (2018) and Li et al. (2021) in modeling the transformation of grain sizes via collisional encounters, and refer to Narayanan et al. (2023) for the relevant species-dependent equations. We employ the default threshold velocities for collisional encounters from Narayanan et al. (2023), noting that these have resulted in both model extinction laws in Milky Way analogs comparable to the Cardelli et al. (1989) observational constraints (Li et al. 2021), as well as polycyclic aromatic hydrocarbon (PAH) mass fractions comparable to constraints from the Local Group (Smith et al. 2007; Narayanan et al. 2023).  \nDust can be destroyed via thermal sputtering in the ISM, and via SNe shocks in star-forming regions. The former process is driven by erosion by hot electrons, where the sputtering timescale is adopted from the analytic approximation of Tsai & Mathews (1995), and is linearly proportional to the grain radius, and inversely proportional to the gas density and temperature. SNe shocks destroy grains via a sputtering process as well. Here, we model the evolution of the grain size distribution near SNe following the models of dust destruction in blastwaves by Nozawa et al. (2006) and Asano et al. (2013) (see Narayanan et al. 2023, for the relevant equations).', '2.3. Zoom-In Technique': 'With the aforementioned /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc+dust physics enabled in the /a.pc/r.pc/e.pc/p.pc/o.pc moving mesh hydrodynamics and gravity solver, we proceed to simulate the evolution of 5 galaxies using the well-established cosmological zoom-in technique (see Vogelsberger et al. 2020, for a review of this technique). We first run a ( 50 / ℎ ) 3 cMpc 3 dark matter only simulation with 256 3 dark matter particles initialized at 𝑧 = 99 down to 𝑧 = 0. The initial conditions are computed with /m.pc/u.pc/s.pc/i.pc/c.pc (Hahn & Abel 2011). We then select halos of interest in the 𝑧 = 0 snapshot, and construct an ellipsoidal mask around all particles within 2 . 5 × the radius of the maximum distance dark-matter particle in the halo with the /c.pc/a.pc/e.pc/s.pc/a.pc/r.pc galaxy analysis package (Thompson et al. 2014). The region encapsulated by this ellipsoidal mask is the Lagrangian high-resolution region that will be re-simulated at higher resolution, and with baryon physics included. We split the particles within the mask until the effective particle count in the zoom-in region is 2048 3 for a baryon (star and gas) mass resolution in the high-resolution region of 3 × 10 5 𝑀 ⊙ . Dust particle masses are a factor ∼ 100 lower mass on average. These resolution limits are set based on convergence studies for the main physical and mock observable properties that we model here. We then re-simulate these halos from the initial conditions ( 𝑧 = 99) down to 𝑧 = 6.  \nWe summarize the main physical properties of our simulated galaxies in § 3 and Table 2.', '2.4. Mock Observations': 'Wegenerate mock UV SEDs from our model galaxies with the publicly available /p.pc/o.pc/w.pc/d.pc/e.pc/r.pc/d.pc/a.pc/y.pc dust radiative transfer package (Narayanan et al. 2021). Functionally, /p.pc/o.pc/w.pc/d.pc/e.pc/r.pc/d.pc/a.pc/y.pc wraps /y.pc/t.pc, /f.pc/s.pc/p.pc/s.pc and /h.pc/y.pc/p.pc/e.pc/r.pc/i.pc/o.pc/n.pc for grid generation, stellar population synthesis, and Monte Carlo radiative transfer, respectively (Turk et al. 2011; Conroy et al. 2010; Robitaille 2011). For each galaxy snapshot, we compute the stellar SED for every star particle based on its age and metallicity with /f.pc/s.pc/p.pc/s.pc. We assume /m.pc/i.pc/s.pc/t.pc stellar isochrones, and a stellar IMF consistent with the hydrodynamic simulations.  \nThe stellar light is emitted isotropically, and is absorbed, scattered, and re-emitted by dust in the galaxy. Here, the radiative transfer occurs on a Voronoi mesh built around the dust particles with the active dust model. Informed by the grain size distribution in every cell, we compute the local extinction law on a cell-by-cell basis following Li et al. (2021). In detail, the wavelength-dependent optical depth seen by every stellar photon can be computed in terms of an extinction efficiency:  \n𝜏 ( 𝑎, 𝜆 ) = ∫ LOS 𝜋𝑎 2 𝑄 ext ( 𝑎, 𝜆 ) 𝑛 𝑑 ( 𝑎 ) d 𝑎 d 𝑠 (2)  \nwhere 𝑄 ext is the extinction efficiency, d 𝑠 is the path length of the ray, and 𝑛 d ( 𝑎 ) is the number density of grains with sizes [ 𝑎, 𝑎 + da ] . We assume extinction efficiencies from Draine & Lee (1984) and Laor & Draine (1993) for silicate and carbonaceous grains, respectively.  \nThedirection and frequency of photons are randomly drawn from uniform distributions, and the photons are propagated through the dusty ISM until they either escape the galaxy, or reach a randomly drawn optical depth drawn from an exponential distribution. This process is iterated until the equilibrium dust temperature has converged (in practice, we iterate until the energy absorbed by 99% of the cells has changed by less than 1%.) At this point, the aggregate SED from the model galaxy is computed via ray tracing. When reporting 𝛽 values, we report the viewing-angle averaged values, noting that there may be some intrinsic dispersion owing to sightline as the galaxies become dustier (e.g. Lovell et al. 2022; Cochrane et al. 2024).', '2.5. Nebular Emission Model': 'For a subset of models, we investigate the impact of nebular continuum emission from star-forming regions. To do this, we implement the subresolution model of Garg et al. (2022, 2023), and refer the reader to these works for details. We summarize the important points here.  \nFor all star particles with ages 𝑇 age < 10 Myr, we compute the shape of the stellar SED from /f.pc/s.pc/p.pc/s.pc, using the stellar age  \nFigure 2. Model star formation histories for galaxies listed in Table 2 From top to bottom, the SFHs are for galaxies h10,h15,h17,h20 &h25 respectively.  \n<!-- image -->  \nFigure 1. Observable and physical evolution of model galaxies withcosmictime. Top: 𝑀 UV vs 𝑧 for model galaxies (circles), colorcoded by their total dust mass. In comparison are the observational results from Austin et al. (2023); Cullen et al. (2024); Tacchella et al. (2022) and Topping et al. (2024). The model galaxy UV luminosities are similar to moderate-area JWST surveys (e.g. CEERS; Finkelstein et al. 2023a), though brighter than deep fields (e.g. NGDEEP; Bagley et al. 2024). In general, our models are comparable to those from the aforementioned observational samples at 𝑧 > 10. This suggests that our model galaxies may serve as reasonable analogs for those detected at 𝑧 > 10 with HST and JWST, which is the primary galaxy sample of interest for this paper. Bottom: Stellar mass evolution as a function of redshift, color-coded by dust mass.  \n<!-- image -->  \nand metallicity as before. We assume this SED is incident on a spherical HII region, and use /c.pc/l.pc/o.pc/u.pc/d.pc/y.pc /v.pc17.00 (Ferland et al. 2017) to compute the nebular emission properties. We fix the inner radius of the HII region as 𝑅 inner = 10 17 cm, following the motivation presented in Garg et al. (2023). The hydrogen density /three.sup of the nebular region is assumed 𝑛 HII = 10 2 cm -2 , and we model the region as dust-free (in order to focus on the impact of diffuse dust on 𝛽 in this paper), and the escape fraction is assumed 𝑓 esc = 0. Using this setup, we compute the emergent nebular spectrum from the HII region, and then replace the stellar SED with this SED for propagation through the dusty ISM with our /p.pc/o.pc/w.pc/d.pc/e.pc/r.pc/d.pc/a.pc/y.pc calculations.', '3. THE PHYSICAL AND OBSERVABLE PROPERTIES OUR MODEL EARLY UNIVERSE GALAXIES': 'In this section, we aim to orient the reader around the main physical and (mock) observable properties of our model  \ngalaxies. We summarize their basic physical properties in Table 2.  \nWebegin our tour of our model galaxies by examining their selectability in current JWST surveys. We present results below 𝑧 < 12 as the dust radiative transfer is not well-resolved at higher redshifts, owing to the relatively few simulation dust particles (which of course represent many real dust particles) constructing the Voronoi mesh. In the top panel Figure 1, we compare the UV absolute magnitude of our model galaxies (computed at a rest-frame wavelength of 𝜆 = 1500 ˚ A) to the observed samples of Tacchella et al. (2022); Austin et al. (2023); Cullen et al. (2024); Morales et al. (2024b) and Topping et al. (2024). The simulated points are color-coded by their total dust mass. It is clear from the top panel of Figure 1 that the general overlap between the observed data and our simulated galaxy mock photometry is reasonable at redshifts 𝑧 > 10, with typical 𝑀 UV ≈ -18 →-20. At lower redshifts ( 𝑧 < 10), our galaxies are brighter than those selected by Topping et al. (2024) and Morales et al. (2024b), owing to the fact that we study individual galaxies in evolution, which get brighter with cosmic time at these redshifts, and the relatively small volumes probed by those surveys.  \nIn the bottom panel of Figure 1, we show the stellar masses of our model galaxies as a function of redshift. In the redshift interval of interest ( 𝑧 ≈ 10 -12), the stellar masses range from 𝑀 ∗ ≈ 10 7 -10 8 𝑀 ⊙ , comparable to inferences by Morales et al. (2024b). In Figure 2, we show the model star formation histories for our 5 model galaxies down to 𝑧 = 6. The /s.pc/m.pc/u.pc/g.pc/g.pc/l.pc/e.pc stellar feedback model drives a bursty star formation history as is seen in other explicit feedback models (e.g. /f.pc/i.pc/r.pc/e.pc; Sparre et al. 2017).  \nFigure 3. Comparison of 𝑀 dust vs 𝑀 ∗ for our model galaxies with observations demonstrates both reasonable dust properties in our model galaxies, as well as an evolutionary connection between 𝑧 > 10 UV-bright galaxies, and 𝑧 = 6 massive, infraredluminous systems (c.f. Figure 1). The circles are our simulation results at 𝑧 > 6, color-coded by their redshift, while the red triangles with errorbars are observational data taken from a range of individual detections and surveys (Cooray et al. 2014; Watson et al. 2015; Knudsen et al. 2017; Laporte et al. 2017; Marrone et al. 2018; Hashimoto et al. 2019; Bakx et al. 2021; Fudamoto et al. 2021; Endsley et al. 2022; Dayal et al. 2022; Ferrara et al. 2022; Topping et al. 2022; Sommovigo et al. 2022).  \n<!-- image -->  \nThe end result for our model galaxies, which may be detectable in the rest-frame UV at 𝑧 = 10 -12 with JWST, are massive, and dust-rich systems at 𝑧 = 6. The individual galaxy stellar masses of our model galaxies rise from a typical value of 𝑀 ∗ ≈ 10 7 𝑀 ⊙ at 𝑧 = 10 to 𝑀 ∗ ≈ 10 9 𝑀 ⊙ at 𝑧 = 6. Concomitant is the rise in dust mass, buoyed first by production processes, followed a sharp rise in mass owing to metal accretion in the ISM. The rapid growth of these galaxies renders them comparable to those detected in major ALMA surveys. As an example, Figure 3 shows the relationship between stellar mass and total dust content of our model galaxies at 𝑧 = 6 as compared to dust and stellar mass inferences from a range of high𝑧 observations (Cooray et al. 2014; Watson et al. 2015; Knudsen et al. 2017; Laporte et al. 2017; Marrone et al. 2018; Hashimoto et al. 2019; Bakx et al. 2021; Fudamoto et al. 2021; Endsley et al. 2022; Dayal et al. 2022; Ferrara et al. 2022; Topping et al. 2022; Sommovigo et al. 2022). It is evident that as galaxies in our modeled mass range approach 𝑧 ≲ 8, their dust-to-stellar mass ratios rise significantly (which is a topic we will return to in § 5). Figure 3 not only highlights model viability by demonstrating a reasonable correspondence between the late time ( 𝑧 = 6) descendants of our 𝑧 = 10 galaxy analogs, but also demon-  \nbright systems detected by JWST at 𝑧 > 10, and massive, infraredluminous systems at 𝑧 = 6.', '4. THE ROLE OF BURSTY STAR FORMATION ON INTRINSIC ULTRAVIOLET COLORS': 'Wenowturn to our analysis of the UV slope 𝛽 in early Universe galaxies. In this section, our key goal is to understand the potential range of UV slopes of unreddened galaxies, and in specific, we demonstrate significant variability in the UV slopes of unreddened galaxies owing entirely to bursty stellar populations. Periods of ongoing star formation naturally have relatively blue UV slopes, whereas periods of star formation dormancy have redder UV slopes.  \nWeshow this point explicitly in Figure 4, where we plot the UVslope ( 𝛽 ) for unattenuated /four.sup stellar populations vs redshift for an arbitrary model galaxy in our simulation sample (h25; we choose this galaxy because its complicated SFH provides an instructive example, though show the equivalent plot for all model galaxies in the appendix). We defer discussion of the evolution of 𝛽 with redshift when dust reddening is included to § 5.  \nGenerally, the unattenuated UV slopes are at their bluest during periods of intense star formation, and then redden during periods of low star formation rate. This is seen explicitly with the pink line in Figure 4, which corresponds to the right hand axis, quantifying the galaxy star formation rate. As the galaxy star formation history progresses, however, the ability for a starburst to blue-en /five.sup the UV slope is diminished, owing to the contribution of the underlying older stellar population to the UV colors. This results in an inverse relationship between the UV slope and the specific star formation rate (sSFR ≡ SFR/ 𝑀 ∗ ) of a galaxy. We quantify this in Figure 5, where we compare 𝛽 0 to the galaxy sSFR for all galaxy snapshots in our simulation campaign. Evident in Figure 5 are both a trend with 𝛽 0 and sSFR, as well as trends at constant sSFR (i.e., vertical slices) where 𝛽 0 becomes redder at later times, reflecting the buildup of an older stellar population.  \nThe upshot from Figures 4 and 5 is that significant reddening of UV slopes in the early Universe can occur, even prior to the onset of significant dust obscuration. In this particular model, we find a range from roughly 𝛽 ≈ -3 → -2 . 2 is attainable even without the inclusion of dust or nebular continuum . Similar effects were noted by both Cullen et al. (2024) and Topping et al. (2024) in stellar population synthesis models. We further note that a consequence of using  \nFigure 4. UV slopes from unattenuated stellar populations can show a diverse range of values, owing to bursty star formation histories in low-mass galaxies in the early Universe. Unreddened 𝛽 0 (blue dots) vs 𝑧 for example galaxy h25. Aging stellar populations redden the UV slope from 𝑧 = 12 → 9 in this example galaxy, until bursts of star formation act to blue-en 𝛽 . We see a range of 𝛽 ≈ [-3 , -2 . 2 ] from intrinsic stellar populations alone (not including any contribution from dust or nebular continuum). This said, as the stellar mass of a galaxy increases, increasingly large bursts of star formation are necessary to blue-en 𝛽 , tamping down the impact of star formation below 𝑧 = 7 in this example (see also Figure 5).  \n<!-- image -->  \nan explicit feedback model with bursty star formation histories is that the resulting galaxies are expected to experience substantial dynamic swings in their age-driven reddening. In contrast, galaxies modeled with smoother star formation histories - such as those that may be expected from cosmological simulations with an artificial pressurization of the ISM (e.g. Narayanan et al. 2024) - will experience less rapid/dramatic changes in the age distribution of their stellar populations, and therefore less rapid/dynamic changes to the age-driven reddening.', '5. THE RISE OF DUST OBSCURATION AT ULTRA-HIGH 𝑍': 'Havingestablished the typical dynamic range in 𝛽 0 owing to bursty star formation histories ( § 4), we now turn our attention to the impact of dust obscuration on the UV slope (i.e., 𝛽 dust). In Figure 6, we show in the top panel the UV slope 𝛽 dust vs 𝑧 for all of our model galaxies, and in the bottom panel  \nΔ 𝛽 dust = 𝛽 dust -𝛽 0. Δ 𝛽 dust quantifies the direct impact of dust reddening on the intrinsic stellar UV slope, 𝛽 0. We arbitrarily define Δ 𝛽 dust < 0 . 05 as the region where dust has relatively minimal impact (this is ∼ 5% of the dynamic range seen in 𝛽 from stellar populations alone), and denote this region in the bottom half of Figure 6.  \nGenerally, at 𝑧 ≈ 8 -9, the galaxy population that we model transitions from having relatively minimal impact of dust on the UV slope to more of an impact, with approximately half of the galaxies transitioning above Δ 𝛽 dust > 0 . 05 by 𝑧 = 8. By 𝑧 = 6 -8, dust begins to impact a substantial fraction of galaxy UV slopes. This corresponds to an overall trend in 𝛽 dust to rise with decreasing redshift, with comparable UV slopes as the observational comparison sets of Cullen et al. (2024) and Topping et al. (2024). In general, between 𝑧 = 12 → 6, we see a typical rise from our dust-free solution of 𝛽 0 ≈ -3 to as red as 𝛽 dust ≈ -2. This overall rise in 𝛽 with 𝑧 is not necessarily monotonic, however, owing to bursts  \nFigure 5. Relationship of galaxy unreddened UV slope from stellar continuum alone ( 𝛽 0 ) vs specific star formation rate (sSFR ≡ SFR/ 𝑀 ∗ ) for all model galaxies. Generally, aging stellar populations tend to redden UV slopes, while bursts of star formation drive bluer UV slopes (c.f. Figure 4). This said, as the stellar mass of a galaxy increases, increasingly large bursts of star formation are necessary to blue-en the UV slope, resulting in an anti-correlation between 𝛽 0 and galaxy sSFR.  \n<!-- image -->  \nof star formation that can somewhat blue-en 𝛽 dust (c.f. § 4). We note that the impact of these bursts on blue-ening the UV SED may be mitigated if nebular continuum emission plays a significant role in the UV SEDs, and we discuss this further in § 6.  \nThe increase in dust reddening at 𝑧 ≈ 8 -9 owes to a sudden increase in the dust masses at this time. In the top panel of Figure 7, we plot the dust mass evolution for all of our model galaxies (this information was otherwise available in Figure 3 as well, and is just distilled slightly differently here). At 𝑧 ≈ 10, we see rapid growth in galaxy dust masses for the mass range considered in our simulation sample, with an increase by a factor 100 -1000 by 𝑧 ≈ 6. This allows the relatively blue galaxies at high𝑧 to transition into dusty star-forming galaxies with dust masses comparable to ALMAdetections at 𝑧 ≈ 6 (Figures 1 and 3). This transition to an epoch of dust mass growth owes to the production of small grains in grain-grain shattering events enabled by a turbulent ISM associated with star formation.  \nIn the bottom panel of Figure 7, we show the average velocity between any two simulation dust particles for every galaxy modeled in our sample, with error bars denoting the 1 𝜎 dispersion in each galaxy snapshot. The pink shaded region shows the imposed velocity thresholds (that are speciesdependent) for grain collisions to result in a shattering event: collisions above this threshold velocity result in a grain-grain shattering, and the increased production of small grains. As  \nthe galaxy star formation rates rise with time, the average interstellar velocity dispersion increases, and subsequently so does the grain-grain shattering rate. The surface areas per mass are larger for small dust grains than large ones, and the growth rates via metal accretion are faster /six.sup (c.f. Equation 1). This is quantified even further in Figure 8, where we show the normalized dust grain size distribution evolution with redshift for example galaxy h25; the size distributions are normalized at 𝑎 = 0 . 06 𝜇 m, the canonical dividing line between small and large grains (though this choice is arbitrary; Hirashita 2015). It is clear that the small-to-large ratio increases dramatically over the redshift range considered. As a result, the epoch of significant dust growth is associated with the epoch of significant rises in galaxy star formation histories for a given galaxy. These growth rates eventually plateau as the dust spatial distribution becomes more diffuse as the galaxy grows, and the metal densities decrease. For the range of masses modeled here, the epoch of peak dust growth during the redshift range considered is approximately 𝑧 ≈ 8 -10.', '6. ON THE ROLE OF NEBULAR CONTINUUM ON THE UV SLOPES OF EARLY GALAXIES': 'HII regions surrounding massive stars can redden the UV continuumowingtofree-free and free-bound processes (Byler et al. 2017; Katz et al. 2024). Indeed, Stanway et al. (2016); Topping et al. (2022, 2024) and Cullen et al. (2024) all find that intrinsic stellar populations with nebular continuum (but no dust) can achieve minimum colors of 𝛽 ≈ -2 . 6. In this section, we investigate the impact of nebular continuum on our model galaxies.  \nIn the top panel of Figure 9, we show a histogram of Δ 𝛽 neb (i.e., isolating the impact of nebular reddening over that of dust) for all of our model galaxies, across redshift. The distribution is broad, but peaks in a region Δ 𝛽 neb ≈ 0 . 2 -0 . 4, quantifying the typical reddening we expect in our models from nebular continuum emission.  \nThis said, we argue that the contribution of nebular continuum to the UV continuum slope cannot be uniform across redshift. In the bottom panel of Figure 9, we show the evolution of 𝛽 neb (i.e. 𝛽 including dust and nebular continuum) with redshift, and as a reference, in the light colored blue circles, show 𝛽 dust (i.e. the UV slope suffering from dust reddening alone). At the earliest times ( 𝑧 > 10), the model with nebular continuum cannot produce UV slopes sufficiently blue to match the extremely blue ( 𝛽 < -2 . 7) observed slopes. At the same time, the model with only dust obscuration can result in sufficiently blue colors at 𝑧 > 11 to compare well against  \nFigure 6. UV slope vs redshift for all model galaxies where we include the impact of dust (but not nebular continuum; 𝛽 dust ). The top panel shows this UV slope, where the red-purple points are our model results color-coded by the total dust mass, while the stars are observations from Austin et al. (2023); Cullen et al. (2024); Morales et al. (2024b); Tacchella et al. (2022) and Topping et al. (2024). The thin grey line denotes the bluest colors seen in our dust-free models. The bottom panel shows the difference in the intrinsic (unreddened) UV slope from the UV slope that includes dust ( Δ 𝛽 dust ), quantifying the impact of dust on the UV slope. Generally, the increase in dust obscuration at lower redshifts drives redder UV 𝛽 slopes.  \n<!-- image -->  \nobservational constraints, owing to the relative lack of dust and the dominance of intrinsic stellar populations in the UV SED. One possible inference from this result is an increasing escape fraction of ionizing photons with increasing redshift. We expand on this possibility further in § 7.', '7.1. Comparison with Observational Studies': "JWST has allowed for a significant number of detections of galaxies at ever-higher redshifts, allowing for a direct comparison between our model results and observations at 𝑧 > 10. While individual studies do not (yet) converge on the magnitude of the trend of 𝛽 dust with redshift, in aggregate there are clear observational trends demonstrating the reddening of the UV slopes of galaxies between 𝑧 = 12 → 6. Of course any such trend would be expected to be mass dependent: more massive halos collapse earlier, begin star formation earlier, and begin to form dust earlier. This said, because our model  \nsample of galaxies overlaps with the rough UV luminosity range of those observed (Figure 1), as well as at least some stellar mass constraints (Morales et al. 2024b), our modeled trend of 𝛽 dust vs 𝑧 appears to be in reasonable correspondence with that observed (Figure 4; also Tacchella et al. 2022; Austin et al. 2023; Cullen et al. 2024; Morales et al. 2024b; Topping et al. 2024). As observations push to deeper limits (or, for example, probe lower intrinsic masses via the aid of gravitational lensing), one may expect bluer UV slopes to persist even below 𝑧 < 10.  \nWhat is less clear, however, is whether the correspondence between our model 𝛽 dust and 𝑧 and observed galaxies owes in reality to the buildup of galactic dust, as we have suggested. Wehavepainted a circumstantial picture that links the buildup of dust in the current 𝑀 UV ≈ -18 -20 JWST detections at 𝑧 > 10 to ALMA detectable galaxies at 𝑧 ≈ 6, though in reality there is little direct observational evidence for the emergence of cosmic dust in these galaxy populations over this redshift  \nFigure 7. Evolution of dust mass and average grain velocity dispersions as a function of redshift. Top: We show the dust mass for all of our model galaxies as a function of redshift. Bottom: This rapid dust growth at 𝑧 ≈ 10 is due to the production of small dust grains via grain-grain shattering events. These events increase with the velocity of grain collisions. Here, we plot the median relative velocity between two dust grains for each model galaxy in our simulation sample. The error bars show the 1 𝜎 dispersion in these velocity dispersions, and the pink shaded region shows the range of velocity (species-dependent) thresholds employed for shattering to occur in grain-grain collisions.  \n<!-- image -->  \nrange. The stellar age-reddening degeneracy prevents this from being fully deducible from UV continuum observations alone (c.f. § 4). Put directly, in our models, we would attribute at least some of the reddening observed in 𝛽 through 𝑧 = 9 by Cullen et al. (2024) and Topping et al. (2024) as being due to stellar aging. In order to distinguish between this scenario vs early dust buildup, we advocate for two potential observational routes to directly test our bespoke assertions. The first is the pursuit of thermal infrared emission in 𝑧 > 10 galaxies, which may require the assistance of gravitational lensing given current instrument sensitivities. The second is via the inference of cosmic dust via the 2175 ˚ 𝐴 absorption feature (e.g. Witstok et al. 2023; Markov et al. 2024). This UV'bump' potentially owes to ultrasmall carbonaceous dust grains, and may therefore betray the emergence of cosmic dust in galaxies. If this feature is related to PAHs, and if PAHs form in a top-down shattering formation scenario as is suggested by at least some models (e.g. Narayanan et al. 2023), then the inference of the 2175 ˚ 𝐴 bumpmayadditionally provide evidence for a rapid growth scenario that owes to grain-grain shattering, and the associated increased growth rates on small dust grains in the turbulent ISM of high𝑧 galaxies (e.g. Figure 7).  \nFigure 8. The relative fraction of small dust grains increases with decreasing redshift, resulting in rapid dust growth (c.f. Figure 7). Shown is the normalized grain size distribution for model galaxy h25 (as an example), color-coded by redshift. The size distributions are normalized at 𝑎 ≈ 0 . 06 𝜇 m, the canonical dividing line between 'small' and 'large' grains in the literature. Shattering events drive increased fractions of small grains, which enhances accretion of metals due to the large surface area-to-mass ratio of small dust grains.  \n<!-- image -->  \nWenowturn our discussion to the implications of very blue UV SEDs ( 𝛽 dust ≈ -3), in light of the potential reddening effects of nebular continuum from HII regions around massive stars. As we demonstrated in Figure 9, a reasonable model for the nebular continuum contribution to the reddening of UV slopes at high𝑧 results in the complete lack of ability to produce the very blue galaxies observed by JWST at 𝑧 > 10 (e.g. Cullen et al. 2024; Topping et al. 2024). As a reminder, in our model, we assume 𝑓 esc = 0 for all HII regions. Therefore, one possible implication of Figure 9 is an evolution of 𝑓 esc with redshift, such that the emission star-forming regions in galaxies at increasingly high redshift are more easily able to escape their HII regions. Indeed, there may be some observational evidence for such a heuristic model. Topping et al. (2024) performed /c.pc/l.pc/o.pc/u.pc/d.pc/y.pc photoionization modeling of HII regions, and demonstrated that large escape fractions may be necessary to reproduce UV colors as blue as has been observed by both their group, as well as Cullen et al. (2024). This is consistent with the weak emission line spectra observed by Topping et al. (2024) in their sample. Beyond this, Chisholm et al. (2022) demonstrated a strong correlation between 𝛽 and 𝑓 esc for low-redshift ( 𝑧 ≈ 0 . 3) galaxies. By applying the Chisholm et al. (2022) 𝛽 -𝑓 esc relation to their observed sample, Cullen et al. (2024) infer a potential rise of up to a factor ∼ 3 in 𝑓 esc for galaxies binned below and above 𝑧 = 10 . 5 (though see Mu˜noz et al. 2024, which suggests lower 𝑓 esc values may be  \nFigure 9. Nebular Continuum drives a Reddening of the UV slope of 𝛽 ≈ 0 . 2 -0 . 4 , but results in galaxies that are too red at very high redshifts ( 𝑧 > 10 ). Top: Histogram of Δ 𝛽 neb (i.e., the difference in UV slope including dust and nebular continuum with that of just including dust) for all of our model galaxies. Bottom: Evolution of 𝛽 neb (i.e., including dust and nebular continuum) with redshift (colored points), compared to the UV slope just including dust ( 𝛽 dust : empty blue circles). When including nebular continuum, the colors are too red to match the highest𝑧 observed data ( 𝑧 > 10). We expand on possible physical explanations for this discrepancy in § 7.  \n<!-- image -->  \nnecessary to match ionizing photon count constraints). This said, we note that it is possible that observed galaxies with very blue UV slopes may be undergoing small bursts of star formation that drive very blue colors. As demonstrated via population synthesis modeling by Topping et al. (2024), 5 Myr after a burst, the nebular continuum emission is weakened sufficiently that for a short period the UV colors can be extremely blue ( 𝛽 ≈ -2 . 8) even for 𝑓 esc = 0 models. This effect may not be captured in our models, given the cadence of our snapshot output.", '7.2. Relationship with Existing Theoretical Models': "The observation of very blue UV-bright galaxies discovered by JWST has generated a flurry of theoretical activity in recent years. Much of the theoretical effort has framed the problem in terms of the relative non-evolution of the UV luminosity function above 𝑧 > 7 (e.g. Naidu et al. 2022; Bouwens et al. 2023; Donnan et al. 2023a; Finkelstein et al. 2023a; Donnan et al. 2023b; Harikane et al. 2023; P'erez-Gonz'alez et al. 2023). Here, one natural solution is the reduction of dust attenuation at 𝑧 > 7, to offset the natural decline in massive halos with increasing redshift. Ferrara et al. (2023) developed an analytic model for the observed galaxy abundances and UV luminosity functions at high-redshift. An important  \naspect to this model, as well as others that aim to model dust at very high redshift (e.g. Zhao & Furlanetto 2024) is that they are constrained by the boundary conditions of relatively little dust attenuation at 𝑧 > 10, as well as the inferred presence of significant dust reservoirs from ALMA detections of massive galaxies at 𝑧 ∼ 6. Ferrara et al. (2023) demonstrate that a decrease in dust attenuation at 𝑧 > 11 can satisfy observational constraints on the UV luminosity function, which is a coupled concept with producing a population of extremely blue galaxies. These authors are agnostic as to the origin of the lack of dust attenuation, positing possibilities such as dust removal due to early winds, offset star-dust geometries (Ziparo et al. 2023), as well as low dust to stellar mass ratios. Our interpretation, based on the cosmological hydrodynamic zoom-in simulations performed here, supports the latter scenario.  \nThe scenarios proposed above are fundamentally different: in two of the physical models (early dust ejection, as well as complex star-dust geometries), the premise is that the dust is already in place at early times, but is simply low opacity. This is in contrast to a model in which galaxies - comparable to the mass currently being detected by JWST at 𝑧 > 10 grow significant dust reservoirs during this transition epoch. Our model clearly advocates for the latter scenario. By and large, simulations that aim to explicitly model the evolution of dust masses during the Epoch of Reionization agree that significant dust growth occurs in galaxies of order the mass range modeled here between 𝑧 > 10 and 𝑧 ≈ 6 (Graziani et al. 2020; Dayal et al. 2022; Lower et al. 2023; Esmerian & Gnedin 2023; Lewis et al. 2023a; Choban et al. 2024). This said, models differ in detail regarding the origin of this dust (i.e., growth-dominated vs production dominated), the magnitude of dust growth, and when the transition from production dominated to growth dominated occurs. Sorting out the origin of these differences is complicated as there is a degeneracy between differences in the underlying algorithms for dust evolution and the ISM physical conditions in varying hydrodynamic models, and would benefit from targeted numerical experiments and code comparison studies comparable to the /a.pc/g.pc/o.pc/r.pc/a.pc simulation suite (Kim et al. 2014).  \nA handful of studies aimed at probing 𝛽 in the context of galaxy simulations and dust find similar levels of reddening in the UV SED owing to the contribution of dust. Esmerian & Gnedin (2023) and Lewis et al. (2023b) implemented a single grain-size dust model into ART and RAMSES, respectively, and found a similar range of 𝛽 values reddening owing to dust at 𝑧 < 9. Katz et al. (2023) find even redder colors when considering 𝛽 neb, though implement a more simplified dust model with an SMC-type extinction curve.  \nFinally, on the front of the contribution of nebular continuum to 𝛽 , Katz et al. (2024) presented a recent series of numerical experiments, studying the impact of nebular continuum on the UV spectra of early Universe galaxies. Katz  \n- et al. 2024 demonstrate a reddening of Δ 𝛽 neb ≈ 0 . 6 (going from 𝛽 0 ≈ -3 . 1 →-2 . 5 in the most extreme case), similar in magnitude to our more simplified models for nebular emission. Even when considering the impact of IMF variations and the assumed density of nebular regions, Katz et al. (2024) finds a minimum slope of 𝛽 neb ≈ -2 . 7, just slightly bluer than our study.", '8. SUMMARY': 'In this paper, we have combined a suite of cosmological zoom-in simulations designed to represent massive galaxies in the early Universe with dust radiative transfer modeling to study the impact of bursty star formation histories, the production and growth of dust, and emission from nebular regions on the rest-frame UV continuum slope ( 𝛽 ) in galaxies. We designed our study as a series of numerical experiments, by adding more physics in the radiative transfer modeling successively, and studying the results. In particular, we simulated the UV slopes of (i) unreddened populations within simulated high-redshift galaxies, (ii) reddened populations by adding the impact of dust (in radiative transfer), and (iii) reddened populations with the contribution of nebular continuum.  \nOur main results follow.  \n- 1. We simulate galaxies with 𝑧 = 0 parent halo masses ∼ 2 -6 × 10 13 𝑀 ⊙ , which results in galaxies with 𝑧 = 6 stellar masses 𝑀 ∗ ≈ 0 . 3 -2 × 10 10 𝑀 ⊙ . Galaxies of this mass have UV luminosities similar to those being detected in JWST surveys at 𝑧 > 10 (Figure 1). We therefore use these galaxies to study the impact of bursty star formation histories, dust, and nebular continuum on 𝛽 :\n- 2. Unreddened stellar populations exhibit a diverse range of intrinsic UV slopes, with values ranging from a dust-free 𝛽 0 ≈ -3 to values as red as 𝛽 0 ≈ -2 . 2 (Figure 4). The somewhat red values of 𝛽 0, even without the reddening effects of dust or nebular continuum, owe to long delays between bursts in bursty star formation histories (Figure 2). Bursts of star formation counter-act this reddening, though are less effective as the underlying older stellar population grows. As a result, there is an inverse correlation between the intrinsic UV slope and sSFR for early galaxies (Figure 5).\n- 3. When considering the effect of dust , the UV colors galaxies become successively redder between 𝑧 = 12 → 6 (Figure 6), due to the rapid growth of dust mass during this period (Figure 7). This growth in dust is due to grain-grain shattering in a turbulent ISM, and enhanced dust growth rates on small dust grains due to  \ntheir larger surface areas (Figures 7 and 8). 𝑧 ≈ 8 -9 is roughly where galaxies transition from being unaffected by dust to affected by dust in the UV (Figure 6), though this is true only for galaxies within the mass range that we simulate. This suggests that the reddening of UV colors seen in recent JWST observations during this epoch owes to the buildup of the first major dust reservoirs in these galaxies.  \n- 4. The inclusion of nebular continuum reddens the UV slope by a median factor Δ 𝛽 neb ≈ 0 . 2 -0 . 4. This said, whenincluding nebular continuum, our highest redshift galaxies ( 𝑧 ≈ 12) are insufficiently blue compared to observations, and we cannot achieve the dust free-like blue UV colors of 𝛽 ≈ -3 as is observed (Cullen et al. 2024; Topping et al. 2024). This may imply an evolving escape fraction from HII regions with redshift, such that the highest redshift sources have significant escape fractions.', '9. ACKNOWLEDGEMENTS': "Desika, Dan and Steve would all like to express appreciation to Adriano Fontana, Paola Santini, Jim Dunlop, Davide Elbaz, Alice Shapley, and Rachel Somerville for having organized 'The Growth of Galaxies in the Early Universe - IV' at the Sexten Center for Astrophysics in 2024. This project idea came from talks at this conference, followed by conversations had during a night out in Venice on the way home. Desika would like to additionally thank Sidney Lower for having left behind enough legacy code from her PhD thesis at the University of Florida in an easily digestible enough format that the analysis for this paper was fairly painless. Desika is grateful to the Institute for Astronomy at the University of Edinburgh, and the Cosmic Dawn Centre at the Niels Bohr Institute for graciously hosting him while he worked on significant parts of this paper, and thanks Gabe Brammer, Romeel Dav'e and Harley Katz for helpful conversations while preparing this manuscript. DN and PT were funded by NASA ATP grant 80NSSC22K0716, and are grateful to NASA for trusting us enough to fund what probably seemed like a fairly risky science idea of modeling active dust in galaxies. DN is additionally grateful to the Aspen Center for Physics, which is supported by NSF grant PHY-1607611, which is where the original framework for /p.pc/o.pc/w.pc/d.pc/e.pc/r.pc/d.pc/a.pc/y.pc was envisioned. F. Cullen acknowledges support from a UKRI Frontier Research Guarantee Grant (PI Cullen; grant reference EP/X021025/1). FM acknowledges funding by the European Union - NextGenerationEU, in the framework of the HPC project - 'National Centre for HPC, Big Data and Quan- tum Computing' (PNRR - M4C2 - I1.4 - CN00000013 - CUP J33C22001170001).", 'N/a.pc/r.pc/a.pc/y.pc/a.pc/n.pc/a.pc/n.pc /e.pc/t.pc /a.pc/l.pc.': 'Figure 11. Unreddened UV slope ( 𝛽 0 ) vs 𝑧 for galaxy h15. See Figure 4, and the associated text, for details.  \n<!-- image -->  \nFigure 12. Unreddened UV slope ( 𝛽 0 ) vs 𝑧 for galaxy h17. See Figure 4, and the associated text, for details.  \n<!-- image -->  \nFigure 13. Unreddened UV slope ( 𝛽 0 ) vs 𝑧 for galaxy h20. See Figure 4, and the associated text, for details.  \n<!-- image -->', 'APPENDIX': 'Here we show the relationship between 𝛽 0, redshift, and SFR for all galaxies in our model sample (aside from h25, which is already shown in Figure 4).'}