The Sphinx Public Data Release: Forward Modelling High-Redshift JWST Observations with Cosmological Radiation Hydrodynamics Simulations

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 Sphinx$^{20}$ public data release. Sphinx$^{20}$ 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 $\sim14,000$ mock images and spectra of the stellar continuum, nebular continuum, and 52 nebular emission lines, including Ly$\alpha$, for each galaxy in Sphinx$^{20}$ with a star formation rate $\geq0.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 Sphinx$^{20}$ data release can be downloaded at the following URL: https://github.com/HarleyKatz/SPHINX-20-data


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 highredshift 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, Finkel- * E-mail: harley.katz@physics.ox.ac.uk stein & 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 (ξ 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 MAP-PINGS (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 lowredshift "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-α (Lyα) 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 ≳ 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.They represent a complementary means to lowredshift 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. 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α 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α (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 Sphinx 20 (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 (≈ 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 highdegree 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 opensourced library of ∼ 14, 000 mock images and spectra for galaxies in the redshift range 4.64 ≤ z ≤ 10 with 10 Myr-averaged star formation rates ≥ 0.3 M ⊙ 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.

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.

The Sphinx 20 Simulations
Sphinx is a suite of cosmological radiationhydrodynamic 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 drives the observability of Lyman-alpha emitters at z ≳ 6 (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).
Sphinx 20 , which is exclusively used in this data release, is the largest volume simulation in the Sphinx suite, with a volume of 20 3 co-moving Mpc 3 .It contains thousands of star-forming galaxies at each snapshot between z = 4.64 and z = 10.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.
Sphinx 20 was run over the period from February 2019 to November 2020 on the JUWELS 1 and Joliot-Curie ROME 2 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(Rosdahl et al. , 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×10 5 M ⊙ , allowing halos to be resolved approximately down to the atomic cooling limit, and stellar particles, representing stellar populations, have initial masses of 400 M ⊙ .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 M ⊙ , which is significantly more than found with standard initial mass functions observed in our local environment (but see Katz et al. 2022b).Sphinx does not include black hole formation or AGN feedback due to its small volume, although new data is suggesting that black holes may be more relevant than realised at early epochs (Greene et al. 2023;Matthee et al. 2023b;Kocevski et al. 2023).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.59eV and 24.59 eV-∞), but post-process the simulation snapshots used in this data release with two additional lower-energy radiation groups (5.6 eV-11.2eV 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 onthe-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 1 https://www.fz-juelich.de/en/ias/jsc/systems/supercomputers/juwels 2 https://www-hpc.cea.fr/en/Joliot-Curie.html 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 Sphinx 20 snapshots, spanning z = 10 to z = 4.64, and only galaxies which have star formation rates larger than 0.3 M ⊙ 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.

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 M ⊙ ) and IMF slope (−1.35) as was used to calculate the ionizing luminosities for star particles in the simulations.

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 S < ∆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,  (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.

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).

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α 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 HI + f ion n HII )σ dust Z/Z 0 , where n HI and n 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 ion = 0.01, σ 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α is computed slightly differently to the method described above.Lyα 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α 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 RAS-CAS runs as described above and we assume a fixed D/H abundance of 3 × 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).

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.

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 ≥ 0.3 M ⊙ 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.
To begin, in Figure 1 we show a gallery of mock RGB images of 100 Sphinx galaxies in the redshift interval 6 ≤ z ≤ 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.
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.

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 Sphinx 20 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 M ⊙ at z = 10, 8, 6, and 4.6, respectively.Our data release includes galaxies with SFR 10 ≥ 0.3 M ⊙ yr −1 that are mostly embedded in halos with M vir ∼ 10 9 M ⊙ at z = 10 or halos with M vir ∼ 10 10 M ⊙ at z = 4.6, but a small fraction of them is also found in halos with 10 8 < M vir < 10 9 M ⊙ , 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).

