# SPICE: the connection between cosmic reionisation and stellar feedback in
the first galaxies

Aniket Bhagwat1, Tiago Costa1,2, Benedetta Ciardi1, Rüdiger Pakmor1 & Enrico
Garaldi1
1Max Planck Institut für Astrophysik, Karl Schwarzschild Straße 1, D-85741
Garching, Germany
2School of Mathematics, Statistics and Physics, Newcastle University, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present SPICE, a new suite of radiation-hydrodynamic, cosmological
simulations targeting the epoch of reionisation. The goal of these simulations
is to systematically probe a variety of stellar feedback models, including
"bursty" and "smooth" forms of supernova energy injection, as well as poorly-
explored physical scenarios such as hypernova explosions and radiation
pressure on dust. We show that even subtle differences in the behaviour of
supernova feedback drive profound differences in reionisation histories, with
burstier forms of feedback causing earlier reionisation. However, we also find
that some global galaxy properties, such as the dust-attenuated luminosity
functions and star formation main sequence, remain degenerate between models.
In particular, we show that stellar feedback and its strength determine the
morphological mix of galaxies emerging by $z\,=\,5$ and that the reionisation
history is inextricably connected to intrinsic properties such as galaxy
kinematics and morphology. While star-forming, massive disks are prevalent if
supernova feedback is "smooth", "bursty" feedback preferentially generates
dispersion-dominated systems. Different modes of feedback produce different
strengths of outflows, altering the interstellar/circumgalactic medium in
different ways, and in turn strongly affecting the escape of Lyman continuum
(LyC) photons. We establish a correlation between galaxy morphology and LyC
escape fraction, revealing that dispersion-dominated systems have escape
fractions 10-50 times higher than their rotation-dominated counterparts at all
redshifts. At the same intrinsic luminosity, dispersion-dominated systems
should thus preferentially generate large HII regions as compared to their
rotation-dominated counterparts. Since dispersion-dominated systems are more
prevalent if stellar feedback is more explosive, reionisation occurs earlier
in our simulation with burstier feedback. We argue that statistical samples of
post-reionisation galaxy morphologies (using both stellar and gaseous
components) probed with telescopes such as JWST, ALMA and MUSE can constrain
stellar feedback at $z>5$ and models of cosmic reionisation.

###### keywords:

Galaxies : – reionisation, galaxy evolution, ISM, feedback – numerical :
radiative transfer – early universe

††pubyear: 2015††pagerange: SPICE: the connection between cosmic reionisation
and stellar feedback in the first galaxies–SPICE: the connection between
cosmic reionisation and stellar feedback in the first galaxies

##  1 Introduction

About $300,000\,\rm yr$ after the Big Bang, protons and free electrons
coalesce for the first time, leading to the formation of hydrogen and helium
atoms. At this time, the Universe is mostly neutral. By $z\,=\,5$ (e.g.
Gaikwad et al., 2023), however, the vast majority of hydrogen has transitioned
to an ionised state. This $\sim 1\,\rm Gyr$ period is known as the epoch of
reionisation (EoR). Reionisation is brought about by the ionising flux
generated by stellar populations in the first galaxies (Shapiro & Giroux,
1987; Madau et al., 1999; Gnedin, 2000; Robertson et al., 2010, 2015; Eide et
al., 2020), with quasars now thought to provide only a minor contribution
(Madau & Haardt, 2015; Mason et al., 2019b). The ionising photon budget is set
by the abundance of young stars, while the ability of photons to escape from
galaxies is governed by the structure of the interstellar medium (ISM), which,
in turn, is shaped by a multitude of ‘feedback’ processes associated to star
formation. Within the $\Lambda$CDM paradigm of galaxy formation, such
processes must be invoked in order to prevent the overproduction of stars even
at $z>6$ (Hopkins et al., 2014; Costa et al., 2014; Somerville & Davé, 2015).

These feedback processes include massive stellar winds (Weaver et al., 1977;
Mackey et al., 2015; Geen et al., 2020; Lancaster et al., 2021; Guszejnov et
al., 2022), photo-ionisation and photo-heating (Geen et al., 2015; Peters et
al., 2017a; Kim & Ostriker, 2018), radiation pressure (Murray et al., 2005;
Menon et al., 2022) and supernova explosions (Dekel & Silk, 1986). The impact
of stellar feedback on the phase structure of the ISM is well established
(McKee & Ostriker, 1977; Chevalier & Clegg, 1985; Kimm & Cen, 2014; Martizzi
et al., 2016). Feedback in the form of ionising radiation suppresses radiative
cooling (Efstathiou, 1992; Gnedin, 2000; Kim et al., 2017). Strong radiation
fields can also reduce the baryon fractions of galaxies via photo-heating,
further suppressing growth (Wise & Abel, 2008; Okamoto et al., 2008; Hasegawa
& Semelin, 2012; Wise et al., 2012; Pawlik et al., 2013). Supernova explosions
heat the ISM (McKee & Ostriker, 1977) and launch galactic outflows, expelling
ISM material that would otherwise form stars (e.g. Dubois & Teyssier, 2008;
Puchwein & Springel, 2013; Fielding et al., 2018; Orr et al., 2019; Martizzi,
2019). Theoretical studies have for decades emphasised the importance of
turbulence in shaping the ISM. Turbulence can be driven by processes such as
gravitational instability (Klessen & Ballesteros-Paredes, 2004) and feedback
from supernovae (Larson, 1981; Solomon et al., 1987; Heyer & Brunt, 2004;
Federrath, 2016). Feedback-driven turbulence in the ISM can both inhibit and
drive star formation by either providing support against gravitational
collapse, or by creating strong density contrasts allowing for rapid collapse,
hence strongly regulating star formation efficiency (Krumholz & McKee, 2005;
Ostriker et al., 2010; Ostriker & Shetty, 2011; Faucher-Giguère et al., 2013;
Gatto et al., 2016).

By influencing the structure of the ISM, stellar feedback plays a direct role
in cosmic reionisation. Recent studies highlight the importance of supernova-
driven outflows in facilitating Lyman continuum (LyC, photons with energy
$>13.6$eV) escape (Trebitsch et al., 2017, 2018; Rosdahl et al., 2022a)
through the creation of low-density channels. The impact of stellar radiation
itself on the ISM, however, remains less clear. While studies agree that
photo-ionisation feedback suppresses star formation bursts by counteracting
the formation of high-density peaks in the ISM (Rosdahl et al., 2015; Peters
et al., 2017b; Haid et al., 2018), the role of radiation pressure is less
settled. Analytical arguments (Murray et al., 2005; Thompson et al., 2015)
suggest that radiation pressure on dust should launch galactic winds if
systems are sufficiently bright. Radiation-hydrodynamic (RHD) simulations,
however, generally find that early stellar feedback (such as radiation
pressure) acts to suppress outflows through a reduction in star formation and
in supernova clustering (Kimm et al., 2018; Costa et al., 2019; Agertz et al.,
2020; Smith et al., 2021). Recent observational evidence points to an
‘effective Eddington limit’ in star-forming galaxies at $z>6.5$, observed
through an absence of systems with high star formation rates and high optical
depths (Fiore et al., 2023). Though possibly not a unique interpretation, such
trends may suggest a link between dust obscuration and the strength of
galactic outflows that is not captured by current models.

Besides questions surrounding the key feedback driving mechanisms, there are
significant numerical uncertainties in the modelling of feedback in galaxy
evolution simulations. An ab initio treatment of supernova feedback remains
challenging in cosmological simulations due to prohibitive resolution
requirements. Cosmological simulations thus have to resort to ‘subgrid’ models
attempting to capture the impact of supernova feedback at the resolution
scale, i.e. $\gtrsim 20\,\rm pc$111In this work, we denote comoving
coordinates as ”cpc, ckpc” etc and physical coordinates as ”pc, kpc” etc..
While such subgrid models have become more sophisticated in the last decade
(Rosdahl et al., 2017; Hopkins et al., 2018), uncertainties persist. For
instance, while commonly adopted strategies account for unresolved ‘PdV’ work
through a momentum boost (e.g. Kimm et al., 2015), current models may
underproduce hot gas, essential for launching galactic outflows, if supernova
explosions are not resolved (Hu, 2019). Even when enforcing the correct
terminal momentum, many existing models have to boost supernova feedback in
order to reproduce realistic galaxy star formation rates and masses. Missing
physics, such as cosmic ray injection and transport (Diesing & Caprioli,
2018), have been proposed as possible ways to further strengthen stellar
feedback (e.g. Martin-Alvarez et al., 2022).

Another possible approach to overcome limitations caused by insufficient
resolution and decouple the impact of stellar radiation and supernova feedback
from the host’s ISM, is through the adoption of an effective model for
galactic winds (Springel & Hernquist, 2003; Pillepich et al., 2017), with a
prescribed mass outflow rate and velocity. While this approach promotes good
numerical convergence properties, it limits the simulations’ insight into the
detailed interaction between supernovae, stellar radiation and the ISM and
properties of host galaxies, restricting its effects to the scale of the
intergalactic medium (IGM). It also introduces difficulties in the generation
of mock observables, such as Ly$\alpha$ or X-ray emission, to which the ISM
contributes significantly.

The recent influx of observations from JWST has begun to provide statistical
samples of galaxies at unprecedented resolution deep into reionisation. Probes
such as spectral energy distribution (SED) fitting, emission line ratios and
UV/H$\alpha$ SFR indicators help constrain star formation properties of high
redshift galaxies. While a large number of studies find that bursty star
formation (Atek et al., 2023; Dressler et al., 2023a; Tacchella et al., 2023;
Langeroodi & Hjorth, 2023; Endsley et al., 2023a; Asada et al., 2023a)
dominates the histories of $M_{*}\leq 10^{9}$ M⊙ galaxies, there is also
evidence of a smoother star formation channel or even a combination of the two
(Ciesla et al., 2023; Dressler et al., 2023b). Recently, a bursty star
formation history (SFH) has been invoked to alleviate (Sun et al., 2023a; Sun
et al., 2023b; Steinhardt et al., 2023) the so-called "too many too bright"
galaxies problem (Ferrara et al., 2023; Boylan-Kolchin, 2023). Other studies
find that galaxies with bursty SFH could be the main drivers of reionisation
(Simmonds et al., 2023a; Endsley et al., 2023a; Endsley et al., 2023b; Atek et
al., 2023). Therefore, understanding the implications of different SFH in high
redshifts galaxies is key, and careful theoretical modelling will help to
interpret various observations. Theoretical studies such as Hartley & Ricotti
(2016) and Furlanetto & Mirocha (2022) use semi-analytical models to
understand the effect of a bursty SFH on reionisation, finding that the sizes
of HII regions are strongly modulated by bursts.

Over the last decade, RHD simulations of reionisation such as CROC (Gnedin,
2014), Technicolor Dawn (Finlator et al., 2018), SPHINX (Rosdahl et al., 2018;
Katz et al., 2018; Rosdahl et al., 2022a), CoDa (Ocvirk et al., 2016, 2020;
Lewis et al., 2022) and THESAN (Kannan et al., 2022; Garaldi et al., 2022;
Smith et al., 2022) have used variations of the stellar feedback prescriptions
previously described to simulate the high redshift universe (see Gnedin &
Madau 2022 for a detailed review). Different simulations adopt varied
approaches toward their modelling, with simulations like THESAN, CoDa and CROC
focusing on the large scale reionisation process with fairly large volumes
(though unable to resolve galaxy scale heights), whereas SPHINX focuses on
resolving the internal structure of galaxies while compromising on volume. Due
to the extreme cost of RHD simulations (tens of millions of cpu hours), most
are performed for only a single fiducial model at the target resolution. The
effect of variations of input baryonic physics has remained comparatively
unexplored.

Available RHD simulations provide either a statistical sample of galaxies
modeled within the same feedback prescription, or a variety of models for
single, zoom-in galaxies (Pallottini & Ferrara, 2023; Katz et al., 2022). A
systematic and statistical study of the effect of different feedback models is
thus missing. The approach taken in this paper is to model galaxy formation
and reionisation such that we maximize how well we resolve the multi-phase ISM
(down to $\approx 28\rm pc$ at $z=5$) in a cosmological volume ($L_{\rm
box}\approx 14.8\rm cMpc$). The emphasis of these simulations is on the
variations in supernova feedback prescriptions in order to quantify the impact
of numerical uncertainties and to uncover observationally-testable connections
between stellar feedback and the properties of the galaxies that drive
reionisation. We introduce SPICE, a suite of simulations which include
mechanical feedback from supernovae, stellar feedback in the form of ionising
radiation (with on-the-fly radiation transport) and radiation pressure. Our
simulations explore different models for supernova feedback with variations in
explosion timing, supernova energies and the presence of hypernovae, predicted
to exist in the pristine conditions of the early universe (Kobayashi et al.,
2006). Within the SPICE framework, we aim to understand how the effect of
different modes of feedback manifests in observable properties of galaxies.
Indeed, by varying only the feedback prescription while keeping everything
else constant, we are able to quantify _relative differences_ between the
feedback models.

The paper is organised as follows. In Section 2 we describe the SPICE suite of
simulations setup along with the physics prescriptions that are included. In
Section 3 we examine our results on the effect of stellar feedback on star
formation, reionisation, UV luminosity functions, galaxy morphologies and LyC
escape fractions. In Sections 4 and 5 we discuss and summarize our findings,
respectively.

##  2 Simulations

In this section, we describe the simulations performed and post-processed for
this study. We perform a total of three flagship simulations to study the
effect of variations in supernova feedback. Relevant global simulation
parameters are summarized in Table 1, while individual variations in feedback
processes are briefly described in Table 2.

All simulations are performed with RAMSES-RT (Rosdahl et al., 2013; Rosdahl &
Teyssier, 2015), which is a radiation-hydrodynamics extension of the Eulerian
adaptive mesh refinement code RAMSES (Teyssier, 2002). RAMSES-RT solves the
coupled gas hydrodynamics, radiative transfer of stellar radiation along with
the self-gravity and N-body dynamics of gas, dark matter and stars. RAMSES-RT
uses a first-order Godunov method with the M1 closure for the Eddington tensor
(Levermore, 1984) to solve the radiation transport equation. We employ the
’Global–Lax–Friedrich’ (GLF) Riemann solver to advect radiation between cells.

The adaptive mesh nature of the code allows for the grid to be dynamically
refined in order to obtain higher numerical resolution within sub-regions of
the simulation domain. We summarize the models used for our simulations below,
with sections 2.1-2.5 describing global setup within RAMSES-RT and sections
2.6-2.11 describing the novel combination of galaxy formation and feedback
models we include in our simulations.

###  2.1 Initial conditions

The initial conditions (ICs) for the simulations are generated at $z=30$ using
MUSIC2-monofonIC222https://bitbucket.org/ohahn/monofonic/ (Hahn et al., 2020;
Michaux et al., 2020) with the 2PPT/2LPT approach. A $\rm\Lambda CDM$
cosmological model is adopted, with parameters $\Omega_{\Lambda}=0.6901$,
$\Omega_{\rm m}=0.3099$, $\Omega_{\rm b}=0.0489$, $H_{0}=67.74$ km/s/Mpc,
$n_{\rm s}=0.9682$, $\sigma_{8}=0.8159$ (Planck Collaboration et al., 2016).
We assume primordial gas mass fraction $X=0.7545$ for Hydrogen (H) and
$Y=0.2455$ for Helium (He), along with a metallicity floor (see 2.3). All
simulation boxes have a length $L_{\rm box}=10$ cMpc/$h$, with $512^{3}$ dark
matter particles, i.e. a mean mass $m_{\rm dm}=6.38\times 10^{5}$ M⊙.
MUSIC2-monofonIC treats dark matter and baryons within the two-fluid
approximation, using the novel 2nd order propagator perturbation theory (2PPT)
(Uhlemann et al., 2019; Rampf et al., 2020). The 2PPT approach involves
perturbing particle masses with a distribution centered at $m_{\rm dm}$ (see
Hahn et al. 2020) which suppresses discreteness errors. RAMSES-RT evolves
baryons from density and momentum fields which are disretized at fixed
locations in Eulerian space. Traditionally, these Eulerian fields are
generated from Lagrangian displacements and interpolated back onto Eulerian
grids, thus introducing interpolation erros. The 2PPT approach directly yields
these fields without ad-hoc interpolation (Porqueres et al., 2020), therefore
providing accurate ICs for codes like RAMSES-RT.

