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

We present SPICE, a new suite of RHD 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 scenarios such as hypernova explosions and radiation pressure. Subtle differences in the behaviour of supernova feedback drive profound differences in reionisation histories, with burstier forms of feedback causing earlier reionisation. We also find that some global galaxy properties, such as the dust-attenuated luminosity functions and star formation main sequence, remain degenerate between models. 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 ISM/CGM in different ways, and in turn strongly affecting the escape of 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. 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. Statistical samples of post-reionisation galaxy morphologies probed with JWST, ALMA and MUSE can constrain stellar feedback and models of cosmic reionisation.


INTRODUCTION
About 300, 000 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  = 5 (e.g.Gaikwad et al. 2023), however, the vast majority of hydrogen has transitioned to an ionised state.This ∼ 1 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;Robertson et al. 2015;Eide et al. 2020), with quasars now thought to provide only a minor contribution (Kulkarni et al. 2019b;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 ΛCDM paradigm of galaxy formation, such processes must be invoked in order to ★ E-mail: abhagwat@mpa-garching.mpg.deprevent the overproduction of stars even at  > 6 (Hopkins et al. 2014;Costa et al. 2014;Somerville & Davé 2015).
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.6eV) escape (Wise & Cen 2009;Kimm & Cen 2014;Trebitsch et al. 2017Trebitsch et al. , 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  > 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. ≳ 20 pc1 .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 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 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  * ≤ 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 (Shen et al. 2023;Sun et al. 2023a,b;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,b;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), Renaissance (O'Shea et al. 2015), Aurora (Pawlik et al. 2017), Technicolor Dawn (Finlator et al. 2018), SPHINX (Rosdahl et al. 2018;Katz et al. 2018;Rosdahl et al. 2022a), CoDa (Ocvirk et al. 2016(Ocvirk et al. , 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 ≈ 28pc at  = 5) in a cosmological volume ( box ≈ 14.8cMpc).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.

SIMULATIONS
In this section, we describe the simulations performed and postprocessed 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.11describing the novel combination of galaxy formation and feedback models we include in our simulations.

Initial conditions
The initial conditions (ICs) for the simulations are generated at  = 30 using MUSIC2-monofonIC2 (Hahn et al. 2020;Michaux et al. 2020) with the 2PPT/2LPT approach.A ΛCDM cosmological model is adopted, with parameters Ω Λ = 0.6901, Ω m = 0.3099, Ω b = 0.0489,  0 = 67.74km/s/Mpc,  s = 0.9682,  8 = 0.8159 (Planck Collaboration et al. 2016).We assume primordial gas mass fraction  = 0.7545 for Hydrogen (H) and  = 0.2455 for Helium (He), along with a metallicity floor (see 2.3).All simulation boxes have a length  box = 10 cMpc/ℎ, with 512 3 dark matter particles, i.e. a mean mass  dm = 6.38 × 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  (Porqueres et al. 2020), therefore providing accurate ICs for codes like RAMSES-RT.We also use the "fixing" technique of Angulo & Pontzen (2016), which fixes the amplitude of each density mode to the expectation of the input power spectrum, therefore ensuring that the ICs generated are representative of the average Universe.We have produced halo mass functions for our ICs and verified that they closely agree with expectations for an average region.

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  = 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 (Δ ≈ 100pc), while at more refined levels we adopt a conjugate gradient solver to improve the computational speed.

Cooling and gas thermochemistry
Gas cooling and heating is calculated as described in Rosdahl et al. (2013) 3 .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  > 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  ≤ 10 4 K we adopt the fine structure cooling rates from Rosen & Bregman (1995), allowing the gas to cool radiatively down to ≈ 300K.We also assume a homogeneous initial gas metallicity floor of  = 3.2 × 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  = 10 4 K.We also adopt the onthe-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.

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 ℓ determines the width of each gas cell as Δ ℓ = 0.5 ℓ  box .We allow refinement levels in the range ℓ = 9 − 16.At the coarsest level (ℓ min = 9) we have a minimum physical resolution of 4.8 kpc at  = 5, while at the finest level (ℓ max = 16) the maximum physical resolution is ≈ 28 pc at  = 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: (i) if  dm + Ω  Ω  b > 8 ×  dm (quasi-Lagrangian refinement criterion), where  dm and  b are the total dark matter and baryonic (i.e.gas + stars) mass within the cell, and  dm is the mean dark matter particle mass.
(ii) The cell width is larger than 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 ℓ max at all redshifts.This allows for the maximum physical resolution to be the highest at very high redshifts (≈10.72 pc at  = 20).

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  vir and mass  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:  SPH = 32,  HOP = 16,  th = 80, and  poisson = 4.We also identify galaxies with AdaptaHOP.Here we require at least 100 stellar particles per galaxy.The galaxy-finder parameters are:  SPH = 10,  HOP = 10,  th = 10 3 , and  poisson = 4.

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  H ≥ 10cm −3 , we adopt a standard Schmidt law, for which the star formation density is given as, where  is the gas density,  ff = √︁ 3/(32 ) is the freefall time, and  ff is the star formation efficiency per freefall time.Typically, in galaxy formation simulations the value of  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)  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)  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 where  = ln(/),   is the variance of ,  = −1/2 2  is the mean logarithmic density, and  is the mean density of the cell.
The local star formation efficiency  ff , as computed in Hennebelle & Chabrier (2011);Federrath & Klessen (2012), is where we assume that gas with density larger than a critical value,  crit , is converted into stars.The parameter  = 1/0.49takes into account the uncertainty in free-fall timescales for gas with different densities and  = 0.5 accounts for the fact that not all gas above  crit is converted into stars (both quantities are derived in Federrath & Klessen 2012).The critical density for star formation ( 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 ( vir < 1) and M < 1.The critical density accounting for both supersonic and subsonic regimes is Here, is the local virial parameter, which is an indicator of local stability.Δ is the cell size,  is the gravitational constant,  is the 1D turbulent velocity (see 2.7) and M = /  is the local Mach number, where   is the speed of sound.
This model therefore allows for two star formation channels.In the first channel, for which  vir < 1 (independent of local M), the entire computational cell is gravitationally unstable.In the second channel, collapse owing to large density fluctuations caused by supersonic turbulence, i.e.M ≫ 1 (see Figure 1 in Kretschmer & Teyssier 2020), becomes possible even if  vir > 1.For reference, typical conditions for star-forming regions within the Milky way are  vir ≃ 3 − 10 with M ≃ 10−20 (Spilker et al. 2022), and correspond to  ff ≃ 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 ⊙ (≈ 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).

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 andSemenov et al. 2018).This is given by where the turbulent kinetic energy  T = 3/2 2 T , is related to the turbulent pressure by  T = 2/3 T , with  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 (  and   , respectively), which can be written as and where    is the mean flow viscous stress tensor (this term is evaluated as in Schmidt & Federrath 2011), ṽ is the mean flow variable (defined as ṽ = /) and the sum is calculated over the nearest neighbours of the computational cell in question.This model has two important parameters, the turbulent viscosity   and the dissipation time-scale  diss , which are related to the cells' size by 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. 2017Trebitsch et al. , 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  = √︁ 2 T /.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).

