# EDGE-INFERNO: Simulating Every Observable Star in Faint Dwarf Galaxies and Their Consequences for Resolved-star Photometric Surveys

Eric P Andersson, Martin P Rey, Andrew Pontzen, Corentin Cadiou, Oscar Agertz, Justin I Read, Nicolas F Martin

Agertz, Justin I Read, Nicolas F Martin

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

## To cite this version:

Eric P Andersson, Martin P Rey, Andrew Pontzen, Corentin Cadiou, Oscar Agertz, et al., EDGE-INFERNO: Simulating Every Observable Star in Faint Dwarf Galaxies and Their Consequences for Resolved-star Photometric Surveys. The Astrophysical Journal, 2025, 978, 10.3847/1538-4357/ad99d6. insu-04918555

HAL Id: insu-04918555

https://insu.hal.science/insu-04918555v1

Submitted on 29 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# EDGE-INFERNO: Simulating Every Observable Star in Faint Dwarf Galaxies and Their Consequences for Resolved-star Photometric Surveys

Eric P. Andersson<fnref n="1" />, Martin P. Rey<fnref n="2" />,<fnref n="3" />, Andrew Pontzen<fnref n="4" />, Corentin Cadiou<fnref n="5" />, Oscar Agertz<fnref n="5" />, Justin I. Read<fnref n="6" />, and Nicolas F. Martin<fnref n="7" />,<fnref n="8" />

Received 2024 September 14; revised 2024 November 15; accepted 2024 December 01; published 2025 January 3

## Abstract

Interpretation of data from faint dwarf galaxies is made challenging by observations limited to only the brightest stars. We present a major improvement to tackle this challenge by undertaking zoomed cosmological simulations that resolve the evolution of all individual stars more massive than 0.5 $M_\odot$, thereby explicitly tracking all observable stars for the Hubble time. For the first time, we predict observable color–magnitude diagrams and the spatial distribution of ≈100,000 stars within four faint ($M_* \approx 10^5 M_\odot$) dwarf galaxies directly from their cosmological initial conditions. In all cases, simulations predict complex light profiles with multiple components, implying that typical observational measures of structural parameters can make the total V-band magnitudes appear up to 0.5 mag dimmer compared to estimates from simulations. Furthermore, when only small (≈100) numbers of stars are observable, shot noise from realizations of the color–magnitude diagram introduces uncertainties comparable to the population scatter in, e.g., the total magnitude, half-light radius, and mean iron abundance measurements. Estimating these uncertainties with fully self-consistent mass growth, star formation, and chemical enrichment histories paves the way for more robust interpretation of dwarf galaxy data.

*Unified Astronomy Thesaurus concepts: Galaxy formation (595); Dwarf galaxies (416); Galaxy properties (615)*

## 1. Introduction

The advent of deeper, wider photometric surveys in the last two decades (e.g., the Sloan Digital Sky Survey (SDSS) and the Dark Energy Survey (DES)) has revealed many, new, ever-f fainter dwarf galaxies around the Milky Way (e.g., V. Belokurov et al. 2007; A. Drlica-Wagner et al. 2015; S. E. Koposov et al. 2015; see J. D. Simon 2019 for a review).

These tiny galaxies ($M_* \lesssim 10^6 M_\odot$) likely inhabit low-mass dark matter halos ($\approx 10^8 M_\odot$; P. Jethwa et al. 2018; E. O. Nadler et al. 2020). Their mere existence provides sensitive constraints on dark matter models that suppress the growth of small-scale structures (e.g., warm dark matter and axions; E. O. Nadler et al. 2019). Furthermore, their low metallicities ([Fe/H] ≈ −2) provide a laboratory to study the first generations of stars and their chemical enrichment (A. Frebel & J. E. Norris 2015; A. P. Ji et al. 2015). Understanding how faint dwarf galaxies form and assemble in a Lambda cold dark matter (ΛCDM) universe is thus of broad importance and urgency, as the forthcoming Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) is expected to vastly expand the census of faint dwarf galaxies in the coming years (e.g., O. Newton et al. 2018; B. Mutlu-Pakdil et al. 2021).

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

CC BY

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

In recent years, numerical models of low-mass galaxies have improved significantly (see, e.g., C.-Y. Hu et al. 2016; C. Wheeler et al. 2019; O. Agertz et al. 2020; M. C. Smith 2021), even unlocking the ability to account for individual stars within a whole galaxy (e.g., C.-Y. Hu et al. 2017; A. Emerick et al. 2018; E. P. Andersson et al. 2020, 2021, 2023; N. Lahén et al. 2020, 2023, 2024a; F. Calura et al. 2022; T. A. Gutcke et al. 2021, 2022a, 2022b; J. M. Hislop et al. 2022; U. P. Steinwandel et al. 2023, 2024; Y. Deng et al. 2024), as opposed to aggregating stars into a single stellar population particle. By explicitly drawing from the initial mass function (IMF) at stellar birth and modeling the stellar evolution tracks of each star, these models provide a much more detailed and robust account of how, when, and where feedback and metals are deposited into the interstellar medium (ISM). Furthermore, this novel model approach opens up the door for incorporating aspects of collisional dynamics present in star clusters that potentially affect both the ISM, galactic winds, and circumgalactic medium (e.g., runaway stars and clustering of feedback; see T. Naab & J. P. Ostriker 2017 for a review). While high-order direct N-body integration is prohibitively computationally expensive for all stars in a galaxy over a Hubble time, its effects can be treated with subgrid models (e.g., E. P. Andersson et al. 2020, 2023; U. P. Steinwandel et al. 2023) using calibrations from highly detailed models applied in idealized setups (M. Y. Grudić et al. 2021, 2022; N. Lahén et al. 2024b; B. Polak et al. 2024).

However, the computational expense has constrained their applications to idealized, noncosmological settings. This limits their ability to interpret redshift $z=0$ data as the cosmological mass-growth history of a dwarf galaxy is key to shaping its

<footnote n="1">¹Department of Astrophysics, American Museum of Natural History, 200 Central Park West, New York, NY 10024, USA; andersson@amnh.org</footnote>

<footnote n="2">²Department of Physics, University of Bath, Claverton Down, Bath BA2 7AY, UK</footnote>

<footnote n="3">³Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="4">⁴Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5">⁵Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden</footnote>

