# The Sphinx public data release: Forward modelling high-redshift JWST
observations with cosmological radiation hydrodynamics simulations

Harley Katz1∗ Joki Rosdahl2 Taysun Kimm3 Jeremy Blaizot2 Nicholas Choustikov1
Marion Farcy4 Thibault Garel5 Martin G. Haehnelt6 Leo Michel-Dansac2 and
Pierre Ocvirk7 1Sub-department of Astrophysics, University of Oxford, Keble
Road, Oxford OX1 3RH, United Kingdom  2Centre de Recherche Astrophysique de
Lyon UMR5574, Univ Lyon1, ENS de Lyon, CNRS, F-69230 Saint-Genis-Laval, France
3Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul
03722, Republic of Korea  4Institute for Physics, Laboratory for Galaxy
Evolution and Spectral Modelling, EPFL, Observatoire de Sauverny, Chemin
Pegasi 51, 1290 Versoix, Switzerland  5Observatoire de Genève, Université de
Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland  6Kavli Institute for
Cosmology and Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK
7Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS UMR
7550, 11 rue de l’Université, 67000 Strasbourg, France

###### Abstract

The recent launch of JWST has ushered in a new era of high-redshift astronomy
by providing detailed insights into the gas and stellar populations of
galaxies in the epoch of reionization. Interpreting these observations and
translating them into constraints on the physics of early galaxy formation is
a complex challenge that requires sophisticated models of star formation and
the interstellar medium (ISM) in high-redshift galaxies. To this end, we
present Version 1 of the Sphinx20 public data release. Sphinx20 is a full box
cosmological radiation hydrodynamics simulation that simultaneously models the
large-scale process of cosmic reionization and the detailed physics of a
multiphase ISM, providing a statistical sample of galaxies akin to those
currently being observed by JWST. The data set contains $\sim 14,000$ mock
images and spectra of the stellar continuum, nebular continuum, and 52 nebular
emission lines, including Ly$\alpha$, for each galaxy in Sphinx20 with a star
formation rate $\geq 0.3\ {\rm M_{\odot}\ yr^{-1}}$. All galaxy emission has
been processed with dust radiative transfer and/or resonant line radiative
transfer, and data is provided for ten viewing angles for each galaxy.
Additionally, we provide a comprehensive set of intrinsic galaxy properties,
including halo masses, stellar masses, star formation histories, and ISM
characteristics (e.g., metallicity, ISM gas densities, LyC escape fractions).
This paper outlines the data generation methods, presents a comparative
analysis with JWST ERS and Cycle 1 observations, and addresses data set
limitations. The Sphinx20 data release can be downloaded at the following URL:
https://github.com/HarleyKatz/SPHINX-20-data.

###### keywords:

high-redshift galaxies, ISM, galaxy spectra, reionization, galaxy formation

††thanks: ∗E-mail<EMAIL_ADDRESS>
,

##  1 Introduction

Elucidating the underlying physics that governs galaxy formation at high-
redshift is one of the primary goals of extragalactic spectroscopic and
imaging surveys with JWST (Gardner et al., 2006). However, understanding how
to map the features of a galaxy image or spectra to the underlying properties
of the stellar populations and gas as well as the physical processes driving
the evolution of the interstellar and circumgalactic medium (ISM, CGM)
represents a key theoretical challenge.

Prior to JWST, high-redshift space-based observations were limited to
primarily photometric surveys with the Hubble Space Telescope (HST) and
Spitzer. The occasional follow-up observation with the WFC grism on HST or
ground-based telescopes (e.g. ALMA, MOSFIRE on Keck, or MUSE on the VLT)
provided additional spectroscopic redshift confirmation for some of the high-
redshift candidates (e.g. Stark et al., 2017; Hashimoto et al., 2018; Jiang et
al., 2021; Inami et al., 2017). These telescopes have been used predominantly
to constrain high-redshift population statistics such as the UV luminosity
function (Bouwens et al., 2015; Livermore, Finkelstein & Lotz, 2017), the
global growth of star formation or UV luminosity density as a function of
redshift (Ellis et al., 2013; Oesch et al., 2018), and galaxy morphology (e.g.
Kawamata et al., 2018; Bouwens et al., 2022). The photometry of individual
galaxies has been used to measure early star formation histories (SFHs) with
spectral energy distribution (SED) fitting codes (e.g. Chevallard & Charlot,
2016; Carnall et al., 2018; Johnson et al., 2021) or to infer the presence of
strong emission lines (e.g. by IR excesses, Roberts-Borsani et al. 2016),
which can be used to measure quantities such as the ionizing photon production
efficiency ($\xi_{\rm ion}$, e.g. Bouwens et al., 2016; Stefanon et al.,
2022). JWST significantly improves upon HST+Spitzer data for numerous reasons.
Not only is JWST more sensitive and has higher spatial resolution compared to
its predecessors, but it also has significantly more filters that probe a
wavelength range overlapping and between HST and Spitzer that is ideal for
capturing rest-frame UV and optical photons at high-redshift. Hence the
accuracy by which one can constrain quantities such as the SFH or stellar mass
of a galaxy, purely from photometry is greatly improved with JWST. Likewise,
the additional spatial resolution and sensitivity allows for spatially
resolved constraints on these properties, deeper in rest-frame UV/optical
luminosity that were impossible before.

Possibly the most important advancement with JWST is its spectroscopic
capabilities. Not only can photometric redshifts be readily confirmed (or not)
with NIRSpec or MIRI (e.g. Arrabal Haro et al., 2023; Roberts-Borsani et al.,
2023), but emission and absorption lines can also be detected which provide an
unprecedented view into the ISM and CGM of early galaxies (e.g. Cameron et
al., 2023; Sanders et al., 2023). The first three NIRSpec observations of
$z>7$ galaxies with JWST released as part of the ERO data have already shown
that the important heating mechanisms in the ISM at high-redshift are
potentially significantly different compared to what is observed in the local
Universe (Katz et al., 2023a). There are now numerous examples of high-
redshift JWST galaxies where the spectrum appears nothing like the galaxies
observed at low-redshift (e.g. Bunker et al., 2023b; Isobe et al., 2023).

Despite the transformative nature of JWST, interpreting the observational data
remains challenging. For example, SED fitting codes are the primary tools by
which photometric redshifts, stellar masses, and SFHs are constrained from
photometry (e.g. Chevallard & Charlot, 2016; Carnall et al., 2018; Johnson et
al., 2021). Outputs from different SED fitting codes can vary significantly in
their estimates of stellar population properties (e.g. Pacifici et al., 2023).
This is particularly problematic at high-redshift where on-going star
formation can dominate the SED leaving the SFH and stellar mass estimates very
sensitive to the chosen prior (e.g. Narayanan et al., 2023). Furthermore,
there exists significant disagreement between spatially resolved and galaxy-
integrated stellar mass estimates (Giménez-Arteaga et al., 2023).

Historically, 1D photoionization or shock models, e.g. those run with CLOUDY
(Ferland et al., 2017) or MAPPINGS (Sutherland & Dopita, 2017), overwhelmingly
dominate the interpretation of galaxy spectra, in particular, for measuring
quantities such as the ISM electron density, temperature, metallicity, and
ionization parameter. These models have been extremely successful in
explaining the physics behind commonly used strong-line diagnostics at low and
high redshift (e.g. Kewley et al., 2001; Strom et al., 2017; Nakajima et al.,
2022). Because of their 1D nature, such models are computationally inexpensive
such that large grids of models with varying physical complexities can be
employed to map the expected parameter space of galaxy properties from the
local Universe to the epoch of reionization (e.g. Morisset, Delgado-Inglada &
Flores-Fajardo, 2015).

However, the simplicity of photoionization and shock models also represents a
limitation. The real ISM and CGM of galaxies are multiphase, exhibit turbulent
geometries — that can deviate significantly from the typically assumed gas
slab or spherical cloud — and include highly dynamic processes such as stellar
winds, radiation pressure, and supernova (SN) feedback, that are much more
difficult, if not impossible, to capture in 1D models. For this reason, there
has been substantial recent effort to expand the capabilities of
photionization models to three dimensions (e.g. Ercolano et al., 2008; Gray &
Scannapieco, 2017; Jin, Kewley & Sutherland, 2022). Moreover, in certain
situations, the underlying assumptions in photoionization models can highly
bias inferences of galaxy properties. This has been shown to be the case for
direct method metallicity measurements, which are widely regarded as the gold
standard, but systematically underestimate metallicity in the presence of
temperature fluctuations (e.g. Cameron, Katz & Rey, 2022; Méndez-Delgado et
al., 2023b). A similar bias also impacts electron density estimates (Méndez-
Delgado et al., 2023a).

One possible solution to the challenges associated with applying SED fitting
codes and photoionization models to high-redshift data is to calibrate them
with low-redshift “analog” galaxies, for which there is ample data. This
technique has been extensively applied to extreme emission line galaxies (e.g.
Izotov et al., 2016; Amorín et al., 2017; Flury et al., 2022) to deduce
properties like the Lyman-$\alpha$ (Ly$\alpha$) or Lyman continuum (LyC)
escape fractions. However, there is some debate about the similarity of the
low-redshift analogue population to those galaxies that form at $z\gtrsim 6$
(Katz et al., 2023a; Schaerer et al., 2022). Furthermore, the issue remains
that the true underlying galaxy properties are still unknown in analogs.

For this reason, cosmological simulations are a valuable tool for
understanding galaxy formation at early times. Such simulations are important
in at least two ways:

  1. 1.

They represent a complementary means to low-redshift analogs in testing the
ability of more simplistic models to recover the true intrinsic parameters of
a galaxy. Unlike the low-redshift analogs, with a cosmological simulation, all
intrinsic galaxy properties are known. While all simulations are inevitably
subject to subgrid modelling, the advantage of simulations is that, with
sufficient resolution, they capture the complex interplay between gas
accretion and cooling, star formation, and feedback in 3D, all of which can
modify the morphology and spectrum of a galaxy, for example, through the star
formation history, production and distribution of dust, and the thermo-
chemical state of the ISM. Hence if the techniques used to infer galaxy
properties do not work on simulated noise-free mock data, it is doubtful that
they would generalize to real observations. Simulations are an ideal
laboratory for exploring the accuracy and underlying systematic biases of SED
codes and photoionization model grids.

  2. 2.

Because simulations attempt to model much of the physics of galaxy formation
from first principles, comparisons between simulations and observations can be
used to test our underlying theories of galaxy formation and to deduce which
physical processes become important in regulating galaxy growth and evolution
as a function of redshift.

In recent years there has been a substantial increase in the number of
simulations specifically targeted towards high-redshift galaxy formation.
These include Aurora (Pawlik et al., 2017), BlueTides (Waters et al., 2016),
CoDa (Ocvirk et al., 2016), CROC (Gnedin, 2014), Flares (Lovell et al., 2021),
Renaissance (O’Shea et al., 2015), Serra (Pallottini et al., 2019), Sphinx
(Rosdahl et al., 2018), Thesan (Kannan et al., 2022a), as well as many others.
In order to 1) employ simulations as a means of evaluating the accuracy of SED
fitting codes and photoionization model grids in inferring underlying galaxy
physical properties and 2) to test our theories of galaxy formation by
comparing simulated outputs with real observations, it is necessary to
completely forward-model mock observations from simulations. While this is
immediately true for the first proposed utility, for the latter it is less-
obviously the case. All simulations require subgrid models that are often
calibrated to reproduce some expected behaviour (see discussion in Section 2.1
of Schaye et al. 2015), for example, the globally inferred star formation rate
density (Madau & Dickinson, 2014), the stellar mass function (Li & White,
2009; Baldry et al., 2012), the stellar mass-halo mass relation (Moster, Naab
& White, 2018; Behroozi et al., 2019), the mass-metallicity relation (Tremonti
et al., 2004), or the UV luminosity function (Bouwens et al., 2015). The
crucial aspect of these calibrations is that the target quantities are often
inferred from observations rather than being the directly observed.

Strong systematic biases often exist between simulations and these
observational inferences. For example, the calibration between star formation
rate (SFR) and H$\alpha$ luminosity is sensitive to both the chosen stellar
initial mass function (IMF), underlying single stellar population (SSP) model,
and metallicity (e.g. Kennicutt & Evans, 2012; Eldridge et al., 2017). Stellar
masses are subject to the prior on SFR (e.g. Narayanan et al., 2023), while
metallicity is highly method-dependent (Andrews & Martini, 2013) and often
biased towards that of H II regions. Something as straightforward as a half-
light radius is subject to how the image is segmented (which sets the outer
radius), whether the quantity is circularized, if a Sérsic fit (Sérsic, 1963)
is used, whether the Sérsic fit is done in 1D or 2D, and whether the slope of
the Sérsic profile is fixed to a specific value (e.g. one). These definitions
not only differ between simulations and observers, but the method used to
report the same quantity in different observational papers often differs. In
order to be able to apply the same measurement techniques to both simulated
and observed data and make valid comparisons, it is paramount that simulations
be forward modelled into observational space.

To forward model the stellar continuum, dust radiative transfer is a common
technique (e.g. Jonsson, 2006; Baes et al., 2011), and this has widely been
applied to large simulations (e.g. Camps et al., 2018; Feldmann et al., 2023;
Yi et al., 2023). However, at high-redshift, nebular emission becomes an
important part of a galaxy SED (e.g. Schaerer & de Barros, 2010). While SSP
libraries are often shipped with emission line estimates (e.g. Byler et al.,
2017; Xiao, Stanway & Eldridge, 2018), the underlying assumptions of these
models (e.g. gas density, ionization parameter) are not obviously
representative of the physics of the high-redshift ISM. Hence, simply applying
these to simulated star particles may have limited fidelity for high-redshift
galaxies. A more accurate technique would be to try to simulate the high-
redshift ISM and its ionization state from first principles. This is
particularly important for simulations that model resonant emission lines such
as Ly$\alpha$ (Blaizot et al., 2023) or low ionization state absorption
features (Gazagnes et al., 2023). While radiative transfer is now included in
many simulations (e.g. Pawlik et al., 2017; Kannan et al., 2022a; Ocvirk et
al., 2016), the vast majority of full-box radiative hydrodynamics simulations
lack a model for the ISM.

Sphinx (Rosdahl et al., 2018) and its successor Sphinx20 (Rosdahl et al.,
2022) are notable exceptions as they simulate a statistically large sample of
galaxies in a full cosmological volume with radiative hydrodynamics while
simultaneously reaching resolutions high enough ($\approx 10$ physical pc) to
begin resolving and modelling the physics of a multiphase ISM. Moreover, the
volume of the simulation is large enough such that many of the galaxies are
bright enough to be observed with telescopes such as JWST, ALMA, and HST.
Hence, Sphinx is a unique resource where the full SEDs, including stellar
continuum, nebular continuum, and nebular emission lines, can be predicted in
a spatially-resolved manner and to a high-degree of accuracy for a large
sample of high-redshift galaxies. This mock data can be used to directly
compare with JWST data and as a laboratory for inferring galaxy properties
from mock observations.