Supernova feedback
We adopt the mechanical feedback scheme introduced by Kimm & Cen (2014) and implemented as in Kimm et al. (2015) (see Hopkins et al. 2014 for a similar setup).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  SN which depends on whether the adiabatic (Sedov-Taylor) phase is resolved.The model introduces a parameter , which gives the ratio between the swept-up mass  swept and the ejected mass  ej for each neighbouring cell (number of neighbouring cells  nbor = 48) as This is compared to a threshold value of where  51 is the explosion energy of an individual SN in units of 10 51 erg,  H is the hydrogen number density in cm −3 and and Here  host is gas density of the SN host cell, and  sn determines what fraction of the gas mass ( ej +  host Δ 3 ) is re-distributed to the host cell of a SN (see Fig. 15 of Kimm & Cen 2014).The parameter  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  >  tr the adiabatic phase is not resolved and the momentum during the snowplow phase ( SN,snow ) is injected to the neighbouring cells.Otherwise the momentum during the adiabatic phase ( SN,ad ) is injected, i.e.
where  e ( ) = 1 − −1 3(  tr −1) smoothly connects the two regimes.The input momentum during the snowplow phase is  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.
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  S is calculated and compared with the cell width Δ.If the Strömgen radius is not well resolved (Δ >>  s ), we adopt a final input momentum of where  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  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 −1 ⊙ , which results in 15 SN per star particle for our minimum stellar mass of 970 M ⊙ .

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.

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

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 & Kruĳssen 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 × 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(Rosdahl et al. , 2022b)).

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. 2006Kobayashi 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 ∼ 1% in the local universe (Podsiadlowski et al. 2004).Simulations of starclusters 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 ∼ 1 Gyr, hence making them relevant to study during reionisation.
In hyper-sn, we assign each stellar particle a metallicity dependant HN fraction  HN given as where  * is the metallicity of the stellar particle. HN represents the fraction of all SN events a stellar particle undergoes that are classified as HN.We assume an initial HN fraction of  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 × 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×10 51 ergs, respectively (Sukhbold et al. 2016;Díaz-Rodríguez et al. 2018, 2021).

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).
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 c = 0.1.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 c, in the range (0.01 − 0.2), while finding convergence with values as low as c = 0.1 (Wu et al. 2019).Therefore, we use results from these studies and adopt c = 0.1.

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 ergs −1 M −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 luminosityweighted 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  d ≡ (/ ⊙ )    H , where  is the local metallicity,  d = 1 −  HII is the neutral fraction of gas that holds dust, and  H is the hydrogen number density.This ensures that the dust density depends not only on enrichment but also on the local thermal state.While dust is generally not expected to survive in photo-ionised gas (Finkelman et al. 2012;Kannan et al. 2020), dust-to-ionised-gas ratios as high as 0.01 have been observed (Harper & Low 1971;Contini & Contini 2003;Laursen et al. 2009) in ionised media.The  d parameter in the previous equation allows for non-zero (albeit small) dust effects even in ionised media.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).

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

