# EDGE: The emergence of dwarf galaxy scaling relations from cosmological radiation-hydrodynamics simulations

Martin P Rey, Ethan Taylor, Emily I Gray, Stacy Y Kim, Eric P Andersson, Andrew Pontzen, Oscar Agertz, Justin I Read, Corentin Cadiou, Robert M Yates, et al.

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

## To cite this version:

Martin P Rey, Ethan Taylor, Emily I Gray, Stacy Y Kim, Eric P Andersson, et al. EDGE: The emergence of dwarf galaxy scaling relations from cosmological radiation-hydrodynamics simulations. Monthly Notices of the Royal Astronomical Society, 2025, 541 (2), pp.1195-1217. 10.1093/mn-ras/staf1058. hal-05005200

HAL Id: hal-05005200

https://hal.science/hal-05005200v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>
# EDGE: the emergence of dwarf galaxy scaling relations from cosmological radiation-hydrodynamics simulations

Martin P. Rey<sup>1,2*</sup> Ethan Taylor,<sup>3</sup> Emily I. Gray<sup>3</sup>, Stacy Y. Kim,<sup>4</sup> Eric P. Andersson<sup>5</sup>

Andrew Pontzen,<sup>6</sup> Oscar Agertz<sup>6</sup>,<sup>7</sup> Justin I. Read<sup>6</sup>,<sup>3</sup> Corentin Cadiou<sup>7,8</sup>, Robert M. Yates<sup>6</sup>,<sup>9</sup>

Matthew D. A. Orkney<sup>10,11</sup> Dirk Scholte,<sup>12</sup> Amélie Saintonge<sup>13,14</sup>, Joseph Breneman,<sup>15</sup>

Kristen B. W. McQuinn,<sup>15,16</sup> Claudia Muni<sup>13</sup> and Payel Das<sup>3</sup>

Accepted 2025 June 26. Received 2025 June 25; in original form 2025 March 13

## ABSTRACT

