# The nature of high [O iii] 88 µm/[C ii] 158 µm galaxies in the epoch of reionization: Low carbon abundance and a top-heavy IMF?

Harley Katz, Joakim Rosdahl, Taysun Kimm, Thibault Garel, Joakim A Rosdahl, Martin Haehnelt, Léo Michel-Dansac, Sergio Martin-Alvarez, Julien Devriendt, Adrianne Slyz, et al.

► To cite this version:

Harley Katz, Joakim Rosdahl, Taysun Kimm, Thibault Garel, Joakim A Rosdahl, et al. The nature of high [O iii] 88 µm/[C ii] 158 µm galaxies in the epoch of reionization: Low carbon abundance and a top-heavy IMF?. Monthly Notices of the Royal Astronomical Society, 2022, 510 (4), pp.5603-5622. 10.1093/mnras/stac028. hal-03541987

HAL Id: hal-03541987

https://hal.science/hal-03541987

Submitted on 19 Mar 2023

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
# The nature of high [O III]88 µm/[C II]158 µm galaxies in the epoch of reionization: Low carbon abundance and a top-heavy IMF?

Harley Katz,<fnref n="1" />★ Joakim Rosdahl,<fnref n="2" /> Taysun Kimm,<fnref n="3" /> Thibault Garel,<fnref n="4" /> Jérémy Blaizot,<fnref n="2" />

Martin G. Haehnelt,<fnref n="5" /> Léo Michel-Dansac,<fnref n="2" /> Sergio Martin-Alvarez,<fnref n="5" />,<fnref n="6" /> Julien Devriendt,<fnref n="1" /> Adrianne Slyz,<fnref n="1" /> Romain Teyssier,<fnref n="6" /> Pierre Ocvirk,<fnref n="7" /> Nicolas Laporte,<fnref n="5" /> and Richard Ellis<fnref n="8" />

Accepted 2021 December 23. Received 2021 December 22; in original form 2021 July 26

## ABSTRACT

ALMA observations of $z > 6$ galaxies have revealed abnormally high [O III]88 µm/[C II]158 µm ratios and [C II]158 µm deficits compared to local galaxies. The origin of this behaviour is unknown. Numerous solutions have been proposed including differences in C and O abundance ratios, observational bias, and differences in ISM properties, including ionization parameter, gas density, or photodissociation region (PDR) covering fraction. In order to elucidate the underlying physics that drives this high-redshift phenomenon, we employ SPHINX²⁰, a state-of-the-art, cosmological radiation–hydrodynamics simulation, that resolves detailed ISM properties of thousands of galaxies in the epoch of reionization which has been post-processed with CLOUDY to predict emission lines. We find that the observed $z > 6$ [O III]88 µm–SFR and [C II]158 µm–SFR relations can only be reproduced when the C/O abundance ratio is ~8 × lower than Solar and the total metal production is ~4 × higher than that of a Kroupa IMF. This implies that high-redshift galaxies are potentially primarily enriched by low-metallicity core–collapse supernovae with a more top-heavy IMF. As AGB stars and type-Ia supernovae begin to contribute to the galaxy metallicity, both the [C II]158 µm–SFR and [C II]158 µm luminosity functions are predicted to converge to observed values at $z \sim 4.5$. While we demonstrate that ionization parameter, LyC escape fraction, ISM gas density, and CMB attenuation all drive galaxies towards higher [O III]88 µm/[C II]158 µm, observed values at $z > 6$ can only be reproduced with substantially lower C/O abundances compared to Solar. The combination of [C II]158 µm and [O III]88 µm can be used to predict the values of ionization parameter, ISM gas density, and LyC escape fraction and we provide estimates of these quantities for nine observed $z > 6$ galaxies. Finally, we demonstrate that [O I]63 µm can be used as a replacement for [C II]158 µm in high-redshift galaxies where [C II]158 µm is unobserved and argue that more observation time should be used to target [O I]63 µm at $z > 6$. Future simulations will be needed to self-consistently address the numerous uncertainties surrounding a varying IMF at high redshift and the associated metal returns.

**Key words:** ISM: kinematics and dynamics – ISM: lines and bands – galaxies: high-redshift – galaxies: ISM – galaxies: star formation – dark ages, reionization, first stars.

## 1 INTRODUCTION

Characterizing the interstellar medium (ISM) of high-redshift galaxies deep into the epoch of reionization (EoR) represents one of the most interesting frontiers in the observational study of galaxy formation. Deep ALMA observations of forbidden infrared (IR) emission lines have now been confirmed out to $z > 9$ (e.g. Hashimoto et al. 2018; Laporte et al. 2019), providing insight into both the ISM properties and star formation histories of these primordial galaxies.

As the sample of high-redshift systems with confirmed emission line detections increases, certain differences have been discovered between reionization-epoch galaxies and those in the local Universe.

The [C II]158 µm line (hereafter [C II]), is one of the dominant coolants of the ISM and hence one of the brightest IR emission lines in the local Universe (e.g. Spitzer 1978; Stacey et al. 1991). With an ionization potential of 11.3 eV, [C II] emission can be stimulated in both ionized and neutral gas, often complicating its interpretation. Despite the wide range of ISM conditions that can emit [C II], it has been empirically shown to strongly correlate with star formation in local galaxies (e.g. Boselli et al. 2002; De Looze et al. 2011, 2014; Herrera-Camus et al. 2015). Furthermore, observations have shown

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

<footnote n="2">²Univ Lyon, Univ Lyon 1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="3">³Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="4">⁴Observatoire de Geneve, Universite de Geneve, 51 Ch. des Maillettes, CH-1290 Versoix, Switzerland</footnote>

<footnote n="5">⁵Kavli Institute for Cosmology and Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="6">⁶Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland</footnote>

<footnote n="7">⁷Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

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

<footnote>* E-mail: harley.katz@physics.ox.ac.uk</footnote>
that the [C II] line is especially bright in low-mass metal-poor dwarf galaxies (e.g. Poglitsch et al. 1995; Cormier et al. 2010), making it a particularly interesting target for studying the high-redshift Universe.

Earlier IR observations of $z > 6$ galaxies focused primarily on this line due to its predicted high luminosity and observability with ALMA. There are now more than 20 galaxies at $z > 6$ with [C II] detections as well as many more with non-detections. At slightly lower redshift (e.g. $4 \le z \le 6$), even larger samples of more than 100 galaxies have been surveyed (Béthermin et al. 2020; Le Fèvre et al. 2020), of which 75 have been detected in [C II] emission (Schaerer et al. 2020). The non-detections are of particular interest because this places many of the EoR galaxies far below the local relation between [C II] and the star formation rate (SFR) (e.g. Ota et al. 2014; Maiolino et al. 2015; Schaerer et al. 2015; Knudsen et al. 2016; Pentericci et al. 2016; Bradač et al. 2017; Laporte et al. 2019). The exact origin of these deficits is currently unknown although some have proposed that the deficits are due to bursts of star formation (e.g. Ferrara et al. 2019). However recent observations of extended, low surface brightness [C II] emission in the early Universe (Fujimoto et al. 2019) has helped us relieve some of the tension between the high-redshift galaxies and the low-redshift relations (Carniani et al. 2020) in terms of absolute luminosity, but not surface brightness.

The [O III]$_{88 \mu m}$ line (hereafter [O III]) is another strong IR emission line that has also been empirically established to correlate with star formation (e.g. De Looze et al. 2014). With an ionization potential of 35.1eV, [O III] only traces ionized gas, in contrast to [C II]. Interestingly, in metal-poor local Universe dwarf galaxies, [O III] often outshines [C II] (e.g. Cormier et al. 2010; Madden et al. 2013), representing up to 2 per cent of the total IR luminosity of the galaxy (Cormier et al. 2015). Since low-redshift metal-poor dwarf galaxies are often considered analogues of high-redshift systems, [O III] emission has been a particularly interesting target for the high-redshift Universe (e.g. Inoue et al. 2016) and has already been detected at $z > 9$ (Hashimoto et al. 2018).

While sample sizes of [O III] emitting galaxies remain small at $z > 6$, there have now been ~10 detections, nearly all of which result in exceptionally high ratios of [O III]/[C II], even compared to local metal-poor dwarf galaxies (e.g. Laporte et al. 2019; Carniani et al. 2020; Harikane et al. 2020). Even when correcting for extended [C II] emission (Carniani et al. 2020), the [O III]/[C II] ratios can reach as high as 8, which represents the upper tail of the distribution of low-redshift metal-poor galaxies (Madden et al. 2013). Similar to the [C II] deficits observed in high-redshift galaxies, the exact origin of the high [O III]/[C II] ratios is still debated as the ionization parameter, photodissociation region (PDR) covering fraction, enrichment patterns, metallicity, ISM density, and many other physical properties of a galaxy can impact the [O III]/[C II] ratio (Harikane et al. 2020).

Theoretical models, either analytic or numerical, are an invaluable tool for interpreting the origin of emission lines and physical properties of high-redshift galaxies. While historically, the understanding of emission lines, emission line ratios, and various strong-line diagnostics has relied upon idealized, 1D photoionization models (e.g. Kewley & Dopita 2002; Byler et al. 2017), using codes such as CLOUDY (Ferland et al. 2017) or MAPPINGS (Sutherland et al. 2018), these models often lack some of the complex physical processes (e.g. SN feedback, magnetic fields, shocks, density inhomogeneities, etc.) that affect the physics of the ISM. Although they are significantly more computationally expensive and it is difficult to run 'controlled' experiments as one can do with 1D photoionization models, cosmological simulations provide the best insight into the complex physical nature of the high-redshift ISM.

Due to their computational expense, most recent studies of IR emission lines in cosmological simulations have focused on small samples or individual galaxies (e.g. Vallini et al. 2015; Katz et al. 2017, 2019b; Olsen et al. 2017; Pallottini et al. 2017, 2019a; Lupi et al. 2019; Arata et al. 2020), while only a few have attempted to model the emission in larger cosmological volumes (e.g. Moriwaki et al. 2018; Leung et al. 2020). Obtaining larger samples of galaxies from simulations often relies on semi-analytic models (e.g. Popping et al. 2016; Lagache 2018; Popping et al. 2019). It is not clear whether the results from small sample sizes can be generalized to the larger galaxy population while semi-analytic models do not resolve much of the physics that dictates ISM properties. Furthermore, previous attempts at large volumes containing many galaxies often compromise on important physics such as radiation transfer or sacrifice resolution due to computational expense.

Not all results from simulations agree on the expected trends between [O III] and [C II] emission and various galaxy properties, likely due to differences in spatial and mass resolution, included physics, subgrid models, and the method for modelling emission lines. For example, the [O III]–SFR relations found in Moriwaki et al. (2018), Arata et al. (2020) point towards larger [O III] luminosities per unit SFR compared to Olsen et al. (2017), Katz et al. (2019b). Likewise, the few galaxies simulated in Katz et al. (2019b) seem to fall closer to the local [C II]–SFR relation than those in Pallottini et al. (2017, 2019a). The origin of these discrepancies can in part be due to galaxy-to-galaxy variance, in addition to differences in modelling techniques, which are well established to be far from straightforward (Olsen et al. 2018; Lupi et al. 2020). Nevertheless, there have been many situations where simulations have agreed on emission-line properties of high-redshift galaxies and provide deep insights into or predictions of the early Universe. For example, the necessary but insufficient condition that LyC leakers have high O32 (Barrow et al. 2020; Katz et al. 2020b) as confirmed at $z \sim 2–3$ (Nakajima et al. 2020) or the strong [O I]$_{63 \mu m}$ (hereafter [O I]) emission line (Olsen et al. 2017; Katz et al. 2019b), that was recently observed at $z > 6$ (Rybak et al. 2020).

While systematic uncertainties due to differences in modelling will continue to persist, issues due to galaxy-to-galaxy variance can be solved by simulating larger cosmological volumes. Until now there has been no full-box cosmological radiation hydrodynamics simulation with high enough resolution to semi-resolve the detailed structure of the ISM where these IR emission lines originate and with a large enough computational volume so that both the statistical properties of high-redshift galaxies can be studied and massive enough galaxies are probed to compare directly with observations. The primary goal of this work is to simulate the [C II] and [O III] emission from such a state-of-the-art simulation to better understand the nature of high [O III]/[C II] ratios at high redshift.

This paper is organized as follows: In Section 2, we briefly describe the SPHINX²⁰ simulation (Rosdahl et al., in preparation) as well as our modelling technique for line emission. In Section 3, we present the general [O III] and [C II] properties of high-redshift galaxies at $z > 4.68$ and discuss the physics that drives high [O III]/[C II] ratios. Finally, in Section 5, we present our discussion and conclusions.

# 2 NUMERICAL METHODS

## 2.1 Cosmological simulations

This work makes use of the SPHINX²⁰ cosmological radiation hydrodynamics simulation (Rosdahl et al., in preparation). SPHINX²⁰ is the largest of all simulations in the SPHINX suite of cosmological radiation
(Rosdahl et al. 2018; Katz et al. 2020a) and magnetoradiation (Katz et al. 2021) hydrodynamics simulations and is designed to model both the inhomogeneous process of reionization and the escape of LyC photons from the ISM. SPHINX²⁰ is ideal for studying the emission line properties of high-redshift galaxies because for more than 20 000 galaxies at $z \ge 4.68$, spanning a wide range in halo mass, stellar mass, and star-formation rate (SFR) in each simulation snapshot, we can probe the detailed properties of the ISM. This includes the inhomogeneous radiation field, as well as the detailed density structure, temperature, and metal enrichment, which represent key underlying quantities that impact emission line luminosities.

