The A
The first galaxies formed 200-300 million years after the big bang, hosting many of the first stars and seeding the creation of every galaxy in the Universe today (V. Bromm & N. Yoshida 2011). The smallest of these early galaxies formed stars for 1-2 Gyr before being quenched by cosmic reionization and stellar feedback (T. M. Brown et al. 2014), so those that survived until now are composed of ancient stars from 13 Gyr ago. Stars from these ultrafaint dwarf galaxies (L < 10 5 L e ) are therefore relics from the era of the first stars and galaxies, preserving signatures of early chemical enrichment (M. S. Bovill & M. Ricotti 2009;J. D. Simon 2019).
Early galaxies played pivotal roles in the cosmic story. Over time, they merged to form the galaxies we observe today, including the Milky Way. Stars stripped from these dwarf galaxies likely compose the metal-poor component of the Milky Way's stellar halo (A. Frebel & J. E. Norris 2015). These galaxies are also tied to reionization both as victims of quenching and as important sources of ionizing photons (T. M. Brown et al. 2014;J. H. Wise et al. 2014). Due to their age and limited star formation histories, they offer invaluable insights into the earliest stages of star formation and chemical enrichment (e.g., A. P. Ji et al. 2015). In the past 20 yr, dozens of these ultrafaint dwarf galaxies have been discovered in the Local Group (e.g., B. Willman et al. 2005; K. Bechtol et al. 2015;J. D. Simon 2019). High-resolution spectroscopy has provided chemical abundances of over 15 of these dwarfs, with many more to come in the next decade. Simultaneously, chemical abundances of millions of stars in the Milky Way stellar halo are being obtained through wide-field spectroscopic programs such as GALAH (G. M. De Silva et al. 2015), H3 (C. Conroy et al. 2019), APOGEE (S. R. Majewski et al. 2017), RAVE (M. Steinmetz et al. 2006), SEGUE (B. Yanny et al. 2009), DESI (C. Allende Prieto et al. 2020), and LAMOST (X.-Q. Cui et al. 2012), with more to come from WEAVE (G. Dalton et al. 2014), 4MOST (R. S. de Jong et al. 2019), SDSS-V (J. A. Kollmeier et al. 2017), and the Rubin Observatory (Ž. Ivezić et al. 2019).
A key pathway to understanding early star and galaxy formation is through the study of stellar chemical abundances in ancient stars, particularly those found in ultrafaint dwarfs and the Milky Way's stellar halo (e.g., K. Brauer et al. 2019). These chemical abundances serve as "fingerprints," preserving crucial information about prior Population III (Pop III) and Population II (Pop II) star formation, hierarchical galaxy assembly, and nucleosynthetic processes. This is especially significant for stars that formed in the smallest, earliest galaxies, as their kinematic signatures are rapidly lost during mergers (e.g., K. Brauer et al. 2022), but their chemical compositions remain detectable. However, the complexity of observed abundance patterns in dwarf galaxies has surpassed the capabilities of many current theoretical models of chemical evolution (e.g., A. P. Ji et al. 2020), many of which still rely on the one-zone models of B. M. Tinsley (1980).
Chemical evolution models that rely on simple parametric approaches and assume homogeneous mixing are capable of reproducing mean trends but incapable of modeling scatter (e.g., B. H. Andrews et al. 2017;B. Côté & C. Ritter 2018). When chemical abundance observations of dwarf galaxies were sparse, these models were sufficient to explain observed trends and differences between dwarf galaxies as a function of luminosity (e.g., E. Tolstoy et al. 2009; E. N. Kirby et al. 2011). Recent observations, however, have revealed abundance scatter in many elements both within and between dwarf galaxies (e.g., V. Hill et al. 2019;A. P. Ji et al. 2020;J. Mead et al. 2024). Such variations from mean trends imply the presence of complex galaxy formation processes like sourceand time-dependent metal mixing, hierarchical galaxy merging, and bursty star formation (e.g., A. Emerick et al. 2020).
Over the past few decades, significant work has been done studying the chemodynamical evolution of galaxies using cosmological hydrodynamics simulations (e.g., B. D. Oppenheimer & R. Davé 2008; R. P. C. Wiersma et al. 2009; J. Schaye et al. 2010Schaye et al. , 2015;;S. Shen et al. 2010; C. M. Simpson et al. 2013;P. F. Hopkins et al. 2014;O. Agertz et al. 2021;F. Renaud et al. 2021;S. Buder et al. 2024). These simulations, coupled with increased attention to feedback processes, have made substantial strides in reproducing global galaxy trends, such as the evolution of the mass-metallicity relationship (e.g., A. Obreja et al. 2014;X. Ma et al. 2016; R. Davé et al. 2017), and more detailed quantities like metallicity distribution functions and the evolution of individual species abundances (A. Marcolini et al. 2008;Y. Revaz et al. 2009;T. Sawala et al. 2010; Y. Revaz & P. Jablonka 2012; Y. Hirai & T. R. Saitoh 2017; M. Jeon et al. 2017). Due to simple chemical enrichment models and resolution, however, these simulations are all unable to explain the observed chemical abundance scatter in Milky Way dwarf galaxies.
More recently, hydrodynamic galaxy simulations have begun to reach parsec resolution and include star particles representing individual stars (e.g., A. Emerick et al. 2019;N. Lahén et al. 2020;T. A. Gutcke et al. 2021;Y. Hirai et al. 2021;F. Calura et al. 2022;J. M. Hislop et al. 2022;E. P. Andersson et al. 2023; U. P. Steinwandel et al. 2023;Y. Deng et al. 2024). This level of detail allows simulations to capture variations in when, where, and how individual stars inject energy and metals. Such detail in a simulation allows investigation into the origins of the full chemical abundance distributions observed in dwarf galaxies beyond just the mean elemental trends. Y. Revaz et al. (2016) showed that below mass resolution of 10 3 M e stellar populations are not representative of the full initial mass function (IMF). This leads to oversimplifications in enrichment and the inability to properly model abundance scatter. This was similarly argued by M. C. Smith (2021) for mass resolution below 500 M e , finding that stochastic variations in stellar populations become important when modeling radiation and galactic wind driving. In the smallest galaxies, this becomes especially important, as each individual feedback event affects the evolution of the galaxy. To understand the origin of these distributions and interpret both present and future chemical abundance data, simulations need to effectively combine extremely high resolution star particles and gas mixing with detailed element yields and metal tracing techniques.
We are developing a series of hydrodynamic simulations to model early galaxy formation in a cosmological context, called the AEOS project, that focus on tracing individual stars throughout their life cycles and modeling their chemical evolution in unprecedented detail. In this paper, we present our initial simulation. This is a 1 cMpc cubical domain containing a variety of early galaxies. In subsequent simulations, we will perform zoom-ins of small early galaxies that can be regarded as early analogs of the surviving ultrafaint dwarf galaxies. We use the methodology of A. Emerick et al. (2019, hereafter E19), with updated models and stellar element yields described in Section 2. These simulations aim to provide insights into the following:
1. metal mixing within the interstellar medium (ISM) at high redshift; 2. the origins of the spreads in the observed stellar chemical abundances; 3. how physically motivated prescriptions for inhomogeneous mixing and abundance scatter can be developed; 4. the identification of signatures associated with stars formed from Population III enrichment; 5. the role of α-and s-process elements in constraining star formation timescales and stellar ages; and 6. the impact of adopting a star-by-star approach to star formation in simulations.
These simulations aim to shed light on the intricate interplay between galaxy formation processes, early nucleosynthetic events, and the stellar chemical abundances we can observe today. This paper introduces the simulation and the novel handling of star-by-star chemical enrichment. In Section 2, we describe the methods implemented in the simulations, focusing on the new or updated physics models. We present an overview of the fiducial simulation in Section 3 and the metal evolution analysis in Section 4. We conclude in Section 5.
Our simulations use the adaptive mesh refinement technique implemented in the code ENZO, a grid-based hybrid code that solves the equations of hydrodynamics and gravity (G. L. Bryan et al. 2014;C. Brummel-Smith et al. 2019). We achieve a physical resolution of 1 pc in regions of high gas density, allowing us to resolve the small-scale structure of the ISM. The computational grid solves for hydrodynamic quantities, such as gravitational potential from gas, stars, and dark matter, and evolves star particles with an adaptive particlemesh N-body solver (Section 2.1.1). Radiative cooling is modeled with a nine-species chemical network (Section 2.1.2), and we implement stellar feedback through mass and energy injection from stellar winds and supernovae (SNe; Section 2.3). We model star formation for Population III and Population II stars independently (Sections 2.2.1 and 2.2.2, respectively), implement their radiation feedback (Section 2.3.1), and model both stellar winds and SNe from each population (Sections 2.3.2 and 2.3.3,respectively). We capture chemical evolution driven by nucleosynthesis from different yield channels, including Population III CCSNe, Population II CCSNe, AGB winds, massive stellar winds, and Type Ia SNe (SNe Ia; Section 2.4). We describe the initial conditions in Section 2.5.
These methods expand on those used in E19 with the addition of a Population III star formation and stellar feedback model, which is an updated version of the one used in J. H. Wise et al. (2012b), and a new, carefully chosen set of stellar yields.
Within ENZO, we employ the piecewise parabolic method (PPM) for hydrodynamics (P. Colella & P. R. Woodward 1984; G. L. Bryan et al. 1995) and a two-shock approximate Riemann solver. In cases where higher-order methods result in negative densities or energies, the solver progressively falls back to more diffusive Riemann solvers to maintain stability. The total gravitational potential is calculated from gas self-gravity, star particles, and dark matter. We compute self-gravity using a multigrid Poisson solver, while the collisionless star particles are evolved with an adaptive particle-mesh N-body solver at an effective force resolution of approximately 2Δx, where Δx is the local cell size.
Mesh refinement is employed to ensure that the thermal Jeans length is resolved by a minimum of eight cells. Refinement continues in any given region until this criterion is met or the region reaches the maximum resolution of 1 pc. At this maximum resolution, the Jeans length may become underresolved, potentially leading to artificial numerical fragmentation. To mitigate this, we follow the guideline from J. K. Truelove et al. (1997), which requires the Jeans length to be resolved by at least four cells to suppress such fragmentation.
Unlike in E19, we do not use an artificial pressure floor. Instead, our simulations operate at sufficiently high resolution, allowing us to accurately resolve dense gas regions. This high resolution enables the efficient conversion of dense gas into stars and ensures that stellar feedback processes, such as SN explosions and stellar winds, effectively disperse the remaining gas.
We use a slightly modified version of GRACKLE (B. D. Smith et al. 2017) to follow both the nine-species nonequilibrium chemistry network including H, H + , He, He + , He ++ , e -, H -, H 2 , and H + 2 and the effects of radiative cooling and heating from metal lines and an ultraviolet (UV) background. We include the effects of H 2 formation on dust using a broken power-law dust-to-gas ratio from A. Rémy-Ruyer et al. (2014) and a F. Haardt & P. Madau (2012) UV metagalactic background accounting for self-shielding in the H I band and propagating its impact self-consistently to other photoionization and photodissociation reaction rates and the metal line cooling (see E19 for more details). We include a photoelectric heating model for far-UV (FUV) radiation from both individual stars using the same dust-to-gas ratio scaling mentioned above and a local attenuation approximation.
New in the AEOS simulations are (1) the additional contribution of the UV background to the FUV band, (2) the use in the photoelectric heating rate of a constant efficiency parameter instead of one that depends on local gas density, (3) the effects that both FUV and infrared (IR) radiation have on H 2 and H -reaction rates, and (4) a UV background that has been extended to high redshift in all bands by adopting rates assuming a blackbody radiation spectrum at 3 × 10 4 K that turns on at z = 50 and is scaled to be continuous with F. Haardt & P. Madau (2012) at z = 10. The UV background is included to account for radiation from stars beyond the edges of the simulated volume. In addition, we use an updated Lyman-Werner (LW) background model from both E19 and J. H. Wise et al. (2012a), adopting the rates at high redshift from Y. Qin et al. (2020).
For the photoelectric heating rate, we estimate the efficiency as ò = 0.05. Calculating the photoelectric heating efficiency requires a determination of the electron number density n e , which is challenging in dense and neutral regions. In these environments, n e is largely influenced by the ionization of carbon, dust grains, and polycyclic aromatic hydrocarbons. Our existing chemical network only accounts for electrons contributed by hydrogen, helium, and molecular hydrogen. To address this gap, we considered a power-law relationship for the heating efficiency ò based on the hydrogen number density n H , as derived from the M. G. Wolfire et al. (2003) model of Γ Pe in the solar neighborhood, as also done in E19: = n 0.0148 H 0.235 . For densities ranging from 1 to 10,000 cm -3 , this ranges from ò = 0.015 to 0.12, with typical values around 0.05. We thus take 0.05 as our efficiency.
Stellar winds and SNe are the two sources of both mass and energy feedback included in these simulations. We consider two distinct types of stellar winds: asymptotic giant branch (AGB) winds and winds from massive stars, as well as both core-collapse SNe (CCSNe) and SNe Ia. The exact mass, metal abundances, and energy deposited by each of these are discussed in Section 2.3. Here we discuss how each is deposited onto the computational grid.
We have sufficiently high resolution (1 pc) in these simulations to reliably resolve the Sedov-Taylor phase of a majority of our SNe for the typical gas densities in which they explode (see E19; see also M. C. Smith et al. 2018;C.-Y. Hu 2019). For this reason, we include only the thermal energy deposition from these events. Each particle deposits mass and energy feedback over a 2 pc radius spherical region centered on each particle. We use a Monte Carlo volume overlap calculation to compute the fractional deposition of mass and energy to grid cells that sit on the boundary of the spherical region.
Similarly, stellar winds are deposited in the same 2 pc radius region. The mass-loss rates for both AGB and massive-star winds are adopted from stellar evolution models but are assumed to have fixed velocities and constant loss rates over their AGB phase or lifetime.
The computational expense of fully resolving fast (10 3 km s -1 ), hot (10 6 K) stellar winds that are continually injected onto the grid is too onerous for long-timescale, galaxyscale simulations, as discussed in E19. For that reason, we fix the wind velocity for all massive stars to a maximum value of 100 km s -1 and fully thermalize the kinetic energy before injection. While this model would reduce the dynamical impact of winds on the evolution of our galaxies, we expect that stellar winds are subdominant to both stellar radiation feedback and SNe, particularly at low metallicity, so we argue that this is a reasonable approximation.
Stars in our simulation form stochastically in cold, dense gas that exhibits a converging flow ∇ • v < 0, assuming that the local star formation rate is proportional to an efficiency per freefall time e ff = 2%. In this work, we allow star formation below T thresh = 500 K and adopt a high density threshold n thresh = 10 4 cm -3 . Stars are formed when at least 100 M e of gas is available that meets star formation conditions. The IMF is sampled, stochastically forming stars until the gas reservoir is depleted. We distinguish between Population III and Population II star formation based on total gas metallicity. Gas enriched with Z > 10 -5 Z e forms Population II stars (R. Schneider et al. 2012;A. P. Ji et al. 2014;G. Chiaki et al. 2015), while gas below this threshold forms Population III stars (J. Tumlinson 2006). For Population III star formation, we place an additional constraint that the molecular hydrogen fraction
, which is consistent with the f H 2 derived from high-resolution simulations of Population III star formation at our adopted threshold number density n thresh (H. Susa et al. 2014;M. Kulkarni et al. 2021). We describe the behavior for each population in more detail below.
Our cosmological simulations are star by star for both Population III and Population II stars. This means that in each star formation event stellar masses are sampled from an adopted IMF and assigned to individual, distinct stellar particles. The exception to this is Population II stars below 2 M e ; these stars are aggregated into a single particle during a star formation event because they do not have significant feedback or enrichment on the timescales of these simulations.
The IMF for Population III stars remains uncertain (e.g., V. Bromm 2013; M. A. Latif et al. 2022; R. S. Klessen & S. C. O. Glover 2023). We adopt the same IMF used in J. H. Wise et al. (2012b), which behaves as the E. E. Salpeter (1955) IMF with power-law slope α = -1.3 above a characteristic mass M char and has an exponential cutoff below M char . Motivated by S. Hirano & V. Bromm (2017), V. Bromm (2013), and N. Yoshida (2006) and by results from our own trials with varying parameter choices, we adopt M char = 10 M e , with a range of Population III stellar masses of 1-100 M e . All Population III star particles represent individual stars over this mass range. We also run an additional simulation with M char = 20 M e that we discuss in a separate paper (K. Brauer et al. 2025, in preparation). Population III stars are assigned lifetimes from D. Schaerer (2002).
Above metallicity Z > 10 -5 Z e , gas is considered to be metal-rich enough to form Population II stars as sampled from a P. Kroupa (2001) IMF with a mass range of 0.08-120 M e . Due to computational constraints, we restrict which stars over this mass range are followed individually to those with M > 2 M e . All stars below this threshold in a star formation event are aggregated together into a single particle. These stars can be combined because they do not have significant feedback or metal enrichment on the timescale of these simulations. Since the low-mass stars with M 1 M e are the only stars that will live to the present day, though, they are key tracers of stellar abundances in present-day low-mass dwarf galaxies. We use the zero-age main-sequence properties from the PARSEC (A. Bressan et al. 2012;J. Tang et al. 2014) stellar evolution data set to assign stellar radii, effective temperature, surface gravity, lifetimes, and the length of the AGB phase, when relevant.
We model detailed multichannel stellar feedback from each of our stars. Our feedback channels include CCSNe from Population III stars, CCSNe from Population II stars, SNe Ia, AGB winds, and massive-star winds from Population II stars. We also include stellar radiation followed in three optically thin bands (IR, FUV, and LW) and H I, He I, and He II ionizing radiation followed with an adaptive ray-tracing radiative transfer method including radiation pressure on H I. These methods are discussed in greater detail below. The yields for each of these events are given in Section 2.4.
In addition to the UV background, we follow the star-by-star radiation in six bands, separated by photon energy E ph . Due to computational constraints, we limit the number of radiation sources-while capturing the vast majority of the photon energy budget of our stars-by restricting radiation to massive stars with M * > 8 M e .
We follow the H I (E ph > 13.6 eV), He I (E ph > 24.6 eV), and He II (E ph > 54.4 eV) ionizing photons using the ENZO +MORAY adaptive ray-tracing radiative transfer model described in detail in J. H. Wise & T. Abel (2011) and G. L. Bryan et al. (2014). Briefly, this method integrates the full equations of radiative transfer, propagating photons mapped onto a HEALPIX grid and adaptively refining once the separation angle between photon packages becomes large.
In addition, we track the stellar IR (0.76 eV < E ph < 5.6 eV), FUV ( 5.6 eV < E ph < 11.2 eV), and LW (11.2 eV < E ph < 13.6 eV) radiation using an optically thin approximation. This allows us to follow local variations in the H 2 (LW), + H 2 (IR, FUV, and LW), and H -(IR) photodissociation rates from each band, in addition to the localized photoelectric heating from stellar FUV radiation.
We use the table of binned photon counts from A. Heger & S. E. Woosley (2010) with the lifetimes in D. Schaerer (2002) to compute the constant photon fluxes for our Population III stars in each radiation bin (IR, FUV, LW, and H I, He I, and He II ionizing radiation) as a function of stellar mass. In practice, this is implemented using a piecewise polynomial fit to these tables.
Population III CCSNe with 10 M e < M * < 100 M e explode with a fixed energy of 10 51 erg.
We use the PARSEC (A. Bressan et al. 2012) grid of stellar evolution tracks to set the lifetime of each star and the start time and length of the AGB phase, if present. This is also used to set the stellar effective temperature, surface gravity, and radiuseach of which remains fixed at its zero-age main-sequence value-which are in turn used to set the radiation properties of each star. Photon fluxes in each radiation band are determined using the OSTAR2002 (T. Lanz & I. Hubeny 2003) grid of O-type stellar models.
However, this table does not have complete coverage over all possible stellar properties encountered in these simulations, particularly for stars below about 15 M e and very massive stars with subsolar metallicity. For stars off of the grid, we adopt a blackbody spectrum with rates scaled to be continuous with the OSTAR2002 grid (see Appendix B of E19). Ionizing photon energies are taken to be the average ionizing photon energy for the corresponding blackbody spectrum of each star. Stellar wind velocities are fixed to 20 km s -1 for AGB stars (M * < 8 M e ) and 100 km s -1 (our wind velocity ceiling; see Section 2.3) for massive stars (M * > 8 M e ). CCSNe occur for stars between 8 and 25 M e with an energy of 10 51 erg, and we assume that stars more massive than 25 M e directly collapse into black holes with no mass or energy feedback (M. Limongi & A. Chieffi 2018).
In E19, we used a power law to describe the delay time distribution (DTD) for the occurrence of SNe Ia assuming a single formation channel. We update our prescription by adopting the standard DTD from A. J. Ruiter et al. (2011, their model A1), which provides the total SN Ia DTD as the sum of four different channels: (1) double-degenerate scenario, (2) single-degenerate scenario, (3) helium-rich donor scenario, and (4) a sub-Chandrasekhar-mass scenario.
Our SN Ia prescription utilizes the initial mass-final mass relation of J. D. Cummings et al. (2019) to assign masses to the corresponding white dwarf particles once stars below 8 M e have reached the end of their lives. Following E19, we assume that stars with initial masses of 3-8 M e form white dwarfs capable of exploding as SNe Ia. Given this, the DTD, and our IMF, the fraction of stars capable of forming SN Ia progenitors that will explode in a Hubble time is 0.1508 (see Equation (2) in E19). We pre-tabulate the cumulative probability distribution for both the total DTD and each underlying DTD. When a white dwarf forms, we use a random number draw over the total DTD to set the time (if any) that each SN Ia candidate will explode and make a separate random number draw to decide which type it will be. For simplicity, we treat the total energy output for each SN Ia as being the same 10 51 erg and differentiate them only by their yields (see Section 2.4.3).
We pay careful attention to capturing the detailed chemical evolution driven by nucleosynthesis from distinct yield channels in both Population III and Population II stars, as detailed below. For a figure showing all yields as a function of progenitor mass, see the Appendix.
In total, we track 10 individual metal abundances (in addition to H, He, and the total metallicity): C, N, O, Na, Mg, Ca, Mn, Fe, Sr, and Ba. This well samples elements from each nucleosynthetic channel, in addition to capturing elements with different mass and metallicity dependence in each channel. Oxygen (O), magnesium (Mg), and calcium (Ca) are produced predominately in CCSNe and show a noticeable evolution with SN progenitor mass, tracing short-timescale (10 Myr) chemical evolution. Iron (Fe) is produced in both CCSNe and SNe Ia, and the relative abundances of O, Mg, and Ca to Fe trace the evolution between these two nucleosynthetic sources on timescales of 0.1-1 Gyr. Nitrogen (N), strontium (Sr), and barium (Ba) trace s-process enrichment in low-mass AGB stars on timescales of 0.1-1 Gyr. N and Ba trace the most massive (4-8 M e ) AGB stars, while Sr traces the less massive (<4 M e ) ones. Sodium (Na) is produced in AGB stars as part of the NeNa cycle, mostly in intermediate-mass (~4 M e ) AGB stars. Carbon (C) is significantly produced in both low-mass AGB stars and CCSNe, and it is also an important tracer of early Population III enrichment (in the form of carbon-to-iron ratios). Manganese (Mn) is predominantly formed in SNe Ia. These elements are readily observed in stellar spectra, with the exception of N and O, with O being the primary tracer of gasphase abundances.
In addition, we follow tracers tracking the total metal mass in each cell from each yield source in our chemical evolution model: Population III CCSNe, AGB winds, massive-star (M > 8 M e ) winds, Population II CCSNe, and SNe Ia. Our SN Ia prescription (Section 2.4.3) includes four metal tracers for different SN Ia progenitor types. We additionally include an r-process yield tracer to allow post-processing of r-process abundances. In total-counting the total metallicity tracer-we follow 20 metal tracers.
For the CCSNe from Population III stars (10 M e < M * < 100 M e ), we adopt the yields from A. Heger & S. E. Woosley (2010) with standard 0.1 mixing. While the exact fate of Population III stars in the range 70-120 M e is uncertain-with some possibly exploding as CCSNe and others undergoing direct collapse with no yield return-it is reasonable to approximate that all of these stars, at least up to 100 M e , end their life in a CCSN event (S. E. Woosley 2017).
For AGB winds (M * < 8 M e ), we adopt the yields of S. Cristallo et al. (2015), with eight grid points in M * ä [1.3, 6.0] M e and 10 grid points in Z * ä [1 × 10 -4 , 0.02]. For the winds and CCSN yields of massive stars, we adopt M. Limongi & A. Chieffi (2018), with nine grid points over M * ä [13, 120] M e and four grid points in Z ä [3.236 × 10 -5 , 0.01345]. Stellar yields are interpolated linearly between mass and metallicity grid points in each of the tables. For stars with masses outside the mass range sampled by the yield tables, we adopt the abundance ratios of the nearest grid point and scale the yield mass linearly with stellar mass. Yields for stars with metallicities outside the covered range are taken to be the same as the yield of the closest grid point with no extrapolation in Z.
The yield models from M. Limongi & A. Chieffi (2018) are presented for three different stellar rotations. Rather than accounting for these differences live in our simulations, we adopt a precomputed mixture model representing a populationaveraged yield set using the metallicity-dependent stellar rotation population fractions from N. Prantzos et al. (2018).
In order to fully sample the variations in N. Prantzos et al. (2018) with metallicity, we precompute an interpolated mixture model using an additional three evenly log-spaced metallicities in between the four existing grid points for a total of 13 metallicities. We adopted our particular set of elemental yields based on comparing the results of a one-zone galactic chemical evolution (GCE) model as applied to a Milky Way-mass galaxy. While these yields generally produced reasonable agreement in this model as compared to observations in [X/Fe] versus [Fe/H] space, Mg is noticeably underproduced in [Mg/Fe] at all [Fe/H], while other α-elements tend to agree well with observed Milky Way chemical abundances. Given the importance of accurate observational comparisons in our simulations, we applied a uniform factor of 2.2 to the Mg yields from all massive stars, which brings the model into agreement with observed Mg abundances. This is discussed in the Appendix.
As discussed in 2.3.3, we use a combined DTD from four different sources of SNe Ia, each with potentially unique abundance signatures. However, given the uncertainty in yields from each of these sources, we opt instead to make the assumption that there is one yield pattern for each source, which allows us to post-process the abundance patterns from each channel separately. Live in the simulation, we assume a single abundance pattern for all SNe Ia from F.-K. Thielemann et al. (1986) and track the contribution of each SN Ia type to the total metallicity as a separate passive scalar tracer field. Given this, and knowing the total number of each SN Ia type that has occurred in the simulation, one can arbitrarily rescale the abundance patterns for each SN Ia type. Doing so implicitly assumes that the yields for each SN Ia do not affect the dynamical evolution of our galaxies (which, in reality, it may, for example, by influencing cooling through the Fe atomic line cooling), since we do not account for local cooling variations due to individual elemental abundances in our simulations.
For the simulations presented here and in our companion paper (K. Brauer et al. 2025, in
Our initial simulation has a comoving (1 Mpc) 3 volume, simulated from redshift z = 130 to z = 14.5 (~300 Myr after the big bang). It has a root-grid resolution of 256 3 , a dark matter resolution of 1840 M e , and a physical 1 pc resolution of the gas at the finest scales. All Population II stars with masses greater than 2 M e and all Population III stars are represented by singlestar particles. The stars and gas have 20 metal tracer fields tracing 10 individual metal abundances and several yield sources (see Section 2.4). Projections of the fiducial volume can be seen in Figure 1.
Because computational limitations do not allow us to simulate the full domain beyond redshift z ~14, our initial analysis focuses on the first 300 Myr of the Universe. We also run an additional AEOS simulation with a different Population III IMF and several comparison simulations without individual stellar feedback; these simulations are discussed and compared in our companion paper (K. Brauer et al. 2025, in preparation). Future work will include zoom-in simulations of ultrafaint dwarfs that will run until reionization.
At redshift z = 14.5, the full 1 Mpc volume contains 91 starforming halos of at least M halo > 2 × 10 5 M e (our resolution limit) with a total of about 250,000 star and stellar remnant particles (see Figure 1). Population III star formation begins at z ~28 (110 Myr after the big bang), and Population II star formation begins a bit after z ~22 (160 Myr after the big bang). The mass of Population II stars, defined as stars with total metals Z 10 -5 Z e , overtakes the mass of Population III stars around z ~17. Figure 2 shows the stellar masses and halo masses of every star-forming halo in the simulations at z = 14.5 and the stellar mass-halo mass relation. Some galaxies in our simulations exhibit unexpectedly high stellar-to-halo mass ratios (such as M * = 10 5 M e for M halo = 10 6 M e ). Abundance matching relations are very poorly calibrated for low-mass halos, but they predict closer to M * = 10 2 -10 3 M e or even lower for a 10 6 M e halo (e.g., P. S. Behroozi et al. 2013). Most of our galaxies are indeed low in stellar mass, and we have significant scatter in stellar mass, which is expected with lower halo masses (S. Garrison-Kimmel et al. 2017). For the massive outliers, their high stellar-to-halo mass ratios may be due in part to environmental factors. Two of the outlier galaxies are the externally enriched galaxies near the central galaxy (see Section 3.1), and their total halo mass is affected by their orbit around a more massive galaxy, allowing for dark matter stripping (a lesser version of the stripping processes described by J. Moreno et al. 2022). The unexpectedly high stellar-to-halo mass ratios may also be due to our star-by-star IMF sampling approach. Recently, M. Jeon & M. Ko (2024) found that star-by-star sampling leads to higher stellar masses in low-mass halos (e.g., 10 8 M e ) compared to traditional single stellar population (SSP) particles, which inject feedback more burstily and suppress star formation more effectively. M. Jeon & M. Ko (2024) demonstrated that the SSP method results in burstier, stronger feedback, while our star-by-star method, with its weaker feedback, allows for higher stellar masses. They also argue that star-by-star sampling more accurately matches observations and accounts for the high scatter in the M halo -M * relation at low halo masses. Similarly, E. P. Andersson et al. (2025) report a systematic increase in stellar mass for a given halo mass compared to classical galaxy models (see Figure 1 of E. P. Andersson et al. 2025). In our other recent AEOS paper (J. Mead et al. 2025), we observe that a single SN disrupts star formation in most galaxies, but with this effect diminishing around 10 7 M e , where our outliers lie.
Of the 91 star-forming halos, the vast majority are tiny. A total of 60 of the halos contain stellar masses of M * 100 M e , while only one galaxy has M * > 10 6 M e . A total of 17 of the galaxies have begun forming Population II stars by z ~14.5. For each of these galaxies, the distribution of Population II versus Population III stellar mass is shown in Figure 3. Galaxies typically undergo multiple episodes of Population III star formation before transitioning to Population II stars. The extent of Population III star formation required for this transition is strongly influenced by the halo mass, which affects the galaxy's ability to retain metals, as well as the star formation activity of neighboring galaxies. For example, in Halo 1 of Figure 3, approximately 1000 solar masses of Population III stars formed over an 80 Myr period before the onset of Population II star formation. Other galaxies exhibited different amounts of Population III star formation, ranging from none or only a few hundred solar masses-particularly in cases of external enrichment, more massive halos, or proximity to other star-forming galaxies-to up to about 1000 solar masses in galaxies that began star formation earlier than their surroundings or were relatively isolated. This variation highlights the significant role that both intrinsic halo properties and environmental factors play in determining the transition from Population III to Population II star formation.
Two galaxies, Halo 3 and Halo 5, were externally enriched in metals by feedback-driven outflows from a larger galaxy, Halo 1. This caused their star formation to begin with Population II stars. This is seen in Figure 3 and shown in more detail in Figure 4. Both Halo 3 and Halo 5 sit just outside the virial radius of Halo 1 and will merge in the future, but even before this future merger they all share gas and metals.
Due to the shallow potentials of their small halos, the AEOS galaxies struggle to retain their gas and metals (J. Mead et al. 2025). Every SN explosion, in the smallest halos, blows out the galaxy's gas, causing neighboring galaxies to freely share metals and gas between themselves. This leads to behavior like the external enrichment shown in Figure 4.
In these simulations, the halos of early galaxies (z 14) are all low in mass (generally M dm 10 7 M e ). Any individual galaxy experiencing an SN explosion thus strongly interacts with any galaxies around it (within ~5 physical kpc). This raises the question of how to treat early galaxies like Halos 1, 3, and 5: as individuals or as systems. Based on the communal sharing of gas and metals between the galaxies in these simulations, we suggest that the earliest galaxies be thought of not as individuals but as systems that contribute their metals to a common reservoir and evolve together.
Examples of external enrichment to the degree of a galaxy beginning its star formation with Population II stars are still rare, however, as evidenced by Figure 3. In this simulation, there are two examples of galaxies that begin star formation with Population II out of 17 galaxies that have begun forming Population II stars. This depends strongly on the amount of clustering. The two externally enriched galaxies exist within a few physical kiloparsecs of the largest galaxy in the simulation and will likely merge with the central galaxy in time. Hence, while clusters of small galaxies within ~5 physical kpc may be considered systems rather than individuals, this does not hold for more isolated early halos.
We also look at how the metallicities differ for stars formed via external enrichment. When looking at first-generation Population II stars in the galaxies shown in Figure 4, the metallicities differ between the externally enriched halos and the central halo, but both of these can be higher or lower in metallicity. The first Population II stars in Halo 1 form at about [Fe/H] = -5, while the first stars in Halo 3 form at about [Fe/H] = -4.5. When SNe blow out the metals from Halo 1, it actually results in Halo 3 having a higher metal density. On the other hand, Halo 5 is about twice as far away (~3 kpc), and when it starts forming stars later at 260 Myr, its stars have [Fe/H] = -5.3. In this case, the metals became diffused during the transfer.
We track 10 individual metal abundances, as described in Section 2.4. These elements trace different nucleosynthetic channels with stars of different mass and metallicity within any given channel. The most important novelty of the AEOS simulations is the ability to trace detailed metal enrichment from individual stars in small galaxies-we do this in a cosmological context with radiative transfer and include enrichment from Population III stars. This star-by-star resolution is necessary to uncover scatter and structure in stellar chemical abundance space.
As a proof of concept, consider the merger of three galaxies shown in Figure 5. In simulation, three early galaxies merge at z ~15. At the time of merger, these galaxies contain only 6000, 2000, and 15,000 M e of stellar mass, respectively, before merging, and after merging they experience bursty star formation until the merged halo has a stellar mass of 10 7 M e at the end of the simulation (Halo 1 at z = 14.5). This is a clustered region of the simulation.
Figure 5 shows the individual star particles of these galaxies, colored according to their progenitor galaxy. As the galaxies merge, star formation increases and the merged galaxy continues to form stars (shown in orange). Here "merged" means that the star-forming regions are indistinguishable from one another at the physical scale of 0.1 kpc, but note that the regions are still oscillating and the halos could still move out and back into Halo 1 before becoming permanently indistinguishable. The projections and stellar birth times shown in Figure 5 can only ever be known in a simulation, never in any observed data.
The stellar chemical abundances shown in Figure 6, however, could be measured in observed stars with spectroscopy. This minimal contribution from AGB winds. In the last ~15 Myr of the simulation, Population II CCSNe take over as the dominant source of metals in these galaxies. For the elements shown in Figure 6, carbon has the greatest AGB contribution, but the mass of carbon from AGB stars is still almost four orders of magnitude lower than the mass from CCSNe (see Section 4.2). The spread in strontium relative to iron is due to strontium's more significant production in Population II AGB winds and CCSNe compared to Population III CCSNe (see Figure A1), exacerbated by the overall extremely low amounts of strontium in the gas (~10 -4 M e of strontium in the entire system). Note that strontium is primarily a tracer of AGB winds but forms in small amounts in CCSNe owing to the weak s-process, and the metal content of Population II stars provides seed nuclei that Population III stars lack. As time continues to pass, winds and SNe Ia will become significant and contribute to more scatter, even within individual progenitor galaxies (such as the spread in strontium seen in Figure 6). There is also a spread seen in [N/Fe] in the different progenitor galaxies, particularly Progenitor 2. This is due to yield differences from differentmass CCSNe, producing an intrinsic source of scatter.
Figure 6 demonstrates that structure exists in chemical abundance space and can correspond to the progenitor origins of the stars. In this simple example, the stars from different progenitor galaxies visibly exist in different regions of chemical abundance space. The multidimensional abundance space contains rich information about the environments in which the stars formed. In this way, stellar chemical abundances encode information about hierarchical galaxy formation.
In this example, the abundance scatter corresponds to differences in early galaxy formation, but scatter can also be due to differences in nucleosynthetic yields (as we have already alluded to when discussing the spreads of strontium). The current simulations have only run to z = 14.5, but the next iteration of the AEOS simulations will run to reionization. This will capture a rich history of different nucleosynthetic events including CCSNe, SNe Ia, s-process elements from AGB stars, and r-process events that will be included in post-processing. With the level of detail shown in Figures 5 and 6, we will be able to identify how the spreads of abundance distributions differ with variations in galaxy evolution (e.g., merger history, star formation history) and also with variations in the different nucleosynthetic events experienced by each galaxy (e.g., different amounts of barium produced in s-versus r-process events or stochastic differences in yields from single sources).
Carbon-enhanced metal-poor (CEMP) stars are a unique and crucial population in the study of stellar evolution and galactic archeology. Observationally, these stars are characterized by their unusually high carbon-to-iron ratios ([C/Fe] > +1.0) in comparison to typical metal-poor stars. CEMP stars are of particular interest because they offer significant insights into the early universe, the processes of nucleosynthesis, and the formation of the first stars and galaxies (T. C. Beers & N. Christlieb 2005).
In this system of galaxies (see Figures 5 and 6), firstgeneration Population II CEMP stars are observed exclusively at extremely low metallicities. The highest stellar metallicity for these CEMP stars is [Fe/H] = -4.3, while the median metallicity of CEMP stars in this halo is [Fe/H] = -5.
More broadly in the simulation, about half of low-metallicity stars that would survive to present day exhibit a CEMP signature. Below [Fe/H] = -4, 54% of the low-mass stars (sub-solar-mass stars that will survive to present day) in the simulation present a CEMP signature. The fraction of stars with a CEMP signature decreases to 29% of the low-mass stars below [Fe/H] = -2.5 and 2% of the low-mass stars below [Fe/H] = -1.5. This is generally in agreement with observations of metal-poor stars that suggest that a significant fraction of metal-poor stars are carbon-rich (J. E. Norris et al. 2013), a fraction that increases with decreasing metallicity and may be as high as ~80% for stars below [Fe/H] = -4 (V. M. Placco et al. 2014).
The plot in Figure 7 . Left: around 160 Myr after the big bang, three small galaxies merge in the simulation. Each progenitor enters with its own stars, with colors labeled in the right panel, and the merged galaxy continues to form stars. Right: we show the birth time and metallicity of each star in this system; the galaxies formed stars at slightly different metallicities prior to merging. Each star is colored by its galaxy of origin. We investigate more chemical abundances in Figure 6 to show that structure in chemical abundance space can correspond to the origins of the stars.
In Figure 7, we also compare to the observed frequencies of CEMP stars in the Milky Way (V. M. Placco et al. 2014). At the lowest metallicities, we underpredict CEMP stars with [C/Fe] > 1.0 when compared to Milky Way observations. This could indicate the need for a Population III IMF that favors higher-mass Population III, because higher-mass CCSNe produce ejecta with a higher ratio of carbon to iron (see Figure A2 and A Because AEOS directly models Population III enrichment, we can determine the median chemical abundances of gas in Figure 6. Example of structure in chemical abundance space. For the system of merging galaxies shown in Figure 5, we plot the chemical abundances of stars (only showing low-mass stars that will still be alive at present day), colored by whether the star formed in Progenitor 1 (blue), Progenitor 2 (purple), Progenitor 3 (green), or later in the fully merged galaxy (orange). The data are binned on expected observational uncertainty for stars in the Milky Way halo today (the gray-white grid). The size of the circle in each bin is proportional to the log of the number of stars in that bin. In this simple example, the location of stars in multidimensional chemical abundance space clearly corresponds to which halo the stars were born in. our Population III halos and the scatter about that median. Figure 8 shows the [X/H] chemical abundances of gas in all pure Population III halos at five different snapshots (equally spaced in time from about 150 to 300 Myr) with 16th-84th percentile scatter. Even with increasing halo mass, the floors are fairly flat, though scatter is significant. We also find that the floors are fairly flat with time, with the most scatter at early times (<175 Myr). These floors depend on Population III yields (A. Heger & S. E. Woosley 2010) and critical metallicity of the transition from Population III to Population II stars (Z < 10 -5 Z e ).
As expected, most of the metal mass in the gas is oxygen, which is the most abundantly produced element in Population III CCSNe (see Figure A1). We find a best-fit metallicity floor of about [O/H] = -4.0, which is lower than that assumed in EDGE, but they found that changing the floor to Z = 10 -4 Z e did not significantly affect their results (O. Agertz et al. 2020). For all of the best-fit metallicity floors, see the red lines in Figure 8.
Our simulation supports a Population III metallicity floor of Z = 10 -4 Z e that generally holds for halo masses of 10 5 -10 7 M e . At the highest halo masses, the metallicity of the gas increases, but there are only four halos above 5 × 10 6 M e , so this may or may not be robust. The best-fit metallicity floor of each element is approximately as follows:
, and [Fe/H] = -5.1. There is significant scatter, however, due to inhomogeneous mixing that follows localized enrichment from different SNe and explosive outflows. We exclude Sr and Ba, as they are not produced in significant amounts by Population III.
We also look at the stellar chemical abundances of the first low-mass Population II stars in each galaxy to see how the gas enrichment translates into stellar enrichment for the firstgeneration Population II stars. A total of 17 galaxies have begun Population II star formation (see Figure 3), and we plot the chemical abundances of the first-generation low-mass Population II stars from each Population II galaxy in Figure 9. These abundances are similar to the gas metallicity floors in Figure 8, as we would expect, with some differences (e.g., some of the median
In this paper, we introduce the AEOS project, a series of simulations designed to trace early galaxy formation processes and model the chemical enrichment of individual stars in unprecedented detail while including Population III enrichment, radiative transfer, and galaxy formation in a cosmological context. Our simulations aim to shed light on the spread of observed stellar chemical abundance patterns in ultrafaint dwarf galaxies, metal mixing in the ISM, the impact of Population III stars, and more. We describe the methods and novel handling of star-by-star chemical enrichment.
By modeling individual stars with their chemical yields, we can capture the intricate interplay between galaxy evolution and different nucleosynthetic processes. As a proof of concept, we demonstrate that the detailed element tracking in the AEOS simulations allows us to reproduce scatter and structure in stellar chemical abundance space that can be directly linked to the progenitors of the stars. This star-by-star approach thus allows us to decode valuable information about hierarchical galaxy formation and early nucleosynthetic events.
We also explore the concept of how to treat small early galaxies: as individuals or as systems. We demonstrate examples of external enrichment in the simulation, where a more massive galaxy enriches the gas of nearby galaxies through outflows. The earliest galaxies live in low-mass halos that lose significant amounts of gas and metals, freely sharing with the galaxies near them. We consider these galaxies not to be individuals but to be systems that functionally evolve together.
Additionally, we estimate metallicity floors for different metals (e.g., [O/H] = -4) as a result of Population III enrichment. These metallicity floors are fairly flat even with variations in halo mass and time, though with significant scatter. There is also over a dex of scatter in the median metallicity of first-generation Population II stars from different galaxies, showing that significant variation exists between first-generation Population II stars. We also investigate CEMP signatures in the first generation of Population II stars and compare to the frequencies of CEMP stars observed in the Milky Way.
The current simulations have been run from redshift z = 130 to z = 14.5, and future work will extend this work to the epoch of reionization for a suite of zoom-ins on a select sample of small galaxies that should be early analogs of the surviving ultrafaint dwarf galaxies. This will enable us to capture a comprehensive history of galaxy growth and different nucleosynthetic events that drove the associated processes. By exploring variations and interdependence of galaxy evolution and nucleosynthetic yields, we will quantify the origins of observed abundance scatter, the metal retention of small bursty galaxies, and the impact of metal mixing in the early ISM.
The Astrophysical Journal, 980:41 (16pp), 2025 February 10 Brauer et al.