Stellar masses
In the right panel of Figure 3, we show the distributions of stellar masses formed (M * ,i ) in Sphinx 20 galaxies.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.
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-toone (primarily because of dust attenuation); however, we cannot determine the UV luminosity until after post-Fig.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.Intrinsic for all emission lines listed in Table 2, dust attenuated along ten sight lines for Hα, Hβ, Hγ, Hδ,  processing (which is the computationally limiting step).Our cut excludes a large fraction of galaxies with M * ,i ≲ 10 8 M ⊙ .For example, our galaxy samples are 90% complete in stellar mass down to log 10 (M * ,i /[M ⊙ ]) = 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.

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. (2018Rosdahl et al. ( , 2022) ) and is available for starforming 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.
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 re-solved stellar mass estimates at high-redshift often result in stellar masses that are ∼ 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 highredshift 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).

Star formation rates
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 Hα recombination line results from gas ionised by young stars, and hence gives a hint of the current SFR on ≲ 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 SFR 100 , SFR 5 , and SFR 3 ) as a function of 10 Myr-averaged SFR (SFR 10 ).All the galaxies in the data release are selected to have SFR 10 ≥ 0.3 M ⊙ yr −1 , but on different time scales, the SFR can drop well below this value.For example, on average, SFR 100 tends to fall well below SFR 10 , indicating that the 10 Myr-averaged SFR is more representative of a burst rather than the continuous star formation history.In contrast, SFR 5 and SFR 3 are more evenly scattered around the one-to-one relation.However, there are a few galaxies with very low SFR 5 and SFR 3 , which is indicative of a short-timescale quenching event.

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 Myraveraged (bottom) SFRs as a function of stellar mass.The SFRs in the Sphinx data release reach as high as 100 M ⊙ 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. Fig. 5.-Star formation rates averaged over 100 Myr (SFR 100 , top), 5 Myr (SFR 5 , middle) and 3 Myr (SFR 3 , bottom) as a function of 10 Myr-averaged SFR (SFR 10 ).Individual points represent Sphinx galaxies at different redshifts, selected to have SFR 10 ≥ 0.3 M ⊙ yr −1 .The black diagonal lines show the one-toone relations.
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 Myraveraged 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 sur-  vey where the fainter galaxies at higher-redshifts must fall consistently higher above the main-sequence to be observed.

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 SFR ≥ 0.3 M ⊙ 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.
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 SFR ≥ 0.3 M ⊙ 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 (ρ 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 Å, M UV,1500 , and Fig. 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 SFR 10 ≥ 0.3 M ⊙ 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 M UV (yellow) or dust obscured 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.
10 Myr-averaged SFR from our data set such that log 10 SFR 10 M ⊙ yr −1 = −0.35M UV,1500 − 6.49.(1) We then select all galaxies that have an intrinsic M UV ≤ −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 M UV 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 UV ≤ −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 ρ UV for spectroscopically confirmed high-redshift JWST galaxies with M UV ≤ −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) 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.