The details of SPHINX²⁰ are very similar to those presented in Rosdahl et al. (2018). We briefly describe the simulations with emphasis on the differences between the simulations used in this work compared to that of Rosdahl et al. (2018). The simulations are run with RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), which is a radiation hydrodynamics version of the open-source, adaptive-mesh refinement RAMSES code (Teyssier 2002). Initial conditions are generated with MUSIC (Hahn & Abel 2011) in a cosmological box with a comoving side length of 20 Mpc assuming the following cosmology: $\Omega_{\Lambda} = 0.6825$, $\Omega_{m} = 0.3175$, $\Omega_{b} = 0.049$, $h = 0.6711$, and $\sigma_{8} = 0.83$. Due to the small volume of the simulations, the initial conditions were selected from multiple realizations to produce an average DM halo mass function at $z = 6$ (see Rosdahl et al. 2018). The initial composition of the gas is set to be 76 percent H and 24 percent He by mass and have an initial metallicity of $3.2 \times 10^{-4}Z_{\odot}$. The dark matter particle mass is $2.5 \times 10^{5}M_{\odot}$ and we have decreased the stellar particle mass from $1000M_{\odot}$ in SPHINX¹⁰ to $400M_{\odot}$ in SPHINX²⁰. The motivation for this change is to increase the stellar mass resolution to better capture star formation in lower mass galaxies. We can confirm that the stellar content of high mass galaxies is not significantly changed. The AMR grid is adaptively refined to resolve the local Jeans length by at least 4 cells or if the gas or DM mass of a cell is more than 8 times its initial value, resulting in a maximum physical resolution of $7.3 \text{ pc h}^{-1}$ at $z = 6$.

Stars are allowed to form in the simulation at the maximum level of refinement following a thermo-turbulent prescription (Kimm et al. 2017; Trebitsch et al. 2017; Rosdahl et al. 2018). Stellar feedback is included in the form of supernova (SN) explosions following a Kroupa stellar IMF (Kroupa 2001) and Lyman-continuum radiation. Star particles with $3 \text{ Myr} < \text{age} < 50 \text{ Myr}$ can explode via individual SNe by sampling from a delay-time distribution. When an SN event occurs, momentum is injected into the host and surrounding cells such that the amount of momentum injected is dependent on the phase of the SN that is resolved (Kimm et al. 2015). Star particles also continuously inject ionizing photons into their host cells depending on their mass, age, and metallicity. In contrast to SPHINX¹⁰, we only include two radiation energy bins to reduce the memory requirements of the simulation, set at the ionization potentials of H I and He I. The spectral energy distribution (SED) model employed in the simulation is BPASS (Eldridge, Izzard & Tout 2008; Stanway, Eldridge & Becker 2016). We have updated the SED from SPHINX¹⁰ which uses v2.0 to v2.2.1 (Stanway & Eldridge 2018) and adopt a stellar IMF with a maximum mass of $100M_{\odot}$, a minimum mass of $0.1M_{\odot}$, that exhibits a broken power-law shape with a high mass slope of $-2.35$, a low mass slope of $-1.3$, and a break-mass of $0.5M_{\odot}$ (very similar to a Kroupa IMF). Non-equilibrium chemistry is evolved locally for H, $H^{+}$, $e^{-}$, He, $He^{+}$, and $He^{++}$ as are the associated cooling processes (see the Appendix of Rosdahl et al. 2013). Cooling from metal lines is included at $T > 10^{4} \text{ K}$ following tables computed with CLOUDY

(Ferland 1996) as well as at $T \le 10^4 \text{ K}$ using the analytic model from Rosen & Bregman (1995).

For all simulation snapshots, we compute halo properties using the ADAPTAHOP halo finder (Aubert, Pichon & Colombi 2004; Tweed et al. 2009) in the most massive submaxima mode as described in Rosdahl et al. (2018). Furthermore, LyC escape fractions are also computed for each halo following Rosdahl et al. (2018), using the open-source RASCAS code (Michel-Dansac et al. 2020).

## 2.2 Emission-line modelling

Because our goal is to measure the [C II] and [O III] emission line properties of high-redshift galaxies, we must compute the [C II] and [O III] emission lines for each cell in the simulation. It is currently computationally intractable with our available resources to follow the full non-equilibrium chemistry and radiation coupling for C and O directly in the simulation so we resort to post-processing the simulations in order to measure the [C II] and [O III] luminosities. We select seven simulation snapshots at $z = 10, 9, 8, 7, 6, 5, \& 4.68$, the final redshift of our simulation for post-processing. Our method is similar in spirit to that used in Katz et al. (2019b, a, 2020b).

For each cell in the simulation, we know the density, temperature, metallicity, and ionization parameter (the ratio of ionizing photons to hydrogen atoms), as well as the cell size and redshift. For computational efficiency, the on-the-fly radiative transfer in the SPHINX simulations only follows H and He ionizing photons (i.e. those with energies >13.6 eV). To obtain the radiation field at energies below 13.6 eV, we have post-processed the simulation snapshots with radiation transfer in the optically thin limit. This is done by restarting the simulation assuming that the gas, dark matter, and stars are static and using RAMSES-RT to evolve the radiation in two additional energy bins (5.6 eV–11.2 eV and 11.2 eV–13.6 eV) for at least the light-crossing time<fnref n="1" /> of the largest halo in the volume at that redshift. We then use CLOUDY (Ferland et al. 2017) to compute the [O III] and [C II] emission lines from the cells based on their density, temperature, metallicity, redshift, and radiation field in the four<fnref n="2" /> energy bins, i.e. 8 local parameters.

Because there are billions of cells across the seven snapshots, it is impractical to compute a CLOUDY model for each cell. Furthermore, computing a grid across our eight parameters would mean that some parameters are only sparsely sampled in the grid which could result in errors in interpolation. Hence we adopt a different approach as follows:

(i) For each snapshot we randomly sample between 1.5 and 4.5 million simulation cells from galaxies across the full halo mass range.

(ii) A CLOUDY model is run for the selected cells assuming an open geometry gas slab that has the same density, temperature, and metallicity as the simulation cell. The width of the slab is set to be the cell length and the shape of the radiation field within each radiation bin is set to follow that of the Milky Way (i.e. the total energy within the wavelength range of the radiation bin matches that of the simulation but the shape within the bin is kept constant for all cells). The CLOUDY models are all run in a constant temperature mode and a CMB background at the relevant redshift is also included.

<footnote n="1">¹Note that because the simulations use the variable speed of light approximation, we use a conservative estimate for the light crossing time by using the lowest value of the speed of light used in the simulation.</footnote>

<footnote n="2">²This includes the two bins computed in post-processing and the two ionizing radiation bins computed on-the-fly in the simulation.</footnote>
(iii) After the photoionization models have been computed, for each snapshot, we randomly split the results into a training, cross-validation, and test set such that 25 per cent of cells are reserved for testing, 22.5 per cent are used for cross-validation, and 52.5 per cent of cells are used for training. For each line, we train two machine learning models. The first is a classifier that is used to predict whether the emission from the cell is >10⁻³L⊙ while the latter is used to predict the log of the total luminosity of the cell. In contrast to Katz et al. (2019b), we predict the log of the luminosity rather than luminosity to better the accuracy in the low luminosity regions which is important for extended emission. For this work, we use the CATBOOST<fnref n="3" /> (Prokhorenkova et al. 2018; Veronika Dorogush, Ershov & Gulin 2018) implementation of gradient boosted decision trees, which we find to be an improvement in terms of accuracy and well as training and prediction speed compared to the Random Forests used in Katz et al. (2019b). Models are trained for a maximum of 10 000 iterations and are stopped early if the loss-function does not improve for a period of 500 consecutive epochs. For the classifier, we use a binary cross-entropy loss function while for the regressor, we aim to maximize the R² score. For [C II], we find on the test set that the mean accuracy of the classifier is 99.75 ± 0.08 per cent, the mean R² score of the regressor is 0.996 ± 0.001, and that the total luminosity is conserved to within 2.89 per cent ± 2.84 per cent. For [O III], we find on the test set that the mean accuracy of the classifier is 99.91 ± 0.06 per cent, the mean R² score of the regressor is 0.997 ± 0.002, and that the total luminosity is conserved to within 11.7 per cent ± 6.53 per cent. Our values of total luminosity conversion are well within the systematic errors of modelling the emission lines.

(iv) Finally, we apply the models to every cell in each halo across all seven snapshots.

To run the CLOUDY models, we must also set the abundances of each metal. For our fiducial model, we use the GASS command that sets the abundance radios to be consistent with those of the Sun (Grevesse et al. 2010). Such abundances may not be representative of the early Universe due to different chemical enrichment scenarios (e.g. Maiolino & Mannucci 2019). For example, high-redshift galaxies may be predominantly enriched by core-collapse SN compared to a mix of core-collapse SN, type-Ia SN, and winds from AGB stars. For this reason, we also run an additional set of models where the chemical abundance ratios are consistent with the low metallicity (i.e. 1/20Z⊙)<fnref n="4" /> SN yields of Nomoto et al. (2006).<fnref n="5" /> This is particularly important for the [O III]/[C II] ratio because the C/O abundance ratio from Nomoto et al. (2006) is eight times less than that of Solar abundance patterns. Note that the SED for all CLOUDY models is kept fixed for all three metal abundance models. This is further discussed in Section 4.3. Finally, we also test models including molecules but find that our results are not substantially different when these are included and so for the results presented in this work, molecules have been excluded.

Our emission line modelling is computed for every halo in each snapshot with a halo mass >10⁸ M⊙. There are an average of 20 265 haloes of this mass across the seven snapshots.

## 3 RESULTS

In this section, we present our results on [C II] and [O III] emission from high-redshift galaxies, with a particular focus on the origin of high [O III]/[C II] ratios at z > 6 and how these two emission lines can be used to constrain the properties of high-redshift galaxies.

We have calculated the [C II] and [O III] emission for thousands of galaxies in each simulation snapshot between z = 10 and z = 4.68. In Fig. 1 we show density, [C II], and [O III] maps of the three most massive haloes at z = 6, with halo masses of 1.31 × 10¹¹M⊙, 1.29 × 10¹¹M⊙, and 9.97 × 10¹⁰M⊙, respectively. These galaxies have 10 Myr averaged star formation rates (SFRs) of 17.6 M⊙yr⁻¹, 20.7 M⊙yr⁻¹, and 39.7 M⊙yr⁻¹, well within the range of galaxies that have been observed at z > 6 (e.g. Maiolino et al. 2015).

These galaxies exhibit very different morphologies. The first is a turbulent, and partially disrupted rotating disc, the second is a triple major merger, and the third, while more similar to the first, possesses two strongly star forming regions and clear evidence of dynamical disturbance. Such properties are rather characteristic of high-redshift galaxies where strong dynamical interactions are more common than in lower redshift galaxies.

Consistent with our earlier work (Katz et al. 2019b), we find that the [C II] emission follows the densest gas structures in the galaxy while [O III] exhibits a different morphology, more coincident with the young star forming regions and mildly lower density gas (see also Pallottini et al. 2019b; Lupi et al. 2020). With a higher ionization potential and lower critical density compared to C⁺, such behaviour is expected for O⁺⁺.

### 3.1 General properties

With such a large sample of galaxies, we can begin to overcome issues related to sample variance and make predictions regarding the relations between [C II] and [O III] luminosity and galaxy properties.

#### 3.1.1 [C II] – Mvir and [O III] – Mvir relations

In Fig. 2 we show the relation between [C II] and [O III] luminosity and halo virial mass as a function of redshift for the models with Solar abundance ratios and core collapse SN abundance ratios. As halo mass increases, so do the [C II] and [O III] luminosities. Our predicted relation between L_{[C II]} and Mvir using Solar abundance ratios is in good agreement with the semi-analytic model of Silva et al. (2015) at Mvir ≳ 9.5M⊙. The two models deviate from each other at the low mass end. Our models with the same abundance ratios exhibit almost no redshift dependence. Similarly, the relation between L_{[O III]} and Mvir shows no evidence of redshift dependence.

Constraining the [C II]–Mvir and [O III]–Mvir relations is of particular importance for forecasting intensity mapping results (e.g. Lagache 2018; Yue & Ferrara 2019) and we have tables of these values available online.<fnref n="6" />

<footnote n="3">³https://github.com/catboost/catboost</footnote>

<footnote n="4">⁴Note that this value is higher than most of the galaxy metallicities in the simulation.</footnote>

<footnote n="5">⁵The stellar IMF used in Nomoto et al. (2006) has an upper mass threshold of 50 M⊙ rather than 100 M⊙ and a slope of -2.35 across the entire mass range of 0.07M⊙-50M⊙, dissimilar to the IMFs used for stellar feedback in the simulation. Because of this difference, we have opted to keep the total metallicity fixed between the models using Solar and Nomoto et al. (2006) abundances patterns and only change the relative proportion of individual metals.</footnote>

<footnote n="6">⁶https://github.com/HarleyKatz/SPHINX-CII-OIII</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Hydrogen column density ($N_H$), (left-hand panel), [C II] surface brightness (centre), and [O III] surface brightness (right-hand panel) maps of the three most massive galaxies in SPHINX$^{20}$. These haloes have masses of $1.31 \times 10^{11}$ M$_\odot$, $1.29 \times 10^{11}$ M$_\odot$, and $9.97 \times 10^{10}$ M$_\odot$, and 10 Myr averaged star formation rates (SFRs) of $17.6$ M$_\odot$yr$^{-1}$, $20.7$ M$_\odot$yr$^{-1}$, and $39.7$ M$_\odot$yr$^{-1}$, respectively.</figcaption></figure>