Name | Value | Description
---|---|---
$L_{\rm box}$ | 10 [cMpc/h] | Box size
$m_{\rm dm}$ |  $6.38\times 10^{5}$ [$\rm\rm M_{\odot}$] | Mean mass of dark matter particles
$N_{\rm dm}$ | $512^{3}$ | Number of dark matter particles
$m_{\ast}$ | 975 [$\rm\rm M_{\odot}$] | Minimum mass of stellar particles
$Z_{\rm floor}$ |  $3.2\times 10^{-4}$ $Z_{\odot}$ | Metallicity floor at initial redshift
$f_{\rm esc}$ | 1.0 | Stellar birth cloud escape fraction
Table 1: Basic properties of the simulated cosmological boxes, including the variable name, its adopted value and a short description. Name | $\rm t_{SNe}$ |  $\rm E_{SN/HN}$ (erg) | Description
---|---|---|---
bursty-sn | 10 Myr | 2$\times 10^{51}$ | All SN in a single stellar particle explode at the same time and have the same energy
smooth-sn | 3-40 Myr | 2$\times 10^{51}$ | The time delay of SN is stochastically sampled, but all SN have the same energy
hyper-sn | 3-40 Myr | $10^{50}-2\times 10^{51}(\rm SN)+10^{52}(\rm HN)$ | The time delay and energy of SN is stochastically sampled + metallicity-dependent HN rate
Table 2: Description of various simulations and supernova feedback model
variations. The columns list (from left to right) the name of the simulation,
delay time for supernova events, energy per supernova event, brief description
of the simulation.

###  2.2 Gravity and hydrodynamics

The Euler hydrodynamic equations are solved employing a second order Godonov
scheme based on a MUSCL-Hancock method. We use the "HLLC" Riemann solver (Toro
et al., 1994) to evaluate fluxes across all interfaces, and the MinMod slope
limiter to construct gas variables at cell interfaces from their cell-centred
values. To close the relation between gas pressure and internal energy, we
adopt an adiabatic index $\gamma$ = 5/3 (for an ideal monatomic gas). The
dynamics of collisionless dark matter and stellar particles are computed using
the Poisson equation with a particle-mesh solver. The dark matter and stellar
particles are projected onto a grid with a cloud-in-cell interpolation
(Guillet & Teyssier, 2011). We employ a multigrid solver (Guillet & Teyssier,
2011) to solve the Poisson equation up to a refinement level of 14 ($\Delta
x\approx 100\rm pc$), while at more refined levels we adopt a conjugate
gradient solver to improve the computational speed.

###  2.3 Cooling and gas thermochemistry

Gas cooling and heating is calculated as described in Rosdahl et al.
(2013)333Metal cooling is not currently modelled self-consistently with the
local radiation, but see Katz (2022) for a development in this direction.. We
evaluate the non-equilibrium ionisation states of H and He (HII , HeII ,
HeIII) in each computational cell, and advect them with the gas, as passive
scalars (details of the quasi-implicit method can be found in Rosdahl et al.
2013). They are fully coupled to the local ionising radiation field and for
these primordial species, we calculate cooling and heating due to
Bremsstrahlung, photoionisation, collisional ionisation, collisional
excitation, Compton cooling off the cosmic microwave background, and di-
electronic recombination.

The cooling contribution from metals at $T>10^{4}$ K is computed using tables
generated with CLOUDY (Ferland et al. 1998, version 6.02), assuming
photoionisation equilibrium with a Haardt & Madau (1995) UV background.
Instead, for $T\leq 10^{4}$ K we adopt the fine structure cooling rates from
Rosen & Bregman (1995), allowing the gas to cool radiatively down to $\approx
300$K. We also assume a homogeneous initial gas metallicity floor of
$Z=3.2\times 10^{-4}$ Z⊙, which is used to mimic missing molecular hydrogen
cooling channels and metal enrichment from Pop-III stars (Wise et al., 2011),
and allows the gas to cool below $T=10^{4}$ K. We also adopt the on-the-spot
approximation, where emission of ionising photons from recombining gas is
ignored, i.e. we assume that it is all absorbed locally, within the same cell.

###  2.4 Refinement strategy

The adaptive mesh refinement nature of RAMSES-RT allows for each parent gas
cell to be split into 8 cells when certain conditions are satisfied (see
below). The cell refinement level $\ell$ determines the width of each gas cell
as $\Delta x_{\ell}\ =\ 0.5^{\ell}L_{\rm box}$. We allow refinement levels in
the range $\ell=9-16$. At the coarsest level ($\ell_{\rm min}=9$) we have a
minimum physical resolution of 4.8 kpc at $z=5$, while at the finest level
($\ell_{\rm max}=16$) the maximum physical resolution is $\approx 28$ pc at
$z=5$. Cells start at the coarsest level and are adaptively refined to higher
levels to increase the numerical resolution of the simulation. A cell is
refined if the following criteria are met:

  1. 1.

if $M_{\rm dm}+\frac{\Omega_{m}}{\Omega{b}}\ M_{\rm b}>8\times m_{\rm dm}$
(quasi-Lagrangian refinement criterion), where $M_{\rm dm}$ and $M_{\rm b}$
are the total dark matter and baryonic (i.e. gas + stars) mass within the
cell, and $m_{\rm dm}$ is the mean dark matter particle mass.

  2. 2.

The cell width is larger than $\frac{1}{8}$ of the local Jeans length.

We adopt a constant comoving resolution throughout time, meaning that we allow
for cells to be refined to $\ell_{\rm max}$ at all redshifts. This allows for
the maximum physical resolution to be the highest at very high redshifts
($\approx$10.72 pc at $z\,=\,20$).

###  2.5 Halo and galaxy finding

To identify dark matter haloes and galaxies, we use the AdaptaHOP halo finder
(Aubert et al., 2004; Tweed et al., 2009) in the most massive submaxima mode.
We fit a triaxial ellipsoid to each (sub-)halo and check that the virial
theorem is satisfied within this ellipsoid, with the center corresponding to
the location of the densest particle. If this condition is not satisfied, we
iteratively decrease its volume until we reach an inner virialized region.
From the volume of this largest ellipsoidal virialized region, we define the
virial radius $R_{\rm vir}$ and mass $M_{\rm vir}$. For the halo finder, we
require a minimum of 100 particles per halo. Halo finder parameters (as
defined in Appendix B of Aubert et al. 2004) used are: $N_{\rm SPH}=32$,
$N_{\rm HOP}=16$, $\rho_{\rm th}=80$, and $f_{\rm poisson}=4$. We also
identify galaxies with AdaptaHOP. Here we require at least 100 stellar
particles per galaxy. The galaxy-finder parameters are: $N_{\rm SPH}=10$,
$N_{\rm HOP}=10$, $\rho_{\rm th}=10^{3}$ , and $f_{\rm poisson}$ = 4.

###  2.6 Star formation

We employ the multi-freefall star formation model implemented by Kretschmer &
Teyssier (2020). This includes a model for subgrid turbulence (as described in
section 2.7), yielding a variable star formation efficiency that depends on
local conditions described in terms of the virial parameter and the turbulent
Mach number (see below). In each gas cell individually where the gas density
is above a threshold of $n_{\rm H}\geq 10\rm cm^{-3}$, we adopt a standard
Schmidt law, for which the star formation density is given as,

$\centering\dot{\rho_{*}}=\epsilon_{\rm ff}\frac{\rho}{t_{\rm
ff}},\@add@centering$ (1)

where $\rho$ is the gas density, $t_{\rm ff}=\sqrt{3\pi/(32G\rho)}$ is the
freefall time, and $\epsilon_{\rm ff}$ is the star formation efficiency per
freefall time. Typically, in galaxy formation simulations the value of
$\epsilon_{\rm ff}$ is assumed to be a constant in the range 1-3 percent
(Agertz et al., 2011), motivated by observations of inefficient star formation
on galactic scales (Bigiel et al., 2008), as well as in Milky Way giant
molecular clouds (GMCs) (Krumholz & Tan, 2007; Murray, 2011; Lee et al., 2016;
Grisdale et al., 2019). This model can produce local (i.e. within a single
computational cell) $\epsilon_{\rm ff}$ values as low as 0.1$\%$ likely
suppressed by feedback (Ostriker & Shetty, 2011; Kim & Ostriker, 2015a), while
also reaching values $>100\%$ (i.e. star formation faster than freefall time).
Feedback modulates star formation in different parts of the galaxy such that
the global (averaged over the entire galaxy) $\epsilon_{\rm ff}$ values
produced are in agreement with observations (see Fig. 6 in Kretschmer &
Teyssier 2020).

We assume that the gas density within a supersonic turbulent medium such as
the ISM is described by a log-normal probability distribution function (PDF)
given as

$p(s)=\frac{1}{\sqrt{2\pi\sigma_{s}^{2}}}{\rm
exp}\frac{(s-\overline{s})^{2}}{2\sigma_{s}^{2}},$ (2)

where $s=\rm ln(\rho/\overline{\rho})$, $\sigma_{s}$ is the variance of $s$
and $\overline{\rho}$ is the mean density of the cell.

The local star formation efficiency $\epsilon_{\rm ff}$, as computed in
Hennebelle & Chabrier (2011); Federrath & Klessen (2012), is

$\displaystyle\epsilon_{\rm ff}$
$\displaystyle=\frac{\epsilon}{\phi}\int_{s_{\rm crit}}^{\infty}\frac{t_{\rm
ff}(\overline{\rho})}{t_{\rm ff}(\rho)}\frac{\rho}{\overline{\rho}}p(s)ds$ (3)
$\displaystyle=\frac{\epsilon}{2\phi}{\rm
exp}\left(\frac{3}{8}\sigma_{s}^{2}\right)\left[1\ +\ {\rm
erf}\left(\frac{\sigma_{s}^{2}-s_{\rm
crit}}{\sqrt{2\sigma_{s}^{2}}}\right)\right],$ (4)

where we assume that gas with density larger than a critical value, $s_{\rm
crit}$, is converted into stars. The parameter $\phi=1/0.49$ takes into
account the uncertainty in free-fall timescales for gas with different
densities and $\epsilon=0.5$ accounts for the fact that not all gas above
$s_{\rm crit}$ is converted into stars (both quantities are derived in
Federrath & Klessen 2012). The critical density for star formation ($s_{\rm
crit}$) is calculated using the model from Krumholz & McKee (2005). In order
to extend the model to the subsonic regime, we allow star formation when a
cell is gravitationally unstable ($\alpha_{\rm vir}<1$) and $\mathcal{M}<1$.
The critical density accounting for both supersonic and subsonic regimes is

$s_{\rm crit}={\rm ln}\left[\alpha_{\rm
vir}\left(1+\frac{2\mathcal{M}^{4}}{1+\mathcal{M}^{2}}\right)\right]$ (5)

Here,

$\alpha_{\rm vir}=\frac{5\sigma^{2}}{3G\rho\Delta x^{2}}$ (6)

is the local virial parameter, which is an indicator of local stability.
$\Delta x$ is the cell size, $G$ is the gravitational constant, $\sigma$ is
the 1D turbulent velocity (see 2.7) and $\mathcal{M}\ =\ \sigma/c_{s}$ is the
local Mach number, where $c_{s}$ is the speed of sound.

This model therefore allows for two star formation channels. In the first
channel, for which $\alpha_{\rm vir}<1$ (independent of local $\mathcal{M}$),
the entire computational cell is gravitationally unstable. In the second
channel, collapse owing to large density fluctuations caused by supersonic
turbulence, i.e. $\mathcal{M}\gg 1$ (see Figure 1 in Kretschmer & Teyssier
2020), becomes possible even if $\alpha_{\rm vir}>1$. For reference, typical
conditions for star-forming regions within the Milky way are $\alpha_{\rm
vir}\simeq 3-10$ with $\mathcal{M}\simeq 10-20$ (Spilker et al., 2022), and
correspond to $\epsilon_{\rm ff}\simeq 0.01-0.02$ (Murray, 2011; Grisdale et
al., 2019).

Using the local star formation efficiency and the gas mass within the cell, we
compute the number of new stellar particles in a given timestep by sampling a
Poisson distribution. The initial stellar mass of new stars is always an
integer multiple of 970 M⊙ ($\approx 95-97\%$ of stars in a simulation have a
mass of 970 M⊙), but it is capped such that we do not deplete more than 90% of
the mass of the star-forming cell. We also assume that each stellar population
represents a fully sampled Chabrier (2003) initial mass function (IMF).

###  2.7 Subgrid turbulence

Turbulence is modelled using a sub-grid scale approach (Schmidt et al., 2006).
In the implementation of Kretschmer & Teyssier (2020), an additional equation
for the turbulent kinetic energy is adopted to account for advection of
turbulent energy and work done by turbulent pressure (Schmidt 2014 and Semenov
et al. 2018). This is given by

$\frac{\partial}{\partial t}K_{\rm T}+\frac{\partial}{\partial x_{j}}(K_{\rm
T}\tilde{v}_{j})+P_{\rm T}\frac{\partial\tilde{v}_{j}}{\partial x_{j}}=C_{\rm
T}-D_{\rm T},$ (7)

where the turbulent kinetic energy $K_{\rm T}=3/2\rho\sigma^{2}_{\rm T}$, is
related to the turbulent pressure by $P_{\rm T}=2/3K_{\rm T}$, with
$\sigma_{\rm T}$ representing the 1D turbulent velocity dispersion. The left-
hand side of the equation includes terms describing the time evolution,
advection and compression. The right-hand side contains creation and
destruction terms ($C_{T}$ and $D_{T}$, respectively), which can be written as

$C_{T}=2\mu_{T}\sum_{ij}\Bigg{[}\frac{1}{2}\Biggl{(}\frac{\partial\tilde{v}_{i}}{\partial
x_{j}}+\frac{\partial\tilde{v}_{j}}{\partial
x_{i}}\Biggl{)}-\frac{1}{3}\bigl{(}\nabla\cdot\tilde{v}\bigl{)}\delta_{ij}\Bigg{]}^{2}\hskip
1.49994pt=\hskip 1.49994pt\frac{1}{2}\mu_{T}\big{|}S_{ij}\big{|}^{2},$ (8)

and

$D_{T}=\frac{K_{T}}{\tau_{\rm diss}},$ (9)

where $\big{|}S_{ij}\big{|}$ is the mean flow viscous stress tensor (this term
is evaluated as in Schmidt & Federrath 2011), $\tilde{v}$ is the mean flow
variable (defined as $\tilde{v}=\overline{\rho v}/\overline{\rho}$) and the
sum is calculated over the nearest neighbours of the computational cell in
question. This model has two important parameters, the turbulent viscosity
$\mu_{T}$ and the dissipation time-scale $\tau_{\rm diss}$ , which are related
to the cells’ size by

$\mu_{T}=\overline{\rho}\Delta x\sigma\hskip 10.00002pt{\rm and}\hskip
10.00002ptt_{\rm diss}=\frac{\Delta x}{\sigma}.$ (10)

Previous implementations of turbulent star formation sub-grid models consider
an in-situ calculation (Kimm et al., 2017) of the turbulent velocity
dispersion (Perret et al., 2015; Trebitsch et al., 2017, 2018; Hopkins et al.,
2018), which can be thought of as the stationary limit of the model used in
this work (i.e. the creation and destruction terms are considered to be
equal). In our case, the model accounts for the non-equilibrium dissipation of
the turbulent kinetic energy using the density and velocity fields without
modifying the hydrodynamic solver. Turbulent velocity dispersion is estimated
by solving Eq. 7 such that $\sigma=\sqrt{2K_{\rm T}/\overline{\rho}}$. The
resulting velocity dispersion is adopted in Eq. 6 to evaluate the Mach number
and the virial parameter entering the star formation efficiency computation
(Eq. 4).

###  2.8 Supernova feedback