<footnote n="6">⁶Department of Physics, University of Surrey, Guildford GU2 7XH, UK</footnote>

<footnote n="7">⁷Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="8">⁸Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>
observables (e.g., A. Fitts et al. 2017; M. P. Rey et al. 2019, 2020, 2022; A. Benitez-Llambay & M. Fumagalli 2021).

In the rare cases when individual-star models have included a full cosmological environment (e.g., T. A. Gutcke et al. 2022a, 2022b; F. Calura et al. 2022), only massive (and thus less numerous) stars ($\gtrsim 4 M_\odot$) are individually sampled to ease the computational load. But slow-moving stellar mass loss from low-mass evolved stars has characteristic velocities ($\approx 10$ km s$^{-1}$; S. Höfner et al. 2018) comparable to the depth of gravitational potential wells of faint dwarf galaxies. As a result, winds from, e.g., asymptotic giant branch stars (AGBs) likely play a role in chemical enrichment and possibly delay the restart of star formation in galaxies quenched by reionization (e.g., M. P. Rey et al. 2020). Furthermore, the lowest-mass galaxies are star formation quenched during the epoch of reionization and hence their mass budget in massive stars by redshift zero is fully depleted. Instead dwarf galaxy properties are characterized through the imaging and spectroscopy of the light from their old, evolved low-mass stars, motivating a complete account of the stellar population.

In this paper, we combine all these advances for the first time to take a major step toward interpreting resolved-star data from dwarf galaxies. We perform high-resolution cosmological simulations within the Engineering Dwarfs at Galaxy formation's Edge (EDGE) framework, combining them with the star-by-star Individual stars with Feedback, Enrichment, and Realistic Natal mOtions (INFERNO) model (E. P. Andersson et al. 2023). With this, we achieve a high numerical resolution (3 pc at $z=0$) that allows us to capture the evolution of supernova (SN) explosions self-consistently over the full cosmological formation of the dwarf, while individually sampling all stars down to $0.5 M_\odot$ (main-sequence lifetime $> 14$ Gyr). This ensures that all major sources of stellar feedback are accounted for on a star-by-star basis, while self-consistently populating the observable spatial and color-magnitude distribution at $z=0$, at which point stars born before reionization and with mass $\sim 0.8 M_\odot$ evolve to the giant branch (thus become bright enough to be observable).

We describe these simulations and methods in Sections 2 and 3 and, in Section 4, present the observational implications of directly predicting resolved-star observables from cosmological initial conditions.

## 2. Star-by-star Cosmological Simulations

We perform zoomed, cosmological simulations of a low-mass dark matter halo ($M_{200} = 1.5 \times 10^9 M_\odot$ at $z=0$). Using the genetic modification technique (N. Roth et al. 2016; M. P. Rey & A. Pontzen 2018; S. Stopyra et al. 2021), we modify the initial conditions of this halo to engineer four different early growth histories for this object, while maintaining the same halo mass at $z=0$ and large-scale cosmological environment (see M. P. Rey et al. 2019 for further details of the genetically modified initial conditions).

Figure 1 shows the halo mass-growth history (top panel) from our new suite of simulations (filled) matching the original EDGE simulations (dotted M. P. Rey et al. 2019). We define halo mass as $M_{200}$, defined as the mass enclosed with the radius $r_{200}$ where the density is 200 times the critical density of the Universe. Genetic modifications of this halo are designed to scan through the 68%-95% population of early mass-growth histories for dark matter halos of this mass (gray contours; see details of the calculation in M. P. Rey et al. 2019), at fixed halo

mass today. Compared to randomly sampling halos, this genetically modified setup allows us to efficiently probe the scatter in observable properties while fixing the main control variable ($M_{200}$).

For this work, we reevolve these initial conditions to $z=0$ using the adaptive mesh refinement code RAMSES (R. Teyssier 2002) and the star-by-star model INFERNO described in depth in E. P. Andersson et al. (2023). The hydrodynamical, gravity, and thermochemistry setup are the same as in E. P. Andersson et al. (2023), while the refinement strategy is identical to M. P. Rey et al. (2019) and O. Agertz et al. (2020). The median resolution is 3 pc within the ISM of the galaxy (dark matter particle mass of $980 M_\odot$).

Star formation follows O. Agertz et al. (2020), assuming stars are born in quantities of $500 M_\odot$ in cells with gas density $> 300$ cm$^{-3}$ and temperature $< 10^3$ K. Spawning stars in groups of $500 M_\odot$ minimizes issues related to incomplete sampling of the IMF (M. C. Smith 2021). In addition to the gas velocity, stars are stirred with a velocity dispersion of $0.01$ km s$^{-1}$ at birth (E. P. Andersson et al. 2023). This velocity dispersion is below the resolution of the gravitational softening (cell size) and serves as a numerical approximation to allow particles with identical initial positions to end up on trajectories that eventually diverge due to small but repeated dynamical perturbation.

The key novelty of this work is the introduction of individual-star particles once a star formation has been flagged. We draw masses from a P. Kroupa (2001) IMF using the sampling technique from M. C. Sormani et al. (2017). Critically, all stars with masses $\ge 0.5 M_\odot$ are explicitly sampled and inject energy, momentum, mass, C, N, O, Mg, Al, Si, and Fe. Stars with masses $\le 0.5 M_\odot$ are expected to remain on their main sequence for the Hubble time and are thus aggregated into a low-mass, passive particle. Although important, stellar multiples (e.g., binaries) are not accounted for in our model (see further discussion in E. P. Andersson et al. 2023).

Stellar feedback in our models includes SNe (core collapse and Type Ia) and stellar winds (from OB type and AGB stars) implemented as in E. P. Andersson et al. (2023). While winds are a subdominant component of the feedback budget (in particular at low metallicity; see, e.g., J. S. Vink et al. 2001; M. R. Mokiem et al. 2007; J. S. Vink 2011), they are an important contributor to metal enrichment and can provide a prompt source of energy injection affecting local star formation (E. P. Andersson et al. 2024). Radiation is not included in our suite of simulations (or in the simulations from M. P. Rey et al. 2019 we compare to) and is expected to further regulate star formation, systematically decreasing the total stellar masses of our galaxies (O. Agertz et al. 2020; M. P. Rey et al. 2024, in preparation). We read the chemical yields and mass loss of each star as a function of mass and metallicity using tabulated values from NUGRID (M. Pignatari et al. 2016; C. Ritter et al. 2018) for core-collapse SNe and stellar winds, and from I. R. Seitenzahl et al. (2013) for Type Ia SNe. In all cases, we calculate the main-sequence lifetime of a star following the fitting function from C. M. Raiteri et al. (1996). Type Ia SNe follow the time-delay distribution for field galaxies from D. Maoz & O. Graur (2017) assuming a rate of $I_{\text{ta}} = 2.6 \times 10^{-13} \text{ yr}^{-1} M_\odot^{-1}$ starting 38 Myr after birth.