### 3.1.2 [C II] and [O III] luminosity functions

Since $M_{vir}$ is not an observed quantity, we do not have any direct observational constraints between [O III], [C II], and $M_{vir}$. However, observational constraints are available for the [C II] luminosity function and in the top panel of Fig. 3, we compare our predicted [C II] luminosity functions as a function of redshift to various observational constraints (Swinbank et al. 2012; Capak et al. 2015; Aravena et al. 2016; Cooke et al. 2018; Uzgil et al. 2021) as well as

other semi-analytic models (Popping et al. 2016; Lagache 2018). We observe a strong redshift dependence of the [C II] luminosity function such that as redshift decreases, the number of [C II] emitters at all luminosities increases. Our simulations have very few galaxies with [C II] luminosities $>10^{42}$ erg s$^{-1}$ which is approximately the lower limit of where observational constraints begin (e.g. at $z = 4.68$ there are five such galaxies while there are none at $z = 6$). Nevertheless, we find that the $z \sim 5.5$ constraint from Capak et al. (2015) falls in
<figure><img src="image_3.png" /><figcaption>**Figure 2.** [C II]–$M_{vir}$ (top) and [O III]–$M_{vir}$ (bottom) relations for galaxies in SPHINX$^{20}$ at $z = 4.68–10$. The solid and dashed lines represent the results for the Solar abundance and core collapse SN abundance models, respectively. The shaded regions represent the $1\sigma$ scatter about the relation. The dotted lines in the top panel represent the semi-analytic models from Silva et al. (2015). Little evolution is seen in these relations as a function of redshift whereas systematic differences occur due to choice of abundance model.</figcaption></figure>

between our $z = 5$ and $z = 6$ predictions. Our models predict far fewer [C II] emitters than the estimates from Aravena et al. (2016) at $z = 6–8$, and are more in line with the $z = 6–8$ upper limit from Uzgil et al. (2021). Compared to results from commonly used semi-analytic models, we predict significantly more [C II] emitters at $L_{[CII]} > 10^8 L_{\odot}$ than both Popping et al. (2016), Lagache (2018). The model that uses core collapse SN abundance ratios predicts a luminosity function much lower than observations.

Unfortunately, there are no observational constraints on the [O III] luminosity function at high redshift and it is not often discussed in the literature (although see Moriwaki et al. 2018). Once again, we expect less redshift evolution in the [O III] luminosity function compared to [C II] due to the fact that O is predominantly produced in core collapse SN; hence, at very high redshifts (e.g. $z \gtrsim 9$), [O III] may be the most likely FIR line to be observed.

### 3.1.3 [C II]–SFR and [O III]–SFR relations

There has been a considerable amount of debate in the literature, both from an observational and theoretical perspective on how the high-

<figure><img src="image_4.png" /><figcaption>**Figure 3.** [C II] (top) and [O III] (bottom) luminosity functions at different redshifts for galaxies in SPHINX$^{20}$ versus observations. The solid and dashed lines represent the results for Solar abundance and core-collapse SN abundance models, respectively. The shaded region represents the $1\sigma$ scatter about the relation. Observational constraints for the [C II] luminosity function are shown for Swinbank et al. (2012), Capak et al. (2015), Aravena et al. (2016), Cooke et al. (2018), Uzgil et al. (2021), Loiacono et al. (2021), and the dotted and dashed lines represent results from the semi-analytic models of Lagache (2018) and Popping et al. (2016), respectively. No observational constraints are available for the [O III] luminosity function at high redshift.</figcaption></figure>

redshift [C II]–SFR relation compares to that in the local Universe. Early observations suggested a large deficit of [C II] emission in $z > 6$ galaxies compared to those at $z = 0$ (e.g. Maiolino et al. 2015; Pentericci et al. 2016; Laporte et al. 2019). However, more recent observations that account for extended [C II] emission reduce some of this difference (Carniani et al. 2020) and by $z \sim 4.5$, the observed relation is nearly identical to that in the local Universe (Schaerer et al. 2020).

The origin of [C II] deficits is currently unknown as simulations also do not agree. For example, the zoom-in simulations of individual high-redshift galaxies in Katz et al. (2019b) show that haloes tend to fall very close to the observed low-redshift relation, and when [C II] deficits do occur in their simulations, it is after a burst of star formation, when the galaxy is also at its brightest, perhaps suggesting an observational bias. The simulations of Arata et al. (2020) produce a [C II]–SFR relation steeper than that observed in the local Universe, likely due to different enrichment channels as well as different ISM properties, although their highly star-forming galaxies are consistent
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mean [C II]–SFR relations (top) and [O III]–SFR relations (bottom) for galaxies in the SPHINX$^{20}$ simulation at $z = 5–9$ compared to observations. The left, centre, and right rows represent models where we assume Solar abundances, core-collapse SN abundances, or core-collapse SN abundances with four times enhanced metallicity, respectively. The shaded regions represent the $1\sigma$ scatter about the relation. The black data points with errors represent observed galaxies at $z = 5–9$ whereas coloured data points represent the brightest [C II] and [O III] emitters at each redshift. We show high-redshift fits from Schaerer et al. (2020), Harikane et al. (2020), Carniani et al. (2020) as well as $z = 0$ relations for dwarfs and starbursts (SBs) from De Looze et al. (2014).</figcaption></figure>

with $z = 0$ observations. In contrast, Pallottini et al. (2017, 2019a) find that for similar halo masses to Katz et al. (2019b), their haloes have both higher star formation rates and significantly larger [C II] deficits.

This disagreement is also present between different simulations regarding the high-redshift [O III]–SFR relation where Katz et al. (2019b), Olsen et al. (2017) find [O III] deficits while Arata et al. (2020), Moriwaki et al. (2018), find an [O III]–SFR relation consistent with $z = 0$. Interestingly, while Arata et al. (2020) cannot produce lower luminosity [O III] emitters, Katz et al. (2019b) struggles to produce high [O III] emitters. In summary, the current lack of agreement in the literature, both in observations and simulations, prohibits robust interpretations of observed high-redshift galaxies and deserves further study.

In the top row of Fig. 4, we show the [C II]–SFR relations for SPHINX$^{20}$ galaxies at different redshifts compared with observations. We show models using Solar abundances, core-collapse SN abundances, and core-collapse SN abundances where we have artificially and uniformly enhanced the metallicities of the gas cells in each simulated galaxy by a factor of 4. Regardless of the abundance model, we see a strong evolution with redshift in that at fixed SFR, the [C II] luminosity increases with decreasing redshift. To first order, the line luminosity scales with the absolute abundance of the metal which is why the Solar abundance model produces higher [C II] for fixed

SFR compared to the core-collapse abundance model. In Fig. 5, we plot the mean ISM metallicity as a function of SFR for various redshifts and see that as redshift decreases, the ISM metallicity at fixed SFR increases. Hence we would expect an increase in [C II] luminosity. The [C II] luminosity is also sensitive to the amount of low temperature neutral gas. We expect more of this gas is present in higher halo masses. In Fig. 6 we show the evolution of the $M_{vir}$–SFR relation as a function of redshift and we also see here that at fixed SFR, halo mass increases with decreasing redshift. Hence this effect should also drive higher [C II] luminosities both due to the increased amount of gas and the decreased susceptibility of the galaxy to SN feedback.

Since [O III] is related to the efficiency of the feedback in the galaxy, we would expect the opposite effect compared to [C II]. The larger halo mass at fixed SFR would reduce the amount of [O III] emission and counteract the increase in metallicity with decreasing redshift. The bottom row of Fig. 4 shows the [O III]–SFR relations as a function of redshift compared with observations. Regardless of the choice of abundance ratios, there is almost no redshift evolution. While it is expected that this relation does not evolve as strongly with redshift as the [C II]–SFR relation, it is an interesting coincidence that there is no redshift evolution at all.

The redshift evolution in the [C II]–SFR relation may be even more pronounced than shown in Fig. 4 due to the expected evolution in
<figure><img src="image_6.png" /><figcaption>**Figure 5.** SFR–Z relation as a function of redshift. Z represents the mass-weighted gas-phase metallicity. The shaded region represents the 1σ standard deviation about the relation. As redshift decreases, galaxies at fixed SFR have higher metallicities.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** SFR–$M_{vir}$ relation as a function of redshift. The shaded region represents the 1σ standard deviation about the relation. As redshift decreases, galaxies at fixed SFR have higher halo masses or alternatively, galaxies at fixed halo mass have lower SFR.</figcaption></figure>

C abundance. While O is predominantly formed in core-collapse SNe soon after the onset of star formation, the formation time-scale for C is significantly longer due to the fact that it is expelled into galaxies via both core-collapse SNe and type-Ia SNe as well as in winds from asymptotic giant branch (AGB) stars (e.g. Maiolino & Mannucci 2019). Hence the amount of C per unit metallicity is expected to increase with decreasing redshift and therefore enhance the redshift evolution in the relation. The top left and top centre panels of Fig. 4 show the extreme cases of Solar abundances and core-collapse SN abundances, respectively. Since we do not explicitly model the varying time-scales of enrichment of different elements in the simulation, we have modulated the abundance ratios in post-processing to explore how this physics may affect the [C II]–SFR relation. Interestingly, when we assume Solar abundance ratios, we find that the simulated galaxies are in much better agreement with observations than when we assume core-collapse SN abundance ratios. At $z \ge 6$, our Solar abundance relations are very consistent

with the fitted relation from Harikane et al. (2020), and by $z = 4.68$, our relation has converged to that of Schaerer et al. (2015). At SFR < $1 M_{\odot}yr^{-1}$, our relations seem to converge more towards those of Carniani et al. (2020) than that of Harikane et al. (2020), although at higher SFR, these two relations are relatively indistinguishable.

Compared to the local relation (shown as the yellow and brown lines in Fig. 4), our $z \ge 6$ galaxies tend to produce small [C II] deficits even for Solar abundance ratios possibly due to the lower metallicity and lower halo mass per fixed SFR. The most [C II] bright galaxies in SPHINX²⁰ are however consistent with the local relation when assuming Solar abundance ratios. This is consistent with the results found in Katz et al. (2019b) which used similar subgrid physics and spatial/mass resolution and demonstrates the need for larger samples of galaxies.

When assuming core-collapse SN abundance ratios, we find very strong [C II] deficits that appear inconsistent with the high-redshift galaxy population due to the drastically reduced C abundance. Taken at face-value, it would appear that our simulations prefer galaxies to exhibit abundance ratios closer to Solar values. However, we argue below that this is not necessarily the case.

One can see in the bottom row of Fig. 4 that regardless whether we assume Solar abundances or core-collapse SN abundances, our [O III]–SFR relation underpredicts individual high-redshift observations as well as the fitted relation from Harikane et al. (2020). There is little difference between the [O III]–SFR relation for the two abundance assumptions because the primary O enrichment channel is via core-collapse SN and thus the O abundance is very similar between the two models. Hence while a Solar abundance ratios model can reproduce the [C II]–SFR relation at high redshift, it struggles with the [O III]–SFR relation, as was the case in our earlier work (Katz et al. 2019b). How then can we reconcile both relations with observations?

Changing the ISM properties of our high redshift galaxies is unlikely to relieve the tension. Because [O III] and [C II] originate in different regions of temperature–density phase-space, attempting to create more [O III]–emitting gas would reduce the amount of [C II] emitting gas. Hence for the Solar abundance model, the reduction in [C II] would disrupt the current agreement in favour of improving the [O III]–SFR relation. For this reason, this is not a preferred solution. If the galaxies were significantly gas richer at fixed halo mass in reality compared to what is predicted by the simulation, then this solution could be an option. However, Katz et al. (2020b) showed that for the SPHINX suite of simulations, these high-mass galaxies accrete gas relatively efficiently. The gas mass could only increase by a factor of about two as to not supersede the cosmic baryon fraction, which is not nearly enough to overcome the near order of magnitude deficit in the [O III]–SFR relation.