We present a new suite of EDGE ('Engineering Dwarfs at Galaxy formation's Edge') cosmological zoom simulations. The suite includes 15 radiation-hydrodynamical dwarf galaxies covering the ultrafaint to the dwarf irregular regime ($10^4 \le M_*(z=0) \le 10^8 M_\odot$) to enable comparisons with observed scaling relations. Each object in the suite is evolved at high resolution ($\approx 3$ pc) and includes stellar radiation, winds, and supernova feedback channels. We compare with previous EDGE simulations without radiation, finding that radiative feedback results in significantly weaker galactic outflows. This generalizes our previous findings to a wide mass range, and reveals that the effect is most significant at low $M_\star$. Despite this difference, stellar masses stay within a factor of two of each other, and key scaling relations of dwarf galaxies (size–mass, neutral gas–stellar mass, and gas-phase mass–metallicity) emerge correctly in both simulation suites. Only the stellar mass–stellar metallicity relation is strongly sensitive to the change in feedback. This highlights how obtaining statistical samples of dwarf galaxy stellar abundances with next-generation spectrographs will be key to probing and constraining the baryon cycle of dwarf galaxies.

**Key words:** methods: numerical – galaxies: dwarf – galaxies: evolution – galaxies: ISM – galaxies: structure.

## 1 INTRODUCTION

Small 'dwarf' galaxies are sensitive probes of galaxy formation and dark matter physics. Within their shallow gravitational potential wells, energetic stellar processes [e.g. supernovae (SNe), stellar winds and radiation, and collectively, 'feedback'] efficiently drive galactic outflows (see Collins & Read 2022 for a review), making dwarf galaxies an ideal laboratory to study how feedback processes regulate the growth of galaxies over cosmic time (e.g. Naab &

Ostriker 2017 for a review). Furthermore, the existence of their low-mass dark matter haloes offer leading constraints on the 'coldness' of dark matter (e.g. Nadler et al. 2020), while their internal dynamics directly probe if dark matter is self-interacting (see Pontzen & Governato 2014; Bullock & Boylan-Kolchin 2017; Sales, Wetzel & Fattahi 2022 for reviews).

The same low mass that makes dwarf galaxies sensitive to important physical processes has historically also made them challenging to observe and characterize. But the advent of wide-field, deep photometric surveys has now revealed an ever-growing number of classical and ultrafaint dwarf galaxies (Simon 2019 for a review). This has in turn enabled dedicated programs to study the stellar

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

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

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

<footnote n="4"><sup>4</sup>Carnegie Theoretical Astrophysics Center, Carnegie Observatories, 813 Santa Barbara St, Pasadena, CA 91106, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Astrophysics, American Museum of Natural History, 200 Central Park West, New York, NY 10024, USA</footnote>

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

<footnote n="7"><sup>7</sup>Department of Physics, Lund Observatory, Division of Astrophysics, Lund University, Box 43, SE-22100 Lund, Sweden</footnote>

<footnote n="8"><sup>8</sup>Institut d'Astrophysique de Paris, Sorbonne Universités, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="9"><sup>9</sup>Centre for Astrophysics Research, University of Hertfordshire, Hatfield AL10 9AB, UK</footnote>

<footnote n="10"><sup>10</sup>Institut de Ciencies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain</footnote>

<footnote n="11"><sup>11</sup>Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="12"><sup>12</sup>Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="13"><sup>13</sup>Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK</footnote>

<footnote n="14"><sup>14</sup>Max-Planck-Institut für Radioastronomie (MPIfR), Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="15"><sup>15</sup>Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Rd, Piscataway, NJ 08854, USA</footnote>

<footnote n="16"><sup>16</sup>Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote>*E-mail: mpr47@bath.ac.uk</footnote>
chemistry and kinematics of these systems (e.g. Ji et al. 2016a, b, c, 2020; Longeard et al. 2018; Buttry et al. 2022; Simon et al. 2023; Bruce et al. 2023; Hansen et al. 2024). Furthermore, the new generation of imaging can now be combined with upgrades in radio capabilities (e.g. with the MeerKAT Telescope) and multi-object spectroscopy (e.g. with the Dark Energy Spectroscopic instrument; DESI). This has fundamentally changed our ability to characterize the gas contents and gas-phase metallicities of faint dwarf galaxies, pushing characterization to the very faintest end (e.g. McQuinn et al. 2021) and vastly extending statistical samples (e.g. Scholte et al. 2025).

This forward trend in observational capabilities and discoveries will continue in the next 5 yr, spearheaded by experiments such as the Vera C. Rubin Observatory which is expected to provide a near-complete census of faint ($M_V \approx -6$) dwarf galaxies within 5 Mpc around the Milky Way (Mutlu-Pakdil et al. 2021) and the *Nancy Roman Space Telescope*. These new photometric data sets will be complemented by forthcoming spectroscopic information, with for example the 4MOST-Dwarf survey expected to obtain chemical abundances for $> 50\,000$ dwarf galaxy stars (Skúladóttir et al. 2023).

Interpreting these new data sets requires us to develop detailed models of dwarf galaxy formation that can match the data's new statistical power while retaining enough physical fidelity to make robust predictions. This is a challenging task. The small sizes of dwarf galaxies ($\approx 100$ pc) require high ($\approx 10$ pc) numerical resolution to resolve their gas reservoirs and star-forming regions over cosmic time. But, at the same time, modelling the full cosmological history of a dwarf galaxy is also essential, as each specific formation scenario plays a key role in setting the $z=0$ properties and observables (e.g. Benítez-Llambay et al. 2015; Benitez-Llambay & Fumagalli 2021; Fitts et al. 2017; Rey et al. 2019b, 2020; Wright et al. 2019; Katz et al. 2020; Tarumi, Yoshida & Frebel 2021; Herzog, Benítez-Llambay & Fumagalli 2023). Simulating many objects, in a cosmological context, and with a high-resolution place conflicting demands on the available computing power.

None the less, the field has made great progress in addressing these challenges in recent years. Improvements in code efficiency and computational power now allow zoomed cosmological simulations of field faint dwarf galaxies with $\approx$ pc resolution over the full *Hubble* time (e.g. Wheeler et al. 2019; Agertz et al. 2020; Gutcke et al. 2022; Go et al. 2025). These gains go beyond an incremental improvement in resolution, as they allow us to resolve explicitly the cooling radius of SNe explosions, in turn enabling a robust modelling of the emerging momentum and its coupling to the gas (e.g. Kimm et al. 2015; Kim & Ostriker 2015; Martizzi, Faucher-Giguère & Quataert 2015; Ohlin, Renaud & Agertz 2019).

In parallel, model improvements are significantly increasing physical fidelity. For example, radiation-hydrodynamics simulations can now account for stellar radiative heating over the cosmological history of a dwarf galaxy (Agertz et al. 2020). Unlike previous subgrid implementation that were attempted to capture radiative effects locally (e.g. Agertz et al. 2013; Hopkins et al. 2018), explicit radiative transfer gives a physical account of how photons ionize, heat, and inject momentum in the gas both locally, and non-locally across the interstellar medium (ISM) and circumgalactic medium (CGM). Beyond radiative transfer, cosmological simulations of dwarf galaxies have also made strides to explicitly sample the initial mass function (IMF) with individual stars (e.g. Gutcke et al. 2022; Andersson et al. 2025), unlocking direct comparisons with resolved-star observations of dwarf galaxies in the Local Volume. And high-resolution cosmological simulations have also started to incorporate detailed models of cosmic ray (e.g. Martin-Alvarez et al. 2023)

and black hole feedback (e.g. Koudmani, Henden & Sijacki 2021; Koudmani, Sijacki & Smith 2022; Arjona-Galvez, Di Cintio & Grand 2024), starting to quantify the importance of these processes at the faint end of galaxy formation at $z=0$.

Furthermore, lessons learned from these detailed numerical simulations directly inform the development of more realistic semi-analytical and semi-empirical models of dwarf galaxy formation (e.g. Benitez-Llambay & Frenk 2020; Wang et al. 2021; Kravtsov & Manwadkar 2022; O'Leary et al. 2023; Kim et al. 2024; Monzon, van den Bosch & Mitra 2024). The statistical power of these models will be invaluable to interpret the populations of dwarf galaxies from next-generation surveys.

The EDGE collaboration has been a key contributor to these advances. At the core of our approach is the undertaking of suites of highly detailed, high-resolution cosmological simulations (e.g. Agertz et al. 2020; Rey et al. 2020; Prgomet et al. 2022; Andersson et al. 2025). These simulations are then leveraged to pinpoint the key factors shaping dwarf galaxy observables (Rey et al. 2019b, 2022, 2024b; Goater et al. 2024; Gray et al. 2025), and dark matter properties (Orkney et al. 2021, 2022, 2023). Learnings from these simulations are then encapsulated into a semi-analytical model to enable the modelling of statistical populations (Kim et al. 2024).

In this paper, we present the new generation of EDGE cosmological zoomed simulations that will be the cornerstone of our future interpretation efforts. This new suite extends the state of the art in multiple ways. First, in physical fidelity, we maintain the characteristic high resolution of the original EDGE simulations ($\Delta x \approx 3$ pc; $m_{DM} \approx 950\,M_{\odot}$) to resolve individual SNe explosions, and combine it with a systematic use of radiative transfer for an explicit account of radiative feedback from stars. This improves the physical realism of the ISM over many more objects than the single dwarf galaxy presented in Agertz et al. (2020).

Second, we extend the scale of the simulation suite. We simulate $> 15$ different dwarf galaxies from ultrafaint dwarf galaxies ($M_{200} \approx 10^9\,M_{\odot}$; $M_* \approx 10^4\,M_{\odot}$ at $z=0$) to dwarf irregulars ($M_{200} \approx 10^{10}\,M_{\odot}$; $M_* \approx 10^7\,M_{\odot}$). All simulations share the same numerical resolution and feedback physics. Bridging to higher mass objects without compromising on the physical fidelity opens new regimes of comparisons with observational data, as well as providing a larger range of data points to calibrate future semi-analytical models.

This combination of scale and fidelity is unique. Furthermore, the ability to compare between two fundamentally different feedback models (without and with photoionization feedback in EDGE1 and EDGE2, respectively) also allows for a new understanding of the theoretical uncertainties in our results. We describe the EDGE2 model and its updates in Section 2 and show how it affects the ISM and outflows of dwarf galaxies in Section 3. Section 4 shows how dwarf galaxy scaling relations naturally emerge from EDGE modelling and their robustness to a large change in subgrid physics modelling. Section 5 discusses the new numerical challenges associated with explicitly modelling stellar radiation, and we conclude in Section 6.

## 2 EDGE2 NUMERICAL SETUP

In this section, we describe the EDGE2 numerical methods and the suite of zoomed cosmological simulations. We briefly summarize the main differences between EDGE1 and EDGE2 in the bullet points below, and refer the reader to Rey et al. (2020) for a more in-depth description of the EDGE1 model.
(i) EDGE2 uses non-equilibrium primordial and molecular chemistry via the implementation of Rosdahl et al. (2013) and Nickerson, Teyssier & Rosdahl (2018), as opposed to EDGE1's equilibrium cooling from Courty & Alimi (2004).

(ii) EDGE2 accounts for stellar photoionization and photoheating using radiative transfer following the setup described in Agertz et al. (2020).

(iii) EDGE2 updates the feedback budget to be normalized to a Kroupa (2001) IMF integrated between 0.1 and 100 M☉ to be consistent with spectral energy distribution (SED) libraries. EDGE1 used a Chabrier (2003) IMF integrated between 0.5 and 100 M☉.

(iv) EDGE2 tracks the enrichment of eight individual elements (C, N, O, Mg, Al, Si, Fe, and Eu) using the NUGRID yields for core-collapse supernovae (CCSNe; Pignatari et al. 2016; Ritter et al. 2018). EDGE1 tracked only O and Fe and used the Woosley & Heger (2007) CCSNe yields for EDGE1.

(v) EDGE2 updates the Type-Ia SN (SNeIa) model to use a delay-time distribution from Maoz, Mannucci & Brandt (2012) and yields from Seitenzahl et al. (2013). EDGE1 used the binary mass function from Raiteri, Villata & Navarro (1996) and yields originally from Thielemann, Nomoto & Yokoi (1986).

(vi) EDGE2 updates the UV background (UVB) from a modified Haardt & Madau (1996) to the Faucher-Giguère (2020) photoionization and photoheating rates.

(vii) EDGE2 uses a new high-cadence infrastructure to allow the dynamical tracking of gas and stars on short time-scales with tracer particles (Cadiou, Dubois & Pichon 2019).

(viii) EDGE2 adds five new objects at higher masses ($M_{200} \approx 10^{10} M_{\odot}$) to the original suite of initial conditions to grow our sample size (Muni et al. 2025). Table A1 in Appendix A provides all the properties of each object in the EDGE1 and EDGE2 models.

Fig. 1 summarizes how these changes affect the energy injection and metal production budget integrated over 14 billion years in a 300 M☉, 0.01 Z☉ stellar population. We describe these changes further in the next sections.

## 2.1 Initial conditions and resolution

All of our simulated galaxies are evolved to $z=0$ from cosmological zoomed initial conditions constructed with the GENETIC software (Stopyra et al. 2021a). All initial conditions use a flat Lambda-cold dark matter (ΛCDM) transfer function generated with CAMB (Lewis, Challinor & Lasenby 2000). Original EDGE1 initial conditions assume a Planck Collaboration XVI (2014) cosmology with parameters $\Omega_m = 0.3086$, $\Omega_\Lambda = 0.6914$, $h = 0.6777$, $\sigma_8 = 0.8288$, and $n_s = 0.9611$. The new initial conditions we introduce below use a Planck Collaboration VI (2020) cosmology with parameters $\Omega_m = 0.3158$, $\Omega_\Lambda = 0.6842$, $h = 0.6732$, $\sigma_8 = 0.8117$, and $n_s = 0.9660$. These cosmologies are compatible within 68 per cent confidence intervals and do not lead to major differences. Initial conditions are evolved analytically using first-order perturbation theory (Zel'dovich 1970) until $z=99$, at which point numerical integration starts.

To build the original EDGE1 zoomed regions, we first simulate a random cosmological volume (50 Mpc at 512³ resolution) only accounting for gravity and dark matter. We then open a cubic zoomed region of 11.5 Mpc (resolution equivalent to 2048³, $m_{DM} = 3.8 \times 10^6 M_{\odot}$) centred on the largest void in the simulation. Having resimulated this first zoom level, we identify haloes within the $z=0$ void using the HOP halo finder (as in Eisenstein & Hut 1998) and keep only isolated central haloes with no neighbours more massive than them within $5r_{200}$. Here, $r_{200}$ is the radius encompassing

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Energy (top), total oxygen (middle), and iron (bottom) injected by a 300 M☉ stellar particle over its lifetime in the EDGE1 and EDGE2 models. The total is the sum of each individual subcomponent and does not include the radiation budget to enable one-to-one comparison between EDGE2 and EDGE1. Due to the change in IMF, yields and SNeIa modelling, EDGE2 injects $\approx 50$ per cent less CCSN energy, half the oxygen, and twice the iron per stellar population.</figcaption></figure>

200 times the critical density of the Universe. We select six dark matter haloes spanning a wide window in present-day halo mass ($M_{200} = 1.5 \times 10^9$ to $7 \times 10^9 M_{\odot}$, where $M_{200}$ is the mass enclosed in a sphere of radius $r_{200}$), track them back to the initial conditions again, and generate initial conditions at our final zoomed resolution ($m_{DM} = 940 M_{\odot}$, equivalent to $16384^3$).

New initial conditions presented in this paper follow the same procedure as previously followed with small operational differences (see also Muni et al. 2025). The original volume is 100 Mpc, first simulated using $256^3$ particles. We centre on the biggest void identified by running a 'paired' simulation (see Stopyra, Peiris & Pontzen 2021b) and select a $\approx 11.5$ Mpc cubic subvolume in which the resolution reaches $m_{DM} = 3.9 \times 10^6 M_{\odot}$. We select three haloes with present-day halo masses from $M_{200} = 7 \times 10^9$ to $1 \times 10^{10} M_{\odot}$ and increase the dark matter resolution in their Lagrangian region to $m_{DM} = 953 M_{\odot}$. Resolution between the original and updated initial conditions match to $\approx 1$ per cent, with slight shift arising from the slight differences in cosmological parameters.

We also use the 'genetic modification' approach to introduce controlled changes in cosmological initial conditions (see Roth, Pontzen & Peiris 2016; Rey & Pontzen 2018; Stopyra et al. 2021a; Cadiou, Pontzen & Peiris 2021 for more information about the method and the range of possible modifications). Using this approach, all initial conditions are genetically modified to ensure
that the Lagrangian region of our dwarf galaxies is (almost) at rest with respect to the cosmological box (Pontzen et al. 2021). This modification minimizes the streaming of the dwarf galaxy through the simulated volume, mitigating advection errors during integration without affecting its mass growth and local environment.

Furthermore, four of our haloes have genetically modified initial conditions that craft a targeted change to their mass growth history. One low-mass halo ('Halo 1459') is engineered to form systematically earlier and later, at fixed halo mass today (see Rey et al. 2019b). Two intermediate-mass haloes ('Halo 624' and 'Halo 600') are modified to respectively be more massive in halo mass overall and to form later at fixed $z = 0$ halo mass (see Rey et al. 2020). Lastly, a high-mass halo ('Halo 383') is modified to form earlier and later at fixed halo mass today (see Gray et al. 2025). Thanks to the nature of the genetic modification algorithm, all other untargeted aspects of each formation scenario is maximally reproduced (e.g. the large-scale environment and cosmic web topology; see e.g. Pontzen et al. 2017; Rey, Pontzen & Saintonge 2019a for visuals). All haloes labelled 'GM' in Table A1 have had their mass accretion histories genetically modified. Future work will perform controlled modifications to additional haloes in the EDGE2 suite.

## 2.2 Hydrodynamics, radiative transfer, and refinement

We follow the evolution of dark matter, stars, gas, and radiation using the adaptive mesh refinement hydrodynamics code RAMSES-RT (Teyssier 2002; Rosdahl et al. 2013). The dynamics of collisionless particles (dark matter and stars) are computed using a multiscale particle-mesh solver estimating densities through a cloud-in-cell approximation (Guillet & Teyssier 2011). Fluid dynamics are computed using an Harten–Lax–van Leer contact Riemann solver (Toro, Spruce & Speares 1994) with the fluid equations closed by assuming an ideal gas equation of state with adiabatic index $\gamma = 5/3$.

A key addition to this suite is the explicit treatment of the local sources of radiation. We solve the dynamics of the radiation field using the M1 method (Rosdahl et al. 2013; Rosdahl & Teyssier 2015) discretizing the light spectrum in six energy bins from the infrared to the UV (same as in Agertz et al. 2020, table 1). Our energy bins are chosen to track radiation that (i) exerts radiation pressure through dust multiscattering (from 0.1 to 1 eV), (ii) exerts direct radiation pressure (from 1 to 12 eV), (iii) dissociates H₂ (from 12 to 13.6 eV), (iv) ionizes H₁ (13.6 to 24.59 eV), (v) ionizes He I (24.59 to 54.42 eV), and (vi) ionizes He II (> 54.42 eV). We update the average energies and cross-sections in each band every 10 coarse time-steps by computing the luminosity-weighted average over the spectra of all stellar populations in the simulation volume (see Rosdahl et al. 2013). Furthermore, to mitigate computational costs, local radiation around stellar sources is propagated at a reduced speed of light ($c_{reduced} = c/100$; see discussions in e.g. Gnedin & Abel 2001; Rosdahl et al. 2013). A time-varying, but spatially uniform, UVB also permeates the simulation box after reionization (see Section 2.3). This source is not propagated but contributes to photoionization and photoheating rates.

We use the adaptive nature of RAMSES-RT to focus computational power in the dense centre of our dwarf galaxies during integration. We split cells when they contain eight dark matter particles, and when their baryonic mass (stars and gas) exceeds $8m_{bar}$, where $m_{bar} \approx 150 M_{\odot}$. Refinement is allowed down to a maximum resolution of $\Delta x \approx 3$ physical pc. This maximum resolution is achieved throughout the dwarf galaxy's ISM and rapidly degrades as densities decrease (see e.g. Pontzen et al. 2021, fig. 3). We maintain an approximately constant resolution in physical units by releasing

additional new levels every two-folding of the scale factor (Snaith et al. 2018). Refinement is only allowed inside the zoomed region.

## 2.3 Cooling and thermochemistry

Compared to the original EDGE simulations that assumed equilibrium chemistry and cooling, the new simulations of this paper now follow the non-equilibrium thermochemistry of the primordial plasma fully coupled to radiative transfer.

We track the individual ionization fractions of H I, H II, He I, He II, He III, and H₂. The mass fractions and cooling contributions from atomic species are computed using the semi-implicit solver described in Rosdahl et al. (2013), accounting for photoionization, collisional ionization and excitation, bremsstrahlung, Compton cooling and heating from the cosmic microwave background, and di-electronic recombination. We assume that ionizing radiation from free-bound ground state recombination radiation is absorbed locally (the on-the-spot approximation; see e.g. Nebrin 2023). We follow the formation, advection and destruction of H₂ and its contribution to the cooling rate following the model described in Nickerson et al. (2018). This model accounts for gas-phase, dust-phase, and collisional H₂ formation, destruction through photodissociation, photoionization, and cooling and heating coupled to the local radiation flux. In all runs, we turn-off cosmic ray heating and ionization due to the large uncertainties on the cosmic ray ionization rate in dwarf galaxies. We leave to future work an exploration of its impact on the ISM and outflows of dwarf galaxies (see also Martin-Alvarez et al. 2023).

Metal cooling uses tabulated cooling rates scaled by metallicity. Below $10^4$ K, we use the Rosen & Bregman (1995) fine structure cooling rates, while above $10^4$ K, we use CLOUDY (Ferland et al. 2017) tables assuming our updated UVB (Faucher-Giguère 2020). These rates are scaled by the total metallicity defined as $Z = (2.09 Y_O + 1.06 Y_{Fe})/Z_{\odot}$, where $Z_{\odot} = 0.02$. The numerical coefficients are derived from assuming a solar mixture from Asplund et al. (2009) and $Y_O$ and $Y_{Fe}$ are the mass fractions of oxygen and iron in the gas cell (see also Kim et al. 2014).

Ionization and heating from the UVB permeates the whole volume in an optically- thin approximation, with high-density regions allowed to self-shield from this background. We exponentially damp the UVB photoionization and photoheating rates according to $\exp(-n_H/10^{-2} \text{cm}^{-3})$ to approximate self-shielding above $n_H \approx 10^{-2} \text{cm}^{-3}$ (Aubert & Teyssier 2010; Rosdahl & Blaizot 2012). Local self-shielding from stellar radiation is treated self-consistently by solving of the coupled equations of radiative transfer and thermochemistry. Optically thin UVBs can overheat the early intergalactic medium (Oñorbe, Hennawi & Lukić 2017), so we also exponentially damp the UVB as a function of redshift by $\exp(z_{reion} - z)$. This ensures a continuous ramp up of photoionizing and photoheating rates, rather than the instantaneous switch-on commonly used in UVB tables (e.g. Faucher-Giguère 2020). With our choice of $z_{reion}$, the UVB heating and ionizing rates are at full strength at $z = 6$ (see Rey et al. 2020).

## 2.4 Star formation

The star formation modelling remains unchanged compared to Agertz et al. (2020). We model star formation following a Kennicutt–Schmidt law (Schmidt 1959; Kennicutt 1998):

$$ \dot{\rho}_{*} = \epsilon_{\text{ff}} \frac{\rho_{g}}{t_{\text{ff}}} \quad \text{for gas cells with } \rho_{g} > \rho_{*} \text{ and } T_{g} < T_{*}, \qquad (1) $$
where $\dot{\rho}_*$ is the instantaneous star formation rate (SFR) in a gas cell, $\epsilon_{\text{ff}}$ is the star formation efficiency per free-fall time, $\rho_g$ and $T_g$ are the gas cell density and temperature, and $t_{\text{ff}} = \sqrt{3\pi/32G\bar{\rho}}$ is the local free-fall time. We choose $\rho_* = 300 m_p \text{ cm}^{-3}$, corresponding roughly to the average density of giant molecular clouds and to the density needed to form a $300 M_{\odot}$ stellar population at our resolution (see below). We pick a constant $\epsilon_{\text{ff}} = 10$ per cent which, at $\approx 3$ pc resolution, has been shown to reproduce key star formation observables such as the density and energy power spectra of the ISM in Milky Way-like galaxies (Grisdale et al. 2017), the properties of giant molecular clouds (Grisdale et al. 2018), and the observed star formation efficiency in molecular clouds (Grisdale et al. 2019). This value is also compatible with typical predictions from multi-free-fall models and star formation simulations with similar resolution (e.g. Federrath & Klessen 2012; Padoan, Haugbølle & Nordlund 2012). We also implement a temperature threshold of $T_* = 1000$ K, which is higher than in EDGE1 (100 K) to compensate for the additional heating from radiative feedback. The original motivation for $T_*$ is to ensure that star-forming gas is cold enough to represent (potentially unresolved) molecular gas. As we will see in Section 5.1.4, however, the choice of this parameter leads to unwanted numerical effects due to the coupling between star formation and radiative heating.

For every gas cell satisfying our threshold conditions, we sample equation (1) stochastically with a Poisson process so that the expectation of the local SFR equals $\dot{\rho}_*$ (Rasera & Teyssier 2006). Stellar particles have initial masses of $300 M_{\odot}$ and assume a Kroupa (2001) IMF. This mass is at the boundary of ensuring that stochastic effects are averaged over at our resolution (see e.g. Smith 2021). This is key for radiative feedback which is still modelled as an IMF-averaged process (as opposed to SN feedback for which SNe are individually sampled). We will explore the consequences of these choices in the future with a star-by-star EDGE implementation (Andersson et al. 2025).

## 2.5 Stellar feedback

A key aspect of the EDGE model is the detailed account of feedback from massive stars. We account for CCSNe, SNeIa, fast winds from massive stars, slow winds from asymptotic giant branch (AGB) stars, and now radiative feedback.

Our model tracks the main-sequence lifetime of different progenitors within a stellar particle (Agertz, Teyssier & Moore 2011), ensuring that stars of different masses inject their feedback on their relevant main-sequence time-scale. We model SN explosions as discrete events, computing at each simulation time-step the number of stars exiting the main sequence to turn into CCSNe (equations 6 in Agertz et al. 2013). This IMF-averaged number is then randomly sampled through a Poisson process to obtain a discrete number of explosions (Agertz et al. 2020). To maintain consistency with the IMF used for radiative feedback, we normalize this number to the Kroupa (2001) IMF integrated between 0.1 and $100 M_{\odot}$. This differs from EDGE1 which used a Chabrier (2003) IMF integrated between 0.5 and $100 M_{\odot}$, leading to a $\approx 40$ per cent decrease in the number of CCSNe per stellar population (Fig. 1, top).

Furthermore, to match the explodability assumptions made by our new chemical enrichment model (Section 2.6), we now assume that massive stars between 8 and $30 M_{\odot}$ explode as CCSNe, while stars above directly collapse into black holes without releasing energy. Previous EDGE1 modelling assumed that all stars between 8 and $40 M_{\odot}$ exploded. When convolved with the IMF, this further reduces the energy budget of CCSNe by $\approx 10$ per cent.

Furthermore, SNeIa now follow a delay time distribution according to Agertz et al. (2021). We assume that SNeIa start when their parent stellar particle has an age of 38.4 Myr (the main-sequence lifetime of $8 M_{\odot}$ star to form the first non-exploding degenerate progenitor) and explode with a rate that initially peaks at $2.6 \times 10^{-13} M_{\odot}^{-1}$ and decays as $t^{-1.12}$ following the empirical determination of Maoz et al. (2012) for field galaxies. As for CCSNe, SNeIa are individually sampled from this rate using a Poisson process. In EDGE1, we instead modelled the SNeIa rate by integrating over the IMF of secondary, binary companions between 1 and $8 M_{\odot}$ (Raiteri et al. 1996; see Agertz et al. 2013 for a detailed description). This modernization leads to a higher number of SNeIa per stellar population (Fig. 1, top), which balances the 40 per cent decrease induced by the change in IMF normalization.

For both SN types, we directly inject thermal energy ($E_{\text{SN}} = 10^{51}$ erg constant at all times) when the cooling radius of the individual SN event is resolved by six resolution elements. This allows us to self-consistently follow the build-up of momentum through the Sedov–Taylor phase by solving the hydrodynamics equations rather than relying on a subgrid implementation. When the cooling radius is unresolved, we switch to a momentum injection. This is a binary switch (see also Kim & Ostriker 2015 for a similar implementation and Kimm & Cen 2014; Kimm et al. 2015 for an interpolated approach). Following Kim & Ostriker (2015), we inject a Sedov–Taylor momentum of $\rho_{\text{ST}} = 2.95 \times 10^5 M_{\odot} \text{ km s}^{-1}$ scaling with density and gas metallicity according to Blondin et al. (1998). We do not account for local velocity dynamics (e.g. Hopkins 2025). Appendix B shows that $> 90$ per cent of CCSNe events and $> 60$ per cent of SNeIa events are resolved by more than $6^3$ ($> 200$) resolution elements, ensuring that we greatly reduce uncertainties associated to the subgrid modelling of SN feedback in an unresolved ISM.

We also account for wind feedback from massive O and B stars ($m \ge 8 M_{\odot}$) during their main sequence, and from stars of masses $0.5-8.0 M_{\odot}$ when they reach their AGB phase. OB winds return mass and momentum to the ISM according to the metallicity-dependent budget described in Agertz et al. (2013). AGB winds continuously release mass and metals over the lifetime of a stellar particle following their IMF-averaged mass loss (equation 17 in Agertz et al. 2013). No changes are made to these models between EDGE1 and EDGE2 and at our characteristic low metallicities, these wind models contribute little to the feedback energy budget (Fig. 1).

In addition to these changes to the SN budget between EDGE1 and EDGE2, a large change in stellar feedback comes from the addition of radiative feedback in EDGE2. In this case, each stellar particle injects radiation according to its age and metallicity following a Bruzual & Charlot (2003) SED. This SED is a ‘soft’ choice for low-metallicity stellar populations, and we explore in Section 5.1.3 the impact of a harder SED accounting for binary stars (Stanway, Eldridge & Becker 2016).

## 2.6 Metal enrichment

Another novelty in this suite is the tracking of new chemical elements. In addition to O and Fe tracked by the previous model (Agertz et al. 2013), we now account for the production and advection of C, N, Mg, Al, Si, and Eu. These elements were chosen to sample main α elements commonly observed in stars within dwarf galaxies, as well as the r-process element Eu.

As with previous EDGE simulations, we do not model the formation, feedback, and metal enrichment from primordial metal-free
stars. Instead, we initialize all simulations with a floor in oxygen metallicity of $Z_o = 10^{-3} Z_⊙$ as an approximation for primordial metal enrichment (e.g. Greif et al. 2010; Jaacks et al. 2018; Visbal, Bryan & Haiman 2020; Brauer et al. 2025). All other elements have vanishing mass fractions initially. Lowering this floor by an order of magnitude leaves the stellar mass and metallicities of our dwarfs unchanged (Agertz et al. 2020).

Once star formation starts, winds from O and B stars, winds from AGB stars, CCSNe, and SNeIa all inject chemical elements on the same time-scale as their feedback (Section 2.5). Once injected, each element is advected passively with the gas.

Yields for CCSNe, OB winds, and AGB winds are interpolated from the tables provided by NUGRID (Pignatari et al. 2016; Ritter et al. 2018). Discrete injection events from CCSNe are interpolated across the table in progenitor mass and metallicity. Since the NUGRID tables are only available for CCSN progenitor masses between 12 and 25 $M_⊙$, we linearly extrapolate to cover the range assumed in our model (8–30 $M_⊙$). Future EDGE simulations will leverage new yield tables explicitly including low progenitor masses (Limongi et al., in preparation) to alleviate this issue, as well as updating our assumptions for the CCSN progenitor mass range to reflect new understanding of the super-AGB phase and electron-capture SNe (Gil-Pons et al. 2018; Limongi et al. 2024).

Discrete SNeIa inject chemical elements according to Seitenzahl et al. (2013) assuming a (constant) solar metallicity (metallicity dependence of these yields is weak). For AGB winds, we compute the IMF-averaged number of AGB stars in a given time-step (Agertz et al. 2013) and inject the corresponding yields from NUGRID for each element. OB winds follow the same procedure, but instead inject according to a time-dependent fitting function that was calibrated in Agertz et al. (2013).

Combined, these updates have a significant impact on the metal production budget per stellar population. Fig. 1 shows that oxygen production is roughly halved in EDGE2 (middle panel), primarily driven by a strong decrease in CCSNe yield. Part of this decrease is due to the change in IMF normalization reducing overall CCSNe numbers by ≈ 40 per cent, but is strongly driven by the change in yield tables. EDGE1 tables (Woosley & Heger 2007) had an exponential scaling of oxygen production with progenitor mass (see e.g. Kim & Ostriker 2015, equation 5) with high-mass stars producing copious amounts of metals. This was known to significantly overproduce oxygen and alpha elements compared to Milky Way observations (e.g. Agertz et al. 2021). The scaling with progenitor mass is much weaker in the NUGRID tables. Importantly, iron production is increased by a factor ≈ 2.5 in EDGE2 (bottom panel), again largely driven by the change in CCSN yield table. This reflects the recent trends in yield computations with newer models all producing significantly more iron than previous iterations (e.g. Limongi & Chieffi 2018).

Finally, for r-process production of Eu, we include an effective model for neutron-star mergers (NSNS) inspired by Naiman et al. (2018). We assume a constant relative fraction $4.6 \times 10^{-2}$ between NSNS and SNeIa consistent with the NSNS rates observed in the local Universe (LIGO and Virgo Scientific Collaboration 2017). NSNS events are sampled discretely using the same mechanism as SNeIa with each NSNS injecting a europium yield of $10^{-5} M_⊙$ (Côté et al. 2018) but no other elements. NSNS events do not inject energy, mass, or momentum.

## 2.7 High-cadence outputs

The Eulerian nature of our code prevents us from efficiently tracking the Lagrangian history of gas flows. This is particularly problematic

in dwarf galaxies where dynamical times and cooling times in the centre are much shorter than the cadence at which we can save simulation outputs (see e.g. Rey et al. 2022, 2024b for examples of these limitations).

To remedy to this, EDGE2 uses a Monte-Carlo particle tracer algorithm (Cadiou et al. 2019). Tracers are designed to statistically track gas flows and exchange mass with stellar tracers to track the full baryon cycle of gas as it is accreted and recycled through star formation. Tracer dynamics is solved using the same physical solvers as the rest of the simulation, but they do not source or contribute towards hydrodynamical, gravity, or radiative forces. We spawn five tracers per high-resolution gas cell, leading to a tracer mass of 223 $M_⊙$ comparable to the stellar particle mass. Our smallest dwarf galaxies have ≈ 2 million tracers, while our most massive objects can have up to 20 million tracers.

The position and velocity of tracers is stored on disc every 4.5 Myr (compared to every 100 Myr for full simulation output), along with the density, pressure, and gravitational potential of the gas at the location of the tracer. Additionally, the same information with the same cadence is dumped for a selection of 10 000 dark matter particles, selected randomly from the main progenitor of the dwarf at $z=2$.

## 2.8 Data processing and analysis

We process each EDGE2 simulation with the ADAPTAHOP halo and subhalo finder (Aubert, Pichon & Colombi 2004; Tweed et al. 2009) retaining dark matter structures with more than 100 particles. EDGE1 simulations were processed with the HOP halo finder (Eisenstein & Hut 1998). We match haloes and subhaloes between simulation snapshots to build merger trees using the PYNBODY (Pontzen et al. 2013) and TANGOS (Pontzen & Tremmel 2018) libraries. Halo centres are identified using the shrinking-sphere algorithm (Power et al. 2003).

# 3 THE IMPACT OF RADIATIVE FEEDBACK ON DWARF GALAXIES

## 3.1 The stellar mass–halo mass relation

Fig. 2 shows how the integrated stellar masses, $M_*$, respond to the change from the original EDGE1 (blue circles) to the updated EDGE2 (red diamonds) model. $M_*$ is computed by summing all stellar particles within $r_{200}$ and lines connect dwarf galaxies sharing the same cosmological initial conditions and formation scenarios (see the legend for individual names).

Comparing our results with empirical determinations of the $M_* - M_{200}$ relation from observed dwarf galaxies (left-hand and middle panels) shows that EDGE dwarf galaxies are within observational uncertainties. In both the EDGE1 and EDGE2 models, low-mass dwarf galaxies ($M_* \approx 10^4 - 10^6 M_⊙$; $M_{200} \approx 10^9 M_⊙$) scatter around the median inferred from Milky Way satellites (blue line showing the peak halo mass from Nadler et al. 2020) and span the breadth of the 16–84 confidence interval at fixed $M_{200}$ (blue contours). In our suite, this scattering around the median is largely driven by our systematic exploration of different formation histories at fixed halo mass (Rey et al. 2019b; Gray et al. 2025).

An agreement between Milky-Way satellites and simulated field dwarfs might appear surprising at first. But in the low-mass regime, the shutdown of gas accretion in the smallest haloes due to cosmic reionization ('reionization feedback') is the dominant mechanism regulating the stellar mass and gas content (e.g. Efstathiou 1992; Gnedin 2000; Hoeft et al. 2006; Okamoto, Gao &
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Stellar and halo masses of dwarf galaxies simulated with the EDGE1 and EDGE2 models (circles and diamonds, respectively). Both suites are broadly compatible with the stellar-mass-halo-mass relation inferred from Milky Way satellites at the low-mass end (left-hand panel, contours; Nadler et al. 2020). Higher mass objects ($M_{\star} \approx 10^{7} M_{\odot}$) in EDGE2 closely match the empirical data measured in isolated dwarf irregulars (left-hand panel, grey points; Read et al. 2017). This match to observational data further extends to stellar sizes, gas contents, and stellar and gas-phase metallicities (Figs 6–9). The average difference in $M_{\star}$ between EDGE1 and EDGE2 is $\approx 20$ per cent, i.e. well converged compared to the scatter in $M_{\star}$ at given $M_{200}$ between different galaxy formation models (right-hand panel, grey points). This is primarily driven by a change in ISM structure (Fig. 3) and galactic outflow strength (Fig. 4). Section 3.1 discusses the response of individual formation histories (symbols linked with a line and caption).</figcaption></figure>

Theuns 2008; Noh & McQuinn 2014; Benitez-Llambay & Frenk 2020). Furthermore, few faint and ultrafaint observed dwarf galaxies have orbital parameters that lead to strong tidal interactions with the Milky Way (Simon 2018; McConnachie & Venn 2020; McConnachie et al. 2021). As a result, the faint-end of the $M_{\star} - M_{200}$ relation is likely to be shaped by their pre-reionization ($z \ge 6$) evolution, rather than limited by the specific environment of our Galaxy.

As dwarf galaxies grow in mass, ($M_{\star} \ge 10^{6} M_{\odot}$; $M_{200} \ge 5 \times 10^{9} M_{\odot}$), reionization feedback becomes less dominant. This enables field dwarf galaxies to accrete late-time gas and enable dark matter halo mass measurements through rotation curves (Read et al. 2017; grey points in left-hand panel). At high $M_{200}$ all the way to $M_{200} \approx 10^{10} M_{\odot}$, our simulated dwarf galaxies closely match the individual measurements of field, gas-rich dwarf irregulars. But, importantly, the two lowest mass objects in the Read et al. (2017) sample lack a simulated counterpart. Rather than a failure of the model, this more likely reflect the rarity of finding such low-mass, gas-rich galaxies on which to perform rotation curve measurements. At halo masses $M_{200} \approx 2 \times 10^{9} M_{\odot}$, only very specific formation histories that form especially late will allow gas accretion and star formation (Benitez-Llambay & Fumagalli 2021). Such objects are not present either in the EDGE1 and EDGE2 suites, but calibrating

a semi-analytical on EDGE1 results to create a large population naturally recovers them (Kim et al. 2024).

On average, $M_{\star}$ in EDGE1 and EDGE2 differ by 34 per cent across the suite, and are within a factor of two of each other for all galaxies except for three specific formation histories that we discuss further below. The magnitude of these shifts is encouragingly tight given the large changes in cooling and heating physics between models. To emphasize that such shifts are well within theoretical uncertainties, the right-hand panel of Fig. 2 shows a compilation of simulated field dwarf galaxies (Benitez-Llambay & Fumagalli 2021; Herzog et al. 2023, plusses, ‘Swift’; Revaz & Jablonka 2018, crosses, ‘Gear’; Wang et al. 2015; Tollet et al. 2016, dots, ‘Nihao’; Fitts et al. 2017, Wheeler et al. 2019, stars, ‘FIRE-2’; Munshi et al. 2019, 2021, pentagons, ‘Changa’]. At given $M_{200}$, predictions from different simulation groups can differ by over an order of magnitude in $M_{\star}$ (e.g. pentagons against stars around $M_{200} \approx 10^{10} M_{\odot}$), highlighting the small-in-comparison shifts between EDGE1 and EDGE2.

Before turning to comparing our dwarf galaxies to observations (Section 4), we highlight several trends in Fig. 2 that will help us establish differences between our two numerical models:

(i) EDGE2 systematically suppresses stellar masses in low-mass systems ($M_{\star} \le 10^{6} M_{\odot}$). This trend is reversed in higher mass
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Temperature–density diagrams of the ISM of the dwarf galaxy ['Halo 383 (early)'; $M_{\star} \approx 10^{7} M_{\odot}$] with the EDGE1 (left) and EDGE2 model (right). The diagrams are 2D mass-weighted probability density function (PDF) averaged over the last 4 billion years of evolution. Many features are shared between models (annotated). The addition of radiative feedback leads to warm ($T \approx 10^{4} \text{ K}$) and dense ($n_{H} \ge 10 \text{ cm}^{-3}$) ISM gas. With this change, the ISM is also significantly denser overall (top panels) and the strength of galactic outflows is reduced (Fig. 4).</figcaption></figure>

dwarfs ($M_{\star} \ge 10^{6} M_{\odot}$), for which $M_{\star}$ is systematically increased. Sections 3.2 and 3.3 show that the inclusion of radiative feedback leads to a fundamentally different structure of the ISM and reduces the efficiency of galactic outflows. This suppression in outflow loading factors is more pronounced at lower masses, driving $M_{\star}$ down, and less at higher masses, driving $M_{\star}$ up.

(ii) One of our simulated dwarfs ('Halo 383: GM higher mass', purple in the right-hand panel) significantly overshoots the stellar-mass–halo-mass relation, with data for this simulation shown at $z = 0.9$ when we stopped the simulation due to the numerical cost incurred by the high $M_{\star}$. Little $M_{200}$ growth is expected after this time. Section 5 shows that, at high $M_{\star}$, the reduced strength of galactic outflows stems from numerical issues leading to increasingly difficult regulation. We stopped 'Halo 153' (pink) at $z = 0.6$ for the same reason.

(iii) Two low-mass galaxies show a large $M_{\star}$ suppression with EDGE2 (right-hand panel, 'Halo 600: GM later' in orange and 'Halo 1459: GM latest' in turquoise). These two objects share the same characteristic of assembling late for their halo masses, from multiple small building blocks during the reionization era each with $M_{200}(z = 6) \le 10^{8} M_{\odot}$ that later come together in dry mergers (Rey et al. 2019b, 2020). Each of these small building blocks sees their $M_{\star}$ suppressed by radiative feedback, compounding the effect on the total $M_{\star}$ compared to a comparatively massive object at the time of reionization.

## 3.2 The structure of the ISM

Fig. 3 shows the difference in ISM structure between models using SN-only feedback (left) and including radiative feedback

(right). For illustrative purposes, we plot the 2D temperature–density distributions averaged over the ≈ 100 snapshots along the formation history of 'Halo 383 (early)' ($M_{\star} \approx 10^{7} M_{\odot}$). We pick this galaxy as it has a stable and gas-rich ISM at all times, in both models. But all points highlighted below generalize across the suite.

In both models, we observe similar features in the temperature–density diagram (labelled on plot). Namely, (i) the upper tail of the IGM temperature–density relation (diagonal track at low densities; see McQuinn 2016 for a review); (ii) the distinctive thermal equilibrium curve around $10^{4} \text{ K}$ where radiative cooling balances photoheating from the external UVB (see e.g. Smith et al. 2017, fig. 8); (iii) the break of this thermal equilibrium when gas starts self-shielding against the UVB; and (iv) high-temperature ($T \ge 10^{5} \text{ K}$) gas at low densities resulting from SN-driven outflows.

The first notable difference is the warm ($T \ge 10^{4} \text{ K}$) and dense ($\rho \ge 10 \text{ m}_{\text{p}} \text{ cm}^{-3}$) gas present in the EDGE2 model (right-hand panel) and absent in EDGE1 (left-hand panel). This is the direct result of the explicit modelling of radiative feedback in H II regions around stars, ensuring that EDGE2 galaxies now capture a key gas phase for the overall structure of the ISM and its emission lines. The 'doubled' horizontal track results from the mixing of photoionized gas with the surrounding ISM – the upper track is fully ionized, while the lower track has significant neutral fraction, driving a difference in molecular weight and thus a ≈ 0.5 shift in temperature. Part of this mixing is physical, at the edges of H II regions, and part is numerical (Section 5.2).

The second difference is that the ISM is overall denser with non-equilibrium cooling and photoionization feedback (top panels show marginalized density PDFs). Surprisingly, this trend extends beyond the density threshold for star formation (dashed lines), with the EDGE2
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Mass-loading factors measured through a spherical shell at 0.25 $r_{200}$ as a function of galaxy stellar mass. Error bars show the 16–84 confidence interval around the median over the cosmological history of each individual dwarf galaxy. Radiative feedback (diamonds) reduces outflow loading factors for every galaxy compared to SN-only feedback (circles). This reduction is up to 1.5 dex for the least massive dwarf galaxies.</figcaption></figure>

model maintaining a significant amount of gas at densities above the star formation threshold (≥ 300 $m_p$ cm$^{-3}$; grey dashed) that contrasts with the plummeting PDF in EDGE1. This reflects our numerical choices in the star formation algorithm, which we explore further in Section 5. As we show next, the much denser ISM in EDGE2 correlates with strongly reduced galactic outflows.

## 3.3 The strength of galactic outflows

Fig. 4 shows the gas mass-loading factor, $\eta_M$, over the cosmological history of each galaxy in the suite as a function of their $M_\star$. We define $\eta_M = M_{\text{out}}/\text{SFR}_{10\,\text{Myr}}$, where the mass outflow rate $\dot{M}_{\text{out}}$ is measured through a spherical shell centred on the galaxy that spans a radial range between 0.2 $r_{200}$ and 0.3 $r_{200}$ and only includes outflowing gas (see also Rey et al. 2024a, for further details). SFR$_{10\,\text{Myr}}$ is the star formation rate averaged over 10 Myr. Error bars show the 16–84 confidence interval of the distribution of $\eta_M$ values over time, with the median shown as a symbol. Some galaxies for which no saved snapshots have SFR$_{10\,\text{Myr}} > 0$ are missing from the plot.

Fig. 4 demonstrates that radiative feedback systematically decreases the strength of galactic outflows, reducing their ability to remove mass from the central galaxy. Scatter as a function of time is large, reflecting the stochastic nature of the star formation-feedback cycle, but the suppression is systematic and close to an order of magnitude. Since $M_\star$ remains close between EDGE1 and EDGE2 (Fig. 2), this highlights a fundamental change in the way star formation is regulated. With radiative feedback, gas is prevented from forming stars by gentle heating, rather than being mechanically removed from the centre in blastwaves when considering only SN feedback. We further show in Appendix C that the star formation histories (SFHs) are on average less bursty as a result, aligning with this picture.

A reduction in the strength of galactic outflows and of the burstiness of star formation due to radiative feedback is well established in isolated non-cosmological dwarf galaxies (e.g. Emerick, Bryan & Mac Low 2018; Smith et al. 2021; Deng et al. 2024b; Andersson et al. 2024) and was previously noted for a single low-mass cosmological dwarf in Agertz et al. (2020). Our simulations extend these findings to a much wider range of masses and highlight trends with host $M_\star$.

In particular, the median $\eta_M$ reduces with increasing $M_\star$ in the EDGE1 model, as expected from the scaling of SN-driven outflows with $M_\star$ (e.g. Muratov et al. 2015; Christensen et al. 2016; Nelson et al. 2019; Mitchell et al. 2020; Pandya et al. 2021). In contrast, EDGE2 sees $\eta_M$ stay relatively flat over a wide range of $M_\star$, even slightly increasing with mass. This is a direct consequence of the different feedback mechanisms at play in the two models.

In EDGE1, there is no rapid feedback mechanism to regulate star formation, the clustering of SNe across all $M_\star$ is high, and the efficiency of mechanical outflows decreases as the gravitational potential wells get deeper.

In contrast, EDGE2 sees some of the weakest outflows and the strongest outflow suppression compared to EDGE1 in the lowest mass objects. This reflects the unique regime of very low-mass galaxies that formed all their stars before reionization.

In this case, stellar radiation efficiently escapes the ISM of the small building blocks that will make the $z=0$ galaxy but remain spatially distinct at $z>10$ (Fig. 5, left-hand panel). The radiation propagates outwards much more rapidly than a mechanical outflow, heating up gas to $T \ge 10^3$ K well beyond the main progenitor's $r_{200}$ (right-hand panel).<fnref n="1" /> This in turn suppresses gas inflows and correlated star formation across a much larger volume. This regulation mechanism is fundamentally different to the mechanical outflows driven by SNe in EDGE1 (middle panel). As the galaxies get larger in mass, this effect occurs at ever higher redshift and contributes less to shaping the overall $M_\star$ of the galaxy.

Lastly, we note that observational values for $\eta_M$ at $M_\star \approx 10^7 M_\odot$ have varied from $\approx 10$ (e.g. Chisholm et al. 2017), to $\approx 0.1-1$ (e.g. McQuinn, van Zee & Skillman 2019; Xu et al. 2022), to $10^{-2}$ (e.g. Marasco et al. 2023), bracketing the values reported here. There are clear limitations in directly comparing the values in Fig. 4 to observations. Here, $\eta_M$ is computed including all gas phases rather than those that are bright in emission lines used to make the measurement. We also use a radius well outside that probed by observations which would likely affect the measured $\eta_M$ (e.g. Muratov et al. 2015). Finally, the values in Fig. 4 are averaged over the full cosmological history, rather than at a time of star formation when the dwarf galaxy's gas is observable. A more careful comparison is thus warranted to establish the realism of either EDGE1 and EDGE2 outflows. This will be the focus of a future paper.

# 4 THE EMERGENCE OF DWARF GALAXY SCALING RELATIONS AND THEIR DIFFERENTIATING POWER

## 4.1 The stellar size–luminosity relation

Fig. 6 shows the absolute V-band magnitude, $M_V$, and projected half-light radius, $r_{1/2,V}$, of simulated dwarf galaxies. To obtain these

<footnote n="1">¹Note that the temperatures in Fig. 5 are such that most of the hydrogen in the volume remains neutral, thereby not impacting the cosmological timing of hydrogen reionization.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Gas density integrated along the line of sight (left) and thin slices of the temperature around the same low-mass dwarf galaxy ($M_* \approx 10^5 M_\odot$) at $z = 10$, in EDGE1 (middle) and EDGE2 (right). Accounting for radiative feedback has a dramatic impact beyond $r_{200}$ (right-hand panel). Radiation efficiently escapes from the dwarf's small building blocks, heating the volume to $T \ge 10^3$ K and suppressing inflows and correlated star formation across a large volume. Hot, SN-driven outflows eventually achieve the same effect (middle), but on a time-scale longer than radiative feedback.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Absolute V-band total magnitude of simulated dwarf galaxies against their V-band half-light radius. With both models, simulated dwarf galaxies populate the observational scatter observed around the Local Volume (grey points; error bars showing the dispersion around the median) showcasing the limited power of stellar sizes in discriminating between galaxy formation models. Three exceptions lack an observational counterpart (labelled), but do not pose significant challenges to either model (see the text discussion).</figcaption></figure>

quantities, we compute the luminosity of all star particles within $r_{200}$ as a function of their mass, age, and metallicity according to the single stellar population model derived from PARSEC isochrones (Bressan et al. 2012; Nguyen et al. 2022). $r_{1/2,V}$ is then obtained along a random line of sight. Simulated properties are compared to observational data from the compilations of McConnachie (2012), Kirby et al. (2013, 2014), and Simon (2019) augmented with individual candidates and detections from Torrealba et al. (2016, 2018, 2019), Homma et al. (2019, 2024), Mau et al. (2020), Bennet et al. (2022), Richstein et al. (2022), Sand et al. (2022), Cerny et al. (2023b, c, a), Jones et al. (2023, 2024), McQuinn et al. (2023), Collins et al. (2024), Li et al. (2024),

Martinez-Delgado et al. (2024), McNanna et al. (2024), Smith et al. (2023, 2024), and Tan et al. (2025). We curate the observational data to use the most recent reference when multiple observations of the same system are available.

Irrespective of the model, nearly all dwarf galaxies are within the observational scatter of $M_V-r_{1/2,V}$ and there is no clear trend or systematic offset between models. On average, EDGE2 dwarfs are marginally larger than their EDGE1 counterparts, and high-mass EDGE2 dwarfs are systematically brighter in line with their increased $M_*$ (Fig. 2).

Two simulated dwarfs lack direct observed analogues (labelled on plot). In the bottom right, the same dwarf galaxy in both models is an ultrafaint ($M_V \ge -6$) that is also extremely extended ($r_{1/2,V} \ge 1$ kpc). This large $r_{1/2,V}$ is driven by the specific assembly history of this galaxy that sees its stellar component built from the mergers of small building blocks that deposit stars on wide orbits (see Rey et al. 2019b for details). As discussed in Section 3.1, photoionization feedback is particularly efficient in these small building blocks, and this object becomes even more diffuse in the EDGE2 model. In both cases, the lack of observational counterpart can be explained by their central surface brightnesses being beyond currently observable capabilities ($\mu_{0,V} > 31$ mag arcsec$^{-2}$). Such dwarfs should be revealed by forthcoming surveys such as on the Vera Rubin C. Observatory.

In the left-hand corner, one EDGE1 dwarf (blue point) has a particular formation history that leads to a compact nuclear star cluster that dominates the light of the galaxy. This galaxy shares observational characteristics with ultracompact dwarfs and nuclear star clusters not included in Fig. 6 (see Gray et al. 2025 for a discussion). Due to the reduced outflow strength and less clustered star formation, this galaxy does not qualify as a nuclear star cluster in EDGE2.

Our results highlight the natural emergence of the $M_V-r_{1/2,V}$ scaling relation, in both median and scatter, from our cosmological modelling. This emergence is robust to a large change in galaxy formation modelling from EDGE1 to EDGE2. Such convergence hints that, if stellar masses are reasonably predicted at a given $M_{200}$, the $M_V-r_{1/2,V}$ relation and its scatter follow from the diversity of cosmological assemblies. This in turn brings confidence that this relation can be well predicted with physics-informed semi-analytical arguments applied to dark-matter-only simulations (Nigudkar et al., in preparation).
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Absolute total V-band magnitude of our simulated galaxies against the average iron content of their stars. Both models track the slope of the mass–metallicity relation, but are offset with one another due to the changes in metal production tables and feedback physics.</figcaption></figure>

## 4.2 The stellar mass–metallicity relation

Fig. 7 shows the $M_V$–$\langle[\text{Fe/H}]\rangle$ relation where $\langle[\text{Fe/H}]\rangle$ is the average stellar iron abundance. We derive [Fe/H] for each stellar particle within $r_{200}$ from their iron mass fractions and compute $\langle[\text{Fe/H}]\rangle$ as the mass-weighted average (see Escala et al. 2018, equations 3 and 4). Observed data are compiled from Kirby et al. (2013, 2014) and Simon (2019) augmented and updated with data from Kirby et al. (2017, 2020), Li et al. (2017, 2018), Longeard et al. (2018), Fritz et al. (2019), Ji et al. (2019, 2021), Collins et al. (2020, 2021), Pace et al. (2020), Taibi et al. (2020), Wojno et al. (2020), Jenkins et al. (2021), Chiti et al. (2021, 2022), Bruce et al. (2023), Charles et al. (2023), Cerny et al. (2023a), Smith et al. (2023), Hansen et al. (2024), Heiger et al. (2024), Kvasova, Kirby & Beaton (2024), and Tan et al. (2025). Data error bars show the dispersion around $\langle[\text{Fe/H}]\rangle$ (when available) rather than measurement errors on the mean.

Both EDGE1 and EDGE2 models accurately track the slope of the $M_V$–$\langle[\text{Fe/H}]\rangle$ relation and are within the observational scatter. However, EDGE1 dwarfs populate the lower end of observed points over the range of $M_V$, while EDGE2 dwarfs populate the upper end. The origin of the ≈ 0.5 dex offset between models is twofold. First, the EDGE2 model produces ≈ 2.5 more iron per stellar mass formed than EDGE1 (Fig. 1). And second, the inclusion of radiative feedback leads to a more gentle regulation of star formation with weaker galactic outflows (Section 3), and thus metal retention in the ISM. These two effects compound one another, with both iron production and retention increased in the EDGE2 ISM, in turn leading to higher $\langle[\text{Fe/H}]\rangle$.

These results cement the $M_V$–$\langle[\text{Fe/H}]\rangle$ relation as a sensitive probe of star formation and stellar evolution physics at low metallicities (see also Agertz et al. 2020; Prgomet et al. 2022; Sanati et al. 2023). The normalization at a given $M_V$ is in particular directly related to the strength of galactic outflows in dwarf galaxies, while the slope emerges from the cosmological relation between stellar and

<figure><img src="image_11.png" /><figcaption>**Figure 8.** H1 and stellar masses of our simulated dwarf galaxies in each model. The structure of the $M_\star – M_{H1}$ is robust to our change of model, with a bimodal structure at the low-mass end and significant time variability in H1 content due to stellar feedback (violins show the distribution over the last 4 Gyr). Both models are within the observed scatter at low masses (grey circles; see the text for the compilation of individual detections), with EDGE2 approaching the $M_\star – M_{H1}$ relation of a stacked mass-complete sample at higher masses (grey diamonds; Scholte et al. 2025).</figcaption></figure>

halo masses. A way to more strongly discriminate between feedback models is thus to increase the precision of $\langle[\text{Fe/H}]\rangle$ measurements in dwarf galaxies. This will occur in the forthcoming years as more stellar abundances per dwarf galaxy become available (e.g. Skúladóttir et al. 2023).

Another promising route is to combine $M_V$–$\langle[\text{Fe/H}]\rangle$ with more observables of the star formation and feedback cycle. Given the sensitivity of iron abundances, further ratios beyond chemical elements are likely to provide complementary constraints on metal production and retention. We will use the extended number of chemical elements tracked in EDGE2 to explore this in future work. Next, we focus on the gas contents and gas-phase metallicity of our simulated dwarf galaxies.

## 4.3 The gas content of dwarf galaxies

Fig. 8 shows the $M_\star – M_{H1}$ relation for our simulated dwarf galaxies, where $M_{H1}$ is the total H1 mass within $r_{200}$ (see Rey et al. 2022 for how we derive $M_{H1}$ in EDGE1, while natively tracks out-of-equilibrium H1 fractions). Violins in Fig. 8 show the distribution of $M_{H1}$ over the last 4 billion years (symbols and lines showing the medians and 16–84 confidence intervals, respectively). This is long enough to average over the time variability due to stellar feedback, but short enough that $M_\star$ does not vary significantly (e.g. Rey et al. 2022).

Both models predict a similar structure for the $M_\star – M_{H1}$ relation, with the most distinct feature being a bimodality between gas-deficient dwarfs with vanishing H1 contents (upper limits at the bottom) overlapping in $M_\star$ with gas-rich systems. This bimodality arises from the interplay between the UVB and stellar feedback which regulate the gas and H1 contents, and the diversity of possible mass-growth histories in a ΛCDM universe (Rey et al. 2022; see also
e.g. Benitez-Llambay & Frenk 2020; Benitez-Llambay & Fumagalli 2021; Kim et al. 2024).

For dwarf galaxies with measurable H I contents, both models produce dwarf galaxies within the observational scatter from a compilation of individual gas-rich field dwarfs (McConnachie 2012; Cole et al. 2014; McQuinn et al. 2015, 2020, 2021; Sand et al. 2015; Adams & Oosterloo 2018; Brunker et al. 2019; Janesh et al. 2019; Hargis et al. 2020; Bennet et al. 2022; Xu et al. 2023; grey circles). In the case of the EDGE2 model, this agreement is striking for $M_{\star} \approx 10^7 M_{\odot}$ where simulated dwarfs cluster around the $M_{\star} - M_{H1}$ relation from Scholte et al. (2025) that combines DESI and ALFALFA (grey diamonds) to stack a mass-complete sample. Note that the error bars on the stack show the error on the median, not the expected population scatter which can be read from individual measurements.

Across the suites, EDGE2 galaxies are on average more H I-rich than EDGE1 galaxies at similar $M_{\star}$, as expected since radiative feedback dampens galactic outflows (Fig. 4). Similarly, both models predict substantial variability (violins), but $M_{H1}$ is more stable over time in the EDGE2 model (the extent of the 16–84 confidence interval is 53 per cent smaller on average). And even more notable difference is the lack of H I-deficient, higher mass dwarf galaxies ($M_{\star} \ge 10^6 M_{\odot}$) in EDGE2. Runaway star formation and highly clustered SN feedback during mergers can vacate the whole H I reservoir and self-quench a dwarf galaxy, while a more gentle, radiative regulation of star formation allows them to retain more gas.

We conclude that the exact shape of the bimodality at low $M_{\star}$, the median $M_{H1}$ at a given $M_{\star}$ and the scatter around this median due to time variability are all sensitive to the efficiency with which dwarf galaxies drive galactic outflows. This is promising for future comparisons with forthcoming radio surveys (e.g. Wallaby, Koribalski et al. 2020; Apertif-Medium deep, van Cappellen et al. 2022). But the fact that the structure and the key features of the $M_{\star} - M_{H1}$ relation robustly emerge from our modelling, despite large changes in the cooling and heating physics of the EDGE model, highlight that constraining power will only be unlocked by leveraging large populations of dwarf galaxies. This requires developing a new generation of semi-analytical models building on the insights presented here to generate large statistical samples of H I dwarfs, which we will present in future work (Hutton et al. in preparation).

## 4.4 The gas-phase mass–metallicity relation

Fig. 9 shows the $M_{\star} - \langle 12 + \log(O/H) \rangle$ relation at $z = 0$, where $\langle 12 + \log(O/H) \rangle$ is the mass-weighted average oxygen metallicity in the gas phase. To compute $\langle 12 + \log(O/H) \rangle$, we select gas within $2r_{1/2,v}$ of dwarf galaxies that have formed new stars in the last 4 billion years, as measured data of oxygen metallicity almost invariably rely on ionized emission lines in H II regions. We do not account for dust depletion which should be small at our low metallicities.

Observed data (grey circles) shows measurements of individual $z = 0$ dwarf galaxies compiled by Yates et al. (2020) and Breneman et al. (in preparation). All these points measure the electron temperature $T_e$ from auroral lines to derive $\langle 12 + \log(O/H) \rangle$. We omit their error bars for visual clarity. We also show the stack of emission-line-selected dwarf galaxies from DESI (Scholte et al. 2025; pentagons) with metallicities derived from strong-line calibration (Nakajima et al. 2022).

Both EDGE1 and EDGE2 are broadly within the scatter of individual measurements of gas-phase metallicities and track the slope of the relation. On average, EDGE2 dwarfs are slightly more oxygen rich ($\approx 0.1$ dex) at similar stellar masses than EDGE1 dwarfs. This much smaller upshift than for stellar $\langle [Fe/H] \rangle$ ($\approx 0.5$ dex) is readily

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Stellar mass versus gas-phase oxygen metallicity for simulated dwarf galaxies that have formed new stars within the last 4 Myr. Both models track the slope of the mass–metallicity relation and are broadly compatible with individual measurements of metal-poor, low-mass dwarfs (grey). EDGE2 dwarfs are slightly more metal-rich at a given stellar mass and populate the upper end of the observational scatter. But these effects are smaller than uncertainties arising from measuring gas-phase metallicities with different calibration methods (strong-line calibration as diamonds, direct method as circles).</figcaption></figure>

explained by the difference in chemical elements used for the observable. While the weaker galactic outflows in EDGE2 increase oxygen retention in the ISM (driving $\langle 12 + \log(O/H) \rangle$ higher), oxygen production per stellar mass formed is more than halved with our updated CCSNe yields (Fig. 1), driving $\langle 12 + \log(O/H) \rangle$ lower. When combined, this leads to a small change in the overall $\langle 12 + \log(O/H) \rangle$ of our galaxies.

The EDGE2 model predicts one example of a very low-mass $M_{\star} \le 10^5 M_{\odot}$, very metal-poor, star-forming dwarf galaxy (bottom-left corner). This object ('Halo 600: GM Later') is quenched by cosmic reionization early on, and achieves a high-enough dynamical mass to re-ignite star formation particularly late ($z = 0.2$, Rey et al. 2020). The delayed growth drives down $M_{\star}$ and $\langle 12 + \log(O/H) \rangle$ along the slope of the mass–metallicity relation. This highlights the possibility for low-mass dwarf galaxies even more oxygen-deficient than currently known forming through such rare assembly histories.

We conclude that, like for the stellar-phase mass metallicity relation (Fig. 7), the normalization of the gas-phase mass metallicity is sensitive to the galaxy formation assumptions made. However, given the implementation choices between EDGE1 and EDGE2, the shifts observed ($\approx 0.1$ dex) are smaller than systematic uncertainties in metallicity measurements. This is illustrated by the $\approx 0.3$ dex offset between individual and stacked measurements (grey circles and diamonds). Such offset is of the magnitude expected given their distinct metallicity calibration (strong line method validated with ‘semidirect’ observations, versus ‘semidirect’ and ‘direct’ auroral measurements; see e.g. Kewley & Ellison 2008; Yates et al. 2020 for a discussion). And shifts in $\langle 12 + \log(O/H) \rangle$ of $\approx 0.1$ dex are also well within the uncertainties of comparing direct-method measurements to simulation values (e.g. Cameron, Katz & Rey
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Stellar mass and metallicity growth over cosmic time (left) and at $z = 0$ for multiple versions of the same dwarf galaxy ('Halo 605') varying input parameters. The increased iron production between EDGE1 and EDGE2 (see Fig. 1) drives a ≈ 0.4 dex increase in $\langle[\text{Fe/H}]\rangle$ at $z = 0$. Further adding radiative feedback increases $\langle[\text{Fe/H}]\rangle$ by another ≈ 0.5 dex by enhancing iron retention in the ISM. Varying the SED and parameters of the star formation algorithm has a smaller overall effect, but modulates the clustering of star formation and the ability to drive more or less powerful outflows. Intrinsic stochasticity due to chaotic noise in the simulations is small compared to these shifts (envelopes).</figcaption></figure>

2023). A more detailed exploration is thus warranted to establish the precise constraining power of the gas-phase mass–metallicity relation on feedback models in the dwarf galaxy regime.

# 5 NUMERICAL LIMITATIONS AND REMAINING UNCERTAINTIES IN THE EDGE2 MODEL

The comparison between new and previous generation of EDGE simulations has revealed impressive convergence in observables (Section 4), but fundamentally different ISM and outflow structures (Section 3). In this section, we discuss the sensitivity of these comparisons to key uncertainties, discuss numerical limitations of our radiative simulations, and highlight key systematics yet to be quantified when modelling small dwarf galaxies.

## 5.1 Direct numerical tests

To explore uncertainties in our modelling, we resimulate multiple times a dwarf galaxy midway through the mass range of our sample ('Halo 605'; $M_* \approx 10^6 M_\odot$). This dwarf survives reionization and maintains star formation for a *Hubble* time, providing a long baseline to compare between model inputs. Fig. 10 shows how $M_*$ and $\langle[\text{Fe/H}]\rangle$ respond to variations in key numerical and physical parameters that we detail next.

### 5.1.1 Chaotic noise and stochasticity

We start by quantifying the intrinsic noise level and stochasticity in our model. The chaotic nature of galaxy formation combined with finite numerical precision can lead to divergent evolution from the

same initial conditions (see e.g. discussion in Keller et al. 2019; Genel et al. 2019).

To quantify the magnitude of this noise term, we resimulate our test galaxy twice, each time seeding a different truncation error. We achieve this by resimulating the same initial condition on a different number of cores, which changes the order of arithmetic operations and their truncation errors when communicating across a supercomputer’s network.

The envelopes (left-hand panels) around the fiducial EDGE2 run (red) in Fig. 10 show the extent of the differences in $M_*$ and $\langle[\text{Fe/H}]\rangle$ spanned by these two additional stochastic resimulations. At early times, when few stellar particles are present, stochasticity effects can be substantial, but the different simulations have converged with one another by $z \approx 6$.

Averaging over the three simulations at $z = 0$, we find $M_* = 2.97^{+0.07}_{-0.09} \times 10^6 M_\odot$ and $\langle[\text{Fe/H}]\rangle = -1.41^{+0.04}_{-0.11}$. Such differences are small compared to the factor-of-a-few changes going from the EDGE1 to the EDGE2 model. They also remain small compared to other input variations that we explore next.

Admittedly, three simulations is a limited sample size to make statistical statements. The numerical costs of these simulations make a wider exploration impractical, but we note that the magnitude of the scatter in $M_*$ is comparable to previous findings in the EDGE1 context (Pontzen et al. 2021). We conclude that chaotic stochasticity, while present, plays a small role in the trends and shifts in observables discussed in this paper.

### 5.1.2 Radiative feedback versus other inputs

In addition to radiative feedback, the EDGE2 model makes significant updates to the chemical enrichment modelling and SN feedback budget (Fig. 1). To isolate the effect of these changes on $\langle[\text{Fe/H}]\rangle$,
we resimulate our test dwarf galaxy in the EDGE2 model turning off non-equilibrium cooling and radiative feedback.

Comparing this new run to the original EDGE1 version in Fig. 10 (blue against purple, respectively), we find a ≈ 2× increase in $M_*$ and ≈ 0.4 dex increase in $\langle[\text{Fe/H}]\rangle$. Conversely, contrasting runs with and without radiative feedback (purple versus red, respectively), $M_*$ is decreased by ≈ 20 per cent, while $\langle[\text{Fe/H}]\rangle$ further increases by another ≈ 0.5 dex. All of these shifts are large compared to chaotic stochasticity (Section 5.1.1).

To summarize, our updates to the IMF and to the chemical enrichment modelling between EDGE1 and EDGE2 drives roughly half of the total increase in $\langle[\text{Fe/H}]\rangle$ at $z=0$. The other half is due to radiative feedback promoting metal retention through weaker outflows. It is thus clear that both of these inputs play an important role in setting the normalization of the stellar mass-metallicity relation. In the future, we will provide a much more detailed exploration of CCSNe yield models and further variations of the SNeIa inputs (Andersson et al., in preparation).

### 5.1.3 The spectral energy distribution

Another key input to radiative feedback is the choice of the SED for stellar populations. Fig. 10 shows the response of our test dwarf galaxy when swapping from our fiducial SED (Bruzual & Charlot 2003, red) to a harder and more ionizing BPASS2.2 SED (Stanway et al. 2016; orange) that takes into account binary populations of massive stars.

As expected with increasing the hardness and number of ionizing photons per stellar population, the dwarf galaxy is ≈ 30 per cent fainter overall and ≈ 0.25 dex more metal poor. These shifts are more modest than the other numerical choices we have explored, but remain larger than the intrinsic stochasticity. We conclude that the choice of SED is a subdominant parameter in setting the stellar masses and metallicities. At present, we therefore settle on the fiducial SED for the EDGE2 suite as a whole.

### 5.1.4 The coupling between photoionization feedback and the star formation algorithm

Fig. 3 shows that the EDGE2 ISM has a longer high-density tail, with significant gas staying above our density threshold for star formation. This stems from the combined conditions that star-forming gas has to be dense and cold, with our fiducial setup restricting star-forming gas to have $T < T_* = 1000$ K (Section 2.4).

Fig. 10 shows two more resimulations of our test dwarf galaxy lowering and raising this temperature threshold, while keeping the density threshold the same at $300\ M_p\ \text{cm}^{-3}$. The input SED in all cases is taken to be BPASS2.2 (Section 5.1.3).

Lowering $T_*$ to 100 K (gold) increases the stellar mass by ≈ 60 per cent compared to the fiducial value with the same SED (orange), despite having made the conditions to form stars more restrictive. This is combined with an increase in $\langle[\text{Fe/H}]\rangle$ by ≈ 0.2 dex, indicating increased metal retention and weaker outflows. Effectively removing $T_*$ by increasing it to $10^7$ K (brown) lowers $M_*$ and $\langle[\text{Fe/H}]\rangle$.

We show in Appendix D that $T_*$ in fact controls the clustering of star formation and SNe events. Following a star formation event, the dense gas surrounding the newborn stellar particle is immediately heated to $T \ge T_* = 10^3$ K by radiative feedback. Any gas in the immediate vicinity of a star-forming event is thus prevented from forming new stars, introducing an effective limit on the clustering

of star formation. Fig. D2 shows how lowering $T_*$ compounds this effect, leading to star formation and SN events forming at ever-increasing densities and in turn reducing the efficiency of galactic outflows at clearing this high-density gas.

This behaviour reflects the difficulties of modelling stellar photoheating, even at the high resolution of these simulations. If all processes around H II regions are captured, a star formation site would remain a mixture of cold gas below $T_*$, allowed to form stars, and of warm photoheated gas above $T_*$ prevented from star formation. If, however, resolution is limited, this multiphase star formation site is numerically mixed into a partially neutral, lukewarm phase. Restricting star formation to cold gas then unphysically shuts down correlated star formation. As shown in Appendix B, our simulations are in this regime of partially resolved radiative feedback, for which removing $T_*$ provides a way to rebuild correlation between SF events that have been spuriously erased.

### 5.1.5 Convergence with dark matter resolution

Using a similar EDGE setup with non-equilibrium cooling and radiative transfer, Agertz et al. (2020) demonstrated that our fiducial dark matter resolution ($m_{\text{DM}} \approx 950\ M_{\odot}$) resolves well the cooling and ISM of small dark matter haloes. Resimulating a low-mass object with eight times more dark matter particles while keeping the spatial resolution fixed, they found well converged stellar masses, sizes, and iron metallicities (compare ‘Fiducial+RT’ to ‘Hires + RT’ in their table 2).

We repeat this experiment with the updated EDGE2 model, resimulating one of our low-mass object, ‘Halo 1459 (earlier)’, improving $m_{\text{DM}}$ from 940 to $118\ M_{\odot}$ and keeping $\Delta x = 3$ pc. At $z=0$, $M_*$ changes from $7.4 \times 10^4$ to $7.9 \times 10^4\ M_{\odot}$, $r_{1/2,\text{V}}$ from 230 to 190 pc, and $\langle[\text{Fe/H}]\rangle$ from −2.03 to −2.11. These shifts are comparable in magnitude to those induced by chaotic stochasticity (Section 5.1.1) and we conclude that the observables presented in this work are robust to the choice of dark matter resolution.

## 5.2 The impact of partially resolved Strömgren spheres

In Appendix B, we show that at our resolution ($\Delta x = 3$ pc) the Strömgren radii, $R_S$, of photoionized H II regions from radiative feedback events are only partially resolved. This is a key numerical limitation of our simulations, which has a number of consequences we discuss here.

First, an unresolved Strömgren radius around a stellar cluster will mix two gas regions, one warm and ionized, one cold and neutral, into a single element that is lukewarm and partially ionized. Such out-of-equilibrium dense, warm, and mostly-neutral gas cools very efficiently through collisional excitation, generating a source of numerically induced cooling that is yet to be fully understood (e.g. Deng et al. 2024a).

Furthermore, this numerically mixed gas from unresolved Strömgren radii interacts with the temperature threshold for star formation (Section 5.1.4), preventing star formation in larger regions of space than should physically be allowed. This then suppresses the clustering of star formation and subsequent CCSNe feedback. In turn, this may mean that EDGE2 outflows are somewhat too weak, but we stress that the impressive agreement with observed scaling relations (Section 4) makes it unclear whether this is a significant issue. In future work, we will quantitatively compare the characteristics of our outflows to observations to obtain an orthogonal constraint (e.g. McQuinn et al. 2019; Xu et al. 2022; Marasco et al. 2023).
How to best address the numerical losses from under-resolved Strömgren spheres is a key open question for the next generation of radiation-hydrodynamics dwarf galaxy simulations.

An obvious solution is to simply increase numerical resolution until $R_S$ is resolved, providing a robust test of the physical-versus-numerical reduction in the strength of galactic outflows. Such an endeavour will be computationally demanding, but could already be achieved using super-Lagrangian strategies around star formation events similar to those implemented for black hole accretion (e.g. Curtis & Sijacki 2015; Beckmann, Slyz & Devriendt 2018; Anglés-Alcázar et al. 2021). However, given that numerical resolution is already at ≈ 3 pc across the ISM, going further will likely require algorithmic improvements to, for example explicitly sample the IMF in star formation events and inject radiative feedback star by star rather than as a population average as is done here.

An alternative approach could be to keep the resolution unchanged, but revise subgrid algorithms. For example, one can probabilistically draw the ionization fraction of H II regions when they are under-resolved (Hopkins et al. 2022), or allow stars to form at either lower densities or in more massive particles to boost their ionizing power and thus $R_S$ (e.g. Katz et al. 2024). Future studies quantifying the effects of these assumptions on dwarf galaxy properties will be crucial.

Lastly, under-resolved Strömgren spheres also have an important impact when comparing simulated and observed dwarf galaxies. Since cooling is efficient in the spurious warm-dense, partially neutral ISM phase, it leads to much brighter emission lines luminosities compared to what should be the true luminosity of the dwarf galaxy. This requires corrections to obtain accurate estimates (e.g. Katz et al. 2023; Ejdetjärn et al. 2024) and will be a key problem to solve when comparing EDGE2 data with observed emission lines.

## 5.3 Untested theoretical uncertainties

Beyond the numerical limitations and parameter exploration discussed above, there remain several astrophysical inputs needed for cosmological zoom simulations whose impact we are yet to explore.

For example, the choice of the IMF, of the UVB, and of the chemical yields all have consequences for the cooling, heating, and metallicity balance of the ISM. We have not yet been able to test all of these independently within the new EDGE2 model, but many of these uncertainties have been explored in previous iterations of the EDGE model. In particular, Rey et al. (2020) discuss the consequences of varying the redshift at which cosmic reionization happens (see also Benitez-Llambay & Frenk 2020; Katz et al. 2020), Prgomet et al. (2022) investigate the impact of a metallicity-dependent IMF, while Andersson et al. (2025) show convergence tests against a new feedback budget sourced by individual stars. Future work will also present a thorough exploration of the impact of yield inputs on EDGE galaxy chemistry (Andersson et al., in preparation).

And despite the milestone advance of including non-local stellar radiative input, other galaxy formation processes remain missing from this breed of EDGE dwarf galaxies. Exotic stellar evolution tracks (e.g. hypernovae, pair-instability SNe, variable-energy SNe, and metal-free primordial stars) can have dramatic effects on such small dwarf galaxies (e.g. Jeon, Besla & Bromm 2017; Jeon et al. 2021; Gutcke et al. 2022; Sanati et al. 2023) and leave distinct traces in their chemistry. The growing abundance of evidence of active galactic nuclei in dwarf galaxies (see e.g. Reines, Greene & Geha 2013; Reines et al. 2020; Burke et al. 2022; Davis et al. 2022; Mezcuá & Sánchez 2024 and summary in Wasleske & Baldassare 2024) and their associated outflows (e.g. Liu et al. 2020, 2024) calls

for the inclusion of black hole processes at this mass scale, although different implementations do not converge on their importance over a Hubble time (e.g. Koudmani et al. 2021, 2022; Sharma et al. 2023; Arjona-Galvez et al. 2024; De Almeida et al. 2024). Similarly, cosmic ray feedback provides a different avenue to regulate star formation, with its efficiency in faint dwarf galaxies now starting to be quantified (e.g. Martin-Alvarez et al. 2023). Future extensions of the EDGE project will explore these avenues.

## 6 SUMMARY AND CONCLUSION

We present new EDGE radiation-hydrodynamics cosmological zoomed simulations covering dwarf galaxy formation from the ultrafaint to the dwarf irregular regime ($10^9 \le M_{200} \le 10^{10} M_\odot$; $10^4 \le M_* \le 10^8 M_\odot$ at $z=0$). Our combination of a uniformly high resolution (≈ 3 pc, $m_{DM} \approx 950 M_\odot$), a large sample size (15 galaxies) and a detailed stellar feedback modelling (resolved SN feedback and explicit radiative feedback) is unprecedented. Leveraging previous-generation EDGE simulations with similar resolution but without radiative feedback, we systematically compare the response of dwarf galaxy observables to this change in physical modelling.

The addition of radiative feedback leads to a fundamentally different ISM structure (Fig. 3) in which star formation and gas accretion are increasingly regulated by radiative heating from massive stars, rather than mechanical removal from the halo centre. This change in regulation mode leads to a strong suppression of the mass outflow rates in dwarf galaxies over their cosmological history (Fig. 4). These findings reaffirm that radiative feedback reduces galactic outflows in dwarf galaxies, extending previous results (Agertz et al. 2020 for a single cosmological dwarf and e.g. Emerick et al. 2018; Smith et al. 2021; Deng et al. 2024b; Andersson et al. 2024 for isolated non-cosmological examples) to a 3-dex range in $M_*$. The extended mass range also shows that the suppression is the most pronounced in low-mass objects ($M_* \le 10^6 M_\odot$) for which shallow gravitational potential wells increases sensitivity to radiative feedback (Fig. 5).

Despite the significant change in outflow behaviours, stellar masses at $z=0$ are converged to ≈ 40 per cent on average across the suite, a shift well within uncertainties (Fig. 2). In fact, most scaling relations of dwarf galaxy observables cannot distinguish between these two models. The $M_V-r_{1/2,V}$ relation (Fig. 6), the $M_* - M_{H1}$ relation (Fig. 8), and the $M_*-\langle 12 + \log(O/H) \rangle$ relation (Fig. 9) all showcase small shifts compared to the width of the observational scatter and its error bars.

Only the $M_V-\langle [Fe/H] \rangle$ presents a strong response to our model changes, with both models tracking the observed slope of the mass-metallicity relation but respectively scattering at the lower and upper ends of the observed scatter (Fig. 7). Higher $\langle [Fe/H] \rangle$ in the new simulations are driven by the combined update to the chemical modelling – yield updates increasing iron production – and the weaker galactic outflows – that increase iron retention in the ISM (Fig. 10).

These results confirm the expectation from Agertz et al. (2020) that the $M_V-\langle [Fe/H] \rangle$ relation is the most sensitive probe of stellar and ISM physics in dwarf galaxies. There, we used a single dwarf galaxy ran with and without radiative transfer. This paper extends these results to a much wider range of stellar and halo masses, showing that the slope of the $M_V-\langle [Fe/H] \rangle$ readily emerges in both models but that the normalization of the relation directly reflects the combination of modelling choices that control iron production (e.g. CCSNe yields) and retention (e.g. outflow strength).

A promising way to distinguish between these two effects will be to leverage the wealth of new chemical data from spectroscopic surveys
(e.g. > 50 000 chemical abundances of dwarf galaxy stars with 4MOST; Skúladóttir et al. 2023). Unpicking the relative abundances of each element and their scatter will provide stronger constraints on the chemical yields of stars at low metallicities. In parallel, direct observations of dwarf galaxies' outflows, their velocities, and their mass-loading factors is now possible thanks to the advancements in integral-field-unit technology (e.g. McQuinn et al. 2019; Xu et al. 2022; Marasco et al. 2023). A careful comparison of these emission-line-derived observations to simulations like the EDGE2 suite should allow us to constrain the physical mechanisms driving galactic outflows in dwarf galaxies. Furthermore, additional dynamical scaling relations (e.g. Vandenbroucke, Verbeke & De Rijcke 2016; Di Paolo, Salucci & Fontaine 2019; Romeo 2020) and new observables such as the distribution of ionized metal absorbers in the CGM of dwarf irregulars (e.g. Zheng et al. 2019, 2020, 2024) can inform us on the baryon cycle of dwarf galaxies and offer orthogonal constraints on the accretion-outflow cycle to be explored in the coming years.

The robust emergence of several dwarf galaxy scaling relations and the relative convergence of our results signal that cosmological simulations of dwarf galaxies have entered an era of precision. Even with a committed theoretical variation (the inclusion of a new feedback channel) which changes by order of magnitudes the outflow behaviour of simulated dwarf galaxies, much of their general characteristics remain broadly compatible with observational data. None the less, there remains clear numerical limitations to our simulations. While SN feedback is now well resolved and converged, radiative feedback remains less accurately captured (Section 5). Moreover, other theoretical uncertainties that impact the strength of galactic outflows in dwarf galaxies remain to be quantified, for example the importance of active galactic nucleus or cosmic rays in this regime. Addressing these challenges will be the focus of future works from the EDGE collaboration. But ultimately, the convergence of basic scaling relations makes for an exciting moment, providing us with a solid foundation to differentiate between these mechanisms and interpret new dwarf galaxy data coming in the next decade.

## AUTHOR CONTRIBUTIONS

The main roles of the authors were, using the CRediT (Contribution Roles Taxonomy) system:<fnref n="2" />

MR: Conceptualization; Data curation; Formal analysis; Funding acquisition; Investigation; Writing – original draft. ET: Data curation; Investigation; Writing – review and editing. EG: Data curation; Investigation; Writing – review and editing. SK: Data curation; Investigation; Writing – review and editing. EA: Investigation; Methodology; Software; Writing – review and editing. AP: Conceptualization; Funding acquisition; Software; Writing – review and editing. OA: Conceptualization; Methodology; Software; Writing – review and editing. JR: Conceptualization; Project administration; Resources; Writing – review and editing. CC: Methodology; Software; Writing – review and editing. RY: Writing – review and editing. MO: Data Curation; Investigation; Writing – review and editing. DS: Data curation; Writing – review and editing. AS: Conceptualization; Writing – review and editing. JB: Data curation; Writing – review and editing. KQ: Writing – review and editing. CM: Writing – review and editing. PD: Writing – review and editing.

## ACKNOWLEDGEMENTS

We thank the referee for constructive criticism that improved the quality of the paper. MR thanks Astha Agarwal, Harley Katz, Taysun Kimm, Joakim Rosdahl, and Mahsa Sanati for useful discussions and comments during the preparation of this manuscript. MR acknowledges support from the Beecroft Fellowship funded by Adrian Beecroft. ET acknowledges the UKRI Science and Technology Facilities Council (STFC) for support (grant ST/V50712X/1). OA acknowledges support from the Knut and Alice Wallenberg Foundation, the Swedish Research Council (grant 2019-04659), the Swedish National Space Agency (SNSA Dnr 2023-00164), and the LMK Foundation. JIR would like to thank the STFC for support from grants ST/Y002865/1 and ST/Y002857/1. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 818085 GMGalaxies. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. The authors acknowledge the use of the UCL Grace High Performance Computing Facility, the Surrey Eureka supercomputer facility, and their associated support services. This work was partially supported by the UCL Cosmoparticle Initiative.

We thank the developers and maintainers of PYNBODY (Pontzen et al. 2013), TANGOS (Pontzen & Tremmel 2018), NUMPY (van der Walt, Colbert & Varoquaux 2011; Harris et al. 2020), SciPy (Virtanen et al. 2020), JUPYTER (Ragan-Kelley et al. 2014), MATPLOTLIB (Hunter 2007), the Astrophysics Data Service, and the arXiv preprint repository for providing open-source softwares that were used extensively in this work.

## DATA AVAILABILITY

The data underlying this article will be shared upon reasonable request to the corresponding author.

## REFERENCES

Adams E. A. K., Oosterloo T. A., 2018, A&A, 612, A26

Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391

Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

Agertz O. et al., 2020, MNRAS, 491, 1656

Agertz O. et al., 2021, MNRAS, 503, 5826

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

Andersson E. P., Rey M. P., Pontzen A., Cadiou C., Agertz O., Read J. I., Martin N. F., 2025, ApJ, 978, 129

Anglés-Alcázar D. et al., 2021, ApJ, 917, 53

Arjona-Galvez E., Di Cintio A., Grand R. J. J., 2024, A&A, 690, A286

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Aubert D., Teyssier R., 2010, ApJ, 724, 244

Beckmann R. S., Slyz A., Devriendt J., 2018, MNRAS, 478, 995

Benítez-Llambay A., Navarro J. F., Abadi M. G., Gottlöber S., Yepes G., Hoffman Y., Steinmetz M., 2015, MNRAS, 450, 4207

Benitez-Llambay A., Frenk C., 2020, MNRAS, 498, 4887

Benitez-Llambay A., Fumagalli M., 2021, ApJ, 921, L9

Bennet P. et al., 2022, ApJ, 924, 98

Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, ApJ, 500, 342

Brauer K. et al., 2025, ApJ, 980, 41

<footnote n="2">²https://authorservices.wiley.com/author-resources/Journal-Authors/open-access/credit.html</footnote>
Bressan A., Marigo P., Girardi Léo., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

Bruce J., Li T. S., Pace A. B., Heiger M., Song Y.-Y., Simon J. D., 2023, ApJ, 950, 167

Brunker S. W. et al., 2019, AJ, 157, 76

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

Burke C. J. et al., 2022, MNRAS, 516, 2736

Buttry R. et al., 2022, MNRAS, 514, 1706

Cadiou C., Dubois Y., Pichon C., 2019, A&A, 621, A96

Cadiou C., Pontzen A., Peiris H. V., 2021, MNRAS, 502, 5480

Cameron A. J., Katz H., Rey M. P., 2023, MNRAS, 552, L89

Cerny W. et al., 2023a, ApJ, 942, 111

Cerny W. et al., 2023b, ApJ, 953, 1

Cerny W. et al., 2023c, ApJ, 953, L21

Chabrier G., 2003, PASP, 115, 763

Charles E. J. E. et al., 2023, MNRAS, 521, 3527

Chisholm J., Tremonti C. A., Leitherer C., Chen Y., 2017, MNRAS, 469, 4831

Chiti A. et al., 2021, Nat. Ast., 5, 392

Chiti A., Simon J. D., Frebel A., Pace A. B., Ji A. P., Li T. S., 2022, ApJ, 939, 41

Christensen C. R., Davé R., Governato F., Pontzen A., Brooks A., Munshi F., Quinn T., Wadsley J., 2016, ApJ, 824, 57

Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252

Cole A. A., Weisz D. R., Dolphin A. E., Skillman E. D., McConnachie A. W., Brooks A. M., Leaman R., 2014, ApJ, 795, 54

Collins M. L. M., Tollerud E. J., Rich R. M., Ibata R. A., Martin N. F., Chapman S. C., Gilbert K. M., Preston J., 2020, MNRAS, 491, 3496

Collins M. L. M. et al., 2021, MNRAS, 505, 5686

Collins M. L. M. et al., 2024, MNRAS, 528, 2614

Collins M. L. M., Read J. I., 2022, Nat. Ast., 6, 647

Côté B. et al., 2018, ApJ, 855, 99

Courty S., Alimi J. M., 2004, A&A, 416, 875

Curtis M., Sijacki D., 2015, MNRAS, 454, 3445

Davis F. et al., 2022, MNRAS, 511, 4109

De Almeida A. P., Mamon G. A., Dekel A., Lima Neto G. B., 2024, A&A, 687, A131

Deng Y., Li H., Kannan R., Smith A., Vogelsberger M., Bryan G. L., 2024a, MNRAS, 527, 478

Deng Y., Li H., Liu B., Kannan R., Smith A., Bryan G. L., 2024b, A&A, 691, A231

Di Paolo C., Salucci P., Fontaine J. P., 2019, ApJ, 873, 106

Efstathiou G., 1992, MNRAS, 256, 43P

Eisenstein D. J., Hut P., 1998, ApJ, 498, 137

Ejdetjäm T., Agertz O., Östlin G., Rey M. P., Renaud F., 2024, MNRAS, 534, 135

Emerick A., Bryan G. L., Mac Low M.-M., 2018, ApJ, 865, L22

Escala I. et al., 2018, MNRAS, 474, 2194

Faucher-Giguère C.-A., 2020, MNRAS, 493, 1614

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Ferland G. J. et al., 2017, Rev. Mex. Astron. Astrofis., 53, 385

Fitts A. et al., 2017, MNRAS, 471, 3547

Fritz T. K., Carrera R., Battaglia G., Taibi S., 2019, A&A, 623, A129

Genel S. et al., 2019, ApJ, 871, 21

Gil-Pons P., Doherty C. L., Gutiérrez J. L., Siess L., Campbell S. W., Lau H. B., Lattanzio J. C., 2018, Publ. Astron. Soc. Aust., 35, e038

Gnedin N. Y., 2000, ApJ, 542, 535

Gnedin N. Y., Abel T., 2001, New Astron., 6, 437

Go M. et al., 2025, ApJ, 986, 214

Goater A. et al., 2024, MNRAS, 527, 2403

Gray E. I. et al., 2025, MNRAS, 539, 1167

Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 2010, ApJ, 716, 510

Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093

Grisdale K., Agertz O., Renaud F., Romeo A. B., 2018, MNRAS, 479, 3167

Grisdale K., Agertz O., Renaud F., Romeo A. B., Devriendt J., Slyz A., 2019, MNRAS, 486, 5482

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

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

Haardt F., Madau P., 1996, ApJ, 461, 20

Hansen T. T., Simon J. D., Li T. S., Sharkey D., Ji A. P., Thompson I. B., Reggiani H. M., Galarza J. Y., 2024, ApJ, 968, 21

Hargis J. R. et al., 2020, ApJ, 888, 31

Harris C. R. et al., 2020, Nature, 585, 357

Heiger M. E. et al., 2024, ApJ, 961, 234

Herzog G., Benítez-Llambay A., Fumagalli M., 2023, MNRAS, 518, 6305

Hoeft M., Yepes G., Gottlöber S., Springel V., 2006, MNRAS, 371, 401

Homma D. et al., 2019, PASJ, 71, 94

Homma D. et al., 2024, PASJ, 76, 733

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Hopkins P. F. et al., 2022, MNRAS, 519, 3154

Hopkins P. F., 2025, Open J. Astrophys., 8, 44

Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Jaacks J., Thompson R., Finkelstein S. L., Bromm V., 2018, MNRAS, 475, 4396

Janesh W., Rhode K. L., Salzer J. J., Janowiecki S., Adams E. A. K., Haynes M. P., Giovanelli R., Cannon J. M., 2019, AJ, 157, 183

Jenkins S. A., Li T. S., Pace A. B., Ji A. P., Koposov S. E., Mutlu-Pakdil B., 2021, ApJ, 920, 92

Jeon M., Besla G., Bromm V., 2017, ApJ, 848, 85

Jeon M., Bromm V., Besla G., Yoon J., Choi Y., 2021, MNRAS, 502, 1

Ji A. P., Frebel A., Simon J. D., Geha M., 2016a, ApJ, 817, 41

Ji A. P., Frebel A., Simon J. D., Chiti A., 2016b, ApJ, 830, 93

Ji A. P., Frebel A., Ezzeddine R., Casey A. R., 2016c, ApJ, 832, L3

Ji A. P., Simon J. D., Frebel A., Venn K. A., Hansen T. T., 2019, ApJ, 870, 83

Ji A. P. et al., 2020, ApJ, 889, 27

Ji A. P. et al., 2021, ApJ, 921, 32

Jones M. G. et al., 2023, ApJ, 957, L5

Jones M. G. et al., 2024, ApJ, 971, L37

Katz H. et al., 2020, MNRAS, 494, 2200

Katz H. et al., 2023, Open J. Astrophys., 6, 44

Katz H., Rey M. P., Cadiou C., Kimm T., Agertz O., 2024, preprint (arXiv:2411.07282)

Keller B. W., Wadsley J. W., Wang L., Kruijssen J. M. D., 2019, MNRAS, 482, 2244

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

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

Kim J.-h. et al., 2014, ApJS, 210, 14

Kim S. Y. et al., 2024, preprint (arXiv:2408.15214)

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kirby E. N., Cohen J. G., Guhathakurta P., Cheng L., Bullock J. S., Gallazzi A., 2013, ApJ, 779, 102

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

Kirby E. N., Rizzi L., Held E. V., Cohen J. G., Cole A. A., Manning E. M., Skillman E. D., Weisz D. R., 2017, ApJ, 834, 9

Kirby E. N., Gilbert K. M., Escala I., Wojno J., Guhathakurta P., Majewski S. R., Beaton R. L., 2020, AJ, 159, 46

Koribalski B. S. et al., 2020, Astrophys. Space Sci., 365, 118

Koudmani S., Henden N. A., Sijacki D., 2021, MNRAS, 503, 3568

Koudmani S., Sijacki D., Smith M. C., 2022, MNRAS, 516, 2112

Kravtsov A., Manwadkar V., 2022, MNRAS, 514, 2667

Kroupa P., 2001, MNRAS, 322, 231

Kvasova K., Kirby E. N., Beaton R. L., 2024, ApJ, 972, 180

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

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

Li T. S. et al., 2017, ApJ, 838, 8

Li T. S. et al., 2018, ApJ, 857, 145

LIGO and Virgo Scientific Collaboration 2017, Phys. Rev. Lett., 119, 161101

Limongi M., Chieffi A., 2018, ApJS, 237, 13

Limongi M., Roberti L., Chieffi A., Nomoto K., 2024, ApJS, 270, 28
Liu W., Veilleux S., Canalizo G., Rupke D. S. N., Manzano-King C. M., Bohn T. U V., 2020, ApJ, 905, 166

Liu W. et al., 2024, ApJ, 965, 152

Longeard N. et al., 2018, MNRAS, 480, 2609

Maoz D., Mannucci F., Brandt T. D., 2012, MNRAS, 426, 3282

Marasco A. et al., 2023, A&A, 670, A92

Martin-Alvarez S., Sijacki D., Haehnelt M. G., Farcy M., Dubois Y., Belokurov V., Rosdahl J., Lopez-Rodriguez E., 2023, MNRAS, 525, 3806

Martinez-Delgado D., Stein M., Pawlowski M. S., Makarov D., Makarova L., Donatiello G., Lang D., 2024, preprint (arXiv:2405.03769)

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

Mau S. et al., 2020, ApJ, 890, 136

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

McConnachie A. W., Venn K. A., 2020, AJ, 160, 124

McConnachie A. W., Higgs C. R., Thomas G. F., Venn K. A., Côté P., Battaglia G., Lewis G. F., 2021, MNRAS, 501, 2363

McNanna M. et al., 2024, ApJ, 961, 126

McQuinn K. B. W. et al., 2015, ApJ, 812, 158

McQuinn K. B. W. et al., 2021, ApJ, 918, 23

McQuinn K. B. W., Mao Y.-Y., Buckley M. R., Shih D., Cohen R. E., Dolphin A. E., 2023, ApJ, 944, 14

McQuinn Kristen. B. W., van Zee L., Skillman E. D., 2019, ApJ, 886, 74

McQuinn Kristen. B. W. et al., 2020, ApJ, 891, 181

McQuinn M., 2016, ARA&A, 54, 313

Mezcuá M., Sánchez H. D., 2024, MNRAS, 528, 5252

Mitchell P. D., Schaye J., Bower R. G., Crain R. A., 2020, MNRAS, 494, 3971

Monzon J. S., van den Bosch F. C., Mitra K., 2024, ApJ, 976, 197

Muni C., Pontzen A., Read J. I., Agertz O., Rey M. P., Taylor E., 2025, MNRAS, 536, 314

Munshi F., Brooks A. M., Christensen C., Applebaum E., Holley-Bockelmann K., Quinn T. R., Wadsley J., 2019, ApJ, 874, 40

Munshi F., Brooks A., Applebaum E., Christensen C., Sligh J. P., Quinn T., 2021, ApJ, 923, 35

Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691

Mutlu-Pakdil B. et al., 2021, ApJ, 918, 88

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

Nadler E. O. et al., 2020, ApJ, 893, 48

Naiman J. P. et al., 2018, MNRAS, 477, 1206

Nakajima K. et al., 2022, ApJS, 262, 3

Nebrin O., 2023, Res. Astron. Astrophys., 7, 90

Nelson D. et al., 2019, MNRAS, 490, 3234

Nguyen C. T. et al., 2022, A&A, 665, A126

Nickerson S., Teyssier R., Rosdahl J., 2018, MNRAS, 479, 3206

Noh Y., McQuinn M., 2014, MNRAS, 444, 503

O'Leary J. A., Steinwandel U. P., Moster B. P., Martin N., Naab T., 2023, MNRAS, 520, 897

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

Okamoto T., Gao L., Theuns T., 2008, MNRAS, 390, 920

Oñorbe J., Hennawi J. F., Lukić Z., 2017, ApJ, 837, 106

Orkney M. D. A. et al., 2021, MNRAS, 504, 3509

Orkney M. D. A. et al., 2022, MNRAS, 515, 185

Orkney M. D. A., Taylor E., Read J. I., Rey M. P., Pontzen A., Agertz O., Kim S. Y., Delorme M., 2023, MNRAS, 525, 3516

Pace A. B. et al., 2020, MNRAS, 495, 3022

Padoan P., Haugbølle T., Nordlund Å., 2012, ApJ, 759, L27

Pandya V. et al., 2021, MNRAS, 508, 2979

Pignatari M. et al., 2016, ApJS, 225, 24

Planck Collaboration XVI, 2014, A&A, 571, A16

Planck Collaboration VI, 2020, A&A, 641, A6

Pontzen A., Governato F., 2014, Nature, 506, 171

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Pontzen A., Tremmel M., Roth N., Peiris H. V., Saintonge A., Volonteri M., Quinn T., Governato F., 2017, MNRAS, 465, 547

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

Pontzen A., Rey M. P., Cadiou C., Agertz O., Teyssier R., Read J., Orkney M. D. A., 2021, MNRAS, 501, 1755

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Prgomet M., Rey M. P., Andersson E. P., Segovia Otero A., Agertz O., Renaud F., Pontzen A., Read J. I., 2022, MNRAS, 513, 2326

Ragan-Kelley M., Perez F., Granger B., Kluyver T., Ivanov P., Frederic J., Bussonnier M., 2014, Am. Geophys. Un., 2014, H44D

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

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Read J. I., Iorio G., Agertz O., Fraternali F., 2017, MNRAS, 467, 2019

Reines A. E., Greene J. E., Geha M., 2013, ApJ, 775, 116

Reines A. E., Condon J. J., Darling J., Greene J. E., 2020, ApJ, 888, 36

Revaz Y., Jablonka P., 2018, A&A, 616, A96

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

Rey M. P., Pontzen A., Saintonge A., 2019a, MNRAS, 485, 1906

Rey M. P., Pontzen A., Agertz O., Orkney M. D. A., Read J. I., Saintonge A., Pedersen C., 2019b, ApJ, 886, L3

Rey M. P., Pontzen A., Agertz O., Orkney M. D. A., Read J. I., Rosdahl J., 2020, MNRAS, 497, 1508

Rey M. P., Pontzen A., Agertz O., Orkney M. D. A., Read J. I., Saintonge A., Kim S. Y., Das P., 2022, MNRAS, 511, 5672

Rey M. P., Katz H. B., Cameron A. J., Devriendt J., Slyz A., 2024a, MNRAS, 528, 5412

Rey M. P. et al., 2024b, MNRAS, 529, 2379

Richstein H. et al., 2022, ApJ, 933, 217

Ritter C., Herwig F., Jones S., Pignatari M., Fryer C., Hirschi R., 2018, MNRAS, 480, 538

Romeo A. B., 2020, MNRAS, 491, 4843

Rosdahl J., Blaizot J., 2012, MNRAS, 423, 344

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J. et al., 2018, MNRAS, 479, 994

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

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

Sales L. V., Wetzel A., Fattahi A., 2022, Nat. Ast., 6, 897

Sanati M., Jeanquartier F., Revaz Y., Jablonka P., 2023, A&A, 669, A94

Sand D. J., Spekkens K., Crnojević D., Hargis J. R., Willman B., Strader J., Grillmair C. J., 2015, ApJ, 812, L13

Sand D. J. et al., 2022, ApJ, 935, L17

Schmidt M., 1959, ApJ, 129, 243

Scholte D. et al., 2025, MNRAS, 535, 2341

Seitenzahl I. R. et al., 2013, MNRAS, 429, 1156

Sharma R. S., Brooks A. M., Tremmel M., Bellovary J., Quinn T. R., 2023, ApJ, 957, 16

Simon J. D., 2018, ApJ, 863, 89

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

Simon J. D. et al., 2023, ApJ, 944, 43

Skúladóttir Á. et al., 2023, The Messenger, 190, 19

Smith B. D. et al., 2017, MNRAS, 466, 2217

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

Smith M. C., Bryan G. L., Somerville R. S., Hu C.-Y., Teyssier R., Burkhart B., Hernquist L., 2021, MNRAS, 506, 3882

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

Smith S. E. T. et al., 2024, ApJ, 961, 92

Snaith O. N., Park C., Kim J., Rosdahl J., 2018, MNRAS, 477, 983

Stanway E. R., Eldridge J. J., Becker G. D., 2016, MNRAS, 456, 485

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

Stopyra S., Peiris H. V., Pontzen A., 2021b, MNRAS, 500, 4173

Strömgren B., 1939, ApJ, 89, 526

Taibi S., Battaglia G., Rejkuba M., Leaman R., Kacharov N., Iorio G., Jablonka P., Zoccali M., 2020, A&A, 635, A152

Tan C. Y. et al., 2025, ApJ, 979, 176

Tarumi Y., Yoshida N., Frebel A., 2021, ApJ, 914, L10

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

Thielemann F. K., Nomoto K., Yokoi K., 1986, A&A, 158, 17

Thornton K., Gaudlitz M., Janka H. T., Steinmetz M., 1998, ApJ, 500, 95
Tollet E. et al., 2016, MNRAS, 456, 3542

Toro E. F., Spruce M., Speares W., 1994, Shock Waves, 4, 25

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

Torrealba G. et al., 2018, MNRAS, 475, 5085

Torrealba G. et al., 2019, MNRAS, 488, 2743

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

van Cappellen W. A. et al., 2022, A&A, 658, A146

van der Walt S., Colbert S. C., Varoquaux G., 2011, Comput. Sci. Eng., 13, 22

Vandenbroucke B., Verbeke R., De Rijcke S., 2016, MNRAS, 458, 912

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Visbal E., Bryan G. L., Haiman Z., 2020, ApJ, 897, 95

Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83

Wang Y., Nadler E. O., Mao Y.-Y., Adhikari S., Wechsler R. H., Behroozi P., 2021, ApJ, 915, 116

Wasleske E. J., Baldassare V. F., 2024, ApJ, 971, 68

Wheeler C. et al., 2019, MNRAS, 490, 4447

Wojno J., Gilbert K. M., Kirby E. N., Escala I., Beaton R. L., Tollerud E. J., Majewski S. R., Guhathakurta P., 2020, ApJ, 895, 78

Woosley S. E., Heger A., 2007, Phys. Rep., 442, 269

Wright A. C., Brooks A. M., Weisz D. R., Christensen C. R., 2019, MNRAS, 482, 1176

Xu J.-L., Zhu M., Yu N., Zhang C.-P., Liu X.-L., Ai M., Jiang P., 2023, ApJ, 944, L40

Xu Y. et al., 2022, ApJ, 929, 134

Yates R. M., Schady P., Chen T.-W., Schweyer T., Wiseman P., 2020, A&A, 634, A107

Zel'dovich Ya. B., 1970, A&A, 5, 84

Zheng Y. et al., 2019, MNRAS, 490, 467

Zheng Y. et al., 2024, ApJ, 960, 21

Zheng Y., Emerick A., Putman M. E., Werk J. K., Kirby E. N., Peek J. E. G., 2020, ApJ, 2020, 133

## APPENDIX A: SUMMARY OF EDGE1 AND EDGE2 PROPERTIES

Table A1 summarizes the properties of the simulated dwarf galaxies with the EDGE1 and EDGE2 models at $z = 0$.

<table><caption>**Table A1.** Properties of simulated dwarf galaxies with the EDGE1 and EDGE2 models (E1 and E2 columns, respectively) at $z = 0$ (see Section 4 for details on how each property is computed). Haloes with only an EDGE2 column are new initial conditions that do not have an EDGE1 counterpart (see Section 2.1). Two haloes were stopped at higher redshifts due to their computational cost (see Section 3.1) – their properties are indicated at the labelled redshift.</caption><thead><tr><th rowspan="2">Object</th><th colspan="2">M<sub>200</sub><br>(M<sub>☉</sub>)</th><th colspan="2">M<sub>⋆</sub> (M<sub>☉</sub>)</th><th colspan="2">M<sub>V</sub></th><th colspan="2">r<sub>1/2,V</sub> (kpc)</th><th colspan="2">M<sub>H1</sub> (M<sub>☉</sub>)</th><th colspan="2">([Fe/H])</th><th colspan="2">⟨12 + log(O/H)⟩</th></tr><tr><th>E1</th><th>E2</th><th>E1</th><th>E2</th><th>E1</th><th>E2</th><th>E1</th><th>E2</th><th>E1</th><th>E2</th><th>E1</th><th>E2</th><th>E1</th><th>E2</th></tr></thead><tbody><tr><td>Halo 1445</td><td>1.2 × 10<sup>9</sup></td><td>7.9 × 10<sup>4</sup></td><td>5.8 × 10<sup>4</sup></td><td>-6.35</td><td>-5.36</td><td>0.20</td><td>0.29</td><td>1.2 × 10<sup>1</sup></td><td>4.9 × 10<sup>1</sup></td><td>-2.58</td><td>-2.24</td><td>6.31</td><td>7.25</td></tr><tr><td>Halo 1459: GM earlier</td><td>1.4 × 10<sup>9</sup></td><td>4.9 × 10<sup>5</sup></td><td>4.4 × 10<sup>5</sup></td><td>-8.29</td><td>-7.50</td><td>0.16</td><td>0.22</td><td>1.3 × 10<sup>1</sup></td><td>4.3 × 10<sup>1</sup></td><td>-2.34</td><td>-1.78</td><td>6.63</td><td>8.39</td></tr><tr><td>Halo 1459</td><td>1.4 × 10<sup>9</sup></td><td>2.6 × 10<sup>5</sup></td><td>3.6 × 10<sup>5</sup></td><td>-7.62</td><td>-7.29</td><td>0.23</td><td>0.25</td><td>1.9 × 10<sup>1</sup></td><td>4.9 × 10<sup>1</sup></td><td>-2.40</td><td>-2.04</td><td>6.73</td><td>8.04</td></tr><tr><td>Halo 1459: GM later</td><td>1.4 × 10<sup>9</sup></td><td>1.2 × 10<sup>5</sup></td><td>7.5 × 10<sup>4</sup></td><td>-6.84</td><td>-5.60</td><td>0.31</td><td>0.28</td><td>1.7 × 10<sup>1</sup></td><td>7.4 × 10<sup>1</sup></td><td>-2.58</td><td>-2.03</td><td>6.55</td><td>7.25</td></tr><tr><td>Halo 1459: GM latest</td><td>1.4 × 10<sup>9</sup></td><td>4.5 × 10<sup>4</sup></td><td>1.5 × 10<sup>4</sup></td><td>-5.72</td><td>-3.94</td><td>1.00</td><td>2.63</td><td>2.0 × 10<sup>1</sup></td><td>1.0 × 10<sup>2</sup></td><td>-2.93</td><td>-2.53</td><td>6.44</td><td>6.35</td></tr><tr><td>Halo 600</td><td>3.4 × 10<sup>9</sup></td><td>5.1 × 10<sup>5</sup></td><td>6.2 × 10<sup>5</sup></td><td>-8.77</td><td>-9.83</td><td>0.23</td><td>0.30</td><td>5.7 × 10<sup>5</sup></td><td>1.8 × 10<sup>6</sup></td><td>-2.43</td><td>-1.49</td><td>7.35</td><td>7.37</td></tr><tr><td>Halo 600: GM later</td><td>3.2 × 10<sup>9</sup></td><td>3.9 × 10<sup>5</sup></td><td>6.4 × 10<sup>4</sup></td><td>-8.18</td><td>-7.00</td><td>0.38</td><td>0.10</td><td>2.1 × 10<sup>5</sup></td><td>3.9 × 10<sup>5</sup></td><td>-2.50</td><td>-2.70</td><td>6.64</td><td>6.80</td></tr><tr><td>Halo 605</td><td>3.3 × 10<sup>9</sup></td><td>1.7 × 10<sup>6</sup></td><td>2.9 × 10<sup>6</sup></td><td>-9.73</td><td>-9.91</td><td>0.26</td><td>0.32</td><td>3.2 × 10<sup>5</sup></td><td>7.5 × 10<sup>5</sup></td><td>-2.21</td><td>-1.35</td><td>7.30</td><td>7.73</td></tr><tr><td>Halo 624</td><td>2.5 × 10<sup>9</sup></td><td>5.9 × 10<sup>5</sup></td><td>4.5 × 10<sup>5</sup></td><td>-8.93</td><td>-8.12</td><td>0.38</td><td>0.63</td><td>2.3 × 10<sup>5</sup></td><td>3.1 × 10<sup>5</sup></td><td>-2.30</td><td>-1.98</td><td>6.68</td><td>7.27</td></tr><tr><td>Halo 624: GM higher<br>mass</td><td>3.8 × 10<sup>9</sup></td><td>1.5 × 10<sup>6</sup></td><td>1.9 × 10<sup>6</sup></td><td>-9.63</td><td>-9.83</td><td>0.31</td><td>0.25</td><td>2.3 × 10<sup>5</sup></td><td>1.6 × 10<sup>6</sup></td><td>-2.17</td><td>-1.60</td><td>7.42</td><td>7.50</td></tr><tr><td>Halo 383: GM earlier</td><td>5.8 × 10<sup>9</sup></td><td>4.4 × 10<sup>6</sup></td><td>2.1 × 10<sup>7</sup></td><td>-10.88</td><td>-12.19</td><td>0.33</td><td>0.49</td><td>3.8 × 10<sup>5</sup></td><td>3.6 × 10<sup>6</sup></td><td>-1.99</td><td>-0.96</td><td>7.70</td><td>8.03</td></tr><tr><td>Halo 383</td><td>5.7 × 10<sup>9</sup></td><td>3.2 × 10<sup>6</sup></td><td>9.4 × 10<sup>6</sup></td><td>-10.49</td><td>-11.85</td><td>0.42</td><td>0.57</td><td>5.8 × 10<sup>2</sup></td><td>7.7 × 10<sup>6</sup></td><td>-2.02</td><td>-1.10</td><td>7.10</td><td>7.74</td></tr><tr><td>Halo 383: stochastic</td><td>5.9 × 10<sup>9</sup></td><td>9.6 × 10<sup>6</sup></td><td>1.1 × 10<sup>7</sup></td><td>-11.76</td><td>-12.23</td><td>0.17</td><td>0.35</td><td>5.5 × 10<sup>1</sup></td><td>4.5 × 10<sup>6</sup></td><td>-1.87</td><td>-1.05</td><td>7.38</td><td>7.84</td></tr><tr><td>Halo 383: GM later</td><td>5.8 × 10<sup>9</sup></td><td>3.3 × 10<sup>6</sup></td><td>3.9 × 10<sup>6</sup></td><td>-10.85</td><td>-11.53</td><td>0.06</td><td>0.78</td><td>5.9 × 10<sup>2</sup></td><td>1.8 × 10<sup>5</sup></td><td>-2.07</td><td>-1.26</td><td>7.78</td><td>7.38</td></tr><tr><td>Halo 383: GM higher<br>mass (z = 0.9)</td><td>7.5 × 10<sup>9</sup></td><td>1.1 × 10<sup>7</sup></td><td>7.4 × 10<sup>7</sup></td><td>-12.45</td><td>-14.00</td><td>0.81</td><td>0.55</td><td>7.8 × 10<sup>6</sup></td><td>3.5 × 10<sup>7</sup></td><td>-1.68</td><td>-1.60</td><td>7.57</td><td>7.12</td></tr><tr><td>Halo 339</td><td>5.2 × 10<sup>9</sup></td><td>-</td><td>5.8 × 10<sup>6</sup></td><td>-</td><td>-11.10</td><td>-</td><td>0.56</td><td>-</td><td>3.5 × 10<sup>6</sup></td><td>-</td><td>-1.05</td><td>-</td><td>7.83</td></tr><tr><td>Halo 261: GM earlier</td><td>6.9 × 10<sup>9</sup></td><td>-</td><td>1.9 × 10<sup>7</sup></td><td>-</td><td>-12.40</td><td>-</td><td>0.37</td><td>-</td><td>5.1 × 10<sup>6</sup></td><td>-</td><td>-1.02</td><td>-</td><td>7.90</td></tr><tr><td>Halo 261</td><td>6.6 × 10<sup>9</sup></td><td>-</td><td>1.5 × 10<sup>7</sup></td><td>-</td><td>-12.20</td><td>-</td><td>0.51</td><td>-</td><td>7.5 × 10<sup>6</sup></td><td>-</td><td>-1.06</td><td>-</td><td>7.82</td></tr><tr><td>Halo 261: GM later</td><td>6.6 × 10<sup>9</sup></td><td>-</td><td>9.7 × 10<sup>6</sup></td><td>-</td><td>-12.65</td><td>-</td><td>0.46</td><td>-</td><td>1.5 × 10<sup>7</sup></td><td>-</td><td>-1.15</td><td>-</td><td>7.59</td></tr><tr><td>Halo 153 (z = 0.6)</td><td>1.1 × 10<sup>10</sup></td><td>-</td><td>4.2 × 10<sup>7</sup></td><td>-</td><td>-14.69</td><td>-</td><td>0.56</td><td>-</td><td>4.3 × 10<sup>7</sup></td><td>-</td><td>-1.02</td><td>-</td><td>7.89</td></tr></tbody></table>
# APPENDIX B: RESOLVED FEEDBACK IN EDGE2 SIMULATIONS

The strength of the EDGE approach lies in its uncompromisingly high resolution to directly resolve key stellar feedback processes in the ISM of our dwarf galaxies. We quantify this statement in this appendix.

## B1 SN feedback

Resolving the cooling radius of SNe remnants is a key modelling milestone to accurately capture the momentum-build up during the Sedov–Taylor phase (e.g. Kim & Ostriker 2015). Our injection scheme computes the cooling radius as

$$R_{\text{cool}} = 30.0 \text{ pc} \left( \frac{n_{\text{H}}}{\text{cm}^{-3}} \right)^{-3/7} \left( \frac{E}{10^{51} \text{ erg}} \right)^{2/7} \left( \frac{Z}{Z_{\odot}} + 0.01 \right)^{1/7}, \quad (\text{B1})$$

where $n_H$ is the hydrogen density in the gas cell, $E$ is the total energy injected (this can be $\ge 10^{51}$ erg if more than one SN are exploding in the same timestep), and $Z$ is the gas metallicity (see Hopkins et al. 2018, appendix D for a derivation of this equation and e.g. Cioffi, McKee & Bertschinger 1988; Thornton et al. 1998; Kim & Ostriker 2015 for similar scalings).

Fig. B1 shows the total number of SNe events in each simulated dwarf galaxy over their whole evolution (bar pair), divided between CCSNe (purple) and SNeIa (orange). The filling of the bar indicates ‘resolved’ events where $R_{\text{cool}} \ge 6 \Delta x$ for which we directly inject thermal energy rather than momentum.

All dwarf galaxies have more than $\ge 87$ per cent of CCSNe events resolved, even as galaxy masses and ISM pressures increase (towards the right). A larger fraction of SNeIa are unresolved, most likely because these events can occur in lowerdensity environments, long after star formation when the adaptive resolution has been degraded. None the less, since CCSNe provide the overwhelming majority of the SN feedback budget, such statistics represent a significant achievement in the robustness of SN feedback modelling.

## B2 Radiative feedback

Even if SN feedback is accurately captured as we have shown above, the coupling of this feedback with the surrounding ISM and the subsequent ability to drive galactic outflows depend on the gas conditions in which the explosions occur. These conditions in turn depend on other, pre-SN feedback channels, especially radiative feedback.

For radiative feedback, the key length-scale to resolve is the Strömgren radius associated to the H II region around the star formation event. Requiring equilibrium between ionization and recombination of a hydrogen sphere, we obtain the Strömgren radius (Strömgren 1939) as

$$R_S = \left( \frac{3 Q}{4 \pi n_H^2 \alpha_B} \right)^{1/3}. \quad (\text{B2})$$

Here, $Q$ is the production rate of hydrogen-ionizing (Lyman-Compton) photons and $\alpha_B$ is the case-B recombination rate of hydrogen. Note that this equation assumes spherical symmetry, a homogeneous gas distribution and neglects the respective contributions of helium and metallic ions to the ionization and recombination balance. None the less, this is enough to gain order-of-magnitude estimates of whether ionization fronts from photoionization feedback are resolved, since gas density is the strongest driver in equation (B2).

<figure><img src="image_14.png" /><figcaption>**Figure B1.** Fractions of resolved CCSNe and SNeIa explosions over the full cosmological history of each galaxy. Bars show the total number of SNe, with filling highlighting those for which we resolve the cooling radius (see the text). Each bar pair correspond to an individual dwarf galaxy, ordered with growing halo mass as in Table A1. Irrespective of the mass scale, $\ge 85$ per cent of CCSNe, which dominate the SN feedback budget, are resolved.</figcaption></figure>

We can start by estimating an order of magnitude for $R_S$ when our star formation algorithm spawns a stellar particle (see Section 2.4 for details). The density at which stellar particles are spawned peaks around $\approx 500 m_p \text{ cm}^{-3}$. Accounting for the depletion of $300 M_{\odot}$ into a stellar particle, and assuming a primordial mixture, we input $n_H = 200 \text{ cm}^{-3}$ into equation (B2). Given our choice of SED (Bruzual & Charlot 2003), a zero-age $300 M_{\odot}$ stellar population has an ionizing rate of $Q \approx 10^{49} s^{-1}$ (see e.g. Rosdahl et al. 2018, appendix D for a visualization). Assuming the temperature-dependent $\alpha_B$ from Hui & Gnedin (1997) evaluated at $T = 2 \times 10^4$ K (see their appendix A), we obtain $R_S \approx 2$ pc. This is already comparable to the maximum resolution of the simulation, $\Delta x \approx 3$ pc. Pushing to $n_H = 10^3 \text{ cm}^{-3}$ which is common for star formation events in our simulations, we find $R_S \approx 0.7$ pc.

The choice of the SED (through $Q$), of the recombination rate ($\alpha_B$) or the assumptions of a primordial plasma will slightly modify these findings, but all these parameters enter equation (B2) with weaker scalings than density. As a result, the Strömgren sphere associated to the birth of a stellar particle is likely to be poorly resolved. This shortcoming is intrinsic to our star formation algorithm and cannot be avoided without modifying the very nature of the recipe (see e.g. Katz et al. 2024 for such a change).

None the less, even if the Strömgren sphere is unresolved at birth, the dynamics of the surrounding ISM could quickly alleviate this issue. Lowering the density to $n_H = 10 \text{ cm}^{-3}$, for example, gives $R_S \approx 14$ pc and a much better-resolved H II region. To test this, we identify young stellar particles with ages $\le 5$ Myr for each saved snapshot in the history of our individual galaxies. For each young stellar particle, we identify the density of their host gas cell and their ionizing output using our assumed SED at the particle’s stellar age and metallicity. We then use equation (B2) to compute $R_S$ (using the same $\alpha_B$ as previously) and consider a Strömgren sphere to be resolved if $R_S \ge 6 \Delta x$. In this case, a spherical H II region would be
<figure><img src="image_15.png" /></figure>

**Figure B2.** Fractions of resolved Strömgren spheres for each dwarf galaxy along the cosmological history of each dwarf, with filling highlighting those for which we resolve the Strömgren radius with six and four resolution elements (>200 and 60 elements per sphere, respectively). On average, ≈ 40 per cent of Strömgren spheres are resolved by 6Δx across the suite, with fractions as low as ≈ 20 per cent for the most massive dwarf galaxies (towards the right). Two dwarf galaxies are absent as their saved snapshots do not sample young stars and H II regions.

captured by > 200 resolution elements, which is enough to accurately capture its ionization and temperature structure (e.g. Deng et al. 2024a). We also compute the fraction of resolved Strömgren spheres with $R_S \ge 4 \Delta x$, which is a more marginal case with ≈ 60 resolution elements per sphere.

Fig. B2 shows the obtained statistics, with bars showing the resolved fraction of Strömgren spheres over the history of our dwarfs. Some galaxies are missing from the plot as their (short) formation histories and snapshot spacing does not allow sampling of < 5 Myr stars. Across the whole suite, we resolve very well 42 per cent of Strömgren spheres on average, with the smallest fraction being 23 per cent for ‘Halo 383: GM early’. Assuming a less stringent requirement of $R_S \ge 4 \Delta x$; this average fraction climb to 52 per cent with a minimum of 31 per cent. Similarly, shifting the age cut for ‘young’ stars from 5 to 10 Myr leads to a 57 per cent average resolved fraction with a minimum 40 per cent fraction.

Overall, unlike the situation for SNe for which the feedback modelling is nearly all resolved (Fig. B1), only a rough half of radiative feedback events are numerically well captured at our resolution. Furthermore, the fraction of resolved events falls rapidly as galaxies become more massive and their ISM denser. This is a clear area of improvement for radiation-hydrodynamics simulations aiming to explicitly model radiative feedback (see Section 5.2 for a discussion).

## APPENDIX C: BURSTINESS OF STAR FORMATION

Fig. C1 shows the distribution of SFRs for each dwarf galaxy over its full cosmological assembly history. Symbols show the median SFR, while the lines show the 16–84 percentiles of the distribution.

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

**Figure C1.** Median and 16–84 percentiles of the SFR over the full assembly history of each dwarf between z = 99 and 0. The *Hubble*-time-average SFRs are well converged between the EDGE1 and EDGE2 models. With radiative feedback, EDGE2 shows a less bursty star formation mode with a reduced scatter around the median.

We compute these statistics by computing the SFR averaged over 10 Myr for the *Hubble* time and eliminate vanishing SFRs. Only dwarf galaxies which have at least 10 snapshots with non-vanishing SFRs are shown on Fig. C1. (Some of our low-mass dwarf galaxies can have very short SFHs and billions of years of vanishing SFRs, leading to heavily biased median and 16–84 estimates.)

The median SFRs are well converged between the EDGE1 and EDGE2 models (blue points and red diamonds, respectively) in line with the aligned $M_*$ reported in Section 3.1. Computing the average extent of the 16–84 interval normalized by its median across each suite, EDGE2 shows a reduction in scatter of ≈ 2 compared to EDGE1. Such reduced burstiness for star formation is expected as a by-product of radiative feedback and aligns with the much calmer regulation mode in EDGE2 (Section 3).

## APPENDIX D: STAR FORMATION AND SUPERNOVA STATISTICS

To visualize the consequences of the much denser ISM of the EDGE2 model highlighted in Section 3, Fig. D1 shows the densities at which stellar particles are formed (left), and at which the CCSNe from a stellar particle explode (i.e. SNe younger than 8 Myr which corresponds to the main-sequence lifetime of a ≈ 20 $M_\odot$ massive star). We show the histograms of all SF and SNe events across four EDGE2 galaxies that span the full range of $M_*$ (We verified that conclusions are unchanged if using other galaxies.) Unfortunately, star formation and SNe statistics were not recorded in the EDGE1 model, so a direct comparison is not possible.

As expected, no star formation event (left) occurs below our density threshold for SF (≥ 300 $m_\odot$ cm$^{-3}$). For low-$M_*$ objects (black and brown), there is a clear peak around ≈ $10^3$ $m_\odot$ cm$^{-3}$ that quickly tapers off. As dwarf's galaxies increase in $M_*$ (green and purple),
<figure><img src="image_17.png" /><figcaption>**Figure D1.** Density distribution of star formation (left) and SNe (right) events across the history of four EDGE2 dwarf galaxies. As objects become more massive, star formation and SN events at increasingly high densities start to appear.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure D2.** Same as Fig. D1, but for the four resimulations of ‘Halo 605’ that vary input parameters as described in Section 5. For the same galaxy, forbidding star formation in warm gas ($T_{\star} = 10^2$ K) recreates a tail of high-density star formation events by reducing the clustering of star formation and subsequent SN feedback.</figcaption></figure>

however, the tail of star-forming events extends to higher densities reaching ($\ge 10^6 m_p \text{ cm}^{-3}$).

The direct consequence is that SNe explosions occur in increasingly dense gas (right-hand panel). The amount of momentum emerg-

ing from the Sedov–Taylor phase of an SN remnant only weakly depends on the background density (e.g. Kim & Ostriker 2015). But since our spatial scales are ultimately limited by numerical resolution, explosions in an increasing dense medium couple this momentum to
an ever-increasing mass, thus yielding smaller velocities and weaker outflows.<fnref n="3" />

Part of these effects are driven by our algorithms, which we illustrate now. In particular, our star formation algorithm requires gas to be both sufficiently dense and sufficiently cold to trigger star formation (Section 2.4). Fig. D2 shows the same statistics as Fig. D1 for the runs described in Section 5, in which the temperature criteria are raised and lowered. In particular, for a fixed object, a lower $T_*$ drives an extended tail of star formation events and SNe at high densities (gold), while a higher $T_*$ truncates both to lower densities.

These results should be attributed to interactions between radiative feedback and our star formation prescription. Following a star formation event, the dense gas surrounding the newborn stellar particle is immediately heated to $T \ge T_* = 10^3$ K by radiative feedback. Any gas in the immediate vicinity of a star-forming event is thus prevented from forming new stars, staying warm and dense until the first SNe can clear it ($\approx 5$ Myr later). This suppresses the clustering of star formation and SNe, making it harder to generate powerful outflows.

Even though we isolated the $T_*$ parameter here, we stress that this parameter is not the only driver of this behaviour. Rather, the combination of $T_*$, the heating power of stellar radiation, and the surrounding ISM densities driving gas cooling and collapse all contribute towards setting the distribution of star-forming densities. For example, a harder SED will raise more gas over a fixed $T_*$ and drive a slightly longer tail of star formation at higher densities (Fig. D2, orange versus red). And, at fixed SED and $T_*$, a more massive galaxy with a deeper gravitational potential well and higher ISM densities will retain high densities more effectively and compound this effect (Fig. D1, purple versus brown).

Combined with under-resolved radiative feedback events (Appendix B), these findings highlight some key uncertainties remaining in our radiation-hydrodynamics simulations. None the less, while these effects drive large changes on the statistics of star formation and SNe events, they have a much more moderate impact on the integrated $M_*$ and other observables (Sections 4 and 5). Observables such as rest-frame optical emission lines are likely more sensitive to changes in the structure of the ISM and its star-forming regions. In future iterations of the EDGE model, we will continue to improve our implementations and compare with these vital observational clues.

This paper has been typeset from a T\_EX/L\_T\_EX file prepared by the author.

<footnote n="3">³For a typical Sedov–Taylor momentum ($\rho_{SN} = 5 \times 10^5$ M$_\odot$ km s$^{-1}$) exploding in a 3 pc-cell with $\rho = 10^5$ m$_p$ cm$^{-3}$, the generated velocity is 7 km s$^{-1}$. This is insufficient to escape the gravitational potential well of even our small dwarfs ($v_{circ} \ge 10$ km s$^{-1}$).</footnote>