# Mg ii in the JWST era: a probe of Lyman continuum escape?

Harley Katz, Thibault Garel, Joakim Rosdahl, Valentin Mauerhofer, Taysun Kimm, Joakim A Rosdahl, Léo Michel-Dansac, Julien Devriendt, Adrianne Slyz, Martin Haehnelt

► To cite this version:

Harley Katz, Thibault Garel, Joakim Rosdahl, Valentin Mauerhofer, Taysun Kimm, et al. Mg ii in the JWST era: a probe of Lyman continuum escape?. Monthly Notices of the Royal Astronomical Society, 2022, 515 (3), pp.4265-4286. 10.1093/mnras/stac1437. hal-03692989

HAL Id: hal-03692989

https://hal.science/hal-03692989

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
# Mg II in the JWST era: a probe of Lyman continuum escape?

Harley Katz,<fnref n="1" />★ Thibault Garel,<fnref n="2" /> Joakim Rosdahl,<fnref n="3" /> Valentin Mauerhofer,<fnref n="2" />,<fnref n="3" /> Taysun Kimm,<fnref n="5" />,⁶
Jérémy Blaizot,<fnref n="3" /> Léo Michel-Dansac,<fnref n="3" /> Julien Devriendt,<fnref n="1" /> Adrianne Slyz<fnref n="1" /> and Martin Haehnelt<fnref n="5" />

Accepted 2022 May 17. Received 2022 March 26; in original form 2021 November 29

## ABSTRACT

Limited constraints on the evolution of the Lyman Continuum (LyC) escape fraction represent one of the primary uncertainties in the theoretical determination of the reionization history. Due to the intervening intergalactic medium, the possibility of observing LyC photons directly in the epoch of reionization is highly unlikely. For this reason, multiple indirect probes of LyC escape have been identified, some of which are used to identify low-redshift LyC leakers (e.g. O32), while others are primarily useful at $z > 6$ (e.g. [O III]/[C III] far-infrared emission). The flux ratio of the resonant Mg II doublet emission at 2796Å and 2803Å as well as the Mg II optical depth have recently been proposed as ideal diagnostics of LyC leakage that can be employed at $z > 6$ with JWST. Using state-of-the-art cosmological radiation hydrodynamics simulations post-processed with CLOUDY and resonant-line radiative transfer, we test whether Mg II is indeed a useful probe of LyC leakage. Our simulations indicate that the majority of bright, star-forming galaxies with high LyC escape fractions are expected to be Mg II emitters rather than absorbers at $z = 6$. However, we find that the Mg II doublet flux ratio is a more sensitive indicator of dust rather than neutral hydrogen, limiting its use as a LyC leakage indicator to only galaxies in the optically thin regime. Given its resonant nature, we show that Mg II will be an exciting probe of the complex kinematics in high-redshift galaxies in upcoming JWST observations.

**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

The Lyman continuum (LyC) escape fraction ($f_{esc}$) represents one of the fundamental uncertainties in constraining the history of reionization. Due to the intervening intergalactic medium, direct observations of escaping LyC photons are essentially impossible at $z > 6$. For this reason, leaking LyC radiation is often studied either via observations of lower redshift galaxies – e.g. at $z \sim 2-3$ (Siana et al. 2015; Steidel et al. 2018; Fletcher et al. 2019; Nakajima et al. 2020) or at $z \lesssim 1$ (e.g. Borthakur et al. 2014; Izotov et al. 2016a, 2018b; Leitherer et al. 2016) – that are expected to be analogues of similar systems at high-redshift, or via numerical simulations that aim to model the physics that enables LyC photons to escape from galaxies (e.g. Paardekooper, Khochfar & Dalla Vecchia 2015; Xu et al. 2016; Rosdahl et al. 2018).

Although direct measurements of LyC radiation cannot be made in the epoch of reionization, numerous proxies exist that can be used to constrain $f_{esc}$ in individual high-redshift galaxies. For example, Zackrisson, Inoue & Jensen (2013), Zackrisson et al. (2017) showed that LyC leakers segregate from non-leakers on the H β equivalent width (EW) and UV slope (β) plane. Katz et al. (2020b)

demonstrated that LyC leakers could be identified based on high ratios of [O III]$_{88\mu m}$/[C III]$_{158\mu m}$ (see also Inoue et al. 2016). Other diagnostics such as high O32 (e.g. Nakajima & Ouchi 2014; Izotov et al. 2016a, 2018b; Nakajima et al. 2020), S II deficits (e.g. Wang et al. 2019, 2021), and Lyα peak separation (e.g. Verhamme et al. 2015, 2017; Izotov et al. 2020) and equivalent width (EW) (e.g. Steidel et al. 2018) have all been successfully employed in the low-redshift Universe. Furthermore, as an alternative to emission, UV absorption lines can be used to either rule out LyC leakage or select high-$f_{esc}$ galaxies (Mauerhofer et al. 2021). Similarly, one can correlate the positions of Lyman break galaxies (LBGs) and the Lyα spectra of background quasars to obtain a statistical measure of $f_{esc}$ directly in the epoch of reionization (Kakiichi et al. 2018).

Recently, Mg II emission has been suggested as a possible probe of $f_{esc}$ (Henry et al. 2018; Chisholm et al. 2020). Mg II can exhibit both emission and absorption lines at 2796.35Å and 2803.53Å due to electrons transitioning from the first excited state to the ground state. The ionization energies of Mg I is ~7.6 eV, well below that of neutral hydrogen, and for Mg II it is ~15 eV. Therefore, Mg II is expected to primarily trace neutral gas. As a resonant line, Mg II has many analogous properties to Lyα emission. However, one key difference is that Mg II emission is not strongly affected by IGM absorption. Thus, if Mg II can be used as a proxy for $f_{esc}$, with the upcoming

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

<footnote n="2">²Observatoire de Genève, Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

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

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

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

<footnote>* E-mail: harley.katz@physics.ox.ac.uk</footnote>
launch of JWST, it represents a promising method that can be used on galaxies at $z > 6$ (Chisholm et al. 2021).

Chisholm et al. (2020) demonstrated that Mg II can be used to constrain $f_{esc}$ in two ways. If one can measure the escape fraction of Mg II emission, then the Mg II optical depth ($\tau_{Mg\text{II}}$) can be converted into a neutral hydrogen optical depth based on the metallicity of the galaxy, which can subsequently be combined with the estimated dust properties of the system (e.g. from the Balmer decrement) to measure the neutral hydrogen column density and the LyC escape fraction. The first method follows Henry et al. (2018) who showed that the intrinsic Mg II emission can be calculated using $[O\,III] 5007\AA$ and $[O\,II] 3727\AA$ doublet emission. The key to this calculation is measuring $\tau_{Mg\text{II}}$. For the second method, Chisholm et al. (2020) also demonstrated that $\tau_{Mg\text{II}}$ can be calculated by using the ratio of the Mg II doublet emission lines. Since the oscillator strength of the 2796$\AA$ line is twice that of the 2803$\AA$ line, the absorption cross-section is twice as large for the former compared to the latter. In a simple model where the intrinsic strengths of the two lines are set by collisional excitation, in the absence of dust, if one places a screen of Mg II gas in front of the source, the emergent ratio of the line fluxes will be sensitive to $\tau_{Mg\text{II}}$ because photons of the shorter wavelength line are expected to scatter more often out of the observer's line of sight. Thus Chisholm et al. (2020) show that the line ratio method can also successfully predict the LyC $f_{esc}$ in observed low-redshift galaxies.

While these methods involving Mg II have shown promise, there are numerous factors that either make them difficult to use in practice or that introduce systematic uncertainties.

(i) Henry et al. (2018) used 1D photoionization models to calculate how the intrinsic Mg II emission scales with various oxygen lines. Whether the calibration based on these simple models applies to galaxies that exhibit a wide range of ISM properties remains to be determined.

(ii) The Chisholm et al. (2020) calculation hinges on the geometry of a system being that of a screen of Mg II in front of a source. The ratio of the two emission lines only changes as a function of optical depth precisely because photons are scattered out of the observer's line of sight. If one assumes a different geometry, for example a shell, regardless of the optical depth, for pure Mg II gas, the intrinsic line ratio will be preserved. The presence of dust will further complicate the calculation. Since the shorter wavelength line will scatter more often than the longer wavelength line, in the optically thick regime the probability of being absorbed by dust is higher. In this case, the line ratio is sensitive to both the dust properties and the Mg II optical depth, therefore making it more difficult to calculate the LyC $f_{esc}$.

(iii) Even if Mg II emission was the perfect tracer of escaping LyC radiation, not all galaxies are Mg II emitters. Feltre et al. (2018) studied a sample of 381 galaxies at $0.7 < z < 2.34$ and showed that among the galaxies with a Mg II detection, ~50 per cent were emitters while the others were either absorbers or exhibited P-Cygni profiles in their spectra. The fraction of galaxies at $z > 6$ that are Mg II emitters versus absorbers is currently unknown. As such, any biases in using $f_{esc}$ measurements from Mg II emitters at high redshift as a representative value for the whole galaxy population must be elucidated before it can be trusted as a diagnostic.

In this work, we study Mg II in the epoch of reionization using a state-of-the-art cosmological radiation hydrodynamics simulation, SPHINX²⁰ (Rosdahl et al., in preparation; Katz et al. 2022). Our goals are

(i) to determine whether the Mg II escape fraction is correlated with the LyC escape fraction,

(ii) to test whether the Henry et al. (2018) method for measuring the intrinsic Mg II emission applies to high-redshift galaxies,

(iii) to test whether the Chisholm et al. (2020) model of using the line ratio of Mg II doublet emission provides an accurate estimate of the LyC escape fraction in high-redshift galaxies, and

(iv) to determine the fraction of high-redshift galaxies that are Mg II emitters, absorbers, or exhibit a P-Cygni profile.

This paper is organized as follows. In Section 2 we introduce the SPHINX²⁰ simulations and describe the methods for emission line modelling and resonant line radiative transfer. In Section 3, we compare the simulated galaxies with low-redshift green peas and blueberries and discuss the utility of Mg II as a tracer of LyC escape. Finally, in Sections 4 and 5 we present our caveats and conclusions.

# 2 NUMERICAL METHODS

## 2.1 Cosmological simulations

This work makes use of the SPHINX²⁰, 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 that are part of the SPHINX project. The details of SPHINX²⁰ are well described in Rosdahl et al. (2018), Katz et al. (2022), and Rosdahl et al. (in preparation).

Briefly, the simulations are run with RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015), a radiation hydrodynamics extension of the RAMSES code (Teyssier 2002). The simulation volume has a comoving side length of 20 Mpc. Initial conditions for $1024^3$ dark matter particles and gas cells were generated with MUSIC (Hahn & Abel 2011) assuming the following cosmology: $\Omega_\Lambda = 0.6825$, $\Omega_m = 0.3175$, $\Omega_b = 0.049$, $h = 0.6711$, and $\sigma_8 = 0.83$. The dark matter particle mass is $2.5 \times 10^5$ M$_\odot$ and gas cells are allowed to adaptively refine up to a maximum physical resolution of $7.3\text{pc} h^{-1}$ at $z = 6$. The initial composition of the gas is set to be 76 per cent H and 24 per cent He by mass, with a small initial metallicity of $3.2 \times 10^{-4} Z_\odot$.

Star formation is modelled following a thermoturbulent prescription (Kimm et al. 2017; Trebitsch et al. 2017; Rosdahl et al. 2018). Star particles are allowed to form by drawing from a Poisson distribution in multiples of 400 M$_\odot$. Star particles impact the gas via gravity, supernova (SN) feedback (Kimm et al. 2015), and LyC radiation feedback (photoionization, photoheating, radiation pressure). Star particles inject ionizing photons into their host cells as a function of their mass, age, and metallicity assuming the BPASS (Eldridge, Izzard & Tout 2008; Stanway, Eldridge & Becker 2016) spectral energy distribution (SED) model.

Non-equilibrium chemistry is followed locally for H I, H II, $e^-$, He I, He II, and He III. Gas cooling is calculated for primordial channels (see the appendix of Rosdahl et al. 2013) as well as for metal lines (Rosen & Bregman 1995; Ferland 1996).

Since SPHINX²⁰ only includes radiative transfer for photons with energies $\ge 13.6$ eV, the simulation is post-processed in the optically thin limit at $z = 10, 9, 8, 7, 6, 5$, and $4.68$ with RAMSES-RT including two additional photon energy bins (5.6–11.2 eV and 11.2–13.6 eV). These radiation bins are not included in the main simulation because the computational cost of the RT scales with the number of radiation bins. This is necessary to be able to model the ionization of metals with ionization energies below 13.6 eV (e.g. the transition from Mg I to Mg II).