In this work, we present Version 1 of the Sphinx Public Data Release (SPDRv1).
The data set is an open-sourced library of $\sim 14,000$ mock images and
spectra for galaxies in the redshift range $4.64\leq z\leq 10$ with 10 Myr-
averaged star formation rates $\geq 0.3\ {\rm M_{\odot}\ yr^{-1}}$,
specifically designed to be compared with JWST observations. We have
supplemented the data set with relevant intrinsic galaxy properties and made
the repository publicly accessible at the following URL:
https://github.com/HarleyKatz/SPHINX-20-data. Below we describe the methods
used to produce SPDRv1 and give a tour of what is available in the data set by
comparing the results with some of the early ERS and Cycle 1 JWST
observations.

##  2 Methods

In order to make direct comparisons between Sphinx and spectroscopic and
photometric observations of galaxies in the early Universe, we must forward
model both the intrinsic emission of the gas and stars in the simulated
galaxies and the propagation of these photons through the intervening ISM and
CGM. We consider both stellar and nebular (lines and continuum) emission as
well as attenuation and scattering by dust and/or H I. Below we describe the
details for generating each component of the mock spectra and images for each
Sphinx galaxy.

###  2.1 The Sphinx20 Simulations

Sphinx is a suite of cosmological radiation-hydrodynamic simulations with
state-of-the-art physical models and spatial/mass resolution designed to
simultaneously capture a multiphase ISM on the scale of galaxies, and the
reionization of the intergalactic medium (IGM) on large scales. The
simulations have been used to discern the physics of early galaxy formation
and the process of reionization, illustrating the sensitivity of the escape of
LyC photons from galaxies to stellar evolution models (Rosdahl et al., 2018),
the creation/amplification and role of magnetic fields (Attia et al., 2021;
Katz et al., 2021), how reionization suppresses the growth of dwarf galaxies
(Katz et al., 2020) and enhances the abundance of Lyman-alpha emitters (Garel
et al., 2021), the redshift-evolution of metal-line emission (Katz et al.,
2022b), how the LyC escape from galaxies correlates with observable galaxy
properties (Maji et al., 2022; Katz et al., 2022a; Choustikov et al., 2023),
and which galaxies predominantly power reionization (Rosdahl et al., 2022;
Katz et al., 2023b).

Sphinx20, which is exclusively used in this data release, is the largest
volume simulation in the Sphinx suite, with a volume of 203 co-moving Mpc3. It
contains thousands of star-forming galaxies at $z=6$. Here, we summarize the
most relevant characteristics of the simulation, but refer to Rosdahl et al.
(2022) and Rosdahl et al. (2018) for full details.

Sphinx20 was run over the period from February 2019 to November 2020 on the
JUWELS111https://www.fz-juelich.de/en/ias/jsc/systems/supercomputers/juwels
and Joliot-Curie ROME222https://www-hpc.cea.fr/en/Joliot-Curie.html
supercomputers and used 63 million CPU-hours to reach a final redshift of
4.64, when the volume is completely reionized. As for all Sphinx volumes, we
use the Ramses adaptive mesh refinement code (Teyssier, 2002) with full
radiation-hydrodynamics to model the propagation of LyC photons and their
interactions with hydrogen and helium (Rosdahl et al., 2018, 2015). Gas cells
are refined on mass and Jeans length, with a minimum cell width of 76 co-
moving pc, corresponding to 11 physical pc at $z=6$. The dark matter particle
mass is $2.5\times 10^{5}~{}{\rm M_{\odot}}$, allowing halos to be resolved
approximately down to the atomic cooling limit, and stellar particles,
representing stellar populations, have initial masses of $400~{}{\rm
M_{\odot}}$. Stars are formed out of the gas with a variable local efficiency
that depends on the thermo-turbulent properties of the gas. Supernova feedback
is performed using the mechanical model described in Kimm & Cen (2014); Kimm
et al. (2015), that adapts to the local conditions in order to generate a
consistent final momentum from the Sedov-Taylor phase. To obtain strong enough
suppression in star formation, we assume 4 SN explosions per 100 $~{}{\rm
M_{\odot}}$, which is significantly more than found with standard initial mass
functions observed in our local environment (but see Katz et al., 2022b). Gas
and stellar metallicity is tracked with a single scalar, and we describe below
how we extract individual metal species and ionization states of the gas in
post-processing. We use the BPASS version 2.2.1 SED model (Stanway & Eldridge,
2018) to calculate the injected LyC luminosities of stellar populations as a
function of their ages and metallicities. We use two LyC radiation groups on-
the-fly (13.6 eV-24.59 eV and 24.59 eV-$\infty$), but post-process the
simulation snapshots used in this data release with two additional lower-
energy radiation groups (5.6 eV-11.2 eV and 11.2 eV-13.6 eV) to more
accurately predict metal ionization states with ionization potentials lower
than H I. Non-equilibrium thermochemistry of hydrogen and helium is performed
on-the-fly, with the contribution of metals to cooling tabulated using Cloudy
for high temperatures and approximated with the fit from Rosen & Bregman
(1995) for low temperature. The cosmological initial conditions are generated
using MUSIC (Hahn & Abel, 2011) and have been selected to produce a typical
patch of the Universe in terms of the halo mass function. We use $\Lambda$CDM
cosmological parameters compatible with the Planck Collaboration et al.
(2014): $\Omega_{\rm m}=0.3175$, $\Omega_{\Lambda}=0.6825$, $\Omega_{\rm
b}=0.049$, $\sigma_{8}=0.83$, $H_{0}=67.11\ \ {\rm km\ s^{-1}\ Mpc^{-1}}$.

AdaptaHop (Tweed et al., 2009) is used for dark matter halo identification,
with halo finder parameters as described in Rosdahl et al. (2022). For each
halo, we designate its occupant galaxy to consist of all gas and stellar
particles enclosed within its virial radius.

In order to maintain manageable data sizes, we limit the current data release
to seven Sphinx20 snapshots, spanning $z=10$ to $z=4.64$, and only galaxies
which have star formation rates larger than $0.3\ \mbox{$\rm M_{\odot}\,{\rm
yr^{-1}}$}$, i.e. omitting the underlying population of galaxies likely too
dim to be observed with current state-of-the-art telescopes. Furthermore, we
only consider main haloes as all subhaloes will be contained within the virial
radius. We list in Table 1 the snapshots and the number of selected galaxies
in each, totaling to 1,380 galaxies.

Table 1: Statistics of the galaxies in the Sphinx data release. Each column shows the number of galaxies with ${\rm SFR}\geq 0.3\,\mbox{$\rm M_{\odot}\,{\rm yr^{-1}}$}$, their maximum and median stellar mass, and the maximum and median virial mass of their host dark matter halo, respectively. Redshift | $N_{\rm gal}$ | $\log\mbox{$M_{*}$}$ | $\log\mbox{$M_{*}$}$ | $\log M_{\rm vir}$ | $\log M_{\rm vir}$
---|---|---|---|---|---
|  | Max | Median | Max | Median
10.0 | 49 | 8.56 | 7.38 | 10.20 | 9.30
9.0 | 66 | 8.72 | 7.52 | 10.48 | 9.44
8.0 | 128 | 9.25 | 7.62 | 10.75 | 9.55
7.0 | 177 | 9.65 | 7.87 | 10.91 | 9.69
6.0 | 276 | 9.93 | 8.07 | 11.12 | 9.89
5.0 | 317 | 10.46 | 8.27 | 11.64 | 10.00
4.6 | 367 | 10.63 | 8.40 | 11.70 | 10.07

###  2.2 Stellar continuum generation

The intrinsic stellar continuum is generated for each halo by summing the
contributions of each star particle inside the virial radius of the halo. For
each star particle, we assign a spectrum based on its age, mass, and
metallicity by interpolating the model SEDs from BPASS v2.2.1 (Stanway &
Eldridge, 2018). We assume the same upper mass cutoff (100 ${\rm M_{\odot}}$)
and IMF slope ($-1.35$) as was used to calculate the ionizing luminosities for
star particles in the simulations.

###  2.3 Emission line generation

Intrinsic emission line luminosities are calculated for 52 emission lines (see
Table 2) for each gas cell within the virial radius of the halo following a
modification to the method described in Katz et al. (2022a) as detailed in
Choustikov et al. (2023). We first perform a search for all gas cells that
host star particles where the Stromgren sphere is unresolved (i.e. $R_{\rm
S}<\Delta x/2$). For these cells, we calculate emission line luminosities by
interpolating tables of CLOUDY models that were generated with v17.03 (Ferland
et al., 2017). The table contains models varying stellar age (computed as the
ionizing luminosity-weighted mean of all star particles within the gas cell),
metallicity, gas density, and total ionizing luminosity. These models assume a
spherical geometry (since the H II region is fully embedded within the host
gas cell). Because Sphinx does not follow the enrichment of individual
elements, metal abundances are assumed to follow Grevesse et al. (2010) scaled
by the metallicity of the cell. The models are iterated to convergence and
stopped at an electron fraction of 1%. For all other gas cells we assume that
the radiation and electrons are well mixed. We tabulate line emissivities from
a grid of one-zone slab models (again using CLOUDY) varying the gas density,
metallicity, ionization parameter, and electron fraction. All models are
iterated to convergence and the shape of the SED varies with metallicity, but
is assumed to have an age of $10$ Myr. Furthermore, for cells with resolved
ionization fractions, emission lines coming from H and He use the non-
equilibrium ionization fractions directly from the simulation to compute
luminosities. The total intrinsic emission of each halo is then the sum of of
the intrinsic line luminosities of all cells within the virial radius.

Table 2: Emission lines included in the Sphinx data release. Species are defined by their element symbol and ionization state. States with a suffix of “R” or “C” represent the recombination or charge exchange contribution to a collisionally excited line. Emission lines shown in magenta have been propagated through dust and/or H I radiative transfer. Wavelengths are the default values in cloudy v17. Species | State | Wavelength | Species | State | Wavelength
---|---|---|---|---|---
H | 1 | 1215.67 Å | O | 3 | 51.80 $\mu$m
H | 1 | 6562.80 Å | O | 3 | 88.33 $\mu$m
H | 1 | 4861.32 Å | Ne | 3 | 3868.76 Å
H | 1 | 4340.46 Å | Ne | 3 | 3967.47 Å
H | 1 | 4101.73 Å | C | 2 | 157.64 $\mu$m
He | 2 | 1640.41 Å | C | 3 | 1906.68 Å
He | 2 | 4685.68 Å | C | 3 | 1908.73 Å
O | 1 | 6300.30 Å | C | 4 | 1548.19 Å
O | 1 | 6363.78 Å | C | 4 | 1550.78 Å
O | 2 | 3726.03 Å | N | 2 | 5754.61 Å
O | 2 | 3728.81 Å | N | 2R | 5755.00 Å
O | 2R | 3726.00 Å | N | 2 | 6548.05 Å
O | 2R | 3729.00 Å | N | 2 | 6583.45 Å
O | 2 | 7318.92 Å | N | 2R | 6584.00 Å
O | 2 | 7319.99 Å | N | 3 | 1748.65 Å
O | 2 | 7329.67 Å | N | 3 | 1753.99 Å
O | 2 | 7330.73 Å | N | 3 | 1746.82 Å
O | 2R | 7332.00 Å | N | 3 | 1752.16 Å
O | 2R | 7323.00 Å | N | 3 | 1749.67 Å
O | 3 | 1660.81 Å | S | 2 | 6716.44 Å
O | 3 | 1666.15 Å | S | 2 | 6730.82 Å
O | 3 | 4363.21 Å | S | 2 | 4076.35 Å
O | 3R | 4363.00 Å | S | 2 | 4068.60 Å
O | 3C | 4363.00 Å | S | 3 | 6312.06 Å
O | 3 | 4958.91 Å | S | 3 | 9068.62 Å
O | 3 | 5006.84 Å | S | 3 | 9530.62 Å

###  2.4 Nebular continuum generation

The final component of the spectrum is the nebular continuum. Here we again
split cells based on whether or not they host an unresolved Stromgren sphere.
In the case where there is an unresolved Stromgren sphere, we again
interpolate tabulated CLOUDY data for the nebular continuum (free-free, free-
bound, two-photon) of only H and He. For all other cells we compute the free-
free, free-bound, and two-photon emission for H and He individually using the
non-equilibrium ionization fractions and electron densities in the cell. Free-
free emission is computed following the method described in Schirmer (2016)
using Gaunt factors from van Hoof et al. (2014). Free-bound emission is
computed with CHIANTI (Dere et al., 2019). Finally, two-photon emission is
also computed following Schirmer (2016).

###  2.5 Dust and Resonant Line Radiative Transfer

All line and continuum radiation is subject to absorption and scattering by
dust. The Sphinx public data release also contains bright resonant lines such
as Ly$\alpha$ that must diffuse both spatially and in frequency space to
escape the galaxy and be viewed by an observer. All of this physics is
modelled using the Monte-Carlo radiative transfer code Rascas (Michel-Dansac
et al., 2020) similarly to the method described in Katz et al. (2022a). Dust
is assigned to all gas cells using the Small Magellanic Cloud (SMC) dust model
described in Laursen, Sommer-Larsen & Andersen (2009). This phenomenological
dust model is normalized such that the dust absorption coefficient in each gas
cell is given by $(n_{\rm HI}+f_{\rm ion}n_{\rm HII})\sigma_{\rm
dust}Z/Z_{0}$, where $n_{\rm HI}$ and $n_{\rm HII}$ are the number densities
of neutral and ionized hydrogen, respectively, the fraction of ionized gas
contributing to the dust density is assumed to be $f_{\rm ion}=0.01$,
$\sigma_{\rm dust}$ is the dust cross section, $Z$ is the metal mass fraction,
and $Z_{0}=0.005$. The dust asymmetry and albedo are adopted from the SMC dust
model of Weingartner & Draine (2001).

For the emission lines shown in magenta in Table 2, 200,000 photon packets are
probabilistically distributed to all gas cells of a halo based on their
intrinsic luminosities. The initial wavelengths of the photon packets are
placed at line-centre of the transition, thermally broadened, and shifted to
the observer’s frame based on the bulk velocity of the cell. Escape fractions
and images for all other emission lines are interpolated from the closest ion
for which we run the radiative transfer due to computational constraints. For
the nebular and stellar continuum, we sample the spectrum at 20 different
wavelengths between 1300 Å and 6583 Å. The wavelengths are chosen to be
consistent with the locations of strong emission lines (so that equivalent
widths can be accurately measured) and consistently throughout the UV where
dust attenuation is the strongest. Another 200,000 photon packets are
probabilistically distributed to each star particle based on the luminosity of
the stellar continuum at each wavelength. Likewise, the same number of packets
are probabilistically distributed for each wavelength to gas cells based on
the nebular continuum luminosity. The full spectrum of each galaxy can thus be
constructed by computing the escape fractions of each component and
interpolating across wavelength.

