# Megatron: disentangling physical processes and observational bias in the multi-phase ism of high-redshift galaxies

Nicholas Choustikov, Harley Katz, Alex Cameron, Aayush Saxena, Julien

Devriendt, Adrianne Slyz, Martin Rey, Corentin Cadiou, Jeremy Blaizot, Taysun Kimm, et al.

► To cite this version:

Nicholas Choustikov, Harley Katz, Alex Cameron, Aayush Saxena, Julien Devriendt, et al. Megatron: disentangling physical processes and observational bias in the multi-phase ism of high-redshift galaxies. 2025. hal-05352217

HAL Id: hal-05352217

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

Preprint submitted on 6 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MEGATRON: DISENTANGLING PHYSICAL PROCESSES AND OBSERVATIONAL BIAS IN THE MULTI-PHASE ISM OF HIGH-REDSHIFT GALAXIES

NICHOLAS CHOUSTIKOV¹*, HARLEY KATZ²,³, ALEX J. CAMERON¹, AAYUSH SAXENA¹,⁴,

JULIEN DEVRIENDT¹, ADRIANNE SLYZ¹, MARTIN P. REY⁵, CORENTIN CADIOU⁶, JEREMY BLAIZOT⁷, TAYSUN KIMM⁸,
ISAAC LASETER⁹, KOSEI MATSUMOTO¹⁰, AND JOKI ROSDAHL⁷

¹Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, United Kingdom

²Department of Astronomy & Astrophysics, University of Chicago, 5640 S Ellis Avenue, Chicago, IL 60637, USA

³Kavli Institute for Cosmological Physics, University of Chicago, Chicago IL 60637, USA

⁴Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, United Kingdom

⁵University of Bath, Department of Physics, Claverton Down, Bath, BA2 7AY, UK

⁶Institut d'Astrophysique de Paris, Sorbonne Universites, CNRS, UMR 7095, 98 bis bd Arago, 75014 Paris, France

⁷Université Claude Bernard Lyon 1, CRAL UMR5574, ENS de Lyon, CNRS, Villeurbanne, F-69622, France

⁸Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea

⁹Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA and

¹⁰Sterrenkundig Observatorium Department of Physics and Astronomy Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium

Version October 9, 2025

ABSTRACT

Now detected out to redshifts of *z* ~ 14.5, the rest-frame ultraviolet and optical spectra of galaxies encode numerous physical properties of the interstellar medium (ISM). Accurately extracting these properties from spectra remains a key challenge that numerical simulations are uniquely suited to address. We present a study of the observed ISM of galaxies in MEGATRON: a suite of cosmological radiation hydrodynamics simulations coupled to on-the-fly non-equilibrium thermochemistry, with multiple prescriptions for star formation/feedback and parsec-scale resolution; capable of directly predicting spectroscopic properties of early galaxies. We find that irrespective of feedback physics used, the ISM of high-redshift galaxies is denser, less metal enriched, and subject to higher ionization parameters and radiation fields compared to similar mass galaxies in the local Universe – in agreement with interpretations of JWST observations. Using common observational techniques to infer bulk galaxy properties, we find that ISM gas density controls the slope of the mass-metallicity relation. Similarly, at the densities reached in some high-redshift galaxies, O32 becomes a density tracer rather than one of ionization parameter. This motivates the use of other line ratios like C43 and N43 to infer the ionization state of the gas. Finally, various feedback models populate different regions of strong-line diagnostic diagrams as the line ratios are sensitive to the feedback-modulated density-temperature structure of the ISM. Therefore, observed strong-line diagnostics can provide a strong constraint on the underlying physics of star formation and feedback in the high-redshift Universe.

Subject headings: high-redshift galaxies, ISM, galaxy formation

1. INTRODUCTION

The launch of JWST (Gardner et al. 2023), has revo-
lutionized our understanding of the interstellar medium
(ISM) at high-redshift. In particular, the Near-Infrared
Spectrograph (NIRSpec, Böker et al. 2023) provides the
capability to directly measure the spectra of galaxies in
the Epoch of Reionization in the rest-frame UV-optical.

Bright emission lines in the rest-frame UV and opti-
cal can be used to understand the properties of the ISM,
such as the electron density, temperature, pressure, ion-
ization parameter, and chemical abundances (see Kew-
ley, Nicholls & Sutherland 2019, for a review). Likewise,
the shape of the continuum emission contains informa-
tion about the age and properties of the stellar popula-
tions or accreting black holes (e.g. Iyer et al. 2025), the
presence of dense neutral gas reservoirs around galaxies
(e.g. Heintz et al. 2025), and electron temperatures in the
ISM (e.g. Laseter et al. 2024; Morishita et al. 2024; Hayes
et al. 2025; Cataldi et al. 2025; Pollock et al. 2025). The
spectra can be used to infer the possible presence of ex-
otic stellar populations (e.g. Cameron et al. 2024; Cullen

et al. 2025), reddening due to dust attenuation or neb-
ular emission (e.g. Saxena et al. 2024b; Topping et al.
2024a; Cullen et al. 2024), or to discover the existence of
entirely new classes of objects, such as ‘Little Red Dots’
(e.g. Matthee et al. 2024).

Spectroscopic surveys of high-redshift galaxies (e.g. Cameron et al. 2023b; Sanders et al. 2023; Hu et al. 2024; Roberts-Borsani et al. 2024, 2025; Hayes et al. 2025; Tang et al. 2025a; Shapley et al. 2025; Pollock et al. 2025) have identified four key characteristics that differentiate the ISM at high redshift from that in the local Universe:

**Harder ionizing spectra:** A large number of high ionization lines (e.g. He II λ1640, N IV] λλ1483, 1486 and C IV λλ1548, 1550) have now been observed at high-redshift (e.g. Bunker et al. 2023; Castellano et al. 2024; Topping et al. 2024b, 2025b; Naidu et al. 2025). While these lines are often attributed to AGN (e.g. Feltre, Charlot & Gutkin 2016; Laporte et al. 2017), this emission can also be produced by low-metallicity massive stars (e.g. Stark et al. 2015). In either case, the detection of these lines in high-redshift galaxies implies the presence of sources with hard ionizing photons (Berg et al. 2021).

**Lower gas-phase metallicities:** Measurements of

<footnote>*E-mail: nicholas.choustikov@physics.ox.ac.uk</footnote>
gas-phase oxygen abundance have found that galaxies of a given mass at high-redshift tend to be systematically less enriched, resulting in redshift-dependent normalizations and slopes in the mass-metallicity relation (MZR, e.g. Heintz et al. 2023; Nakajima et al. 2023; Curti et al. 2023, 2024; Morishita et al. 2024; Chemerynska et al. 2024a; Sarkar et al. 2025; Pollock et al. 2025). Furthermore, these works have found evidence for departures from the low-redshift ‘Fundamental Metallicity Relation’, indicating that galaxies of a given stellar mass and star formation rate are also less enriched at high-redshift.