Overview
Primordial gas cools to form the first stellar particle at  = 22.4.Feedback from massive stars via photo-ionisation begins immediately, and supernova feedback kicks in as early as 3 Myr (10 Myr) after a stellar particle is born in smooth-sn and hyper-sn (bursty-sn).Galaxies continue to grow and the first  * > 10 7 M ⊙ dwarf galaxies form as early as  = 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.
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 B shows the infrared radial flux emerging from galaxies, which includes intrinsic emission from stars, dust emission and multi-scattered IR radiation.We see a network of infrared-bright galaxies with a well established background by  = 5.
In panel C we zoom further into smaller scales, showing a single galaxy.We plot the radial flux of hydrogen-ionising radiation (ℎ > 13.6 eV) emerging from a  * = 2 × 10 10 M ⊙ 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 4 ×  vir (≈ 120kpc) centred on a ≈ 2.7 × 10 11 M ⊙ halo at  = 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  = 5.The projections The right-hand panels of Figure 2 show metallicity projections for the various feedback models at  = 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 ( > 1 ⊙ ) at the very centers of halos.hyper-sn produces a distribution similar to smooth-sn, where the enriched regions are compact and with higher maximum metallicities.(the gas-phase metallicity of the most massive halo at ∼ 0.9 ⊙ , ∼ 0.49 ⊙ for bursty-sn,smooth-sn and hyper-sn, respectively).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).