Due to the H I opacity, Ly$\alpha$ is computed slightly differently to the
method described above. Ly$\alpha$ photon packets are probabilistically
assigned to gas cells based on their intrinsic luminosity, which is the sum of
the recombination and collisional components. We use $10^{6}$ photon packets
to better capture the spectrum when the escape fraction is low. The Ly$\alpha$
photon packets are propagated in a medium that consists of H I, D, and dust.
The dust model is the same as used for all of our other RASCAS runs as
described above and we assume a fixed D/H abundance of $3\times 10^{-5}$.

To produce mock observations along individual sight lines, we employ the
peeling algorithm (Yusef-Zadeh, Morris & White, 1984; Zheng & Miralda-Escudé,
2002; Dijkstra, 2017) to make images and spectra along ten directions. The
viewing angles are kept fixed for every galaxy. As our data set contains 1,380
galaxies, viewing each along ten sight lines results in 13,800 total spectra
(not including the 1,380 intrinsic spectra).

###  2.6 Available Data

In Table 3 we list the full set of galaxy properties, mock observations, and
derived quantities that can be downloaded from the public repository. The data
set contains numerous intrinsic properties of the galaxies such as stellar
mass, halo mass, SFR averaged over multiple time scales, SFHs, etc. as well as
the full spectra described above. We additionally provide certain derived
quantities such as UV continuum slopes, UV magnitudes, and JWST filter
magnitudes. The breadth of data should encapsulate the vast majority of JWST
observations available at present.

Table 3: Details of the galaxy properties available as part of the Sphinx public data release. Quantity |  Units |  Notes
---|---|---
Halo ID |  |
Redshift |  |
Halo mass |  ${\rm\log_{10}(M/M_{\odot})}$ |
Stellar mass |  ${\rm\log_{10}(M/M_{\odot})}$ |  This value is the total stellar mass formed (i.e. the integral of the star formation history) and is not adjusted for mass loss due to stellar feedback
$R_{\rm vir}$ |  20 cMpc |
$x,\ y,\ z$ position |  20 cMpc |  3D position of the halo in the simulation volume
Star formation rate |  ${\rm M_{\odot}\ yr^{-1}}$ |  Provided as an average over 3, 5, 10, and 100 Myr and can be recomputed for any other interval from the star formation history
Star formation history |  ${\rm M_{\odot}\ yr^{-1}}$ |  Provided for every galaxy on a 1 Myr time cadence
Stellar ages |  Myr |  Mass-weighted and LyC luminosity-weighted
Stellar metallicity |  Absolute |  Mass-weighted and LyC luminosity-weighted over all stars
Stellar metallicity history |  Absolute |  Mass-weighted stellar metallicity of all star particles that formed in bins of 1 Myr
Ionizing luminosity |  photons ${\rm s^{-1}}$ |
LyC escape fraction |  |  Angle-averaged (for all photons with $E>13.6\ {\rm eV}$) and along ten sight lines (for photons with a wavelength of 900 Å)
ISM gas density |  $\log_{10}({\rm n_{H}/cm^{-3}})$ |  Weighted by intrinsic [O II] $\lambda\lambda$3727 or [C III] $\lambda\lambda$1908
Gas metallicity |  $\log_{10}(Z/Z_{\odot})$ |  Mass-weighted as well as [O II] $\lambda\lambda$3727, [O III] $\lambda$5007, [N II] $\lambda$6583, and H$\beta$ weighted
Emission line luminosities |  ${\rm erg\ s^{-1}}$ |  Intrinsic for all emission lines listed in Table 2, dust attenuated along ten sight lines for H$\alpha$, H$\beta$, H$\gamma$, H$\delta$, [O II] $\lambda\lambda$3727, [Ne III] $\lambda$3869, [O III] $\lambda$4363, [O III] $\lambda$4959, [O III] $\lambda$5007, [N II] $\lambda$6583
Stellar continuum luminosities |  ${\rm erg\ s^{-1}\ \text{\AA}^{-1}}$ |  Intrinsic & dust attenuated along ten sight lines for 20 wavelengths (1300Å, 1400Å, 1500Å, 1600Å, 1700Å, 1800Å, 1900Å, 2000Å, 2500Å, 3000Å, 3727Å, 3869Å, 4102Å, 4341Å, 4363Å, 4861Å, 4959Å, 5008Å, 6563Å, 6583Å)
Nebular continuum luminosities |  ${\rm erg\ s^{-1}\ \text{\AA}^{-1}}$ |  Intrinsic & dust attenuated along ten sight lines for 20 wavelengths (1300Å, 1400Å, 1500Å, 1600Å, 1700Å, 1800Å, 1900Å, 2000Å, 2500Å, 3000Å, 3727Å, 3869Å, 4102Å, 4341Å, 4363Å, 4861Å, 4959Å, 5008Å, 6563Å, 6583Å)
Full SEDs |  ${\rm erg\ s^{-1}\ Hz^{-1}\ cm^{-2}}$ |  Intrinsic & dust attenuated along ten sight lines and redshifted to the relevant $z$. Spectra are computed at 1 Å resolution by interpolating the escape fractions at the 20 continuum wavelengths and for each emission line. SED files provide the total SED as well as the three separate components
E(B$-$V) |  |  Along ten sight lines. Computed from the Balmer decrement (H$\alpha$ and H$\beta$)
Effective radii ($R_{\rm eff}$) |  pc |  Measured in each of the JWST filters along each line of sight for the largest segment after our image segmentation procedure. We provide the corresponding flux density (nJy) of the segment in addition to its circularized radius
UV continuum slopes ($\beta$) |  |  Intrinsic & dust attenuated along ten sight lines. Measured from the full SED (stellar + nebular continuum) as well as only the stellar continuum. Additional values can be measured from the photometry with the inclusion of emission lines
UV magnitudes |  AB |  Intrinsic & dust attenuated along ten sight lines. Measured at 1500 Å from the stellar and nebular continuum
JWST filter magnitudes |  AB |  Dust attenuated along ten sight lines. Computed for all NIRCam wide and medium filters (F070W, F090W, F115W, F140M, F150W, F162M, F182M, F200W, F210M, F250M, F277W, F300M, F335M, F356W, F360M, F410M, F430M, F444W, F460M, F480M)
Ly$\alpha$ and H$\alpha$ spectra |  ${\rm erg\ s^{-1}}$ |  Dust attenuated along ten sight lines. Spectral resolution of 0.1 Å. Values should be divided by the wavelength bins to obtain appropriate units
Ly$\alpha$ and H$\alpha$ surface brightness profiles |  ${\rm erg\ s^{-1}}$ |  Dust attenuated along ten sight lines. Spatial resolution of $R_{\rm vir}/250$. Values should be divided by the pixel size to obtain surface brightness
Galaxy images |  ${\rm nJy\ pixel}$ |  Dust attenuated along ten sight lines for each JWST filter. Due to data size, these are made available upon request for any emission line or continuum (nebular or stellar) wavelength. Example RGB images combining multiple filters are shown in Figure 1

##  3 Results

The primary purpose of SPDRv1 is to help better understand the physical
properties of high-redshift galaxies that are (potentially) observable with
existing facilities. Thus, we have limited our analysis to actively star-
forming galaxies with a 10 Myr-averaged SFR $\geq 0.3\,\mbox{$\rm
M_{\odot}\,{\rm yr^{-1}}$}$. Below we provide a tour of data set, organized
by:

  * •

Intrinsic galaxy properties (e.g. stellar and halo masses, metallicities,
etc.),

  * •

Photometric properties (i.e. those that can be compared to observational data
sets with imaging), and

  * •

Spectroscopic properties (e.g. full SEDs and emission lines).

Our aim is to demonstrate the diversity of what is available as part of the
data release as well as what can be computed. We provide a few sample
comparisons with real observations to contextualize the data release among
state-of-the-art observational data.

Figure 1: RGB images of the Sphinx galaxies with the highest escaping LyC
luminosity (subject to a stellar mass threshold) at each redshift. F365W,
F444W, and F410M are used for the red channel, F200W, and F277W are used for
green, and F115W and F150W are used for blue, identical to the UNCOVER mosaic
(Bezanson et al., 2022). As redshift changes, the filters sample a different
wavelength range of the SEDs of the galaxies and hence they appear to change
colour.

To begin, in Figure 1 we show a gallery of mock RGB images of 100 Sphinx
galaxies in the redshift interval $6\leq z\leq 10$. Following the UNCOVER
mosaic (Bezanson et al., 2022), we combine F365W, F444W, and F410M in the red
channel, F200W, and F277W in green, and F115W and F150W in blue. To highlight
the lower surface brightness features, we scale the pixel intensities by an
arcsinh filter. Galaxies are selected as those with the highest escaping LyC
luminosities (i.e. the galaxies that are currently reionizing the simulated
volume) subject to a stellar mass threshold. As redshift changes, the filters
sample a different wavelength range of the SEDs of the galaxies and hence they
appear to change colour. The morphological diversity exemplifies how the
galaxies leaking LyC radiation can look very different. This is true not only
between different galaxies, but also for the same galaxy along different lines
of sight. However, in the vast majority of galaxies, the ISM is clearly
disrupted and there are strong outflows, allowing for the escape of LyC
radiation.

Figure 2: Example spectrum of a massive $z=6$ galaxy in the Sphinx public data
release. We show the intrinsic stellar continuum (olive), intrinsic nebular
continuum (red), and total intrinsic spectrum (black). NIRCam filter
magnitudes are shown as the orange data points. The different blue lines
represent the dust-attenuated spectra along the ten different sight lines in
the data set.

In Figure 2 we show the full intrinsic spectrum (black) broken down into
stellar continuum, nebular continuum, and nebular emission lines, as well as
the dust-attenuated spectra along the ten sight lines for the most massive
$z=6$ galaxy. Data points in this image represent the JWST filter magnitudes
for the 20 wide and medium filters on NIRCam. Data such as that shown in
Figures 1 and 2 is available for every galaxy in SPDRv1.

Figure 3: (Left) Dark matter halo mass functions for haloes that host star-
forming galaxies with $\rm SFR_{\rm 10}\geq 0.3\,M_{\odot}/yr$ in Sphinx for
different redshifts (solid lines). For comparison, the mass functions for the
entire sample of main halo population are shown as dashed lines, while the
theoretical estimates by Despali et al. (2016) are shown as dotted lines.
(Right) Stellar mass distributions of our main sample with $\rm SFR_{\rm
10}\geq 0.3\,M_{\odot}/yr$ for different redshifts (solid lines). The stellar
mass distributions for all Sphinx galaxies are shown as dashed lines. Note
that ${\rm M_{*}}$ corresponds to the total stellar mass formed, so does not
take into account the mass loss due to stellar evolution.

###  3.1 Intrinsic Galaxy Properties

####  3.1.1 Halo masses

The initial conditions for Sphinx were chosen to simulate an average region of
the Universe in which the total ionizing luminosity is closest to its mean
value (e.g., see Figure 1 of Rosdahl et al., 2018). The mass function of dark
matter halos in Sphinx20 is shown in Figure 3, at different redshifts in the
interval $4.6<z<10$. It is in good agreement with the theoretical expectations
from dark matter-only N-body simulations (e.g., Despali et al., 2016) despite
the limited volume of the simulation.

Sphinx has 278, 1,069, 2,789, 4,508 DM halos with masses above $10^{9}~{}{\rm
M_{\odot}}$ at $z=10$, $8$, $6$, and $4.6$, respectively. Our data release
includes galaxies with $\mbox{${\rm SFR}_{\rm 10}$}\geq 0.3\,\mbox{$\rm
M_{\odot}\,{\rm yr^{-1}}$}$ that are mostly embedded in halos with
$\mbox{$M_{\rm vir}$}\sim 10^{9}~{}{\rm M_{\odot}}$ at $z=10$ or halos with
$\mbox{$M_{\rm vir}$}\sim 10^{10}~{}{\rm M_{\odot}}$ at $z=4.6$, but a small
fraction of them is also found in halos with $10^{8}<\mbox{$M_{\rm
vir}$}<10^{9}~{}{\rm M_{\odot}}$, as star formation can be bursty and
temporarily enhanced. This effect is important as analytic models often assume
that SFR scales with the accretion rate of dark matter and the scatter may be
important for explaining the number of bright high-redshift galaxies (e.g.
Shen et al., 2023; Sun et al., 2023).

####  3.1.2 Stellar masses

In the right panel of Figure 3, we show the distributions of stellar masses
formed (${\rm M_{*,i}}$) in Sphinx20galaxies. ${\rm M_{*,i}}$ is computed as
the sum of the initial mass of all star particles within the virial radius
(i.e. including subhaloes). This is equivalent to computing the integral of
the SFH of each galaxy and neglects mass loss due to stellar evolution.
Because we consider all star particles within the virial radius, the data
release excludes subhaloes to avoid double counting.

Within the $20^{3}\,{\rm cMpc}^{3}$ volume, Sphinx20 has 2,417, 6,099, 11,747,
and 15,455 galaxies with ${\rm M_{\rm*,i}>10^{6}\ M_{\odot}}$ at $z=10$, $8$,
$6$, and $4.64$, respectively. Of these, only 49, 128, 276, and 367 galaxies
with ${\rm SFR}_{\rm 10}\geq 0.3\,\mbox{$\rm M_{\odot}\,{\rm yr^{-1}}$}$ are
selected to simulate emission lines (see Table 1 for more details). The
typical galactic stellar masses we probe with the SFR criterion range between
$7\la\log_{10}({\rm M_{\rm*,i}}/{\rm M_{\odot}})\la 9$ at $z>7$, which are
comparable to those from the early JWST results (e.g., Santini et al., 2023).
More massive galaxies with stellar masses up to $4\times 10^{10}\ {\rm
M_{\odot}}$ appear in later epochs ($z=4.64$), providing the opportunity to
learn more about the physical nature of Ly$\alpha$ emitters, such as those
observed with MUSE (e.g. Bacon et al., 2017).

Our choice of SFR cut is motivated by the recent analysis from JWST
observations (Leethochawalit et al., 2023; Fujimoto et al., 2023; Shapley et
al., 2023) that are subject to a flux limit. This of course ignores the fact
that UV luminosity and SFR do not necessarily map one-to-one (primarily
because of dust attenuation); however, we cannot determine the UV luminosity
until after post-processing (which is the computationally limiting step). Our
cut excludes a large fraction of galaxies with ${\rm M_{*,i}}\la 10^{8}\ {\rm
M_{\odot}}$. For example, our galaxy samples are 90% complete in stellar mass
down to $\log_{10}({\rm M_{*,i}}/[{\rm M_{\odot}}])=8.0$, $8.5$, $8.9$, and
$8.8$ at $z=10$, $8$, $6$, and $4.64$, respectively, and thus are biased (by
construction) towards star-forming galaxies. This should be kept in mind in
the following sections when we compare Sphinx data to observations.

####  3.1.3 Stellar mass - Halo mass relation

Although neither quantity is actually observable, the stellar mass-halo mass
relation provides important information on the galaxy-halo connection at all
redshifts and constraints from abundance matching are often used to tune
subgrid feedback models in simulations. This relation has been presented for
Sphinx data for all galaxies in Rosdahl et al. (2018, 2022) and is available
for star-forming galaxies as part of the Sphinx data release. We show our
results in Figure 4. For comparison we also plot the stellar mass-halo mass
relation from the model of Tacchella et al. (2018) as well as from abundance
matching (Behroozi et al., 2019), both of which are estimated at $z=6$.