**Higher densities:** Numerous galaxies with direct determinations of electron densities have been observed (e.g. Isobe et al. 2023a; Reddy et al. 2023; Abdurro'uf et al. 2024; Li et al. 2025a; Topping et al. 2025a). In most cases, densities are found to be higher in the early Universe. Furthermore, some objects have now been found with extreme densities $n_e > 10^5$ cm$^{-3}$ (Senchyna et al. 2024; Topping et al. 2024b), rarely seen in the local Universe except in extreme star-forming galaxies (Mingozzi et al. 2022). Finally, recent work combining rest-frame optical and IR emission have necessitated multiple H II regions of different density to explain the observed emission (Harikane et al. 2025; Usui et al. 2025). Altogether, this implies that the ISM of high-redshift galaxies is typically highly complex and inhomogeneous.

**Higher ionization parameters:** Many high-redshift galaxies are found to have line ratios consistent with large ionization parameters (e.g. Cameron et al. 2023b; Mascia et al. 2023; Nakajima et al. 2023; Topping et al. 2025b; Tang et al. 2025b; Glazer et al. 2025; Hayes et al. 2025), defined as the number of ionizing photons per hydrogen atom and typically traced by O32 ([O III]λ5007/[O II]λ3726, 3728). In particular, stacks have been used to demonstrate a correlation between redshift and O32 (Sanders et al. 2023; Roberts-Borsani et al. 2024), showing that these conditions become more extreme in the early Universe. Spectroscopic measurements of high-redshift galaxies reveal high emission line equivalent widths (Boyett et al. 2024; Endsley et al. 2024) and high ionizing photon production efficiencies (Atek et al. 2024; Saxena et al. 2024a; Simmonds et al. 2024a,b; Laseter et al. 2025; Choustikov et al. 2025). Together, these imply that the ISM of galaxies at high-redshift had significantly higher ionization parameters.

Finally, JWST observations reveal an over-abundance of bright galaxies in the early Universe (e.g. Finkelstein et al. 2023; Donnan et al. 2023; Leung et al. 2023; Harikane et al. 2024; Chemerynska et al. 2024b) compared to what was predicted by most pre-JWST galaxy formation models (see review by Robertson 2022). Possible solutions include increasing the star-formation efficiency to produce more stars (Dekel et al. 2023, c.f. Ferrara, Manzoni & Ntormousi 2025), allowing burstier star formation histories to increase scatter in the $L_{UV}-M_{halo}$ relation (Mason, Trenti & Treu 2023; Shen et al. 2023; Sun et al. 2023; Kravtsov & Belokurov 2024), or modifying the SFR-$L_{UV}$ relation by including the nebular continuum (Katz et al. 2025), varying the stellar initial mass function (e.g. Yung et al. 2024; Cueto et al. 2024; Trinca et al. 2024), modifying the physical properties of stars (Katz et al. 2024a), and including contributions from AGN (e.g. Hegde, Wyatt & Furlanetto 2024).

While all of these models have been explored in terms of their effects on the UV luminosity function, their impact on the ISM remains less understood.

The most robust approach for inferring physical properties from spectra relies on using theoretical line ratios, set by quantum mechanics that are sensitive to quantities such as electron temperature or density. When this is not possible (e.g. ionization parameter), one typically relies on equilibrium photoionization models that make simplistic assumptions on star formation history and geometry (see Kewley, Nicholls & Sutherland 2019, for a review). One-dimensional models run with codes such as CLOUDY (Ferland et al. 2017; Chatzikos et al. 2023) and MAPPINGS (Dopita & Sutherland 1996; Sutherland et al. 2018) are either used in isolation (e.g. Feltre, Charlot & Gutkin 2016; Gutkin, Charlot & Bruzual 2016), assembled from weighted combinations of models (e.g. Lebouteiller & Ramambason 2022; Marconi et al. 2024) or applied as a post-processing step to cosmological simulations (e.g. Hirschmann et al. 2017; Vijayan et al. 2021, 2025; Hirschmann et al. 2023; Katz et al. 2019a,b, 2022c,a, 2023a; Wilkins et al. 2023; Choustikov et al. 2024a; Nyhagen et al. 2024; Lovell et al. 2024, 2025; Giovinazzo et al. 2025).

While this approach is very computationally efficient, it begins to break down under certain conditions. For example, the presence of turbulent gas in H II regions can have non-negligible effects on emission line ratios, even when compared to one-dimensional models with the same average conditions (Gray & Scannapieco 2017; Jin, Kewley & Sutherland 2022). Secondly, non-equilibrium effects can be important in the presence of non-thermal physics such as cosmic rays (Katz et al. 2022b), or in regions where cooling timescales are sufficiently short that gas can be over-ionized (Oppenheimer & Schaye 2013; Katz et al. 2022b; Richings et al. 2022; Ploeckinger et al. 2025). Finally, galaxy spectra integrated over multiple unresolved H II regions can also lead to unrepresentative emission line ratios (Cameron, Katz & Rey 2023). As a result, we are motivated to explore the use of cosmological simulations with on-the-fly non-equilibrium thermochemistry that can complement photoionization models by overcoming some of these aforementioned issues.

In this work, we aim to understand the physics underpinning observed properties of the high-redshift ISM, in particular focusing on the physical and observational effects which bias our ability to extract information from spectra. To accomplish this, we use the MEGATRON simulations (Katz et al. 2025; Rey et al. 2025; Cadiou et al. 2025). These simulations employ a detailed thermochemistry network of primordial species, metals, and molecules, coupled to on-the-fly multi-frequency radiation transport and a state-of-the-art galaxy formation model, which allows us to self-consistently predict the spectra of the simulated high-redshift galaxies. In particular, we use this unique capability to conduct a direct comparison with JWST-observed high-redshift galaxies, most notably measuring all bulk galaxy properties from the mock spectra in the same way that is typically done for JWST observations.

This paper is organized as follows. In Section 2, we describe the MEGATRON simulations and outline our spectroscopic analysis approach. In Section 3, we investigate each of the four unique ISM properties of high-redshift
galaxies. Finally, in Section 4, we discuss caveats of our approach before concluding in Section 5.

## 2. METHODS

In this section, we provide an overview of our methodology: describing the variety of simulations and observational techniques used in the present work.

### 2.1. Galaxy Evolution with Non-Equilibrium Thermochemistry

For the purpose of this work, we use data from the high-redshift suite of the MEGATRON simulations (Katz et al. 2025). These are four separate high-resolution cosmological radiation hydrodynamical simulations run with RAMSES-RTZ<fnref n="1" /> (Katz 2022) and various models (described in Section 2.1.1) for star formation and feedback. In each simulation, we follow the evolution of an early-forming Milky Way mass halo in a relative over-density to a final redshift of $z = 8.5$. The simulations use a constant-comoving resolution scheme, which refines cells based on a quasi-Lagrangian approach as well as requiring that the Jeans length is resolved by at least 4 cells. As a result, at $z \sim 8.5$, the maximum physical resolution is $\sim 5$ pc $h^{-1}$, compared to $\sim 1$ pc $h^{-1}$ when the first stars form at $z \sim 30$. For full numerical details, we direct the interested reader to Katz et al. (2024b); Katz et al. (2025). Below we summarize the numerical methods relevant to the ISM.

Radiation and hydrodynamics are coupled on-the-fly to a non-equilibrium chemistry solver (Katz 2022) for primordial species, metals, and molecules; following the evolution of H I–II, He I–III, $e^-$, C I–VI, N I–VII, O I–VIII, Ne I–X, Mg I–X, Si I–XI, S I–XI, Fe I–XI, H$_2$, and CO. All metal ionization states not followed are assumed to be in collisional ionization equilibrium. Gas heating and cooling processes include photo-heating, photoelectric heating, H$_2$ formation heating, H$_2$ excitation/dissociation heating, compton heating/cooling, primordial cooling, H$_2$ cooling, CO cooling, dust recombination cooling, dust-gas collisional cooling (or heating), and metal line cooling (Katz et al. 2022b; Katz et al. 2024b, 2025). We assume that the dust-to-gas mass ratio varies with metallicity following Rémy-Ruyer et al. (2014) using the BARE-GRS dust composition model of Zubko, Dwek & Arendt (2004).

Star-formation is modeled based on the local, turbulent properties of the gas (Padoan & Nordlund 2011; Federrath & Klessen 2012; Kimm et al. 2017; Rosdahl et al. 2018). Specifically, stars can form in a gas cell if the following criteria are fulfilled: it has a density $n_H > 10 \text{ cm}^{-3}$, is a local density maximum, if the fluid flow is locally convergent, and if the turbulent Jeans length is smaller than the cell width. If these criteria are fulfilled, star-formation proceeds using a Schmidt law (Schmidt 1959) with an efficiency per free-fall time that depends on the turbulent properties of the gas (Padoan & Nordlund 2011; Federrath & Klessen 2012). For Population II star particles, we adopt a Kroupa IMF (Kroupa 2001) with a maximum mass of 300 $M_{\odot}$. Radiation from

Population II stars<fnref n="2" /> is taken from BPASS v2.2.1 (Eldridge et al. 2017; Stanway & Eldridge 2018). Stellar winds, core-collapse SN, and type Ia SN are modeled following Agertz et al. (2021) with minor modifications (Katz, et al., 2025).

#### 2.1.1. Variations to Sub-Grid Physics

The four high-redshift MEGATRON simulations differ in their subgrid physics in ways that loosely follow proposed solutions to the over-abundance of high-redshift bright galaxies. Here we summarize the key differences between the four models.

* **Efficient Star Formation:** This model represents the fiducial simulation, as described above. As this simulation has the weakest feedback compared to the others, it has the highest conversion efficiency of gas into stars.

* **Bursty Star Formation:** In this simulation, the energy injected by all SNe has been increased by a factor of 5 compared to the efficient star formation model. As a result, stellar masses are regulated and star-formation histories become significantly more stochastic and bursty.

* **Variable IMF:** In this simulation, the high-mass slope of the IMF is allowed to vary according to the local gas density and metallicity according to Marks et al. (2012). Furthermore, a metallicity-dependent fraction of high-mass stars (Kobayashi et al. 2006) explode as hypernovae with a mass-dependent energy (Nomoto et al. 2006). While the IMF is allowed to vary, we cap the maximum mass at 120 $M_{\odot}$ and use SEDs from STARBURST99 (Leitherer et al. 1999). As a result, stellar populations can form with very low mass-to-light ratios.

* **High $\epsilon_{\text{ff}}$, HN:** In this simulation, stars are formed in lower-density gas with an efficiency per free-fall time of 100%. Like the variable IMF model, a fraction of high-mass stars are able to explode as hypernovae (Kobayashi et al. 2006), with variable energy (Nomoto et al. 2006). As a result, this model regulates stellar masses significantly better and produces a less-dense ISM with bursty star-formation histories.

#### 2.1.2. Self-Consistent Synthetic Observations

The key benefit of evolving the full non-equilibrium thermochemistry of the ISM is that we are able to predict the intrinsic spectra of galaxies, accounting for the nebular continuum, nebular line emission, and the stellar continuum. Each mock spectrum considers all gas cells within 25% of the virial radius of a particular halo. Line emissivities and the nebular continuum are computed using PyNeb (Luridiana, Morisset & Shaw 2015) with atomic data for collisional emission from Version 10 of the CHIANTI database (Dere et al. 1997, 2019; Del Zanna et al. 2021). Finally, the stellar continuum is computed for each star particle based on the SEDs used in the simulation.

<footnote n="2">² Population III stars are also modeled explicitly (Storck et al. 2025). However, all of the galaxies considered in this work are already completely dominated by Population II stars.</footnote>

<footnote n="1">¹ A fork of RAMSES (Teyssier 2002) and RAMSES-RT (Rosdahl et al. 2013).</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1.— Projected UV-optical emission line maps normalized by the total Hβ flux for the galaxy with the largest specific star-formation rate in the bursty star formation simulation. Low surface-brightness gas has been brightened. On the bottom, we show the total intrinsic spectrum (in arbitrary log $f_ν$ units), highlighting the relative strengths of each line. This galaxy has a stellar mass of $10^{7.9} M_\odot$ at z = 8.93.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2.— Temperature-density phase diagrams for an example ~$10^{8.5} M_{\odot}$ stellar mass galaxy at $z = 8.6$ in the efficient star formation run, with all cells colored by their contribution to the total luminosity of a given emission line. We include UV, optical, and FIR lines in blue, green, and red respectively. No two emission lines are produced by the same gas, introducing potential biases when emission lines are combined to make measurements of ISM properties.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3.— Distributions of common UV, optical, and FIR emission line ratios for each MEGATRON simulation, highlighting the impacts that different subgrid prescriptions for star-formation/feedback have on the ISM, and therefore on its spectroscopic properties.</figcaption></figure>

Like all numerical simulations, we are limited by finite numerical resolution. At high densities, the Strömgren spheres of star particles are sometimes unresolved (see discussion in Choustikov et al. 2024a). To correct for this, we replace the emission of cells that contain unresolved Strömgren spheres with CLOUDY models<fnref n="3" /> (Ferland et al. 2017) using the local gas conditions and star particle properties (Katz et al. 2025). Finally, for the purpose of this work we ignore the effect of dust attenuation, scattering, and emission on the synthetic spectra of MEGATRON galaxies (see discussion in Section 4).

Figure 1 shows example emission line maps of a highly star-forming galaxy. A 1D spectrum for the same galaxy is also shown. While each emission line primarily orig-

inates from star-forming regions, they are not entirely co-spatial. For example, the Balmer series and lines like [O III] λ5007 can be seen to also trace the large-scale gas distribution around the galaxy, while high-ionization lines like C IV λλ1549 and He II λ1640 and auroral lines like [N II] λ5755 and [O II] λλ7325 primarily trace the H II regions.

Figure 2 shows density-temperature phase diagrams for a typical galaxy, weighted by emission line luminosity. UV lines are labeled in blue, optical lines are in green and FIR lines are given in red. No two emission lines are produced by the same gas in density-temperature space. For example, [O III] λ4363 is produced in preferentially hotter gas than [O III] λ5007. As a result, electron temperature measurements from these lines are typically biased high, skewing the inferred metallicity low (Cameron, Katz & Rey 2023, see also discussion in Kewley, Nicholls

<footnote n="3">³ RAMSES-RTZ has been shown to agree with CLOUDY to within 10% for idealized Strömgren sphere models (Katz et al. 2024b).</footnote>
& Sutherland 2019).

A simple constant-density or isobaric photoionization model cannot sample the full distribution of ISM properties that we find in these simulations. To highlight differences between observed properties for each simulation, Figure 3 presents histograms of common UV, optical, and FIR emission line ratios for each run. We find that in the case of each ratio, every simulation produces a different distribution. In particular, density- and temperature-dependent ratios show the largest variance between each run. These differences are driven by changes to sub-grid star formation/stellar feedback physics, which directly modify the density-temperature structure of the ISM (see discussion in Katz et al. 2024b). While each model has its own limitations, variations to the sub-grid physics allow for an ISM parameter exploration that would be otherwise impossible with a single simulation. These variations demonstrate that it may be possible to use an ensemble of observed emission line ratios to constrain sub-grid physics. We leave this exploration to future work. Altogether, the complex nature of galaxies in the MEGATRON simulations make them a unique tool for understanding the physics that give rise to observed galaxy spectra.

### 2.1.3. LyC Escape Fractions

Lyman continuum (LyC) escape fractions ($f_{esc}$) are computed for each halo resolved by more than 1,000 dark matter particles that hosts at least one star particle. The $f_{esc}$ calculation is performed in post-processing with the Monte Carlo radiation transfer code RASCAS (Michel-Dansac et al. 2020). For each star particle, we propagate monochromatic photon packets at 900 Å until they reach either $0.75R_{vir}$<fnref n="4" />, or they are absorbed by neutral hydrogen or dust. The number of photon packets is calculated as $100 \times$ the number of star particles with a minimum of $10^5$ and maximum of $10^7$ photon packets per halo. The photon packets are randomly distributed to star particles based on a multinomial distribution, weighted by the LyC luminosity of each star particle. The dust model used in post-processing follows exactly that used when running the simulation — it is a combination of the BARE-GR-S composition of Zubko, Dwek & Arendt (2004) with the dust-to-gas mass ratio as a function of metallicity from Rémy-Ruyer et al. (2014). The escape fractions quoted in this work represent angle-averaged values, i.e. the quantity needed for reionization calculations, rather than line-of-sight values, which is what an observer would see. The latter will be considered in future work.

## 2.2. Observational Methods

As our goal is to make the best possible comparison to high-redshift observational data, a unique aspect of our work is that we measure galaxy properties from the mock spectra and photometry, rather than relying on the true value from the simulation. Here we outline the

methods that we use, inspired by typical observational approaches.

### 2.2.1. ISM Properties from Emission Line Diagnostics: PyNeb

We employ a combination of strong emission line ratio diagnostics, defined in Table 1 (see also Kewley, Nicholls & Sutherland 2019), to infer properties of the simulated galaxies and compare with observations. We focus on line ratios from the same element in order to limit uncertainties due to chemical yields that are innate to all numerical simulations (e.g. Buck et al. 2021).

More specifically, rather than using the true electron number density ($n_e$), electron temperature ($T_e$), gas phase metallicity, or ionization, we measure these quantities using emission line ratios. This allows us to make the most direct comparisons with observational inferences. For example $n_e$ is measured in different ionization zones using S2, O2, O3, and N4. Likewise, $T_e$ is estimated directly with auroral lines (e.g. [O III] λ4363) and a corresponding strong line (e.g. [O III] λ5007). In practice, we use the GetTemDen function of PyNeb (Luridiana, Morisset & Shaw 2015) with atomic data for collisional emission from Version 10 of the CHIANTI database (Dere et al. 1997; Del Zanna et al. 2021) to measure these quantities. As temperature and density measurements are rare for the same observed high-redshift galaxy, we follow a relatively common observational procedure of assuming a temperature when measuring a density and similarly assume a density when measuring temperature. We adopt $T_e \sim 1.5 \times 10^4$ K for density measurements and $n_e \sim 300 \text{ cm}^{-3}$ for temperature measurements.

### 2.2.2. Oxygen Abundances from the Direct Method

To compute gas-phase oxygen abundance, we assume that the vast majority of oxygen emission is coming from the low and intermediate ionization zones, giving $\frac{O}{H} = \frac{O^+}{H^+} + \frac{O^{++}}{H^+}$ as follows:

$$ \frac{\mathrm{O}}{\mathrm{H}} = \frac{\epsilon_{\mathrm{H}\beta}(T_{[\mathrm{O}\,\mathrm{II}]}, n_{\mathrm{e},[\mathrm{C}\,\mathrm{III}]})}{\epsilon_{[\mathrm{O}\,\mathrm{II}]}(T_{[\mathrm{O}\,\mathrm{II}]}, n_{\mathrm{e},[\mathrm{C}\,\mathrm{III}]})} \cdot \frac{L([\mathrm{O}\,\mathrm{II}]\lambda\lambda3726, 3728)}{L(\mathrm{H}\beta)} + \frac{\epsilon_{\mathrm{H}\beta}(T_{[\mathrm{O}\,\mathrm{III}]}, n_{\mathrm{e},[\mathrm{C}\,\mathrm{III}]})}{\epsilon_{[\mathrm{O}\,\mathrm{III}]}(T_{[\mathrm{O}\,\mathrm{III}]}, n_{\mathrm{e},[\mathrm{C}\,\mathrm{III}]})} \cdot \frac{L([\mathrm{O}\,\mathrm{III}]\lambda5007)}{L(\mathrm{H}\beta)}. \quad (1) $$

Ionization corrections for metallicity are expected to be small because the neutral states of H and O are highly coupled via charge-exchange reactions and stars typically do not emit significant quantities of photons hard enough to create a substantial amount of $O^{+++}$ (Berg et al. 2021).

Observationally, deriving $T_{[\mathrm{O}\,\mathrm{II}]}$ is difficult due to the fact that the $O^+$ auroral lines at 7320,7330Å drop out of NIRSpec coverage at $z \sim 6.2$ and are extremely weak for a highly ionized ISM. Therefore, while we do measure this temperature directly (see the discussion in Section 3.2.2), we choose to compute $T_{[\mathrm{O}\,\mathrm{II}]}$ from $T_{[\mathrm{O}\,\mathrm{III}]}$ based on an empirical relation. While several such relations have been presented (Izotov et al. 2006; Pilyugin et al. 2012; Cataldi et al. 2025), we use Equation 3 of Cameron, Katz & Rey (2023), which was derived from an isolated galaxy simulation also run with RAMSES-RTZ, and corrects for temperature inhomogeneities in the ISM. Nonetheless,

<footnote n="4">⁴ Note that we have adopted 75% of the virial radius as the escape radius rather than at $R_{vir}$ due to a data format incompatibility between RASCAS and MEGATRON outputs. We have tested this parameter for a few haloes and found our results to be well converged as most absorption occurs in the ISM and inner CGM.</footnote>
<table><caption>TABLE 1
A LIST OF ALL EMISSION LINE DIAGNOSTIC RATIOS USED IN THIS WORK. THROUGHOUT THIS WORK, ALL EMISSION LINE RATIOS ARE CONSIDERED IN LOG BASE 10.</caption><thead><tr><th>Diagnostic Name</th><th>Emission Line Ratio</th><th>Property Traced</th></tr></thead><tbody><tr><td>R3</td><td>[O III]λ5007/Hβ</td><td>Metallicity</td></tr><tr><td>R23</td><td>([O III]λλ4959, 5007 + [O II]λλ3726, 3728)/Hβ</td><td>Excitation</td></tr><tr><td>S2</td><td>[S II]λ6731/[S II]λ6719</td><td>Density</td></tr><tr><td>O2</td><td>[O II]λ3728/[O II]λ3726</td><td>Density</td></tr><tr><td>O2O2</td><td>[O II]λλ3731, 7320, 7330, 7331/[O II]λλ3726, 3728</td><td>Temperature</td></tr><tr><td>O32</td><td>[O III]λ5007/[O II]λλ3727</td><td>Ionization parameter</td></tr><tr><td>O3O3</td><td>[O III]λ4363/[O III]λ5007</td><td>Temperature</td></tr><tr><td>O3</td><td>[O III]λ88μ/[O III]λ52μ</td><td>Density</td></tr><tr><td>C3</td><td>C III]λ1908/C III]λ1906</td><td>Density</td></tr><tr><td>C43</td><td>C IV]λλ1550/C III]λλ1907</td><td>Ionization parameter</td></tr><tr><td>N4</td><td>N IV]λ1486/N IV]λ1483</td><td>Density</td></tr><tr><td>N43</td><td>N IV]λλ1483, 1486/N III]λ1750</td><td>Ionization parameter</td></tr></tbody></table>