Since when assuming Solar abundance ratios, the [C II]–SFR relation is in good agreement with observations, we argue that it is important that the total carbon mass in galaxies is not significantly different from our fiducial model. One solution to our problem could be to increase the abundances of all elements for the Nomoto et al. (2006) abundance patterns. To explore this option, in the right-hand panel of Fig. 4, we plot the relations where we have re-run CLOUDY models where the total metallicity is increased by a factor of 4. Here we can see that the good agreement between the [C II]–SFR relation is once again recovered while the [O III]–SFR relation is in significantly better agreement with the observed high-redshift galaxies. We emphasize that in this model, the C/O relation is the same as in the core-collapse SN abundance model, but the total amount of metals is increased.
While this solution may seem ad hoc, there is precedent for this in high-redshift galaxies. Simulations show that low metallicity galaxies may produce a more top-heavy stellar initial mass function (IMF) compared to what is observed in the local Universe (e.g. Safranek-Shrader, Milosavljević & Bromm 2014; Chon, Omukai & Schneider 2021). More top-heavy IMFs result in an increased amount of SNe which leads to an increase in metal production as well as more high-energy photons which can form more O⁺⁺. Chon et al. (2021) show that the high-mass stellar IMF slope can become significantly flatter at stellar metallicities in the range 0.10Z⊙ < Z < 0.1Z⊙. While the most massive galaxies in our simulation have gas-phase metallicities slightly less than 0.2Z⊙ at z = 6, we highlight the fact that only the nebular metallicity must be enhanced for our proposed solution to work, not the stellar metallicity. Steidel et al. (2016) showed that in order to explain the SEDs of certain z ~ 2–3 galaxies, the nebular metallicities appear 5 × enhanced compared to stellar metallicities. This is partly driven by the way in which stellar and nebular abundances are measured, for example using different elements (e.g. O and Fe) to measure each. Any true enhancement in the nebular metallicity compared to stellar metallicity will support our argument. Furthermore, the work of Chon et al. (2021) and Safranek-Shrader et al. (2014) assume solar abundance patterns. At low temperatures, their cooling functions are dictated by [C II] emission. If one assumes that metal abundances are more consistent with Nomoto et al. (2006), low temperature cooling is likely suppressed which could lead to a top-heavy IMF at even higher metallicities than reported in their work. We emphasize again that given the uncertainties in both the modelling in this work and in determining the stellar IMF from first principles in simulations, both our factor of 4 boost and the discussion in this section are intended to be qualitative rather than quantitative. In SPHINX²⁰, we already artificially boost the energy from SNe by a factor of 4 in order to match the high-redshift UV luminosity function (Rosdahl et al. 2018). This extra energy is thought to be needed in order to overcome numerical overcooling; however it also may be reflecting the necessity for a more top-heavy IMF. Although the energy from SN in the simulation is artificially boosted, the metal yields are still consistent with a standard (Kroupa 2001) IMF. A more top-heavy IMF may produce an ISM that is similar to what we find in this work already due to the boosted feedback although this neglects the fact that the additional metals may enhance cooling. Therefore, we argue that our results are potentially indicative of a top-heavy IMF in the early Universe.

Arata et al. (2020) have explicitly shown with simulations that as metallicity increases, the O/C ratio drops due to enrichment from AGB stars and type Ia SNe. Thus the high-redshift metal-poor galaxies in their work naturally exhibit high [O III]/[C II] ratios. In contrast to our results where our fiducial model has difficulty producing bright [O III] emitters, their model struggles in predicting lower luminosity [O III] and [C II] emitters at high redshift. However since they only simulate a few haloes, the origin of this discrepancy is inconclusive. Arata et al. (2020) choose a slightly more top-heavy IMF (Chabrier) than used in the fiducial model of this work (Kroupa). Their halo at z = 6 is more massive than any of the haloes in SPHINX and thus a deeper potential well may allow the object to retain more of its metals, although their feedback is more efficient at regulating the stellar mass of the object. Perhaps the largest difference between the two simulations is that Arata et al. (2020) does not include cooling below 10⁴ K which is important for a multiphase ISM and emergent line luminosities. This discrepancy makes it difficult to compare the line emission from the two simulations because the [C II] emission is highly sensitive to gas temperature.

### 3.1.4 [O III]/[C II] ratios

In Fig. 7 we plot mean [O III]/[C II] versus SFR relations for SPHINX²⁰ galaxies compared with observations. We show the Solar abundance model (left-hand panel), the core-collapse SN abundance model (centre), and the core-collapse SN abundance model where the metallicity has been enhanced by a factor of 4 (right-hand panel). Note that the latter two models predict very similar [O III]/[C II] ratios. Thus any future discussion of [O III]/[C II] ratios in the core-collapse SN abundance model also applies to that with boosted metallicity.

In contrast to the observed relations, we find that the trend between [O III]/[C II] and SFR is not monotonic. Between −2 < log₁₀(SFR/M☉yr⁻¹) < 0, we find decreasing [O III]/[C II] ratios as SFR decreases however, at log₁₀(SFR/M☉yr⁻¹) > 0, the trend flattens, regardless of the assumption on abundance ratios. Redshift evolution pushes galaxies towards lower [O III]/[C II] with decreasing redshift, especially at z < 6.

As expected, compared to the observed high-redshift, z = 5–9 galaxies, the Solar abundance model produces [O III]/[C II] ratios that are far lower. At SFRs of 10 M☉yr⁻¹ the [O III]/[C II] deficit for the mean galaxy in our simulations is more than an order of magnitude compared to observations (see also Laporte et al. 2019). The [O III]/[C II] ratios are even lower than the metal-poor dwarf galaxies in the local Universe and more consistent with local starburst galaxies. However, when we assume that the O and C abundances are similar to that of core-collapse SNe, we find that the high [O III]/[C II] ratios at high redshift can be reproduced by our simulation. Still, we do not have many galaxies with [O III]/[C II] > 10 as seen in Harikane et al. (2020). When accounting for [C II] surface brightness dimming corrections (Carniani et al. 2020), we find very good agreement between z = 5–9 galaxies and the highly star-forming galaxies in SPHINX²⁰.

Our results indicate that the C/O abundance ratios are a crucial factor for high-redshift [O III]/[C II] ratios. Such an argument was proposed by Steidel et al. (2016) to explain the offsets of z ~ 2–3 galaxies on the BPT diagram. Similarly, Arata et al. (2020) demonstrated with cosmological simulations that explicitly track C and O abundances from different chemical enrichment channels that the C/O ratio is indeed significantly lower in high-redshift galaxies due to the delayed onset of AGB winds. However, Harikane et al. (2020) suggest that while lower C/O ratios may explain some of the enhanced [O III]/[C II] ratios observed at high redshift, the reduced C abundance alone is not enough to explain some of the non-detections of [C II]; hence, the ratio by itself cannot explain the galaxies with >1 dex lower [C II]/SFR. Once again, these [C II] deficits are partially relieved by surface brightness dimming corrections (Carniani et al. 2020); however, it is important to understand what other physical mechanisms can drive high [O III]/[C II] ratios as well as [C II] deficits.

## 3.2 What else drives high [O III]/[C II] ratios?

Harikane et al. (2020) systematically discuss eight different effects that can impact the [O III]/[C II] ratio in high-redshift galaxies. They employ idealized models of the ISM using CLOUDY to better understand the impact of different ISM properties on the [O III]/[C II] ratio. While these simplified photoionization models are useful for isolating individual effects, real high-redshift galaxies exhibit a complicated, multiphase ISM structure. In this section, we explore seven of the eight effects<fnref n="7" /> in the context of the SPHINX²⁰ simulations.

<footnote n="7">⁷The last effect is the lower C/O abundance ratios due to different enrichment channels in the early Universe which has already been discussed above.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** [O III]/[C II] versus SFR for galaxies in our simulation at different redshifts compared to observations for the Solar abundance model (top) and the core-collapse SN abundance model (bottom). The shaded regions represent the 1σ scatter about the relation. The black data points represent $z = 5–9$ galaxies from Harikane et al. (2020) while red data points represent those same galaxies accounting for surface brightness dimming (Carniani et al. 2020). The orange points represent individual low redshift dwarf galaxies from the DGS survey (Madden et al. 2013). The red dotted line is the high-redshift fitted relation from Harikane et al. (2020) while the maroon and orange lines are the $z = 0$ relations for dwarf and starburst galaxies, respectively, from De Looze et al. (2014).</figcaption></figure>

### 3.2.1 Lower PDR covering fractions

[O III] is generally emitted very close to stars from H II regions while [C II] tends to originate further out in PDRs. Therefore if the H II region becomes large enough so that the gas cloud is density-bounded (i.e. the outer radius of the cloud is smaller than the Stromgren radius), we expect there to be very little [C II] emission. Molecular clouds exhibit complex density structures where some regions around the star are expected to remain shielded from the ionizing radiation and exist in a neutral state. Harikane et al. (2020) argue that if the covering fraction of PDRs ($C_{PDR}$) surrounding stars is significantly less than unity, we would expect an increase in [O III]/[C II].

Because the ISM structure in our simulations is highly complex, we cannot define a direct metric for $C_{PDR}$ in our simulations that matches the idealized models of Harikane et al. (2020). Nevertheless, we have measured the LyC escape fraction ($f_{esc}$) for each galaxy in SPHINX²⁰ which provides a proxy for $C_{PDR}$ on a galaxy scale. We expect that galaxies that have higher $f_{esc}$ have lower PDR covering fractions and should thus exhibit higher ratios of [O III]/[C II].

From simulations, it is well established that $f_{esc}$ is a feedback regulated quantity (e.g. Kimm et al. 2017; Trebitsch et al. 2017). Hence it is the SN feedback that dictates the PDR covering fraction in our simulation more than the ionizing photon production. Such an effect can only be realized in numerical simulations.

In the bottom-centre panel of Fig. 8, we show a 2D histogram of [C II]/SFR versus [O III]/SFR coloured by $log_{10}(f_{esc})$. Galaxies in our simulation that exhibit the highest values of $f_{esc}$ and, by deduction, lower PDR covering fractions also tend to exhibit weaker [C II]/SFR and higher [O III]/SFR, or in other words, higher [O III]/[C II] ratios.

As the Universe evolves and becomes more metal enriched, the typical mass of a galaxy also increases while the luminosity-weighted global LyC escape fractions tend to decrease (e.g. Rosdahl et al. 2018). Certain analytical models predict that $f_{esc}$ evolves very strongly with redshift such that $f_{esc} \propto (1+z)^{3.4}$ (Haardt & Madau 2012). Harikane et al. (2020) show that for $C_{PDR} = 0$, [C II] luminosity is decreased by 99 per cent and therefore, the evolving escape fraction

can contribute to the high [O III]/[C II] ratios seen at high redshift. Furthermore, given the correlation between [O III]/[C II] and $f_{esc}$, Inoue et al. (2016), Katz et al. (2020b) argued that one can select LyC leakers from non-leakers based on this ratio alone. This metric was tentatively tested (Katz et al. 2020b) on the few low-redshift galaxies where observations of both lines exist and a common trend for the systems that were predicted to be leakers was that many appeared to have strong signatures of feedback, consistent with the results from our simulations.

To demonstrate this relation more explicitly, in Fig. 9 we plot the mean value and 1σ scatter of $f_{esc}$ versus $log_{10}([O III]/[C II])$ for each redshift interval for models where we assume either Solar abundance ratios or core-collapse SN abundance ratios. Although there is considerable scatter in the relation, there is a very strong trend that the galaxies that have the highest values of [O III]/[C II] also exhibit the highest $f_{esc}$. Assuming core-collapse SN abundances systematically shifts the relation towards higher values of [O III]/[C II]; nevertheless, the trend persists.

Consistent with Harikane et al. (2020), we find that lower PDR covering fractions or equivalently increased LyC escape fractions can drive galaxies towards higher [O III]/[C II].

### 3.2.2 Higher ISM densities

Galaxies that exhibit the highest ISM densities are expected to exhibit the lowest ratios of [C II]/SFR and [O III]/SFR because once the density surpasses the critical density of the line, the emissivity increases as $n$ rather than $n^2$ due to collisional de-excitation. In the top right-hand panel of Fig. 8 we show a 2D histogram of [C II]/SFR versus [O III]/SFR, coloured by the 90th percentile density for gas that emits<fnref n="8" /> [C II] and [O III]. There is a clear gradient indicating that the galaxies with the lowest values of [C II]/SFR and [O III]/SFR also have the highest ISM densities.

<footnote n="8">⁸This gas is defined as having $T < 10^5$ K and $\rho > 1$ cm$^{-3}$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** 2D histograms of [C II]/SFR versus [O III]/SFR for galaxies in the SPHINX²⁰ simulation between $6 \le z \le 10$ coloured by virial mass (top left-hand panel), sSFR (top centre), 90th percentile ISM density (top right-hand panel), ISM metallicity (bottom left-hand panel), LyC escape fraction (bottom centre), and 90th percentile ISM ionization parameter (bottom right-hand panel). We show results for the core-collapse SN abundance model with the 4× enhanced metallicity, $n_H$, Z, and q have are calculated as the mass-weighted quantities in the [C II] and [O III] emitting regions (Z, $n_H$) or H II regions ($\bar{q}$). The cyan points represent individual high-redshift galaxies from Carniani et al. (2020).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** [O III]/[C II] versus $f_{\text{esc}}$ for galaxies in our simulation versus observations. The coloured lines and shaded regions indicate the mean relation and 1σ scatter for different redshift bins. The solid and dashed lines represent the results when we assume Solar abundances and core-collapse SN abundances, respectively. We only include galaxies in this plot that have an SFR > $10^{-2} \text{M}_\odot \text{yr}^{-1}$.</figcaption></figure>

The critical density of [C II] for interactions with hydrogen is ~2800 cm⁻³ while the critical density of [O III] for interactions with electrons is ~510 cm⁻³ (Carilli & Walter 2013). Hence we expect that galaxies with higher ISM densities will exhibit lower ratios of [O III]/[C II] because as density increases, we reach the critical density of [O III] before [C II]. In Fig. 10, we plot the 90th percentile ISM density of the [O III] and [C II] emitting gas in the simulated galaxies versus [O III]/[C II] for the cases where we assume both Solar abundances as well as core-collapse SN abundances. For lower values of [O III]/[C II], we see a decrease in ISM density with decreasing redshift. Consistent with expectations, galaxies with lower ISM densities exhibit higher ratios of [O III]/[C II].