Figure 4: Stellar mass-halo mass relation for galaxies in the Sphinx data
release compared to the $z=6$ model of Tacchella et al. (2018) and $z=6$
estimates from abundance matching Behroozi et al. (2019).

We note a few important points regarding high-redshift stellar mass. Despite
its ubiquity as a constraint on feedback, stellar mass at these redshifts is
not an observed quantity; hence, constraints from observational data are
model-dependent. This is particularly true at high-redshift where the JWST
photometry primarily probes the rest-frame UV and shorter-wavelength optical
and HST data probes only rest-frame UV. These regions of the spectrum are
dominated by young stellar populations, especially if the galaxies exhibit
bursty star formation histories and thus the prior on the star formation
history is a key parameter in determining the final stellar mass (in addition
to stellar IMF and various other quantities). This is reflected by the fact
that spatially resolved stellar mass estimates at high-redshift often result
in stellar masses that are $\sim 0.5-1$ dex higher than the galaxy-integrated
photometry (Giménez-Arteaga et al., 2023). This offset is consistent with the
difference between the Sphinx prediction and that from abundance matching. The
Sphinx data release can perhaps play a role in improving the priors used in
SED fitting codes (see also Narayanan et al. 2023). Finally the models of
Tacchella et al. (2018); Behroozi et al. (2019), which have been used as
benchmarks for theoretical models of high-redshift galaxy formation, do not
reproduce the cumulative number counts of high-redshift galaxies from JWST
(see e.g. Leung et al. 2023 and Section 4.1.3), neither do the simulations
that have been tuned to reproduce them (e.g. Kannan et al., 2022a).

####  3.1.4 Star formation rates

Figure 5: Star formation rates averaged over 100 Myr ($\rm SFR_{100}$, top), 5
Myr ($\rm SFR_{5}$, middle) and 3 Myr ($\rm SFR_{3}$, bottom) as a function of
10 Myr-averaged SFR ($\rm SFR_{10}$). Individual points represent Sphinx
galaxies at different redshifts, selected to have $\rm SFR_{\rm 10}\geq
0.3\,M_{\odot}\ yr^{-1}$. The black diagonal lines show the one-to-one
relations.

Observationally, different tracers exist for determining the SFR of a galaxy
such as UV and IR luminosity as well as various emission lines across the
wavelength spectrum (e.g. Kennicutt, 1998; Kennicutt & Evans, 2012; De Looze
et al., 2014). Different SFR indicators probe star formation on various
timescales. For instance, the ${\rm H}\alpha$ recombination line results from
gas ionised by young stars, and hence gives a hint of the current SFR on
$\lesssim 10$ Myr time scales. In contrast, FUV continuum emission traces a
longer-time scale SFR, up to 100 Myr. Therefore, the SFR from a galaxy
spectrum may differ depending on the SFR indicator that is used. In order to
illustrate the consequences of the star formation history variability, we
provide in SPDRv1 SFRs measured over multiple time scales.

Figure 5 shows SFR averaged over 100, 5, and 3 Myr (respectively $\rm
SFR_{100}$, $\rm SFR_{5}$, and $\rm SFR_{3}$) as a function of 10 Myr-averaged
SFR ($\rm SFR_{10}$). All the galaxies in the data release are selected to
have $\rm SFR_{\rm 10}\geq 0.3\ M_{\odot}\ yr^{-1}$, but on different time
scales, the SFR can drop well below this value. For example, on average, $\rm
SFR_{100}$ tends to fall well below $\rm SFR_{10}$, indicating that the 10
Myr-averaged SFR is more representative of a burst rather than the continuous
star formation history. In contrast, $\rm SFR_{5}$ and $\rm SFR_{3}$ are more
evenly scattered around the one-to-one relation. However, there are a few
galaxies with very low $\rm SFR_{5}$ and $\rm SFR_{3}$, which is indicative of
a short-timescale quenching event.

####  3.1.5 Star formation main sequence

There is a well established relationship between the stellar mass of a galaxy
and its SFR (e.g. Brinchmann et al., 2004; Salim et al., 2007) that extends to
the epoch of reionization (e.g. Popesso et al., 2023). Sphinx galaxies
similarly exhibit such a relationship. In Figure 6 we show both the 10 Myr-
averaged (top) and 100 Myr-averaged (bottom) SFRs as a function of stellar
mass. The SFRs in the Sphinx data release reach as high as 100 ${\rm
M_{\odot}\ yr^{-1}}$ and span more than three dex in SFR on longer timescales.
On 100 Myr timescales there is a tight correlation such that specific SFRs
(sSFRs) typically fall between $10^{-8}-10^{-9}$ yr-1 while on 10 Myr
timescales sSFRs can go as high as $10^{-7}$ yr-1. Such high values are
consistent with some of the most highly star-forming extreme emission line
galaxies at lower redshift (e.g. Amorín et al., 2017) as well as tentative
estimates from early JWST observations (e.g. Shapley et al., 2023; Curti et
al., 2023). The high sSFRs are clearly not sustained over long timescales due
to the fact that the 100 Myr-averaged SFRs are in general lower than the 10
Myr-averaged values for our most highly star-forming galaxies.

We note that the flattening in the main sequence at the low-mass end of the 10
Myr-averaged figure is entirely due to the threshold SFR used to create the
Sphinx public data set. If all galaxies were included, our relations would
trend towards lower SFR, albeit the relation would still exhibit significantly
more scatter than the 100 Myr-averaged main-sequence. This is an important
point to consider, and a recurring theme throughout this work that higher
redshift galaxies in SPDRv1 tend to probe galaxies further above the main
sequence. In many ways this is reflective of the behaviour of a flux-limited
survey where the fainter galaxies at higher-redshifts must fall consistently
higher above the main-sequence to be observed.

Figure 6: Star formation main sequence in Sphinx for different redshifts for
10 Myr and 100 Myr averaged SFRs (upper and lower panel, respectively). The
points represent individual galaxies and the solid lines represent median star
formation rates. The alignment of points at low SFR is due to the formation of
integer numbers of star particles.

####  3.1.6 Star formation histories

The Sphinx data release contains the star formation histories (SFHs) of all
sampled halos, with a time resolution of 1 Myr. We show in Figure 7 the SFHs
of five example galaxies in the $z=4.64$ snapshot, selected to span the range
of galaxy masses in our sample, as indicated in the plot legend. As is typical
for all Sphinx galaxies, the SFHs shown are very stochastic and bursty,
increasingly so for decreasing galaxy stellar masses. The plot also
demonstrates that the lower-mass galaxies have spent most of their lives with
star formation rates below the selection limit of the data release sample,
i.e. with ${\rm SFR\geq 0.3\ M_{\odot}\ yr^{-1}}$, and therefore would
approximately phase in and out of being observable with the JWST. One
important property of these SFHs is that they represent star formation in all
of the progenitor haloes, rather than the main branch along the merger tree.
This is important for SED modelling as the total star formation in all
progenitors is the quantity constrained by the integrated SED.

Figure 7: Star formation histories (top) and stellar metallicity histories
(bottom) of five Sphinx galaxies, with stellar masses at $z=4.64$ as indicated
in the legend.

Figure 8: (Top) SFRD of the Sphinx volume, using the galaxies selected for the
data release, i.e. traced back in time for all galaxies with ${\rm
SFR_{10}\geq 0.3\ \mbox{$\rm M_{\odot}\,{\rm yr^{-1}}$}}$ from Sphinx
snapshots at the redshifts indicated in the legend. For comparison we show the
dust-corrected estimates for the SFRD from Bouwens et al. (2015). (Bottom)
SFRD of the Sphinx volume measured from the UV magnitude of each Sphinx galaxy
at each redshift using either the intrinsic ${\rm M_{UV}}$ (yellow) or dust
obscured ${\rm M_{UV}}$ (red) down to a limiting magnitude of $-18$. For
comparison we show observational constraints from Harikane et al. (2023);
Bouwens et al. (2023) that have been integrated to the same limiting
magnitude. We have converted their reported UV luminosity density into an SFRD
using our Equation 1.

We can also extract the SFR density of the early universe from the Sphinx
data. This is shown in the top panel of Figure 8, where we plot the SFR
density traced back in time for all sampled galaxies (i.e. with ${\rm SFR\geq
0.3\ \mbox{$\rm M_{\odot}\,{\rm yr^{-1}}$}}$) at the given redshifts. Note how
the total SFRD for galaxies in our $z=4.64$ catalogue always tracks above
those in the higher redshift catalogues. This is due to the fact that the low
mass progenitors that do not reach our SFR threshold merge into larger haloes
by $z=4.64$. Thus their star formation is included in the history of the
massive galaxies in the lower-redshift snapshots but not at high redshift.

This intrinsic SFRD is, however, not directly comparable to observations,
which instead measure a UV luminosity density ($\rho_{\rm SFR}$) above some
magnitude threshold that is then converted into an SFRD. To demonstrate this
comparison, we first derive an empirical relation between absolute magnitude
at 1500 Å, MUV,1500, and 10 Myr-averaged SFR from our data set such that

$\log_{10}\left(\frac{\rm SFR_{10}}{\rm M_{\odot}\ yr^{-1}}\right)=-0.35\,{\rm
M_{UV,1500}}-6.49.$ (1)

We then select all galaxies that have an intrinsic M${}_{\rm UV}\leq-18$ and
calculate their total SFRD at each snapshot. This is shown as the yellow line
in the bottom panel of Figure 8. We then use the observed MUV values along
each of our ten sight lines (i.e. those uncorrected for dust) and convert to
an SFRD using Equation 1. The red line and shaded region in the bottom panel
of Figure 8 represent the SFRD that an observer would derive for our sample at
each snapshot for galaxies with observed M${}_{\rm UV}\leq-18$. At $z=4.64$,
the observed and true SFRDs differ by an order of magnitude due to obscuration
by dust. For comparison, we have converted $\rho_{\rm UV}$ for
spectroscopically confirmed high-redshift JWST galaxies with M${}_{\rm
UV}\leq-18$ from Harikane et al. (2023) to an SFRD using our Equation 1 (for
consistency). Because these measurements do not correct for dust and are only
spectroscopically confirmed galaxies, they constitute a firm lower-limit on
the SFRD. These are shown at the light blue arrows in the bottom panel of
Figure 8 and they are surprisingly consistent with our mock observations. We
emphasize that this agreement is to some extent serendipitous because Sphinx
does not probe galaxies as bright as seen in Harikane et al. (2023). For
example, the brightest dust-obscured $z=10$ galaxy in Sphinx has a UV
magnitude of -19.6 while Harikane et al. (2023) probe galaxies brighter than
$-21$ at this redshift. Selecting only fainter galaxies from the Harikane et
al. (2023) sample (i.e. those with ${\rm M_{UV}>-19.6}$) inevitably decreases
the SFRD as shown by the dark blue arrows. We have performed the same
experiment using data from Bouwens et al. (2023) which has a similar magnitude
limit and the results differ compared to the Harikane et al. (2023) data. This
may simply be due to cosmic variance or different estimates of survey volume.
Nevertheless, the key feature is that even at such high redshifts, dust plays
an important role in obscuring star formation and the Sphinx data set provides
a means to directly compare SFRD estimates between simulations and
observations.

####  3.1.7 Stellar ages

The stellar age of a galaxy is an important parameter that encapsulates its
mass growth history. Mass-weighted stellar ages are provided in SPDRv1 and can
likewise be recomputed from the SFHs. Figure 9 shows the mass-weighted stellar
ages as a function of stellar mass. Interestingly, at ${\rm M_{*}\gtrsim
10^{7.5}\ M_{\odot}}$ the trend is relatively flat at a fixed redshift such
that galaxies, independent of their stellar mass, exhibit similar stellar
ages. However, as redshift increases, ages of galaxies decrease. At lower
stellar masses, the galaxies tend to exhibit younger stellar ages but indeed
this is due to the SFR threshold in the database. At such low stellar masses,
a 10 Myr-averaged SFR of $0.3\ {\rm M_{\odot}\ yr^{-1}}$ represents a stellar
mass of $3\times 10^{6}\ {\rm M_{\odot}}$, which is a significant fraction of
the total stellar mass of objects with ${\rm M_{*}<10^{7.5}\ M_{\odot}}$

Figure 9: Mass-weighted stellar age in Sphinx for different redshifts as a
function of galaxy stellar mass. The solid lines and large data points
indicate the mean values while smaller data points show individual galaxies.

Stellar age need not be defined as a mass-weighted quantity. At high-redshift,
JWST is primarily sensitive to the rest-frame UV and optical parts of a
spectrum. Furthermore, the emission from nebular regions cares only about the
stars that emit LyC photons. For this reason, we also provide stellar ages
weighted by LyC emission. A comparison between the mass-weighted and LyC-
weighted ages is shown in Figure 10. While the mass-weighted ages reach values
of hundreds of Myr and represent a significant fraction of the Hubble time at
each particular redshift, the LyC-weighted ages are only sensitive to the
young stellar populations and hence probe the recent star formation in the
past $\sim 10$ Myr.

Figure 10: LyC-weighted stellar age versus mass-weighted stellar age for
Sphinx galaxies in the data release.

####  3.1.8 Stellar metallicity histories

In addition to the SFHs, we also provide stellar metallicity histories for
each galaxy in the Sphinx data release. These are computed as the mass-
weighted metallicity of all star particles that form in the same 1 Myr bins as
were used to compute the SFHs. The bottom panel of Figure 7 shows the stellar
metallicity evolution for the same five galaxies as in the top panel. Note
that the stellar metallicity functions are not monotonic. Pristine inflows can
dilute the gas metallicity of the galaxies which is then imprinted on the
stars. Similarly, the metallicity evolution contains the star formation events
in all progenitors with different metallicity evolution and SFHs. Similar to
the SFHs, the metallicity evolution is a key component of the SED modelling,
and thus the Sphinx data can be used as a prior when analyzing real
observations.

####  3.1.9 Mass-metallicity relation

The stellar mass/gas-phase metallicity relation (MZR) is one of the primary
galaxy scaling relations that encapsulates information on the baryon cycle.
Hence, constraining this relation is one of the primary goals of high-redshift
JWST observations (e.g. Nakajima et al., 2023; Curti et al., 2023). The Sphinx
data release contains all of the information needed to compare simulated
metallicity to that which can be probed with observations.

In the top panel of Figure 11 we show the MZR where metallicity is weighted by
gas mass within each halo. There is little evolution with redshift among our
star-forming sample. This lack of evolution is partially driven by the fact
that the SFR cut means that the higher redshift galaxies tend to fall further
above the star formation main-sequence. However we emphasize that while the
mass-weighted metallicity is the value predicted by most cosmological
simulations, this is not a quantity that is accessible to observers. Rather,
observations are only sensitive to metallicity in the bright, line-emitting H
II regions. In the middle panel we show the metallicity weighted by the [O
III] $\lambda$5007, [O II] $\lambda\lambda$3727, and H$\beta$ luminosities of
the gas cells, which better represents the metallicity of H II regions. These
values are significantly higher than the mass-weighted metallicity of the
galaxy demonstrating the early self-enrichment of galaxies.