Photon group |  $\epsilon_{0}$ [eV] |  $\epsilon_{1}$ [eV] |  $\langle\epsilon\rangle$ [eV] |  $\rm\sigma_{HI}$ [$\rm cm^{2}$] |  $\rm\sigma_{HeI}$ [$\rm cm^{2}$] |  $\rm\sigma_{HeII}$ [$\rm cm^{2}$] |  $\rm\kappa_{abs}^{0}$ [$\rm g\ cm^{-2}$] |  $\rm\kappa_{scat}^{0}$ [$\rm g\ cm^{-2}$]
---|---|---|---|---|---|---|---|---
Infrared | 0.1 | 1 | 0.61 | 0 | 0 | 0 | 0 | 10
Optical | 1 | 13.6 | 5.51 | 0 | 0 | 0 | $10^{3}$ | 0
UV I | 13.6 | 24.59 | 17.9 | $3.19\times 10^{-18}$ | 0 | 0 | $10^{3}$ | 0
UV II | 24.59 | 54.42 | 32.32 | $6.27\times 10^{-19}$ | $4.78\times 10^{-18}$ | 0 | $10^{3}$ | 0
UV III | 54.42 | $\infty$ | 67.77 | $7.37\times 10^{-20}$ | $1.12\times 10^{-18}$ | $8.71\times 10^{-19}$ | $10^{3}$ | 0
Table 3: Columns left to right: Photon groups, the lower and upper energies
defining their energy interval, mean photon group energies , ionisation cross-
sections to HI, HeI and HeII, absorption and scattering opacities to dust for
each photon group.

We adopt the mechanical feedback scheme introduced by Kim & Ostriker (2015b)
and implemented as in Kimm et al. (2015) (see also Kimm et al. 2016; Kimm et
al. 2017; Rosdahl et al. 2017; Dubois et al. 2021, and a similar setup in
Hopkins et al. 2014). The key idea behind this model is to correctly capture
the terminal momentum associated to the snowplow phase of a SN remnant, by
injecting into the surrounding cells a radial momentum $p_{\rm SN}$ which
depends on whether the adiabatic (Sedov-Taylor) phase is resolved. The model
introduces a parameter $\chi$, which gives the ratio between the swept-up mass
$M_{\rm swept}$ and the ejected mass $M_{\rm ej}$ for each neighbouring cell
(number of neighbouring cells $N_{\rm nbor}$ = 48) as

$\chi\equiv dM_{\rm swept}/dM_{\rm ej}.$ (11)

This is compared to a threshold value of

$\chi_{\rm tr}\equiv 69.58\ n_{\rm H}^{-4/17}E_{\rm
51}^{-2/17}Z^{\prime-0.28},$ (12)

where $E_{51}$ is the explosion energy of an individual SN in units of
$10^{51}$ erg, $n_{\rm H}$ is the hydrogen number density in $\rm cm^{-3}$ and
$Z^{\prime}={\rm{max}}[0.01,Z/{\rm Z}_{\odot}]$. Additionally,

$dM_{\rm ej}=(1-\beta_{\rm sn})M_{\rm ej}/N_{\rm nbor},$ (13)

and

$dM_{\rm swept}=\rho_{\rm nbor}\biggl{(}\frac{\Delta
x}{2}\biggl{)}^{2}+\frac{(1-\beta_{\rm sn})\rho_{\rm host}\Delta x^{3}}{N_{\rm
nbor}}+dM_{\rm ej}.$ (14)

Here $\rho_{\rm host}$ is gas density of the SN host cell, and $\beta_{\rm
sn}$ determines what fraction of the gas mass ($M_{\rm ej}\ +\ \rho_{\rm
host}\Delta x^{3}$) is re-distributed to the host cell of a SN (see Fig. 15 of
Kimm & Cen 2014). The parameter $\beta_{\rm sn}$ is set to 4/52 to distribute
the mass as evenly as possible to the host and neighbouring cells when the
cells are on the same level of refinement.

If $\chi>\chi_{\rm tr}$ the adiabatic phase is not resolved and the momentum
during the snowplow phase ($p_{\rm SN,snow}$) is injected to the neighbouring
cells. Otherwise the momentum during the adiabatic phase ($p_{\rm SN,ad}$) is
injected, i.e.

$p_{\rm SN}=\begin{cases}p_{\rm SN,ad}=\sqrt{2\chi M_{\rm ej}f_{\rm
e}({\rm\chi})E_{\rm SN}}&\hskip 20.00003pt\text{$\chi<\chi_{\rm tr}$\,
\rm,}\\\ p_{\rm SN,snow}&\hskip
20.00003pt\text{$\chi\geq\chi_{tr}\,,$}\end{cases}$ (15)

where $f_{\rm e}(\rm\chi)\ =\ 1-\frac{\chi-1}{3(\chi_{\rm tr}-1)}$ smoothly
connects the two regimes. The input momentum during the snowplow phase is

$p_{\rm{SN},\rm{snow}}=3\times 10^{5}\ {\rm{km\ s^{-1}\rm M_{\odot}}}\
E^{16/17}_{51}\ n_{\rm{H}}^{-2/17}\ Z^{\prime-0.14},$ (16)

following Blondin et al. (1998); Thornton et al. (1998); Geen et al. (2015);
Kim & Ostriker (2015b). Geen et al. (2015) have shown that the final radial
momentum from SN in the snowplow phase can be augmented by photoionisation of
their environments by massive stars. Following their approach, for each SN
event, the local Strömgen radius $r_{\rm S}$ is calculated and compared with
the cell width $\Delta x$. If the Strömgen radius is not well resolved
($\Delta x>>r_{\rm s}$), we adopt a final input momentum of

$p_{\rm{SN,snow}}=p_{\rm SN,snow}{\rm exp}\biggl{(}-\frac{\Delta x}{r_{\rm
s}}\biggl{)}\ +\ p_{\rm SN+PH}\biggl{(}1-{\rm exp}\biggl{[}-\frac{\Delta
x}{r_{\rm s}}\biggl{]}\biggl{)}$ (17)

where $p_{\rm SN+PH}$ is obtained using a fit to the results of Geen et al.
(2015).

Our assumption of a Chabrier (2003) IMF implies that a fraction $\eta_{\rm
sn}=0.31$ of the initial mass of the stellar particle is recycled back into
the ISM. Out of the recycled mass, a fraction of $0.05$ is in form of metals
heavier than He (Chabrier, 2003). Due to computational cost, we do not store
individual element abundances, but rather track the ejected metal mass as a
metallicity per cell. Our choice of the IMF gives a fixed supernova rate of
0.01639 SN M${}_{\odot}^{-1}$, which results in $15$ SN per star particle for
our minimum stellar mass of 970 M⊙.

###  2.9 Supernova feedback variations

Given significant uncertainties in numerical models for supernova feedback
(see Section 1), we explore a number of variations in the timing and energy
injected per supernova event, keeping the basic implementation described in
Section 2.8 fixed. The aim of these model variations is to quantify how the
behaviour of stellar feedback affects reionisation and galaxy properties at
high redshift.

####  2.9.1 Bursty supernova feedback

In our first model, labeled bursty-sn due to the bursty star formation
behaviour it drives, energy from all SN events associated to a given stellar
particle is injected simultaneously. Thus, we assume that each stellar
particle hosts one single SN explosion event at 10 Myr, i.e. equivalent to the
mean time at which SN occur for a Chabrier (2003) IMF. We adopt a progenitor
mass of 19.1 M⊙, and assume that each individual SN injects 2 × $10^{51}$ ergs
into its neighbouring cells. Since large amounts of energy are deposited
simultaneously at the same position, supernova feedback in bursty-sn is
particularly efficient. A similar feedback scheme has been previously adopted
by the NewHorizonAGN and Obelisk simulations (Dubois et al., 2021; Trebitsch
et al., 2021).

####  2.9.2 Smooth supernova feedback

The second model we explore is referred to as smooth-sn, because, unlike
bursty-sn, it produces a smooth star formation history. Previous studies (Kimm
et al., 2015; Su et al., 2018; Keller et al., 2022; Keller & Kruijssen, 2022)
have shown that discretising SN injections in time is crucial. Indeed, for a
given stellar population, SN with progenitor masses $>$ 20 M⊙ can explode as
early as 3 Myr, whereas SN with a progenitor mass of 8 M⊙ can explode as late
as 40 Myr. In smooth-sn, the SN delay times are accounted for by randomly
sampling the lifetimes of individual stars using the inverse sampling method
of a polynomial fit for the integrated SN rate from STARBURST99 (see Fig 2. in
Kimm et al. 2015). As before, for each individual SN we assume an energy
injection of $2\times 10^{51}$ ergs and a progenitor mass of 19.1 M⊙. This
setup is a step toward a more physically-motivated feedback model as compared
to bursty-sn, and a variation has been previously adopted by the SPHINX
simulations (Rosdahl et al., 2018, 2022b).

####  2.9.3 Variable supernovae and hypernovae feedback

In this model, called hyper-sn, we vary the injected supernova energy per
explosion event and, in addition, include a contribution from hypernova
events. Indeed, a (metallicity-dependent) fraction of type-II SN is theorized
to be associated to very large energies (Kobayashi et al., 2006), in the range
$10^{51}-10^{53}$ erg. Such events are referred to as hypernovae (HN), and are
expected to be more frequent in low metallicity environments (Kitayama &
Yoshida, 2005; Kobayashi et al., 2006; Kobayashi et al., 2007; Smidt et al.,
2014). Indeed, to explain chemical compositions observed in metal-poor stars,
Grimmett et al. (2020) argue that the hypernova fraction (fraction of all
type-II SN events) in the early universe needs to be $>50\%$, while it
decreases with increasing metallicity to reach the estimated rate of $\sim
1\%$ in the local universe (Podsiadlowski et al., 2004). Simulations of star-
clusters and idealized galaxies have tested the effect of hypernova on small
scales (Su et al., 2018; Brown & Gnedin, 2022), finding that HN feedback can
quench dwarf galaxies for up to $\sim 1$ Gyr, hence making them relevant to
study during reionisation.

In hyper-sn, we assign each stellar particle a metallicity dependant HN
fraction $f_{\rm HN}$ given as

$f_{\rm HN}={\rm max}\biggl{[}0.5\ {\rm
exp}\biggl{(}-\frac{Z_{\ast}}{0.001}\biggl{)},0.01\biggl{]},$ (18)

where $Z_{\ast}$ is the metallicity of the stellar particle. $f_{\rm HN}$
represents the fraction of all SN events a stellar particle undergoes that are
classified as HN. We assume an initial HN fraction of $f_{\rm HN}=50\%$ to
explore the extreme case of strong feedback in metal-free environments. In
case of an HN event, we inject an explosion energy of $10^{52}$ erg. A
constant progenitor mass of 30 M⊙ is assumed for all HN events (Grimmett et
al., 2020).

For SN explosions, we assume that progenitors with different masses (between
8-40 M⊙, see Díaz-Rodríguez et al. 2018) explode at different times and inject
different energies and masses into the surrounding medium. Thus, in addition
to sampling the explosion times as described in the smooth-sn model, we also
stochastically assign different energies to each SN event following a normal
distribution centered at $1.2\times 10^{51}$ erg (best fit distribution using
the Z9.6+W18 model as in Sukhbold et al. 2016), with minimum and maximum
energies of $10^{50}$ and $2\times 10^{51}$ ergs, respectively (Sukhbold et
al., 2016; Díaz-Rodríguez et al., 2018, 2021).

###  2.10 Radiative transfer

RAMSES-RT solves radiation transport by taking the first two moments of the
radiative transfer equation and obtaining a system of conservation laws which
is closed with the M1 closure of the Eddington tensor (Levermore, 1984).
Further details of the methods used for the injection, propagation, and
interaction of the ionising radiation with hydrogen and helium are described
in Rosdahl et al. (2013), while the diffusion of multi-scattering IR radiation
is followed with the trapped/streaming photon scheme presented in Rosdahl &
Teyssier (2015). Each photon group, defined by a frequency interval, is
described in each grid cell by the radiation energy density and the bulk
radiation flux, which corresponds approximately to the radiation intensity
integrated over all solid angles. RAMSES-RT solves the non-equilibrium
evolution of the ionisation fractions of hydrogen and helium, along with
photon fluxes and the gas temperature in each grid cell.

Since radiation is advected on the grid with an explicit solver, the
simulation time-step is subject to the Courant condition. Since the RT time-
step can become significantly smaller than the hydrodynamic time-step, we
subcycle the RT on each AMR level, with a maximum of 500 RT steps performed
after each hydro step (if the projected number of RT steps exceeds 500, the
hydro time-step length is decreased accordingly). During the subcycling on
each level, radiation is prevented from propagating to other levels and the
radiation flux across boundaries is treated with Dirichlet boundary
conditions. Details of the subcycling scheme in the context of flux-limited
diffusion are given in Commerçon et al. (2014).

Figure 1:  The SPICE simulations and their scope and dynamic range. In panel A
we plot the volume-weighted gas density $z=5$, showing the cosmic web on the
largest scales captured by the simulations. Panel B shows a square region
projected in two different quantities: the acceleration on gas produced by
radiation pressure on dust (upper triangle), and infrared flux emergent from
the galaxies (lower triangle). Panel C shows the flux of ionising radiation
escaping from a disc galaxy. We see channels of high (bright orange) LyC
escape within the spiral arms of the galaxy and regions of inefficient escape
(blue). SPICE traces LyC radiation escape on scales of $\Delta x\approx 28\rm
pc$. In panel D we show the surface brightness of the same galaxy in [CII]. In
panel E, we zoom on a region of size $2\times R_{\rm vir}$ centred on a
$\approx 2\times 10^{11}\rm M_{\odot}$ halo. The central galaxy is fed by
large-scale cold inflows. Sub-panels E1,E2,E3 highlight the very different gas
density distributions emerging for this galaxy across our different
simulations. In some simulations (e.g. smooth-sn, hyper-sn), gas has settled
into a disc, while in others (bursty-sn), the gas component is strongly
disrupted by strong outflows. Note that quantities shown in panels A,B,C,D and
E are taken from smooth-sn and can vary between models.

To prevent prohibitively small time-steps and a large number of RT subcycles,
we adopt the ’reduced speed of light approximation’ (RSLA) (Gnedin & Abel,
2001) with the global reduced speed of light set to $\tilde{c}=0.1c$. Note
that such approximation is used to advect the radiation field, while processes
such as radiation pressure are treated with the full speed of light (Rosdahl &
Teyssier, 2015). Studies performed by Costa et al. (2018a, b) demonstrate a
mild effect of the RSLA on the spatial extend of outflows driven by radiation
pressure. When adopted in cosmological simulations targeting reionisation
(Gnedin, 2016; Ocvirk et al., 2016; Deparis et al., 2019), the RSLA is found
to affect mainly the post-reionisation neutral hydrogen fractions, which tend
to be overestimated for a lower speed of light, while the differences are
minimal during the overlap phase. Previous large scale simulations of
reionisation (Rosdahl et al., 2018; Ocvirk et al., 2020; Kannan et al., 2021;
Rosdahl et al., 2022b) adopt different values for $\tilde{c}$, in the range
$(0.01-0.2)c$, while finding convergence with values as low as
$\tilde{c}=0.1c$ (Wu et al., 2019). Therefore, we use results from these
studies and adopt $\tilde{c}=0.1c$.

###  2.11 Stellar radiative feedback

We sample radiation fluxes from stellar populations in five photon frequency
groups, i.e. infrared, optical and three bands of ionising ultraviolet
photons, whose widths are determined by the ionising potentials of HI, HeI and
HeII. The frequency ranges, characteristic energies, ionisation cross-sections
and dust opacities are listed in Table 3. Each of these radiation groups plays
a key role in terms of radiative feedback on galaxy formation. For each
stellar particle, the mass, age and metallicity-dependent stellar specific
luminosities are extracted on-the-fly from the SED model of BPASSv2.2.1
(Eldridge et al., 2017; Stanway & Eldridge, 2018) assuming a Chabrier
(Chabrier, 2003) IMF. The SED spectra (in units of $\rm ergs^{-1}\ \rm
M_{\odot}^{-1}$ Å-1) are tabulated according to age- and metallicity-dependent
luminosities for each radiation group by integrating over the energy
intervals. We use this tabulated SED (single table for each simulation) to
extract the number of photons emitted by each stellar particle. The photons in
each group are directly injected in the host cell of the stellar particle. We
update the photo-ionisation cross sections and energies at fixed intervals for
the radiation groups to represent the luminosity-weighted average of all
emitting sources (see Rosdahl et al. 2013).