this choice will not significantly impact our results because $O^{++}$ typically dominates the ionization state of oxygen in galaxies with direct [O III] λ4363 detections (Andrews & Martini 2013; Curti et al. 2017, 2020, 2023; Laseter et al. 2022).

Finally, when inferring metallicity, we use electron number densities computed from the $C^{++}$ doublet as these can be detected at high redshift (Topping et al. 2025a) and better trace the gas phases where [O III] λ5007 originates (see Figure 2). In contrast, we find that using $S^+$ densities often give us incorrect results, due to the fact that this density measurement is tracing a different ionization zone (see discussion in Berg et al. 2021). The efficacy of gas-phase metallicity measurements is presented in Appendix A.

### 2.2.3. SED Fitting with BAGPIPES

As our goal is to make a direct comparison with typical observational approaches, we derive bulk stellar population properties (i.e. total stellar mass) by performing SED fitting on the photometry of all MEGATRON galaxies with $M_{UV} \le -15$ with the public version of BAGPIPES (Carnall et al. 2018), using an approach similar to that outlined in Choustikov et al. (2025). To accomplish this, we first compute photometric magnitudes in numerous standard JWST filters (F115W, F150W, F200W, F277W, F335M, F356W, F410M, and F444W) by integrating over each intrinsic SED with SEDPY (Johnson 2019). Next, we follow the ‘best case’ scenario by fixing the object redshift to the known value and ignoring the presence of dust (see also Narayanan et al. 2024; Cochrane et al. 2025).

We use the 2016 stellar population synthesis models from the BCO3 templates<fnref n="5" /> (Bruzual & Charlot 2003; Chevallard & Charlot 2016) with a Kroupa (2001) IMF. Nebular emission (both continuum and line) from these stellar populations is modeled using the 2017 version of CLOUDY (Ferland et al. 2017). Next, we use a non-parametric star formation history (SFH), following the continuity prior approach presented in Leja et al. (2019) with time bins set at $t_{bins} = [0, 5, 10, 25, 50, 100, \min(250, t_z)]$ Myr where $t_z$ is the age of the Universe at redshift $z$. Furthermore, we allow the non-parametric model to better reconstruct a more

bursty SFH by using $\sigma = 1$ and $\nu = 2$ in the Student t-distribution priors (Tacchella et al. 2022). This SFH prescription is chosen as it was recently shown to reproduce the stellar masses of SPHINX²⁰ galaxies well (Cochrane et al. 2025). Our remaining model priors are as follows: we use uniform priors for mass and metallicity of $log_{10}(M_*/M_\odot) \in [0, 10]$ and $log_{10}(Z/Z_\odot) \in [-3, 1]$ and for simplicity take a fixed ionization parameter of $log_{10}U = -2$. This simplification will not severely impact our results, as this approach was also shown to recover the intrinsic ionizing luminosity of SPHINX²⁰ galaxies well (Choustikov et al. 2025).

We present the efficacy of stellar mass recovery for MEGATRON galaxies in Appendix B. In general, we find that below a redshift of $z \sim 10.5$, the masses of galaxies in the efficient SF, bursty SF, and high $\epsilon_{ff}$ are well recovered, though with a slight tendency to over-predict the masses of the lowest mass objects (Choe et al. 2025). Next, we find that the stellar masses of the variable IMF run are typically over-predicted by a factor of $\sim 2$, in accordance with the fact that a more standard IMF was used for the SED fitting, whereas variable IMF models typically exhibit lower mass-to-light ratios for young stellar populations. However, at redshifts above $z \sim 10.5$, the scatter in stellar mass increases significantly, because the Balmer break redshifts beyond the F444W filter. Here, stellar masses are typically incorrect by an order of magnitude, highlighting the need for MIRI data in studies attempting to apply SED fitting to very high-redshift galaxies (e.g. Leung et al. 2024). In all cases, we find that the residuals are correlated with the amount of burstiness in the recent SFH, with galaxies with higher SFR$_{10}$/SFR$_{10}$ having up-scattered stellar masses. While outshining can in principle be important (Narayanan et al. 2024; Harvey et al. 2025), these galaxies are not sufficiently massive for it to matter in this case.

### 2.2.4. Ionizing Photon Production Efficiency

To compute ionizing photon production efficiencies ($\xi_{ion}$), we follow a common approach and assume Case B recombination with $T_e = 10^4$ K, $n_e = 100 \text{ cm}^{-3}$, $f_{esc} = 0$ so that:

$$ \xi_{ion} / [\text{Hz/erg}] = \frac{7.3 \times 10^{11} L(H\alpha)}{L_{1500,\text{int}}}, \quad (2) $$

where $L(H\alpha)$ is the flux of the $H\alpha$ line and $L_{1500,\text{int}}$ is the intrinsic UV luminosity density (e.g. Maseda et al.

<footnote n="5">⁵ While using consistent SPS models such as BPASS would potentially provide more accurate results, we chose to use the default setup of BAGPIPES to mimic uncertainties inherent to inferring stellar masses from observed data.</footnote>
<figure><img src="image_5.png" /></figure>

FIG. 4.— Smoothed BPASS spectral energy distributions for a 1 Myr, $10^6$ M$_\odot$ stellar cluster for a variety of metallicities, normalized at 1500Å (Eldridge et al. 2017; Stanway & Eldridge 2018). This illustrates the fact that lower-metallicity stellar populations show characteristically harder and more dominant ionizing spectra. These differences are important for higher-ionization species (e.g. N IV, C IV, O III, N III, C III, etc.) which are often observed in the high-redshift Universe. In the top panel, we show the ionization energies for species of elements considered in this work.

2020). When comparing to observations at $z \gtrsim 7$ where $H\alpha$ drops out, we instead use $H\beta$ with the case B result that $L(H\alpha) = 2.863L(H\beta)$. This efficiency should be corrected for the escape fraction of LyC radiation. However, to best compare with observational inferences, we ignore this effect (see discussion in Saxena et al. 2024a).

## 3. RESULTS

In this section, we explore the physics that governs the ISM conditions in the MEGATRON simulations, exploring each of the four key characteristics of the high-redshift ISM. To best compare with observations, we consider only galaxies that have $M_{UV} < -15$.

### 3.1. Hard Ionizing Spectra

In the early Universe, the ISM is generally expected to be denser and less metal enriched than in the local Universe. These conditions limit the degree of fragmentation (Safranek-Shrader, Milosavljević & Bromm 2014). Such environments are particularly conducive to the formation of stellar binaries (Moe, Kratter & Badenes 2019) which have been shown to boost the ionizing production efficiency of a given stellar population (Stanway, Eldridge & Becker 2016; Stanway & Eldridge 2019; Götberg et al. 2019). Furthermore, at fixed angular momentum, metal poor stars are predicted to have smaller radii, rotate more rapidly (e.g. Yoon & Langer 2005; Ekström et al. 2012), and live longer (Sana et al. 2012), all of which sustain the production of ionizing photons (Topping & Shull 2015; Xiao, Stanway & Eldridge 2018). While the ionizing SEDs used in MEGATRON are an assumption of the model rather than a prediction, nearly all SED models predict that lower metallicity stars produce more ionizing photons. The shape of low-metallicity stellar SEDs is instructive to discuss as context for the rest of the paper.