Cosmic star formation history
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  = 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  = 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 multiwavelength 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  ≳ 10, bursty-sn at  ≲ 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  = 10, across the mass range, bursty-sn shows the largest scatter (≈ 0.8 dex) followed by hyper-sn (≈ 0.5 dex) and finally smooth-sn (≈ 0.3 dex).By  = 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  * < 10 7 M ⊙ .While smooth-sn and hyper-sn show minimal evolution in scatter between  = 7 and  = 5, the scatter for bursty-sn drops to ≈ 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 ±0.25 around the redshifts shown.At  = 10, the vast majority of observations lie in a mass range not probed well by SPICE, however, smooth-sn shows the best agreement with observations.At  = 7 and  = 5, the median relations from all models is in excellent agreement with data for galaxies with  * > 10 7 M ⊙ .
Figure 5 shows the ratio between stellar mass and halo mass, which we view as a proxy for global star formation efficiency (SFE),  * /( halo   ), as a function of halo mass  halo , where 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  ≈ 9, increases, and by  = 7 becomes comparable to the one of the bursty-sn model.The HN rate in haloes with  halo > 10 10 M ⊙ reaches 1% (by  ≈ 9), therefore, the energetic component of feedback becomes negligible 4 .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 (≈ 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 =7 observations from (Stefanon et al. 2021).We note that the trend reversal at low halo masses is due to these halos undergoing their first starburst events (average stellar age ∼ 5 − 15 Myr), which temporarily boosts their SFEs.
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  = 5.Finally, the hyper-sn model emphasizes the importance of striking a balance between bursty and smooth feedback mechanisms.

Reionisation histories
In Figure 6 we show thin-slices (depth of 100 ckpc) of ionised fractions ( HII ) in the three feedback models representing the progress of 4 The range  = 9 − 10 also marks the onset of galaxies settling into disks.
We plan to quantify this morphological evolution in a follow up work.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  = 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  = 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  = 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 8 × 8 grid and mark the location of the most luminous galaxies within each grid sub-division.The galaxies marked within each ionised bubble account for > 95% of the escaping ionising photon budget, effectively representing the galaxies that are responsible for the formation and growth of such bubble.Each galaxy is coded by its morphology (stars and spirals represent dispersion-and rotation-dominated galaxies, respectively; see section 3.5 for details), LyC escape fraction (different colors) and escaping ionising luminosity (marker sizes).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 reionises5 the Universe by  = 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, x HI , dropping to ≈ 0.9 already at  ≈ 10.From  ≈ 7, 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. 2019a;Keating Figure Thin-slice projections showing HII fraction in the bursty-sn, smooth-sn and hyper-sn models (from top to bottom) at  = 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, with the sizes of markers representing escaping ionising photon luminosity.Already at  = 10 ionised regions show differences in sizes, with bursty-sn and hyper-sn producing the largest and smallest, respectively.By  = 5, reionisation is complete in bursty-sn, whereas significant neutral patches persist in the other two models.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  = 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 ≳ 100 cMpc/ℎ 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  = 5 for the smooth-sn, bursty-sn and hyper-sn model are ≈ 5 × 10 10 M ⊙ , ≈ 2 × 10 11 M ⊙ and ≈ 1 × 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 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  < 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.

Luminosity functions at 1500 Å
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,a;Leung et al. 2023).
From the top panel of Figure 8 we see that at  = 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  = 7. Similarly to  = 10, smooth-sn shows a larger number (by 0.3-0.4dex) of galaxies at all magnitudes compared to the other models.
By  = 5, smooth-sn produces extremely UV-bright ( 1500 ∼ −23) objects and an abundance of galaxies at all magnitudes.Also hyper-sn results in a similar abundance of objects at  1500 ∼ −23, while showing less dimmer galaxies (−20 <  1500 < −12).Meanwhile, bursty-sn lies in between the other two models.
At  = 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  1500 < −17.Dust attenuation becomes more at  = 7 for smooth-sn and hyper-sn, when an effect is visible at the bright end of the LF.Indeed, it becomes relevant at  1500 < −16 for smooth-sn, and at  1500 < −18 for hyper-sn, while in bursty-sn dust attenuation is minor even at the brightest end of the LF.At  = 5, dust attenuation for smooth-sn is significant (> 1 dex) at  1500 < −18, while in hyper-sn it is present only at  1500 < −20.Finally, for bursty-sn it is minimal (≈ 0.2dex) for  1500 < −17.At  = 10, while bursty-sn and smooth-sn agree very well with data, hyper-sn shows a deficit of bright galaxies.However, both at  = 7 and  = 5, dust-attenuated LFs from all three models are in excellent agreement with observations for  1500 < −16.
Overall, in this section we show that despite the intrinsic LFs being extremely different, the dust-attenuated ones are very similar below  1500 < −16, suggesting that UVLFs cannot be directly used to probe stellar feedback at high redshifts.

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 ≈ 50 M ⊙ yr −1 .Within the same haloes, bursty-sn generates systems with 10 times lower star formation rates and (0.8 − 1) 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 ≈ 50 M ⊙ 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 ≈ 107 M ⊙ (see Figure 10).
While in Figure 9 we show 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 / (Dubois et al. 2016;Pillepich et al. 2019) as a proxy for morphology, where  is the rotational velocity and  is the 3D velocity dispersion.This ratio classifies the degree to which a galaxy is supported by dispersion or rotation.We calculate / for both stars and gas for each galaxy with  * > 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 galaxy6 .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  vir adopting an H 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 volume emissivity for each gas cell as (Hummer & Storey 1987): where  B = 2.753 × 10 −14  1.5 /(1.0 + (/2.74) 0.407 ) 2.242 cm 3 s −1 with  = 315614/ is the case-B recombination coefficient (Hui & Gnedin 1997),  is the gas temperature,  e and  p are the number densities of electrons and protons, and  B () is the conversion probability per recombination event, which is ≈ 0.451 for  = 10 4 K.
To calculate gas kinematics we select gas cells with  H > 10 −6 erg s −1 cm −3 , noting that this choice does not affect our qualitative results 7 .Following the process used in various theoretical and observational studies (Pillepich et al. 2019;Rizzo et al. 2020), / is calculated using  at the peak of the rotation curve and  as the mean velocity dispersion.We classify galaxies with / ≥ 1 (for both stars and gas) as rotationally-supported, otherwise as dispersionsupported.
In Figure 10 we show the median / 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 ≈ 18% at  = 7 to ≈ 6%  = 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 ≈ 31% at  = 10 to ≈ 42% at  = 5), smooth-sn has a similar trend but the increase is more drastic (from ≈ 27% to ≈ 53%), while hyper-sn shows a mild increase in this fraction between  = 10 and 7, followed by a drop in the range  = 7 − 5.
At  = 10, both stellar and gas components are indistinguishable between models.Similarly at  = 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  = 5 for both stellar and gas components.smooth-sn shows a population Figure 10.Median / values for stars (top row) and gas (bottom row) as a function of stellar mass for galaxies in the various feedback models at  = 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  = 5, with smooth-sn, bursty-sn and hyper-sn producing mainly galaxies which are rotationally-supported, dispersion-supported, and mixed, respectively.
of highly rotationally-supported galaxies (/ > 2) and very large scatter at  * ≳ 10 8 M ⊙ , while bursty-sn produces a dispersionsupported galaxy population with a small scatter around the median.
The hyper-sn model exhibits a transition at  * ≈ 10 8 M ⊙ : 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.We note that the gaseous component shows a larger degree of rotational support as compared to the stellar component (i.e. star / star <  gas / gas ), a trend particularly prominent in the smooth-sn and hyper-sn models.The number of galaxies classified as rotationally-vs dispersionsupported is also widely different for the three models, especially at  = 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 / as an indicator at  > 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.

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 (  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  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  esc to each galaxy, and use  gas / gas (as calculated in section 3.5) to sort galaxies into rotationally-and dispersionsupported categories.In Figure 11 we show the luminosity-weighted mean 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  esc of rotation-and dispersiondominated systems, with the latter exhibiting ≈ 10 − 50× 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  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.The relative difference between the  esc for rotationvs 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 com- pared to the other two models, where rotation-dominated systems are preferentially produced, especially at the highest masses.In Figure 6 we show the connection between galaxy morphologies, escape fraction and their resulting ionised bubbles.At  = 10 smooth-sn and hyper-sn exhibit a mix of star symbols, i.e dispersion-dominated systems, and spiral symbols, i.e rotation-dominated systems, while the latter dominate by  = 5.Differently, bursty-sn is dominated by star symbols (i.e.dispersion-dominated systems) at all redshifts.The morphological trends are also reflected in the  esc , where smooth-sn and hyper-sn largely show small escape fractions (pink/orange symbols), whereas in bursty-sn larger values are more prominent (yellow/orange symbols).We note that particularly at  = 7, the largest ionised bubble contains a network of yellow star symbols (i.e.strong LyC leakers which are dispersion-dominated) as in bursty-sn, whereas bubbles in which there are more pink spiral symbols (i.e.rotation-dominated systems with low escape fractions) are smaller in size.By  = 5, the strong LyC leakers (star symbols) in bursty-sn are able to completely reionise their surroundings, whereas this is not the case for the rotation-dominated galaxies in smooth-sn and hyper-sn (which tend to be massive, see Figure 10) We next quantify differences in reionisation histories by looking at the net intrinsic and escaping emissivities from all galaxies, defined as  esc  ion , where  ion is the production rate of photons with   > 13.6eV calculated by integrating the stellar SED, and  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 ≈ 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  * > 10 7 M ⊙ to both intrinsic and escaping emissivities.Differently, the intrinsic emissivity in smooth-sn and hyper-sn is largely dominated by galaxies with  * ⪆ 10 9 M ⊙ (≈ 50% and ≈ 65%, respectively).While the escaping emissivity shows a distribution similar to the intrinsic one (i.e.dominated by galaxies with  * ⪆ 10 9 M ⊙ ), hyper-sn exhibits an even contribution from galaxies with  * > 10 7 M ⊙ .Therefore, in bursty-sn and hyper-sn, galaxies with  * > 10 7 M ⊙ contribute evenly to reionisation, whereas in smooth-sn ≈ 46% of the escaping emissivity is produced by  * > 10 9 M ⊙ galaxies.
As measuring  ion at high redshifts is challenging, other quantities are used to evaluate the ionising photon production, such as the ionising photon production efficiency  ion,0 =  ion / 1500 , where  1500 is the intrinsic luminosity of the galaxy at 1500 Å.  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  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 ≈ (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,b;Saxena et al. 2023a) find median values of log( ion,0 ) ≈ 24.8 − 25.7 for galaxies in the range  ∼ 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  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 rotationor 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  esc (≈ 4× higher than smooth-sn) and reionises the simulation volume by  = 5.1.

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.

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,  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  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 (≈ 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  > 5, and estimate a disk (pure disk + disks with bulges) fraction in the range ∼ 30 − 40% at  ∼ 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 equivalent width and H emission line luminosities (Endsley et al. 2023a,b;Simmonds et al. 2023a) to characterise modes of star formation in JWST galaxies, finding that a bursty SFH is prevalent in galaxies with masses  * ∼ 10 7−10 M ⊙ in the range  = 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, supporting 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  = 5), while galaxies with a non-bursty or smooth SFH are mainly rotationsupported 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 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 /) 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  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  esc , and to constrain stellar feedback.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 and [O III] in low mass galaxies ( * < 10 10 M ⊙ ) in the range  ∼ 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  * ≈ 10 8 M ⊙ galaxy at  ∼ 8.5 is likely associated to a starburst which consequently promotes strong ionising photon escape.It is important to note, however, that in obscured systems (such as SMG/ULIRGs) ionising photon escape into the IGM can be strongly suppressed due to optical depth effects despite incidence of outflows (Cen 2020).Moreover, nonbursty (but high) star-formation (as seen in the smooth-sn model) fails to drive strong outflows, which allows the gas to settle within the disk boosting the gas fractions of galaxies and consequently suppressing  esc (see Yoo et al. 2020 for similar findings).'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.