As in E. P. Andersson et al. (2023), massive stars ($8 \le m_* \le 60 M_\odot$) inject fast (1000 km s$^{-1}$) winds while on the main sequence and end their life with core collapse directly injecting material and $10^{51}$ erg if $m_* \le 30 M_\odot$ (assuming direct
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Halo mass ($M_{200}$, top panel) and stellar mass (bottom panel) as a function of redshift (bottom axis) and age of the Universe (top axis) showing the four galaxies simulated star by star (filled lines) and with star particles (300 $M_{\odot}$ resolution, dotted lines; M. P. Rey et al. 2019). The dark (light) gray-shaded region encloses 95% (68%) of halo growth histories in the entire simulation box (50 Mpc$^3$) from which our initial conditions were selected. The right y-axis of the top panel shows the fraction of mass relative to the halo mass at redshift zero.</figcaption></figure>

collapse to black hole otherwise). While the high resolution of our simulations allows direct injection of SN energy to obtain accurate capture of the final momentum (L. Ohlin et al. 2019) in almost all cases (>96%), a switch to terminal momentum injection (C.-G. Kim & E. C. Ostriker 2015) is applied in the rare cases when the resolution is not sufficient (see O. Agertz et al. 2020; E. P. Andersson et al. 2023 for further details).

Low-mass stars ($<8 M_{\odot}$) eject slow (10 km s$^{-1}$) winds with a constant mass-loss rate ($10^{-5} M_{\odot} \text{ yr}^{-1}$; K. Eriksson et al. 2014; S. Höfner et al. 2018) during their AGB phase. Wind injection stops when their stellar mass follows the initial-to-final mass relation from J. D. Cummings et al. (2016).

### 3. Stellar Mass Growth

The lower panel of Figure 1 shows the stellar mass growth of our dwarf galaxies. With $M_{*}(z=0) = 10^{5}-10^{6} M_{\odot}$ and $M_{200} = 1.5 \times 10^{9} M_{\odot}$, our galaxies are within the $2\sigma$ scatter inferred for local Milky Way satellites (e.g., E. O. Nadler et al. 2020, Figure 6).

Compared with the traditional stellar-particle approach (dotted line; M. P. Rey et al. 2019), swapping to individual-star feedback systematically increases $M_{*}$ at all redshifts, by at most a factor of 2. Such variations are well within theoretical uncertainties at this mass scale (see, e.g., the discussion in

O. Agertz et al. 2020; review in L. V. Sales et al. 2022), which should be expected given the extensive updates to the chemical enrichment channels (and thus gas cooling) and updates to the stellar feedback implementation.<fnref n="9" />

Whether using individual stars or stellar particles, earlier-forming dwarfs systematically assemble more stellar mass by $z=0$. This confirms that relative effects between mass-growth histories are robust to model changes, even if astrophysical assumptions leave residual uncertainties in the absolute normalization of the stellar mass (see also the discussion in M. P. Rey et al. 2019). Combined, these consistency checks allow us to confidently turn to harnessing the new explicit connection between assembly and resolved-star observables.

### 4. Direct Predictions of Dwarf Galaxy Resolved-star Observables

#### 4.1. Color–Magnitude Diagrams

Figure 2 highlights the vastly improved stellar sampling of our simulated galaxies, showing the surface brightness (top row) as if the galaxies are observed at 400 kpc with the resolution of DES. Each star is assigned stellar properties (e.g., luminosity) derived by interpolation of PARSEC (A. Bressan et al. 2012) isochrones (see Appendix A for details) accounting for stellar mass, age, and metallicity. In the top panels of Figure 2, low-mass stellar aggregates are accounted for by sampling the IMF in the range 0.08 and 0.5 $M_{\odot}$. Bright individual stars are visible at this distance.

The bottom row of Figure 2 shows each galaxy’s respective CMD for r- and g-band magnitudes scaled with a $D = 400$ kpc distance modulus. For each CMD, we add a scatter similar to the typical noise associated with the equivalent observation (left points in each panel; D. Homma et al. 2018). The dashed red line visualizes the typical point-source limit accessible with the DES (J. D. Simon 2019). As the stellar mass decreases (left to right), the number of stars above this observable threshold decreases (top-left number). Such observable constructs are a natural by-product of our simulations, as is predicting the spatial location and chemistry of the brightest, observable giant stars. This enables us to derive dwarf galaxy properties directly with observational tools, to which we turn now. We specifically chose $D = 400$ kpc to match the tens to hundreds of stars typically observed in dwarfs with a DES-like setup. We vary these assumptions in Appendix C.

All our dwarf galaxies have diffuse light profiles, with hierarchical merging leading to distant bound members up to $\ge 5 r_h$ similar to those reported in recent observations (e.g., A. Chiti et al. 2021; F. Sestito et al. 2023). As we will see, these multicomponent spatial distributions can bias the recovery of faint dwarf structural parameters, particularly when few stars are bright enough to be observed. Such outer stellar halo for faint dwarf galaxies was already pointed out in simulations with stellar particles (e.g., Y. Tarumi et al. 2021; A. Goater et al. 2024). In theory, postprocessing techniques sampling the IMF could be used to transform such stellar particles into individually resolved stars (e.g., M. R. Krumholz et al. 2015; R. J. J. Grand et al. 2018; R. E. Sanderson et al. 2020; E. I. Gray et al. 2024) and test observational pipelines. However, up-sampling techniques can