SPHINX²⁰ is evolved to a final redshift of $z = 4.64$ and semirevolves the detailed ISM and CGM structure of tens of thousands
<table><caption>**Table 1.** Emission lines considered in this work. From left to right we list the element and ionization state or common name for the line, the rest-frame wavelength of the line, if the intrinsic luminosity was computed with CLOUDY, if the intrinsic luminosity was computed with an analytic model, the possible absorbers for photons at the wavelength of the line, whether the line is resonant, and the reason why the line was included in this work.</caption><thead><tr><th>Line</th><th>Rest frame<br>wavelength (Å)</th><th>CLOUDY</th><th>Analytic</th><th>Absorbers</th><th>Resonant</th><th>Reason</th></tr></thead><tbody><tr><td>Lyα</td><td>1215.67</td><td>✗</td><td>✓</td><td>H I, Dust</td><td>✓</td><td>For comparison as an alternative f<sub>esc</sub> diagnostic</td></tr><tr><td>[C III]</td><td>1906.68</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>Alternative for O3 to calibrate intrinsic Mg II emission.</td></tr><tr><td>C III]</td><td>1908.73</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>Alternative for O3 to calibrate intrinsic Mg II emission.</td></tr><tr><td>Mg II</td><td>2796.35</td><td>✓</td><td>✗</td><td>Mg II, Dust</td><td>✓</td><td>–</td></tr><tr><td>Mg II</td><td>2803.53</td><td>✓</td><td>✗</td><td>Mg II, Dust</td><td>✓</td><td>–</td></tr><tr><td>[O II]</td><td>3726.03</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>For calibrating intrinsic Mg II emission.</td></tr><tr><td>[O II]</td><td>3728.81</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>For calibrating intrinsic Mg II emission.</td></tr><tr><td>Hβ</td><td>4861.32</td><td>✗</td><td>✓</td><td>Dust</td><td>✗</td><td>Balmer decrement needed for dust measurement</td></tr><tr><td>[O III]</td><td>4958.91</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>For calibrating intrinsic Mg II emission.</td></tr><tr><td>[O III]</td><td>5006.84</td><td>✓</td><td>✗</td><td>Dust</td><td>✗</td><td>For calibrating intrinsic Mg II emission.</td></tr><tr><td>Hα</td><td>6562.80</td><td>✗</td><td>✓</td><td>Dust</td><td>✗</td><td>Balmer decrement needed for dust measurement</td></tr></tbody></table>

of galaxies, making it an ideal tool for studying emission lines in the epoch of reionization. In this work, we will focus only on the more massive galaxies in the $z = 6$ snapshot. Haloes are found using the ADAPTAHOP halo finder (Aubert, Pichon & Colombi 2004; Tweed et al. 2009) in the most massive submaxima (MSM) mode (see Rosdahl et al. 2018 for more details). For this work, we consider all gas and star particles inside the virial radius when computing emission lines and continuum emission (see below) and do not separate subhaloes.

## 2.2 Emission-line, stellar continuum, and escape fraction modelling

In this work, we consider both resonant and non-resonant emission lines, as well as continuum radiation from stars. The emergent luminosities of each emission line, the emergent stellar continuum, and the escape fractions are computed by post-processing SPHINX²⁰ with the Monte Carlo radiative transfer code RASCAS (Michel-Dansac et al. 2020). Intrinsic emission line luminosities are either computed analytically or with CLOUDY (Ferland et al. 2017), while the intrinsic stellar continuum is given by the BPASS SED used in the cosmological simulation.

We consider 11 different emission lines in this work as listed in order of increasing rest-frame wavelength in Table 1. Intrinsic luminosities from each gas cell for primordial species (i.e. H and He) are computed analytically while intrinsic metal-line luminosities are numerically computed with the combination of CLOUDY and the machine learning method described in Katz et al. (2019, 2022).

In order to test whether Mg II is a useful tracer of the LyC escape fraction at high redshift, we must model the resonant line radiation transfer of the 2796.35Å and 2803.53Å lines. This is a crucial step because CLOUDY does not model the radiative transfer of the Mg II photons through the ISM and CGM of the galaxies. CLOUDY only predicts the intrinsic emission from each gas cell in the simulation. Regardless of the method used to measure the intrinsic Mg II emission, for the Chisholm et al. (2020) method, one must also know the dust attenuation in order to measure $f_{esc}$. Hence we also model H α and H β, which can be used to measure the Balmer decrement. Since the Henry et al. (2018) method for estimating the intrinsic Mg II emission relies on the O32 diagnostic, we also model the [O II] 3726.03Å, 3728.81Å doublet as well as the [O III] 5006.84Å and [O III] 4958.91Å lines. Due to their longer wavelengths, the [O III] lines will drop out of the JWST filters at lower redshift compared to the [O II] and Mg II lines. For this reason, we also model C III] 1908Å

doublet emission to determine whether these lines can be used as a replacement for [O III] at high-redshift to measure the intrinsic Mg II emission. Finally, we also model Ly α emission and radiative transfer to determine which of the two (Ly α or Mg II) provides a more faithful representation of $f_{esc}$. This is useful for deciding where to spend telescope time, especially in the low-redshift Universe where Ly α is less subject to IGM attenuation. Full details of our method for computing intrinsic luminosities, and the nearby stellar continuum are provided in Appendix A.

Once the intrinsic luminosities and stellar continuum have been calculated for each cell and star particle in every halo, the emergent luminosities are measured by using Monte Carlo radiative transfer. We employ different approaches depending on the line.

For each line, we first generate a set of initial conditions consisting of positions, initial directions, and frequencies. Initial positions are sampled from a multinomial distribution across all cells or star particles where the weights on each cell or star particle are computed as the fraction of the total intrinsic luminosity that each individual cell or star particle represents. Initial directions are randomized across an isotropic sphere. For cells, frequencies are drawn from a Gaussian profile with line-width dependent on the temperature of the cell and the mass of the ion (Michel-Dansac et al. 2020) and further modulated based on the bulk gas velocity of the cell.

When modelling the emission lines from gas cells we use $10^6$ photon packets for the entire galaxy (distributed by randomly assigning photon packets to cells weighted by luminosity) while for the continuum, resonant regions of the spectrum are sampled with $5 \times 10^6$ photon packets while non-resonant regions are sampled with $10^7$ photon packets. These values were selected to be the maximum number of photons feasible for the computational time and storage space we had available for this project.

After the initial conditions are defined, we propagate the photon packets through different media out to the virial radius of the halo. For non-resonant lines (e.g. [O III] 5007Å), we ignore any self-absorption and consider only intervening dust for absorption and scattering. We employ the phenomenological dust model of Laursen, Sommer-Larsen & Andersen (2009), normalized for the SMC such that the dust absorption coefficient is given by $(n_{H1} + f_{ion}n_{HII})\sigma_{dust}(\lambda)Z/Z_0$. Here $f_{ion} = 0.01$ and $Z_0 = 0.005$. The albedo and dust asymmetry parameters needed to compute the phase function for scattering off of dust are calculated as a function of wavelength from Weingartner & Draine (2001). We note that our choice of dust model is not unique and can impact the results. For example, setting $f_{ion}$ to a higher or lower value is expected to decrease and increase the Mg II escape
fraction, respectively. This scaling is non-linear since the escape is proportional to the exponential of the optical depth. Likewise we have assumed that the effective dust cross-section scales linearly with metallicity where observations (e.g. Rémy-Ruyer et al. 2014) show that this may not be the case. As we do not track dust self-consistently in the simulation, we have adopted this commonly used approximation.

For Ly α, we propagate photon packets through a mixture of hydrogen, deuterium, and dust. The neutral hydrogen density is taken directly from our RAMSES-RT simulation and we assume a fixed D/H abundance of $3 \times 10^{-5}$. Finally, for Mg II, we propagate photon packets through a mixture of Mg II and dust. In order to do this, we must know the Mg II fraction of each cell. The Mg abundance in each cell is calculated assuming GASS solar ratios (Grevesse et al. 2010) and scaled by the metallicity of the cell. In other words, we assume a fixed Mg/H abundance ratio as a function of metallicity and we ignore depletion of Mg on to dust.<fnref n="1" /> We then use CLOUDY to calculate the fraction of Mg in each cell that is in the form of Mg II as a function of density, temperature, metallicity, and local radiation field. Unsurprisingly, the Mg II fraction is most sensitive to temperature. The relation is well fit by a logistic function such that

$$f_{\text{Mg II}} = \frac{1}{1 + e^{-[179.04 - 43.40 \log_{10}(\frac{I}{K})]}} \quad (1)$$

We use this simple approximation to compute the Mg II fraction for each cell in our RASCAS simulations and assume that there is no Mg I present due to photoionization. We note that this approximation may fail in cells in the simulation that have low densities and temperatures (such as in the neutral IGM); however, the cells we consider in this work are all inside the virial radii of haloes and thus the approximation suffices for our purposes.

Finally, we also compute the escape fraction for each line as well as the LyC escape fraction at 912Å. Compared to our previous work (Rosdahl et al. 2018), rather than sampling 500 randomly directed rays from each star particle in the halo, we use $10^7$ photon packets with initial positions sampled from a multinomial, based on the locations and ionizing emissivities of each star particle. Note that the LyC escape fractions in this work refer specifically to those at 912Å, more akin to observational values, in contrast to the LyC escape fractions quoted in Rosdahl et al. (2018) which represent a luminosity-weighted sum across the ionizing spectrum.

For most of this work, the escape fractions and spectra discussed will be the ‘angle-averaged’ (or global) quantities. However, since observations are made along individual sightlines, we use a standard peeling algorithm (Yusef-Zadeh, Morris & White 1984; Zheng & Miralda-Escudé 2002; Dijkstra 2017) as implemented in RASCAS to compute images and spectra along the three principal axes of the cosmological simulation. Furthermore, we have also computed the escape fractions along the same lines of sight. We will explicitly indicate when the results along individual sightlines are used.

## 3 RESULTS

In this section, we will analyse whether Mg II can be used as an indicator for $f_{\text{esc}}$ in the epoch of reionization by studying the

emission line properties of 694 galaxies with $M_{\text{halo}} > 3 \times 10^9 M_{\odot}$ in SPHINX$^{20}$ at $z=6$. This mass selection corresponds to stellar masses of $\sim 10^7 M_{\odot} - 10^{10} M_{\odot}$ and was implemented for numerous reasons. Resonant line radiative transfer is sensitive to how well the gas distribution of a galaxy is resolved. By limiting our analysis to higher mass galaxies, the impact of resolution will be less severe compared to lower mass systems. Furthermore, the higher mass galaxies are much more metal enriched (see Katz et al. 2022) which means that the likelihood of being observable scales with halo mass. However, by placing a mass cut at $M_{\text{halo}} > 3 \times 10^9 M_{\odot}$, we are excluding a significant number of galaxies in SPHINX$^{20}$ that contribute to the reionization of the volume.

### 3.1 General properties

Before determining whether Mg II emission can be used as a probe of the LyC escape fraction, it is first important to determine how the properties of the simulated galaxies compare with the samples of lower redshift objects that are thought to be analogues of galaxies in the epoch of reionization. More specifically, we compare with Green Peas, Blueberries, as well as with suspected analogues at $z \sim 3$. Furthermore, we show specific comparisons with known LyC leakers at low redshift.

We emphasize here that our goal is not to model low-redshift analogues. Rather our simulations are designed to model $z \ge 6$ galaxies from first principles. This comparison only serves to highlight any biases that may arise from attempting to apply empirical indirect relations for the LyC escape fraction derived on low-redshift galaxies to simulated $z \ge 6$ galaxies. Observed low-redshift galaxies may not be perfect analogues of high-redshift systems and similarly our simulated galaxies are unlikely to be perfect representations. Nevertheless, by comparing these two complementary approaches and understanding the differences, we aim to improve our ability to elucidate the physics of reionization epoch galaxies from future JWST observations.

#### 3.1.1 Comparison of SPHINX$^{20}$ galaxies with green peas, blueberries, $z \sim 3$ analogues, and $z \sim 1$ MUSE galaxies

In Fig. 1 we compare numerous properties of low-redshift Green Pea galaxies as compiled by Yang et al. (2017a), Blueberry galaxies<fnref n="2" /> Yang et al. (2017b), $z = 2.5-3.5$ high-redshift analogues from Amorín et al. (2017), and $z \sim 1$ MUSE galaxies (Feltre et al. 2018) with SPHINX$^{20}$ galaxies at $z=6$. We emphasize here that the goal of this exercise is not to demonstrate that SPHINX$^{20}$ galaxies are direct analogues of Green Peas or Blueberries (or any other low-redshift galaxy population) as there is no reason that they should be. Rather, the purpose is to highlight and understand any differences between the galaxy populations which should be kept in mind when extrapolating our Mg II results to low-redshift.

In general, Green Pea galaxies exhibit SFRs towards the higher end of the SPHINX$^{20}$ sample at fixed stellar mass (i.e. higher specific star formation rates, sSFRs). This is particularly true at stellar masses of $M_* < 10^9 M_{\odot}$. Similarly we find that in this same mass interval, Green Pea galaxies fall towards the higher end of the metallicity distribution of $z=6$ SPHINX$^{20}$ galaxies. Blueberry galaxies have lower SFRs than Green Peas at fixed stellar mass which are more comparable to what we find in SPHINX$^{20}$. At the lowest stellar masses (i.e. $M_* \sim 10^7 M_{\odot}$), the SFRs of Blueberry galaxies tend to be higher than

<footnote n="1">¹Similar to our discussion on choice of dust model, adopting different depletion fractions can impact the Mg II escape fraction non-linearly. Since the simulated galaxies are at low metallicity, the impact is expected to be mild, especially for the lower stellar mass objects.</footnote>