Numerical uncertainties and physical limitations
The range of physics included in SPICE, along with a physical resolution of ∼ 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 ∼ (28 − 84)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 in-gredient in the feedback modelling (e.g.runaway stars Andersson et al. 2020;Steinwandel et al. 2023), 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 Kimm & Cen 2014).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(Götberg et al. , 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.1.We note that at the tail end of reionisation, the I-fronts can travel as fast as 0.1 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 ≳ 100cMpc/ℎ are required (Iliev et al. 2014;Gnedin & Madau 2022).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  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. 2021Koudmani et al. , 2022;;Sharma et al. 2022).At the same time, observational evidence is mounting supporting the presence of a population of ∼ (10 5 − 10 8 ) M ⊙ 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.

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 ≈ 28 pc at  = 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  = 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  > 9, akin to bursty-sn.As the hypernova rate declines (due to metal enrichment), below  ∼ 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  = 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  = 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  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 ( 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  = 10, by  = 7 all models exhibit noticeable levels of attenuation, which becomes even more pronounced at  = 5.All three models show an excellent agreement with observations, especially at  = 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 / (section 3.5), where we classify galaxies with / ≥ 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)  * ≈ 10 8 M ⊙ typically show strong dispersion (rotational) support.
• The differences in morphology translate into variations in Lyman continuum (LyC) escape fractions (  esc ; section 3.6 and Figure 11).We find that  esc of dispersion-dominated systems is enhanced by a factor of ≈ 10 − 50 in comparison to that of rotationdominated 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  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  esc .The connection between morphologies and  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 escaping emissivity (  esc  ion , dashed lines Figure 12) in the three models is dominated by galaxies in different mass ranges.In smooth-sn the largest contribution (≈ 46%) comes from  * > 10 9 M ⊙ galaxies, whereas bursty-sn and hyper-sn show roughly uniform contribution from galaxies in the range  * ≈ 10 7−10 M ⊙ .Upon comparing the escaping and intrinsic emissivities, we find a global (i.e. from the entire galaxy populations at all redshifts) escape fraction of ≈ 6.8, 1.8 and 1.5% for bursty-sn, smooth-sn and hyper-sn, respectively.
• The ionising photon production efficiency ( ion,0 =  ion / 1500 ) evaluated over the full duration of the simulation shows comparable median trends (see Figure 13) for the three models, with differences of ≈ 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  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  esc to connect feedback to reionisation using kinematics derived from stars and H bright gas.However, imprints of feedback will manifest in a wide range of observables.The multiphase ISM model, along with the additional radiation transport in IR and optical, allows us to investigate a variety of emission lines (such as H, [C II], [O III] and Ly).Indeed, all these probes offer a unique opportunity to make predictions for and comparisons to observations across a wide range of wavelengths, from state-ofthe-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.