The photons in the UV groups not only interact with gas via photo-ionisation
and photo-heating, but also via radiation pressure from photo-ionisation and
dust. In contrast, the infrared and optical groups do not ionise hydrogen or
helium, but interact with gas via radiation pressure on dust. In SPICE we do
not model dust as an active ingredient but assume that the dust number density
scales with the gas-phase metallicity and the gas thermal state. Following
Nickerson et al. (2018), we assume a local dust number density $n_{\rm
d}\equiv(Z/Z_{\odot})f_{d}n_{\rm H}$, where $Z$ is the local metallicity,
$f_{\rm d}=1-x_{\rm HII}$ is the neutral fraction of gas that holds dust, and
$n_{\rm H}$ is the hydrogen number density. This ensures that the dust density
depends not only on enrichment but also on the local thermal state, since dust
is not expected to survive in photo-ionised gas (Finkelman et al., 2012;
Kannan et al., 2020). To each photon group, we assign dust absorption and
scattering opacities as given in Table 3. If absorbed by dust, the photon flux
from the group is added to the infrared group, where it will interact with gas
via multi-scattered radiation pressure (Costa et al., 2018a).

Figure 2: Temperature (left column) and gas-phase metallicity (right column)
of the various simulations at $z=5$. The projections shown cover a region
spanning $L_{\rm box}\times L_{\rm box}/3\times L_{\rm box}$. Consequences of
varying the stellar feedback are imprinted in the thermal state and enrichment
of the large-scale structure. Figure 3: Cosmic star formation rate density as
a function of redshift for the bursty-sn (red line), smooth-sn (blue) and
hyper-sn (turquoise) models. Observational constraints (see text) are marked
with grey circles. hyper-sn most closely reproduces the observed star
formation density at $z\gtrsim 10$, while bursy-sn is closest to observational
constraints at $z\lesssim 8$.

##  3 Results

In this section we present our main results, starting with a qualitative
overview of the simulations.

###  3.1 Overview

Primordial gas cools to form the first stellar particle at $z=22.4$. Feedback
from massive stars via photo-ionisation begins immediately, and supernova
feedback kicks in as early as $3\,\rm Myr$ ($10\ \rm Myr$) after a stellar
particle is born in smooth-sn and hyper-sn (bursty-sn). Galaxies continue to
grow and the first $M_{*}>10^{7}\rm M_{\odot}$ dwarf galaxies form as early as
$z=12$. In Figure 1, we illustrate the dynamic range of the SPICE simulations.
We note that the quantities shown in this figure are largely taken from
smooth-sn, but they can vary dramatically between models.

Panel A shows the dark matter distribution on Mpc scales, illustrating the
cosmic web at $z=5$. At this redshift, SPICE has produced a population of
$\approx 13,000,10,000,7,000$ resolved (see Section 2.5) galaxies with
$M_{*}>10^{5}\rm M_{\odot}$ for the smooth-sn, bursty-sn and hyper-sn models,
respectively. The most massive halo at $z=5$ in each simulation has a virial
mass of $2.7\times 10^{11}\rm M_{\odot}$ and hosts a stellar mass of
$1.9\times 10^{10}\,\rm M_{\odot}$, $2.1\times 10^{9}\,\rm M_{\odot}$ and
$1.7\times 10^{10}\,\rm M_{\odot}$ for smooth-sn, bursty-sn and hyper-sn,
respectively.

In panel B, we zoom onto an over-density in the simulations and highlight the
quantities that describe the effect of novel physics included in SPICE, i.e,
infrared radiation transport and consequently, radiation pressure on dust. The
upper triangle in panel B shows the acceleration imparted to gas via radiation
pressure on dust, with the highest radiation pressure observed very close to
the centers of dusty galaxies. The lower triangle in panel shows the infrared
radial flux emerging from galaxies, which includes intrinsic emission from
stars, reprocessed radiation and dust emission. We see a network of infrared-
bright galaxies with a well established background by $z=5$.

In panel C we zoom further into smaller scales, showing a single galaxy. We
plot the radial flux of hydrogen-ionising radiation ($h\nu>13.6\,\rm eV$)
emerging from a $M_{*}=2\times 10^{10}\rm M_{\odot}$ galaxy. We can see the
escape channels of LyC radiation, with bright orange/yellow colors
representing HII regions through of which ionising radiation escapes most
efficiently, and blue colours represent regions of inefficient escape. We
investigate the connection between galaxy morphology and escape fractions of
galaxies in Section 3.6.

Panel D shows the surface brightness of the same galaxy in [CII] (a line which
efficiently traces the cold phase of the ISM), calculated using a subgrid
model which will be presented in detail in a follow-up study (Bhagwat et al,
in prep). As [CII] is sensitive to the thermal state of the ISM, which is, in
turn, shaped by feedback, this line can be potentially used to constrain
feedback models.

Finally, panel E is a zoom onto a region of size $2\times R_{\rm vir}$ centred
on a $\approx 2\times 10^{11}\rm M_{\odot}$ halo at $z=5$. We can see this
massive halo being fed by infalling gas filaments. We find further galaxies
along the filaments, some with low mass companions. In E1, E2 and E3, we show
the gas distribution of the central galaxies for our three different feedback
models. The bursty-sn (E1) produces a highly irregular and disturbed gas
distribution, in contrast to the smooth-sn (E2) model, which shows a well-
formed disk galaxy. Finally, hyper-sn (E3) produces a compact disk-like
density distribution. These distinct gas morphologies provide us with a first
hint of the widely different effects of our feedback models on the morphology
of the first galaxies.

In the left-hand panel of Figure 2, we show the gas temperature as obtained in
the three feedback models at $z=5$. The projections are calculated in a region
of volume $L_{\rm box}\times L_{\rm box}/3\times L_{\rm box}$. Bright
red/maroon areas trace recent hot ($T>10^{6}\,\rm K$) outflows driven by SN
explosions, while the lighter brown regions ($T\approx 10^{5}$ K) trace older
adiabatically-cooled ionised gas bubbles produced by previous star formation
episodes. While these giant HII regions are present and co-spatial in all
simulations, they are hottest and largest in bursty-sn, and least pronounced
in smooth-sn. Beyond the SN-driven ionised gas lies a volume-filling, diffuse
component with $T\sim 10^{4}\,\rm K$ gas, shown by the bright white regions.
This component consists of photoionised material which has been irradiated by
stellar radiation at distances beyond those reached by feedback-powered
ionised gas. The spatial scale of this diffuse component is set by the
strength of SN feedback, as SN explosions create channels through which
ionising photons can escape into the low density inter-galactic medium.

Figure 4: Star formation rate averaged over 10 Myr as a function of stellar
mass for the bursty-sn (red line), smooth-sn (blue) and hyper-sn (turquoise)
at $z=10$, $7$ and $5$. Solid lines represent the median SFR, with the shaded
regions covering the 16th and 84th percentiles. Observations from various JWST
programs (see text) are shown with grey circles. The dashed line refers to the
best fit relation for the main sequence at $z=6$ (Iyer et al., 2018).

A low volume-filling fraction of $T>10^{4}$ K gas could arise from either
inefficient feedback or a smaller ionising photon budget. However, we show in
Sections 3.3 and 3.6 that the smaller volume-filling fractions seen in smooth-
sn (middle-left panel) and hyper-sn (lower-left panel) arise in spite of a
_higher_ photon budget and are caused by low escape fractions. We see the
presence of cold voids in both smooth-sn and hyper-sn, with the latter being
colder and more extended. The difference between these two models is solely
the SN feedback, highlighting the major role this plays in reionisation.

The right-hand panels of Figure 2 show metallicity projections for the various
feedback models at $z=5$. Simulation bursty-sn produces extended metal-
enriched regions, again illustrating the stronger feedback present in this
simulation. For instance, smooth-sn shows a larger number of compact enriched
gas with typically higher maximum metallicities ($Z>1Z_{\odot}$) at the very
centers of haloes. hyper-sn produces a distribution similar to smooth-sn,
where the enriched regions are compact and with higher maximum metallicities.
However, hyper-sn contains visibly less structure than smooth-sn owing to the
fact that low mass galaxies are much more effectively quenched (see section
3.2 and 3.4).

###  3.2 Cosmic star formation history

Figure 5: Star formation efficiency as a function of halo mass for bursty-sn
(red line), smooth-sn (blue) and hyper-sn (turquoise) at redshift $z=$ 10
(left panel), 7 (middle) and 5 (right). Solid lines represent the median SFE
in each mass bin, with the shaded regions covering the 16th and 84th
percentiles. Grey dots are observations from Stefanon et al. (2021), grey
triangles abundance matching estimates from Tacchella et al. (2018), and the
grey dashed line refers to the best fit from Behroozi et al. (2019).

In Figure 3 we show the time evolution of the star formation rate density
(SFRD) in the various feedback models calculated over the full simulation
volume. We see that the SFRD increases with decreasing redshift as massive
galaxies assemble. The SFRD is consistently highest in the smooth-sn model,
while it is more effectively suppressed in bursty-sn, particularly below
$z\,=\,8$. The SFRD is initially lowest in hyper-sn compared other models,
owing to early strong feedback from HN. However, as the gas metallicity
increases, the HN rate in the massive haloes is reduced exponentially (see
Equation 18), and the weakening feedback results in a rapid increase in star
formation at $z=9-10$. We also notice that the burstiness of star formation in
the different models varies as mentioned in section 2.9. Both bursty-sn and
hyper-sn have strong isolated bursts of star formation followed by relatively
quiescent phases, whereas smooth-sn shows a comparatively smooth star
formation history. The strength of the bursts increases with decreasing
redshift for bursty-sn, while it decreases for hyper-sn and smooth-sn.

By comparing our models to observational constraints from multi-wavelength
studies (Madau & Dickinson, 2014; Rowan-Robinson et al., 2016; Khusanova et
al., 2021; Donnan et al., 2022; Asada & Ohta, 2022; Harikane et al., 2023;
Algera et al., 2023), we see that hyper-sn shows the best agreement at
$z\gtrsim 10$, bursty-sn at $z\lesssim 8$, while smooth-sn is consistently
higher than observations.

Figure 4 shows the SFR of galaxies (averaged over the previous 10 Myr) as a
function of stellar mass for the three feedback models and redshifts. We see
that the median relations from all models roughly lie on a locus, forming a
star formation main sequence (SFMS). Although the median relation does not
vary significantly across models, the scatter differs. We quantify the scatter
by fitting a Gaussian to the distributions of SFRs within bins of stellar
mass. At $z=10$, across the mass range, bursty-sn shows the largest scatter
($\approx 0.8$ dex) followed by hyper-sn ($\approx 0.5$ dex) and finally
smooth-sn ($\approx 0.3$ dex). By $z=7$, the scatter in smooth-sn and hyper-sn
reduces to about $(0.2-0.3)$ dex, in agreement with observations from Speagle
et al. (2014). However, bursty-sn continues to show a large scatter of
$(0.8-0.9)$ dex with it being the largest below $M_{*}<10^{7}\rm M_{\odot}$.
While smooth-sn and hyper-sn show minimal evolution in scatter between $z=7$
and $z=5$, the scatter for bursty-sn drops to $\approx 0.4$ dex.

We compare the SFMS from the three models with with observations from recent
JWST programs (Fujimoto et al., 2023; Haro et al., 2023a; Jung et al., 2023a;
Long et al., 2023; Haro et al., 2023b; Leethochawalit et al., 2023; Robertson
et al., 2023; Heintz et al., 2023a; Jin et al., 2023; Helton et al., 2023;
Atek et al., 2023; Treu et al., 2023; Heintz et al., 2023b; Asada et al.,
2023b; Bouwens et al., 2023b; Looser et al., 2023; Papovich et al., 2023).
Data is collected in bins of $\pm 0.25$ around the redshifts shown. At $z=10$,
the vast majority of observations lie in a mass range not probled well by
SPICE, however, smooth-sn shows the best agreement with observations. At $z=7$
and $z=5$, the median relations from all models is in excellent agreement with
data for galaxies with $M_{*}>10^{7}\rm M_{\odot}$.

Figure 5 shows the ratio between stellar mass and halo mass, which we view as
a proxy for global star formation efficiency (SFE), $M_{*}/(M_{\rm
halo}f_{b})$, as a function of halo mass $M_{\rm halo}$, where $f_{\rm
b}=\Omega_{\rm b}/\Omega_{\rm m}$. We note that all the stellar mass within
the virial radius of a halo is included in the calculation. The smooth-sn
model produces the highest SFE at all times as compared to the other models,
especially at the high-mass end. The bursty-sn model has only a weakly time-
dependent SFE, while the SFE in hyper-sn starts at very low values at early
times, but, following a starburst episode at $z\approx 9$, increases, and by
$z=7$ becomes comparable to the one of the bursty-sn model. The HN rate in
haloes with $M_{\rm halo}>10^{10}\,\rm M_{\odot}$ reaches 1% (by $z\approx
9$), therefore, the energetic component of feedback becomes negligible.
Consequently, the SFE in these haloes increases rapidly due to inefficient
feedback. Lower-mass halos have small SFEs, likely because the high HN rate
combined with the shallower potential wells lead to more effective quenching.
All three models show consistently higher ($\approx 0.6-1$ dex) SFEs as
compared to best fit estimate from Behroozi et al. (2019). They are however in
agreement with abundance matching estimates from Tacchella et al. (2018), and
$z$=7 observations from (Stefanon et al., 2021).

Overall, the results shown in this section highlight how some observables (e.g
SFRD, SFE) vary strongly with implementation of feedback, while others (e.g
SFMS) are less affected. The explosiveness of the bursty-sn model highlights
the need for sustained strong feedback to maintain a sufficiently low star
formation rate. In contrast, the smooth-sn model, characterized by less
energetic feedback, struggles to effectively regulate star formation, thus
allowing for the formation of massive star-forming galaxies by $z=5$. Finally,
the hyper-sn model emphasizes the importance of striking a balance between
bursty and smooth feedback mechanisms.

###  3.3 Reionisation histories

Figure 6: Projections showing HII fraction in the bursty-sn, smooth-sn and
hyper-sn models (from top to bottom) at $z=10$, $7$ and $5$ (from left to
right). Locations of galaxies are over-plotted (see text for selection
criteria), with spiral (star) symbols representing rotationally- (dispersion-)
dominated galaxies (gas morphologies; see section 3.5). Each symbol is color
coded with the LyC escape fraction of the respective galaxy. Already at
$z\,=\,10$ ionised regions show differences in sizes, with bursty-sn and
hyper-sn producing the largest and smallest, respectively. By $z\,=\,5$,
reionisation is complete in bursty-sn, whereas significant neutral patches
persist in the other two models. Figure 7: Volume-weighted neutral hydrogen
fraction as a function of redshift for the three feedback models.
Observational data shown in grey points are taken from Fan et al. (2006);
McGreer et al. (2011); Ono et al. (2012); Schroeder et al. (2012); McGreer et
al. (2015); Greig et al. (2017); Mason et al. (2018, 2019a); Greig et al.
(2019); Hoag et al. (2019); Yang et al. (2020); Wang et al. (2020); Lu et al.
(2020); Jung et al. (2020); Choudhury et al. (2021); Bosman et al. (2022); Zhu
et al. (2022); Gaikwad et al. (2023). bursty-sn reionises within observational
constraints whereas smooth-sn and hyper-sn reionise late.

In Figure 6 we show the ionised fractions ($x_{\rm HII}$) in the three
feedback models representing the progress of reionsation, i.e, formation of
ionised regions at high redshifts, overlap and post-overlap phase. The
difference in their sizes and abundance is already evident at $z=10$, when
bursty-sn produces generally larger HII regions, while the smooth-sn and
hyper-sn models show a combination of extended and localised regions, with
hyper-sn exhibiting a larger number of small ones. The differences are
amplified at $z=7$, where we observe that the bubbles in bursty-sn are already
deep into the overlap phase, with islands of neutral hydrogen localised in
voids with small volume-filling fractions, while both smooth-sn and hyper-sn
still show mostly isolated HII regions with large neutral voids in between. By
$z=5$, the overlap phase is completed in bursty-sn, with neutral hydrogen
present only in dense self-shielded structures, while in both smooth-sn and
hyper-sn small neutral islands still persist.