Stellar ages
The stellar age of a galaxy is an important parameter that encapsulates its mass growth history.Massweighted 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 M * ≳ 10 7.5 M ⊙ 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 M ⊙ yr −1 represents a stellar mass of 3 × 10 6 M ⊙ , which is a significant fraction of the total stellar mass of objects with M * < 10 7.5 M ⊙ 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 sen- sitive to the young stellar populations and hence probe the recent star formation in the past ∼ 10 Myr.

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 highredshift 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 starforming 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] λ5007, [O II] λλ3727, and Hβ 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 2009) conversion between t O3 and t O23 .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).

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] λλ3727 (top) or [C III] λλ1908 (bottom) emission.We highlight two interesting features.First, [C III] λλ1908 emission is generally probing higher gas densities than [O II] λλ3727 in our galaxies.Second, the [O II] λλ3727 gas density has little redshift evolution while the [C III] λλ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).
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 fractions 4 .LyC escape fractions from Sphinx have been discussed at-length in previous works (Rosdahl et al. 2018(Rosdahl et al. , 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 esc (which includes all photons with E > 13.6 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 esc of star-forming galaxies to decrease with decreasing redshift, as is discussed at length in Rosdahl et al. (2022).
3.1.12.Ionizing photon production efficiency 4 These predictions are valid on the scale of the spatial resolution of the simulation (i.e.∼ 10 pc).Similarly important for reionization as f esc is the ionizing photon production efficiency, ξ ion ≡ Q LyC /L UV , where Q LyC is the production rate of LyC photons and L 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 ξ ion and their evolution are non-trivial.Furthermore, comparing the observed ξ ion with the intrinsic ξ 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 ξ 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).

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.

Photometry
The Sphinx data can be used to make mock photometric catalogs in any filter that covers the rest-frame UV and optical5 .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 lowredshift interlopers.At z = 9, the Lyα 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] λ5007 and Hβ 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] λ5007+Hβ (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α 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.

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.

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 ≥ 7. The top panel shows the results using the default EAZY FSPS templates while the bottom panel uses the templates designed for high-redshift galaxies from Larson et al. ( 2022).For the fits, we have assumed a 5% error on the flux in each filter.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 photometric redshift estimation is also very good with a few systems preferring a higher redshift solution which is the upper limit of our prior.Interestingly, when the default templates are used, 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.This confusion does decrease significantly when adopting the templates from Larson et al. (2022).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.

UV spectral slopes
Prior to JWST, UV spectral slope, β, represented one of the few diagnostics available that could be directly measured from observations (e.g.Dunlop et al. 2012;Finkelstein et al. 2012;Bouwens et al. 2014;Wilkins et al. 2016a;Bhatawdekar & Conselice 2021).Because β depends on both age and metallicity of the underlying stellar population as well as the dust content in galaxies, it traces the recent 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, β 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 β 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 β is bluer than that from the spectra.There are two effects that can cause this behaviour.Emission lines are present in the observed photometry and can thus bias the beta measurement.To assess the importance of emission lines, we have generated the photometry for each galaxy without emission lines and fit the UV slope.In doing this, we find that the median absolute difference drops to 0.08 and thus emission lines are not the dominant source of error.The remaining and dominant source of error is thus due to the shape of the filter response curves.Both biases can partially be removed with SED fitting codes.Nevertheless, naive measurements of β that do not correct for these effects 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.We have also tested using only the wide filters and find similar results.Furthermore, we have assumed an IGM that fully attenuates Lyα and the bias may be stronger for Lyα-emitting galaxies that live in ionized bubbles.
The UV 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 β 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 β.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 β 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 β − 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 SFR < 0.3 M ⊙ yr −1 .

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 (dustattenuated 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 ≈ 5 to 10.The overlap between observations and data at z = 9 and 10 only spans ∼ 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 M UV,1500 Å ≈ −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).

M * -MUV
In Section 3.1.6,we connected UV magnitudes (M UV ) with SFR.Since SFR and stellar mass are correlated via the star formation main sequence, 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.

M * -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  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.should be lower than between stellar mass and UV magnitude.
In the bottom panel Figure 21  tance, which will decrease the observed flux, and 3) the [O III] λ5007+Hβ 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 Ends-   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] λ5007 luminosity for Sphinx galaxies, both intrinsic (top) and dust attenuated (bottom).The JADES galaxies have a typical [O III] λ5007 luminosity of 10 42 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  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.

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 segment6 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 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.
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.

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 ≥ 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 simulations7 .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.

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.

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α/Hβ) in Figure 25.As redshift decreases, the typical E(B−V) of star-forming galaxies increases from < 0.1 at z ≥ 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.

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 β 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, β is computed from only the intrinsic stellar continuum.Here we find a clear correlation between age and β, with β ∼ −3 for the youngest galaxies and ∼ −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 β 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 β.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 β with age.
However, when dust is included, the intrinsic trend between β and age can be further altered.In the bottom panel, we compute β 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 β values, even positive ones, and no correlation between β 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.

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 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] λ5007 luminosity function is in very good agreement with results from the Eiger survey (Kashino et al. 2023;Matthee et al. 2023a) but underpredicts estimates from the four serendipitous line emitters detected in JWST commissioning data (Sun et al. 2022).This likely explains why the Sphinx Hα luminosity function under-predicts the same commissioning data.We note that the luminosity function from the commissioning data is based on only four galaxies.

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), CON-CERTO (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α, [O III] λ5007, and [C II] 158µ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µ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.

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] λ4959,5007+Hβ 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 mainsequence.One important feature is that the data release contains nearly 400 sight lines with EW > 1000 Å, which are typically rare in simulations.The Sphinx data re- lease is among the first along with the FLARES simulation Wilkins et al. (2023b) to produce galaxies with such extreme emission lines.However, neither simulation produces the long tail of high EW galaxies observed at high-redshift (e.g.De Barros et al. 2019).The origin of this discrepancy is currently unknown.
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 highredshift 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 dereddened 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β.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 gas8 .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.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.