Because the data release contains all of the necessary emission lines to apply
the direct method, which is widely considered the gold-standard metallicity
estimator, in the bottom panel of Figure 11 we show the MZR where metallicity
is computed along each sight line. Specifically, we use the reddening
corrected values of [O III] $\lambda$5007, [O III] $\lambda$4363, [O II]
$\lambda\lambda$3727, and H$\beta$ and use the Pilyugin et al. (2009)
conversion between tO3 and tO2333tO3 and tO2 correspond to the electron
temperature in the high and low ionization zones, respectively.. We reproduce
the result from Cameron, Katz & Rey (2022) where H II region metallicity is
underpredicted by the direct method due to temperature fluctuations between H
II regions across the simulated galaxies. The Sphinx data release contains
auroral lines for many elements and ionization states so that the direct
method metallicities can be compared between Sphinx and observations for
numerous JWST programs (e.g. ID 2953 - CECILIA - PIs Strom, Rudie, ID 1914 -
AURORA - PIs Shapley, Sanders, ID 1879 - PI Curti).

Figure 11: Stellar mass/gas-phase metallicity relation for SPHINX20 galaxies
computed in three different ways. The top, middle, and bottom panels show the
gas mass-weighted metallicity, the intrinsic [O III] $\lambda$5007-, [O II]
$\lambda\lambda$3727-, and H$\beta$-weighted metallicity, and the metallicity
measured using the direct method along each viewing angle, respectively. The
different colours represent galaxies at different redshifts and the solid
lines on all panels represent running median fits to the mass-weighted
metallicity shown in the top panel. In the bottom panel, we show seven high-
redshift galaxies with [O III] $\lambda$4363 detections from Nakajima et al.
(2023) as black circles where we have recomputed metallicities to be
consistent with the atomic data used in Sphinx. We also show the stacked data
with [O III] $\lambda$4363 detections from Matthee et al. (2023) based on JWST
NIRCam grism data.

####  3.1.10 ISM conditions

NIRSpec access to emission lines at high-redshift means that certain ISM
properties such as gas density, metallicity, ionization parameter,
temperature, etc. can be inferred. As we showed above, this is paramount for
constraining the mass-metallicity relation; although, the emission lines
provide a biased view of the metallicity of the galaxy. The Sphinx data
release also contains various other ISM properties including gas density. In
Figure 12 we show histograms of ISM density in each redshift bin weighted by
intrinsic [O II] $\lambda\lambda$3727 (top) or [C III] $\lambda\lambda$1908
(bottom) emission. We highlight two interesting features. First, [C III]
$\lambda\lambda$1908 emission is generally probing higher gas densities than
[O II] $\lambda\lambda$3727 in our galaxies. Second, the [O II]
$\lambda\lambda$3727 gas density has little redshift evolution while the [C
III] $\lambda\lambda$1908 seems to increase with increasing redshift. This is
due to the fact that the different ionization potentials of O+ and C++ mean
that the emission lines originate from different regions of the nebula. These
points are both crucial for interpreting the evolution of ISM properties with
multi-line tracers and hence the Sphinx data release can be a valuable tool in
this regard as standard photoionization models cannot capture this behaviour
(see discussion in Katz et al. 2022b).

Figure 12: Histograms of ISM gas densities weighted by intrinsic [O II]
$\lambda\lambda$3727 (top) or [C III] $\lambda\lambda$1908 (bottom) emission
as a function of redshift.

####  3.1.11 LyC escape fractions

Constraining the sources that reionized the Universe is one of the primary
science goals of JWST. Because of its spatial resolution and multiphase ISM,
Sphinx is one of the few full-box radiation hydrodynamics simulations that can
predict escape fractions444These predictions are valid on the scale of the
spatial resolution of the simulation (i.e. $\sim 10$ pc).. LyC escape
fractions from Sphinx have been discussed at-length in previous works (Rosdahl
et al., 2018, 2022; Choustikov et al., 2023) and hence we refrain from further
discussion here. Methods for computing LyC escape fractions can be found in
the aforementioned references. The data set contains both angle-averaged
$f_{\rm esc}$ (which includes all photons with $E>13.6\ {\rm eV}$) as well as
those along the fiducial ten sight lines (for photons specifically at 900 Å).
In Figure 13 we show histograms of the angle-averaged LyC escape fraction
(i.e. the value that is important for reionization) for each redshift. The
values range from nearly 100% to zero. We highlight the tendency for the
typical $f_{\rm esc}$ of star-forming galaxies to decrease with decreasing
redshift, as is discussed at length in Rosdahl et al. (2022).

Figure 13: Angle-averaged LyC escape fractions for Sphinx galaxies as a
function of redshift.

####  3.1.12 Ionizing photon production efficiency

Similarly important for reionization as $f_{\rm esc}$ is the ionizing photon
production efficiency, $\xi_{\rm ion}\equiv Q_{\rm LyC}/L_{\rm UV}$, where
$Q_{\rm LyC}$ is the production rate of LyC photons and $L_{\rm UV}$ is the
monochromatic luminosity at 1500 Å. Early JWST observations have attempted to
constrain this value for star-forming galaxies during the epoch of
reionization (e.g. Simmonds et al., 2023). While the ionizing photon
production efficiency of each star particle is given by its age and
metallicity from our chosen BPASS SED, the complex SFHs of our galaxies mean
that distributions of $\xi_{\rm ion}$ and their evolution are non-trivial.
Furthermore, comparing the observed $\xi_{\rm ion}$ with the intrinsic
$\xi_{\rm ion}$ is key for interpreting how well observations constrain this
quantity. We refrain from that exercise here and show in Figure 14 the
intrinsic distribution of $\xi_{\rm ion}$ for galaxies in the Sphinx data
release. The typical values increase with increasing redshift as we sample
lower metallicity galaxies further above the main-sequence at higher redshift.
Due to our adopted SED, our values do not strongly deviate from those
typically assumed in reionization models (e.g. Robertson et al., 2015).

Figure 14: Histograms of intrinsic $\xi_{\rm ion}$ for the different redshift
snapshots available in the Sphinx data release.

###  3.2 Photometry and Imaging

In the previous section we primarily focused on the intrinsic galaxy
properties that are distributed as part of SPDRv1 as well as how they relate
to certain observable quantities. In this section, we continue our
demonstration of the database by highlighting data that is relevant to compare
with imaging and aperture photometry.

####  3.2.1 Photometry

The Sphinx data can be used to make mock photometric catalogs in any filter
that covers the rest-frame UV and optical555All photometry released in SPDRv1
is computed from the spectra with sedpy (https://github.com/bd-j/sedpy).. We
demonstrate this in Figure 15 where the top and middle panels show the filter
magnitudes of Sphinx galaxies for all NIRCam wide filters plotted against each
other for intrinsic and dust attenuated emission, respectively. For
comparison, the bottom panel shows the same relations for JADES galaxies from
the publicly available catalog (Bunker et al., 2023a; Eisenstein et al., 2023;
Hainline et al., 2023; Rieke & the JADES Collaboration, 2023) selected by
photometric redshift (or spectroscopic when available) and binned by similar
redshifts as the Sphinx data. We only consider galaxies with a signal-to-noise
ratio $>3$ in each of the filter combinations.

The dust-attenuated Sphinx galaxies cover a very similar magnitude
distribution as what is observed in JADES, making it an ideal comparison
sample. Nevertheless, there are some differences. First, JADES contains some
brighter galaxies that are not present in the Sphinx data release due to the
finite volume of the simulation. Similarly, although we implement an SFR
threshold, the Sphinx data release contains galaxies below the magnitude limit
of JADES. The most notable difference is that the observed galaxies are bluer
than those in Sphinx. For example, comparing F115W versus F444W, many JADES
galaxies fall below the one-to-one relation indicating they are very blue.

By comparing the top and middle panels of Figure 15, we can gain insight into
the origin of this discrepancy. Without dust, many of the Sphinx galaxies fall
below the one-to-one relation and are thus more consistent with JADES.
However, without dust the scatter in the Sphinx data is likely too small
compared to observations. One way to solve this it to adopt either a dust
attenuation law that is less steep with wavelength compared to the SMC. For
example the Reddy et al. (2015) and Calzetti et al. (2000) curves exhibit such
behaviour. One could also adopt a dust-to-gas mass ratio that falls off as a
power-law with metallicity rather than linearly, which is perhaps more
supported by observations (e.g. Rémy-Ruyer et al., 2014). A more top-heavy IMF
can also make the galaxies appear bluer (e.g. Stanway & Eldridge, 2023), and
simultaneously solves various other issues related to the number of observed
bright galaxies (e.g. Yung et al., 2023) and ALMA line emission (e.g. Katz et
al., 2022b). Finally, a reduction in total stellar mass (e.g. stronger
feedback) would preferentially reduce the flux at longer wavelengths.

From the observational perspective, the galaxies with the highest sSFRs are
likely to be most easily detected so there might be a small bias towards
seeing bluer galaxies in observations, or some of the galaxies may be low-
redshift interlopers. At $z=9$, the Ly$\alpha$ break falls in the F115W
filter, so one might expect that the pink points in Figure 15, representing
the $z=9$ galaxies, would be slightly offset to the left, e.g. from the $z=8$
sample, if there is no strong spectral evolution. This is not true for the
JADES data where the $z=9$ photometrically selected sample overlaps with the
other, lower-redshift galaxies. Hence the observed sample must be increasingly
blue with redshift as reported in Topping et al. (2023).

Finally, one should note the importance of emission lines. At $z=9$, [O III]
$\lambda$5007 and H$\beta$ fall in the F444W filter. If the emission lines in
Sphinx are too strong, our galaxies may appear too red. However, simulations
generally struggle to produce the observed tail of high equivalent widths of
[O III] $\lambda$5007$+$H$\beta$ (e.g. Wilkins et al., 2023b) so we believe
this is unlikely to cause the discrepancy. Furthermore, Sphinx galaxies appear
too red compared to JADES galaxies in other filters besides F444W and F115W.
If many of the JADES galaxies are Ly$\alpha$ emitters, this could similarly
boost the F115W flux but not solve the issue of Sphinx galaxies being too red
in F150W, unless other UV emission lines are similarly as bright. However such
bright UV lines are simply not observed in either Sphinx or spectroscopically
confirmed JADES galaxies (Curtis-Lake et al., 2022). Further analysis that is
beyond the scope of this paper will be required to elucidate the origin of the
differences between Sphinx and JADES.

Figure 15: Distribution of NIRCam wide filter magnitudes for Sphinx galaxies
using their intrinsic magnitudes (Top), dust-attenuated magnitudes (Middle),
and for galaxies observed as part of the JADES GTO program (Bottom). JADES
galaxies are selected based on photometric (or spectroscopic where available)
redshift. The dashed black lines represent the one-to-one relation.

####  3.2.2 Colour-colour selection

Colour selections designed to capture the Lyman break represent one of the
primary techniques for photometrically selecting samples of high-redshift
galaxies (e.g. Steidel et al., 1996; Giavalisco, 2002). Understanding how well
certain colour selections work in terms of false positive and false negative
rates are key for accurate determinations of luminosity functions and for
follow-up spectroscopic observations. The Sphinx public data release contains
magnitudes for all JWST filters along each sight line and filter magnitudes
for any other telescope can be trivially computed from the full spectra
provided (see below).

In Figure 16 we show an example colour-colour diagram and plot F115W$-$F150W
versus F150W$-$F277W for each galaxy at each redshift. The black line
demarcates the colour-colour selection from Harikane et al. (2023) for $z=9$
galaxies. It is clear that the vast majority of Sphinx $z=9$ galaxies are
selected by this criteria; however, there is some contamination from $z=8$
systems and a handful of $z=9$ galaxies fall outside the selection boundary.
For comparison, we also show the locations of JADES galaxies in the plane and
highlight the ones with photometric or spectroscopic redshifts $>9$. The
realistic spectra in the Sphinx public data release can be used to optimize
these criteria in combination with observations of potential lower-redshift
contaminants.

Figure 16: Example colour-colour diagram of JWST filters F115W$-$F150W versus
F150W$-$F277W. Sphinx galaxies are shown as coloured data points. Small black
points represent galaxies from JADES and cyan points are JADES galaxies with
photometric or spectroscopic redshifts of $z>9$. For comparison, the region
marked in black shows the $z=9$ selection criteria from Harikane et al.
(2023). Note that all Sphinx galaxies at $z=10$ are above the y-range of the
plot.

####  3.2.3 Photometric redshifts

Beyond colour-colour diagrams, one can attempt to estimate the exact redshifts
of the galaxies from their photometry. As an example, in Figure 17, we show a
violin plot of the distribution of photometric redshifts estimated with EAZY
(Brammer, van Dokkum & Coppi, 2008) using all of the wide and medium band
filters in JWST for all galaxies in the Sphinx data release at $z\geq 7$666For
the fits, we have assumed a 5% error on the flux in each filter and use the
default FSPS templates.. For galaxies at $z=7$ and $z=9$, the vast majority of
galaxies are assigned the correct redshift, with only a few galaxies being
placed at low-redshift. At $z=10$, the photomoetric redshift estimation is
also very good with a few systems preferring a higher redshift solution which
is the upper limit of our prior. Interestingly, most of the galaxies in the
database at $z=8$ are assigned a low-redshift solution where the Lyman break
is confused with a Balmer break. This is likely due to the fact that the Lyman
break occurs towards the edge of the F115W filter so the bluest band has much
lower flux than the second bluest (F140M), thus mimicking a weak break (as can
be seen in Figure 16). Hence, by working with only the maximum likelihood
redshifts, there could be a significant population of $z=8$ galaxies that are
incorrectly assigned a low-redshift solution in real JWST data. Here we have
used every JWST filter; however, the vast majority of imaging surveys use far
fewer filters, which would make photometric redshift estimation even more
uncertain than in this experiment.

Figure 17: Violin plots of photometric redshifts estimated with EAZY versus
the true redshift of the galaxies in the Sphinx database at $z\geq 7$. The
data points represent the median redshift and the vertical line represents the
interquartile range (which is only visible for the $z=8$ distribution).

####  3.2.4 UV spectral slopes

Because the UV spectral slope, $\beta$, depends on both age and metallicity of
the underlying stellar population as well as the dust content in galaxies, it
is a sensitive tracer to the SFH of a galaxy. Particular emphasis is placed on
finding galaxies with extremely blue UV slopes, which would be indicative of
very metal-poor or even metal-free stellar populations (e.g. Topping et al.,
2022). In principle, $\beta$ can be measured for large samples of galaxies
from photometry, making it one of the primary constraints from early JWST
data.

Using both the photometry and spectra in SPDRv1, in Figure 18 we show $\beta$
measured from the photometry for each Sphinx galaxy along each sight line
versus the true continuum slope along that direction (i.e. the stellar $+$
nebular continuum attenuated by dust). In general, most of the points fall
close to the one-to-one relation (dashed black line) and the median absolute
difference is 0.12 with a bias such that photometrically measured $\beta$ is
bluer than that from the spectra. The origin of this discrepancy is the
emission lines which have been removed from the spectra but not from the
photometry. This bias can partially be removed with SED fitting codes that
model UV emission lines. There are a few points that scatter above the one-to-
one relation and this is likely due to the fact that the transmission curves
are not perfectly uniform. This bias can once again be removed with SED
fitting and an accurate photometric redshift. Nevertheless, naive measurements
of $\beta$ that do not correct for UV emission lines can exhibit a typical
bias of 0.12. Our measurement was made for galaxies using all of the wide and
medium band filters on NIRCam and results may be different for different
subsets of filters. Furthermore, we have assumed an IGM that fully attenuates
Ly$\alpha$ and the bias may be stronger for Ly$\alpha$-emitting galaxies that
live in ionized bubbles.

Figure 18: Comparison between the UV slope measured photometrically and that
from the spectra for Sphinx galaxies at various redshifts. The dashed black
line represents the one-to-one relation. In general, UV slopes measured from
photometry tend to be steeper than slopes measured from the spectra.

Spectral slope is often compared to the UV magnitude of a galaxy. In Figure 19
we compare Sphinx galaxies to high-redshift photometrically-selected JWST
$\beta$ measurements from Cullen et al. (2023); Topping et al. (2023). The
Cullen et al. (2023) measurements represent individual galaxies that exhibit
significant scatter, often to unphysically blue values of $\beta$. In the
magnitude regime where the measurements overlap, the scatter is consistent
with what we see in Sphinx. The points from Topping et al. (2023) represent
the mean relation for samples of galaxies photometrically selected at
different redshifts. At the faint-end, there is little observed evolution in
$\beta$ between $z=6$ and $z=9$, consistent with Sphinx. The typical JADES
galaxy in Topping et al. (2023) is slightly bluer than what we find in Sphinx,
consistent with our earlier comparison to JADES galaxies. Among the many
effects that can cause this are the choice of dust attenuation curve in Sphinx
and the flux limit constraints in JADES (i.e. redder galaxies at the faint-end
are less likely to meet signal-to-noise thresholds).

At UV magnitudes fainter than $-17$ there is an upturn in the $\beta-{\rm
M_{UV}}$ relation in Sphinx. This is due to the SFR threshold used to create
SPDRv1 and would go away if we include the more numerous and fainter galaxies
with ${\rm SFR}<0.3\ {\rm M_{\odot}\ yr^{-1}}$.

Figure 19: $\beta$ measured photometrically versus ${\rm M_{UV}}$ for Sphinx
galaxies compared to individual high-redshift JWST measurements from (Cullen
et al., 2023) as well as mean relations of high-redshift JADES galaxies from
(Topping et al., 2023). The colour of the points indicates the redshift.

####  3.2.5 UV luminosity function

One of the earliest results from JWST was constraints on the high-redshift UV
luminosity function Bouwens et al. (2023); Harikane et al. (2023); Donnan et
al. (2023). Results for the Sphinx UV luminosity function (dust-attenuated
from an angle-averaged perspective) were compared to HST data in Rosdahl et
al. (2022). The Sphinx data release now contains the UV magnitudes along ten
different lines of sight. In Figure 20 we show that the luminosity functions
in the Sphinx data release (averaged over the ten directions) are in good
agreement with both photometric constraints from HST and JWST from $z\approx
5$ to 10. The overlap between observations and data at $z=9$ and 10 only spans
$\sim 2$ magnitudes; however, this will improve as observations go deeper and
gravitational lensing probes fainter galaxies. We note that the turnover in
our luminosity function at ${\rm M_{UV,1500\AA}}\approx-17$ is due to our SFR
threshold and not due to dwarf galaxy quenching from reionization or any other
process. It is clear that dust attenuation remains a necessary ingredient for
agreement with observations and the UV luminosity function is strongly
sensitive to the chosen dust model. However, we stress that our assumed dust
model was in no way designed to guarantee such an agreement and perhaps even
makes faint galaxies too red (see above discussion).

Figure 20: UV luminosity function for Sphinx galaxies at different redshifts
compared to high-redshift photometric constraints. The curves show the mean UV
luminosity functions over the ten viewing angles and the shaded areas
correspond to the standard deviation. Squares are observational data from HST
at $z\leq 8$ Bouwens et al. (2021), while circles, triangles, diamonds, stars,
cross, and hexagons show JWST measurements at $z\geq 8$ from Donnan et al.
(2023), Harikane et al. (2023) and Bouwens et al. (2023), Leung et al. (2023),
Franco et al. (2023), respectively. The cutoff in Sphinx data at high bright
magnitudes is due to the finite volume of our simulation box while the
downturn at faint magnitudes is a result of our SFR threshold.

####  3.2.6 ${\rm M_{*}}$ \- $\rm M_{UV}$

In Section 3.1.6, we connected UV magnitudes ($\rm M_{UV}$) with SFR. Since
SFR and stellar mass are correlated via the star formation main sequence, $\rm
M_{UV}$ also correlates with stellar mass. In the top panel of Figure 21 we
show stellar mass as a function of intrinsic and attenuated UV magnitudes at
1500 Å for galaxies at different redshifts in SPDRv1. We show intrinsic
magnitudes as dots and in horizontal lines the range of attenuated magnitudes
along the ten different sight lines. We also show fits to observational
inferences from Song et al. (2016) and Kikuchihara et al. (2020). The Sphinx
galaxies compare well against these observations only if dust attenuation is
ignored. At face value, this can be interpreted as Sphinx galaxies being too
massive and/or having too strong dust attenuation for a fixed UV magnitude,
both plausible outcomes if star formation is happening too rapidly very early
in the evolution of our simulated galaxies. This could also be caused by an
attenuation curve that is too strong at rest-frame UV wavelengths or any of
the other possibilities discussed in Section 3.2.1. However, we emphasize that
there remains significant uncertainty on the accuracy of stellar mass
estimated from photometry.

Figure 21: (Top) Galaxy stellar mass as a function of UV magnitude at 1500 Å
at various redshifts. The points represent intrinsic magnitudes for Sphinx
galaxies, while horizontal lines represent the range of attenuated magnitudes
along the 10 directions. For comparison we show observational inferences from
HST $+$ Spitzer data (Song et al., 2016; Kikuchihara et al., 2020). (Bottom)
Galaxy stellar mass as a function of F444W flux in $\mu$Jy. Over-plotted are
stellar masses inferred from JWST by CEERS (magenta, $z\sim 7-8$, Endsley et
al., 2023), GLASS (black, $z\sim 5-7$, Dressler et al., 2023), PEARLS (cyan,
$z\sim 7-12$, Trussler et al., 2023), COSMOS-Web (yellow, $z\sim 9-11$, Franco
et al., 2023), and NGDEEP (green, $z\sim 8-16$ Austin et al., 2023).

####  3.2.7 ${\rm M_{*}}$ \- ${\rm F444W}$

Rather than compare stellar masses with UV magnitude, one can alternatively
compare stellar mass with F444W flux. The advantages of doing so are 1) F444W
flux is much less sensitive to dust than UV magnitude so using this filter
avoids some of the uncertainty in dust modelling for high-redshift galaxies,
2) it samples a longer wavelength part of the spectrum which is more sensitive
to radiation from older stellar populations than UV magnitude and should thus
better probe stellar mass, and 3) it is a directly observable quantity and
therefore does not require modelling. Hence the scatter between stellar mass
and observed F444W flux at fixed redshift should be lower than between stellar
mass and UV magnitude.