<footnote n="9">⁹ Even if the total feedback budget and injection mechanisms are fixed, the location and clustering of individual stars can be significantly different and affect galactic outflows (e.g., E. P. Andersson et al. 2020, 2023; M. C. Smith 2021).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** V-band surface brightness maps for each galaxy (top row) at 400 kpc distance at a resolution of 0.3 pixel⁻¹ (equivalent to DES). Compared with traditional stellar population approaches, modeling every star with $≥0.5 M_⊙$ greatly increases stellar sampling throughout the galaxy, both in the center and the outskirts. It also directly predicts the simulated dwarf's color–magnitude diagram (CMD; here scaled by a 400 kpc distance modulus in *r*-* and *g*-band magnitudes, error bars reflecting the observationally motivated applied scatter) and spatial distribution, enabling the estimation of structural properties as an observer would (red ellipses show the elliptical half-light radius; see text for details). As galaxies become less massive (left to right), fewer stars stay above a given magnitude limit (e.g., red dashed line for a DES-like survey), in turn introducing greater scatter in their derived properties (Figure 3).</figcaption></figure>

create strong phase-space artifacts, especially when stellar sampling is limited (see, e.g., the discussion in S. H. Lim et al. 2024). The vastly improved stellar sampling of our resolved-star simulations—≈100,000 stars per galaxy compared to ≈500 stellar particles in M. P. Rey et al. (2019)—significantly improves the robustness in predicting the kinematic and chemical properties of the simulated galaxy. We turn next to quantifying how this impacts the observational recovery of structural parameters and will explore radial trends and properties of distant members in future work.

## 4.2. Structural Properties

We derive the half-light radius, total magnitude, and average iron metallicity ($r_h$, $M_V$, and $\langle[\text{Fe/H}]\rangle$, respectively) as an observer would. Using the projected positions ($\alpha_i, \beta_i$) for each star $i$ brighter than the assumed limit ($r=23$ mag), we fit exponential profiles to the projected stellar number density (see N. F. Martin et al. 2008; N. F. Martin et al. 2016 and Appendix B for further details). Red ellipses in Figure 2 show the best-fit elliptical half-light radius $r_h$ for each galaxy, with their stellar density profiles and best-fit parameters shown in Figure 4. $\langle[\text{Fe/H}]\rangle$ is the average of stars above the magnitude limit.

Figure 3 shows the derived $M_V$, $r_h$, and $\langle[\text{Fe/H}]\rangle$ with these procedures (colored crosses) compared to observational data. We also compare to “simulated” values (colored stars), determined by counting all stars within $r_{200}$ for $M_V$ and $[\text{Fe/H}]$, and interpolating the spherical half-light radius for $r_h$. These “simulated” definitions are intrinsically different from the observational quantities, but are commonly used across the literature to compare models to data (e.g., C. Wheeler et al. 2019; O. Agertz et al. 2020; M. P. Rey et al. 2020) and provide a complete account of each galaxy’s stars. For all galaxies, the total magnitudes are underestimated using the observational approach,

by up to 0.3 mag for our brightest objects. This should be expected because a magnitude-limited observation will miss stars available to the simulation. We emphasize that the lower magnitude is not due to missing the light from stars below the observed limit, which is accounted for by sampling luminosity functions (see Appendix B). Instead, the missing light is from stars not accounted for by the fitted exponential profile. Furthermore, brighter dwarfs (with more observed stars) remain biased in $M_V$ and also show less accurate $r_h$ determinations (up to ≈1.5 × offset).

All our galaxies exhibit surface density profiles with multiple components including a roughly exponential center and an additional extended halo (see Figure 5 and Appendix B for details). When fitted to a single exponential, the extended population is not accounted for, biasing $M_V$ to fainter values and $r_h$ to more compact values. ΛCDM structure formation predicts that even faint dwarf galaxies exhibit extended halos and distant members (Figure 2; see also, e.g., Y. Tarumi et al. 2021; A. Goater et al. 2024), so this inadequacy should be expected. Furthermore, the physical extended halo of the dwarf galaxy is largely indistinguishable from the typical background assumed for Milky Way observations ($\Sigma \approx 10^{-1} \text{ arcmin}^{-2}$; see Equation (B1); N. F. Martin et al. 2008, 2016). This gives the appearance of a good fit of the stellar population on the sky, but misses physical members of the dwarf galaxy (see Figure 5) resulting in further systematic offsets in $M_V$ compared to the values estimated directly from the simulations. Distinguishing background stars from members will remain a challenge for real observations in the Local Group, but our simulations provide a first estimate of the biases when comparing models and observations. Furthermore, our simulated dwarfs evolve in the field and are therefore less prone to tidal stripping compared with the satellite galaxies observed in the Local Group. We
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Half-light radius (top panel) and average iron abundance (bottom panel) as a function of total absolute V-band magnitude of our simulations, compared to Local Group data (gray). Obtaining $r_h$ and $M_V$ via mock observation of resolved stars (crosses) provides estimates of the total V-band magnitude that are systematically lower than the magnitude derived from the simulations (star symbols). This apparent dimming arises from light profiles emerging from cosmological formation histories that are more complex than typically assumed for observed dwarfs. Stochastically resampling the CMD for the same history (contours around crosses) introduces a noise term, as different realizations of the light and metallicity distributions get over the magnitude limit. This noise becomes comparable to the population scatter for our fainter objects, with only a few tens of observable stars. The observed data are from the compilations of A. W. McConnachie (2012), E. N. Kirby et al. (2013, 2014), and J. D. Simon (2019) augmented with individual candidates and detections from G. Torrealba et al. (2016, 2018, 2019), D. Homma et al. (2018, 2024), S. Mau et al. (2020), P. Bennet et al. (2022), H. Richstein et al. (2022), D. J. Sand et al. (2022), W. Cerny et al. (2023b), W. Cerny et al. (2023a, 2023c), M. G. Jones et al. (2023, 2024), K. B. W. McQuinn et al. (2023), M. L. M. Collins et al. (2024), J. Li et al. (2024), D. Martinez-Delgado et al. (2024), M. McNanna et al. (2024), S. E. T. Smith et al. (2023, 2024), and C. Y. Tan et al. (2024).</figcaption></figure>

leave testing how our galaxies evolve in an external potential field to future work.

Another source of uncertainty when estimating dwarf galaxy properties comes from the inherent stochasticity associated with which stars reside above the magnitude limit at $z=0$, and are thus bright enough to be observed (CMD shot noise). Observed stars only provide a subset of the true spatial distribution, light profile, and metallicity distribution of the whole dwarf galaxy.

The exact CMD realization will depend on the initial mass of each low-mass star, its precise metallicity, and age. All these quantities are deterministic in our simulations, providing us with one example realization of the CMD for a given dwarf galaxy’s history.