Figure 1 .
Figure 1.The SPICE simulations and their scope and dynamic range.In panel A we plot the volume-weighted gas density  = 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 Δ ≈ 28pc.In panel D we show the surface brightness of the same galaxy in [CII].In panel E, we zoom on a region of size 4 ×  vir (≈ 120kpc) centred on a ≈ 2.7 × 10 11 M ⊙ 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.

Figure 2 .Figure 3 .
Figure 2. Temperature (left column) and gas-phase metallicity (right column) of the various simulations at  = 5.The projections shown cover a region spanning  box ×  box /3 ×  box .Consequences of varying the stellar feedback are imprinted in the thermal state and enrichment of the large-scale structure.

Figure 4 .
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  = 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  = 6 (Iyer et al. 2018).

Figure 5 .
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  = 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 fromTacchella et al. (2018), and the grey dashed line refers to the best fit fromBehroozi et al. (2019).

Figure 8 .
Figure8.1500 Å luminosity functions (LF) for the various feedback models at  = 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.

Figure 9 .
Figure9.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  = 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 / 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.

Figure 11 .
Figure11.Luminosity-weighted mean escape fractions of SPICE galaxies as a function of stellar mass.Solid and dashed lines refers to dispersion-and rotation-dominated systems ( gas / gas ), respectively.Solid (hollow) stars represent single dispersion-(rotation-) dominated galaxies in bins with fewer than 3 galaxies.

Figure 12 .
Figure12.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.

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

Table 1 .
Hahn et al. 2020f et al. 2020)ted cosmological boxes, including the variable name, its adopted value and a short description.etal.2019;Rampfetal. 2020).The 2PPT approach involves perturbing particle masses with a distribution centered at  dm (seeHahn et al. 2020) which suppresses discreteness errors.RAMSES-RT evolves baryons from density and momentum fields which are discretized 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

Table 2 .
The time delay of SN is stochastically sampled, but all SN have the same energy hyper-sn 3-40 Myr 10 50 − 2 × 10 51 (SN) + 10 52 (HN)The time delay and energy of SN is stochastically sampled + metallicity-dependent HN rate 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.
3 Metal cooling is not currently modelled self-consistently with the local radiation, but seeKatz (2022)for a development in this direction.

Table
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 characteristics of galaxies in the three feedback models.