In Figure 4, we show BPASS$^6$ SEDs normalized at 1500 Å for a $10^6$ M$_\odot$ stellar population with an age of 1 Myr for a variety of metallicities. While SED shape is broadly consistent in the rest-frame UV (up to metal absorption features), they are different in the Lyman Continuum (LyC: $\lambda \le 912$ Å). At these wavelengths, stellar populations with sub-solar metallicities (0.1Z$_\odot$) produce significantly more ionizing photons, especially at higher energies (i.e. harder photons).

The hardness of the radiation field is crucial for determining the size of high-ionization zones of the ISM (see discussion in Berg et al. 2021). This is demonstrated with the top panel of Figure 4, where we show the wavelength ranges occupied by each ionization species of the elements focused on in this work (H, He, C, N, O, S). Certain high ionization species (e.g. N IV/N III, C IV/C III, O III/O II) fall within the range of the BPASS SEDs most dependent on metallicity. Therefore, line ratios such as C43, N43, O32 can be particularly affected by the SEDs of high-redshift stellar populations (e.g. Bunker et al. 2023; Katz et al. 2023b; Topping et al. 2024b; Hayes et al. 2025; Pollock et al. 2025). While higher-ionization ratios do exist (e.g. N54), these require ionizing photons with wavelengths $\lambda \lesssim 200$ Å, whose contributions from stellar populations are considered minimal (c.f. Sander et al. 2025) and likely require an AGN component (Feltre, Charlot & Gutkin 2016).

### 3.2. Low Metallicities

Next, we explore the second key physical aspect of the high-redshift ISM: investigating the metallicity properties of MEGATRON galaxies.

#### 3.2.1. The Mass-Metallicity Relation

The mass-metallicity relation (MZR: Lequeux et al. 1979; Tremonti et al. 2004; Lee et al. 2006) and fundamental metallicity relation (FMR: Ellison et al. 2008; Mannucci et al. 2010; Yates, Kauffmann & Guo 2012) have strong diagnostic power over key galaxy formation physics such as inflows, outflows, and enrichment (see review by Maiolino & Mannucci 2019). With JWST, we can now infer the MZR and FMR at high redshifts directly from observations.

Initially, Heintz et al. (2023) and Nakajima et al. (2023) both found evidence for evolution in the normalization and slope of the low-mass MZR towards higher redshifts, producing an MZR that was offset from low-redshift galaxies (Curti et al. 2020). Since then, other studies have confirmed the apparent offset in the MZR normalization, though much debate still exists about the slope of the MZR (Heintz et al. 2023; Nakajima et al. 2023; Curti et al. 2024; Morishita et al. 2024; Chemerynska et al. 2024a; Sarkar et al. 2025). This is of particular interest, as simulations (e.g. Ma et al. 2016; Torrey et al. 2019; Langan, Ceverino & Finlator 2020; Dubois et al. 2021; Pallottini et al. 2022; Ucci et al. 2023; D'Silva et al. 2023; Marszewski et al. 2024) and chemical evolution models (e.g. Kobayashi & Taylor 2023) often struggle with the normalization and slope of the low-mass end of the MZR (e.g. Curti et al. 2024).

<footnote>$^6$ Similar trends hold for STARBURST99 (Leitherer et al. 1999), which is used for the model with a variable IMF.</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 5.— Mass-metallicity relations for each of the high-redshift MEGATRON simulations. Unless specified as intrinsic, coloured points and lines represent masses and metallicities inferred using observational techniques (i.e. the direct method for metallicities and BAGPIPES SED fitting for stellar masses). Left: Observed gas-phase metallicity as a function of observed stellar mass, coloured by the Hα-derived star formation rate for all galaxies in MEGATRON with $M_{UV} \le -15$. We show comparisons with observational data of galaxies with direct [O III] λ4363 measurements above redshift $z = 5$ (Curti et al. 2020, 2023; Arellano-Córdova et al. 2022, 2025; Cameron et al. 2023a; Morishita et al. 2024; Pollock et al. 2025), running means (colored lines), and contours showing the distribution of intrinsic values (i.e. mass-weighed metallicity and true stellar mass) for each simulation (pale blue). Observational measurements of these quantities tend to imprint a slight bias on the MZR, by flattening the slope of the relation. Right: We show comparisons between the running mean of each simulation with the same observational data. We also indicate the impact of key physical processes, highlighting intrinsic biases (e.g. different feedback models) in gray and observational biases (e.g. temperature/density fluctuations) in blue.</figcaption></figure>

We show the observed mass-metallicity relation for each of the MEGATRON simulations, using direct-method metallicities (see Section 2.2.2) and stellar masses from BAGPIPES fitting (see Section 2.2.3) in Figure 5. On the left of Figure 5, we present the full dataset of mock-observed MEGATRON galaxies with $M_{UV} \le -15$, showing direct-method metallicities as a function of BAGPIPES stellar mass, colored by Hα-derived star formation rates (Kennicutt 1998). We also provide contours representing the distribution of intrinsic measurements of MEGATRON galaxies (i.e. mass-weighted average metallicity and simulation-derived stellar masses), running means for each simulation (colored lines) and observational data for galaxies at $z \ge 5$ with direct-method metallicities as black diamonds (Curti et al. 2020, 2023; Arellano-Córdova et al. 2022, 2025; Cameron et al. 2023a; Morishita et al. 2024; Pollock et al. 2025). On the right, we show how running means from

each simulation compare to observational inferences. We also provide indications as to the effect of various intrinsic (e.g. feedback approach, gray) and observational (e.g. temperature fluctuations, blue) processes which can affect the observed MZR. We find that each of the simulations produce a different overall mass-metallicity relation, suggesting that stellar/feedback physics has a strong impact on the shape of the MZR (see discussion in Maiolino & Mannucci 2019). Nevertheless, all simulations are able to reproduce the full scatter at the low-mass end of the observed MZR.

Furthermore, we compare the different simulations to understand the impact of each physical prescription choice on the shape of the MZR. First, we note that the efficient SF, bursty SF, and variable IMF runs all show characteristically shallower MZRs than the high $\epsilon_{ff}$ simulation. This is likely due to the fact that compared to the other three models, the high $\epsilon_{ff}$ run preferentially
forms stars in lower-density gas. This means that stellar feedback in the high $\epsilon_{ff}$ is better coupled to the ISM, driving more powerful outflows that are capable of ejecting enriched gas further out of low-mass galaxies (and limiting subsequent star formation), thus producing a steeper MZR. Next, considering the bursty SF run, we find that a stronger feedback prescription for the same stellar populations tends to reduce the overall normalization of the MZR, by strengthening outflows and removing enriched gas – impacting both the mass and abundance of a galaxy at later times. In contrast, the variable IMF simulation produces more massive stars for a given stellar population. While these massive stars can increase the amount of feedback, this is not guaranteed for a very shallow IMF slope (see Appendix C) and the dominant effect is to reduce the final mass of the galaxy (measured after mass-loss due to winds/SNe) while also increasing the gas-phase metallicity of each galaxy since massive stars yield more metals. The presence of hypernovae can further exacerbate the impact of massive stars, by injecting comparatively more oxygen than supernovae (e.g. Nomoto et al. 2006), while also increasing the strength of feedback and thus reducing the final stellar masses of galaxies.

Finally, the turbulent structure of the ISM will make H II region properties vary across a galaxy, impacting the measured oxygen abundance (e.g. Cameron, Katz & Rey 2023; Harikane et al. 2025; Usui et al. 2025). Specifically, temperature fluctuations tend to bias the [O III] $\lambda$4363-derived electron temperature high (see Figure 2), thus driving the measured metallicity low (Kewley, Nicholls & Sutherland 2019). Similarly, density fluctuations tend to bias density measurements towards the lower-density gas$^7$, therefore skewing metallicity measurements low.

Comparing the observed scatter and intrinsic distribution of data from each MEGATRON simulation, we find that the observational approaches do not introduce large biases into the MZR when averaged across many galaxies. However, we find that the observed MZR tends to be slightly shallower in slope than the intrinsic relations. Tests of metallicity and stellar mass measurements (see Appendices A and B respectively) show that galaxy stellar masses are typically over-predicted at the low mass end. Furthermore, while the direct method recovers the line-weighted gas-phase metallicity well, it struggles in simulations with dense gas to accurately infer the mass-weighted metallicity, preferentially over-predicting the abundance of each galaxy. Overall, both of these processes will tend to produce a shallower MZR slope. may help to explain the apparent discrepancy in MZR slopes between simulations and observational results.

In summary, we find that all four simulations in the high-redshift suite produce galaxies that fall near the observed low-mass end of the MZR. While it is clear that different SF/feedback physics can have a strong impact on the slope and normalization of the MZR, we find that simulations with dense and turbulent gas (e.g. efficient and bursty star formation) are capable of producing the flat slope which has been observed at low stellar masses (Curti et al. 2024).

### 3.2.2. Temperature Structure of the ISM

In an H II region, the temperatures of different ionization zones can be inferred by various auroral line ratios. However, the relation between the temperatures in each of these zones in high-redshift galaxies is largely unknown.

In Figure 6, we show observed $O^+$ and $O^{++}$ electron temperatures for all galaxies in each run with $M_{UV} \le -15$. On the left, we show the $T_{[O III]}$ as a function of $T_{[O III]}$, colored by electron density as measured by the C3 ratio. Finally, we provide comparisons with observations of H II regions in low redshift galaxies (Mingozzi et al. 2022), local low-metallicity galaxies (Rickards Vaught et al. 2024), intermediate-redshift galaxies (Cataldi et al. 2025) and high-redshift galaxies (Sanders et al. 2025b). We also provide several modeled $T_{[O III]}-T_{[O III]}$ relations from observations (Pilyugin et al. 2009; Cataldi et al. 2025) as well as from an isolated galaxy simulation also run with RAMSES-RTZ (Cameron, Katz & Rey 2023). On the right panel, we show contours describing the loci of MEGATRON galaxies.

We find that there is significant scatter in the $T_{[O III]} - T_{[O III]}$ relation. While low-density, low-ionization galaxies seem to follow the relation from Cameron, Katz & Rey (2023) (suggesting that this is a lower-limit irrespective of feedback prescription), galaxies with higher gas densities tend to scatter up to higher values of $T_{[O III]}$. Furthermore, galaxies with larger gas-phase metallicities tend to shift to lower values of $T_{[O III]}$ and $T_{[O III]}$, consistent with the relative normalizations of each simulation on the MZR in Figure 5 and the expectations of stronger cooling. Galaxies with dense gas are measured to have very high values of $T_{[O III]}$, often greater than $T_{[O III]}$. This is due to the fact that the emission line ratio used to determine this temperature ($O2O2$) is also density dependent, due to the fact that both lines in the [O II]$\lambda$3727 doublet have low critical densities$^8 \lesssim 10^{3.5}$ cm$^{-3}$ (see discussion in Section 3.4). The density range within which this doublet is a suitable density diagnostic will also produce larger values of $O2O2$ for a given temperature - thus over-estimating the true temperature (see Figure 2 of Méndez-Delgado et al. 2023). This underscores the need to carefully fold in density diagnostics when using emission line ratios to infer properties of the ISM at high redshift.

Altogether, we find that the temperature structure of MEGATRON galaxies with low densities is entirely consistent with low-metallicity objects measured at low-redshift and cosmic noon (see discussion in Cataldi et al. 2025). In contrast, the majority of galaxies in MEGATRON have a sufficiently dense ISM that $O2O2$ becomes an unreliable temperature diagnostic without having a density indicator (Méndez-Delgado et al. 2023). Nevertheless, the fact that oxygen abundances at high-redshift are typically dominated by $O^{++}$ (Laseter et al. 2022; Curti et al. 2023) means that any uncertainties due to density are generally subdominant compared to those due to temperature fluctuations (e.g. Cameron, Katz & Rey 2023).

<footnote>$^7$ As gas in very dense regions begins to be collisionally de-excited, emission from lower-density regions begins to dominate.</footnote>

<footnote>$^8$ The same effect will impact the $O3O3$ ratio, though at significantly larger densities.</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 6.— Observed O$^+$ and O$^{++}$ electron temperatures as measured using auroral line ratios ([O II] λλ7325/[O II] λλ3727 and [O III] λ4363/[O III] λ5007 respectively) for all galaxies in MEGATRON with $M_{UV} \le -15$. Left: $T_{[O II]}$ measurements as a function of $T_{[O III]}$ for all galaxies in each MEGATRON simulation, colored by the electron density measured from the C3 ratio. Broadly, we find that for low-density, lower-ionization gas, $T_{[O II]}$ is fairly constant, growing rapidly with increased densities. We provide comparisons with similarly derived observational data of low-redshift analogue galaxies (lime: Mingozzi et al. 2022), local H II regions (green: Rickards Vaught et al. 2024), and intermediate redshift (cyan: Cataldi et al. 2025) and high redshift (black: Sanders et al. 2025b) JWST samples. We also provide modeled relations (green: Pilyugin et al. 2009, blue: Cataldi et al. 2025, salmon: Cameron, Katz & Rey 2023) and the one-to-one relation (black dashed). Right: We show comparisons between contours for each simulation, highlighting the fact that each simulation typically falls in a slightly different $T_{[O III]} - T_{[O II]}$ space.</figcaption></figure>

### 3.3. High Densities

In this section, we explore the third key physical aspect of the high-redshift ISM: investigating the gas density properties of MEGATRON galaxies.

Traditionally, low-ionization state optical lines (e.g. [O II]λλ3726, 3728 and [S II]λλ6719, 6731) have been used to measure the densities of star-forming regions due to their strength and doublet separation. However, the information that these lines provide is inherently incomplete, as this low-ionization gas is not representative of the ISM as a whole, particularly tracing the edge of H II regions and diffuse ionized gas. In contrast, the C III]λλ1906, 1908 and N IV]λλ1483, 1486 doublets are sensitive to increasingly higher energies, tracing denser gas closer to the young, massive stars (e.g. James et al. 2014; Berg et al. 2021; Mingozzi et al. 2022).