Nonetheless, we can estimate the magnitude of this noise term by perturbing the initial mass of each of our stars. In other words, the true spatial, age, and metallicity distributions are unchanged, but we resample which stars ascend the giant branch. We achieve this by resampling the IMF 100 times, in very fine bins (0.25 $M_\odot$) to minimize the impact that a different mass has on the dynamical evolution of any given star throughout the simulation. We then refit all properties as previously and show the 68%-95% confidence intervals as contours in Figure 3.

This noise has limited impact on bright galaxies with a well-sampled CMD (purple and violet), where it is far outweighed by the inaccurate $r_h$ and $M_V$ from fitting. For fainter galaxies (cyan and blue) with only a few tens of visible stars, this shot noise introduce significant uncertainty, with $2\sigma$ shifts of up to 0.4 mag, 230 pc, and 0.3 dex for $M_V$, $r_h$, and $\langle[\text{Fe/H}]\rangle$, respectively.

These numbers are highly significant to the interpretation of dwarf galaxy data—they are comparable to the population scatter in the $\langle[\text{Fe/H}]\rangle - M_V$ plane when galaxies only have a few tens of observable stars. The number of detected member stars, and thus the extent of this noise, will depend on the specifics of each observation. But, even with modern facilities, faint dwarf galaxies often have <50 candidate member stars (e.g., W. Cerny et al. 2023b; D. Homma et al. 2024; S. E. T. Smith et al. 2024), hinting that some of the observed scatter could be explained through CMD shot noise. Furthermore, this effect will be compounded by the need for spectroscopic follow-up. Here, all observable stars above a magnitude limit are considered when estimating a mean, while in reality, only a subset of stars can be followed up with high-resolution spectroscopy to determine [Fe/H] and metal abundances. We emphasize that when the metallicity distribution is sampled by a few tens of stars, a mean value for, e.g., iron abundance is a poor constraint on models. We provide a detailed analysis of specific abundance features, their observability, and the uncertainty of varying element yields in E. P. Andersson et al. (2024, in preparation).

## 5. Discussion and Conclusion

We present a novel set of cosmological simulations of small dwarf galaxies with unprecedented detail, tracking and evolving all stars observable throughout the Hubble time on a star-by-star basis. With this, we self-consistently predict how $z=0$ resolved-star observables (e.g., observable spatial distribution and CMD; Figure 2) arise from the cosmological mass-growth, star formation, and chemical enrichment history of each dwarf galaxy. This allows us to directly apply observational parameter estimation to test the accuracy and precision of these pipelines on realistic objects (Figure 3).

In all four of dwarfs, hierarchical assembly creates light profiles with at least two different components. This leads to systematic offsets when recovering structural parameters with methods assuming a single exponential for profiles (e.g., N. F. Martin et al. 2016), with the total V-band magnitude being shifted by up to 0.5 mag compared to the value derived from simulations with commonly used methods. This outer
component is faint, indistinguishable from the background, and likely affected by tidal stripping in galaxies that reside in group environments. In future work, we will investigate how this can be accounted for and if methods including different density models (e.g., H. C. Plummer 1911; J. L. Sérsic 1963; I. R. King 1966) including multiple components can give better estimates of the light profiles for observed dwarfs.

In realistic cases of observed faint dwarf galaxies, few member stars are bright enough to be confidently assigned to the galaxy. This introduces an intrinsic uncertainty to recovering the dwarf's galaxy underlying light profile and metallicity distribution from observationally estimated properties. We show that when $\le 100$ members stars are observable, CMD shot noise introduces significant uncertainty. This uncertainty makes interpretation particularly challenging in the $\langle [Fe/H] \rangle - M_V$ plane, where it is comparable to the population scatter.

This work marks a significant step toward bridging simulations and observations of dwarf galaxies. The differences in quantities estimated using techniques commonly applied to survey pipelines compared to those commonly applied to simulations uncovered in this work have direct implications on how we map observables such as $M_V$ (and thus $M_*$) and $\langle [Fe/H] \rangle$ onto cosmological dark matter halos, and in turn constrain the nature of dark matter (e.g., E. O. Nadler et al. 2019) or the strength of galactic outflows (e.g., O. Agertz et al. 2020) with faint dwarf galaxies. While the magnitude of these differences is left to be quantified across a wide range of surveys for galaxies at all distances, this work underscores the importance of self-consistently accounting for the full cosmological history rather than assuming idealized objects when testing observational pipelines (but see, e.g., B. Mutlu-Pakdil et al. 2021). The advent of individual-star models such as INFERNO, and their application to cosmological simulations of galaxies down to low-mass stars, will enable their quantification in the near future and unlock much more robust interpretation of forthcoming data from, e.g., the LSST.

## Acknowledgments

We thank the anonymous referee for comments and suggestions that greatly improved the quality of this work. We acknowledge Mordecai-Mark Mac Low for helpful comments improving the quality of this manuscript.

E.A. acknowledges support from NASA ATP grant 80NSSC24K0935 and NSF grant AST23-07950. M.R. is supported by the Beecroft Fellowship funded by Adrian Beecroft. O.A. acknowledges support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019-04659), and the Swedish National Space Agency (SNSA Dnr 2023-00164). A.P. was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 818085 GMGalaxies). J.I.R. would like to acknowledge support from STFC grants ST/Y002865/1 and ST/Y002857/1.

This work used the DiRAC Data Intensive service (DIaL2/ DiLaL) at the University of Leicester, managed by the University of Leicester Research Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Leicester was funded by BEIS, UKRI, and STFC capital funding and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure.

We thank the Astrophysics Data Service funded by NASA under Cooperative Agreement 80NSSC21M00561 and the

arXiv preprint repository for providing services that were used extensively in this work.

Software: GENETIC (S. Stopyra et al. 2021), RAMSES (R. Teyssier 2002), PYNBODY (A. Pontzen et al. 2013, 2022), TANGOS (A. Pontzen & M. Tremmel 2018), NUMPY (S. van der Walt et al. 2011; C. R. Harris et al. 2020), SCIPY (P. Virtanen et al. 2020), EMCEE (D. Foreman-Mackey et al. 2013), MATPLOTLIB (J. D. Hunter 2007), and JUPYTER (M. Ragan-Kelley et al. 2014).