<footnote n="2">²Blueberry galaxies are the $z < 0.05$ counterparts of Green Peas.</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1.** Comparison of $z = 6$ SPHINX$^{20}$ galaxies (red circles) with low-redshift Green Peas from Yang et al. (2017a) (green pentagons), $z < 0.05$ blueberry galaxies from Yang et al. (2017b) (blue circles), $z = 2.5–3.5$ high-redshift analogues from Amorín et al. (2017) (downward black triangles), and $z \sim 1$ MUSE galaxies (Feltre et al. 2018) (beige triangles). From the top left, we show SFR, metallicity ($12 + \log (O/H)$) as calculated in the H II regions of SPHINX$^{20}$ galaxies, colour excess ($E(B-V)$), Ly α luminosity, Ly α escape fraction, H α luminosity, H α equivalent width, H β luminosity, [O III] 5007Å luminosity, [O III] 4959Å luminosity, [O II] 3727Å doublet luminosity, and O32 versus the stellar mass of the galaxy. For all luminosities, we compare the dust attenuated ('escaped') value with observations while for O32, we compare the intrinsic dust corrected value. Colour excess values have been computed in two ways, either assuming an intrinsic ratio of 2.86 as is common in the literature (red) or using the true intrinsic H α/H β ratio in the simulation (beige). For additional comparison, we show [O III] 5007Å luminosities for intense star-forming galaxies from the MOSEL sample at $z \sim 3–4$ (Tran et al. 2020). In general, Green Pea galaxies exhibit SFRs at the higher end of the SPHINX$^{20}$ sample at fixed stellar mass (i.e. higher sSFR) and have significantly higher nebular oxygen emission. Blueberry galaxies have more similar properties to $z = 6$ SPHINX$^{20}$ galaxies in nearly all aspects considered. Note that observational data are only included on each panel when available.
SPHINX$^{20}$ galaxies. Furthermore, the Blueberry galaxy population seems to follow the same stellar mass–metallicity trend as Green Peas, which overlaps with SPHINX$^{20}$ galaxies at $M_{*} \gtrsim 10^{7.5}$ M$_{\odot}$.

Dust content can be probed by measuring the colour excess ($E(B-V)$). Green pea galaxies exhibit a significant amount of scatter in this quantity, independent of stellar mass; however, most galaxies have $E(B-V) \lesssim 0.2$. $E(B-V)$ was computed for Green Pea galaxies assuming a Calzetti et al. (2000) dust extinction law and an intrinsic $H\alpha/H\beta$ ratio of 2.86 which is appropriate for Case B recombination at a temperature of $10^4$ K (Osterbrock 1989). Since the temperatures in the ISM of SPHINX$^{20}$ galaxies are not all $10^4$ K and similarly, since we also include the collisional contribution, we compute $E(B-V)$ in two ways to compare with observations. The first method (denoted as ‘True’) uses the true intrinsic $H\alpha/H\beta$ ratio (shown in beige), while the second method ignores the true intrinsic ratio and assumes it to be 2.86 (shown in red and denoted as 2.86). The colour excess measurements using the true method are generally lower than the 2.86 method because the intrinsic values tend to be higher. Similar to Green Peas, SPHINX$^{20}$ galaxies exhibit a large scatter in colour excess and many of the simulated galaxies scatter to $E(B-V)$ values that are nearly twice as high as observations.

Green Pea galaxies are categorized by their extreme emission line equivalent widths. Beginning with $Ly\alpha$, in Fig. 1 we see that at stellar masses below $10^9$ M$_{\odot}$ Green Peas tend to have higher $Ly\alpha$ luminosities than SPHINX$^{20}$ galaxies. At higher stellar masses, the less luminous Green Pea $Ly\alpha$ emitters are more similar to the simulated systems. Although there is considerable scatter, the $Ly\alpha$ escape fractions of SPHINX$^{20}$ galaxies tend to decrease with increasing stellar mass. The same holds true for the sample of ~40 galaxies in the Yang et al. (2017a) catalogue. However, Green Peas at $M_{*} > 10^9$ M$_{\odot}$ exhibit $Ly\alpha$ escape fractions up to ~20 per cent while these values tend to be lower than ~10 per cent in SPHINX$^{20}$.

Similar to $Ly\alpha$, the $H\alpha$ and $H\beta$ luminosities of Green Peas also fall towards the higher end of the SPHINX$^{20}$ distribution. Blueberry galaxies have lower $H\beta$ luminosities at fixed stellar mass compared to Green Peas which is more consistent with the $H\beta$ luminosities of SPHINX$^{20}$ galaxies. In terms of equivalent widths (EWs), Green Peas tend to have higher $H\alpha$ EWs than SPHINX$^{20}$ galaxies at $M_{*} \gtrsim 10^{8.5}$ M$_{\odot}$ while at lower stellar masses, the EWs are comparable.

The ‘green’ aspect of Green Peas is due to intense oxygen emission lines that appear in the *r* band of SDSS images. In Fig. 1, we compare the line luminosities of [O III] 5007Å, [O III] 4959Å, and [O II] 3727Å doublet with SPHINX$^{20}$ galaxies. We see a general trend that the oxygen emission lines from these galaxies are nearly an order of magnitude brighter in Green Peas than in the simulated galaxies. Such extreme oxygen line luminosities are also characteristic of *z* ~ 3–4 galaxies that sit well above the star formation main sequence (e.g. Tran et al. 2020). Blueberry galaxies have also been selected to have intense oxygen emission lines; however, similar to $H\beta$, Blueberries, have lower oxygen emission line luminosities at fixed stellar mass compared to Green Peas. The oxygen line luminosities still fall towards the upper end or even above the emission line luminosities found in SPHINX$^{20}$; however, we highlight that the oxygen abundance ratios used in this work are Solar (Grevesse et al. 2010) values scaled by the metallicity of the gas and deviations from non-solar ratios, such as those described in Katz et al. (2022) may enhance the oxygen emission line luminosities. We note that both Green Pea and Blueberry galaxies are star bursting systems, while this is not necessarily the case for all SPHINX$^{20}$ galaxies. Hence many of

the SPHINX$^{20}$ galaxies exhibit SFRs that are multiple orders of magnitude lower than what is measured for the Green Pea and Blueberry populations.

Despite the weaker oxygen emission line luminosities, the O32 ratios for Green Pea and Blueberry galaxies are very similar to those of SPHINX$^{20}$ galaxies. At $M_{*} > 10^9$ M$_{\odot}$, Green Pea galaxies seem to have significantly lower O32 ratios compared to the simulated galaxies. The observed trend points towards mildly decreasing O32 with increasing stellar mass while SPHINX$^{20}$ galaxies exhibit a flat or slightly increasing trend of O32 with stellar mass. We note that there are only a few galaxies in SPHINX$^{20}$ with $M_{*} > 10^9$ M$_{\odot}$ so a larger sample of galaxies at higher stellar masses will be needed to confirm this trend. Nevertheless, the O32 ratios of Green Peas, Blueberries, and SPHINX$^{20}$ galaxies are all indicative of high ionization parameters compared to more typical galaxies at lower redshift.

In summary, Green Peas and our high-redshift simulated galaxies are not perfect analogues of each other. Green Peas tend to exhibit higher sSFRs as well as brighter hydrogen and oxygen emission line luminosities. The SPHINX$^{20}$ galaxies are more akin to Blueberries, the *z* < 0.05 counterparts of Green Pea galaxies, in almost all metrics considered.

### 3.1.2 Comparison with known low-redshift LyC leakers

As our primary aim is to understand the utility of Mg II as a probe of LyC escape, we continue our comparison by showing properties of the simulated high-redshift galaxies with a sample known low-redshift LyC leakers that have Mg II detections. Data for these galaxies was compiled from Izotov et al. (2016a,b, 2018a,b, 2021), Guseva et al. (2020).

In Fig. 2 we show stellar mass, SFR, sSFR, metallicity (12 + log (O/H)) as calculated in the H II regions of SPHINX$^{20}$ galaxies, $H\alpha$ luminosity, $H\beta$ luminosity, [O II] 3727Å doublet luminosity, [O III] 5007Å luminosity, and O32 versus $f_{\text{esc}}$ for SPHINX$^{20}$ galaxies and local LyC leakers. A detailed study on how many of these properties correlate with $f_{\text{esc}}$ will be presented in Rosdahl et al. (in preparation) and here we aim to highlight the similarities and differences between simulated $z=6$ galaxies and low-redshift leakers.

While low-redshift leakers overlap in stellar mass (albeit on the higher end) of our simulated galaxy population, low-redshift leakers tend to exhibit higher metallicities as well as significantly higher SFRs and sSFRs. Compared to the bulk of the low-redshift galaxy population, candidate leakers are often selected based on extreme emission line properties (possibly driven by the high sSFRs) and these galaxies are also more extreme than the simulated galaxy population. This can be observed in the plots in Fig. 2 that show the luminosities of $H\alpha$, $H\beta$, [O II] 3727Å and [O III] 5007Å against the LyC escape fraction. Local leakers exhibit significantly higher emission line strengths compared to the simulated $z=6$ leakers. Interestingly, the O32 emission line ratios are similar (although perhaps slightly higher) for low-redshift leakers compared to the simulated galaxies. This will be important later as this ratio has been used to predict the intrinsic Mg II luminosity.

We emphasize again that the SPHINX$^{20}$ simulations were designed to predict the properties of $z=6$ galaxies (both LyC leakers and non-leakers) rather than low-redshift LyC leakers. However, since indirect probes of LyC leakage are often developed empirically from observations, as is the case for Mg II, it is important to understand the differences between the two galaxy populations. With these caveats
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of $z = 6$ SPHINX$^{20}$ galaxies (red circles) with low-redshift LyC leakers that have Mg II detections (dark blue circles) from Izotov et al. (2016a,b, 2018a,b, 2021), Guseva et al. (2020). Local galaxies with only upper limits on $f_{esc}$ are indicated by left-pointing triangles. From the top left, we show stellar mass, SFR, sSFR, metallicity ($12 + \log(O/H)$) as calculated in the H II regions of SPHINX$^{20}$ galaxies, H$\alpha$ luminosity, H$\beta$ luminosity, [O II] 3727Å doublet luminosity, [O III] 5007Å luminosity, and O32 versus $f_{esc}$.</figcaption></figure>

in mind, we consider how SPHINX$^{20}$ galaxies compare with the low-redshift Mg II emitters.

### 3.1.3 Comparison of Mg II emission

In Fig. 3 we plot the Mg II luminosity (i.e. sum of the doublet emission) as a function of stellar mass for SPHINX$^{20}$ galaxies at $z = 6$ compared with observations (Rubin et al. 2011; Martin et al. 2013; Izotov et al. 2016a,b, 2018a,b, 2021; Feltre et al. 2018; Henry et al. 2018; Guseva et al. 2020; Zabl et al. 2021) as well as other cosmological simulations (Nelson et al. 2021). We show both the intrinsic emission (red) as well as the escaped flux (green). Consistent with our earlier results, we find that the Mg II emission from low-redshift Green Pea galaxies and low-redshift LyC leakers is significantly higher at fixed stellar mass compared to the predicted Mg II luminosities of $z = 6$ galaxies from SPHINX$^{20}$. The Green Peas/local leakers exhibit higher Mg II luminosities compared to other observed low-redshift galaxies that have significantly higher stellar masses. These other observations exist at higher stellar masses

than are probed by SPHINX$^{20}$; however, they seem to follow the general trend of Mg II luminosity with stellar mass that we see in our simulation.

The $z \sim 1$ galaxy sample from MUSE (Feltre et al. 2018) tends to have Mg II luminosities closer to SPHINX$^{20}$ galaxies than do Green Peas. At higher stellar masses, SPHINX$^{20}$ galaxies tend to exhibit higher Mg II luminosities however the lowest stellar mass MUSE Mg II emitter has an Mg II luminosity at the upper end of what is exhibited by SPHINX$^{20}$ galaxies. This is likely due to the fact that the low stellar mass galaxies from MUSE tend to have high SFRs compared to the galaxy main sequence while the higher stellar mass galaxies fall either on or below the $z \sim 1-2$ galaxy main sequence (Whitaker et al. 2014).

Compared with IllustrisTNG galaxies at $z = 2$ (Nelson et al. 2021), SPHINX$^{20}$ galaxies at $z = 6$ exhibit approximately an order of magnitude higher Mg II luminosity. We emphasize that there is no reason that galaxies at fixed stellar mass should have the same Mg II luminosities due to the different redshift probed. Furthermore, simulation differences such as mass and spatial resolution, feedback
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Stellar mass versus Mg II luminosity for SPHINX²⁰ galaxies compared to observations at $z < 1$. We show both the intrinsic (red) and escaped (light green) Mg II luminosity for our $z = 6$ galaxies. The solid lines of the same colours represent the running medians. The blue band represents the $z = 2$ simulation results from Nelson et al. (2021). Mg II emission from Green Pea galaxies from Henry et al. (2018) are plotted as green pentagons, for low-redshift LyC leakers from Izotov et al. (2016a,b, 2018a,b, 2021), Guseva et al. (2020) as black pentagons, for MUSE galaxies at $z \sim 1$ from Feltre et al. (2018) are shown as beige triangles, and various other low-redshift observations (Rubin et al. 2011; Martin et al. 2013; Zabl et al. 2021) are shown in black. Empty symbols show observations where only one of the two Mg II lines was detected.</figcaption></figure>