We do however stress that one should exercise caution when interpreting the trend between ISM density and [O III]/[C II] for our simulated high-redshift galaxies. While such effects can be easily tested in controlled idealized models such as those presented in Harikane et al. (2020), in cosmological simulations, other parameters are not held fixed with the changing ISM density. For example, low density may indicate that the galaxy has been susceptible to strong SN feedback which can destroy the [C II] emitting regions of the galaxy, enhance the amount of $O^{++}$, and increase $f_{\text{esc}}$, all of which are important for setting the [O III]/[C II] ratio. Furthermore, it is well established that [O III] is emitted from different regions of temperature-density phase-space than [C II] and thus the gas that
<figure><img src="image_11.png" /><figcaption>**Figure 10.** [O III]/[C II] versus $log_{10}(\bar{\rho}_{90}/\text{cm}^{-3})$ for galaxies in our simulation versus observations. $\rho$ is defined to be the 90th percentile value from the [C II] and [O III] emitting gas in the galaxies. The coloured lines indicate the median relation for different redshift bins. The solid and dashed lines represent the results when we assume Solar abundances and core-collapse SN abundances, respectively. We only include galaxies in this plot that have an SFR > $10^{-2} \text{ M}_\odot\text{yr}^{-1}$.</figcaption></figure>

emits [O III] may be characterized by a different density compared to that which emits [C II].

To demonstrate this, in Fig. 11 we plot the luminosity-weighted mean gas density for [C II] against that of [O III] coloured by the ratio of [O III]/[C II]. While there is considerable scatter, in general, the galaxy regions that emit strongly in [O III] have significantly lower gas densities compared to the regions with high [C II] luminosity. It is often the case that the densities are more than an order of magnitude lower in the [O III] emitting regions compared to [C II]. This suggests that idealized 1D photoionization models with a fixed ISM density poorly mimic the conditions of [O III] and [C II] emitting regions in the ISM of high-redshift galaxies. Rather, our simulations are more akin<fnref n="9" /> to a constant pressure assumption for the 1D models (see e.g. Cormier et al. 2019) that can, in principle, capture these density contrasts. Since, many of the galaxies with high [O III]/[C II] ratios do not have ISM densities above the critical density for [O III] and [C II], despite the large density contrasts between [C II] and [O III] emitting regions, the higher ISM densities do not play a strong role in setting the [O III]/[C II] ratio.

### 3.2.3 Higher ionization parameters

The photoionization models of Harikane et al. (2020) show that models with higher ionization parameters ($q = cn_{\gamma,\text{ion}}/n_{\text{H}}$) tend to exhibit higher ratios of [O III]/[C II]. In our simulations, [C II] is predominantly emitted from higher density neutral regions whereas [O III] originates in lower density, ionized gas with $T \gtrsim 10^4$ K. Additional high-energy photons can ionize some of the $C^+$ into $C^{++}$ which also serves to reduce [C II] emission while additional

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Luminosity-weighted gas density of [C II] regions versus [O III] regions for galaxies at $z=6$ coloured by the ratio of [O III]/[C II]. We only include galaxies that have SFR > $10^{-2} \text{ M}_\odot\text{yr}^{-1}$ as well as [C II] and [O III] luminosities >$10^3 L_\odot$. [O III]/[C II] ratios are shown for the Solar abundance model although the results are not fundamentally different for the core-collapse SN abundance model.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** [O III]/[C II] versus $log_{10}(\bar{q})$ for galaxies in our simulation. $\bar{q}$ is defined to be the 90th percentile value from the H II regions in each galaxy. The coloured lines indicate the median relation for different redshift bins. The solid and dashed lines represent the results when we assume Solar abundances and core-collapse SN abundances, respectively. We only include galaxies in this plot that have an SFR > $10^{-2} \text{ M}_\odot\text{yr}^{-1}$.</figcaption></figure>

hard ionizing photons have the potential to create more $O^{++}$, thereby increasing the ratio of [O III]/[C II].