To demonstrate the varied ability of galaxies to reionise their surroundings,
in Figure 6 we divide each projection in a $4\times 4$ grid and mark the
location of the three most luminous galaxies within each grid sub-division.
Each galaxy is coded by its morphology (stars and spirals represent
dispersion- and rotation-dominated galaxies, respectively; see section 3.5 for
details) and LyC escape fraction (different colors). In Section 3.6, we
discuss in more detail the connection between galaxy morphology and LyC escape
fractions.

In Figure 7 we show the resulting reionisation histories. We find that while
bursty-sn reionises444We consider reionisation to be complete if the volume-
averaged neutral hydrogen fraction of the simulation reaches $10^{-4}$. the
Universe by $z\,=\,5.1$, this is not the case for smooth-sn and hyper-sn
models. In bursty-sn the reionisation process starts early, with the hydrogen
neutral fraction, $\rm x_{HI}$, dropping to $\approx 0.9$ already at $z\approx
10$. From $z\,\approx\,7$, $\rm x_{HI}$ experiences a steep decline in a
series of steps that are time-correlated with three star formation bursts
clearly visible in Figure 3, and resulting in a ‘late reionisation scenario’
(Kulkarni et al., 2019; Keating et al., 2020; Gaikwad et al., 2023), with
$x_{\rm HI}$ becoming lower than $10^{-4}$ at $z\approx 5.1$. Despite
producing more star-forming galaxies (see Section 3.2), smooth-sn and hyper-sn
result in a much later reionisation, failing to bring the volume-weighted
neutral fraction below $30\%$ even by $z\,=\,5$, the time when our simulations
end. However, we note that the simulation volumes are not large enough to
derive a converged mean reionisation history (Iliev et al., 2014; Gnedin &
Madau, 2022), and a box size of $\gtrsim 100$ cMpc/$h$ is required to better
represent an average region of the universe and to account for the typical
sizes of ionised regions in the late stages of reionisation, as these can be
tens of cMpc in size.

The stellar masses integrated over the entire simulation volume at $z=5$ for
the smooth-sn, bursty-sn and hyper-sn model are $\approx 5\times 10^{10}$ M⊙,
$\approx 2\times 10^{11}$ M⊙ and $\approx 1\times 10^{11}$ M⊙, respectively (a
ratio of $1:4:2$). These estimates imply a different number of SN explosions
and intrinsic ionising photon budget, as both increase for larger stellar
masses. Therefore, smooth-sn and hyper-sn have a larger ionising photon budget
as compared to bursty-sn, suggesting that the late reionisation is not a
result of a lower ionising photon budget, but rather it is likely due to
inefficient escape of photons.

Explosiveness of feedback affects star formation which determines the ionising
photon budget and the extent to which gas is disturbed (modulation of escape
fractions). smooth-sn and hyper-sn have a weaker feedback (especially at
$z<8$), which allows for the formation of stable gas configurations (see
Figure 1) and the overproduction of stars (see Figures 3,5 and Section 3.2),
but it is unable to significantly disrupt the gas configuration, hence
suppressing the escape of radiation.

###  3.4 Luminosity functions at 1500 Å

Figure 8: 1500 Å luminosity functions (LF) for the various feedback models at
$z=10$, $7$ and $5$ (from top to bottom). The solid and dashes lines refer to
the intrinsic and dust-attenuated luminosity functions, respectively.
Observations from HST and JWST fields are marked in grey (see text). Stars
represent intrinsic luminosities in bins with fewer than $3$ galaxies. Despite
intrinsic LFs being different, dust-attenuated ones show minimal differences
between models. Therefore, LFs _cannot_ be used to constrain feedback models.

In Figure 8 we show the evolution of the 1500 Å luminosity function, where the
luminosities are calculated in a 10 Å bin around 1500 Å using the stellar
SEDs. Solid and dashed lines refer to the intrinsic and dust-attenuated
luminosity functions, respectively. The latter is calculated with the Monte
Carlo line transfer code RASCAS (Michel-Dansac et al., 2020), by casting 100
rays from each stellar particle within a galaxy to the edge of the halo, and
evaluating the solid angle-averaged attenuation per galaxy. The orientations
of the rays are sampled randomly, and the dust attenuation along each ray is
calculated using the dust model described in Section 2.11. We also show
observational estimates of the luminosity functions from HST legacy fields and
recent JWST programs (Finkelstein et al., 2015; Bouwens et al., 2015; Harikane
et al., 2022; McLeod et al., 2023; Adams et al., 2023; Harikane et al., 2023;
Bouwens et al., 2023b; Bouwens et al., 2023a; Leung et al., 2023).

Figure 9: The columns show JWST RGB composite projections (face-on and edge-
on) in the F200W, F277W and F444W filters for the four most massive galaxies
at $z\,=\,5$ in the bursty-sn (top two rows), smooth-sn (middle two rows), and
hyper-sn (bottom two rows) model. Numbers in each panel indicate the stellar
mass, SFRs and $V/\sigma$ values for the stellar component. The smooth-sn
model produces mainly rotationally supported, massive disk galaxies that are
bulge-heavy; the bursty-sn model dispersion-supported, disturbed systems which
are redder and less-massive; the hyper-sn model a combination of massive
rotatinally-supported spiral galaxies and dispersion-supported ellipsoids.
Stellar feedback processes profoundly affect galaxy kinematics, colours and
morphology.

From the top panel of Figure 8 we see that at $z=10$ hyper-sn produces fewer
galaxies at all magnitudes compared to the other two models. The latter have
similar luminosity functions, with smooth-sn producing the highest number of
(intrinsically) bright objects. Following a decline of the HN rate in massive
haloes in hyper-sn, the increase in star formation boosts the intrinsic LF,
which becomes almost indistinguishable from the LF in the bursty-sn model at
$z=7$. Similarly to $z=10$, smooth-sn shows a larger number (by 0.3-0.4 dex)
of galaxies at all magnitudes compared to the other models.

By $z=5$, smooth-sn produces extremely UV-bright ($M_{1500}\sim-23$) objects
and an abundance of galaxies at all magnitudes. Also hyper-sn results in a
similar abundance of objects at $M_{1500}\sim-23$, while showing less dimmer
galaxies ($-20<M_{1500}<-12)$. Meanwhile, bursty-sn lies in between the other
two models.

At $z=10$ the dust-attenuated LFs are similar to the intrinsic ones for
bursty-sn and hyper-sn, while we observe significant attenuation for smooth-sn
at $M_{1500}<-17$. Dust attenuation becomes more significant at $z=7$ for
smooth-sn and hyper-sn, when an effect is visible at the bright end of the LF.
Indeed, it becomes relevant at $M_{1500}<-16$ for smooth-sn, and at
$M_{1500}<-18$ for hyper-sn, while in bursty-sn dust attenuation is minor even
at the brightest end of the LF. At $z=5$, dust attenuation for smooth-sn is
significant ($>1$ dex) at $M_{1500}<-18$, while in hyper-sn it is present only
at $M_{1500}<-20$. Finally, for bursty-sn it is minimal ($\approx 0.2{\rm
dex}$) for $M_{1500}<-17$. At $z=10$, while bursty-sn and smooth-sn agree very
well with data, hyper-sn shows a deficit of bright galaxies. However, both at
$z=7$ and $z=5$, dust-attenuated LFs from all three models are in excellent
agreement with observations for $M_{1500}<-16$.

Overall, in this section we show that despite the intrinsic LFs being
extremely different, the dust-attenuated ones are very similar below
$M_{1500}<-16$, suggesting that UVLFs cannot be directly used to probe stellar
feedback at high redshifts.

Figure 10: Median $V/\sigma$ values for stars (top row) and gas (bottom row)
as a function of stellar mass for galaxies in the various feedback models at
$z=10,7$ and $5$ (from left to right). The shaded regions show the 16th and
84th percentiles. Percentages in each panel refer to the number of
rotationally-supported galaxies at a given redshift. Large differences emerge
in the kinematics by $z=5$, with smooth-sn, bursty-sn and hyper-sn producing
mainly galaxies which are rotationally-supported, dispersion-supported, and
mixed, respectively.

###  3.5 Galaxy morphologies

Figure 9 shows stellar light projections for the four most massive galaxies
produced in each feedback model, which are matched across simulations using
the IDs of the dark matter particles comprising their parent haloes (which
remain identical across simulations). The images consist of RGB composites in
the F200W(B) + F277W(G) + F444W(R) JWST filters. Each galaxy is shown both
face-on (upper rows) and edge-on (bottom rows), where the total angular
momentum vector was used to define the orientation. Note that dust attenuation
is not taken into account, i.e. the images illustrate intrinsic stellar
emission only.

The most massive galaxies in smooth-sn are blue, spiral galaxies (see Figure
9), which host also bright, red bulges, and have star formation rates reaching
$\approx 50\,\rm\rm M_{\odot}\,yr^{-1}$. Within the same haloes, bursty-sn
generates systems with 10 times lower star formation rates and
$(0.8\\--1)\,\rm dex$ lower stellar mass. The host galaxies now show no
evidence of disks, and appear much more irregular, hosting a number of stellar
clumps and streams. Finally, the hyper-sn model results in a mixture of blue
star-forming spirals with SFRs of $\approx 50\,\rm\rm M_{\odot}\,yr^{-1}$,
similar to those of smooth-sn, and irregular galaxies which also tend to be
star-forming, unlike those in bursty-sn. Variations in supernova feedback have
a particularly striking impact on galaxy morphology, with differences that are
systematic and persist down to stellar masses of $\approx 10^{7}\,\rm\rm
M_{\odot}$ (see Figure 10).

While in Figure 9 we shows the strong qualitative differences caused by
feedback on the four most massive galaxies, we also quantify them for the
entire galaxy population using the ratio $V/\sigma$ (Dubois et al., 2016;
Pillepich et al., 2019) as a proxy for morphology, where $V$ is the rotational
velocity and $\sigma$ is the 3D velocity dispersion. This ratio classifies the
degree to which a galaxy is supported by dispersion or rotation. We calculate
$V/\sigma$ for both stars and gas for each galaxy with $M_{*}>10^{5}$ M⊙. For
stars, we construct a velocity dispersion radial profile, as well as the 3D
rotation curve for all stars within twice the stellar half mass radius of each
galaxy555We bin the rotation and dispersion curves at a fixed bin-width of 0.2
ckpc ($\sim$8 times the maximum resolution of the simulation) to ensure
statistical mean/medians in each bin.. As for gas using all cells inside a
fixed radius can lead to noisy calculations due to inflows and outflows, we
select gas cells within $R_{\rm vir}$ adopting an H$\alpha$ emissivity
threshold. Indeed, recombination lines of atomic hydrogen from the ionised ISM
are often used to measure gas kinematics in a wide redshift range (de Graaff
et al., 2023). We evaluate the case-B recombination surface emissivity for
each gas cell as (Hummer & Storey, 1987):

$\epsilon_{{\rm H}\alpha}=h\nu\ P_{\rm B}(T)\ \alpha_{B}(T)n_{\rm e}n_{\rm
p},$ (19)

where $\alpha_{\rm B}=2.753\times 10^{-14}\
\lambda^{1.5}/(1.0+(\lambda/2.74)^{0.407})^{2.242}$ with $\lambda=315614/T$ is
the case-B recombination coefficient (Hui & Gnedin, 1997), $T$ is the gas
temperature, $n_{\rm e}$ and $n_{\rm p}$ are the number densities of electrons
and protons, and $P_{\rm B}(T)$ is the conversion probability per
recombination event, which is $\approx 0.451$ for $T=10^{4}$ K.

To calculate gas kinematics we select gas cells with $\epsilon_{{\rm
H}\alpha}>\rm 10^{-6}\ erg\ s^{-1}\ cm^{-2}$, noting that this choice does not
affect our qualitative results666This limit is taken using results on
completeness estimates and theoretical studies from Belfiore et al. (2022) and
Tacchella et al. (2022). Following the process used in various theoretical and
observational studies (Pillepich et al., 2019; Rizzo et al., 2020), $V/\sigma$
is calculated using $V$ at the peak of the rotation curve and $\sigma$ as the
mean velocity dispersion. We classify galaxies with $V/\sigma\geq 1$ (for both
stars and gas) as rotationally-supported, otherwise as dispersion-supported.

In Figure 10 we show the median $V/\sigma$ for stars (top row) and gas
(bottom) as a function of stellar mass. The numbers in the various panels give
the fraction of all galaxies that are rotationally-supported systems at a
given redshift. The stellar component in bursty-sn and hyper-sn shows a
consistent decrease with decreasing redshift in the fraction of rotationally-
supported galaxies (hyper-sn exhibits a sudden drop from $\approx 18\%$ at
$z=7$ to $\approx 6\%$ $z=5$), whereas the fraction does not change
significantly in smooth-sn. The gas component behaves differently for all
models. bursty-sn shows a steady increase in the fraction of rotationally-
supported galaxies with decreasing redshift (from $\approx 31\%$ at $z=10$ to
$\approx 42\%$ at $z=5$), smooth-sn has a similar trend but the increase is
more drastic (from $\approx 27\%$ to $\approx 53\%$), while hyper-sn shows a
mild increase in this fraction between $z=10$ and $7$, followed by a drop in
the range $z=7-5$.

At $z=10$, both stellar and gas components are indistinguishable between
models. Similarly at $z=7$, with the only differences appearing at the highest
masses, where smooth-sn and hyper-sn show formation of highly rotationally-
supported systems. The differences between the various models become most
pronounced at $z=5$ for both stellar and gas components. smooth-sn shows a
population of highly rotationally-supported galaxies ($V/\sigma>2$) and very
large scatter at $M_{*}\gtrsim 10^{8}\rm M_{\odot}$, while bursty-sn produces
a dispersion-supported galaxy population with a small scatter around the
median. The hyper-sn model exhibits a transition at $M_{*}\approx 10^{8}\rm
M_{\odot}$: scatter is small and galaxies are mainly dispersion-supported for
masses lower than this value, while more massive galaxies are typically
rotationally-supported and scatter is more significant.

Figure 11: Median escape fraction of SPICE galaxies as a function of stellar
mass. Solid and dashed lines refers to dispersion- and rotation-dominated
systems ($V_{\rm gas}/\sigma_{\rm gas}$), respectively. Solid (hollow) stars
represent single dispersion-(rotation-) dominated galaxies in bins with fewer
than $3$ galaxies.

We note that the gaseous component shows a larger degree of rotational support
as compared to the stellar component (i.e. $V_{\rm star}/\sigma_{\rm
star}<V_{\rm gas}/\sigma_{\rm gas}$), a trend particularly prominent in the
smooth-sn and hyper-sn models. The number of galaxies classified as
rotationally- vs dispersion-supported is also widely different for the three
models, especially at $z=5$, implying the presence of kinematic misalignment
between the stellar and gas components in the galaxies.

Examining the galaxy morphologies, we show that feedback plays a key role in
shaping the morphological mix of galaxies that emerge post-reionisation. The
populations are hard to distinguish using $V/\sigma$ as an indicator at $z>7$,
except for the most massive galaxies. Below this redshift however, the
morphological mixes diverge and hence can act as an indicator to distinguish
feedback models.

###  3.6 Implications for LyC escape : What galaxies drive reionisation?

LyC escape can proceed either through low density channels created by stellar
feedback, or if the ISM is highly ionised and optically thin (Zackrisson et
al., 2013; Katz et al., 2022). These scenarios are not mutually exclusive and
previous studies (Trebitsch et al., 2017; Kimm et al., 2017; Rosdahl et al.,
2018; Barrow et al., 2020; Katz et al., 2022; Yeh et al., 2023) show that LyC
escape fraction ($f_{\rm esc}$) is strongly regulated by feedback through a
complex multiphase ISM. Therefore, the morphology of the ISM holds signatures
that can aid in the identification of the galaxies responsible for
reionisation.