Figure 22: [O III] $\lambda$5007 luminosity as a function of F444W luminosity
for Sphinx galaxies. We show both the intrinsic values (top) and dust
attenuated (bottom). Black points represent high-redshift galaxies from JADES
where the size of the point is indicative of redshift, with larger points
being at higher redshift.

In the bottom panel Figure 21 we plot stellar mass as a function of flux in
the NIRCam F444W filter for Sphinx galaxies as well as a large collection of
observational inferences of stellar mass as a function of F444W from multiple
JWST programs: CEERS (magenta, $z\sim 7-8$, Endsley et al., 2023), GLASS
(black, $z\sim 5-7$, Dressler et al., 2023), PEARLS (cyan, $z\sim 7-12$,
Trussler et al., 2023), COSMOS-Web (yellow, $z\sim 9-11$, Franco et al.,
2023), and NGDEEP (green, $z\sim 8-16$, Austin et al., 2023). These programs
all sample similar redshifts as the Sphinx data. However, despite Sphinx
galaxies showing a rather tight trend, there is significant scatter in the
data reported from the literature.

There are numerous physical effects that can drive scatter in this relation.
1) As redshift increases, the F444W filter samples a different part of the
galaxy SED, 2) As redshift increases, so does the luminosity distance, which
will decrease the observed flux, and 3) the [O III] $\lambda$5007$+$H$\beta$
emission lines enter the F444W filter at redshifts slightly below $z=7$.
Despite these effects, Sphinx galaxies still exhibit a rather tight trend;
although, there is clear evidence of the redshift effect. Nevertheless, these
effects cannot completely explain what is reported in the literature. Most
notably, the Dressler et al. (2023) sample is at lower redshift than the
Endsley et al. (2023) sample and predicts much higher stellar masses, yet both
exhibit similar ranges of F444W fluxes. Hence the Dressler et al. (2023)
galaxies should have lower stellar masses than those in Endsley et al. (2023),
unless [O III] $\lambda$5007$+$H$\beta$ overwhelmingly dominate the filter in
the Endsley et al. (2023) sample. This is unlikely given how wide the filter
is and the discrepancy in stellar mass is up to three orders of magnitude. In
Figure 22 we show the F444W flux against [O III] $\lambda$5007 luminosity for
Sphinx galaxies, both intrinsic (top) and dust attenuated (bottom). The JADES
galaxies have a typical [O III] $\lambda$5007 luminosity of $10^{42}\ {\rm
erg\ s^{-1}}$ and can thus increase the flux in the band by around a factor of
two. This is clearly not enough to explain the stellar mass discrepancy.
However, consistent with our above JADES discussion, this comparison shows
that the dust model in Sphinx is likely too strong at all wavelengths.
Removing dust would shift the Sphinx galaxies to the right on the stellar
mass-F444W relation.

In contrast to the Dressler et al. (2023) and Endsley et al. (2023)
inferences, the results from Trussler et al. (2023); Austin et al. (2023);
Franco et al. (2023) all overlap with Sphinx data. We can speculate the reason
for such different stellar mass estimates is somehow related to the way the
SEDs are being fit, the choice of emission line library and the prior on SFH.
However, based on these results, we emphasize that before they can be used on
constraints for galaxy formation physics, stellar mass estimates and therefore
any relation involving stellar mass (e.g the stellar mass-UV magnitude
relation, stellar mass-halo mass relation) should be made consistent amongst
the various observational campaigns.

####  3.2.8 Galaxy sizes

Galaxy sizes are computed in a similar manner to what is typically done for
real observations. For each mock observation in SPDRv1 we sum the stellar
continuum, nebular continuum, and nebular emission line images for each JWST
filter and use the PHOTUTILS (Bradley et al., 2023) package to segment the
image. We select the segment corresponding to the central galaxy in each halo
and measure the flux of the segment777Note that this flux of the main segment
can differ significantly from the flux of all stars within the virial radius
of the halo. This is particularly true in the case of a merger. as well as the
circularized half-light radius. By doing this, in principle, we can directly
compare galaxy sizes with observations.

In Figure 23 we show the circularized half-light radius in arcseconds as a
function of the segment flux in nJy in the F150W filter. For comparison, we
show recent data from JADES (Bunker et al., 2023a; Eisenstein et al., 2023;
Hainline et al., 2023; Rieke & the JADES Collaboration, 2023) for galaxies
that have a photometric redshift of $z_{\rm phot}>4$. There is significant
overlap between the observational and simulated data sets in terms of segment
fluxes; however, Sphinx galaxies scatter to significantly smaller radii than
what is observed with JWST. This is because we have not convolved our images
with a PSF. The horizontal dotted red lines show the pixel sizes for the short
and long-wavelength modes of NIRCam and thus all Sphinx galaxies with a radius
smaller than this line would likely be observed as a point source with a size
given by the PSF for the particular filter. The sizes of resolved sources in
our data set are very consistent with what is observed in JADES.

Figure 23: Circularized half-light radii in arcseconds for Sphinx galaxies in
the NIRCam F200W filter as a function of their segment flux in nJy. For
comparison, we show JADES galaxies with $z_{\rm phot}>4$ and an ${\rm S/N>3}$
in the F200W filter as black contours. The dotted horizontal red lines show
the pixel sizes of the short and long wavelength channels of NIRCam. Galaxies
smaller than this would appear as point sources convolved with the filter PSF.

For individual galaxies, we find a significant amount of dispersion in their
measured sizes depending on viewing angle. One can interpret that a
significant fraction of the observed dispersion in observations is explained
by viewing geometry, without needing to invoke intrinsic dispersion in SFR,
star formation history or extinction at a given magnitude, although variations
in all these parameters should be expected. In certain cases, the galaxy size
can vary by more than an order of magnitude. Part of this is due to the fact
that the galaxies are very clumpy and frequently undergo mergers. If the
merger is aligned along the line of sight, then the galaxy will have a very
small effective size while if the merger is perpendicular to the sight line,
the flux radius becomes very extended. The problem becomes exacerbated after
convolution with a PSF because these galaxies are often small enough where the
merger becomes blended.

####  3.2.9 Size - luminosity relation

While in the previous section we compared galaxy sizes with their observed
flux density, galaxy size is more commonly compared to the intrinsic galaxy UV
luminosity (e.g. Kawamata et al., 2018; Bouwens et al., 2022). In Figure 24 we
show galaxy size in the JWST F277W filter against the observed UV magnitude
for Sphinx galaxies at $z\geq 7$ compared with JWST observations in the same
filter and redshift interval from Yang et al. (2022); Franco et al. (2023). In
general, fainter galaxies have smaller sizes, but we predict a significant
$>1$ dex scatter at any magnitude. There are simply not enough observed
galaxies that overlap in magnitude with SPDRv1 to make any assessment on the
agreement between observations and our simulations888A further analysis of
JADES public data would likely have much more overlap due to the deeper
observations.. Furthermore, we highlight that despite our attempt to measure
sizes in a similar manner to observers, differences remain, beyond the fact
that we have not convolved with the PSF and added noise. For example, it is
typical to fit observations with a Sérsic profile with a slope of $n=1$.
Empirically, we find that this does not represent the clumpy nature of our
galaxies. Thus the reported definition of radius can vary between different
observational studies and what we have measured for Sphinx galaxies.
Nevertheless, because we provide the raw images with the data release, it is
possible to create a one-to-one mapping by applying the same methods to both
observations and the Sphinx galaxies.