models, and metal enrichment strategies are expected to cause
differences in the Mg II luminosity and a more systematic study is
needed to understand this effect. Note that the luminosities reported
in Nelson et al. (2021) represent the intrinsic luminosities and thus
should be compared with the red data points in Fig. 3. Interestingly,
Nelson et al. (2021) finds very little evolution in the relation between
Mg II luminosity and stellar mass between $z = 2$ and $z = 0$ and their
predictions are in reasonable agreement with MUSE galaxies. If
SPHINX²⁰ and IllustrisTNG are consistent, we would expect certain
ISM physics to change between $z = 6$ and $z = 2$ that would drive
this difference. This may include an evolution in the galaxy main
sequence as well as other physics such as depletion of Mg on to
dust.

### 3.1.4 Observability of Mg II at $z \ge 6$ with JWST

In the left-hand panel of Fig. 4 we show the Mg II 2796Å and Mg II 2803Å luminosity functions at $z = 6$ for SPHINX²⁰ galaxies with halo masses $\ge 3 \times 10^9 M_\odot$. The intrinsic luminosities reach values $> 10^{41} \text{ erg s}^{-1}$. A $10\sigma$ detection of the brightest Mg II emitters in SPHINX²⁰ at $z = 6$ would require $O(10^5 \text{ s})$ integration times with JWST NIRSpec.

One important question is whether SPHINX²⁰ galaxies are rep-
resentative of the high-redshift galaxy population. In the centre
panel of Fig. 4 we compare the Ly α luminosity function at z =
6 with observational constraints from Konno et al. (2018), Drake
et al. (2017) at a similar redshift. We show both the intrinsic
Ly α luminosity for SPHINX²⁰ galaxies (dashed green line) as well
as the values after attenuation by the ISM and CGM. Note that we
do not consider attenuation by the IGM (see e.g. Garel et al. 2021

for the importance of this effect in SPHINX). While the volume
of SPHINX²⁰ is not large enough to capture the bright end of the
luminosity function where the majority of the observations probe,
we find that at luminosities of ~10⁴² erg s⁻¹, the predictions from
SPHINX²⁰ are in good agreement with observations.

The 1500Å UV luminosity function (UVLF) has significantly better constraints at faint magnitude at $z = 6$ compared to the Ly α luminosity function. In the right-hand panel of Fig. 4 we compare the results from SPHINX²⁰ to observations. Similar to Ly α, we find that the intrinsic UVLF (i.e. not accounting for attenuation) overpredicts the number of bright galaxies compared to observations; however, when dust attenuation is accounted for, the simulations are in much better agreement with the observational data (see also Garel et al. 2021). Note that the turnover at faint UV magnitudes is due to only considering haloes with $M_{\text{halo}} > 3 \times 10^9 M_\odot$ in the computation of the UVLF. While these results do not guarantee that our Mg II predictions are correct, they provide confidence that the model employed in SPHINX²⁰ can produce certain realistic aspects of galaxy populations at high-redshift.

## 3.2 Mg II emitters versus absorbers

The spectral profiles of Mg II emission can provide significant insight into various galaxy properties. In general, pure emission would represent a relatively optically thin line of sight, pure absorption indicates a high column density medium, and a P-Cygni or inverse P-Cygni profile can signify outflows and inflows as well as other complex kinematics.

Feltre et al. (2018) demonstrated that among a sample of 123 galaxies with Mg II detections at $0.7 \le z \le 2.34$, about 50 per cent were emitters while the others either exhibited P-Cygni profiles or were Mg II absorbers. No particular redshift dependence was found for the fraction of emitters versus absorbers in Feltre et al. (2018). In order to determine whether Mg II emission is a useful probe of the LyC escape fraction, we must first determine whether high-redshift galaxies are Mg II emitters and if there are any biases in the galaxy population of emitters versus non-emitters.

### 3.2.1 Angle-averaged spectra

For all SPHINX²⁰ galaxies at $z = 6$, we have combined the emission spectra from the gas with the processed continuum radiation from star particles and manually classified the spectra as either a Mg II emitter, a Mg II absorber, as having a P-Cygni profile, or as either having flat or arbitrarily complex spectra. We stress that this classification is qualitative and there is subjectivity in the classification. We further emphasize that we do not consider whether the galaxy is likely to be detected (or the significance of the detection as some emitters or absorbers may be very weak).

We first consider the angle-averaged spectra. In Table 2 we list the number of SPHINX²⁰ galaxies at $z = 6$ as well as the percentage of the total sample that fall into each of the four different classifications. Most of the SPHINX²⁰ galaxies (64 per cent) are emitters, consistent with low-redshift observations (Feltre et al. 2018).<fnref n="3" /> In contrast to the Feltre et al. (2018) sample, we find twice the number of P-Cygni galaxies and half the number of absorbers. Our P-Cygni galaxies often have weaker emission and absorption compared to either the

<footnote n="3">³Note that there are spectral resolution differences between the simulations and observations that can also determine which spectral features are observable. For the simulations we have computed the spectra at 1Å resolution.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Mg II (left-hand panel), Ly α (centre), and 1500Å (right-hand panel) luminosity functions for SPHINX20 galaxies at z = 6. The intrinsic luminosity functions are shown as the dashed lines while the solid lines represent the luminosity functions after attenuation by the ISM and CGM, hence the observed values. The shaded region around our predictions represent the 1σ Poisson uncertainty. Note that the turnover in the luminosity function at the faint end is due to only considering haloes with $M_{\text{halo}} > 3 \times 10^9 M_{\odot}$ (this is also true for the Ly α and Mg II luminosity functions). For Mg II we show results for both resonant lines as indicated in the legend. We compare the Ly α luminosity functions to the z ~ 6 results from Konno et al. (2018), Drake et al. (2017) while the UV luminosity functions are compared with observational constraints from Atek et al. (2018), Finkelstein et al. (2015), Livermore, Finkelstein & Lotz (2017), Bouwens et al. (2017).</figcaption></figure>

<table><caption>**Table 2.** Class distribution of the angle-averaged and line of sight spectra for SPHINX20 galaxies at z = 6. Galaxies are classified as being either emitters, absorbers, having a P-Cygni profile, or as N/A when the spectrum is either flat or arbitrarily complicated. The final column lists the distribution seen at z ~ 1 with MUSE (both the number of galaxies and the percentage of the total) from Feltre et al. (2018) where N/A here corresponds to non-detections.</caption><thead><tr><th>Class</th><th># of Galaxies</th><th>Per cent of sample</th><th>MUSE</th></tr></thead><tbody><tr><td><strong>Angle averaged</strong></td><td></td><td></td><td></td></tr><tr><td>Emitters</td><td>446</td><td>64 per cent</td><td>63 (16 per cent)</td></tr><tr><td>P-Cygni</td><td>207</td><td>30 per cent</td><td>19 (5 per cent)</td></tr><tr><td>Absorbers</td><td>10</td><td>1 per cent</td><td>41 (11 per cent)</td></tr><tr><td>N/A</td><td>31</td><td>5 per cent</td><td>258 (68 per cent)</td></tr><tr><td><strong>Line of sight</strong></td><td></td><td></td><td></td></tr><tr><td>Emitters</td><td>823</td><td>40 per cent</td><td>63 (16 per cent)</td></tr><tr><td>P-Cygni</td><td>1175</td><td>56 per cent</td><td>19 (5 per cent)</td></tr><tr><td>Absorbers</td><td>83</td><td>4 per cent</td><td>41 (11 per cent)</td></tr><tr><td>N/A</td><td>1</td><td>0 per cent</td><td>258 (68 per cent)</td></tr></tbody></table>

emitters or absorbers which may make it hard to achieve a significant detection without long integration times. We once again highlight the subjectivity in the classification as we find that many of the emitters and absorbers exhibit very weak signatures of P-Cygni profiles as well. Furthermore, we note that our mass selection of haloes to analyse is completely different for how galaxies were selected in Feltre et al. (2018) which is likely responsible for the differences in Mg II characteristics.

In Fig. 5 we show five example spectra that fall into each of the four classes. While many emitters exhibit spectra that have nearly Gaussian profiles (as was seen in Chisholm et al. 2020), the second emitter in Fig. 5 shows a double peaked profile while the third exhibits a very weak P-Cygni profile. Mg II absorbers also have very diverse spectra. The first absorber in Fig. 5 has very deep and broad absorption in the continuum and the emission from the gas fills in the absorption to some extent, while the last absorber is much weaker. The third absorber exhibits a complicated profile with very weak emission. Among the five P-Cygni spectra, we show three examples of classic P-Cygni profiles where the absorption

is blueshifted compared to the emission, usually interpreted as a signature of outflows, as well as two inverse P-Cygni spectra where the absorption is redshifted, potentially a signature of inflows. P-Cygni galaxies are often dominated by the continuum; however, we do find examples where the collisional emission from gas is also a significant component of the spectra.

In Fig. 6 we show examples of flux maps for the line emission and continuum for a Mg II emitter, absorber, and P-Cygni galaxy. The three galaxies correspond to those in the top row of Fig. 5. Here one can see that the different galaxies have extremely different morphologies. The Mg II emitter exhibits very extended emission (see Burchett et al. 2021; Zabl et al. 2021 for examples of extended Mg II emission in observations) whereas the flux from the absorber is very centrally concentrated. Flux from the P-Cygni galaxy is dominated by the continuum and has a very clumpy morphology, consistent with the locations of star forming regions and satellite galaxies. The spectral shape results from the complex kinematics of these different components.

In the top row of Fig. 7, we analyse the differences between the galaxy populations that are Mg II emitters versus absorbers. As stellar mass, SFR, and LyC escape fraction increases, so does the percentage of galaxies that are Mg II emitters. In contrast, the fraction of Mg II absorbers and P-Cygni galaxies decrease as these other quantities increase. At $M_* > 10^8 M_{\odot}$, SFR > 0.1 $M_{\odot}yr^{-1}$, and $f_{\text{esc, LyC}} > 0.01$, more than 60 per cent of galaxies are emitters. Because the emitter fraction increases with stellar mass and also with $f_{\text{esc, LyC}}$, this may give the impression that $f_{\text{esc, LyC}}$ increases with increasing stellar mass. However, the total number of galaxies in SPHINX20 is dominated by low stellar mass objects. We find that at fixed stellar mass, galaxies with higher escape fractions tend to be biased towards being Mg II emitters, while low mass galaxies tend to exhibit the highest $f_{\text{esc, LyC}}$.

Interestingly, our trend differs from that of Feltre et al. (2018) who found that Mg II absorbers tend to have higher stellar masses than P-Cygni galaxies, which, in turn, have higher stellar masses than Mg II emitters. The same holds true for SFR. The underlying cause of this difference is not entirely clear as we find that many of the low stellar mass galaxies in our simulations are dominated by the continuum. It is important to consider that the vast majority of galaxies in the Feltre et al. (2018) sample are non-detections.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Example spectra of Mg II emitters (first column), Mg II absorbers (second column), P-Cygni galaxies (third column), and either flat or arbitrarily complicated galaxies (fourth column). All spectra have been normalized to 1.1$I_{\lambda, \max}$. Green, beige, and magenta lines show the processed spectra for the gas, stellar continuum radiation, and the total.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Flux maps of an example Mg II emitter (left-hand panel), absorber (centre), and P-Cygni galaxy (right-hand panel). The top row shows the flux from the line emission while the bottom row shows the continuum. The galaxies in each column correspond to the spectra in the first row of Fig. 5. The width of each image spans the virial diameter of each halo and a 1 physical kpc scale-bar is plotted for each halo.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fraction of SPHINX²⁰ galaxies at a given stellar mass (left-hand panel), SFR (centre), and LyC escape fraction (right-hand panel) that are Mg II emitters, absorbers, or P-Cygni galaxies for both the angle averaged spectra (top row) as well as line of sight spectra (bottom row). The window size for all quantities is set to 0.4 dex. The shaded region represents the 1σ scatter about the relation. Note that the crossover between P-Cygni and emitters at the highest escape fractions in the angle average plot is due to only one galaxy in the highest escape fraction bin being a P-Cygni galaxy.</figcaption></figure>

Furthermore, at lower stellar masses, many of the galaxies populate the regime above the star formation main-sequence, while this is less true for the higher stellar mass galaxies. However, if our simulations are an adequate representation of the high-redshift Universe, selecting galaxies based on their location with respect to the star formation main sequence may not produce a bias in terms of the number of each class of galaxy. This is because in Fig. 8, we show that all different spectral types populate the stellar mass-SFR plane approximately equally. We also note that the stellar masses of SPHINX²⁰ galaxies tend to be at the lower end of the Feltre et al. (2018) galaxy sample. In the stellar mass regime where the two samples overlap, Mg II emitters represent the dominant population. Thus, it is possible that if we simulated a larger volume and probed more massive haloes, the simulated trend may reverse and decrease at higher stellar masses to be more consistent with Feltre et al. (2018).

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Star formation main sequence (i.e. SFR versus stellar mass) for SPHINX²⁰ galaxies (circles) at z = 6 coloured by their spectral class. Observational results at z = 6 from Pearson et al. (2018) and at z = 5.5 from ALPINE (Schaerer et al. 2020) are shown as the black and pink dashed lines while lower redshift galaxies from Feltre et al. (2018) are shown as stars. Solid lines show the running median SFRs as a function of stellar mass for each spectral class. The colours of stars have the same meaning as the solid coloured lines. We find no difference in the star formation main sequence for the different spectral classes.</figcaption></figure>

### 3.2.2 1D Spectra