Observations of H II regions at $z \gtrsim 1$ with optical lines suggest that they are considerably denser than in the local Universe (e.g. Erb et al. 2006; Hainline et al. 2009; Steidel et al. 2014; Sanders et al. 2016; Berg et al. 2022). JWST now provides access to the rest UV-optical spectra of galaxies at higher redshifts, allowing multi-

tracer determinations of electron densities (Isobe et al. 2023a; Reddy et al. 2023; Abdurro'uf et al. 2024; Topping et al. 2025a; Li et al. 2025b). JWST observations have found elevated gas densities in the early Universe resulting and a strong evolution in gas density with redshift. Likewise a few objects with exceptionally high densities ($\gtrsim 10^5$ cm$^{-3}$, Senchyna et al. 2024; Topping et al. 2024b), were discovered.

Electron densities measured from O2, S2, O3, C3, and N4 for all MEGATRON galaxies with $M_{UV} < -15$ are shown in Figure 7. Each density diagnostic traces a different region of the ISM — therefore, relations between different density measurements are not one-to-one (e.g. Kewley, Nicholls & Sutherland 2019; Méndez-Delgado et al. 2023). Focusing on optical density tracers, we find that $n_e([O II]) \sim n_e([S II])$ (e.g. Kewley, Nicholls & Sutherland 2019; Berg et al. 2021). Next, we find that the O3 FIR tracer generally reports higher densities than that of O2. Given the fact that O3 is sensitive to a lower density range than O2, this suggests that O$^{++}$ exists in systematically denser gas than O$^+$. We note that in more massive metal-enriched galaxies than what is modeled in
<figure><img src="image_8.png" /><figcaption>FIG. 7.— Electron number densities of galaxies with $M_{UV} \le -15$ in each MEGATRON simulation measured using a variety of emission line ratios, including those in the optical (O2 and S2), FIR (O3) and UV (C3 and N4). We find that galaxies of a given simulation produce a rich, inhomogeneous density structure, the exact properties of which are highly dependent on the star-formation/feedback prescriptions used.</figcaption></figure>

MEGATRON, observations show that the reverse may be true Harikane et al. (2025).

Considering UV density diagnostics, we find that both C3- and N4-based density measurements are significantly greater than those traced by O2, in line with observational results (Berg et al. 2022; Mingozi et al. 2022). In particular, the $n_e(\text{C III}) - n_e(\text{O II})$ space yields the greatest differences between each simulation. For instance, we find that while the bursty model exhibits similar C3 and N4 densities, the variable IMF, efficient and high $\epsilon_{ff}$ simulations all find that N4 traces systematically denser regions.

Overall, it is clear that the efficient star formation run produces the densest galaxies, due to the weaker feedback. Nevertheless, it can be instructive to highlight differences between the other models. Beginning with the

high $\epsilon_{ff}$ simulation, the key point is that here, stars are encouraged to form in less dense regions. As a result, these stars will illuminate lower density gas (producing emission lines that trace more diffuse gas) and produce SNe which couple more effectively to the gas$^9$, again reducing the density of the ISM. In contrast, while the bursty and variable IMF runs both feature stronger feedback, these events will occur in denser gas, making the momentum coupling to the ISM comparatively weaker (Blondin et al. 1998; Thornton et al. 1998).

### 3.4. High Ionization Parameters

Here, we explore the final key physical aspect of the high-redshift ISM: investigating the observed ionization

<footnote>$^9$ The presence of hypernovae in this run will also aid to better clear the stellar birth clouds.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 8.— O32-R23 ionization-excitation diagnostic diagram for galaxies in each MEGATRON simulation, with each point colored by the O II-derived electron density. A compilation of observational data is provided including high redshift (black: Cameron et al. 2023b; Mascia et al. 2023; Sanders et al. 2023; Topping et al. 2024b) and low-redshift blueberries (cyan: Yang et al. 2017b) and green peas (green: Yang et al. 2017a). Accounting for metallicity and mass biases (see text), MEGATRON galaxies reproduce the observed scatter well, particularly for galaxies with high values of O32. Interestingly, a large number of galaxies with high O32, low R23, and a dense ISM are predicted, though these tend to be UV-faint. We have highlighted the expected region within which O32 behaves as a density tracer.</figcaption></figure>

parameters of MEGATRON galaxies.

### 3.4.1. *Ionization-Excitation with the O32-R23 Diagram*

Numerous galaxies that have been spectroscopically confirmed at high redshift have been found to produce strong emission lines. Of these, the vast majority appear to have line ratios consistent with high ionization parameters (e.g. Cameron et al. 2023b; Mascia et al. 2023; Sanders et al. 2023; Nakajima et al. 2023; Roberts-Borsani et al. 2024; Topping et al. 2025b; Tang et al. 2025a; Glazer et al. 2025; Hayes et al. 2025), as measured from O32 (Kewley, Nicholls & Sutherland 2019). Furthermore, various studies of stacked spectra have suggested a redshift trend, with galaxies in the earlier Universe having ever higher values of O32 (e.g. Sanders et al. 2023; Roberts-Borsani et al. 2024). Finally, these high ionization parameters are also consistent with the high EWs (e.g. Boyett et al. 2024; Endsley et al. 2024), high

ionizing photon production efficiencies (e.g. Atek et al. 2024; Saxena et al. 2024a; Simmonds et al. 2024a,b; Laseter et al. 2025; Choustikov et al. 2025), and high ionization-state lines (e.g. Bunker et al. 2023; Topping et al. 2024b; Calabrò et al. 2024; Tang et al. 2025b).

The most commonly used emission line ratio diagram to diagnose ionization and excitation is the O32-R23 diagram. In Figure 8 we show O32 as a function of R23 for all galaxies with $M_{UV} \le -15$ in each MEGATRON simulation, colored by the O2-measured electron number density. We also provide a compilation of observed $z \ge 5$ galaxies (Cameron et al. 2023b; Mascia et al. 2023; Sanders et al. 2023; Topping et al. 2024b) as well as the population of low-redshift blueberries (Yang et al. 2017b) and green peas (Yang et al. 2017a) to provide a comparison. First, visual inspection suggests that diagnostic diagrams like O32-R23 depend fundamentally
<figure><img src="image_10.png" /><figcaption>FIG. 9.— O32 as a function of O2-weighted gas density for all four simulations in the high-redshift MEGATRON suite and coloured by the LyC escape fraction. For values of O32 ≳ 1.25, we find a strong correlation between the two, owing to the low critical density of the [O II] λλ3727 lines (indicated by a vertical dashed line). Systems with high values of O32 at low densities are LyC leakers.</figcaption></figure>

on star-formation/feedback prescription, suggesting an exciting avenue to constrain stellar physics in the early Universe (Katz et al. 2024b). We find that the efficient and variable IMF runs are best able to reproduce the observed scatter. In contrast, the bursty and high $\epsilon_{\text{ff}}$ simulations produce galaxies with lower apparent metallicities (see discussion in Section 3.2), which has the effect of shifting points to lower R23. In all cases, MEGATRON struggles to produce a large number of galaxies with high R23 and low O32. Recently, Witten et al. (2025) used data from the SPHINX public data release (Katz et al. 2023a) to show that these are likely to be massive post-starburst systems. While MEGATRON does contain post-starburst galaxies, these are not sufficiently massive or metal-enriched (due to the small simulation volume) to fill this parameter space. We find that all four MEGATRON simulations are capable of producing the scatter of high-O32-high-R23 points observed both in the local universe and at high redshift, including the galaxy with the highest observed O32: RXCJ2248-ID (Topping et al. 2024b). Interestingly, all four MEGATRON simulations predict that there is a large population of low-metallicity, UV faint galaxies with low R23 values and extreme values of O32 (~1000). While these are yet to be observed in the high-redshift universe, their unique electron density properties ($n_e$([O II]) ≳ 1000) motivate us to explore them in more detail.

### 3.4.2. Dependence of O32 on Density

O32 is the most common diagnostics of ionization parameter at high-redshift due to the strength of the lines and their rest-frame wavelengths. However, the critical density of the O2 lines is lower than we often find in MEGATRON galaxies and what is sometimes measured at high-redshift. Figure 9 shows O32 as a function of O2-weighted gas density for galaxies in each MEGATRON simulation, coloured by their angle-averaged LyC escape fraction. Broadly, we find two regimes of behavior. For values of O32 ≲ 1.25, there is no strong relationship between O32 and gas density. In contrast, when O32 ≳ 1.25, we find a strong correlation between O32 and $\langle n_{\text{H}} \rangle_{\text{O2}}$. This is a clear manifestation of the low critical density (indicated by a vertical line) of the O2 lines biasing O32 towards higher values. Finally, we find that galaxies with a large value of O32 at low gas densities tend to be strong LyC leakers, in contrast to those with the same values of O32 at high densities. This is a demonstration of the complex relationship between O32 and $f_{\text{esc}}$ (Nakajima & Ouchi 2014; Katz et al. 2020; Barrow et al. 2020; Flury et al. 2022; Choustikov et al. 2024a), highlighting the dependence of $f_{\text{esc}}$ on the density of the ISM (Choustikov et al. 2024a,b).

Another approach to constrain the ionization parameter might be to use C43 and N43, since these ratios are less sensitive to density. To illustrate this, on the left of Figure 10, we show N43 as a function of C43 for galaxies in each MEGATRON simulation. We find that for the most part, these two ratios are directly proportional, suggesting that both are dependent only on the ionization parameter, to first order. Each simulation is consistent with measured line ratios from the lensed galaxy RXCJ2248-ID (Topping et al. 2024b) as well as other extreme sources (Bunker et al. 2023; Castellano et al. 2024; Naidu et al. 2025). In the right panel of Figure 10 we show O32 as a function of C43 colored by the O2-weighted gas density taken from each simulation, compared again to RXCJ2248-ID. We find that strong correlations between O32 and C43 exist only in a particular gas density range, again confirming that in the vast majority of MEGATRON galaxies, high values of O32 are partially driven by high gas density.

Finally, it is important to consider the effects of turbulence in the ISM. Due to the fact that [O III] λ5007 and [O II] λλ3727 are produced by volume and boundary species respectively, the O32 ratio can be significantly impacted by turbulence (see discussions in Gray & Scanapnieco 2017; Jin, Kewley & Sutherland 2022, Xing et al., in prep.). For instance, introducing turbulence into an initially spherically symmetric Strömgren sphere will increase the surface area of the ionized region for a given volume, thereby boosting [O II] λλ3727 flux for a given [O III] λ5007 luminosity and therefore reducing the O32 ratio at fixed (mean) ionization parameter (Jin, Kewley & Sutherland 2022).

### 3.4.3. Ionization-Excitation with the C43-R3 Diagram

Given the impact of density, temperature, and turbulence on O32, we are inspired to consider a different ionization diagnostic for the ionization-excitation diagram. Ideally, for this purpose we would use N43; due to the involved lines having a higher critical density than those used in O32 or C43; however, significantly more strong carbon emitters have been detected at high-redshift (e.g. Roberts-Borsani et al. 2025), despite the high N/O abundances that have also been measured in the early universe (e.g. Cameron et al. 2023a; Isobe et al. 2023b, 2025; Schaerer et al. 2024; Marques-Chaves et al. 2024; Arellano-Córdova et al. 2025). Therefore, we propose the use of the C43-R3 diagram, which we show in Figure 11 for all MEGATRON galaxies considered, colored by the
<figure><img src="image_11.png" /><figcaption>FIG. 10.— Left: N43 as a function of C43 for all galaxies in each MEGATRON simulation. This acts as a calibration, showing that each of these line ratios depends primarily on the ionization parameter, consistent with observational data (Bunker et al. 2023; Topping et al. 2024b; Castellano et al. 2024; Naidu et al. 2025). Right: O32 as a function of C43 for each galaxy in MEGATRON, colored by the O2-weighted gas density taken from the simulation. Strong correlations between O32 and C43 only exist for a given gas density, suggesting that in the high-redshift Universe, the dependence of O32 on ionization parameter is dominated by a dependence on gas density, due to the low critical densities of the [O II] λλ3727 lines.</figcaption></figure>

direct method gas-phase metallicity. Histograms are pro-
vided to show differences between each of the simulations
in each diagnostic. For comparisons, we provide observa-
tional data from the DJA (DAWN JWST Archive<fnref n="10" />: Heintz
et al. 2024) as well as for other specific objects (Marques-
Chaves et al. 2024; Topping et al. 2024b, 2025b) and
stacks (Hayes et al. 2025). We find very good agree-
ment between oxygen-rich MEGATRON galaxies and the
observed data, implying that even without the density ef-
fects discussed above, high-redshift galaxies exhibit gas
with high ionization parameters. This effect is further
compounded with metallicity, where galaxies with low
oxygen abundances tend to have higher values of C43 for
a given value of R3. Moreover, we find that for the lowest
values of R3, two tails appear at C43 ~ 0 and 1. These
correspond to systems with ionizing photon production
efficiencies of <~25 and >~26 respectively. Finally, many
galaxies with the highest ionization parameters tend to
have very low values of R3, due to the fact that these sys-
tems have low metallicities and/or are sufficiently dense
to begin collisionally de-exciting [O III] λ5007 emission.
This may help explain the particularly weak O III and
strong C IV emission of LAP1-B (Nakajima et al. 2025).

However, it is important to note that there remain
complexities with interpreting C43. First, while the crit-
ical densities of the C III doublet are significantly higher
than those of O II, in principle at sufficiently high densi-

ties a similar effect will take place, with C43 increasing as
a function of density. Furthermore, temperature can also
have an impact on C43, boosting C43 by a maximum of
~0.5 dex at the most extreme H II region temperatures.
Moreover, the C IV doublet is a resonant line, meaning
that high optical depths can bias the line ratio in non-
trivial ways (e.g. Hayes et al. 2025) including potentially
preferentially tracing outflows (e.g. Steidel et al. 2010).
Finally, there is uncertainty in the shape of the atten-
uation law in high-redshift galaxies (Reddy et al. 2020,
2025; McKinney et al. 2025; Sanders et al. 2025a), which
may have an impact on C43. Nevertheless, for ISM con-
ditions more representative of those observed in galaxies
at high-redshift, this diagram should prove very helpful
to discern the properties of galaxies below redshift ~9
where [O III] λ5007 redshifts out of the NIRSpec range.
In those cases, it is possible to use similar diagrams with
O III] λλ1663/Hγ, or [O III] λ4363]/Hγ, though the first
introduces systematics due to dust attenuation and the
latter can be difficult to detect.