Figure 24: Size-luminosity relation for Sphinx galaxies at $z\geq 7$ in the
JWST F277W filter. For comparison, we show recent JWST observations in the
same filter from Yang et al. (2022); Franco et al. (2023) as black and red
points, respectively.

###  3.3 Spectra

One of the key advantages of JWST in comparison to previous instruments is its
spectroscopic capabilities at high-redshift. In this section, we both describe
the spectral data products and demonstrate a few example use-cases of the
Sphinx public data release.

The spectra are stored at high resolution (1 Å) and can be degraded to match
any instrument with the appropriate line-spread function.

####  3.3.1 E(B$-$V)

As galaxies chemically evolve so does their dust content which can have a
significant impact on the observed spectrum. We show the impact of dust in the
form if interstellar reddening, E(B$-$V), as measured from the Balmer
decrement (H$\alpha$/H$\beta$) in Figure 25. As redshift decreases, the
typical E(B$-$V) of star-forming galaxies increases from $<0.1$ at $z\geq 9$
to nearly 0.2 at $z=4.64$. Since the database contains galaxies with both high
and low reddening, it can be used, for example, to test our ability to dust-
correct spectra and to predict completeness rates for luminosity functions.

Figure 25: E(B$-$V) distributions for spectra in the Sphinx database as a
function of redshift.

####  3.3.2 UV spectral slopes

We return to UV spectral slopes to demonstrate how the different components of
the spectrum (stellar continuum, nebular continuum, and dust) impact the
observed $\beta$ of a galaxy. We show in Figure 26 the spectral slope versus
galaxy age. The slope is computed between 1400 Å and 2500 Å (rest-frame) and
the galaxy age is the mass-weighted mean age of all star particles in the
galaxy. In the top panel, $\beta$ is computed from only the intrinsic stellar
continuum. Here we find a clear correlation between age and $\beta$, with
$\beta\sim-3$ for the youngest galaxies and $\sim-2.5$ for the oldest ones.
This highlights the expected intrinsic distribution of UV slopes given the
complex SFHs in Sphinx. The middle panel shows the same, but now with added
contributions from the nebular continuum. This modifies the correlation
between $\beta$ and age, with the nebular continuum significantly reddening
the youngest galaxies. When the nebular continuum is added, there is indeed
very little correlation between the mean mass-weighted age of the galaxy and
$\beta$. This is a combination of two effects. First, the nebular continuum
preferentially reddens the youngest and bluest galaxies. Second, even for
galaxies with relatively older ages, the UV part of the spectrum is still
often dominated by the youngest stars which can easily out-shine the older
population. This conspires to produce a very flat trend of $\beta$ with age.

However, when dust is included, the intrinsic trend between $\beta$ and age
can be further altered. In the bottom panel, we compute $\beta$ from the
stellar plus nebular continuua, dust attenuated along each of the ten sight
lines. In this case, the slope is widely scattered to very large $\beta$
values, even positive ones, and no correlation between $\beta$ and galaxy age
remains. It is therefore impossible to determine anything about the ages of
stellar populations in galaxies from observed UV slopes, unless a correction
for dust attenuation and nebular continuum can be accurately made. Many modern
SED fitting codes account for these effects simultaneously and we anticipate
that the Sphinx data can be used to test the accuracy of such methods.

Figure 26: UV slope, $\beta$, as a function of mass-weighted galaxy age,
measured from the intrinsic stellar continuum (top), intrinsic stellar +
nebular continuum (middle), and dust attenuated stellar + nebular continuum
along each of the ten sight lines (bottom) for each Sphinx galaxy. Small
circles represent individual galaxies and the large connected circles
represent mean values per age bin. Colour indicates the redshift of the
galaxy.

####  3.3.3 Emission line luminosity functions

Because Sphinx is a full-box simulation with initial conditions selected to
produce a typical halo mass function at $z=6$, it can be used to predict
emission line luminosity functions in the same way that we predicted the UV
luminosity function. In Figure 27 we show reddening corrected H$\alpha$ (top)
and [O III] $\lambda$5007 (bottom) luminosity functions at multiple redshifts
compared to JWST observational constraints from Matthee et al. (2023); Sun et
al. (2022). For each redshift, we show ten luminosity functions representing
each of the ten mock sight lines. This demonstrates the sample variance in the
Sphinx data set. As redshift decreases, the number density of line emitters
increases. We note that the turnover at the faint-end of the luminosity
function is due to our SFR cut rather than a physical effect, for example,
suppression due to the radiative feedback from reionization.

The Sphinx [O III] $\lambda$5007 luminosity function is in very good agreement
with results from the Eiger survey (Kashino et al., 2023; Matthee et al.,
2023) but under-predicts estimates from the four serendipitous line emitters
detected in JWST commissioning data (Sun et al., 2022). This likely explains
why the Sphinx H$\alpha$ luminosity function under-predicts the same
commissioning data. We note that the luminosity function from the
commissioning data is based on only four galaxies.

Figure 27: Line luminosity functions for H$\alpha$ (top) and [O III]
$\lambda$5007 (bottom) as a function of redshift. For each redshift, we show a
luminosity function for each of the ten different sight lines to demonstrate
the expected variance with viewing angle. For comparison we show results from
$z\sim 6$ JWST observations from Matthee et al. (2023); Sun et al. (2022).

####  3.3.4 Halo mass - Line luminosity relations

Understanding how halo mass connects with emission line luminosities is not
only important for understanding the galaxy-halo connection at high redshift,
but it is a key input for interpreting intensity mapping surveys. For example,
SPHEREx (Doré et al., 2014), CONCERTO (CONCERTO Collaboration et al., 2020),
and TIME (Crites et al., 2014) can conduct high-redshift intensity mapping for
emission lines in the rest-frame UV, optical, and IR up to high redshift.

In Figure 28 we show the relationship between observed H$\alpha$, [O III]
$\lambda$5007, and [C II] 158$\mu$m emission and halo virial mass. While there
is significant scatter in the relationship between nebular lines and halo
mass, there is a much tighter correlation between [C II] 158$\mu$m and halo
mass as this line tends to trace neutral gas. The Sphinx public data release
can be used to build relations between any of the available emission lines and
intrinsic galaxy properties which can be used to populate larger simulations
and make predictions for upcoming intensity mapping experiments.

Figure 28: Observed H$\alpha$ (top), [O III] $\lambda$5007 (middle), and [C
II] 158$\mu$m (bottom) emission versus halo virial mass. Because [C II]
158$\mu$m is not impacted by dust, we show only the intrinsic value of the
emission line as it is the same across all ten sight lines.

####  3.3.5 Equivalent width distributions

Because the Sphinx data set contains both emission line luminosities and the
stellar (and nebular) continuum, it is trivial to compute equivalent widths
(EWs) for any emission line in our sample. In Fig. 29 we show a histogram of
[O III] $\lambda$4959,5007$+$H$\beta$ EWs for each redshift bin in the data
set. As redshift increases, the typical EW in our data set also increases due
to the fact that we are sampling galaxies further above the main-sequence. One
important feature is that the data release contains nearly 400 sight lines
with ${\rm EW}>1000\ \text{\AA}$, which are typically rare in simulations. The
Sphinx data release is among the first to produce galaxies with such extreme
emission lines.

Figure 29: Equivalent width distributions of [O III] $\lambda$5007+H$\beta$
for Sphinx galaxies at each redshift.

####  3.3.6 Diagnostic diagrams

Diagnostic diagrams, for example the BPT and VO diagrams (Baldwin, Phillips &
Terlevich, 1981; Veilleux & Osterbrock, 1987), are important tools for
understanding the excitation mechanisms and conditions of the ISM throughout
the Universe. Understanding how these diagrams are expected to evolve with
changing ISM conditions as a function of redshift is key for interpreting JWST
spectra. For example, it is now well established that high-redshift galaxies
are offset on the BPT diagram compared to their low-redshift counterparts
(e.g. Steidel et al., 2014) and using the BPT to identify AGN at high-redshift
may fail (e.g. Groves, Heckman & Kauffmann, 2006; Feltre, Charlot & Gutkin,
2016).

In Figure 30 we show an example BPT diagram of Sphinx galaxies using the dust-
attenuated and de-reddened emission line luminosities along the ten sight
lines. The data set nicely maps out the region populated by the lowest
metallicity galaxies known in the low-redshift Universe (black triangles,
Isobe et al., 2022; Guseva et al., 2017) as well as JWST observations of $z>5$
star-forming galaxies from CEERS and JADES (pink circles, Sanders et al.,
2023; Cameron et al., 2023).

The upper cutoff in our diagram is driven by the fact that we assume a fixed
N/O ratio as a function of metallicity. This is a limitation of the Sphinx
data as elements are not enriched individually and we thus do not capture the
expected primary to secondary sequence as a function of metallicity (e.g.
Pilyugin et al., 2012). Decreasing N/O at lower metallicity would cause some
galaxies to scatter left and impact our upper sequence. Similarly if C/O also
decreases at lower O/H, the number of cooling channels is reduced which would
increase [O III]/H$\beta$. We also note a significant number of galaxies fall
above the classic Kewley et al. (2001); Kauffmann et al. (2003) demarcations
between star-forming galaxies and AGN. This is consistent with high-redshift
JWST observations (Sanders et al., 2023; Cameron et al., 2023) and is a result
of Sphinx galaxies achieving very high ionization parameters in extremely
dense gas999Part of this is a modelling choice because in gas cells with
unresolved Stromgren spheres, we act as if all constituent star particles are
co-located rather than assuming they are uniformly scattered throughout the
gas cell. In the latter case, the ionization parameter at the inner radius of
the shell is significantly decreased.. Since galaxies in Sphinx were selected
to be star-forming, we do not sample the downturn in the BPT seen for low-
redshift, much more passive galaxies with low-ionization parameters and high-
metallicities.

Figure 30: (Top) N II BPT diagram for Sphinx galaxies compared to SDSS (grey
2D histogram, Thomas et al., 2013), local metal-poor galaxies (black
triangles, Isobe et al., 2022; Guseva et al., 2017), and high-redshift JWST
galaxies (pink circles, Sanders et al., 2023; Cameron et al., 2023). (Bottom)
R23-O32 diagram for Sphinx galaxies compared to SDSS and high-redshift JWST
galaxies.

SPDRv1 is not limited to only the BPT diagram. In the bottom panel of Figure
30 we show an example R23-O32 diagram for Sphinx galaxies compared to SDSS and
$z>5$ JWST galaxies. Similar to the BPT, we find that Sphinx galaxies populate
the higher-excitation and lower-metallicity regions of the diagram, coincident
with high-redshift JWST galaxies. Because this diagram only depends on oxygen
emission lines, it suffers fewer systematics in modelling compared to the BPT.

####  3.3.7 Ly$\alpha$ luminosity function

Similar to the nebular line luminosity functions, direct comparisons can be
made between Sphinx and observations for Ly$\alpha$. For all redshifts we
provide data to make intrinsic Ly$\alpha$ luminosity functions and for
$z=4.64,\ 5,\&\ 6$ we provide Ly$\alpha$ luminosities for ten sight lines that
have been post-processed with Monte Carlo radiation transfer out to the virial
radius of the halo. We emphasize that IGM absorption can be significant at
these redshifts (e.g. Inoue et al., 2014) and our modelling choice is such
that additional IGM attenuation can be applied in post-processing to
understand how the IGM neutral fraction impacts the observability of
Ly$\alpha$ after it escapes the ISM and CGM.

Previous Ly$\alpha$ luminosity functions for Sphinx were presented in Garel et
al. (2021) and Katz et al. (2022a). The differences between those works and
here are that the former only studied the angle-averaged escape rather than
the line-of-sight values that are part of this data release. Furthermore, this
work has updated the method for estimating intrinsic Ly$\alpha$ emission, most
notably the fix for unresolved Stromgren spheres. As this fix primarily occurs
in very dense regions where Ly$\alpha$ photons are readily absorbed, our
results are very similar to our earlier work.

Figure 31: Ly$\alpha$ luminosity functions for Sphinx galaxies at $z=4.64,5$
and 6. Dashed lines show the intrinsic luminosity function while solid lines
correspond to luminosity functions after radiation transfer (i.e. absorption +
scattering, up to the virial radius of each halo) for ten viewing angles.
Circles, diamonds, upward triangles, downward triangles, hexagons, and squares
show observational constraints from Drake et al. (2017); Herenz et al. (2019);
Sobral et al. (2018); Santos, Sobral & Matthee (2016); Cassata et al. (2011);
Konno et al. (2018), respectively. Note that IGM transmission is not included
in these figures.

In Figure 31 we show Ly$\alpha$ luminosity functions for Sphinx galaxies
compared to observational constraints from $z=4.64-6$. Very good agreement is
found between simulations and observations. The only discrepancy occurs at
$z\sim 6$ at $L_{\rm Ly\alpha}<10^{42}\ {\rm erg\ s^{-1}}$ where Sphinx
predicts many more faint emitters than are observed. This may simply be due to
flux limits, incompleteness corrections in the observations, or IGM
attenuation (which is not included in our analysis). Nevertheless, we find
very little evolution in the Ly$\alpha$ luminosity function of star-forming
Sphinx galaxies in the redshift interval $z=4.64-6$. Future data releases may
contain processed Ly$\alpha$ at higher redshifts as it has now been detected
in very UV faint galaxies (Saxena et al., 2023) at up to $z\sim 11$ (Bunker et
al., 2023b).

####  3.3.8 Ly$\alpha$ and H$\alpha$ spectra and radial profiles

In addition to Ly$\alpha$ luminosities and escape fractions, the Sphinx data
release contains the full spectrum and radial profiles of Ly$\alpha$ and
H$\alpha$. This allows for determining how the shape of the spectrum and its
offset from line-centre relates to the intrinsic properties of the galaxy.

In Figure 32 we show the Ly$\alpha$ and H$\alpha$ spectrum of an example $z=6$
galaxy. The top left panel demonstrates the diversity of spectral profiles one
can see for Ly$\alpha$ for the same galaxy (e.g. Blaizot et al., 2023). The
spectrum varies between complete line-centre emission and full absorption. The
H$\alpha$ spectrum exhibits less variance in terms of luminosity; however,
there are strong features in the spectrum that correspond to individual star-
forming regions in the galaxy. Spectra, especially those at high-redshift, are
rarely sampled as well as in the Sphinx data release and are subject to a
line-spread function. To demonstrate the impact of the line-spread function,
we convolved the Ly$\alpha$ and H$\alpha$ spectra using a 1D Gaussian filter
across five spectral pixels and the results are shown in the right panels of
Figure 32. Where previously the H$\alpha$ spectra exhibited a significant
number of features, the convolved spectra appear significantly more Gaussian
as one typically sees with a real spectrograph. Similarly, the features of the
Ly$\alpha$ spectra are less well resolved after the convolution. We note that
the spectra have been computed in the rest-frame of the simulation box rather
than that of the galaxy. Hence when viewing the galaxy along different sight
lines, the peculiar velocities can shift the emission line from line-centre.
This can be corrected (as an observer might do) by shifting the spectrum in
velocity-space so that H$\alpha$ falls on line-centre. Finally, we re-
emphasize that Ly$\alpha$ spectra have only been integrated to the virial
radius of the galaxy and not through the IGM. Hence further IGM attenuation
must be applied to understand how these galaxies may appear in regions with
different neutral fractions.