We have performed a similar exercise on the 1D, line of sight Mg II spectra and manually classified three different directions (i.e. 2,082 spectra) as either emitters, absorbers, as containing both strong emission and absorption features, or as having a rather flat spectrum, for each SPHINX²⁰ halo with $M_{vir} > 3 \times 10^9 M_{\odot}$ at z = 6. For this experiment, we have switched the third class from P-Cygni as many of the 1D spectra have both strong emission and absorption that are not typical of the standard P-Cygni profile.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Example spectra for the same galaxy along three different sightlines. Depending on the viewing angle, Mg II may be observed in emission, absorption, or in both emission and absorption.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Images of Mg II emission for the three spectra shown in Fig. 9. The top row shows the emission from the gas while the bottom row shows the stellar continuum.</figcaption></figure>

In contrast to the angle-averaged spectra, we find that for the 1D spectra, 56 per cent exhibit both strong emission and absorption features, including P-Cygni profiles. 40 per cent of the spectra are Mg II emitters and 4 per cent are Mg II absorbers. This reflects the fact that the structure of the ISM and CGM gas in our simulations is complex and it is not uncommon to have both optically thin and thick channels along the same line of sight.

In Fig. 9 we show an example galaxy where depending on the sightline, the galaxy may be observed as an Mg II emitter, absorber, or as having both emission and absorption. Images of the Mg II emission corresponding to the spectra are shown in Fig. 10. This complicates the use of Mg II as a diagnostic or constraint for individual galaxy properties, such as the LyC escape fraction which is also known to have a significant directional dependence (e.g Cen & Kimm 2015). Large samples of Mg II emitters will be necessary to make inferences about the galaxy population in aggregate.

In the bottom row of Fig. 7 we show the class fraction distribution of the line of sight spectra as a function of stellar mass, SFR, and line of sight LyC escape fraction. The trends along the line of sight are very similar to those seen in the angle averaged spectra. There is a particularly strong dependence of the emitter fraction with the line of sight LyC escape fraction. Interestingly, the emitter fraction peaks at a stellar mass of $10^{8.5}$ M$_{\odot}$ before dropping at higher stellar

mass. This may help reconcile our data with the observations from Feltre et al. (2018). This peak also corresponds to the stellar mass of maximum angle averaged LyC escape fraction in SPHINX$^{20}$ (see Rosdahl et al., in preparation).

One of the interesting features that we find in our simulations is that in certain cases, the continuum radiation from star particles can be reprocessed to appear like a Mg II emission line along a particular line of sight. This can be observed in the top row of Fig. 9 as bumps in the continuum near the line centres of the emission lines (note that more extreme examples exist). A simple example of how this might happen would be to envision a screen of pure Mg II gas behind a source. Photons that are emitted in a direction opposite from the line of sight that are not close to either of the resonances will pass through the screen unaffected. However, photons near the resonance will be scattered/reflected with some probability back into the direction of the observer. Hence in this simple setup, depending on the optical depth of the screen, we would expect an enhancement in the number of photons an observer sees near the wavelength of the resonance. Such a feature would not appear in the angle-averaged spectra. The strongest emission lines that we see in our 1D data set are still dominated by emission from gas. However, we find examples where reflected continuum emission can both contribute to the strength of an observed emission line for strong emission, or dominate the entire observed emission line for the weaker Mg II emitters.

In summary, according to the angle averaged spectra, we predict that a significant fraction (i.e. 64 per cent) of $z = 6$ galaxies with halo masses $>3 \times 10^9$ M$_{\odot}$ are Mg II emitters, which is a promising result in the context of using Mg II emission as a probe of LyC escape. Along individual sightlines, we expect that this value drops to 40 per cent. However, it is important to consider that Mg II emitters are predicted to have higher stellar masses, SFRs, and LyC escape fractions than non-emitters and thus Mg II emitters are unlikely to be completely representative of the high-redshift galaxy population. When analysing the 1D line of sight spectra, we find that the majority of spectra exhibit both significant absorption and emission features. Furthermore, depending on sightline, the same galaxy may exhibit both emission and absorption which may complicate the interpretation of Mg II for individual galaxies at high-redshift. Finally, we find that along individual sightlines, reflected/scattered continuum radiation can contribute to, and in some cases dominate the strength of emission lines.

## 3.3 Predicting the Mg II escape fraction

The idea behind using Mg II as a probe of the LyC escape fraction is that since Mg II likely traces neutral gas, if one can measure the Mg II escape fraction, then this can be used as a proxy to measure the neutral hydrogen optical depth. In the following subsections, we test whether different methods used in the literature to measure the Mg II escape fraction adequately capture the physics of Mg II escape in our simulated, high-redshift galaxies.

### 3.3.1 The Chisholm model

Chisholm et al. (2020) consider three models to describe the escape of Mg II photons from a galaxy. We briefly summarize them here.

(i) *Model 1 - optically thin slab:* In their optically thin model, they consider a slab of dust-free Mg II gas in front of a source. Photons that interact with the slab are preferentially scattered out of the line of sight of the observer. Because there is a factor of two difference in the oscillator strength between the 2796Å and 2803Å
<figure><img src="image_12.png" /><figcaption>**Figure 11.** (Left-hand panel) The global, or angle-averaged Mg II escape fraction versus the line-of-sight escape fraction for SPHINX²⁰ galaxies at z = 6. For each galaxy, we measure the line-of-sight escape fraction along three different directions. We show the escape fractions for the Mg II 2796Å and Mg II 2803Å lines in green and beige, respectively. The dashed black line shows the one-to-one relation. The absolute scatter increases significantly at high global escape fractions. (Centre) Effective line of sight or (Right-hand panel) global/angle-averaged Mg II escape fraction versus the ratio of the 2796Å and 2803Å lines for SPHINX²⁰ galaxies at z = 6. In the centre panel, the ratio is computed as the ratio of the observed flux along the line of sight while in the right-hand panel, the ratio is computed as the sum of the flux leaking along all directions. We highlight the Mg II emitters from non-emitters using green and beige points, respectively. The red line represents the running median of all the Mg II emitters. Because we compute the line of sight values along three sightlines, the centre panel has three times the number of data points as the right-hand panel.</figcaption></figure>

lines, for the same slab of Mg II gas, the optical depth is twice as large for the 2796Å line as it is for the 2803Å line. Thus the observed ratio, $R = F_{2796Å}/F_{2803Å}$, of Mg II emission should decrease with increasing optical depth such that $R = 2e^{-\tau_{2803}}$.

(ii) *Model 2 - picket fence:* In their second model, they consider Mg II escape through a clumpy geometry that consists of optically thin channels (with zero column density) and optically thick (i.e. $\tau \gg 1$) clouds with some covering fraction. Escape along the optically thin channels preserves the intrinsic ratio whereas Mg II photons that interact with the optically thick clouds are either scattered out of the line of sight or are destroyed and thus do not contribute to the observed flux ratio. Hence in this model $R = 2$.

(iii) *Model 3 - optically thick slab with optically thin holes:* Finally, in their third model, they consider a setup of low column density clouds surrounded by high density regions. In this case, the emerging ratio is $R = 2e^{-\tau_{2803,chan}}$, where the optical depth here represents the optical depth in the low column density channels.

In the latter two scenarios, one must use another method (see Section 3.3.2) to measure the escape fraction rather than the ratio as the optically thick regions do not impact the observed ratio but remove some of the flux, thus lowering the overall escape fraction. If one only uses the ratio for these latter two models the true optical depth will be underestimated. Note that in optically thick slabs, it does not matter if the Mg II photons are absorbed by dust or scattered out of the line of sight by Mg II gas.

For the galaxy analysed in Chisholm et al. (2020), the spectral shape does not show any strong features of scattering or absorption as it appears very Gaussian. This is interpreted to mean that there is a near zero covering fraction of optically thick gas in their galaxy. Hence, the first model should apply and we begin by analysing this model in the current section. For this reason, in the remainder of this section, $R$ values will be computed as if the Mg II emission lines can be perfectly separated from the continuum, thus providing the strongest test of the three models.

It is first important to understand how the line-of-sight escape fraction for Mg II emission relates to the global, angle-averaged value as this will give a sense for how accurately the observed Mg II escape fraction can translate to the global value. In the left-hand panel of Fig. 11 we compare the line-of-sight Mg II escape fraction for

three different viewing angles with the global, angle-averaged value for both Mg II emission lines. In general, galaxies with high angle-averaged Mg II escape fractions also exhibit high line-of-sight escape fractions. However, the absolute scatter increases significantly from nearly zero at low global escape fractions to nearly a factor of 2 at high values while the relative scatter remains approximately constant. Note that the line-of-sight escape fraction can be >1, in contrast to what is usually the case for the LyC escape fraction. This is because photons that were emitted out of an observers line of sight can be scattered into it. In such a scenario, the measured $R$ value from observations for the two Mg II lines would actually be >2. Since the ISM in our simulations has a very complex structure, we may expect that all models considered above and various $R$ values less than and greater than 2 might be possible along different lines of sight.

In the centre panel of Fig. 11 we plot $R$ as a function of the effective line of sight escape fraction<fnref n="4" /> ($\bar{f}_{esc,los}$), defined as

$$ \bar{f}_{esc,los} = \frac{2f_{esc,los,2796} + f_{esc,los,2803}}{3} \quad (2) $$

As $\bar{f}_{esc,los}$ tends towards zero, $R$ approaches 1.9 and for $\bar{f}_{esc,los} > 1$, $R > 2$, consistent with the behaviour expected from the optically thin Chisholm et al. (2020) model. Assuming that this model is correct, we can then derive $\tau_{2803}$ and $\tau_{2796}$ as described above using $R$ and then estimate the Mg II escape fractions for each line as $e^{-\tau}$.

In the top panel of Fig. 12 we show the predicted (i.e. that derived from the Chisholm et al. 2020 model) $\bar{f}_{esc,los}$ versus the value calculated from the Monte Carlo radiative transfer simulations. It is immediately obvious that the optically thin model does not predict the correct escape fractions for our simulated galaxies as the trend does not align with the one-to-one relation. Similar to the simple slab model, the ‘picket fence’ is also not a good representation of the data because the $R$ values are not all 2.

The third Chisholm et al. (2020) model is a potential solution in that the $R$ value is set by the optically thin channels while the total

<footnote n="4">⁴We have chosen to use this effective value as it incorporates information from both lines and it represents an intrinsic flux-weighted average of the escape fractions of the two lines. Our results are not significantly different if we choose to use the escape fraction of one of the two Mg II lines.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** (Top) Effective line-of-sight Mg II escape fraction as measured in our simulation versus the predicted value from the optically thin Chisholm et al. (2020) model. The black dashed line shows the one-to-one relation. (Bottom) Estimated Mg II covering fraction assuming the third model of Chisholm et al. (2020).</figcaption></figure>

Mg II emission is substantially reduced by optically thick clouds with some covering fraction <1. Assuming that the optical depth in the thick clouds is ≫1, the covering fraction can be estimated as $1 - \frac{f_{esc, true}}{f_{esc, pred}}$, where $f_{esc, pred}$ is the value predicted by the optically thin model and $f_{esc, true}$ is the value from the simulation. Note that in real observations, $f_{esc, true}$ is unknown and one would need another method to derive this before the covering fraction can be calculated. In the bottom panel of Fig. 12 we show what the covering fractions as a function of the effective line-of-sight Mg II escape fraction would have to be (in the context of the Chisholm et al. 2020 models) in order to recover the correct intrinsic R and Mg II escape fraction.

For the lowest $\bar{f}_{esc,los}$, the covering fractions approach 1 which is needed to reconcile the differences between the $\bar{f}_{esc,los}$ predicted by the optically thin model and the true values from the simulation. We note that the covering fractions go below 0 for $\bar{f}_{esc,los} > 1$ which is, once again, indicative of photons being scattered into the line of sight.

In summary, our experiments show evidence against the first and second models presented in Chisholm et al. (2020) because the optically thin model does not predict the correct line-of-sight escape fractions, nor are our R values all 2. The question remains whether the third model is an accurate representation of the ISM of real galaxies.

We argue that the third model is also likely inadequate for a few reasons. First, because the model is a slab or screen of low column density channels surrounded by high column density regions, it clearly cannot accurately model the galaxies where $\bar{f}_{esc,los} > 1$. If we simply modify the geometry by considering, for example, a shell of low column density channels surrounded by high column density regions, this issue can be remedied. However, for this modified geometry, if we look at the global or angle-averaged value of the Mg II escape fraction, we would expect the value to be one. Since we are considering resonant lines, for a pure Mg II gas, the photons will never be destroyed, rather for optically thick media, they will scatter into the wings of the line profile. Hence, in this case, any Mg II photon that is input into the simulation should escape the galaxy with a different frequency. We also highlight that for any 3D symmetric distribution of matter, as long as the holes in the distribution are small compared to the field of view, the line of sight Mg II escape fraction will also be 1. Thus as long as this condition holds, our arguments above apply to both global and line of sight escape fractions.

In the left-hand panel of Fig. 11, it is clear that the global Mg II escape fractions are not all 1. This is because our simulation is not pure Mg II gas, rather the simulation also contains dust which can both scatter and more importantly absorb the photons, thus reducing the escape fraction. The presence of dust can impact the R value because the shorter wavelength Mg II photons spend more time scattering, thus increasing the probability that they are absorbed by dust.

This effect is further explored in the right-hand panel of Fig. 11 where we plot the angle-averaged effective Mg II escape fraction versus the angle-averaged R value. The galaxies that have high $\bar{f}_{esc}$ nearly all exhibit R ~ 2 while galaxies that have low $\bar{f}_{esc}$ have a range in R from ~1.7–2. In the context of our interpretation, the galaxies with low $\bar{f}_{esc}$ must have ISMs where enough scatterings occur that the cumulative probability of being absorbed by dust is not insignificant, whereas for the galaxies that have R ~ 2, dust has a far less significant role.