4. CAVEATS

As with any simulation, there are numerous caveats
that should be kept in mind when interpreting our re-
sults.

**Numerical Resolution:** MEGATRON is subject to a finite spatial and mass resolution. For the purpose of this work, this can have two main effects. The first is that we are unable to always completely resolve every single

<footnote n="10">¹⁰ https://dawn-cph.github.io/dja/</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 11.— C43-R3 excitation ionization diagnostic diagram. We show data for galaxies with $M_{UV} \le -15$ from all MEGATRON simulations colored by the direct method oxygen abundance. On the top and right, we show histograms of the emission line ratios in question, demonstrating the impact of sub-grid star-formation/feedback prescriptions on these diagnostics. For a comparison, we include all galaxies at $z \ge 6$ in the DJA (Heintz et al. 2024) with signal-to-noise ratios greater than 10 for all lines as well as data from individual objects (Marques-Chaves et al. 2024; Topping et al. 2024b, 2025b) and stacks (Hayes et al. 2025).</figcaption></figure>

Strömgren sphere (and importantly the detailed physics that occurs at the ionization front). These cells will have an incorrect temperature (this can affect predicted emission line luminosities, see discussion in Choustikov et al. 2024a) which is corrected according to the strategy described in Katz et al. (2024b). Secondly, our resolution limits our ability to accurately follow the effects of radiation feedback (see discussions in Kimm et al. 2019, 2022) which may pre-process dense regions of the ISM ahead of SNe. This will be addressed in upcoming work.

**Dust Attenuation:** In producing mock spectroscopic observations of MEGATRON galaxies, we have ignored the effects of dust – in terms of absorption, emission, and scattering. This is likely a good assumption, as we are predominantly focusing on low-mass, low-metallicity objects which are not expected to have a dominant dust component. This will be addressed in future work.

**Missing Physics:** At present, it is important to note that the MEGATRON simulations have not included certain physical processes such as magnetic fields, cosmic rays, and on-the-fly dust formation and destruction – all of which have the potential to impact the properties of the ISM as extra heating/cooling and/or pressure terms. These will be addressed in an upcoming work. Similarly, we have not included a prescription for AGN, which are a source of radiation and feedback. While we

expect these not to dominate at the mass ranges considered in MEGATRON, we do note that some observations have suggested an abundance of over-massive rapidly-accreting black holes (Maiolino et al. 2024, c.f. Li et al. 2025a).

**Sub-Grid Uncertainties:** As with any physical simulation which needs to account for a high dynamic range of scales, we have had to use certain sub-grid prescriptions for stellar populations, formation and feedback. While prescriptions in the efficient star formation run have been shown to match the stellar mass - halo mass relation at $z=0$ (Agertz et al. 2021, see also Section 3 of Katz et al. 2025), the remaining three simulations use prescriptions that are also physically plausible. Given the dependence of ISM properties on low-metallicity sources of ionizing radiation, in future work it will be necessary to include detailed prescriptions of x-ray binaries (e.g. Saxena et al. 2021), Wolf-Rayet stars (e.g. Sander et al. 2025), and AGN (e.g. Kubota & Done 2018; Trebitsch et al. 2021). All of these will contribute by increasing the number of hard ionizing photons, leading to an even more extreme ISM.

A more in-depth discussion of the caveats is provided in Katz et al. 2025. Similarly, an extensive discussion on the impacts of sub-grid modeling on galaxy properties in MEGATRON-like simulations is given in Katz et al.
(2024b).

## 5. CONCLUSIONS

We have presented early results from the MEGATRON simulations (Katz et al. 2025; Rey et al. 2025; Cadiou et al. 2025), investigating the physics behind observed properties of the ISM of galaxies in the high-redshift Universe. The MEGATRON simulations represent an attempt to employ a detailed non-equilibrium thermochemical network of primordial species, metals, and molecules, coupled to on-the-fly radiation transport and galaxy formation physics in a cosmological context. Due to this approach, we are able to predict spectra on a cell-by-cell basis, accounting for the stellar continuum, nebular continuum and nebular emission lines.

We have approached the problem as observers: beginning with the intrinsic spectrum and photometry of each galaxy and measuring galaxy properties using techniques often applied to JWST observations. For example, electron temperatures and densities are measured from emission line ratios using PyNeb (Luridiana, Morisset & Shaw 2015), gas-phase metallicities are measured using the direct (auroral line) method (Cameron, Katz & Rey 2023), and stellar masses are measured by SED-fitting with BAGPIPES (Carnall et al. 2018). This is a unique approach to comparing simulations with observations, that exploits some of the key strengths of the MEGATRON simulations.

Our results can be summarized as follows:

1. The ISM of MEGATRON galaxies is fundamentally different to that of observed low-redshift galaxies. It is denser, less metal-enriched, subject to harder radiation fields, and has a higher ionization parameter. This is in agreement with most JWST observations (e.g. Bunker et al. 2023; Cameron et al. 2023b; Sanders et al. 2023; Mascia et al. 2023; Curti et al. 2024; Hayes et al. 2025; Tang et al. 2025a; Pollock et al. 2025).

2. No two emission lines are produced by gas of the same density and temperature. As a direct consequence of the turbulent, multi-phase structure of high-redshift galaxies, this can significantly bias measurements of ISM properties (e.g. Cameron, Katz & Rey 2023; Harikane et al. 2025; Usui et al. 2025).

3. Simulations that produce a very dense ISM are observed to have a shallow MZR slope, consistent with what is observed at high-redshift (e.g. Heintz et al. 2023; Nakajima et al. 2023; Curti et al. 2023, 2024; Morishita et al. 2024; Chemerynska et al. 2024a; Sarkar et al. 2025). Observational techniques designed to measure galaxy masses and metallicities systematically bias the MZR to be flatter at the low-mass end.

4. ISM densities inferred by UV tracers are systematically larger than those measured by optical lines, in agreement with observations (Isobe et al. 2023a; Reddy et al. 2023; Abdurro'uf et al. 2024; Li et al. 2025b; Topping et al. 2025a). Interestingly, we find that O III FIR lines also trace denser gas than

O II. The exact density structure is sensitive to the star-formation/feedback model assumed via the strength of feedback, IMF choice, and the ambient density during stellar feedback.

5. At the densities reached in the ISM of high-redshift galaxies, [O III] $\lambda5007/[O II]$ $\lambda3727$ becomes a density tracer as [O II] is collisionally de-excited. If not properly taken into account, measurements of ionization parameter from O32 will be artificially high. We instead propose the use of C IV $\lambda1549/C III$ $\lambda1907$ or N IV $\lambda1484/N III$ $\lambda1750$. When considering these diagnostics, MEGATRON galaxies match observations very well, suggesting that the ionization parameter of observed galaxies is indeed very high (e.g. Marques-Chaves et al. 2024; Topping et al. 2024b, 2025b; Hayes et al. 2025).

6. Galaxies simulated with various feedback models populate different regions of strong-line diagnostic diagrams as line ratios are sensitive to the density-temperature structure of the ISM (Katz et al. 2024b). This suggests a possible avenue for constraining high-redshift feedback physics with line ratios directly.

This work has demonstrated the ability of simulations with on-the-fly out-of-equilibrium thermochemistry like MEGATRON to help interpret the diversity of emission line properties observed by JWST. Likewise, when large samples of JWST galaxies with a well-defined selection function become public, comparisons with simulations such as MEGATRON will be ideal for to constraining the physics of star formation as well as the impact of stellar feedback on the ISM of high-redshift galaxies.

## ACKNOWLEDGMENTS

NC acknowledges support from the Science and Technology Facilities Council (STFC) for a PhD studentship. AJC and AS acknowledge funding from the “FirstGalaxies” Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No.789056). TK is supported by the National Research Foundation of Korea (RS-2022-NR070872 and RS-2025-00516961) and by the Yonsei Fellowship, funded by Lee Youn Jae. KM acknowledges the Flemish Fund for Scientific Research (FWO-Vlaanderen), Grant number 1169822N.

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. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC op-
erations grant ST/R000832/1. This work was per-
formed using resources provided by the Cambridge
Service for Data Driven Discovery (CSD3) operated
by the University of Cambridge Research Computing
Service (www.csd3.cam.ac.uk), provided by Dell EMC
and Intel using Tier-2 funding from the Engineering
and Physical Sciences Research Council (capital grant
EP/T022159/1), and DiRAC funding from the Science
and Technology Facilities Council (www.dirac.ac.uk).
DiRAC is part of the National e-Infrastructure. The
authors thank Jonathan Patterson for smoothly run-
ning the Glamdring Cluster hosted by the University
of Oxford, where part of the data processing was per-
formed. The authors also acknowledge financial support
from Oriel College's Research Fund. The material in this