Figure 32: (Left) Ly$\alpha$ and H$\alpha$ spectra along ten sight lines are
shown for an example $z=6$ galaxy at the fiducial spectral resolution. (Right)
The same spectra after being convolved with a Gaussian filter representing an
arbitrary line-spread function.

In Figure 33 we show the Ly$\alpha$ and H$\alpha$ radial profiles of the same
galaxy. For this particular object, H$\alpha$ is much more centrally
concentrated than Ly$\alpha$ as is does not resonantly scatter off of H I and
D. For many sight lines, there is no Ly$\alpha$ emission at all in the central
regions consistent with the spectral profiles that show complete absorption.
The comparison between the two profiles could provide indications into the
state of the ISM and CGM at high redshift and the LyC escape fraction
(Choustikov et al. in prep.).

Figure 33: Ly$\alpha$ and H$\alpha$ radial surface brightness profiles for the
same galaxy as shown in Figure 32.

##  4 Discussion

###  4.1 Comparison with Other Simulations

Although the Sphinx data release is one the most complete data sets in terms
of breadth of observational comparisons that can be made, our work is by no
means the first simulated data set for which mock observations have been
computed. In this section, we contextualize the Sphinx data release with
respect to other simulations where mock observations were published by briefly
highlighting some of the similarities and differences.

####  4.1.1 Comparison with semi-analytic models

Semi-analytic models are a computationally efficient way of producing mock
observational catalogs for a large population of galaxies. Simplistic models
are often applied to dark matter halo merger-trees in order to predict the
physical properties of galaxies. In the context of JWST, numerous such models
have been published (e.g. Williams et al., 2018; Yung et al., 2019). The
benefit to such models is that they can be tuned to match a wide range of
galaxy properties and predict the full SEDs of galaxies as they evolve. Such
models have been readily applied to design JWST observations in early cycles
when little was known about the spectroscopic properties of the high-redshift
galaxy population. However, due to their simplicity, semi-analytic models lack
detailed spatial information of galaxy properties, they cannot self-
consistently predict the evolving ISM properties of galaxies which are crucial
for emission line predictions and dust attenuation, and finally, they often
tie star formation with mass accretion, which conflicts with the bursty star
formation histories often predicted in cosmological simulations at early
epochs (e.g. Ma et al., 2018).

Due to computational expense, Sphinx only contains a single set of subgrid
models and is thus not as ideal for sampling the wide parameter space of
physical recipes in the way that semi-analytic models provide. However, Sphinx
is a significant improvement over such models as it resolves the vast majority
of physics that is important for generating the emission that is observed in a
cosmological volume large enough to make meaningful comparisons with
observations.

####  4.1.2 Comparison with hydrodynamic simulations

Over the past decade, large suites of cosmological hydrodynamics simulations
(e.g. Vogelsberger et al., 2014; Dubois et al., 2014; Schaye et al., 2015)
have been run that can reproduce many of the observational properties and the
diversity of galaxies at low-redshift. Because the subgrid models are tuned to
reproduce the low-redshift galaxy population, if is often the case that the
models differ at high-redshift. Thus JWST observations are a crucial
ingredient for differentiating between galaxy formation models.

Vogelsberger et al. (2020) recently presented UV and H$\alpha$ luminosity
functions for multiple dust models for galaxies in the IllustrisTNG
simulations. In comparison to our work, due to the much higher spatial
resolution and mass resolution of Sphinx compared to IllustrisTNG, as well as
ISM modelling, we find burstier star formation histories. This results in an
increased scatter in the relation between UV magnitude and stellar mass. For
example, at $M_{\rm UV}=-18$, we find a scatter of $\sim 0.5$ dex in stellar
mass compared to Vogelsberger et al. (2020) where the value is typically
$\lesssim 0.2$ dex. Our value is in much better agreement with observations
(Song et al., 2016), even if the normalization of Sphinx is offset in this
relation. Because IllustrisTNG does not model the ISM, H$\alpha$ luminosities
are painted on to star particles in order to predict an H$\alpha$ luminosity
function. Hence their predictions are not properly modulated by the different
ISM conditions (e.g. gas temperature, density, metallicity contrast between
dust and stars, LyC escape fraction, diffuse ionized gas, etc., see Tacchella
et al. 2022) that impact H$\alpha$ emission. Similar to the UV luminosity
function, their model will predict significantly less scatter than ours for a
given star formation rate.

Hirschmann et al. (2022) have performed an alternative post-processing of
IllustrisTNG compared to Vogelsberger et al. (2020) in order to predict the
evolution of emission line luminosities and diagnostic diagrams. The upside to
these models compared to Sphinx is that they sample both star-forming galaxies
and AGN and they have an explicit model for capturing emission from shocks
when it is not always well resolved in Sphinx. However, the major limitation
still remains that IllustrisTNG does not model the ISM. For this reason,
Hirschmann et al. (2022) must assume a gas density (which they choose to be
$10^{2}~{}{\rm cm^{-3}}$) in H II regions; whereas observations show that this
value likely increases with redshift (Isobe et al., 2023). Despite Sphinx
having a much smaller volume compared to IllustrisTNG-50 (20 cMpc vs. 50
cMpc), we find that our emission line luminosities tend to be higher than
those in Hirschmann et al. (2022). For example, the brightest [O III]
$\lambda$5007 emitter at $z=6$ in IllustrisTNG-50 (based on their Figure 11)
has a similar intrinsic luminosity to the one in Sphinx. This once again may
be due to the burstier nature of star formation in Sphinx compared to
IllustrisTNG as well as less regulatory feedback in Sphinx. Other factors that
could impact this result are differences in ISM gas density (Sphinx has much
higher gas densities which can lead to higher intrinsic emission) or mass-
metallicity relation.

The post-processing technique also matters significantly in terms of mock
observations. Shen et al. (2020) also post-processed IllustrisTNG and find
significantly different emission line luminosities compared to Hirschmann et
al. (2022). Hence when the ISM is not resolved, the freedom in post-processing
parameters likely represents the dominant systematic uncertainty in the
ability to make observational predictions with simulations.

BlueTides (Wilkins et al., 2016, 2017; Marshall et al., 2022) is one of the
largest full-box cosmological simulations of high-redshift galaxy formation.
Due to its large volume, it probes galaxies orders of magnitude more massive
than what is in Sphinx and hence it is much better for understanding the
physics dictating the bright-end of the galaxy UV luminosity function. Like
IllustrisTNG, BlueTides does not have an ISM so even though the full SEDs of
galaxies can be predicted, the impact of nebular emission is calculated by
assigning photoionization models with a density of $10^{2}~{}{\rm cm^{-3}}$
and an escape fraction of zero to star particles. This ignores any potential
difference between stellar and nebular metallicity (as we show above can be
large in the case of a very dense ISM) and essentially fixes the ionization
parameter for a given stellar age and metallicity. Similar to Sphinx, Wilkins
et al. (2016) find that the nebular continuum can be an important component
for high-redshift galaxy SEDs; however, the mass ranges of the two simulations
do not overlap well enough to be compared.

The FLARES simulations (Lovell et al., 2021) are a suite of high-redshift zoom
simulations employing the EAGLE model (Crain et al., 2015) for galaxy
formation. FLARES publicly provides numerous galaxy properties as well as SED
fluxes, line luminosities, and EWs, similar to the Sphinx data release. While
not a full volume simulation, the suite is large enough to contain a diversity
of galaxies. For example, FLARES contains high-redshift galaxies with Balmer
breaks (Wilkins et al., 2023a) that are also present in Sphinx (see Steinhardt
et al. 2023) and have been observed at high redshift (Hashimoto et al., 2018),
although c.f. Bradač et al. (2023); Stiavelli et al. (2023). Similar to
BlueTides and IllustrisTNG, FLARES does not have an ISM so nebular emission is
assumed to originate from an ionization-bounded nebula with a maximum
ionization parameter of $-2$, assuming that gas near the star particles have
the same metallicity as the star. One of the notable differences between
Sphinx and FLARES is that many of our galaxies can reach [O III]
$\lambda$5007/H$\beta>10$ which are not seen in FLARES. This may be due to the
fact that in general Sphinx galaxies can exhibit ionization parameters much
greater than $-2$.

FirstLight (Ceverino, Klessen & Glover, 2019) adopts a similar approach to
FLARES by using multiple zoom simulations to capture the diversity of galaxies
at high-redshift. Unlike the previously discussed simulations in this section,
FirstLight attempts to self-consistently model the ISM. Despite this, nebular
regions are modelled using a constant gas density of $10^{2}~{}{\rm cm^{-3}}$
and intrinsic emission line luminosities are adopted from those provided with
BPASS (Xiao, Stanway & Eldridge, 2018).

FirstLight deals only with intrinsic emission rather than dust-attenuated
quantities. Galaxies in FirstLight exhibit bursty star formation histories,
similar to Sphinx and they similarly find that the nebular continuum sets a
lower limit on UV slope (although ours is $\sim-2.7$ whereas it is closer to
$-2.5$ in their work). One considerable difference is the lack of scatter in
the N II BPT diagram in FirstLight compared to Sphinx as well as local metal
poor galaxies (Guseva et al., 2017). Furthermore, Sphinx probes [O III]
emitters with observed, rest-frame equivalent widths $>1000\ \text{\AA}$ that
are not absent from the intrinsic estimates in FirstLight. This is partially
due to the fact that Sphinx probes galaxies with higher sSFR; but, also
possibly due to the high densities reached in the nebular regions in Sphinx.

####  4.1.3 Comparison with radiative-hydrodynamic simulations without a
multiphase ISM

Simulations like THESAN (Kannan et al., 2022a) improve upon their related
predecessor IllustrisTNG due to the addition of on-the-fly radiation
hydrodynamics. Public data is provided for Ly$\alpha$, SEDs, and various other
galaxy properties similar to FLARES and now Sphinx. The volume of THESAN is
significantly larger than Sphinx; hence, it probes much more massive galaxies
and a more diverse IGM.

Like the simulations in the previous section, THESAN does not have a
multiphase ISM which means that Ly$\alpha$ and nebular predictions are
entirely subject to the ISM assumptions made in post-processing. The same
holds true for the LyC escape fraction estimates from THESAN. As with
IllustrisTNG this can lead to order of magnitude systematic uncertainties for
emission predictions from the same simulation. This is particularly true for
the intensity-mapping predictions presented in Kannan et al. (2022b) regarding
[C II] 158$\mu$m emission, which is predicted to mostly come from neutral ISM
gas at high-redshift (e.g. Pallottini et al., 2017), not H II regions with a
fixed ionization parameter of $-2$. We have processed Sphinx with sub-ionizing
radiation specifically to better capture the physics of photodissociation
regions. For similar reasons, without an ISM, THESAN cannot reliably predict
the spectral shapes of Ly$\alpha$ escaping the ISM (Smith et al., 2022) which
are set by the complex kinematics and ionization state of the dense gas near H
II regions (e.g. Blaizot et al., 2023). The simulation is perhaps better
suited for measuring IGM transmission curves, similar to CODA (Park et al.,
2022) due to the large volume compared to Sphinx.

One final notable difference is that Sphinx exhibits a much higher conversion
efficiency of baryons into stars compared to THESAN. While the high efficiency
is less in agreement with pre-JWST abundance matching estimates (e.g. Behroozi
et al., 2019), and may still be too efficient given our discussion comparing
JADES and Sphinx photometry, these older estimates are inconsistent with the
surface densities of bright galaxies at high-redshift (e.g. Finkelstein et
al., 2023; Leung et al., 2023). This is shown in Figure 34 where we compare
the cumulative surface density of sources with F277W magnitudes $<29.5$ from
Sphinx to recent results from NGDEEP (e.g. Leung et al., 2023) and other
simulations101010All curves for this plot besides those from Sphinx were
digitized (Rohatgi, 2022) from Figure 4 in Leung et al. (2023).. Indeed THESAN
and Universe Machine (Behroozi et al., 2019) under predict the results from
NGDEEP while Sphinx is in agreement with the data up to $z\sim 9$, where it
then begins to over predict the galaxy counts.

Figure 34: Cumulative surface density of sources with F277W magnitudes $<29.5$
at redshifts $>z$ for Sphinx galaxies (red) compared to JWST NGDEEP
observations (black, Leung et al., 2023) and other simulations/models. The
thick red line shows the mean across all sight lines in our simulations while
the thin red lines represent individual viewing angles. The grey shaded region
represents the uncertainty on the observations.

####  4.1.4 Comparison with zoom-in radiative-hydrodynamic simulations

The Renaissance simulation (O’Shea et al., 2015) is a suite of high resolution
cosmological radiation hydrodynamics zoom simulations that samples large
regions from under-dense, to over-dense, stopping at progressively higher
redshift. The spatial resolution is higher than Sphinx and the stellar masses
in Renaissance are in good agreement with inferences from high-redshift JWST
observations (McCaffrey et al., 2023). SEDs and JWST photometry for
Renaissance galaxies were presented in Barrow et al. (2017) and due to the
simulation volume, Renaissance samples much lower mass galaxies, and only at
higher redshifts. The post-processing method shares many similarities with the
methods presented in this work. Consistent with Sphinx, they find a large
scatter in the [N II] BPT diagram at low-metallicity and that scattering can
play an important role in modulating the observed spectrum. Similar analysis
techniques were applied to cosmological radiative hydrodynamics simulations of
Barrow et al. (2020) who find similar results to Sphinx (specifically
Choustikov et al. 2023) with regards to how LyC escape trends with nebular
emission.

Numerous other zoom-in cosmological radiation hydrodynamics simulation of
individual high-redshift objects have been presented in the literature (e.g.
Katz et al., 2017, 2019; Pallottini et al., 2017, 2019; Lupi et al., 2020;
Trebitsch et al., 2021). For those that provide emission lines and spectra,
post-processing techniques are similar to those used in this work — i.e. non-
equilibrium emission is computed when possible, otherwise photoionization
models are used. The simulations of Katz et al. (2019) are most similar to
Sphinx in terms of mass, spatial resolution, and feedback, although that work
included more radiation bins (primarily sampling lower energies), the dust
model followed Rémy-Ruyer et al. (2014), and molecular hydrogen was self-
consistently modelled. No significant differences can be reported in terms of
emission line properties, validating the choice of a metallicity floor and the
RT post-processing in the context of the emission lines studied here.

###  4.2 Future Updates

While the Sphinx data release represents the current state-of-the-art in
comparing high-redshift galaxy simulations with observations, we consider the
data set to be an evolving tool rather than a static reference. We describe
below some of the intended upgrades following the initial data release.

  1. 1.

Larger galaxy sample: The current data release contains only $\sim 1,400$
galaxies at seven different redshift bins. This is significantly fewer than
the nearly 30,000 resolved galaxies in each Sphinx snapshot. Data is stored
from the simulation on a 5 Myr cadence; hence, there is significant room to