While it is trivial to set the ionization parameter in a CLOUDY calculation, real galaxies exhibit huge ranges in $q$ depending on location within the galaxy. In order to describe the galaxy by a single value, we select all H II regions defined as gas with $T < 10^{5.5}$ K, $n_H [> 1 \text{ cm}^{-3}$ and $x_{\text{HII}} > 50$ per cent, and measure the 90th percentile mass-weighted ionization parameter. In the bottom right-hand panel of Fig. 8 we show a 2D histogram of [C II]/SFR versus [O III]/SFR coloured by $log_{10}(\bar{q})$. As predicted by Harikane et al. (2020), we find the gradient in $q$ evolves towards the upper left of the figure (i.e. high ionization parameter results in lower [C II]/SFR and higher [O III]/SFR). This can be seen more clearly in Fig. 12.

<footnote n="9">⁹This statement may appear to conflict with the CLOUDY models that were used to post-process the cosmological simulations in this work which assumed constant density and temperature. However, the ISM in our simulations is semi-resolved so each individual model is more similar to an individual zone in a CLOUDY model that captures H II regions, and PDRs in the same slab.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** $\log_{10}(\bar{q})$ versus metallicity for galaxies in our simulation. $\bar{q}$ is defined to be the 90th percentile value from the H II regions in each galaxy whereas metallicity is the mass-weighted mean gas-phase metallicity in the ISM. We only show galaxies with SFR > $1 M_{\odot} yr^{-1}$ for the core-collapse SN abundance model with enhanced metallicity. Data points with error bars represent ranges in ionization parameter and metallicity for SDSS galaxies, $z \sim 1$ galaxies, LyC leakers, Lyman-break analogues (LBAs), green pea (GP) galaxies, $z = 2-3$ Lyman-break galaxies (LBGs), and $z = 2-3$ Ly$\alpha$ emitters (LAEs) from Nakajima & Ouchi (2014).</figcaption></figure>

where we plot $\log_{10}(\bar{q})$ directly against $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ for both the Solar abundance and core-collapse SN abundance models. Ionization parameters of galaxies with $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratios of 10 can reach as high as $\log_{10}(\bar{q}) \gtrsim 9$. Therefore the higher $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratios seen at high redshift could be partially due to higher ionization parameters.

While measurements of the ionization parameter at $z > 6$ are not yet possible, proxies from lower redshift may give an indication as to whether reionization epoch systems truly exhibit higher values of $q$. The O32 diagnostic (i.e. $[\mathrm{O\,III}]_{4960\AA},_{5007\AA}/[\mathrm{O\,II}]_{3727\AA},_{3729\AA}$) is often used as a proxy for ionization parameter (e.g. Penston et al. 1990). Strom et al. (2017) demonstrated using the KBSS survey that star-forming galaxies at $z \sim 2-3$ exhibit considerably higher O32 compared to low-redshift SDSS galaxies, indicative of higher ionization parameters. Similarly, Katz et al. (2019b) showed that massive high-redshift galaxies in the epoch of reionization have very similar strong-line diagnostics compared to these lower redshift systems. Consistent results were also found for LyC leakers in the redshift range $z \sim 2-4$ (Nakajima et al. 2020).

In Fig. 13, we compare SPHINX$^{20}$ galaxies with SFRs > $1 M_{\odot} yr^{-1}$ to different classes of observed galaxies on the ionization parameter–metallicity plane assuming the core-collapse SN abundance model with enhanced metallicity. The vast majority of star-forming SPHINX$^{20}$ galaxies exhibit ionization parameters greater than that of SDSS galaxies, and often equal to or greater than the ionization parameters in $z = 2-3$ Ly$\alpha$ emitters or known LyC leakers. There is a clear trend that lower metallicity star-forming galaxies exhibit higher ionization parameters. Furthermore, it is often the lower metallicity, higher ionization parameter galaxies that exhibit the highest $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratios. However, at fixed $q$, the $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ seems to be relatively constant (this is explored further in the next section), hence we argue that $q$ is more important for setting the $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratio compared to $Z$.

In summary, observational evidence combined with theoretical models suggests that high-redshift galaxies exhibit higher ionization

parameters which can contribute to the high ratios of $[\mathrm{O\,III}]/[\mathrm{C\,II}]$. However, in the case of the solar abundance model, all of the observed high-redshift galaxies would need to exhibit ionization parameters of the order of $\log(q) \sim 9$ in order to reach the $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ observed. The core-collapse SN models can reach the observed ratios with a lower ionization parameter.

### 3.2.4 Lower gas metallicities

To first order, the line luminosity of each simulation cell scales with the amount of metals present in the gas. Harikane et al. (2020) argue that for lower metallicity gas, one expects a drop in $[\mathrm{O\,III}/\mathrm{SFR}]$ but a constant $[\mathrm{C\,II}/\mathrm{SFR}]$. Since the bulk of the $[\mathrm{C\,II}]$ emission arises from PDRs, the PDR column density is proportional to $1/Z$ if dust dominates the self-shielding from UV photons (Kaufman, Wolfire & Hollenbach 2006). For dwarf galaxies with SFR $\ge 1 M_{\odot} yr^{-1}$ in the local Universe, De Looze et al. (2014) find an enhancement in $[\mathrm{C\,II}/\mathrm{SFR}]$ for metal-poor dwarf galaxies compared to other galaxy types.

Similar to ISM density, the $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratio is sensitive to numerous different effects that either affect or are impacted by metallicity. For example, lower metallicity stars have harder SEDs which could result in the production of more $O^{++}$, while the additional ionizing photons may ionize some of the PDRs that would be present otherwise. The abundance ratios are expected to change with metallicity and Arata et al. (2020) predicted from simulations that lower metallicity galaxies exhibit higher ratios of $[\mathrm{O\,III}]/[\mathrm{CII}]$. The final amount of momentum injected during the SN snowplow phase is weakly sensitive to the metallicity of the gas (Blondin et al. 1998; Thornton et al. 1998; Karpov et al. 2020) and similarly the ability for gas to cool is also sensitive to metallicity (e.g. Sutherland & Dopita 1993). Finally, there is a well established mass-metallicity relation in the local Universe (Tremonti et al. 2004) and total halo mass is partially responsible for setting the maximum densities achievable in our simulation.

In the bottom left-hand panel of Fig. 8, we show a 2D histogram of $[\mathrm{C\,II}/\mathrm{SFR}]$ versus $[\mathrm{O\,III}/\mathrm{SFR}]$ for galaxies in our simulation coloured by average gas-phase metallicity in the regions of the galaxies that emit $[\mathrm{O\,III}]$ and $[\mathrm{C\,II}]$. The higher metallicity galaxies evolve towards the top right of the diagram, consistent with the idea that to first order, $[\mathrm{O\,III}]$ and $[\mathrm{C\,II}]$ luminosity scale with total metallicity. Note that our CLOUDY calculations do include dust, but we have chosen a dust-to-metal ratio that scales with metallicity (Rémy-Ruyer et al. 2014) so most of the galaxies are relatively dust free. In fact, our results are very similar when the CLOUDY models are run without dust. Therefore, we do not recover the trend predicted by Harikane et al. (2020). Furthermore, as the abundance ratios in our simulation do not evolve with metallicity, we also do not recover the trend seen in Arata et al. (2020), although we emphasize that their trend is likely more realistic and the evolution is most likely bounded by our two different abundance ratio models. Rather, metallicity in our simulation, in the context of the $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ ratio, is a more sensitive indicator of the ISM physics and star formation history which helps us set the emission line luminosities.

In Fig. 14, we plot the mean relations of $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ versus $\log_{10}(Z/Z_{\odot})$ as a function of redshift. There is a strong evolution with decreasing redshift towards higher values of metallicity simply because as the simulation evolves, the galaxies become more metal enriched. The trend between $[\mathrm{O\,III}]/[\mathrm{C\,II}]$ and metallicity is essentially flat indicating that for dust-poor galaxies, excluding the impact of evolving abundance ratios with metallicity, in general, the metallicity
<figure><img src="image_15.png" /><figcaption>**Figure 14.** [O III]/[C II] versus log₁₀(Z/Z⊙) for galaxies in our simulation. The coloured lines indicate the mean relation for different redshift bins. The solid and dashed lines represent the results when we assume Solar abundances and core-collapse SN abundances, respectively. We only include galaxies in this plot that have an SFR > 10⁻²M⊙yr⁻¹.</figcaption></figure>

of the galaxy does not strongly impact the [O III]/[C II] ratio. We do however point out that for the galaxies with high SFRs, we indeed find a trend that lower metallicity galaxies tend to also exhibit higher ionization parameters and higher [O III]/[C II] ratios (see Fig. 13). However, at fixed q, the [O III]/[C II] ratio seems constant. While in general, metallicity may not impact the relation for the majority of galaxies in our simulation, at higher SFRs, since there is a trend between ionization parameter and metallicity, it may appear that low metallicity correlates with high [O III]/[C II], opposite of what is suggested in Harikane et al. (2020).

### 3.2.5 Cosmic microwave background attenuation

The CMB can have multiple effects on high-redshift infrared line emission because the radiation temperature is significantly higher at $z > 6$ than it is in the local Universe. The CMB acts as both a background against which one has to observe the line and can also act as a heating term which is important for low temperature gas (e.g. Lagache 2018). While the CMB can impact [O III], it is particularly important in the context of [C II], especially for lower density, low-temperature gas (see e.g. fig. 1 of Kohandel et al. 2019). Strong CMB attenuation of the [C II] line can drive higher ratios of [O III]/[C II] and partially explain the elevated ratio at high-redshift compared to the local Universe.

Both Laporte et al. (2019) and Harikane et al. (2020) considered the impact of CMB attenuation and concluded that the effect was not strong enough to be completely responsible for the elevated [O III]/[C II] ratios. Lagache (2018) predicted that the maximum increase in [O III]/[C II] due to CMB attenuation is ~0.5 dex. This value is sensitive to both the temperature and density structure of the gas that emits [C II].

To measure the impact of CMB attenuation, we have run an additional set of CLOUDY models<fnref n="10" /> for our simulated galaxies at $z = 6$ where we have artificially set the CMB background to be consistent

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Stacked [C II] and [O III] surface brightness profiles at $z = 6$ for the three most massive galaxies in the simulation viewed along three different angles compared to an observed stack at $z = 5-7$ from Fujimoto et al. (2019). The [O III] profile is only slightly more concentrated than that of [C II]. The thick lines represent the profiles that have been smoothed with the ALMA beam while the thin lines represent the intrinsic stacked profile. All profiles have been normalized so that the [C II] profile matches the innermost data point from Fujimoto et al. (2019).</figcaption></figure>

with the value at $z = 0$, thus removing the CMB attenuation effects. Considering only galaxies that have [C II] and [O III] luminosities of at least 10³L⊙, for most of the systems at $z = 6$, the enhancement in [O III]/[C II] is only ~0.05 dex and thus we concur with Laporte et al. (2019) and Harikane et al. (2020) that while CMB attenuation enhances the [O III]/[C II] ratio slightly, it cannot be a dominant factor.

### 3.2.6 Extended [C II] emission

Using stacks of $z = 5-7$ galaxies, Fujimoto et al. (2019) identified large haloes of [C II], extending out to ~10 kpc in galaxies with SFRs of ~10–70 M⊙yr⁻¹. Taking into account this extended emission, Carniani et al. (2020) determined that earlier observations may underestimate the [C II] luminosity by up to a factor of ~2. Simulations have struggled in reproducing large [C II] haloes (e.g. Katz et al. 2019b; Pallottini et al. 2019a; Arata et al. 2020) such as those observed in Fujimoto et al. (2019). Nevertheless, if the [C II] surface brightness profile is significantly more extended than that of [O III], the underestimated values [C II] may lead to artificially high values of [O III]/[C II].

In Fig. 15, we plot the [O III] and [C II] surface brightness profiles for a stack of the three most massive galaxies in our simulation at $z = 6$ for three different viewing angles. These are the galaxies shown in Fig. 1 and represent the most similar galaxies in our simulation to those observed in Fujimoto et al. (2019). We show results for when the stacks are smoothed with the ALMA beam (thick lines) as well as the unsmoothed values (thin lines). The surface brightness of the stacks have been renormalized by a constant factor so that the [C II] surface brightness profiles match the innermost data point from Fujimoto et al. (2019). While the [C II] surface brightness profile is slightly more extended than that of [O III], consistent with observations from Carniani et al. (2020), the difference is small.

The smoothed profile is in good agreement with the data from Fujimoto et al. (2019) up to ~4.5 kpc before dropping below the observational data points. Smoothing moves the intrinsic profile into significantly better agreement with observed data, but there is still a

<footnote n="10">¹⁰These CLOUDY models were run without dust and compared to a matched set of CLOUDY models with a $z = 6$ CMB that also neglect dust.</footnote>
discrepancy at large radii. Part of the extended profile that we see is driven by mergers. In the second row of Fig. 1, we can see that the second most massive halo at $z = 6$ is a triple merger of three smaller galaxies. It is not particularly obvious how to choose a centre for stacking the images given the complex morphology and this can also play a role in the shape of the surface brightness profile. In Fig. 15, we have opted to use the centre of the halo.

Because the [O III] surface brightness profile is so similar to that of [C II], it is unlikely that galaxy [O III] and [C II] morphologies are responsible for the high [O III]/[C II] ratios observed at high redshift, consistent with the arguments posed in Harikane et al. (2020) because even when accounting for surface brightness dimming in observations (Carniani et al. 2020), high-redshift galaxies still exhibit high [O III]/[C II] ratios.

### 3.2.7 Inclination effects

Kohandel et al. (2019) demonstrated using hydrodynamics simulations that because the shape of the [C II] spectrum changes with viewing angle of a galaxy, for orientations that lead to large velocity spreads (e.g. edge-on) the [C II] flux may not be high enough to be detected by ALMA. In contrast, when the galaxy is viewed face-on, it is more likely to be detected. Harikane et al. (2020) speculate that this effect is unlikely to change the [O III]/[C II] ratio because if [C II] is undetected due to orientation, [O III] is also likely to be impacted in a similar way. However, this argument implicitly assumes that [C II] and [O III] have similarly shaped spectra or that [O III] is not always more centrally concentrated. Katz et al. (2019b) showed that in certain cases, the first of these two assumptions is incorrect, especially during strong dynamical interactions or feedback events. This once again arises because [O III] is emitted from a different gas phase than [C II].

In Fig. 16 we show the [C II] and [O III] spectra of the three most massive galaxies at $z = 6$ for three different viewing angles. We have normalized each spectrum so the full-width at half-maximum (FWHM) can be more easily compared between the lines that have different luminosities. Consistent with Kohandel et al. (2019), we also see strong variations in the spectra as a function of viewing angle as expected. While the [O III] and [C II] spectra exhibit many similar features, the shapes are not identical. In some cases, the [C II] spectrum is more extended than that of [O III] (see the solid lines in the top panel) which can lead to a [C II] deficit for this galaxy; however, this is not a generic feature. For example, in the second row, the dotted line peak for [O III] is nearly 200 km s$^{-1}$ offset from the [C II] peak towards negative velocity.

While galaxy orientation may explain [C II] deficits in specific high-redshift galaxies, it cannot explain the high [O III]/[C II] ratios in all of the observed high-redshift systems, unless there is an observational bias that preferentially selects concentrated [O III] emission or spatially extended [C II] emission.

## 3.3 Alternative explanations for high [O III]/[C II] ratios

We have discussed eight different effects that can potentially drive high-redshift galaxies towards high [O III]/[C II] ratios; however, there is still physics that is not considered in our model that may play a role in setting the [O III]/[C II] ratios of high redshift galaxies.

At slightly lower redshifts than considered here, it has been shown that the [C II] line is optically thick in some of the brightest submillimetre galaxies (e.g. Neri et al. 2014). We can obtain a rough estimate of [C II] by measuring the optical depth to [C II] in every

<figure><img src="image_17.png" /><figcaption>**Figure 16.** [C II] (cyan) and [O III] (magenta) spectra for three different galaxies at $z = 6$ for three different viewing angles (solid, dashed, and dotted lines). The shape of the spectra change depending on the viewing angle and it is not necessarily the case that the [O III] and [C II] spectra exhibit the same shape. We show the spectra for the Solar abundance model (although the shape will be similar for the core-collapse SN abundance model) and the spectra have been smoothed over bins of 20 km s$^{-1}$. We have normalized each spectrum to its maximum value so that the FWHM can be easily compared.</figcaption></figure>

simulation cell. Following Goldsmith et al. (2012), Lagache (2018), the optical depth ($\tau$) can be written as

$$ \tau = \tau_0 \frac{1 - e^{-T_* / T_{\text{ex}}}}{1 + \left( \frac{g_u}{g_l} \right) e^{-T_* / T_{\text{ex}}}}, \quad (1) $$

where $T_* = 91.25$ K is the equivalent temperature of the transition, $g_u = 4$ and $g_l = 2$ are the statistical weights of the upper and lower levels, $T_{\text{ex}}$ is the excitation temperature of the gas, and

$$ \tau_0 = 7.49 \times 10^{-18} N_{C^+} / \delta v, \quad (2) $$

where $N_{C^+}$ is the C II column density and $\delta v$ is the line width. The line width is estimated based on the temperature of the cell and the microturbulence of the gas such that $\delta v = \sqrt{a_{\text{therm}}^2 + a_{\text{turb}}^2}$. Here, $a_{\text{therm}} = \sqrt{\frac{2k_B T}{\mu_B m_p}}$, where T is the temperature of the cell, $k_B$ is the Boltzmann constant, $m_p$ is the proton mass, and $\mu$ is the mean weight of $^{12}$C$^+$. To estimate the microturbulence, we follow Larson (1981), Franeck et al. (2018), and set $a_{\text{turb}} = 1.1 \times (l/\text{pc})^{0.38}$, where l is the physical length of a cell in the simulation. Based on our CLOUDY models, the fraction of C that is in the form of C$^+$ is well described as a function of temperature such that

$$ f_{C^+} = \frac{1}{1 + e^{-(147.27 - 35.69 \log_{10}(T/K))}}, \quad (3) $$

$N_{C^+}$ can then be calculated for every cell knowing the density, metallicity, $f_{C^+}$, and cell size. We can measure the excitation temperature in each cell using the collisional de-excitation rate coefficients for neutral hydrogen and electrons as given in Goldsmith et al. (2012) and assuming a CMB background at $z = 6$. To calculate the true extinction, one must then ray trace through the simulation accounting also for the bulk velocities of every gas cell. Rather, to obtain a simple estimate for the effect of self-absorption, we assume that the emission can only be absorbed locally and that nearby cells exhibit similar $\tau$.
Under this simple approximation, for the most massive galaxy at $z = 6$, we find that $L_{[CII]}$ is reduced from $5.24 \times 10^{41}$ ergs s$^{-1}$ to $2.62 \times 10^{41}$ ergs s$^{-1}$, which is equal to a 50 per cent decrease in total luminosity. We have completed this exercise for the next two most massive galaxies at $z = 6$ and find similar reduction factors. Hence self-absorption can indeed be important for the most massive objects in the epoch of reionization; however, it is unlikely to be fully responsible for the high [O III]/[C II] ratios seen, especially in the lower mass objects that are less metal enriched and when turbulence is high compared to the thermal velocity.

Another physical effect not considered here is different rates of depletion of C and O on to dust. While C and O are not expected to be heavily depleted like some other elements (e.g. Fe; Jenkins 2009), if for example C is more depleted than O, we would expect that the [O III]/[C II] ratio would be enhanced. Because our simulated galaxies have low metallicity, we expect the effect of depletion of metals on to dust to be small.

We have hinted earlier that we also cannot rule out observational bias and the preferential selection of galaxies at high-redshift that either exhibit [C II] deficits or enhanced [O III] emission. While observational bias does not explain why some high-redshift galaxies exhibit exceptionally high [O III]/[C II] ratios, it may help us explain why all of the nine observed high-redshift galaxies with both [O III] and [C II] detections fall high compared to the local relation. Galaxies are often selected by means of UV dropout. Since UV photons are primarily emitted by young stars, at fixed halo mass, a galaxy is more likely than not to be observed during a burst of star formation. Since both ionizing photons from young stars and SN feedback enhance [O III] emission, we might expect that UV selection biases galaxies towards having high [O III]/[C II]. Vallini et al. (2021), using the model of Ferrara et al. (2019) also found that extreme gas conditions due to starbursts in UV bright galaxies lead to high [O III]/[C II]. Furthermore, galaxies that are more disrupted by feedback are likely to have less dust attenuation which favours brighter magnitudes at fixed SFR. This in turn suggests that galaxies not selected in the UV may not exhibit [C II] deficits. While the sample is still small, we point out that two galaxies at $z > 6$ that were selected based on blue, rest-frame optical colours for ALMA follow-up are more consistent with, or perhaps even above the local [C II]–SFR relation (Smit et al. 2018). We do note that these galaxies may host AGN which have not been explicitly modelled in our work.

Finally, our models only consider stars as the sources of ionizing radiation; however, if there are a substantial number of accreting black holes, either low-mass or obscured, they may evade observational detection while providing additional photons and feedback needed to increase [O III]. Follow-up observations of additional lines, such as [NV] may help constrain this scenario.

# 4 DISCUSSION

## 4.1 Predicting galaxy properties with [O III] and [C II]

As was shown in Section 3, the [O III]/[C II] ratio is sensitive to many galaxy properties including the LyC escape fraction, ionization parameter, ISM density, SFR, etc. [C II] and [O III] are already used as SFR indicators (e.g. De Looze et al. 2014) whereas Katz et al. (2020b) showed that the combination of [O III] and [C II] could be used to identify LyC leakers. Here, we aim to constrain the individual properties of high-redshift galaxies using the correlations found in SPHINX$^{20}$.

In Fig. 17, we show 2D histograms of [C II] versus [O III] where we have coloured each individual pixel by the mean value of $f_{esc}$ (left-hand panel), 90th percentile H II region ionization parameter

(centre), and 90th percentile H II region density (right-hand panel) for simulated galaxies at $6 \le z \le 10$. The trends that were discussed in Section 3 are immediately visible. For example, the bias towards higher $f_{esc}$ and $q$ at high [O III] for fixed [C II]. We show the model for the core-collapse SN abundances with the four times enhanced metallicity as this is the model that best reproduces the [C II]–SFR and [O III]–SFR relations at high redshift. We have fit relations between each parameter and [C II] and [O III] using a generalized linear model such that:

$$ \log_{10}(f_{esc}) = \text{Min} \left[ \text{Max} \left[ -0.91 + 0.95 \log_{10} \left( \frac{L_{[OIII]}}{L_{\odot}} \right) \right. \right. \\ \left. \left. - 1.11 \log_{10} \left( \frac{L_{[CII]}}{L_{\odot}} \right), -10 \right], 0 \right], \quad (4) $$

$$ \log_{10} \left( \frac{q}{\text{cm s}^{-1}} \right) = \text{Min} \left[ \text{Max} \left[ 8.66 + 0.27 \log_{10} \left( \frac{L_{[OIII]}}{L_{\odot}} \right) \right. \right. \\ \left. \left. - 0.29 \log_{10} \left( \frac{L_{[CII]}}{L_{\odot}} \right), 7.5 \right], 9.5 \right], \quad (5) $$

and

$$ \log_{10} \left( \frac{n_H}{\text{cm}^{-3}} \right) = \text{Min} \left[ \text{Max} \left[ 1.63 - 0.18 \log_{10} \left( \frac{L_{[OIII]}}{L_{\odot}} \right) \right. \right. \\ \left. \left. + 0.21 \log_{10} \left( \frac{L_{[CII]}}{L_{\odot}} \right), 0 \right], 4 \right]. \quad (6) $$

We have defined these relations over fixed intervals for each parameter by the Min and Max functions in each equation. We have underlaid each of these relations in Fig. 17 where one can see that the generalized linear models provide an adequate representation of the simulated galaxies. To determine the uncertainty of these relations, we calculate the expected $f_{esc}$, $q$, and $n_H$ for each of the simulated galaxies and measure the difference between the true and predicted values. The 50th percentile errors on the parameters are 0.006 for $f_{esc}$, 0.47 dex for $n_H$, and 0.29 dex for $q$.

Also shown in Fig. 17 are nine $z > 6$ galaxies with [C II] and [O III] detections from Carniani et al. (2020). Six of the nine galaxies have [C II] and [O III] luminosities that fall either outside or very close to the end of the regions of the [C II]–[O III] plane probed by SPHINX$^{20}$, while three galaxies, including MACS1149-JD1, the highest redshift galaxy in the sample, fall well within range of the simulation.

In Table 1, we list the values of $f_{esc}$, $q$, and $n_H$ for each of the nine observed high-redshift galaxies as predicted by our model. Interestingly, these galaxies all have $f_{esc} \lesssim 5$ per cent, which is below the luminosity weighted values at all redshifts in SPHINX$^{20}$ and more consistent with the globally averaged value of $f_{esc}$ (Rosdahl et al., in preparation). Such values are typical of more massive galaxies at $z > 6$ in our simulations. These nine galaxies all have 90th percentile ionization parameters with $\log_{10}(\bar{\bar{q}}/\text{cm s}^{-1}) > 8.5$ and 90th percentile H II region densities of $\log_{10}(n_H/\text{cm}^{-3}) > 1.5$. The LyC escape fractions scatter below that expected from the relation in Fig. 9 due to the mass dependence of this quantity. The relations presented in Section 3 are weighted by number of galaxies rather than luminosity and are thus dominated by lower mass systems compared to what has been observed.

Compared to lower redshift galaxies, we find that our estimates for the ionization parameters of the observed high-redshift galaxies tend to be more typical of $z \sim 2-3$ Ly $\alpha$ emitting galaxies rather than local galaxies from SDSS (see table 3 of Sanders et al. 2015 and Fig. 13). We emphasize that our values of $q$ and $n_H$ may not be directly comparable to those constrained by observations because of the different methods of measuring ionization parameter. Our values represent the mass-weighted values of $q$ and $n_H$ in H II regions
<figure><img src="image_18.png" /><figcaption>**Figure 17.** 2D histograms of [C II] versus [O III] coloured by the mean values of $f_{esc}$ (left-hand panel), 90th percentile ionization parameter (centre), and 90th percentile ISM density (right-hand panel) for galaxies in SPHINX$^{20}$ at $6 \le z \le 10$. The underlaid colour map shows the results from the fitted generalized linear model for each parameter. We have shown the results for the model that assumes the core-collapse SN abundances with the metallicity scaled by a factor of 4 as this best reproduces the high-redshift [C II]–SFR and [O III]–SFR relations. For comparison, we show nine $z > 6$ galaxies from Carniani et al. (2020).</figcaption></figure>

<table><caption>**Table 1.** Estimated ISM and H II region properties for nine $z > 6$ galaxies. Redshifts, [C II], and [O III] luminosities are the surface brightness dimming corrected values from Carniani et al. (2020) while $f_{esc}$, ionization parameter and density and predicted from our models. Ionization parameter and density are estimated for H II regions.</caption><thead><tr><th>Galaxy</th><th>Redshift</th><th>L<sub>[C II]</sub>/10<sup>8</sup>L<sub>☉</sub></th><th>L<sub>[O III]</sub>/10<sup>8</sup>L<sub>☉</sub></th><th>f<sub>esc</sub></th><th>log<sub>10</sub>(q̄/cm s<sup>-1</sup>)</th><th>log<sub>10</sub>(n<sub>H</sub>/cm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>MACS1149-JD1</td><td>9.11</td><td>0.21 ± 0.05</td><td>0.88 ± 0.19</td><td>0.021 ± 0.009</td><td>8.79 ± 0.04</td><td>1.9 ± 0.04</td></tr><tr><td>A2744-YD4</td><td>8.38</td><td>0.28 ± 0.09</td><td>0.88 ± 0.21</td><td>0.017 ± 0.021</td><td>8.77 ± 0.05</td><td>1.93 ± 0.05</td></tr><tr><td>MACS0416-Y1</td><td>8.31</td><td>1.73 ± 0.25</td><td>14.0 ± 4.0</td><td>0.022 ± 0.007</td><td>8.92 ± 0.05</td><td>1.96 ± 0.03</td></tr><tr><td>SXDF-NB1006-2</td><td>7.21</td><td>3.0 ± 0.7</td><td>13.1 ± 2.9</td><td>0.012 ± 0.005</td><td>8.86 ± 0.04</td><td>2.03 ± 0.04</td></tr><tr><td>B14-65666</td><td>7.16</td><td>13.9 ± 1.8</td><td>39.0 ± 5.0</td><td>0.005 ± 0.001</td><td>8.84 ± 0.02</td><td>2.14 ± 0.02</td></tr><tr><td>BDF-3299</td><td>7.11</td><td>1.16 ± 0.16</td><td>2.13 ± 0.24</td><td>0.006 ± 0.001</td><td>8.73 ± 0.02</td><td>2.04 ± 0.02</td></tr><tr><td>J0217-0208</td><td>6.2</td><td>15.6 ± 2.2</td><td>89.3 ± 2.1</td><td>0.01 ± 0.002</td><td>8.94 ± 0.02</td><td>2.1 ± 0.02</td></tr><tr><td>J0235-0532</td><td>6.09</td><td>5.1 ± 0.8</td><td>39.9 ± 3.1</td><td>0.017 ± 0.003</td><td>8.95 ± 0.02</td><td>2.02 ± 0.02</td></tr><tr><td>J1211-0118</td><td>6.03</td><td>15.1 ± 1.1</td><td>51.0 ± 7.0</td><td>0.006 ± 0.001</td><td>8.87 ± 0.02</td><td>2.13 ± 0.01</td></tr></tbody></table>

which is represented by gas with $n_H > 1 \text{ cm}^{-3}$, $T < 10^{5.5}$ K, and $x_{HII} > 50$ per cent.

We highlight that estimates of gas density in these high-redshift galaxies are not always consistent between different methods in the literature. For example for MACS1149_JD1, Yang & Lidz (2020) find an upper limit gas density of $\log_{10}(n_H/cm^{-3}) = 2.4-2.8$ or $1.6-2.0$ depending on how scatter between H II regions is accounted for, whereas Vallini et al. (2021) estimate a value of $0.88^{+0.55}_{-0.19}$. We find a 90th percentile H II region (log) density of 1.69 for the same galaxy. These density estimates are not all measuring the exact same quantity. The Vallini et al. (2021) values represent the average over the [O III] and [C II] emitting regions while the Yang & Lidz (2020) values represent the upper limits in H II regions. Our 90th percentile values are perhaps more comparable to the upper limits quoted in Yang & Lidz (2020). Consistent with Vallini et al. (2021), we find that including both H II regions and PDRs in the calculation increases the density estimate.

## 4.2 The case for [O I] observations

Much of the complexity in interpreting ISM properties in real high-redshift galaxies using [O III] and [C II] is that without knowing the

star formation history and metallicity of the system, the different time-scales for O and C enrichment can impact the determination of galaxy properties. Values of $f_{esc}$, q, and $\rho$ are all sensitive to the C/O abundance ratio assumed. One potential way to circumvent this issue is to focus on line ratios that are not as impacted by the exact metal abundance ratios.

The [O I] 63 µm line has many similar properties to [C II] in that it is predominantly emitted by neutral gas since the ionization energy of O I is 13.6 eV (nearly exactly that of H I). [O I] has a very high critical density (~$5 \times 10^5$ cm$^{-3}$) and is known to strongly correlate with star formation in the local Universe (De Looze et al. 2014). Numerical simulations predict that [O I] can outshine [C II] at high redshift (e.g. Olsen et al. 2017; Katz et al. 2019b) and an [O I]/[C II] ratio >1 was recently detected with ALMA at $z > 6$ (Rybak et al. 2020). [O I] is however impacted by self-absorption (Kraemer, Jackson & Lane 1996; Goldsmith 2019). The primary benefit of observing [O I] is that the comparison with [O III] should not be as impacted by different enrichment channels.

In Fig. 18 we show a 2D histogram of [O I] versus [O III] coloured by the mean value of $f_{esc}$ for galaxies in SPHINX$^{20}$ at $6 \le z \le 10$. Similar to [C II], galaxies with higher [O III] at fixed [O I] exhibit higher escape fractions due to the fact that [O I] originates in
<figure><img src="image_19.png" /><figcaption>**Figure 18.** 2D histogram of [O I] versus [O III] coloured by the mean value of $f_{esc}$ for galaxies in SPHINX$^{20}$ at $6 \le z \le 10$. The underlaid colour map shows results from the fitted generalized linear model. We have shown the results for the model that assumes core-collapse SN abundances with the metallicity scaled by a factor of 4.0.</figcaption></figure>

predominantly neutral gas whereas [O III] is emitted by ionized gas. The relation between [O III], [O I], and $f_{esc}$ is well defined by the generalized linear model:

$$ \log_{10}(f_{esc}) = \text{Min} \left[ \text{Max} \left[ -0.31 + 1.02 \log_{10} \left( \frac{L_{[OIII]}}{L_{\odot}} \right) \right. \right. \\ \left. \left. - 1.07 \log_{10} \left( \frac{L_{[OI]}}{L_{\odot}} \right), -10 \right], 0 \right]. \quad (7) $$

The coefficients in this equation are nearly identical to the relation between [O III], [C II], and $f_{esc}$, while the intercept has increased due to the fact that [O I] is brighter. Given the utility of [O I] as a supplement or replacement for [C II], we argue that more observation time should be used to constrain [O I] emission from $z > 6$ galaxies.

## 4.3 Caveats

As with all numerical simulations, there are numerous caveats we must be aware of when interpreting our results. We have modelled the emission lines in post-processing with CLOUDY and our method for setting up the models can introduce systematic uncertainties (see e.g. Olsen et al. 2018). For example, we have assumed a gas slab being irradiated from one direction whereas one could assume a shell instead. Simulation cells that host star particles are more likely to act like shells whereas gas cells without stars are more similar to slabs. Furthermore, we employ a constant temperature assumption to account for shocks and SN feedback that are present in the simulation rather than allow CLOUDY to compute the photoionization-equilibrium temperature of the gas. Such an assumption can lead to different line luminosities as Lupi et al. (2020) showed that a variable temperature CLOUDY model can increase the [O III]/[C II] ratio by a factor of two. One could also choose a different shape for the input SED than used in this work which can impact both the ionization states of individuals as well as the electron density.

Rather than model a different SED for each gas cell, we have fixed the shape of the SED within each radiation bin to be consistent with the Milky Way ISM and have normalized the photon count in each photon energy bin to match the simulation to conserve energy. The choice of SED within the cosmological simulation itself is also uncertain and it is well established that depending on the choice, the number of ionizing photons as well as the escape fraction can change significantly (e.g. Ma et al. 2016; Rosdahl et al. 2018). The SED will also vary with IMF as well as metallicity. Our results may change if the model with a top-heavy IMF and enhanced metallicity is self-consistently followed in the simulation. During post-processing, we also modelled the sub-ionizing bands in the optically thin limit when in reality, molecules like H$_2$ or dust particles can absorb such photons and perhaps change the C$^+$ content of the galaxy. Finally, additional physics such as cosmic rays, magnetic fields, and subgrid turbulence are not included in the calculation, all of which can impact the emergent luminosities.

Our simulations are limited by both spatial and mass resolution which prevents us from modelling the detailed physics of the ISM on scales that we do not resolve. While our resolution is state-of-the-art in the context of full-box radiation hydrodynamics simulations, we do not resolve the formation of individual stars and their individual H II regions/PDRs where [O III] and [C II] are emitted. Rather, we represent stellar clusters and therefore the structure of the ISM in the simulation may be different than in high-redshift galaxies. Furthermore, because we use a quasi-Lagrangian scheme for refining the AMR grid, most of our resolution tends to be concentrated on regions of high density. The origin of large [C II] haloes at high redshift is currently unknown and one suggestion is that there are small clumps of gas or small satellite galaxies in the outer regions of the haloes that emit [C II] (e.g. Fujimoto et al. 2019). Simulations that add additional resolution to the CGM (e.g. Hummels et al. 2019; van de Voort et al. 2019) tend to resolve significantly more of these clumps.

Finally, we have assumed specific subgrid models for star formation and stellar feedback, all of which impact the stellar mass, metal enrichment, and ionizing photon production. As we have argued, a more top-heavy stellar IMF which leads to enhanced SN feedback, ionizing photon production, and metal enrichment may go a long way to reconciling the high-redshift [C II]–SFR, [O III]–SFR, and [O III]/[C II]–SFR relations. However, abundance ratios change for different mass stars and modifications to the assumed metal abundances ratios will result from a change in the stellar IMF. Enhancements in metallicity will also lead to additional cooling which is not completely self-consistent with our post-processing method. Our goal with this work is not to demonstrate that an exact metallicity increase of a factor of 4 is needed. Rather, we simply mean to demonstrate that increasing the metallicity from what is in our fiducial SPHINX model will push the simulated galaxies in the direction of observations. Future simulations will be needed to study these effects.

## 5 CONCLUSIONS

In this work, we have studied the possible origins of high [O III]/[C II] ratios observed in $z > 6$ galaxies. We have used a new, state-of-the-art, full-box, cosmological radiation hydrodynamics simulation – SPHINX$^{20}$ – to model [C II] and [O III] for more than 20 000 simulated star-forming galaxies at $z > 6$, allowing us to disentangle the physics that drives high [O III]/[C II] ratios. Our main conclusions can be summarized as follows:
(i) We can only reproduce the observed [C II]–SFR, [O III]–SFR, and [O III]/[C II]–SFR relations at $z > 6$ if we assume that high-redshift galaxies are primarily enriched by low-metallicity core-collapse SN and that the high-redshift stellar IMF is significantly more top heavy than the Kroupa IMF. This is consistent with expectations for low-metallicity star formation in the early Universe.

(ii) The [C II]–SFR relation naturally converges to the observed values at $z \sim 4.5$ when Solar abundance ratios with a standard IMF are assumed, indicating that both the stellar IMF and enrichment channels evolve significantly as redshift decreases.

(iii) While the reduced C/O ratios at $z > 6$ may be the primary reason why high-redshift galaxies exhibit high [O III]/[C II] ratios, increased ionization parameters, increased ISM density, increased LyC escape fractions, and to a small extent CMB attenuation also contribute to the high [O III]/[C II] ratios. In contrast, lower gas metallicities, extended [C II] emission, and inclination effects are unlikely to explain the behaviour for all of the nine high-redshift galaxies with observations of both [C II] and [O III].

(iv) [O III] emitting regions have significantly lower ISM densities than regions that emit [C II] and thus 1D models that consider a fixed gas density must be treated with caution when being used to study the ISM of the early Universe.

(v) [O III] and [C II] can be used to constrain ISM properties and LyC escape fractions of high-redshift galaxies. We find that typical escape fractions of the observed systems are ~1 per cent – 4 per cent while gas densities in the [O III] and [C II] emitting regions can reach as high as >1000 cm⁻³.

(vi) [O I] can be used as both an alternative and supplement to [C II] observations in the high-redshift Universe because not only is it predicted to be brighter than [C II], it also originates from similar regions of the ISM meaning that it can be used to constrain similar galaxy characteristics as [C II]. We argue that more observation time should be spent on [O I] observations of $z > 6$ galaxies.

Our model for predicting emission lines from simulations is applicable to numerous other lines across the wavelength spectrum. In future work we aim to compare how the methods for predicting galaxy properties such as metallicity, ionization parameter, and gas density from strong-line diagnostics translate to the early Universe.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for their comments which improved the manuscript. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. MGH acknowledges support from the UKRI Science and Technology Facilities Council (grant numbers ST/N000927/1 and ST/S000623/1). Support by ERC Advanced Grant 320596 ‘The Emergence of Structure during the Epoch of reionization’ is gratefully acknowledged. SMA acknowledges support by ERC Starting Grant 638707 ‘Black holes and their host galaxies: co-evolution across cosmic time’. TK was supported by the National Research Foundation of Korea (NRF-2019K2A9A1A0609137711 and NRF-2020R1C1C100707911). TG acknowledges support from the European Research Council under grant agreement ERC-stg757258 (TRIPLE). The research of AS and JD is supported by Adrian Beecroft and STFC. Computing time for this work was provided by the Partnership for Advanced Computing in Europe (PRACE) as part of the ‘First luminous objects and reionization with SPHINX (cont.’) (2016153539, 2018184362, 2019215124) project. We thank Philipp Otte and Filipe Guimaraes for helpful support throughout the project and for the extra storage

they provided us. We also thank GENCI for providing additional computing resources under GENCI grant A0070410560. RSE was supported by the European Research Council under the European Union Horizon 2020 research and innovation programme (grant agreement no. 669253).

## DATA AVAILABILITY

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

## REFERENCES

Arata S., Yajima H., Nagamine K., Abe M., Khochfar S., 2020, MNRAS, 498, 5541

Aravena M. et al., 2016, ApJ, 833, 71

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

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

Béthermin M. et al., 2020, A&A, 643, A2

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

Boselli A., Gavazzi G., Lequeux J., Pierini D., 2002, A&A, 385, 454

Bradač M. et al., 2017, ApJ, 836, L2

Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44

Capak P. L. et al., 2015, Nature, 522, 455

Carilli C. L., Walter F., 2013, ARA&A, 51, 105

Carniani S. et al., 2020, MNRAS, 499, 5136

Chon S., Omukai K., Schneider R., 2021, MNRAS, 508, 4175

Cooke E. A. et al., 2018, ApJ, 861, 100

Cormier D. et al., 2010, A&A, 518, L57

Cormier D. et al., 2015, A&A, 578, A53

Cormier D. et al., 2019, A&A, 626, A23

De Looze I. et al., 2014, A&A, 568, A62

De Looze I., Baes M., Bendo G. J., Cortese L., Fritz J., 2011, MNRAS, 416, 2712

Eldridge J. J., Izzard R. G., Tout C. A., 2008, MNRAS, 384, 1109

Ferland G. J. et al., 2017, RMXAA, 53, 385

Ferland G. J., 1996, University of Kentucky Internal Report, Hazy, A Brief Introduction to Cloudy 90. Univ. Kentucky, Lexington

Ferrara A., Vallini L., Pallottini A., Gallerani S., Carniani S., Kohandel M., Decataldo D., Behrens C., 2019, MNRAS, 489, 1

Franeck A. et al., 2018, MNRAS, 481, 4277

Fujimoto S. et al., 2019, ApJ, 887, 107

Goldsmith P. F., 2019, ApJ, 887, 54

Goldsmith P. F., Langer W. D., Pineda J. L., Velusamy T., 2012, ApJS, 203, 13

Grevesse N., Asplund M., Sauval A. J., Scott P., 2010, Ap&SS, 328, 179

Haardt F., Madau P., 2012, ApJ, 746, 125

Hahn O., Abel T., 2011, MNRAS, 415, 2101

Harikane Y. et al., 2020, ApJ, 896, 93

Hashimoto T. et al., 2018, Nature, 557, 392

Herrera-Camus R. et al., 2015, ApJ, 800, 1

Hummels C. B. et al., 2019, ApJ, 882, 156

Inoue A. K. et al., 2016, Science, 352, 1559

Jenkins E. B., 2009, ApJ, 700, 1299

Karpov P. I., Martizzi D., Macias P., Ramirez-Ruiz E., Kolborg A. N., Naiman J. P., 2020, ApJ, 896, 66

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

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

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

Katz H. et al., 2021, MNRAS, 507, 1254

Katz H., Kimm T., Sijacki D., Haehnelt M. G., 2017, MNRAS, 468, 4831

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

Kaufman M. J., Wolfire M. G., Hollenbach D. J., 2006, ApJ, 644, 283
Kewley L. J., Dopita M. A., 2002, ApJS, 142, 35

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

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

Knudsen K. K., Richard J., Kneib J.-P., Jauzac M., Clément B., Drouart G., Egami E., Lindroos L., 2016, MNRAS, 462, L6

Kohandel M., Pallottini A., Ferrara A., Zanella A., Behrens C., Carniani S., Gallerani S., Vallini L., 2019, MNRAS, 487, 3007

Kraemer K., Jackson J. M., Lane A. P., 1996, American Astronomical Society Meeting Abstracts, 21.02

Kroupa P., 2001, MNRAS, 322, 231

Lagache G., 2018, in Jelić V., van der Hulst T., eds, Proceedings of the International Astronomical Union, IAU Symposium, Vol. 333, Peering towards Cosmic Dawn. Cambridge Univ. Press, Cambridge UK, p. 228

Laporte N. et al., 2019, MNRAS, 487, L81

Larson R. B., 1981, MNRAS, 194, 809

Le Fèvre O. et al., 2020, A&A, 643, A1

Leung T. K. D., Olsen K. P., Somerville R. S., Davé R., Greve T. R., Hayward C. C., Narayanan D., Popping G., 2020, ApJ, 905, 102

Loiacono F. et al., 2021, A&A, 646, A76

Lupi A., Volonteri M., Decarli R., Bovino S., Silk J., Bergeron J., 2019, MNRAS, 488, 4004

Lupi A., Pallottini A., Ferrara A., Bovino S., Carniani S., Vallini L., 2020, MNRAS, 496, 5160

Ma X., Hopkins P. F., Kasen D., Quataert E., Faucher-Giguère C.-A., Kereš D., Murray N., Strom A., 2016, MNRAS, 459, 3614

Madden S. C. et al., 2013, PASP, 125, 600

Maiolino R. et al., 2015, MNRAS, 452, 54

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

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

Moriwaki K. et al., 2018, MNRAS, 481, L84

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

Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ, 889, 161

Neri R., Downes D., Cox P., Walter F., 2014, A&A, 562, A35

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

Olsen K. et al., 2018, Galaxies, 6, 100

Olsen K., Greve T. R., Narayanan D., Thompson R., Davé R., Niebla Rios L., Stawinski S., 2017, ApJ, 846, 105

Ota K. et al., 2014, ApJ, 792, 34

Pallottini A. et al., 2019a, MNRAS, 487, 1689

Pallottini A. et al., 2019b, MNRAS, 487, 1689

Pallottini A., Ferrara A., Gallerani S., Vallini L., Maiolino R., Salvadori S., 2017, MNRAS, 465, 2540

Penston M. V. et al., 1990, A&A, 236, 53

Pentericci L. et al., 2016, ApJ, 829, L11

Poglitsch A., Krabbe A., Madden S. C., Nikola T., Geis N., Johansson L. E. B., Stacey G. J., Sternberg A., 1995, ApJ, 454, 293

Popping G., van Kampen E., Decarli R., Spaans M., Somerville R. S., Trager S. C., 2016, MNRAS, 461, 93

Popping G., Narayanan D., Somerville R. S., Faisst A. L., Krumholz M. R., 2019, MNRAS, 482, 4906

Prokhorenkova L., Gusev G., Vorobev A., Veronika Dorogush A., Gulin A., 2018, in Bengio S., Wallach H., Larochelle H., Grauman K., Cesa-Bianchi N., Garnett R., eds, CatBoost: Unbiased Boosting with Categorical Features, 32nd Conference on Neural Information Processing Systems (NeurIPS 2018). Curran Associates, Inc., Montréal, Canada

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

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

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

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

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

Rybak M., Zavala J. A., Hodge J. A., Casey C. M., Werf P. v. d., 2020, ApJ, 889, L11

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

Sanders R. L. et al., 2015, ApJ, 799, 138

Schaerer D. et al., 2020, A&A, 643, A3

Schaerer D., Boone F., Zamojski M., Staguhn J., Dessauges-Zavadsky M., Finkelstein S., Combes F., 2015, A&A, 574, A19

Silva M., Santos M. G., Cooray A., Gong Y., 2015, ApJ, 806, 209

Smit R. et al., 2018, Nature, 553, 178

Spitzer L., 1978, Physical Processes in the Interstellar Medium. Wiley, Hoboken, NJ

Stacey G. J., Geis N., Genzel R., Lugten J. B., Poglitsch A., Sternberg A., Townes C. H., 1991, ApJ, 373, 423

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

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

Steidel C. C., Strom A. L., Pettini M., Rudie G. C., Reddy N. A., Trainor R. F., 2016, ApJ, 826, 159

Strom A. L., Steidel C. C., Rudie G. C., Trainor R. F., Pettini M., Reddy N. A., 2017, ApJ, 836, 164

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Sutherland R., Dopita M., Binette L., Groves B., 2018, MAPPINGS V: Astrophysical Plasma Modeling Code, record ascl:1807.005

Swinbank A. M. et al., 2012, MNRAS, 427, 1066

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

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

Trebitsch M., Blaizot J., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 470, 224

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

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

Uzgil B. D. et al., 2021, ApJ, 912, 67

Vallini L., Gallerani S., Ferrara A., Pallottini A., Yue B., 2015, ApJ, 813, 36

Vallini L., Ferrara A., Pallottini A., Carniani S., Gallerani S., 2021, MNRAS, 505, 5543

van de Voort F., Springel V., Mandelker N., van den Bosch F. C., Pakmor R., 2019, MNRAS, 482, L85

Veronika Dorogush A., Ershov V., Gulin A., 2018, preprint (arXiv:1810.1136)

Yang S., Lidz A., 2020, MNRAS, 499, 3417

Yue B., Ferrara A., 2019, MNRAS, 490, 1928

This paper has been typeset from a *TEX/LATEX* file prepared by the author.