Here we connect the morphologies of galaxies as discussed in the previous
section to their $f_{\rm esc}$, which is calculated using RASCAS (Michel-
Dansac et al., 2020). Photon packets are injected at the position of each
stellar particle with a probability proportional to its LyC luminosity,
evaluated using the stellar SEDs. We set the number of photon packets per halo
to be 100 times the number of star particles, with a maximum of $10^{7}$.
Photons are propagated until they are absorbed or reach a distance equal to
the virial radius of the host halo. The escape fraction is then defined as the
ratio between the number of escaping and injected photons.

We assign a LyC $f_{\rm esc}$ to each galaxy, and use $V_{\rm gas}/\sigma_{\rm
gas}$ (as calculated in section 3.5) to sort galaxies into rotationally- and
dispersion-supported categories. In Figure 11 we show the median LyC escape
fractions as a function of stellar mass for dispersion- (solid lines) and
rotationally-supported (dashed) galaxies. We observe a clear difference in the
$f_{\rm esc}$ of rotation- and dispersion-dominated systems, with the latter
exhibiting $\approx 10-50\times$ higher escape fractions as compared to their
rotational counterparts. This trend holds for all three feedback models at all
redshifts in consideration. The differences in $f_{\rm esc}$ imply that at the
same intrinsic luminosity and in the same environment, dispersion-dominated
systems produce ionised bubbles larger than those of their rotation-dominated
counterparts. This effect is observed in Figure 6 and is particularly evident
at $z=10$ (left-column), where ionised regions surrounding dispersion-
dominated galaxies (stars) are the largest (e.g in bursty-sn), while those
around rotation-dominated galaxies (spirals) are the smallest (e.g in smooth-
sn). The relative difference between the $f_{\rm esc}$ for rotation- vs
dispersion- dominated systems is the largest for smooth-sn, followed by
bursty-sn and hyper-sn. However, the morphological mix in the three models is
redshift and stellar mass dependent. Since dispersion-dominated systems are
more prevalent if stellar feedback is more explosive, reionisation occurs
earlier in bursty-sn as compared to the other two models, where rotation-
dominated systems are preferentially produced, especially at the highest
masses.

Figure 12: Top: Intrinsic (solid lines) and escaping (dashed line) ionising
emissivity as a function of stellar mass for bursty-sn (red), smooth-sn (blue)
and hyper-sn (turquoise). Bottom: Cumulative distribution of the above
quantities indicating net contributions from galaxies in different mass bins.
Emissivities are calculated over the full duration of the simulation.

We next quantify differences in reionisation histories by looking at the net
intrinsic and escaping emissivities from all galaxies, defined as $f_{\rm
esc}\dot{N}_{\rm ion}$, where $\dot{N}_{\rm ion}$ is the production rate of
photons with $E_{\gamma}>13.6$eV calculated by integrating the stellar SED,
and $f_{\rm esc}=1$ for the instrinsic emissivity. The top panel of Figure 12
shows the total intrinsic (solid lines) and escaping (dashed lines)
emissivities of the galaxies over the entire duration of the simulation as a
function of stellar mass. As previously discussed, the net intrinsic
emissivity is the highest for smooth-sn, followed by hyper-sn and bursty-sn.
The intrinsic contributions for smooth-sn and hyper-sn are largely dominated
by the most massive galaxies (see also discussion below), whereas bursty-sn
shows a more even distribution across the mass range. However, the net
escaping emissivity is the highest for bursty-sn, followed by smooth-sn and
finally hyper-sn. We can also use the total intrinsic and escaping emissivity
to calculate a global escape fraction for each simulation, finding $\approx
6.8,1.8$ and $1.5\%$ for bursty-sn, smooth-sn and hyper-sn, respectively.

To understand which galaxies contribute the largest to reionisation, we look
at the cumulative distribution of emissivities, which are shown in the bottom
panel of Figure 12 as intrinsic (solid lines) and escaping (dashed lines).
bursty-sn results in a fairly even contribution from galaxies with
$M_{*}>10^{7}\rm M_{\odot}$ to both intrinsic and escaping emissivities.
Differently, the intrinsic emissivity in smooth-sn and hyper-sn is largely
dominated by galaxies with $M_{*}\gtrapprox 10^{9}\rm M_{\odot}$ ($\approx
50\%$ and $\approx 65\%$, respectively). While the escaping emissivity shows a
distribution similar to the intrinsic one (i.e. dominated by galaxies with
$M_{*}\gtrapprox 10^{9}\rm M_{\odot}$), hyper-sn exhibits an even contribution
from galaxies with $M_{*}>10^{7}\rm M_{\odot}$. Therefore, in bursty-sn and
hyper-sn, galaxies with $M_{*}>10^{7}\rm M_{\odot}$ contribute evenly to
reionisation, whereas in smooth-sn $\approx 46\%$ of the escaping emissivity
is produced by $M_{*}>10^{9}\rm M_{\odot}$ galaxies.

Figure 13: The median ionising photon production efficiency $\xi_{\rm ion,0}$
as a function of stellar mass. Scatter around the median represents 16th and
84th percentiles. $\xi_{\rm ion,0}$ is calculated over the full duration of
the simulations. All models show very similar $\xi_{\rm ion,0}$, although
leading to very different reionisation histories suggesting that reionisation
is $f_{\rm esc}$ limited.

As measuring $\dot{N}_{\rm ion}$ at high redshifts is challenging, other
quantities are used to evaluate the ionising photon production, such as the
ionising photon production efficiency $\xi_{\rm ion,0}$ = $\dot{N}_{\rm
ion}/L_{1500}$, where $L_{1500}$ is the intrinsic luminosity of the galaxy at
1500 Å. $\xi_{\rm ion,0}$ can give key insights about the stellar populations
and dust obscuration and is used to study escape fractions of high redshift
galaxies (Simmonds et al., 2023b). In Figure 13, we show $\xi_{\rm ion,0}$ as
a function of stellar mass of galaxies calculated over the full duration of
the simulation. All three models show comparable median trends with minor
differences of $\approx(0.1-0.3)$ dex, however, bursty-sn has a very large
scatter at all masses. Previous estimates from the HST surveys (Robertson et
al., 2013) and more recent JWST programs (Atek et al., 2023; Simmonds et al.,
2023b, a; Endsley et al., 2023a; Endsley et al., 2023b; Saxena et al., 2023a)
find median values of log($\xi_{\rm ion,0}$) $\approx 24.8-25.7$ for galaxies
in the range $z\sim 5-11$, which is in agreement with all three models.

Despite a very similar ionising photon production efficiency, the three models
result in quite different reionisation histories (see section 3.3), suggesting
that reionisation is strongly $f_{\rm esc}$ limited. Indeed, we have seen that
different feedback models produce a different morphological mix (bottom panel
in Figure 10), which in turn affects radiation escape depending on the
relative predominance of rotation- or dispersion-dominated systems, and hence
the escaping ionising emissivity budget. Therefore, models such as smooth-sn
which preferentially produce rotationally-dominated systems (especially the
most massive and luminous galaxies) show very low global escape fractions and
produce very late reionisation histories. Meanwhile, bursty-sn produces a
majority of dispersion-dominated systems (especially the most massive and
luminous galaxies) which implies a high global $f_{\rm esc}$ ($\approx
4\times$ higher than smooth-sn) and reionises the simulation volume by
$z=5.1$.

##  4 Discussion

Here we discuss the implications of feedback variations on galaxy properties,
the prospects of constraining stellar feedback using observations, and the
limitations of our work.

###  4.1 Connecting galaxy morphology and LyC escape: insights for
observations

We have shown that the SN feedback mode hardly affects some observables (e.g
SFMS, UVLF), while strongly altering others (e.g SFRD, SFE, reionisation,
morphologies, $f_{\rm esc}$). A strong indicator to test feedback models we
point to is the morphological mix that emerges post-reionisation, as feedback
alters in a fundamental way stellar and gas morphologies. The gas morphologies
in turn affect the $f_{\rm esc}$ of galaxies. As galactic morphology can be
probed using different tracers (including stellar light and emission line
kinematics), a multi-wavelength study of morphological characteristics
provides a key insight to constrain feedback models (see Figure 6).

Telescopes like JWST and ALMA allow us to map galaxies at high angular
resolutions ($\approx 0.1-0.2$ arcsec) deep into the epoch of reionisation.
Studies using JWST observations have already started to produce galactic
morphology statistics at various redshifts using e.g. stellar surface
brightness profiles, Sersic indices, and Gini indices (Huertas-Company et al.,
2023; Treu et al., 2023; Jacobs et al., 2023; Kartaltepe et al., 2023; Vega-
Ferrero et al., 2023; Tacchella et al., 2023; Sun et al., 2023c). All find a
rich morphological diversity of galaxies well established already at $z>5$,
and estimate a disk (pure disk + disks with bulges) fraction in the range
$\sim 30-40\%$ at $z\sim 6$, in agreement with our busty-sn model. In this
work we just consider dispersion- or rotation-dominated galaxies, finding that
their relative distribution varies strongly depending on the feedback model.
We defer to a future work a more direct comparison to observations, which will
help to better constrain our models.

Recent studies have used SED fitting (Looser et al., 2023; Dressler et al.,
2023a), [OIII]+H$\beta$ equivalent width and H$\alpha$ emission line
luminosities (Endsley et al., 2023a; Endsley et al., 2023b; Simmonds et al.,
2023a) to characterise modes of star formation in JWST galaxies, finding that
a bursty SFH is prevalent in galaxies with masses $M_{*}\sim 10^{7-10}\rm
M_{\odot}$ in the range $z=6-12$. These studies also show that galaxies with a
bursty SFH are likely to drive reionisation. In sections 3.3 and 3.5 we have
discussed the strong influence that a bursty SFH has on galaxy morphology and
on the reionisation history, confirming the claim that galaxies with bursty
SFH are the likely drivers of reionisation. While observations have not yet
connected galaxy morphologies to their ionising efficiencies, in SPICE we find
that galaxies with a bursty SFH are largely dispersion-dominated (especially
evident at $z=5$), while galaxies with a non-bursty or smooth SFH are mainly
rotation-supported and are unlikely to drive reionisation (see section 3.6).
Such prediction can be tested with JWST observations.

The impact of cold streams (Dekel et al., 2009; Bournaud & Elmegreen, 2009; Oh
et al., 2018), gas fractions (Barnes & Hernquist, 1996; Barnes, 2001; Naab et
al., 2006), mergers (Toomre, 1977; Hopkins et al., 2009; Kannan et al., 2015;
Velázquez et al., 2020), tidal effects (Bekki, 1998) and stellar feedback
(Okamoto et al., 2005; Agertz & Kravtsov, 2016) on galaxy morphologies has
been a topic of debate over the last few decades. While it is widely accepted
that disk galaxies form due to gas accretion and transform to spheroidal
galaxies via mergers (Toomre & Toomre, 1972; Toomre, 1977; Naab et al., 2006),
the picture becomes unclear at high redshifts. For instance, the short
dynamical timescales at high redshifts along with the long depletion
timescales can allow for disk galaxies to stabilize even after major mergers
(Robertson et al., 2006). As in all SPICE simulations the initial conditions
are the same, we expect the timing of halo mergers as well as the cold streams
feeding the haloes to be very similar across all models. Therefore, we argue
that the stark differences in galaxy properties are a consequence of supernova
feedback systematically altering the structure of galaxies on ISM/CGM scales.
These differences translate not only to galaxy morphology (see section 3.5),
but also to the LyC escape fractions (see section 3.6).

As a direct measurement of LyC escape fractions from high redshift galaxies is
not possible because of the high IGM opacity, alternative indirect diagnostics
have been suggested, such as metal line ratios (Wang et al., 2019; Chisholm et
al., 2020; Saxena et al., 2022; Schaerer et al., 2022), Ly$\alpha$ peak
separations (Verhamme et al., 2015; Verhamme et al., 2017) and recent SFR
(Calzetti, 2012; Kennicutt & Evans, 2012; Velázquez et al., 2020). While all
these probes are sensitive to dust attenuation, line strengths and orientation
effects, galaxy morphology can be reliably traced with high resolution
observations, and could thus be used as a better proxy for the escape
fraction. Indeed, here we show that a strong correlation exists between the
morphology of galaxies (described in terms of $V/\sigma$) and their escaped
radiation, with dispersion-dominated systems exhibiting the highest LyC escape
fractions at all redshifts (by factor of 10-50). Our results also indicate
that, as reionisation is limited by the escape fraction, it is strongly
dependent on the behavior of stellar feedback, as this needs to be explosive
and bursty (see section 3.3, 3.6).

Due to the differences in LyC $f_{\rm esc}$, at the same intrinsic luminosity,
dispersion-dominated galaxies should preferentially create large expanding HII
regions as compared to their rotation-dominated counterparts (see Figure 6).
The growth of the HII regions is thus connected to the galaxy morphology,
suggesting that concurrent observations of galaxy morphology and ionised
regions could help to establish more firmly the connection between morphology
and $f_{\rm esc}$, and to constrain stellar feedback. Observations of
Ly$\alpha$ emitters (LAE) from HST/Keck-MOSFIRE (Stark et al., 2016; Matthee
et al., 2018) suggest that bright LAEs are preferentially surrounded by large
ionised regions, and also characterised by intense star formation. However,
recent JWST observations of LAEs at $z\sim 7-9$ (Endsley et al., 2023a;
Whitler et al., 2023; Tang et al., 2023; Jung et al., 2023b; Saxena et al.,
2023b) suggest a more diverse HII regions distribution, where not all strong
LAEs lie within the largest ionised regions. Therefore, deep spectroscopic
observations combined with imaging of LAEs will help to better constrain
bubble sizes along with galactic properties. In a companion study (Bhagwat et
al. in prep.) we will investigate Ly$\alpha$ characteristics of galaxies in
the three feedback models.

Understanding the demographics of outflows driven by SN feedback and radiation
pressure (Hayward & Hopkins, 2016; Li et al., 2017; Costa et al., 2018b; Menon
et al., 2023) is key to constrain feedback models at high redshifts. Recent
JWST observations (Zhang et al., 2023; Carniani et al., 2023) investigate the
incidence of ionised outflows using H$\alpha$ and [O III] in low mass galaxies
($M_{*}<10^{10}\rm M_{\odot}$) in the range $z\sim 6-9$. These studies find
that the inferred outflow velocities and mass loading factors are,
respectively, 3 and 100 times larger than those observed for local dwarfs
(Marasco et al., 2023). Multi-wavelength investigations combining JWST and
ALMA data ([O III] and [C II], respectively; e.g. Fujimoto et al. 2022)
suggest that the incidence of an outflow from a $M_{*}\approx 10^{8}\rm
M_{\odot}$ galaxy at $z\sim 8.5$ is likely associated to a starburst. This
would then also promote a strong ionising photon escape. In SPICE we have
shown how feedback affects the outflow characteristics (see Figure 2)
suggesting that outflow properties and incidence rates can be a key tool to
constrain feedback models.

###  4.2 Numerical uncertainties and physical limitations

The range of physics included in SPICE, along with a physical resolution of
$\sim 28$ pc, allows to make robust predictions for a variety of galaxy
observables. However, due to approximations in numerical schemes, subgrid
models and stochasticity, SPICE is affected by uncertainties. Understanding
them is critical to improve models for future simulations, as well as the
interpretation of their predictions.

While SPICE models metal line cooling, we do not include a molecular line
cooling channel, which is relevant in the dense gas phase of the ISM. This can
potentially impact, among others, the presence of cold gas in the outflows
driven by SN feedback (Richings et al., 2014a, b; Biernacki & Teyssier, 2018),
as well as the distribution and efficiency of star formation. While the
refinement strategy adopted in SPICE allows us to resolve the CGM of galaxies
at $\sim(28-84)\rm pc$ scales, the resolution deteriorates further away from
the centers of haloes, with possible consequences on the correct modeling of
the energy, mass loading and multiphase nature of outflows (Rey et al., 2023).
This in turn can have important effects on the modeling of the escape of LyC
radiation.