manuscript is based upon work supported by NASA un-
der award No. 80NSSC25K7009.

The authors thank Romain Teyssier and Leo Michel-
Dansac for both developing and open-sourcing RAMSES
and RASCAS, respectively. We thank the developers and
maintainers of PYNBODY (Pontzen et al. 2013), YT (Turk
et al. 2011), NUMPY (van der Walt, Colbert & Varo-
quaux 2011; Harris et al. 2020), SciPy (Virtanen et
2020), JUPYTER (Ragan-Kelley et al. 2014), MATPLOTLIB
(Hunter 2007), CMASHER (van der Velden 2020), PANDAS
(pandas development team 2020; Wes McKinney 2010),
the Astrophysics Data Service, and the arXiv pre-print
repository for providing open-source software and ser-
vices that were used extensively in this work.

REFERENCES

Abdurro'uf et al., 2024, ApJ, 973, 47

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

Andrews B. H., Martini P., 2013, ApJ, 765, 140

Arellano-Córdova K. Z. et al., 2022, ApJ, 940, L23

Arellano-Córdova K. Z. et al., 2025, MNRAS, 540, 2991

Atek H. et al., 2024, Nature, 626, 975

Barrow K. S. S., Robertson B. E., Ellis R. S., Nakajima K.,
Saxena A., Stark D. P., Tang M., 2020, ApJ, 902, L39

Berg D. A., Chisholm J., Erb D. K., Skillman E. D., Pogge
R. W., Olivier G. M., 2021, ApJ, 922, 170

Berg D. A. et al., 2022, ApJS, 261, 31

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

Böker T. et al., 2023, PASP, 135, 038001

Boyett K. et al., 2024, MNRAS, 535, 1796

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

Buck T., Rybiczki J., Buder S., Obreja A., Macciò A. V.,
Pfrommer C., Steinmetz M., Ness M., 2021, MNRAS, 508, 3365

Bunker A. J. et al., 2023, A&A, 677, A88

Calabrò A. et al., 2024, ApJ, 975, 245

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

Cameron A. J., Katz H., Rey M. P., Saxena A., 2023a, MNRAS,
523, 3516

Cameron A. J., Katz H., Witten C., Saxena A., Laporte N.,
Bunker A. J., 2024, MNRAS, 534, 523

Cameron A. J. et al., 2023b, A&A, 677, A115

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018,
MNRAS, 480, 4379

Castellano M. et al., 2024, ApJ, 972, 143

Cataldi E. et al., 2025, arXiv e-prints, arXiv:2504.03839

Chatzikos M. et al., 2023, Rev. Mexicana Astron. Astrofis., 59,
327

Chemerynska I. et al., 2024a, ApJ, 976, L15

Chemerynska I. et al., 2024b, MNRAS, 531, 2615

Chevallard J., Charlot S., 2016, MNRAS, 462, 1415

Choe J., Kimm T., Katz H., Rey M., Han D., Jang J. K., Rosdahl
J., 2025, arXiv e-prints, arXiv:2510.00400

Choustikov N. et al., 2024a, MNRAS, 529, 3751

Choustikov N. et al., 2024b, MNRAS, 532, 2463

Choustikov N., Stiskalek R., Saxena A., Katz H., Devriendt J.,
Slyz A., 2025, MNRAS, 537, 2273

Cochrane R. K., Katz H., Begley R., Hayward C. C., Best P. N.,
2025, ApJ, 978, L42

Cueto E. R., Hutter A., Dayal P., Gottlöber S., Heintz K. E.,
Mason C., Trebitsch M., Yepes G., 2024, A&A, 686, A138

Cullen F. et al., 2025, MNRAS, 540, 2176

Cullen F. et al., 2024, MNRAS, 531, 997

Curti M., Cresci G., Mannucci F., Marconi A., Maiolino R.,
Esposito S., 2017, MNRAS, 465, 1384

Curti M. et al., 2023, MNRAS, 518, 425

Curti M. et al., 2024, A&A, 684, A75

Curti M., Mannucci F., Cresci G., Maiolino R., 2020, MNRAS,
491, 944

Dekel A., Sarkar K. C., Birnboim Y., Mandelker N., Li Z., 2023,
MNRAS, 523, 3201

Del Zanna G., Dere K. P., Young P. R., Landi E., 2021, ApJ, 909,
38

Dere K. P., Del Zanna G., Young P. R., Landi E., Sutherland
R. S., 2019, ApJS, 241, 22

Dere K. P., Landi E., Mason H. E., Monsignori Fossi B. C.,
Young P. R., 1997, A&AS, 125, 149

Donnan C. T. et al., 2023, MNRAS, 518, 6011

Dopita M. A., Sutherland R. S., 1996, ApJS, 102, 161

D'Silva J. C. J., Lagos C. D. P., Davies L. J. M., Lovell C. C.,
Vijayan A. P., 2023, MNRAS, 518, 456

Dubois Y. et al., 2013, A&A, 651, A109

Ekström S. et al., 2012, A&A, 537, A146

Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S.,
Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, PASA, 34,
e058

Ellison S. L., Patton D. R., Simard L., McConnachie A. W., 2008,
ApJ, 672, L107

Endsley R. et al., 2024, MNRAS, 533, 1111

Erb D. K., Shapley A. E., Pettini M., Steidel C. C., Reddy N. A.,
Adelberger K. L., 2006, ApJ, 644, 813

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

Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354

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

Ferrara A., Manzoni D., Ntormousi E., 2025, arXiv e-prints,
arXiv:2509.02566

Finkelstein S. L. et al., 2023, ApJ, 946, L13

Flury S. R. et al., 2022, ApJ, 930, 126

Gardner J. P. et al., 2023, PASP, 135, 068001

Giovinazzo E. et al., 2025, arXiv e-prints, arXiv:2507.01096

Glazer K. S. et al., 2025, arXiv e-prints, arXiv:2504.21080

Götberg Y., de Mink S. E., Groh J. H., Leitherer C., Norman C.,
2019, A&A, 629, A134

Gray W. J., Scannapieco E., 2017, ApJ, 849, 132

Gutkin J., Charlot S., Bruzual G., 2016, MNRAS, 462, 1757

Hainline K. N., Shapley A. E., Kornei K. A., Pettini M.,
Buckley-Geer E., Allam S. S., Tucker D. L., 2009, ApJ, 701, 52

Harikane Y., Nakajima K., Ouchi M., Umeda H., Isobe Y., Ono
Y., Xu Y., Zhang Y., 2024, ApJ, 960, 56

Harikane Y. et al., 2025, arXiv e-prints, arXiv:2505.09186

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

Harvey T. et al., 2025, arXiv e-prints, arXiv:2504.05244

Hayes M. J., Saldana-Lopez A., Citro A., James B. L., Mingozzi
M., Scarlata C., Martinez Z., Berg D. A., 2025, ApJ, 982, 14

Hegde S., Wyatt M. M., Furlanetto S. R., 2024, J. Cosmology
Astropart. Phys., 2024, 025

Heintz K. E. et al., 2023, Nature Astronomy, 7, 1517

Heintz K. E. et al., 2025, A&A, 693, A60

Heintz K. E. et al., 2024, Science, 384, 890

Hirschmann M. et al., 2023, MNRAS, 526, 3610

Hirschmann M., Charlot S., Feltre A., Naab T., Choi E., Ostriker
J. P., Somerville R. S., 2017, MNRAS, 472, 2468

Hu W. et al., 2024, ApJ, 971, 21

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

Isobe Y. et al., 2025, MNRAS, 541, L71

Isobe Y., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y.,
Zhang Y., Umeda H., 2023a, ApJ, 956, 139

Isobe Y. et al., 2023b, ApJ, 959, 100

Iyer K. G., Pacifici C., Calistro-Rivera G., Lovell C. C., 2025,
arXiv e-prints, arXiv:2502.17680

Izotov Y. I., Stasińska G., Meynet G., Guseva N. G., Thuan
T. X., 2006, A&A, 448, 955

James B. L. et al., 2014, MNRAS, 440, 1794

Jin Y., Kewley L. J., Sutherland R. S., 2022, ApJ, 934, L8

Johnson B. D., 2019, SEDPY: Modules for storing and operating
on astronomical source spectral energy distribution.
Astrophysics Source Code Library, record ascl:1905.026

Katz H., 2022, MNRAS, 512, 348

Katz H. et al., 2025, The Open Journal of Astrophysics, 8, 104

Katz H. et al., 2019a, MNRAS, 487, 5902

Katz H. et al., 2022a, MNRAS, 515, 4265
Katz H., Ji A. P., Telford G., Senchyna P., 2024a, The Open Journal of Astrophysics, 7, 106

Katz H., Laporte N., Ellis R. S., Devriendt J., Slyz A., 2019b, MNRAS, 484, 4054

Katz H. et al., 2022b, arXiv e-prints, arXiv:2211.04626

Katz H., Rey M. P., Cadiou C., Kimm T., Agertz O., 2024b, arXiv e-prints, arXiv:2411.07282

Katz H. et al., 2023a, The Open Journal of Astrophysics, 6, 44

Katz H. et al., 2022c, MNRAS, 510, 5603

Katz H. et al., 2023b, MNRAS, 518, 592

Katz H. et al., 2020, MNRAS, 498, 164

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

Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511

Kimm T., Bieri R., Geen S., Rosdahl J., Blaizot J., Michel-Dansac L., Garel T., 2022, ApJS, 259, 21

Kimm T., Blaizot J., Garel T., Michel-Dansac L., Katz H., Rosdahl J., Verhamme A., Haehnelt M., 2019, MNRAS, 486, 2215

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Kobayashi C., Taylor P., 2023, arXiv e-prints, arXiv:2302.07255

Kobayashi C., Umeda H., Nomoto K., Tominaga N., Ohkubo T., 2006, ApJ, 653, 1145

Kravtsov A., Belokurov V., 2024, arXiv e-prints, arXiv:2405.04578

Kroupa P., 2001, MNRAS, 322, 231

Kubota A., Done C., 2018, MNRAS, 480, 1247

Langan I., Ceverino D., Finlabor K., 2020, MNRAS, 494, 1988

Laporte N., Nakajima K., Ellis R. S., Zitrin A., Stark D. P., Mainali R., Roberts-Borsani G. W., 2017, ApJ, 851, 40

Laseter I. H., Barger A. J., Cowie L. L., Taylor A. J., 2022, ApJ, 935, 150

Laseter I. H. et al., 2024, A&A, 681, A70

Laseter I. H. et al., 2025, ApJ, 988, 73

Lebouteiller V., Ramab践on L., 2022, A&A, 667, A34

Lee H., Skillman E. D., Cannon J. M., Jackson D. C., Gehrz R. D., Polomski E. F., Woodward C. E., 2006, ApJ, 647, 970

Leitherer C. et al., 1999, ApJS, 123, 3

Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3

Lequeux J., Peimbert M., Rayo J. F., Serrano A., Torres-Peimbert S., 1979, A&A, 80, 155

Leung G. C. K. et al., 2023, ApJ, 954, L46

Leung G. C. K. et al., 2024, arXiv e-prints, arXiv:2411.12005

Li J. et al., 2025a, ApJ, 981, 19

Li S. et al., 2025b, ApJ, 979, L13

Lovell C. C., Roper W. J., Vijayan A. P., Wilkins S. M., Newman S., Seeyave L., 2025, arXiv e-prints, arXiv:2508.03888

Lovell C. C. et al., 2024, arXiv e-prints, arXiv:2411.13960

Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42

Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, MNRAS, 456, 2140