Lyα luminosity function
Similar to the nebular line luminosity functions, direct comparisons can be made between Sphinx and observations for Lyα.For all redshifts we provide data to make intrinsic Lyα luminosity functions and for z = 4.64, 5, & 6 we provide Lyα 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α after it escapes the ISM and CGM.
Previous Lyα 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  Konno et al. (2018), respectively.Note that IGM transmission is not included in these figures.
than the line-of-sight values that are part of this data release.Furthermore, this work has updated the method for estimating intrinsic Lyα emission, most notably the fix for unresolved Stromgren spheres.As this fix primarily occurs in very dense regions where Lyα photons are readily absorbed, our results are very similar to our earlier work.
In Figure 31 we show Lyα 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 ∼ 6 at L Lyα < 10 42 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 atten-uation (which is not included in our analysis).Nevertheless, we find very little evolution in the Lyα luminosity function of star-forming Sphinx galaxies in the redshift interval z = 4.64 − 6. Future data releases may contain processed Lyα at higher redshifts as it has now been detected in very UV faint galaxies (Saxena et al. 2023) up to z ∼ 11 (Bunker et al. 2023b). 3.3.8.Lyα and Hα spectra and radial profiles In addition to Lyα luminosities and escape fractions, the Sphinx data release contains the full spectrum and radial profiles of Lyα and Hα.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α and Hα spectrum of an example z = 6 galaxy.The top left panel demonstrates the diversity of spectral profiles one can see for Lyα for the same galaxy (e.g.Blaizot et al. 2023).The spectrum varies between complete line-centre emission and full absorption.The Hα 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 highredshift, 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α and Hα 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α 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α 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α falls on linecentre.Finally, we re-emphasize that Lyα 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.
In Figure 33 we show the Lyα and Hα radial profiles of the same galaxy.For this particular object, Hα is much more centrally concentrated than Lyα as is does not resonantly scatter off of H I and D. For many sight lines, there is no Lyα 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.).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.

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 selfconsistently 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.

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.In general, these hydrodynamic simulations have much larger volumes compared to Sphinx due to the computational demands of on-the-fly radiative transfer.Therefore, they capture the brighter end of the UV luminosity function in a regime where our current work cannot make predictions.Hence they represent a complementary view of the high-redshift Universe.Vogelsberger et al. (2020) recently presented UV and Hα 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 UV = −18, we find a scatter of ∼ 0.5 dex in stellar mass compared to Vogelsberger et al. (2020) where the value is typically ≲ 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α luminosities are painted on to star particles in order to predict an Hα 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α 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 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] λ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. 2016b(Wilkins et al. , 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 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. (2016b) 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] λ5007/Hβ > 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.However, FLARES contains a higher proportion of galaxies with [O III] λ4959,5007+Hβ > 1000 Å in better agreement with observations.This may be due to differences in dust modelling and star formation recipe.
FirstLight (Ceverino, Klessen & Glover 2019) adopts a similar approach to FLARES by using multiple zoom simulations to capture the diversity of galaxies at highredshift.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 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 ∼ −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 Å 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α, 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.Moreover, THESAN includes black holes which are absent from Sphinx and may be relevant at high redshift (e.g.Greene et al. 2023;Matthee et al. 2023b;Kocevski et al. 2023).
Like the simulations in the previous section, THESAN does not have a multiphase ISM which means that Lyα 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 THE-SAN.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µm emission, which is predicted to mostly come from neutral ISM gas at highredshift (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α 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 better suited than Sphinx 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 simulations 9 .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 ∼ 9, where it then 9 All curves for this plot besides those from Sphinx were digitized (Rohatgi 2022) from Figure 4 in Leung et al. (2023).
begins to over predict the galaxy counts.

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. 2017Katz et al. , 2019;;Pallottini et al. 2017Pallottini et al. , 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.

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. Larger galaxy sample: The current data release contains only ∼ 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 increase the sample size of the data release.This may help capture the even rarer, extreme objects that represent the "tip-of-the-iceberg" of what is currently being observed with the first JWST observations.The limiting factors in providing all galaxies in the initial release are due to computation time and data storage and these will be less problematic moving forward.
2. Extensions beyond star-forming galaxies: A star-formation threshold of 0.3 M ⊙ yr −1 was chosen to limit the galaxies that were being postprocessed to only those that are likely observable with JWST.However, JWST observes flux, not SFR.Without post-processing the data, we are unable to know which galaxies have high enough magnitudes to be observed; hence SFR was our proxy.Unfortunately this choice means that certain classes of galaxies, e.g.quiescent or remnant leakers, are absent from the initial data release, despite the fact that they may represent critical phases in galaxy evolution.Future releases would ideally contain Sphinx galaxies that are in between star formation events rather than at a peak.
3. Alternative dust models: Since Sphinx does not explicitly follow the formation and destruction of dust, we have adopted an effective dust model that assumes a particular SMC attenuation curve and connects the dust content with the neutral gas.While this model is successful in that our UV luminosity and Lyα luminosity functions are in good agreement with observations, it is not clear whether this model is unique in that regard.For this reason, we plan to sample more dust models to better understand their effect on mock observations.

Additional emission lines:
The data set currently contains the vast majority of emission lines that have been observed by JWST at high-redshift; however, there are notable exceptions (for example Mg II λ2796,2803 and N IV λ1483,1486).We have also not provided most of the IR lines that can be observed with ALMA.Such additional emission lines may be part of future data releases.

5.
Varying metal abundance patterns: Because Sphinx does not follow enrichment of individual elements, we have scaled solar abundance patterns to the metallicity of each gas cell.However, it is well known that ratios such as C/O and N/O vary with metallicity.This effect is important for both line ratios and cooling in the ISM and future data releases may attempt to address this effect.
6. Spectral data cubes: While the current data release contains line luminosities of the integrated galaxy spectra, which is comparable to the vast majority of JWST spectral data, one of the unique aspects of JWST is its IFU capabilities.RASCAS has the capabilities of outputting spectral data cubes for any emission line and the limiting factor in their public release is purely data storage.These can be currently made available on-request for individual galaxies and will ideally be standard in future releases.

Lyα IGM propagation:
The Lyα spectra that we released are calculated only out to the virial radius of the galaxy.While this is instructive for understanding how the ISM and CGM modulate the intrinsic Lyα emission, it is not directly comparable to high-redshift galaxy observations where the CGM remains important, even in the postreionization epoch (e.g.Inoue et al. 2014).Although the Sphinx volume does not quite probe the cosmological homogeneity scale (∼ 100 Mpc), we are still able to assess the impact of the IGM on our spectra more locally (see e.g.Garel et al. 2021).
8. Environmental Information: Additional environmental information, such as neighbour properties, filamentary properties, and whether the galaxies reside in ionized bubbles may also be useful for comparing with observations and understanding the role of radiation feedback on galaxy properties (see e.g.Katz et al. 2020).For example, the recent observation that high-redshift Lyα emitters tend to have companions (Witten et al. 2023) would not currently be testable with the Sphinx data set as we only include positions of star-forming galaxies.

Caveats
The Sphinx data release represents a unique tool to both compare theory with observation and to help interpret observations; however, like all numerical simulations, there remain important limitations of the data set.Most importantly, Sphinx has finite spatial and mass resolution limiting how much of the ISM turbulent structure is truly resolved.While the resolution of Sphinx is the state-of-the-art for radiative hydrodynamics simulations of this size, in many cases, Stromgren spheres remain unresolved.This is in many ways a result of the way stars form in the simulation (see discussion in Rosdahl et al. 2015), but the modelling choices have a strong impact on our results.For example, in the case where multiple star particles exist in a cell where the net Stromgren sphere is unresolved, we model this as a single source at the centre of the cell.Thus the ionization parameter is much higher than the case where we model these as separate H II regions embedded in the same cell.Choosing the former can result in higher e.g.[O III]/Hβ than the latter method.Finite spatial resolution also has a strong impact on Lyα and LyC escape physics.It is not yet known how well these quantities are converged with resolution.Furthermore, due to the fixed mass and comoving spatial resolution, high-mass galaxies are better resolved than low mass ones at a fixed redshift and the ISM is better resolved in the higher redshift snapshots compared to the lower redshift ones.The important impact of this numerical effect is that the gas can reach higher densities in more massive galaxies and at high-redshift for a fixed halo mass.The gas density is a fundamental parameter that impacts star formation and much of the emission presented in this work.Intrinsic emission will be more strongly affected as in nearly all galaxies, the resolution in Sphinx is high enough such that the densities impacted by this effect are typically optically thick (at short wavelengths).
One may argue that due to the inevitable use of subgrid models and the imperfect match between simulations and all high-redshift observations, Sphinx is not a completely accurate representation of high-redshift galaxy formation.This is not unique to Sphinx and likely true of all galaxy formation simulations.However, photoionization models that are very simplistic have been the dominant means of interpreting high-redshift spectra.Despite the fact that Sphinx is imperfect, it undoubtedly captures a much more complex ISM and SFHs that better reflect reality than photoionization models, especially since a Sphinx galaxy can be thought of as a complex combination of photoionization models with the additional elements of diffuse gas and a non-zero escape fraction.This is consistent with the fact that groups are now adopting combinations of photoionization models to explain emission lines rather than relying on individual models (e.g.Ramambason et al. 2022).Moreover, Sphinx provides a unique test-bed for understanding how well ISM properties can be inferred from an integrated SED.In this case, it does not matter in which galaxy the ISM resides, this is simply a photon counting exercise for a given density, metallicity, and stellar population distribution.For this reason we argue that the Sphinx data release represents a drastic improvement over previous methods.At worst, our results should be interpreted as more complex, 3D photoionziation models.This will only be superseded by simulations where non-equilibrium metal chemistry is fully coupled to radiative transfer (e.g.Katz 2022) and the Stromgren spheres are resolved (e.g.Kimm et al. 2022).
In addition to subgrid models, Sphinx also relies on a series of modelling assumptions.For example, we have chosen a specific stellar IMF (Kroupa-like), a stellar population synthesis library (BPASS) for the stellar spectra, a particular dust model (SMC), etc. Different choices for these parameters can systematically shift galaxy properties such as mass-to-light ratios, ionizing photon luminosities and therefore emission line luminosities, dust attenuation as a function of wavelength and metallictiy and thus E(B − V) values and β slopes.Once again, this is not only true for Sphinx, but is also true for other numerical simulation as well as SED fitting codes.For these reasons, we have tried to compare directly with observations when possible rather than derived or inferred quantities to mitigate these systematic biases between codes.
5. CONCLUSIONS Here we have described Version 1 of the Sphinx data release, which is publicly available to download from https://github.com/HarleyKatz/SPHINX-20-data.Sphinx is currently the largest full box simulation of galaxy formation in the epoch of reionization with a multiphase ISM and thus represents a valuable resource to the wider community, especially given the recent launch of JWST.The data release contains galaxies spanning ∼ 5 dex in stellar mass and star formation rate that exhibit a diverse set of mass growth and star formation histories.The unique aspect of this data set is the focus on forward modelling observational quantities such that the data can be directly compared with observations.We provide full spectra and photometry including self-consistently calculated nebular emission lines (including Lyα) and continuum, accounting for absorption and scattering by dust (and H I). We have provided a tour of all of the different quantities available in the data set and made comparisons with ERS and Cycle 1 JWST data where possible to demonstrate the utility of the data in interpreting high-redshift observations and to show where constraints can be made on the physics of early galaxy formation.
While the results presented here represent a static view of the Sphinx public data release, we emphasize, as indicated in Section 4.2 that the data set will be updated continuously over time to provide both new data and enhance existing data products.We encourage users to quote the git hash of the data set when publishing with Sphinx data for reproducibility purposes.
The data is organized in a series of CSV tables as well as JSON files that have been indexed by halo ID and redshift.These files can be flexibly loaded into any relevant computing software (e.g.Python, IDL, Julia, etc.).In addition to the curated data, we also release a repository of Juptyer notebooks with basic tutorials on how to load and manipulate the data in Python.Our aim is to lower the barrier to entry for interested users of high-redshift galaxy formation simulations and to ease the comparison between simulations and observations.
Fig. 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.

FigFig
Fig.3.-(Left)Dark matter halo mass functions for haloes that host star-forming galaxies with SFR 10 ≥ 0.3 M ⊙ /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 byDespali et al. (2016) are shown as dotted lines.(Right) Stellar mass distributions of our main sample with SFR 10 ≥ 0.3 M ⊙ /yr for different redshifts (solid lines).The stellar mass distributions for all Sphinx galaxies are shown as dashed lines.Note that M * corresponds to the total stellar mass formed, so does not take into account the mass loss due to stellar evolution.

Fig. 6 .
Fig. 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.
Fig. 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.

FigFig
Fig. 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.
we show the MZR where metallicity is computed along each sight line.Specifically, we use the reddening corrected values of [O III] λ5007, [O III] λ4363, [O II] λλ3727, and Hβ and use the Pilyugin et al. (

Fig. 11
Fig. 11.-Stellar mass/gas-phase metallicity relation for SPHINX 20 galaxies computed in three different ways.The top, middle, and bottom panels show the gas mass-weighted metallicity, the intrinsic [O III] λ5007-, [O II] λλ3727-, and Hβ-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 highredshift galaxies with [O III] λ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] λ4363 detections from Matthee et al. (2023a) based on JWST NIRCam grism data.
Fig. 12.-Histograms of ISM gas densities weighted by intrinsic [O II] λλ3727 (top) or [C III] λλ1908 (bottom) emission as a function of redshift.

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

Fig
Fig. 14.-Histograms of intrinsic ξ ion for the different redshift snapshots available in the Sphinx data release.

Fig. 15
Fig. 15.-Distribution of NIRCam wide filter magnitudes for Sphinx galaxies using their intrinsic magnitudes (Top), dustattenuated 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.

Fig. 16
Fig. 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 fromHarikane et al. (2023).Note that all Sphinx galaxies at z = 10 are above the y-range of the plot.
Fig. 17.-Violin plots of photometric redshifts estimated with EAZY using the default templates (top) or the templates from Larson et al. (2022) (bottom) versus the true redshift of the galaxies in the Sphinx database at z ≥ 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).