We attempted to improve the feedback modelling by making it increasingly
physically based. bursty-sn represents the "IMF averaged" model, whereas
hyper-sn includes realistic SN energies and explosion times, along with a
theoretically supported prescription for hypernovae (Sukhbold et al., 2016;
Kobayashi et al., 2006; Grimmett et al., 2020). However, we note that the
latter, which is the most physically motivated model, is unable to complete
reionisation by the end of the simulation, whereas bursty-sn (the least
physically motivated model) is more consistent with observational constraints.
While this could point to issues in the assumed HN rates or SN energy
distributions, it could also indicate a potentially missing ingredient in the
feedback modelling, or the need for an improved numerical scheme for injection
of energy and momentum.

The ISM model used in SPICE is a combination of subgrid prescriptions based on
idealised simulations (e.g turbulence from Federrath 2016 and mechanical SN
feedback from Kim & Ostriker 2015b). These prescriptions include a number of
correction terms that are not predicted ab initio, but rather derived from
smaller-scale simulations. Therefore, central assumptions such as log-normal
density PDF and missing feedback processes at sub-resolution scale are a
source of uncertainties.

Detailed studies of population synthesis models and SEDs have shown that the
SED choice can strongly affect LyC escape fractions and reionisation histories
(Ma et al., 2016; Rosdahl et al., 2018; Ma et al., 2022; Rosdahl et al.,
2022a). Further observational evidence (Götberg et al., 2019, 2020; Götberg et
al., 2023) suggests that e.g. the contribution from stripped binary stars
could boost the ionising photon budget from massive stars. Hence, the SED and
the ingredients of populations synthesis models are important uncertainties
which should be further explored in future work.

SPICE employs the reduced speed of light approximation to advect radiation.
This approximation produces converged reionisation histories as long as the
speed of the ionisation fronts (I-fronts) remains below the value of the
adopted reduced speed of light, which in SPICE is $0.1c$. We note that at the
tail end of reionisation, the I-fronts can travel as fast as $0.1c$ as they
traverse the voids (D’Aloisio et al., 2019), but the approximation is
unavoidable because of the extreme computational costs of running simulations
with the full speed of light.

We also note that the simulation volume is too small to study the global
reionisation on large scales. Indeed, typical sizes of ionised regions can
become tens of cMpc, suggesting that to produce converged reionisation
histories boxes $\gtrsim 100$cMpc/$h$ are required. In follow-up studies we
plan on extending the SPICE models to larger simulation volumes to get better
galaxy statistics and to further improve numerical schemes of stellar
feedback.

SPICE addresses the effect of variations in stellar feedback models while
including some poorly-explored physical processes, like radiation pressure on
dust and dust itself. However, due to the computational costs involved, our
models are far from complete, as processes such as cosmic rays (CRs), AGN
feedback and magnetic fields are still missing. Magnetic fields and CRs can
both affect the structure of the ISM which in turn affects galaxy obervables.
Numerous studies have attempted to quantify the effect of magnetic fields
(Gnedin et al., 2000; Marinacci & Vogelsberger, 2015; Krumholz & Federrath,
2019; McKee et al., 2020; Garaldi et al., 2021; Katz et al., 2021; Martin-
Alvarez et al., 2020) and CRs (Sazonov & Sunyaev, 2015; Leite et al., 2017;
Farcy et al., 2022; Martin-Alvarez et al., 2022) on galaxy formation and
reionisation itself. However, the implications remain unclear.

Dust is relevant for various processes such as cooling, fragmentation,
absorption and scattering of radiation and obscuration of UV photons. SPICE
adopts a simplified dust model in which dust is hydrodynamically coupled to
the gas and its creation and destruction are not accounted for. Additionally,
we assume that dust opacities are independent of temperature, although in
dense cold gas, they can scale as $T^{2}$ (Semenov et al., 2003). Scaling of
opacities can also strongly affect the radiation pressure on dust, which we
include as a feedback channel. Therefore, on-the-fly modelling of dust grain
creation and destruction, as well as more accurate dust opacities, are
important for a better modelling of the ISM structure and of the strength of
radiation pressure driven feedback.

Theoretical studies have argued that AGN feedback can be important to
understand the evolution of dwarf galaxies, even at very high redshifts
(Dashyan et al., 2017; Koudmani et al., 2021, 2022; Sharma et al., 2022). At
the same time, observational evidence is mounting supporting the presence of a
population of $\sim(10^{5}-10^{8})\ \rm M_{\odot}$ black holes at very high
redshifts (Scholtz et al., 2023; Mezcua et al., 2023; Schneider et al., 2023;
Maiolino et al., 2023). Therefore, AGN feedback is a key missing ingredient
worth exploring in the future to investigate when and in which galaxies it is
relevant.

##  5 Conclusions

In this work we introduce SPICE, a suite of radiation hydrodynamical
simulations of galaxy formation and reionisation performed with RAMSES-RT. Our
aim is to systematically study the effects of variations in stellar feedback
on reionisation and properties of the first galaxies. SPICE resolves atomic
cooling haloes and has a spatial resolution of $\approx 28$ pc at $z=5$. The
galaxy formation model includes cooling, non-equilibrium chemistry, a multi-
freefall star formation model, which employs a subgrid prescription for
turbulence to evaluate the star formation efficiency in each computational
cell, and a mechanical feedback scheme to inject energy and mass from
supernovae. We include radiation transport (through gas and dust) in 5
frequency bands, i.e. IR, optical and three UV groups. Our base model remains
identical across our different simulations with the sole exception of SN
feedback, which is modelled in three different way (see Table 2): bursty
(bursty-sn), smooth (smooth-sn) and smooth feedback with a time-evolving
energetic component (hyper-sn).

In this introductory paper, we showcase the global properties of the
simulations in terms of galaxy populations, reionisation histories, and
connection between galaxy morphologies and LyC escape fractions. A summary of
our key findings is presented below:

  * •

Feedback strongly affects the burstiness and amplitude of the star formation
rate density (SFRD; section 3.2). bursty-sn consistently exhibits bursts of
SF, with a more pronounced intensity below $z=8$. Conversely, smooth-sn has a
minimal burstiness, consistently reaching the highest SFRD among the three
models. Finally, in hyper-sn, the SFRD is the lowest at the highest redshifts
because of the strong feedback from hypernovae. This also induces pronounced
burstiness at $z>9$, akin to bursty-sn. As the hypernova rate declines (due to
metal enrichment), below $z\sim 9$ the model exhibits similarities with
smooth-sn and shows a strong upswing in SFRD.

  * •

We observe the emergence of a star formation main sequence (SFMS; section 3.2)
already by $z=10$. The median SFMS remains very similar across models,
however, bursty-sn exhibits a much larger scatter ($0.3-0.5$ dex higher) as
compared to the other two models. All three models show good agreement with
observations from various JWST programs.

  * •

Feedback strongly affects reionisation histories (section 3.3). bursty-sn
completes reionisation by $z=5.1$, consistent with observational constraints.
This is not the case for smooth-sn and hyper-sn, despite both models yielding
an excess of ionising photons in comparison to bursty-sn. We thus confirm that
reionisation is very sensitive to the modulation of $f_{\rm esc}$ due to
feedback, rather than being driven solely by the number of available ionising
photons.

  * •

The impact of feedback is visible also in the intrinsic $1500$ Å luminosity
function (LF; section 3.4). While in smooth-sn and hyper-sn we observe an
abundance of highly UV-bright galaxies ($M_{1500}<-23$), these are absent in
bursty-sn. The hyper-sn (smooth-sn) model consistently yields the lowest
(highest) galaxy counts at all redshifts.

  * •

The evolution of the dust-attenuated LFs is different from that of the
intrinsic ones (section 3.4). Interestingly, while only in smooth-sn we
observe signs of dust attenuation at $z=10$, by $z=7$ all models exhibit
noticeable levels of attenuation, which becomes even more pronounced at $z=5$.
All three models show an excellent agreement with observations, especially at
$z=7$ and 5. Despite intrinsic LFs being different, dust-attenuated ones show
minimal discrepancies across models. Therefore, LFs cannot be directly used to
constrain feedback models.

  * •

We see striking differences in post-reionisation galaxy morphologies,
characterized by the ratio $V/\sigma$ (section 3.5), where we classify
galaxies with $V/\sigma\geq 1$ as rotationally-supported and dispersion-
supported otherwise. bursty-sn exhibits a preference for dispersion-dominated
systems, primarily at the higher mass end, while in smooth-sn most galaxies
are rotationally-dominated at any mass. Finally, in hyper-sn galaxies with
masses below (above) $M_{*}\approx 10^{8}\rm M_{\odot}$ typically show strong
dispersion (rotational) support.

  * •

The differences in morphology translate into variations in Lyman continuum
(LyC) escape fractions ($f_{\rm esc}$; section 3.6 and Figure 11). We find
that $f_{\rm esc}$ of dispersion-dominated systems is enhanced by a factor of
$\approx 10-50$ in comparison to that of rotation-dominated ones at all
redshifts.

  * •

Galaxies that undergo strong bursts of feedback preferably produce dispersion-
dominated systems (especially at the highest mass and luminosities), and hence
higher $f_{\rm esc}$ as compared to galaxies that have smoother feedback,
which allow for the formation of stable, rotationally-supported systems (in
particular at the highest mass and luminosities). Therefore, feedback
determines the global morphological mix of galaxies, and, as a consequence,
the global $f_{\rm esc}$. The connection between morphologies and $f_{\rm
esc}$, therefore, can potentially be an excellent probe not only of feedback
models at high redshifts, but also of the galaxies that drive reionisation.

  * •

The intrinsic contribution to ionising emissivity ($\dot{N}_{\rm ion}$; solid
lines in Figure 12) is dominated by galaxies with $M_{*}\gtrsim 10^{9}\rm
M_{\odot}$ in hyper-sn ($\approx 65\%$) and smooth-sn ($\approx 50\%$),
whereas bursty-sn shows a similar contribution from galaxies of all stellar
masses above $\sim 10^{7}\rm M_{\odot}$.

  * •

The escaping emissivity ($f_{\rm esc}\dot{N}_{\rm ion}$, dashed lines Figure
12) in the three models is dominated by galaxies in different mass ranges. In
smooth-sn the largest contribution ($\approx 46\%$) comes from
$M_{*}>10^{9}\rm M_{\odot}$ galaxies, whereas bursty-sn and hyper-sn show
roughly uniform contribution from galaxies in the range $M_{*}\approx
10^{7-10}\rm M_{\odot}$. Upon comparing the escaping and intrinsic
emissivities, we find a global (i.e. from the entire galaxy populations at all
redshifts) escape fraction of $\approx 6.8,1.8$ and $1.5\%$ for bursty-sn,
smooth-sn and hyper-sn, respectively.

  * •

The ionising photon production efficiency ($\xi_{\rm ion,0}=\dot{N}_{\rm
ion}/L_{\rm 1500}$) evaluated over the full duration of the simulation shows
comparable median trends (see Figure 13) for the three models, with
differences of $\approx 0.1-0.3$ dex. Despite a very similar ionising photon
production efficiency, the three models result in quite different reionisation
histories, suggesting that reionisation is strongly $f_{\rm esc}$ limited.

The feedback variations as well as the high resolution of SPICE allow us to
probe in detail the first sources of reionisation in their stellar, gas and
radiative components. In this work, we focused on the morphology and LyC
$f_{\rm esc}$ to connect feedback to reionisation using kinematics derived
from stars and H$\alpha$ bright gas. However, imprints of feedback will
manifest in a wide range of observables. The multi-phase ISM model, along with
the additional radiation transport in IR and optical, allows us to investigate
a variety of emission lines (such as H$\alpha$, [C II], [O III] and
Ly$\alpha$). Indeed, all these probes offer a unique opportunity to make
predictions for and comparisons to observations across a wide range of
wavelengths, from state-of-the-art facilities such as JWST, ALMA and MUSE, as
well as from planned ones such as SKA and ELT. Therefore, a comparison between
synthetic SPICE observables and real data will help to disentangle and
understand different feedback mechanisms at high redshifts.

## Acknowledgements

Computations were performed on the HPC systems Raven and Cobra at the Max
Planck Computing and Data Facility (MPCDF). AB would like to thank the MPCDF
team for their support with HPC tickets. The authors would like to extend a
great deal of gratitude to Joakim Rosdahl for helping with RAMSES-RT setups
and sharing analysis pipelines employed in this work. The authors are grateful
to Taysun Kimm for developing and sharing the mechanical feedback model
adopted and modified in this work. AB would also like to thank Christian
Partmann, Eileen Herwig, Volker Springel, Chiaki Kobayashi, Thorsten Naab,
Martyna Chruślińska, Ylva G$\ddot{\rm o}$tberg, Thomas Janka and Selma de Mink
for helpful discussions; Jeremy Blaizot and Leo-Michel Dasnac for developing
RASCAS and Romain Teyssier for developing the code RAMSES and making it
public. This work made extensive use of multiple publicly available software
packages: matplotlib (Hunter, 2007), numpy (Van Der Walt et al., 2011), scipy
(Jones et al., 01) and CoReCon (Garaldi, 2023). The authors would like to
thank the community of developers and those maintaining these packages.

## Data Availability

Data will be provided on reasonable request to the authors.