We highlight the fact that there is a bias in R values between Mg II emitters and non-emitters. In the right-hand panel of Fig. 11 we have separated Mg II emitters from non-emitters as green and beige points, respectively. Most of the galaxies at low $\bar{f}_{esc}$ are non-emitters and all but two emitters in our galaxy sample at z = 6 have R > 1.9. The same holds true when analysing the spectra along individual lines of sight where we find a mean value of R = 2.0 for emitters and R = 1.95 for non-emitters. The spectra with the lowest R values are preferentially non-emitters. Since there is even less diversity among the Mg II emitters than the general galaxy population, this must be kept in mind when interpreting future observational results at high redshift.

Chisholm et al. (2020) has successfully applied their model to a set of observed galaxies and have shown that they can reproduce

<footnote n="5">⁵This can be derived from equation (18) in Chisholm et al. (2020). $F_{obs}/F_{int}$ is the true escape fraction while $e^{-\tau_{chan}}$ is the predicted escape fraction one obtains from the R value in the optically thin limit.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Example spectra of three galaxies where the Chisholm et al. (2020) model accurately predicts the Mg II escape fraction along three different sightlines. In most cases, the galaxies are Mg II emitters however, the emission may include significant contributions from the continuum.</figcaption></figure>

escape fractions estimated from other methods. While in general, the Chisholm et al. (2020) models fail to reproduce the true effective line of sight Mg II escape fractions in our experiments, there are some galaxies that fall along the one-to-one relation in the top panel of Fig. 12 and we aim to explore these further. For each SPHINX²⁰ galaxy, we measure the sum of the squares of the distance between the true effective line of sight escape fraction and the value predicted from the Chisholm et al. (2020) models.

Unsurprisingly, we find that the galaxies where the model works best are the most metal poor (and thus less dusty). These galaxies also tend to be emitters or have complex spectra along the sightline. In general, we find that simulated galaxies that are emitters along the line of sight also tend to have higher Mg II escape fractions along the same sightline. The median line of sight Mg II escape fraction is 91 per cent for emitters while only 63 per cent for non-emitters. In Fig. 13 we show the spectra along three sightlines for three example galaxies where the model works well. In many cases, the line emission from gas dominates the spectrum; however, it is often the case that the continuum is non-negligible (see the bottom row). Thus if the Chisholm et al. (2020) models are to be applied, we recommend selecting particularly low metallicity galaxies where there is confidence that the observed emission line profiles are dominated by emission from gas and there is little evidence of radiative transfer effects.

### 3.3.2 The Henry model

For the ‘picket fence’ model and the optically thick clouds surrounded by optically thin channels model presented in Chisholm et al. (2020), the *R* value alone is not enough to constrain the Mg II optical depth and Mg II escape fraction. Rather, another method is needed to constrain the intrinsic Mg II luminosity. Henry et al. (2018) demonstrated using idealized 1D CLOUDY models that there exists a correlation between [O III]/[O II] and Mg II/[O III]. If one can obtain an intrinsic (dust corrected) measurement of

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Mg II 2796Å/[O III] 5007Å versus O32 for galaxies in SPHINX²⁰ at *z* = 6 (scatter points) compared to the predicted relationship from Henry et al. (2018) (dashed black line). The black dotted lines represent the 0.1 dex scatter around the Henry et al. (2018) relationship. Beige data points represent all haloes in our sample while pink points highlight Mg II emitters with 7.8 ≤ 12 + log₁₀(O/H) ≤ 8.3 which is the metallicity regime where the Henry et al. (2018) relation was derived. The underlaid 2D histogram shows Mg II 2796Å/[O III] 5007Å versus O32 for individual CLOUDY models used in our training data set at *z* = 6. The Henry et al. (2018) relation follows the same trend as our CLOUDY models and entire simulated galaxies; however, the scatter in the simulated galaxies is much larger than that predicted by simple 1D photoionization models.</figcaption></figure>

O32 = log₁₀([O III] 5007Å/[O II] 3727Å) then the intrinsic Mg II luminosity can be calculated. While this method works in idealized models (Henry et al. 2018), here we test whether the calibrated relation applies to entire galaxies in the epoch of reionization.

In Fig. 14, we show Mg II 2796Å/[O III] 5007Å versus O32 for individual SPHINX²⁰ galaxies (data points) as well as the CLOUDY models that were used to predict the luminosity (2D histogram) compared with the predictions from Henry et al. (2018). First, it is clear that the relation presented in Henry et al. (2018) is a relatively good match to many of the CLOUDY models in our simulation. This is indeed a very good consistency check because although the cells in our simulation exhibit a wide range in metallicity, density, and ionization parameter, many should overlap with the parameter choices in Henry et al. (2018). Gas cells in our simulation exhibit significantly more scatter and at high O32, our CLOUDY models predict somewhat higher Mg II 2796Å/[O III] 5007Å compared to Henry et al. (2018). Note that the CLOUDY models used in Henry et al. (2018) were at a fixed density and spanned a smaller range of metallicity and ionization parameter than is represented by gas cells in our simulation so it is not surprising that the scatter is large in comparison.

More importantly, in Fig. 14 we see how SPHINX²⁰ galaxies compare with the Henry et al. (2018) relation. We have differentiated galaxies that are both Mg II emitters and within the metallicity range of Henry et al. (2018) in magenta, versus all others in beige. While the Henry et al. (2018) relation captures the general behaviour of the trend exhibited by SPHINX²⁰ galaxies, consistent with the CLOUDY models, we see that the scatter is significant and that most SPHINX²⁰ galaxies tend to fall above the relation. We speculate that this is because simple CLOUDY models cannot account for
the complexity of the ISM seen in simulations. In some cases, the galaxies fall ~1 dex above the published model. This is potentially very problematic because it would mean that individual predictions for the Mg II escape fraction may be uncertain by a factor of 10.

We emphasize that the O32 and Mg II 2796Å/[O III] 5007Å ratios of individual CLOUDY models are often not very good representations of entire galaxies. This is primarily due to the fact that [O II] and Mg II are emitted from different regions of temperature-density phase space compared to [O III]. For example Mg II can be emitted from neutral gas while this is less likely for [O II]. Our CLOUDY models have been computed at fixed density and temperature and, for example, the cells that are brightest in Mg II and [O II] may be weaker in [O III]. This was also demonstrated to be the case for the [C III] and [O III] infrared emission lines presented in Katz et al. (2022). Thus we conclude that the scatter in the Henry et al. (2018) relation between both Mg II 2796Å/[O III] 5007Å and Mg II 2803Å/[O III] 5007Å and O32 is likely significantly larger than the ~0.1 dex indicated in Henry et al. (2018). Measuring this value directly from the SPHINX²⁰ z = 6 galaxy population, we find a scatter of ~0.2 dex at a typical O32 of 0.5.

Although we argue that the scatter in the Henry et al. (2018) relation is perhaps too large to be useful, it is possible that another relation using the same emission lines may result in better intrinsic Mg II 2796Å luminosity predictions. Rather than using the traditional O32 diagnostic, we fit a generalized linear model to SPHINX²⁰ galaxies at z = 6 using the same three oxygen emission lines.<fnref n="6" /> To calculate the intrinsic Mg II 2796Å luminosity from the three oxygen emission lines, the following equation can be used:

$$ \log_{10}(L_{[Mg II],2796\AA}) = 3.22 + 8.49 \log_{10}(L_{[O II],3726\AA}) \\ - 7.86 \log_{10}(L_{[O II],3728\AA}) \\ + 0.32 \log_{10}(L_{[O III],5007\AA}), \quad (3) $$

where luminosities are in units of erg s⁻¹. The intrinsic Mg II 2803Å luminosity can also be easily calculated from equation (3) by simply halving the Mg II 2796Å luminosity.

In Fig. 15 we compare the true intrinsic Mg II 2796Å luminosity with the value predicted by our generalized linear model. The plot demonstrates that the model provides a non-biased estimator whereas the Henry et al. (2018) prediction (shown as beige points) generally falls below the one-to-one relation, indicating that Mg II 2796Å luminosities are consistently underpredicted. Our generalized linear model is still not perfect as we find a median absolute error of 0.071 dex. However, this is far better than the 0.385 dex found for the Henry et al. (2018) relation. We have attempted to fit a generalized linear model using higher order polynomial combinations of the three oxygen emission lines and found that increasing to degree two or three decreases the median absolute error to 0.054 and 0.044 dex, respectively. We advocate for the lower order polynomial given its simplicity and reduced chance of overfitting.

In the case where the [O II] doublet cannot be separated, the following equation can be used:

$$ \log_{10}(L_{[Mg II],2796\AA}) = -2.39 + 0.85 \log_{10}(L_{[O II],3726,3728\AA}) \\ + 0.21 \log_{10}(L_{[O III],5007\AA}). \quad (4) $$

In this case, the median absolute error increases to 0.12 dex.

<figure><img src="image_16.png" /><figcaption>**Figure 15.** (Top) True, intrinsic Mg II 2796Å luminosity of SPHINX²⁰ galaxies as predicted by our CLOUDY models compared to the predicted value from various models. Our generalized linear models using [O III] and C III] as given in equations (3) and (5) are shown in green and red, respectively. Mg II line luminosities predicted from the Henry et al. (2018) model are shown in beige. The dashed black line shows the one-to-one relation. In general our models provide a non-biased estimate of the intrinsic Mg II 2796Å luminosity of a galaxy with significantly less scatter than that presented in Henry et al. (2018). We provide the median absolute errors (MAE) for each model in the legend. (Bottom) Probability distribution function (PDF) of the log difference between the predicted and intrinsic luminosities for each model. The colours are consistent with the top panel.</figcaption></figure>

We note that the scaling in our relations is considerably different from Henry et al. (2018).<fnref n="7" /> In our equation, Mg II 2796Å luminosity scales with both [O III] and [O II]. In contrast, the Henry et al. (2018) relation has almost no scaling with [O III] and weakly scales as [O III]².

### 3.3.3 C III] as an alternative to [O III]

Due to its longer wavelength, [O III] 5007Å will drop out at a lower redshift than [O II] or Mg II. For this reason, we consider an alternative set of lines that may be more practical at high redshifts. C III] has already been detected at z > 6 (Stark et al. 2017) and because it originates from gas in a high ionization state, we test whether C III] emission can be used as a possible replacement for [O III] when it is not available.

<footnote n="6">⁶Because the intrinsic [O III] 4959Å luminosity is a constant fraction of [O III] 5007Å, the use of this additional line will not improve predictive power.</footnote>

<footnote n="7">⁷They find $\log_{10}(L_{[Mg II],2796\AA}) = 0.079 \times O32^2 - 1.04 \times O32 - 0.54 + \log_{10}(L_{[O III],5007\AA})$.</footnote>
In Fig. 15, we show the intrinsic Mg II 2796Å predicted from C III] versus the true Mg II 2796Å luminosity as pink points. Compared to using [O III], the median absolute error is almost identical, marginally increasing from 0.071 to 0.073 dex. Similar to [O III], we also find that C III] can be used as a non-biased estimator of the intrinsic Mg II 2796Å luminosity.

$$ \log_{10}(L_{[\text{Mg II}, 2796\AA]}) = 2.46 + 9.09 \log_{10}(L_{[\text{O III}, 3726\AA]}) \\ - 8.44 \log_{10}(L_{[\text{O III}, 3728\AA]}) \\ + 0.32 \log_{10}(L_{\text{C III], 1906, 1908\AA}}). \quad (5) $$

Note once again that the intrinsic Mg II 2803Å luminosity can also be calculated by halving the intrinsic Mg II 2796Å luminosity.

### 3.3.4 Predicting the Mg II escape fraction with UV and optical emission lines

The final question we aim to address in this section is how well can our newly calibrated relations constrain the Mg II escape fraction in the epoch of reionization? To do this, we assume that the intrinsic UV and optical emission line luminosities can be dust-corrected and are thus known. We then use these lines to measure the intrinsic Mg II 2796Å and Mg II 2803Å luminosities and compute the effective escape fraction from the observed Mg II 2796Å and Mg II 2803Å luminosities.

In the left-hand panel of Fig. 16 we show the predicted values for the effective Mg II escape fraction using the [O III] and C III] relations derived in this work as well as the relation from Henry et al. (2018). While our predictions tend to scatter around the one-to-one line, we find two clear biases. At escape fractions of ~50 per cent we tend to overpredict the true value while at $\tilde{f}_{\text{esc}} \gtrsim 80$ per cent we under predict the true values. Nevertheless, the median absolute error for both of our models is only 0.11. In other words, we expect the typical galaxy will have an absolute error of ~11 per cent in the escape fraction prediction. The Henry et al. (2018) relation significantly overpredicts the escape fractions for nearly all simulated galaxies due to the underprediction of the intrinsic Mg II luminosities as noted in Fig. 15.

The left-hand panel of Fig. 16 shows all galaxies, rather than only the Mg II emitters. However, when isolating only the emitters, we find the exact same biases and trend in the escape fraction predictions as the non-emitters.