Fig. 18
Fig. 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.

Fig. 19 .
Fig. 19.-β measured photometrically versus 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 and the solid lines represent the binned median values.

Fig
Fig. 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 ≤ 8 Bouwens et al. (2021), while circles, triangles, diamonds, stars, cross, and hexagons show JWST measurements at z ≥ 8 from Donnan et al. (2023),Harikane et al. (2023) andBouwens 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.

Fig. 22
Fig. 22.-[O III] λ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.
ley 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] λ5007+Hβ overwhelmingly dominate the filter in the Endsley et al. ( 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. (

Fig. 23
Fig.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 phot > 4 and an 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.

Fig. 24
Fig. 24.-Size-luminosity relation for Sphinx galaxies at z ≥ 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, as well as various highredshift dropout selected galaxies fromMorishita et al. (2023) in magenta.
we show reddening corrected Hα (top) and [O III] λ5007 (bottom) luminosity functions at multiple redshifts compared to JWST observational constraints from Matthee et al. (2023a); Sun et al. (

Fig. 26
Fig. 26.-UV slope, β, 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.
Fig. 27.-Line luminosity functions for Hα (top) and [O III] λ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 ∼ 6 JWST observations from Matthee et al. (2023a); Sun et al. (2022).
Fig. 28.-ObservedHα (top), [O III] λ5007 (middle), and [C II] 158µm (bottom) emission versus halo virial mass.Because [C II] 158µ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.
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

Fig
Fig. 32.-(Left)Lyα and Hα 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.

Fig
Fig.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.

TABLE 1
Statistics of the galaxies in the Sphinx data release.Each column shows the number of galaxies with SFR ≥ 0.3 M ⊙ yr −1 , their maximum and median stellar mass, and the maximum and median virial mass of their host dark matter halo, respectively.Redshift N gal log M * log M * log M vir log M vir

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.

TABLE 3
Details of the galaxy properties available as part of the Sphinx public data release.
Angle-averaged (for all photons with E > 13.6 eV) and along ten sight lines (for photons with a wavelength of 900 Å) ISM gas density log 10 (n H /cm −3 ) Weighted by intrinsic [O II] λλ3727 or [C III] λλ1908 Gas metallicity log 10 (Z/Z ⊙ ) Mass-weighted as well as [O II] λλ3727, [O III] λ5007, [N II] λ6583, and Hβ weighted Emission line luminosities erg s −1