Maiolino R., Mannucci F., 2019, A&A Rev., 27, 3

Maiolino R. et al., 2024, A&A, 691, A145

Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

Marconi A. et al., 2024, A&A, 689, A78

Marks M., Kroupa P., Dabringhausen J., Pawlowski M. S., 2012, MNRAS, 422, 2246

Marques-Chaves R. et al., 2024, A&A, 681, A30

Marszewski A., Sun G., Faucher-Giguère C.-A., Hayward C. C., Feldmann R., 2024, ApJ, 967, L41

Mascia S. et al., 2023, A&A, 672, A155

Maseda M. V. et al., 2020, MNRAS, 493, 5120

Mason C. A., Trenti M., Treu T., 2023, MNRAS, 521, 497

Matthee J. et al., 2024, ApJ, 963, 129

McKinney J., Cooper O. R., Casey C. M., Muñoz J. B., Akins H., Lambrides E., Long A. S., 2025, ApJ, 985, L21

Méndez-Delgado J. E. et al., 2023, MNRAS, 523, 2952

Michel-Dansac L., Blaizot J., Garel T., Verhamme A., Kimm T., Trebitsch M., 2020, A&A, 635, A154

Mingozzi M. et al., 2022, ApJ, 939, 110

Moe M., Kratter K. M., Badenes C., 2019, ApJ, 875, 61

Morishita T. et al., 2024, ApJ, 971, 43

Naidu R. P. et al., 2025, arXiv e-prints, arXiv:2505.11263

Nakajima K., Ouchi M., 2014, MNRAS, 442, 900

Nakajima K. et al., 2025, arXiv e-prints, arXiv:2506.11846

Nakajima K., Ouchi M., Isobe Y., Harikane Y., Zhang Y., Ono Y., Umeda H., Oguri M., 2023, ApJS, 269, 33

Narayanan D. et al., 2024, ApJ, 961, 73

Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nucl. Phys. A, 777, 424

Nyhagen C. T., Schimek A., Cione C., Decataldo D., Shen S., 2024, arXiv e-prints, arXiv:2410.18471

Oppenheimer B. D., Schaye J., 2013, MNRAS, 434, 1043

Padoan P., Nordlund Å., 2011, ApJ, 741, L22

Pallottini A. et al., 2022, MNRAS, 513, 5621

pandas development team T., 2020, pandas-dev/pandas: Pandas Peimbert M., 1967, ApJ, 150, 825

Pilyugin L. S., Mattsson L., Vílchez J. M., Cedrés B., 2009, MNRAS, 398, 485

Pilyugin L. S., Vílchez J. M., Mattsson L., Thuan T. X., 2012, MNRAS, 421, 1624

Ploeckinger S., Richings A. J., Schaye J., Trayford J. W., Schaller M., Chaikin R., 2025, MNRAS

Pollock C. L. et al., 2025, arXiv e-prints, arXiv:2506.15779

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

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

Reddy N. A. et al., 2020, ApJ, 902, 123

Reddy N. A. et al., 2025, arXiv e-prints, arXiv:2506.17396

Reddy N. A., Topping M. W., Sanders R. L., Shapley A. E., Brammer G., 2023, ApJ, 952, 167

Rémy-Ruyer A. et al., 2014, A&A, 563, A31

Richings A. J., Faucher-Giguère C.-A., Gurvich A. B., Schaye J., Hayward C. C., 2022, MNRAS, 517, 1557

Rickards Vaughn R. J. et al., 2024, ApJ, 966, 130

Roberts-Borsani G. et al., 2025, arXiv e-prints, arXiv:2508.21708

Roberts-Borsani G. et al., 2024, ApJ, 976, 193

Robertson B. E., 2022, ARA&A, 60, 121

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

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

Safranek-Shrader C., Milosavljević M., Bromm V., 2014, MNRAS, 438, 1669

Sana H. et al., 2012, Science, 337, 444

Sander A. A. C. et al., 2025, arXiv e-prints, arXiv:2508.18410

Sanders R. L. et al., 2016, ApJ, 816, 23

Sanders R. L. et al., 2025a, ApJ, 989, 209

Sanders R. L. et al., 2025b, arXiv e-prints, arXiv:2508.10099

Sanders R. L., Shapley A. E., Topping M. W., Reddy N. A., Brammer G. B., 2023, ApJ, 955, 54

Sarkar A. et al., 2025, ApJ, 978, 136

Saxena A. et al., 2024a, A&A, 684, A84

Saxena A. et al., 2024b, arXiv e-prints, arXiv:2411.14532

Saxena A. et al., 2021, MNRAS, 505, 4798

Schaerer D., Marques-Chaves R., Xiao M., Korber D., 2024, A&A, 687, L11

Schmidt M., 1959, ApJ, 129, 243

Senchyna P., Plat A., Stark D. P., Rudie G. C., Berg D., Charlot S., James B. L., Mingozzi M., 2024, ApJ, 966, 92

Shapley A. E. et al., 2025, ApJ, 980, 242

Shen X., Vogelsberger M., Boylan-Kolchin M., Tacchella S., Kannan R., 2023, MNRAS, 525, 3254

Simmons C. et al., 2024a, MNRAS, 527, 6139

Simmons C. et al., 2024b, MNRAS, 535, 2998

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stanway E. R., Eldridge J. J., 2019, A&A, 621, A105

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

Stark D. P. et al., 2015, MNRAS, 454, 1393

Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289

Steidel C. C. et al., 2014, ApJ, 795, 165

Sun G., Faucher-Giguère C.-A., Hayward C. C., Shen X., Wetzel A., Cochrane R. K., 2023, ApJ, 955, L35

Sutherland R., Dopita M., Binette L., Groves B., 2018, MAPPINGS V: Astrophysical plasma modeling code. Astrophysics Source Code Library, record ascl:1807.005

Tacchella S. et al., 2022, ApJ, 927, 170

Tang M., Stark D. P., Mason C. A., Gelli V., Chen Z., Topping M. W., 2025a, arXiv e-prints, arXiv:2507.08245

Tang M. et al., 2025b, arXiv e-prints, arXiv:2505.06359

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

Thornton K., Gaudlitz M., Janka H. T., Steinmetz M., 1998, ApJ, 500, 95

Topping M. W. et al., 2025a, arXiv e-prints, arXiv:2502.08712

Topping M. W., Shull J. M., 2015, ApJ, 800, 97

Topping M. W. et al., 2024a, MNRAS, 529, 4087

Topping M. W. et al., 2025b, ApJ, 980, 225

Topping M. W. et al., 2024b, MNRAS, 529, 3301

Torrey P. et al., 2019, MNRAS, 484, 5587

Trebitsch M. et al., 2021, A&A, 653, A154

Tremonti C. A. et al., 2004, ApJ, 613, 898

Trinca A., Schneider R., Valiante R., Graziani L., Ferrotti A., Omukai K., Chon S., 2024, MNRAS, 529, 3563

Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, ApJS, 192, 9

Ucci G. et al., 2023, MNRAS, 518, 3557

Usui M. et al., 2025, arXiv e-prints, arXiv:2507.02054
<figure><img src="image_13.png" /></figure>

FIG. 12.— Recovery test for gas-phase oxygen abundance inferred through the direct method (e.g. Peimbert 1967; Cameron, Katz & Rey 2023) compared to the line-weighted (top) and mass-weighted (bottom) measured for each galaxy in the simulation. We find that generally, the direct method performs well at recovering the line-weighted oxygen abundance. In contrast, in both the efficient and bursty star formation simulations, the direct method tends to over-predict the mass-weighted oxygen abundance of each galaxy. This is primarily due to fluctuations in the temperature, density and abundance structure of the ISM.

van der Velden E., 2020, The Journal of Open Source Software, 5, 2004

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

Vijayan A. P., Lovell C. C., Wilkins S. M., Thomas P. A., Barnes D. J., Irodotou D., Kuusisto J., Roper W. J., 2021, MNRAS, 501, 3289

Vijayan A. P. et al., 2025, arXiv e-prints, arXiv:2507.20190

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

Wes McKinney, 2010, in Proceedings of the 9th Python in Science Conference, Stéfan van der Walt, Jarrod Millman, eds., pp. 56 – 61

Wilkins S. M. et al., 2023, MNRAS, 522, 4014

Witten C. et al., 2025, MNRAS, 537, 112

Xiao L., Stanway E. R., Eldridge J. J., 2018, MNRAS, 477, 904

Yang H. et al., 2017a, ApJ, 844, 171

Yang H., Malhotra S., Rhoads J. E., Wang J., 2017b, ApJ, 847, 38

Yates R. M., Kauffmann G., Guo Q., 2012, MNRAS, 422, 215

Yoon S. C., Langer N., 2005, A&A, 443, 643

Yung L. Y. A., Somerville R. S., Finkelstein S. L., Wilkins S. M., Gardner J. P., 2024, MNRAS, 527, 5929

Zubko V., Dwek E., Arendt R. G., 2004, ApJS, 152, 211

# APPENDIX

## A. GAS-PHASE METALLICITY RECOVERY WITH THE DIRECT METHOD

In Figure 12, we show the gas-phase metallicities inferred for galaxies in each MEGATRON simulation as a function of the line-weighted metallicity (top) and mass-weighted metallicity (bottom) taken from the simulation. We find that the direct method traces the line-weighted values very well, with residuals correlating with the S2-weighted gas density of each galaxy. This is due to the fact that under-estimating the density of the system tends to under-predict the overall abundance of the gas. In contrast, we find that the direct method tends to significantly over-predict oxygen abundances in galaxies with a denser ISM. This is due to the fact that in these systems, emission lines are particularly biased to dense, star-forming regions which may be unrepresentative of the galaxy as a whole.

## B. STELLAR MASS RECOVERY WITH BAGPIPES

In Figure 13, we show recovery tests for the stellar masses of MEGATRON galaxies, as measured using BAGPIPES (Carnall et al. 2018). In particular, we show predicted stellar masses as a function of true (after mass loss) stellar masses for each simulation, colored by SFR<sub>10</sub>/SFR<sub>100</sub> (an indicator of burstiness). We find that below $z \lesssim 10.5$, BAGPIPES performs well, although it tends to over-predict the masses of the least massive galaxies. Above $z \gtrsim 10.5$, it tends to struggle considerably more, with the residuals correlating with burstiness. This is primarily due to the fact that at these redshifts, the Balmer break/jump redshifts out of the red-most filter used – suggesting that the inclusion of MIRI bands may help to alleviate this issue.
<figure><img src="image_14.png" /></figure>

FIG. 13.— Recovery test for stellar masses inferred through SED-fitting by BAGPIPES (Carnall et al. 2018) compared to their values taken from each simulation in two redshift bins. Each point is colored by the burstiness in the recent star formation history given by SFH<sub>10</sub>/SFH<sub>100</sub>. We also provide running medians and one-to-one relations in black and blue, respectively. We find that at lower redshifts, masses are recovered fairly well while at higher redshifts this approach struggles, due to the fact that the Balmer jump has shifted out of the last filter used. In all cases, we find that galaxy with more bursty recent star formation tend to have over-predicted stellar masses.

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

FIG. 14.— Average number of supernovae for a Kroupa initial mass function (Kroupa 2001) as a function of high-mass slope, assuming that either all stars with $M_{*} > 8 M_{\odot}$ explode (blue), or that only stars with $8 M_{\odot} < M_{*} < 25 M_{\odot}$ do (red).

### C. EXPECTED NUMBER OF SUPERNOVAE FOR A VARIABLE IMF

As one of the MEGATRON simulations uses a variable IMF, it is instructive to demonstrate the non-trivial impact that this can have on SN feedback. In Figure 14, we show the expected number of SN as a function of high-mass IMF slope for two assumptions about which stars will explode (as opposed to collapsing into a black hole, releasing little energy). Here, we see that for a sufficiently shallow slope, the expected number of SN will actually begin to decrease, due to the fact that many of the massive stars are either collapsing directly to black holes. Furthermore, the presence
of a massive star will directly impact the potential number of intermediate mass SNe. Therefore, the presence of a variable, or top-heavy IMF does not guarantee that typical SN feedback will be stronger.

This paper was built using the Open Journal of Astrophysics LATEX template. The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike. Learn more at http://astro.theoj.org.