## Appendix A: Isochrones and Luminosity Functions

We use isochrones and luminosity function derived from the PARSEC web-interface<fnref n="10" /> applying PARSEC version 1.2S (Y. Chen et al. 2014, 2015; J. Tang et al. 2014) with COLIBRI version S_37 for thermally pulsing AGB end-phase evolution (P. Marigo et al. 2013; P. Rosenfield et al. 2016; G. Pastorelli et al. 2019, 2020). Circumstellar dust is included using scaling relations from P. Marigo et al. (2008) with 60% silicate and 40% AlO for M stars and 85% AMC and 15% SiC for C stars (M. A. T. Groenewegen 2006). Long-period variability during the red giant branch and AGB is included following M. Trabucchi et al. (2021). For magnitudes, we used photometric systems UVRIJHK (M. S. Bessell 1990; J. Maïz Apellániz 2006) and ugriz from SDSS. We do not consider dust extinction for any of our results.

## Appendix B: Method for Determining Galaxy Luminosity and Structural Properties

We determine the structural properties of the galaxies by the technique presented in N. F. Martin et al. (2008; see also N. F. Martin et al. 2016). Given a distance modulus and magnitude limit, the observable stars $i$ in each galaxy are projected to sky coordinates $(\alpha_i, \beta_i)$. An additional set of stars is added to this group by uniformly sampling a background stellar field around the galaxy. For all examples presented throughout this work we use a background of 0.1 stars arcmin⁻². While the number density of the background is somewhat arbitrary, it is of the order expected in the Local Group and is expected by the fitting model. As in N. F. Martin et al. (2008), the projected number density in units of stars arcmin⁻² is given by

$$n(r) = \frac{1.68^2 N_* \exp(-1.68r/r_h)}{\underbrace{2\pi(1-e)r_h^2}_{\text{galaxy}} + \frac{(N_{\text{tot}} - N_*)/(3600\pi)}{\underbrace{\text{background}}_{\text{background}}}} \quad (B1)$$

where $N_*$ is the number of stars in the galaxy, $e$ is the eccentricity, $r_h$ is the half-density radius in units of arcminutes, $N_{\text{tot}}$ is the total number of stars in the sample, and $r$ is the elliptical radius centered on the galaxy coordinates $(\alpha_0, \beta_0)$ rotated by a north-to-east angle $\theta$ such that the major axis aligns with the position axis. We define ellipticity as one minus the ratio between the scale length along the minor and the major axes. The fitting procedure use applies the Markov Chain

<footnote n="10">¹⁰ http://stev.oapd.inaf.it/cmd</footnote>
Monte Carlo sampler (J. Goodman & J. Weare 2010) implemented in EMCEE (D. Foreman-Mackey et al. 2013) using 64 walkers and 5000 steps to determine the maximum of the log-likelihood function $\mathcal{L}(\alpha_0, \beta_0, \theta, e, r_h, N_*) = \prod_i n(r_i)$.

To determine the total magnitude, we use luminosity functions (see Appendix A) scaled to the assumed distance and sample stars from an IMF (P. Kroupa 2001) until the number of stars with a magnitude brighter than the set limit reaches the best fit for $N_*$. In this case, the luminosity function is derived assuming the age and metallicity of those of the mean of the galaxy. This provides a complete stellar sample used to calculate the total magnitude in a given band. To ensure robustness of these estimates we perform

this exercise 1000 times for each galaxy and compute the mean and standard deviation. We note that in all cases, this standard deviation is negligible compared to other errors.

## Appendix C: Estimating Structure Parameters of Simulated Data

Figure 4 shows a typical example of one fitting procedure for the DES magnitude limit ($r=23$) at 400 kpc used throughout the main text. Note that this is only one out of 100 fits used to determine the distributions presented in Figure 3. The figure includes the distribution of stars (top panel, galaxy members in red points), the observed stellar number density (black points

<figure><img src="image_8.png" /><figcaption>**Figure 4.** The projected surface number density of stars brighter than 23 mag in $r$ band at 400 kpc as a function of projected distance from the galaxy center (black points with error bars indicating Poisson noise). Red lines show profiles (Equation (B1)) fitted by a maximum likelihood search using EMCEE (D. Foreman-Mackey et al. 2013) with parameters and uncertainties from the fitting procedure written for each galaxy in each panel. ($\alpha_0, \beta_0$) are the offsets from the halo center (determined through a shrinking sphere applied to dark matter particles) in arcminutes, $\theta$ is the angle in degrees between the major axis and the projection of the galaxy in the simulation $x-y$ plane, $e$ is the ellipticity defined as one minus the ratio between major and minor axis length, $r_h$ is the half-density/mass-radius, and $N_*$ is the number of stars the fitting procedure assigned as galaxy members. The true number of stars is indicated in the top-left corner of each panel.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The projected surface number density of stars in our simulations. The gray dashed lines denote where the profile visually deviates from an exponential (extended halo stars). Red and blue lines show the exponential profile (Equation (B1)) without the second term) using parameters obtained from the fitting procedure (Appendix B) with and without a population of background stars (see text for details), respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6. Same as Figure 4 but for 22 mag in r band at 200 kpc.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7. Same as Figure 6 but for 24 mag in r band at 1000 kpc.</figcaption></figure>

including Poisson noise), the resulting profile (red line), and
parameters with errors from the fitting procedure (bottom panel).

Figure 5 shows the surface number density of each simulation
calculated using spherical radii. The gray dashed lines denote the
radius where we find a transition into an extended diffuse
halo through visual inspection, highlighting multiple structural

components in our galaxies. This transition results from galaxy
growth in a cosmological environment (in particular through
mergers) and makes fitting the profile in Equation (B1) challenging.
This extended halo cannot be distinguished when we artificially add
a sample of background stars, but affects the parameters derived
through the fitting method.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Each column is same as Figure 3, showcasing *r*-band limiting magnitudes at 22 (200 kpc), 23 (400 kpc), and 24 (1 Mpc), i.e., typical point-source magnitude limits for SDSS, DES, and the LSST, respectively, at distances where the smallest galaxy ("GM3: Latest") is at the limit of observability. Note that the center column is identical to Figure 3.</figcaption></figure>