Observations, of course, can only measure the Mg II flux along a single line of sight. In the centre panel of Fig. 16, we apply our relations to predict the effective line-of-sight Mg II escape fraction for three viewing angles for each SPHINX²⁰ galaxy at $z = 6$. Surprisingly, the median absolute error of our predictions remains the same compared to the angle-averaged predictions, despite the large scatter between the line of sight and angle-averaged Mg II escape fractions for the same galaxy. Furthermore, it seems that some of the bias in the predictions is reduced at intermediate escape fractions; although, our relations still tend to underpredict higher values of $\tilde{f}_{\text{esc,los}}$. As expected, the Henry et al. (2018) relation once again drastically overpredicts the line-of-sight Mg II escape fractions for all galaxies.

There exists a small sample of galaxies at low redshift that we can apply our relations too and try to directly predict the Mg II escape fraction. In the right-hand panel of Fig. 16 we show the observed R values versus the predicted Mg II escape fractions for a subsample of galaxies from Feltre et al. (2018) and Henry et al. (2018). Galaxies were selected based on whether there were observations of the emission lines required to apply equation (3) or (4) and dust

corrections were applied to the oxygen emission line luminosities when necessary. In general, we find that the observed galaxies tend to exhibit Mg II escape fractions below ~30 per cent. While this is not inconsistent with our simulated $z = 6$ galaxies, where escape fractions along individual sightlines approach 0 per cent, if the observed galaxies were representative of the high-redshift population, we might have expected to observe a few with higher Mg II escape fractions. We have shown that there exist important differences between the Green Pea galaxy population and our simulated $z = 6$ galaxies so this must be kept in mind when interpreting the results. Furthermore, depletion of Mg on to dust may be more important for lower redshift galaxies.

In summary, we predict that our newly calibrated models are able to constrain both the global/angle-averaged and line-of-sight Mg II escape fractions with an absolute accuracy of ~11 per cent for the typical galaxy in the epoch of reionization.

### 3.4 LyC $f_{\text{esc}}$ from Mg II

Having demonstrated a model that adequately predicts the Mg II escape fraction for high-redshift galaxies, in this section we test whether Mg II can be used to constrain the LyC escape fraction.

#### 3.4.1 The Chisholm model

Chisholm et al. (2020) show that once the Mg II optical depth is determined, this can be converted into a neutral hydrogen column density as long as the metallicity is known. While Chisholm et al. (2020) compute the optical depth from the Mg II doublet flux ratio, which we have shown is likely an inadequate representation of the physics in our simulated galaxies, an effective optical depth can also be computed if the Mg II escape fraction is known, for example, by using the calibrated relations in the previous section. Here, we test the accuracy of this approach.

First, in the top left-hand panel of Fig. 17, we compare the angle-averaged LyC escape fractions of SPHINX²⁰ galaxies at $z = 6$ with both their Ly$\alpha$ escape fractions and their Mg II escape fractions. We find a stronger trend between the LyC escape fraction and that of Ly$\alpha$ compared to Mg II. While LyC leakers in SPHINX²⁰ also tend to have higher Mg II escape fractions, there is a significant amount of scatter which questions whether the Mg II escape fraction has predictive power for the global LyC escape fraction. Interestingly, the situation is slightly different along individual sightlines. In this case, we find correlations amongst all three quantities as shown in the bottom row of Fig. 17, albeit with considerable scatter.

As we have demonstrated earlier, the angle-averaged Mg II escape fraction is not a probe of the neutral hydrogen column density, but rather a probe of dust. However, this is not necessarily the case for the Mg II escape fraction along individual sightlines. Thus, we test the Chisholm et al. (2020) approach using our simulated line of sight escape fractions. As discussed earlier, the models have the best chance of working when the Mg II spectra appears optically thin along the line of sight. For all 2,082 sightlines, we have manually classified the Mg II spectra as appearing optically thin or not. This subjective determination is made by selecting only galaxies with strong emission peaks in their total spectra where the emission line profiles look relatively Gaussian. 19 per cent of all sightlines are given this classification. We note that this determination is made on the total spectra (gas emission plus stellar continuum), but we continue to assume that the gas emission can be perfectly separated from the continuum, even if the continuum contributes to
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Predicted versus true angle-averaged (left-hand panel) or line of sight (centre) effective Mg II escape fraction for SPHINX²⁰ galaxies at z = 6. The green and pink points show the predictions when using our [O III] and C III] relations as given in equations (3) and (5), respectively. The beige points show the predictions when applying the relation from Henry et al. (2018). (Right-hand panel) Predicted Mg II escape fractions versus observed R for a subsample of low-redshift galaxies from Feltre et al. (2018), Henry et al. (2018), Izotov et al. (2016a,b, 2018a,b, 2020), Guseva et al. (2020). Galaxies were selected based on whether emission line observations existed such that equation (3) or (4) could be applied.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** (Top left-hand panel) Angle-averaged LyC escape fraction versus the angle-averaged Mg II (green) or Lyα (beige) escape fraction for SPHINX²⁰ galaxies at z = 6. We find a much stronger trend between the LyC escape fraction and the Lyα escape fraction compared to LyC and Mg II. (Top Centre) Predicted line of sight LyC escape fraction from the Chisholm et al. (2020) relation either assuming a pure hydrogen gas (green) or a gas and dust mixture (beige) for SPHINX²⁰ galaxies at z = 6. Darker points represent galaxies that have been identified as having spectra that appear optically thin along the line of sight while lighter points represent galaxies impacted by significant radiative transfer effects. The dashed black line shows the one-to-one relation. We use the method presented in Chisholm et al. (2020) to compute the neutral hydrogen column density from the Mg II optical depth and the dust content from the Balmer decrement. Models without dust significantly overpredict the LyC escape fraction. (Top right-hand panel) Predicted line-of-sight LyC escape fraction from the Chisholm et al. (2020) relation assuming a gas and dust mixture using a SMC bar extinction law (Gordon et al. 2003) instead of a Reddy et al. (2016) extinction law as well as the true intrinsic Hα/Hβ ratio for calculating E(B–V) rather than 2.86. (Bottom left-hand panel) Line-of-sight LyC escape fraction versus the line-of-sight Lyα escape fraction. (Bottom centre) Line-of-sight LyC escape fraction versus the line-of-sight Mg II escape fraction. (Bottom right-hand panel) Line-of-sight Mg II escape fraction versus the line-of-sight Lyα escape fraction. Note that escape fractions can be greater than one, especially for resonant lines, due to radiative transfer effects.</figcaption></figure>
the emission line via back scattering. This may be more difficult observationally and will change the perceived R value.

Following the methodology in Chisholm et al. (2020), we compute the Mg II column density as

$$N_{\text{Mg II}} = 3.8 \times 10^{14} \frac{\AA}{\text{cm}^2} \frac{\tau_{2803}}{f\lambda}, \qquad (6)$$

where $\lambda = 2803.53\AA$, $f = 0.303$ is the oscillator strength and $\tau_{2803} = -\ln(f_{\text{esc, 2803}})$, with $f_{\text{esc, 2803}}$ computed using the UV/optical emission line method presented above assuming the intrinsic luminosities are known. The neutral hydrogen column density is then estimated as $N_{\text{H I}} = 46 H/O N_{\text{Mg II}}$ assuming solar abundances from Asplund et al. (2009). Note that this scaling implicitly assumes that all Mg II gas is found in neutral regions. Since Mg II can also exist in ionized regions, there may be additional unaccounted for scatter in this relation. In a dust-free scenario, the LyC escape fraction can then be computed as $f_{\text{esc, LyC}} = e^{-N_{\text{H I}}\sigma}$, where $\sigma = 6.3 \times 10^{-18} \text{cm}^{-2}$ is the hydrogen photoionization cross-section at 912Å. However, in the case where dust is present, the escape fraction is further modulated by a factor of $10^{-0.4E(B-V)k(912)}$, where k(912) is the dust extinction at 912Å.

In the top centre panel of Fig. 17 we show the predicted line of sight LyC escape fraction estimated with (beige) and without (green) dust versus the true LyC escape fraction for three different sightlines for each SPHINX²⁰ galaxy. We find that the pure hydrogen model significantly over-predicts the true LyC escape fraction while the hydrogen and dust model follows the one-to-one relation, albeit with a significant amount of scatter. For the model with dust, we find a median absolute error of the optically thin galaxies is 0.8 dex. This value is slightly smaller (0.67 dex) for galaxies that have a predicted LyC escape fraction > 1 per cent.

Because the true extinction law at high-redshift is unknown as is the intrinsic $H\alpha/H\beta$ ratio, we have followed Chisholm et al. (2020) and adopted a Reddy et al. (2016) attenuation curve and a value of 2.86 (Osterbrock 1989) for the intrinsic $H\alpha/H\beta$ ratio. These assumptions will introduce systematic uncertainties into the calculation as in our simulations, dust is modelled using the synthetic extinction curve for the SMC bar (Gordon & Clayton 1998; Weingartner & Draine 2001; Gordon et al. 2003) and the intrinsic $H\alpha/H\beta$ ratios deviate from 2.86 (due to e.g. contributions from collisional emission and temperature deviations from $10^4$ K). Nevertheless, our goal is to follow as closely as possible to the methods used in the literature. Note that these two assumptions only impact the predictions for the hydrogen and dust mixture.

To better understand how these assumptions impact our predictions, in the top right-hand panel of Fig. 17 we show the predicted LyC escape fractions versus the true values when using the SMC bar extinction law as well as the true intrinsic $H\alpha/H\beta$ ratios. In this case, the scatter is as significant as the previous method and we see a systematic shift towards higher predicted escape fractions. This is because $E(B-V) \times k(912)$ tends to be higher for the Reddy et al. (2016) extinction law assuming an intrinsic $H\alpha/H\beta$ of 2.86 than it is for the SMC bar using the true intrinsic $H\alpha/H\beta$ ratio.

Computing a precise value for the LyC escape fraction using indirect proxies is certainly a difficult problem. The concerning aspect in this section is that the escape fractions are drastically overpredicted when dust is not included in the calculation. Low-mass high-redshift galaxies are expected to have little dust, so in principle, the escape fraction should be predominantly dictated by the neutral hydrogen content (see e.g. fig. 5 of Kimm et al. 2019). The median line-of-sight neutral hydrogen column density computed from the Mg II method is only $2.3 \times 10^{17} \text{cm}^{-2}$ which would correspond to

an escape fraction of 23 per cent. This is much higher than even the luminosity-weighted escape fraction of SPHINX²⁰ galaxies, which is ~ 3 per cent at $z = 6$ (Rosdahl et al., in preparation).

To demonstrate that high-redshift LyC escape fractions are dominated by neutral hydrogen rather than dust, we have repeated our escape fraction calculation but have removed dust. We find a maximum fractional difference (i.e. $|f_{\text{esc, LyC, w/dust}} - f_{\text{esc, LyC, w/dust}}|$) of 3 per cent along any given sightline or angle average. This in itself does not prove that the escape fraction is dominated by neutral hydrogen. If photons primarily escape via very optically thin channels and the optically thick regions are thick to both neutral hydrogen and dust, then the escape fraction should not change with the removal of either neutral hydrogen or dust. However, repeating the calculation a third time without neutral hydrogen results in substantially increased LyC escape fractions, thereby confirming our hypothesis.

Unfortunately, our calculations show that the line-of-sight Mg II escape fractions are like the angle-averaged values and probes of dust rather than neutral hydrogen. Thus the Mg II escape fraction is not an ideal probe of LyC escape at high redshift. Perhaps the shape of the spectra, the peak separation, or deviation from line-centre can provide a more robust diagnostic of the Mg II and neutral hydrogen optical depth as it does for $H\alpha$ (e.g. Verhamme et al. 2015, 2017); however, this is beyond the scope of this work.

## 4 CAVEATS

There are numerous caveats that arise when modelling emission lines from cosmological radiation hydrodynamics simulations. In the context of the simulations used in this work, many are discussed in Katz et al. (2022). Here, we focus on a few issues that are particularly important for resonant lines and Mg.

The emergent luminosities and spectral profiles of resonant lines are very sensitive to the small scale distribution of gas and dust as well as the velocity and turbulent structure of the medium. Like all simulations, SPHINX²⁰ has a finite spatial resolution of ~10 pc at $z = 6$ which is not enough to fully resolve the detailed structure of the ISM or giant molecular clouds (e.g. Kimm et al. 2022). Both the spectral profiles and emergent luminosities are expected to change for a simulation of both lower and higher resolution (e.g. Camps et al. 2021). Correctly modelling the ISM is key for accurately modelling emission lines and thus it is important to consider that any issues with the ISM model in SPHINX will lead us to systematic biases in our results.

We have adopted solar abundance ratios that are scaled by metallicity throughout this work whereas it is well-known that different chemical enrichment pathways can result in non-solar abundance patterns that in many cases vary with metallicity (e.g. Maiolino & Mannucci 2019). Being an $\alpha$ element, there is a trend of higher [Mg/Fe] a low metallicity; however, Mg can also be produced in other phenomena besides core-collapse supernova that cause a decrease in [O/Mg] with increasing metallicity (e.g. Nomoto, Kobayashi & Tominaga 2013). Such effects can only be captured in simulations that specifically model the yields from different types of SNe, stellar winds, and other enrichment phenomena. Because high-redshift galaxies are more likely to be enriched by core-collapse SN, not modelling the varying abundance patterns may under predict the Mg II luminosities.