## References

  * Adams et al. (2023) Adams N. J., et al., 2023, EPOCHS Paper II: The Ultraviolet Luminosity Function from $7.5<z<13.5$ using 110 square arcminutes of deep, blank-field data from the PEARLS Survey and Public Science Programmes (arXiv:2304.13721)
  * Agertz & Kravtsov (2016) Agertz O., Kravtsov A. V., 2016, The Astrophysical Journal, 824, 79
  * Agertz et al. (2011) Agertz O., Teyssier R., Moore B., 2011, Monthly Notices of the Royal Astronomical Society, 410, 1391
  * Agertz et al. (2020) Agertz O., et al., 2020, MNRAS, 491, 1656
  * Algera et al. (2023) Algera H. S. B., et al., 2023, MNRAS, 518, 6142
  * Asada & Ohta (2022) Asada Y., Ohta K., 2022, ApJ, 924, 71
  * Asada et al. (2023a) Asada Y., et al., 2023a, Bursty star formation and galaxy-galaxy interactions in low-mass galaxies 1 Gyr after the Big Bang (arXiv:2310.02314)
  * Asada et al. (2023b) Asada Y., et al., 2023b, Monthly Notices of the Royal Astronomical Society: Letters, 523, L40
  * Atek et al. (2023) Atek H., et al., 2023, JWST UNCOVER: Discovery of $z>9$ Galaxy Candidates Behind the Lensing Cluster Abell 2744 (arXiv:2305.01793)
  * Aubert et al. (2004) Aubert D., Pichon C., Colombi S., 2004, Monthly Notices of the Royal Astronomical Society, 352, 376
  * Barnes (2001) Barnes J. E., 2001, in Hibbard J. E., Rupen M., van Gorkom J. H., eds, Astronomical Society of the Pacific Conference Series Vol. 240, Gas and Galaxy Evolution. p. 135 (arXiv:astro-ph/0010144), doi:10.48550/arXiv.astro-ph/0010144
  * Barnes & Hernquist (1996) Barnes J. E., Hernquist L., 1996, ApJ, 471, 115
  * Barrow et al. (2020) Barrow K. S. S., Robertson B. E., Ellis R. S., Nakajima K., Saxena A., Stark D. P., Tang M., 2020, ApJ, 902, L39
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, Monthly Notices of the Royal Astronomical Society, 488, 3143
  * Bekki (1998) Bekki K., 1998, A&A, 334, 814
  * Belfiore et al. (2022) Belfiore F., et al., 2022, A&A, 659, A26
  * Biernacki & Teyssier (2018) Biernacki P., Teyssier R., 2018, Monthly Notices of the Royal Astronomical Society, 475, 5688
  * Bigiel et al. (2008) Bigiel F., Leroy A., Walter F., Brinks E., De Blok W., Madore B., Thornley M. D., 2008, The Astronomical Journal, 136, 2846
  * Blondin et al. (1998) Blondin J. M., Wright E. B., Borkowski K. J., Reynolds S. P., 1998, The Astrophysical Journal, 500, 342
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Bournaud & Elmegreen (2009) Bournaud F., Elmegreen B. G., 2009, The Astrophysical Journal, 694, L158
  * Bouwens et al. (2015) Bouwens R. J., et al., 2015, The Astrophysical Journal, 803, 34
  * Bouwens et al. (2023a) Bouwens R., Illingworth G., Oesch P., Stefanon M., Naidu R., van Leeuwen I., Magee D., 2023a, Monthly Notices of the Royal Astronomical Society, 523, 1009
  * Bouwens et al. (2023b) Bouwens R. J., et al., 2023b, Monthly Notices of the Royal Astronomical Society, 523, 1036
  * Boylan-Kolchin (2023) Boylan-Kolchin M., 2023, Nature Astronomy, 7, 731
  * Brown & Gnedin (2022) Brown G., Gnedin O. Y., 2022, Monthly Notices of the Royal Astronomical Society, 514, 280
  * Calzetti (2012) Calzetti D., 2012, Star Formation Rate Indicators (arXiv:1208.2997)
  * Carniani et al. (2023) Carniani S., et al., 2023, JADES: The incidence rate and properties of galactic outflows in low-mass galaxies across 3 < z < 9 (arXiv:2306.11801)
  * Chabrier (2003) Chabrier G., 2003, Publications of the Astronomical Society of the Pacific, 115, 763
  * Chevalier & Clegg (1985) Chevalier R. A., Clegg A. W., 1985, Nature, 317, 44
  * Chisholm et al. (2020) Chisholm J., Prochaska J. X., Schaerer D., Gazagnes S., Henry A., 2020, Monthly Notices of the Royal Astronomical Society, 498, 2554
  * Choudhury et al. (2021) Choudhury T. R., Paranjape A., Bosman S. E. I., 2021, MNRAS, 501, 5782
  * Ciesla et al. (2023) Ciesla L., et al., 2023, Identification of a transition from stochastic to secular star formation around $z=9$ with JWST (arXiv:2309.15720)
  * Commerçon et al. (2014) Commerçon B., Debout V., Teyssier R., 2014, A&A, 563, A11
  * Costa et al. (2014) Costa T., Sijacki D., Trenti M., Haehnelt M. G., 2014, MNRAS, 439, 2146
  * Costa et al. (2018a) Costa T., Rosdahl J., Sijacki D., Haehnelt M. G., 2018a, MNRAS, 473, 4197
  * Costa et al. (2018b) Costa T., Rosdahl J., Sijacki D., Haehnelt M. G., 2018b, MNRAS, 479, 2079
  * Costa et al. (2019) Costa T., Rosdahl J., Kimm T., 2019, MNRAS, 489, 5181
  * D’Aloisio et al. (2019) D’Aloisio A., McQuinn M., Maupin O., Davies F. B., Trac H., Fuller S., Upton Sanderbeck P. R., 2019, ApJ, 874, 154
  * Dashyan et al. (2017) Dashyan G., Silk J., Mamon G. A., Dubois Y., Hartwig T., 2017, Monthly Notices of the Royal Astronomical Society, 473, 5698
  * Dekel & Silk (1986) Dekel A., Silk J., 1986, ApJ, 303, 39
  * Dekel et al. (2009) Dekel A., Sari R., Ceverino D., 2009, The Astrophysical Journal, 703, 785
  * Deparis et al. (2019) Deparis N., Aubert D., Ocvirk P., Chardin J., Lewis J., 2019, A&A, 622, A142
  * Díaz-Rodríguez et al. (2018) Díaz-Rodríguez M., Murphy J. W., Rubin D. A., Dolphin A. E., Williams B. F., Dalcanton J. J., 2018, The Astrophysical Journal, 861, 92
  * Díaz-Rodríguez et al. (2021) Díaz-Rodríguez M., Murphy J. W., Williams B. F., Dalcanton J. J., Dolphin A. E., 2021, Monthly Notices of the Royal Astronomical Society, 506, 781
  * Diesing & Caprioli (2018) Diesing R., Caprioli D., 2018, Phys. Rev. Lett., 121, 091101
  * Donnan et al. (2022) Donnan C. T., et al., 2022, Monthly Notices of the Royal Astronomical Society, 518, 6011
  * Dressler et al. (2023a) Dressler A., et al., 2023a, Building the First Galaxies – Chapter 2. Starbursts Dominate The Star Formation Histories of 6 < z <12 Galaxies (arXiv:2306.02469)
  * Dressler et al. (2023b) Dressler A., et al., 2023b, The Astrophysical Journal Letters, 947, L27
  * Dubois & Teyssier (2008) Dubois Y., Teyssier R., 2008, A&A, 477, 79
  * Dubois et al. (2016) Dubois Y., Peirani S., Pichon C., Devriendt J., Gavazzi R., Welker C., Volonteri M., 2016, Monthly Notices of the Royal Astronomical Society, 463, 3948
  * Dubois et al. (2021) Dubois Y., et al., 2021, A&A, 651, A109
  * Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
  * Eide et al. (2020) Eide M. B., Ciardi B., Graziani L., Busch P., Feng Y., Di Matteo T., 2020, MNRAS, 498, 6083
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Endsley et al. (2023a) Endsley R., et al., 2023a, The Star-forming and Ionizing Properties of Dwarf z 6-9 Galaxies in JADES: Insights on Bursty Star Formation and Ionized Bubble Growth (arXiv:2306.05295)
  * Endsley et al. (2023b) Endsley R., Stark D. P., Whitler L., Topping M. W., Chen Z., Plat A., Chisholm J., Charlot S., 2023b, Monthly Notices of the Royal Astronomical Society, 524, 2312
  * Fan et al. (2006) Fan X., et al., 2006, AJ, 132, 117
  * Farcy et al. (2022) Farcy M., Rosdahl J., Dubois Y., Blaizot J., Martin-Alvarez S., 2022, MNRAS, 513, 5000
  * Faucher-Giguère et al. (2013) Faucher-Giguère C.-A., Quataert E., Hopkins P. F., 2013, Monthly Notices of the Royal Astronomical Society, 433, 1970
  * Federrath (2016) Federrath C., 2016, Journal of Physics: Conference Series, 719, 012002
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, The Astrophysical Journal, 761, 156
  * Ferland et al. (1998) Ferland G., Korista K., Verner D., Ferguson J., Kingdon J., Verner E., 1998, Publications of the Astronomical Society of the Pacific, 110, 761
  * Ferrara et al. (2023) Ferrara A., Pallottini A., Dayal P., 2023, Monthly Notices of the Royal Astronomical Society, 522, 3986
  * Fielding et al. (2018) Fielding D., Quataert E., Martizzi D., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3325
  * Finkelman et al. (2012) Finkelman I., Brosch N., Funes J. G., Barway S., Kniazev A., Väisänen P., 2012, Monthly Notices of the Royal Astronomical Society, 422, 1384
  * Finkelstein et al. (2015) Finkelstein S. L., et al., 2015, The Astrophysical Journal, 810, 71
  * Finlator et al. (2018) Finlator K., Keating L., Oppenheimer B. D., Davé R., Zackrisson E., 2018, MNRAS, 480, 2628
  * Fiore et al. (2023) Fiore F., Ferrara A., Bischetti M., Feruglio C., Travascio A., 2023, ApJ, 943, L27
  * Fujimoto et al. (2022) Fujimoto S., et al., 2022, JWST and ALMA Multiple-Line Study in and around a Galaxy at $z=8.496$: Optical to FIR Line Ratios and the Onset of an Outflow Promoting Ionizing Photon Escape (arXiv:2212.06863)
  * Fujimoto et al. (2023) Fujimoto S., et al., 2023, The Astrophysical Journal Letters, 949, L25
  * Furlanetto & Mirocha (2022) Furlanetto S. R., Mirocha J., 2022, Monthly Notices of the Royal Astronomical Society, 511, 3895
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, arXiv e-prints, p. arXiv:2304.02038
  * Garaldi (2023) Garaldi E., 2023, The Journal of Open Source Software, 8, 5407
  * Garaldi et al. (2021) Garaldi E., Pakmor R., Springel V., 2021, MNRAS, 502, 5726
  * Garaldi et al. (2022) Garaldi E., Kannan R., Smith A., Springel V., Pakmor R., Vogelsberger M., Hernquist L., 2022, MNRAS, 512, 4909
  * Gatto et al. (2016) Gatto A., et al., 2016, Monthly Notices of the Royal Astronomical Society, 466, 1903
  * Geen et al. (2015) Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3248
  * Geen et al. (2020) Geen S., Bieri R., Rosdahl J., de Koter A., 2020, Monthly Notices of the Royal Astronomical Society, 501, 1352
  * Gnedin (2000) Gnedin N. Y., 2000, The Astrophysical Journal, 542, 535
  * Gnedin (2014) Gnedin N. Y., 2014, The Astrophysical Journal, 793, 29
  * Gnedin (2016) Gnedin N. Y., 2016, The Astrophysical Journal, 833, 66
  * Gnedin & Abel (2001) Gnedin N. Y., Abel T., 2001, New Astronomy, 6, 437
  * Gnedin & Madau (2022) Gnedin N. Y., Madau P., 2022, Modeling Cosmic Reionization (arXiv:2208.02260)
  * Gnedin et al. (2000) Gnedin N. Y., Ferrara A., Zweibel E. G., 2000, The Astrophysical Journal, 539, 505
  * Götberg et al. (2019) Götberg Y., de Mink S. E., Groh J. H., Leitherer C., Norman C., 2019, A&A, 629, A134
  * Götberg et al. (2020) Götberg Y., de Mink S. E., McQuinn M., Zapartas E., Groh J. H., Norman C., 2020, A&A, 634, A134
  * Götberg et al. (2023) Götberg Y., et al., 2023, arXiv e-prints, p. arXiv:2307.00074
  * Greig et al. (2017) Greig B., Mesinger A., Haiman Z., Simcoe R. A., 2017, MNRAS, 466, 4239
  * Greig et al. (2019) Greig B., Mesinger A., Bañ ados E., 2019, Monthly Notices of the Royal Astronomical Society, 484, 5094
  * Grimmett et al. (2020) Grimmett J., Karakas A. I., Heger A., Müller B., West C., 2020, Monthly Notices of the Royal Astronomical Society, 496, 4987
  * Grisdale et al. (2019) Grisdale K., Agertz O., Renaud F., Romeo A. B., Devriendt J., Slyz A., 2019, Monthly Notices of the Royal Astronomical Society, 486, 5482
  * Guillet & Teyssier (2011) Guillet T., Teyssier R., 2011, Journal of Computational Physics, 230, 4756
  * Guszejnov et al. (2022) Guszejnov D., Grudić M. Y., Offner S. S. R., Faucher-Giguère C.-A., Hopkins P. F., Rosen A. L., 2022, Monthly Notices of the Royal Astronomical Society, 515, 4929
  * Haardt & Madau (1995) Haardt F., Madau P., 1995, arXiv preprint astro-ph/9509093
  * Hahn et al. (2020) Hahn O., Rampf C., Uhlemann C., 2020, Monthly Notices of the Royal Astronomical Society, 503, 426
  * Haid et al. (2018) Haid S., Walch S., Seifried D., Wünsch R., Dinnbier F., Naab T., 2018, Monthly Notices of the Royal Astronomical Society, 482, 4062
  * Harikane et al. (2022) Harikane Y., et al., 2022, The Astrophysical Journal Supplement Series, 259, 20
  * Harikane et al. (2023) Harikane Y., et al., 2023, The Astrophysical Journal Supplement Series, 265, 5
  * Haro et al. (2023a) Haro P. A., et al., 2023a, Spectroscopic confirmation of CEERS NIRCam-selected galaxies at $\boldsymbol{z\simeq 8-10}$ (arXiv:2304.05378)
  * Haro et al. (2023b) Haro P. A., et al., 2023b, Spectroscopic verification of very luminous galaxy candidates in the early universe (arXiv:2303.15431)
  * Hartley & Ricotti (2016) Hartley B., Ricotti M., 2016, Monthly Notices of the Royal Astronomical Society, 462, 1164
  * Hasegawa & Semelin (2012) Hasegawa K., Semelin B., 2012, Monthly Notices of the Royal Astronomical Society, 428, 154
  * Hayward & Hopkins (2016) Hayward C. C., Hopkins P. F., 2016, Monthly Notices of the Royal Astronomical Society, 465, 1682
  * Heintz et al. (2023a) Heintz K. E., et al., 2023a, Dilution of chemical enrichment in galaxies 600 Myr after the Big Bang (arXiv:2212.02890)
  * Heintz et al. (2023b) Heintz K. E., et al., 2023b, Extreme damped Lyman-$\alpha$ absorption in young star-forming galaxies at $z=9-11$ (arXiv:2306.00647)
  * Helton et al. (2023) Helton J. M., et al., 2023, The JWST Advanced Deep Extragalactic Survey: Discovery of an Extreme Galaxy Overdensity at $z=5.4$ with JWST/NIRCam in GOODS-S (arXiv:2302.10217)
  * Hennebelle & Chabrier (2011) Hennebelle P., Chabrier G., 2011, The Astrophysical Journal Letters, 743, L29
  * Heyer & Brunt (2004) Heyer M. H., Brunt C. M., 2004, ApJ, 615, L45
  * Hoag et al. (2019) Hoag A., et al., 2019, The Astrophysical Journal, 878, 12
  * Hopkins et al. (2009) Hopkins P. F., et al., 2009, Monthly Notices of the Royal Astronomical Society, 397, 802
  * Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, Monthly Notices of the Royal Astronomical Society, 445, 581
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1578
  * Hu (2019) Hu C.-Y., 2019, MNRAS, 483, 3363
  * Huertas-Company et al. (2023) Huertas-Company M., et al., 2023, Galaxy Morphology from $z\sim 6$ through the eyes of JWST (arXiv:2305.02478)
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27
  * Hummer & Storey (1987) Hummer D. G., Storey P. J., 1987, MNRAS, 224, 801
  * Hunter (2007) Hunter J. D., 2007, Computing in science & engineering, 9, 90
  * Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, Monthly Notices of the Royal Astronomical Society, 439, 725
  * Iyer et al. (2018) Iyer K., et al., 2018, The Astrophysical Journal, 866, 120
  * Jacobs et al. (2023) Jacobs C., et al., 2023, ApJ, 948, L13
  * Jin et al. (2023) Jin S., et al., 2023, A&A, 670, L11
  * Jones et al. (01 ) Jones E., Oliphant T., Peterson P., et al., 2001–, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Jung et al. (2020) Jung I., et al., 2020, ApJ, 904, 144
  * Jung et al. (2023a) Jung I., et al., 2023a, CEERS: Diversity of Lyman-Alpha Emitters during the Epoch of Reionization (arXiv:2304.05385)
  * Jung et al. (2023b) Jung I., et al., 2023b, arXiv e-prints, p. arXiv:2304.05385
  * Kannan et al. (2015) Kannan R., Macciò A. V., Fontanot F., Moster B. P., Karman W., Somerville R. S., 2015, Monthly Notices of the Royal Astronomical Society, 452, 4347
  * Kannan et al. (2020) Kannan R., Marinacci F., Vogelsberger M., Sales L. V., Torrey P., Springel V., Hernquist L., 2020, Monthly Notices of the Royal Astronomical Society, 499, 5732
  * Kannan et al. (2021) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2021, Monthly Notices of the Royal Astronomical Society, 511, 4005
  * Kannan et al. (2022) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2022, MNRAS, 511, 4005
  * Kartaltepe et al. (2023) Kartaltepe J. S., et al., 2023, The Astrophysical Journal Letters, 946, L15
  * Katz (2022) Katz H., 2022, Monthly Notices of the Royal Astronomical Society, 512, 348
  * Katz et al. (2018) Katz H., Kimm T., Haehnelt M. G., Sijacki D., Rosdahl J., Blaizot J., 2018, Monthly Notices of the Royal Astronomical Society, 483, 1029
  * Katz et al. (2021) Katz H., et al., 2021, Monthly Notices of the Royal Astronomical Society, 507, 1254
  * Katz et al. (2022) Katz H., et al., 2022, Monthly Notices of the Royal Astronomical Society, 518, 270
  * Keating et al. (2020) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020, MNRAS, 491, 1736
  * Keller & Kruijssen (2022) Keller B. W., Kruijssen J. D., 2022, Monthly Notices of the Royal Astronomical Society, 512, 199
  * Keller et al. (2022) Keller B. W., Kruijssen J. D., Chevance M., 2022, Monthly Notices of the Royal Astronomical Society, 514, 5355