Throughout the main text, we presented our results as if observed at 400 kpc with a limiting magnitude in *r* band of 23. In this way, our estimates suffer from the same limitations as if observed with state-of-the-art facilities and therefore directly comparable with the data shown by error bars in Figure 3 (J. D. Simon 2019). In addition, we performed the fitting procedure using SDSS limits (*r*-band limit of 22 mag, Figure 6) and the upcoming LSST (single pointing with *r*-band limit of 24 mag, Figure 7). Similarly to the main text, we place the galaxies at a distance of 200 kpc (1 Mpc) where the smallest galaxy is barely observable with SDSS (LSST) to provide constraints for a range of observational quality.

For completeness, Figure 8 presents the half-light radius and average [Fe/H] as functions of total *V*-band magnitude for all three facilities. Our conclusions remain the same for all three cases, only showing an increase in the scatter in the estimated observables directly related to the number of observed stars.

## ORCID iDs

Eric P. Andersson https://orcid.org/0000-0003-3479-4606

Martin P. Rey https://orcid.org/0000-0002-1515-995X

Andrew Pontzen https://orcid.org/0000-0001-9546-3849

Corentin Cadiou https://orcid.org/0000-0003-2285-0332

Oscar Agertz https://orcid.org/0000-0002-4287-1088

Justin I. Read https://orcid.org/0000-0002-1164-9302

Nicolas F. Martin https://orcid.org/0000-0002-1349-202X

## References

* Agertz, O., Pontzen, A., Read, J. I., et al. 2020, MNRAS, 491, 1656

* Andersson, E. P., Agertz, O., & Renaud, F. 2020, MNRAS, 494, 3328

* Andersson, E. P., Agertz, O., Renaud, F., & Teyssier, R. 2023, MNRAS, 521, 2196

* Andersson, E. P., Mac Low, M.-M., Agertz, O., Renaud, F., & Li, H. 2024, A&A, 681, A28

* Andersson, E. P., Renaud, F., & Agertz, O. 2021, MNRAS, 502, L29

* Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007, ApJ, 654, 897

* Benitez-Llambay, A., & Fumagalli, M. 2021, ApJL, 921, L9

* Bennet, P., Sand, D. J., Crnojević, D., et al. 2022, ApJ, 924, 98

* Bessell, M. S. 1990, PASP, 102, 1181

* Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

* Calura, F., Lupi, A., Rosdahl, J., et al. 2022, MNRAS, 516, 5914

* Cerny, W., Drlica-Wagner, A., Li, T. S., et al. 2023a, ApJL, 953, L21

* Cerny, W., Martínez-Vázquez, C. E., Drlica-Wagner, A., et al. 2023b, ApJ, 953, 1

* Cerny, W., Simon, J. D., Li, T. S., et al. 2023c, ApJ, 942, 111

* Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068

* Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525

* Chiti, A., Frebel, A., Simon, J. D., et al. 2021, NatAs, 5, 392

* Collins, M. L. M., Karim, N., Martínez-Delgado, D., et al. 2024, MNRAS, 528, 2614

* Cummings, J. D., Kalirai, J. S., Tremblay, P. E., & Ramirez-Ruiz, E. 2016, ApJ, 818, 84

* Deng, Y., Li, H., Liu, B., et al. 2024, A&A, 691, A231

* Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al. 2015, ApJ, 813, 109

* Emerick, A., Bryan, G. L., & Mac Low, M.-M. 2018, ApJL, 865, L22

* Eriksson, K., Nowotny, W., Höfner, S., Aringer, B., & Wachter, A. 2014, A&A, 566, A95

* Fitts, A., Boylan-Kolchin, M., Elbert, O. D., et al. 2017, MNRAS, 471, 3547

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Frebel, A., & Norris, J. E. 2015, ARA&A, 53, 631

* Goater, A., Read, J. I., Noël, N. E. D., et al. 2024, MNRAS, 527, 2403

* Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65

* Grand, R. J. J., Helly, J., Fattahi, A., et al. 2018, MNRAS, 481, 1726

* Gray, E. I., Read, J. I., Taylor, E., et al. 2024, arXiv:2405.19286

* Groenewegen, M. A. T. 2006, A&A, 448, 181

* Grudić, M. Y., Guszejnov, D., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2021, MNRAS, 506, 2199

* Grudić, M. Y., Guszejnov, D., Offner, S. S. R., et al. 2022, MNRAS, 512, 216

* Gutcke, T. A., Pakmor, R., Naab, T., & Springel, V. 2021, MNRAS, 501, 5597

* Gutcke, T. A., Pakmor, R., Naab, T., & Springel, V. 2022a, MNRAS, 513, 1372

* Gutcke, T. A., Pfrommer, C., Bryan, G. L., et al. 2022b, ApJ, 941, 120

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

* Hislop, J. M., Naab, T., Steinwandel, U. P., et al. 2022, MNRAS, 509, 5938
Höfner, S., & Olofsson, H. 2018, A&ARv, 26, 1

Homma, D., Chiba, M., Komiyama, Y., et al. 2024, PASJ, 76, 733

Homma, D., Chiba, M., Okamoto, S., et al. 2018, PASJ, 70, S18

Hu, C.-Y., Naab, T., Glover, S. C. O., Walch, S., & Clark, P. C. 2017, MNRAS, 471, 2151

Hu, C.-Y., Naab, T., Walch, S., Glover, S. C. O., & Clark, P. C. 2016, MNRAS, 458, 3528

Hunter, J. D. 2007, CSE, 9, 90

Jethwa, P., Erkal, D., & Belokurov, V. 2018, MNRAS, 473, 2060

Ji, A. P., Frebel, A., & Bromm, V. 2015, MNRAS, 454, 659

Jones, M. G., Mutlu-Pakdil, B., Sand, D. J., et al. 2023, ApJL, 957, L5

Jones, M. G., Sand, D. J., Mutlu-Pakdil, B., et al. 2024, ApJL, 971, L37

Kim, C.-G., & Ostriker, E. C. 2015, ApJ, 802, 99

King, I. R. 1966, AJ, 71, 64

Kirby, E. N., Bullock, J. S., Boylan-Kolchin, M., Kaplinghat, M., & Cohen, J. G. 2014, MNRAS, 439, 1015

Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102

Koposov, S. E., Belokurov, V., Torrealba, G., & Evans, N. W. 2015, ApJ, 805, 130

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., Fumagalli, M., da Silva, R. L., Rendahl, T., & Parra, J. 2015, MNRAS, 452, 1447