We have not accounted for the depletion of Mg on to dust. In certain models, Mg can be heavily depleted by up to 0.7–1.0 dex (e.g. Dopita et al. 2000; Jenkins 2009). Although this effect is expected to be weaker at low metallicity (Guseva et al. 2019) and it has
been shown that Mg is less depleted in the SMC than it is for the Milky Way (Jenkins & Wallerstein 2017), it is unlikely to be completely negligible. Neglecting Mg depletion would likely result in an overprediction of Mg II luminosities, opposite to impact of varying abundances so it is likely that these two effects partially cancel each other out. Taking values from Jenkins & Wallerstein (2017), intrinsic Mg II luminosities may be overpredicted by a factor of ~2.5; however, this is unlikely to impact the intrinsic or emergent flux ratios of the Mg II doublet. A reduction in gas phase Mg II may lead to increased line-of-sight escape fractions in galaxies where the escape fraction is not set by dust.

The dust model used in this work (Laursen et al. 2009) is not a prediction by the simulation but rather a phenomenological model where the dust scales with the neutral hydrogen content of each gas cell. Constraints on both the dust properties and masses of high-redshift galaxies are limited and our results will be sensitive to the choice of model. For example, if we assumed that the gas was completely dust-free, the emergent Mg II line ratios would only be sensitive to the geometry of the gas along the observed sightline which would have the effect of generally increasing the Mg II escape fractions.

As discussed in Mauerhofer et al. (2021), the Doppler parameter which enters the computation of optical depth is sensitive to both the thermal and turbulent velocity of the gas. The simulation self-consistently models the temperature of the gas and thus the thermal velocity is, in principle, calculated accurately. However, there is no unique model for subgrid turbulence. Turbulence is less important for Ly α than it is for Mg II because the thermal velocity is inversely proportional to the atomic mass (which is 24 × higher for Mg than it is for H). Since our fiducial model does not include subgrid turbulence and only considers the thermal and bulk velocity of the gas cells (which already incorporates the turbulence on the scales that we resolve), it is possible that we may be overestimating the optical depth at line centre and underestimating it further away from the line. If we instead assume the microturbulence model following observations from Larson (1981) where $v_{\text{turb}} = 1.1 \Delta x^{0.38} \text{ km s}^{-1}$ and $\Delta x$ is the physical size of the gas cell in pc, we find that the turbulent velocity is approximately equal to the thermal velocity for the most refined gas cells in our simulation at a temperature of $10^4$ K. At temperatures below this value, the turbulent velocity will dominate. We demonstrate the impact of our fiducial assumption compared to this alternative microturbulence model in Appendix B.

We have not included the impact of the nebular continuum on the spectra. While this does not impact our work on the Mg II escape fraction, the flux ratios, or their use as a probe of LyC escape, it may impact the fraction of galaxies that are seen as emitters versus absorbers.

## 5 CONCLUSIONS

In this work, we have used SPHINX²⁰, a state-of-the-art cosmological radiation hydrodynamics simulation of galaxy formation in the epoch of reionization to study the Mg II emission line properties of high-redshift galaxies at $z=6$. The resonant line Mg II doublet has recently been proposed as an indicator of the LyC escape fraction (Henry et al. 2018; Chisholm et al. 2020) and will be targeted in Cycle 1 JWST observations (e.g. Chisholm et al. 2021). Our primary goals were to determine the fraction of high-redshift galaxies that are Mg II emitters, to determine whether Mg II can be used as a probe of LyC escape, and to test the Henry et al. (2018) and Chisholm et al.

(2020) models for estimating the Mg II and LyC escape fractions from Mg II observations. Our main conclusions can be summarized as follows:

(i) When compared with Green Peas and Blueberry galaxies (which are thought to be low-redshift analogues of reionization-epoch galaxies), SPHINX²⁰ galaxies at $z=6$ are less extreme in terms of their emission line luminosities (e.g. Balmer and nebular oxygen lines) despite exhibiting similar stellar masses, star formation rates, and metallicities. This disagreement is exacerbated when comparing to only low-redshift LyC leakers that exhibit significantly higher sSFRs compared to SPHINX²⁰. We highlight that low-redshift galaxies are often selected based on their extreme emission line properties and there is no guarantee that this selection provides direct analogues of the bulk of the $z=6$ galaxy population that we simulate with SPHINX²⁰. Likewise, as we discuss in the Section 4, due to various numerical effects and assumptions, SPHINX²⁰ galaxies may not be perfect replicas of $z=6$ galaxies. These differences should be kept in mind when developing indirect probes of LyC leakage from both simulations and observations.

(ii) High-redshift galaxies are predicted to exhibit a diversity of Mg II spectral shapes including emission, absorption, and P-Cygni profiles. However, the majority of $z>6$ galaxies are predicted to be Mg II emitters or have complex spectra with both emission and absorption features. There is a tendency for galaxies with higher stellar masses, SFRs, and LyC escape fractions to be Mg II emitters.

(iii) We find that the model proposed by Chisholm et al. (2020) that uses the Mg II doublet ratio, R, as an indicator of the Mg II escape fraction does not sufficiently capture the physics that dictates Mg II escape in our simulations except for some of the lowest metallicity systems where emission from gas dominates the spectra. This is due to both the importance of dust in our simulation as well as the geometry of SPHINX²⁰ galaxies not being well represented as a Mg II source behind a screen of Mg II gas.

(iv) The model proposed by Henry et al. (2018) to compute the intrinsic Mg II emission from [O II] and [O III] emission lines tends to underpredict the intrinsic Mg II luminosities in our simulation. We attribute this to the fact that the Henry et al. (2018) model was calibrated on CLOUDY models that do not fully capture the diversity of ISM properties in our simulated galaxies. We provide newly calibrated relations in equations (3), (4), and (5).

(v) The model proposed by Chisholm et al. (2020) to convert a Mg II escape fraction to a LyC escape fraction has a typical error of 1.14 dex. Moreover, when applied to our galaxies, the model predicts that dust plays a stronger role in setting the LyC escape fraction compared to neutral hydrogen. This is opposite what we find in our simulations as removing dust results in a maximum fractional difference of 3 per cent in the LyC escape fraction.

We have shown that the complex physics that governs resonant line emission from high-redshift galaxies can make its interpretation difficult, especially since the same galaxy can exhibit numerous, different spectral profiles along different lines of sight. However, if a large sample of galaxies can be obtained, using the relations derived in this work, there is a prospect of utilizing Mg II emission as an aggregate probe LyC escape at $z>6$ for optically thin systems. This method will be complementary to those in other wavebands, such as the IR, and thus provide an orthogonal measure of LyC leakage in the epoch of reionization. Future work aimed at studying the spectral profiles of Mg II emission may reveal correlations with LyC escape as has been shown for Ly α.
## ACKNOWLEDGEMENTS

TG is supported by the ERC Starting grant 757258 ‘TRIPLE’. 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. TK was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2020R1C1C1007079). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 693024).

## DATA AVAILABILITY

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

## REFERENCES

Amorín R. et al., 2017, Nat. Astron., 1, 0052

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

Atek H., Richard J., Kneib J.-P., Schaerer D., 2018, MNRAS, 479, 5184

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

Borthakur S., Heckman T. M., Leitherer C., Overzier R. A., 2014, Science, 346, 216

Bouwens R. J., Oesch P. A., Illingworth G. D., Ellis R. S., Stefanon M., 2017, ApJ, 843, 129

Burchett J. N., Rubin K. H. R., Prochaska J. X., Coil A. L., Vaught R. R., Hennawi J. F., 2021, ApJ, 909, 151

Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchibergmann T., 2000, ApJ, 533, 682

Camps P., Behrens C., Baes M., Kapoor A. U., Grand R., 2021, ApJ, 916, 39

Cen R., Kimm T., 2015, ApJ, 801, L25

Chisholm J., Prochaska J. X., Schaerer D., Gazagnes S., Henry A., 2020, MNRAS, 498, 2554

Chisholm J. et al., 2021, The First Observations of the Ionizing Luminosity of Galaxies within the Epoch of Reionization, JWST Proposal. Cycle 1

Dijkstra M., 2017, preprint (arXiv:1704.03416)

Dopita M. A., Kewley L. J., Heisler C. A., Sutherland R. S., 2000, ApJ, 542, 224

Drake A. B. et al., 2017, A&A, 608, A6

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

Feltre A. et al., 2018, A&A, 617, A62

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

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

Finkelstein S. L. et al., 2015, ApJ, 810, 71

Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

Garel T., Blaizot J., Rosdahl J., Michel-Dansac L., Haehnelt M. G., Katz H., Kimm T., Verhamme A., 2021, MNRAS, 504, 1902

Gordon K. D., Clayton G. C., 1998, ApJ, 500, 816

Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279

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

Guseva N. G., Izotov Y. I., Fricke K. J., Henkel C., 2019, A&A, 624, A21

Guseva N. G. et al., 2020, MNRAS, 497, 4293

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

Henry A., Berg D. A., Scarlata C., Verhamme A., Erb D., 2018, ApJ, 855, 96

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

Izotov Y. I., Schaerer D., Thuan T. X., Worseck G., Guseva N. G., Orlitová I., Verhamme A., 2016a, MNRAS, 461, 3683

Izotov Y. I., Orlitová I., Schaerer D., Thuan T. X., Verhamme A., Guseva N. G., Worseck G., 2016b, Nature, 529, 178

Izotov Y. I., Schaerer D., Worseck G., Guseva N. G., Thuan T. X., Verhamme A., Orlitová I., Fricke K. J., 2018a, MNRAS, 474, 4514

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018b, MNRAS, 478, 4851

Izotov Y. I., Schaerer D., Worseck G., Verhamme A., Guseva N. G., Thuan T. X., Orlitová I., Fricke K. J., 2020, MNRAS, 491, 468

Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Chisholm J., Thuan T. X., Fricke K. J., Verhamme A., 2021, MNRAS, 503, 1734

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

Jenkins E. B., Wallerstein G., 2017, ApJ, 838, 85

Kakiichi K. et al., 2018, MNRAS, 479, 43

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

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. et al., 2022, MNRAS, 510, 5603

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

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

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

Konno A. et al., 2018, PASJ, 70, S16

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

Laursen P., Sommer-Larsen J., Andersen A. C., 2009, ApJ, 704, 1640

Leitherer C., Hernandez S., Lee J. C., Oey M. S., 2016, ApJ, 823, 64

Livermore R. C., Finkelstein S. L., Lotz J. M., 2017, ApJ, 835, 113

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

Martin C. L., Shapley A. E., Coil A. L., Kornei K. A., Murray N., Pancoast A., 2013, ApJ, 770, 41

Mauerhofer V., Verhamme A., Blaizot J., Garel T., Kimm T., Michel-Dansac L., Rosdahl J., 2021, A&A, 646, A80

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

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

Nelson D., Byrohl C., Peroux C., Rubin K. H. R., Burchett J. N., 2021, MNRAS, 507, 4445

Nomoto K., Kobayashi C., Tominaga N., 2013, ARA&A, 51, 457

Osterbrock D. E., 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, New York

Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544

Pearson W. J. et al., 2018, A&A, 615, A146

Reddy N. A., Steidel C. C., Pettini M., Bogosavljević M., 2016, ApJ, 828, 107

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

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

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

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

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

Rubin K. H. R., Prochaska J. X., Ménard B., Murray N., Kasen D., Koo D. C., Phillips A. C., 2011, ApJ, 728, 55

Schaerer D. et al., 2020, A&A, 643, A3
Siana B. et al., 2015, ApJ, 804, 17

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

Stark D. P. et al., 2017, MNRAS, 464, 469

Steidel C. C., Bogosavljević M., Shapley A. E., Reddy N. A., Rudie G. C., Pettini M., Trainor R. F., Strom A. L., 2018, ApJ, 869, 123

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

Tran K.-V. H. et al., 2020, ApJ, 898, 45

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

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

Verhamme A., Orlitová I., Schaerer D., Hayes M., 2015, A&A, 578, A7

Verhamme A., Orlitová I., Schaerer D., Izotov Y., Worseck G., Thuan T. X., Guseva N., 2017, A&A, 597, A13

Wang B., Heckman T. M., Leitherer C., Alexandroff R., Borthakur S., Overzier R. A., 2019, ApJ, 885, 57

Wang B. et al., 2021, ApJ, 916, 3

Weingartner J. C., Draine B. T., 2001, ApJ, 548, 296

Whitaker K. E. et al., 2014, ApJ, 795, 104

Xu H., Wise J. H., Norman M. L., Ahn K., O'Shea B. W., 2016, ApJ, 833, 84

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

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

Yusef-Zadeh F., Morris M., White R. L., 1984, ApJ, 278, 186

Zabl J. et al., 2021, MNRAS, 507, 4294

Zackrisson E., Inoue A. K., Jensen H., 2013, ApJ, 777, 39

Zackrisson E. et al., 2017, ApJ, 836, 78

Zheng Z., Miralda-Escudé J., 2002, ApJ, 578, 33

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Figure S1.** Emissivity of $Ly\alpha$, $H\alpha$, $H\beta$, and He II as a function of electron temperature.

**Figure S2.** Comparison of the collisional emissivities used in this work (dark green) versus others in the literature for $Ly\alpha$ (left), $H\alpha$ (centre), and $H\beta$ (right).

**Figure S3.** Comparison of Mg II spectra along three different lines of sight for simulations with (solid lines) and without (dotted lines) subgrid turbulence.

**Table S1.** Fitting parameters for collisional emissivities given by equation A3.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a TeX/LA TeX file prepared by the author.