Lahén, N., Naab, T., Johansson, P. H., et al. 2020, ApJ, 891, 2

Lahén, N., Naab, T., & Szécsi, D. 2024a, MNRAS, 530, 645

Lahén, N., Rantala, A., Naab, T., et al. 2024b, arXiv:2410.01891

Lahén, N., Naab, T., Kauffmann, G., et al. 2023, MNRAS, 522, 3092

Li, J., Greene, J. E., Carlsten, S. G., & Danieli, S. 2024, ApJ, 975, L23

Lim, S. H., Raman, K. A., Buckley, M. R., & Shih, D. 2024, MNRAS, 533, 143

Maíz Apellániz, J. 2006, AJ, 131, 1184

Maoz, D., & Graur, O. 2017, ApJ, 848, 25

Marigo, P., Bressan, A., Nanni, A., Girardi, L., & Pumo, M. L. 2013, MNRAS, 434, 488

Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883

Martin, N. F., de Jong, J. T. A., Rix, H. W., et al. 2008, ApJ, 684, 1075

Martin, N. F., Ibata, R. A., Lewis, G. F., et al. 2016, ApJ, 833, 167

Martinez-Delgado, D., Stein, M., Pawlowski, M. S., et al. 2024, arXiv:2405.03769

Mau, S., Cerny, W., Pace, A. B., et al. 2020, ApJ, 890, 136

McConnachie, A. W. 2012, AJ, 144, 4

McNanna, M., Bechtol, K., Mau, S., et al. 2024, ApJ, 961, 126

McQuinn, K. B. W., Mao, Y.-Y., Buckley, M. R., et al. 2023, ApJ, 944, 14

Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603

Mutlu-Pakdil, B., Sand, D. J., Crnojević, D., et al. 2021, ApJ, 918, 88

Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59

Nadler, E. O., Gluscevic, V., Boddy, K. K., & Wechsler, R. H. 2019, ApJL, 878, L32

Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020, ApJ, 893, 48

Newton, O., Cautun, M., Jenkins, A., Frenk, C. S., & Helly, J. C. 2018, MNRAS, 479, 2853

Ohlin, L., Renaud, F., & Agertz, O. 2019, MNRAS, 485, 3887

Pastorelli, G., Marigo, P., Girardi, L., et al. 2019, MNRAS, 485, 5666

Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283

Pignatari, M., Herwig, F., Hirschi, R., et al. 2016, ApJS, 225, 24

Plummer, H. C. 1911, MNRAS, 71, 460

Polak, B., Mac Low, M.-M., Klessen, R. S., et al. 2024, A&A, 690, A207

Pontzen, A., Cadiou, C., Stopyra, S., et al. 2022, synbody/genetIC: v1.5.0, Zenodo, doi:10.5281/zenodo.3874396

Pontzen, A., Roškar, R., Stinson, G., & Woods, R. 2013, synbody: N-Body/SPH Analysis for Python, Astrophysics Source Code Library, ascl:1305.002

Pontzen, A., & Tremmel, M. 2018, ApJS, 237, 23

Ragan-Kelley, M., Perez, F., Granger, B., et al. 2014, AGU Fall Meeting, 2014, H44D

Raiteri, C. M., Villata, M., & Navarro, J. F. 1996, A&A, 315, 105

Rey, M. P., & Pontzen, A. 2018, MNRAS, 474, 45

Rey, M. P., Pontzen, A., Agertz, O., et al. 2019, ApJL, 886, L3

Rey, M. P., Pontzen, A., Agertz, O., et al. 2020, MNRAS, 497, 1508

Rey, M. P., Pontzen, A., Agertz, O., et al. 2022, MNRAS, 511, 5672

Richstein, H., Patel, E., Kallivayalil, N., et al. 2022, ApJ, 933, 217

Ritter, C., Herwig, F., Jones, S., et al. 2018, MNRAS, 480, 538

Rosenfield, P., Marigo, P., Girardi, L., et al. 2016, ApJ, 822, 73

Roth, N., Pontzen, A., & Peiris, H. V. 2016, MNRAS, 455, 974

Sales, L. V., Wetzel, A., & Fattahi, A. 2022, NatAs, 6, 897

Sand, D. J., Mutlu-Pakdil, B., Jones, M. G., et al. 2022, ApJL, 935, L17

Sanderson, R. E., Wetzel, A., Loebman, S., et al. 2020, ApJS, 246, 6

Seitenzahl, I. R., Ciaraldi-Schoolmann, F., Röpke, F. K., et al. 2013, MNRAS, 429, 1156

Sérsic, J. L. 1963, BAAA, 6, 41

Sestito, F., Zaremba, D., Venn, K. A., et al. 2023, MNRAS, 525, 2875

Simon, J. D. 2019, ARA&A, 57, 375

Smith, M. C. 2021, MNRAS, 502, 5417

Smith, S. E. T., Cerny, W., Hayes, C. R., et al. 2024, ApJ, 961, 92

Smith, S. E. T., Jensen, J., Roediger, J., et al. 2023, AJ, 166, 76

Sormani, M. C., Treß, R. G., Klessen, R. S., & Glover, S. C. O. 2017, MNRAS, 466, 407

Steinwandel, U. P., Bryan, G. L., Somerville, R. S., Hayward, C. C., & Burkhart, B. 2023, MNRAS, 526, 1408

Steinwandel, U. P., Kim, C.-G., Bryan, G. L., et al. 2024, ApJ, 960, 100

Stopyra, S., Pontzen, A., Peiris, H., Roth, N., & Rey, M. P. 2021, ApJS, 252, 28

Tan, C. Y., Cerny, W., Drlica-Wagner, A., et al. 2024, arXiv:2408.00865

Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287

Tarumi, Y., Yoshida, N., & Frebel, A. 2021, ApJL, 914, L10

Teyssier, R. 2002, A&A, 385, 337

Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016, MNRAS, 459, 2370

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2018, MNRAS, 475, 5085

Torrealba, G., Belokurov, V., Koposov, S. E., et al. 2019, MNRAS, 488, 2743

Trabucchi, M., Wood, P. R., Mowlavi, N., et al. 2021, MNRAS, 500, 1575

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Vink, J. S. 2011, Ap&SS, 336, 163

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wheeler, C., Hopkins, P. F., Pace, A. B., et al. 2019, MNRAS, 490, 4447