Astraeus I: The interplay between galaxy formation and reionization

We introduce a new self-consistent model of galaxy evolution and reionization, ASTRAEUS (semi-numerical rAdiative tranSfer coupling of galaxy formaTion and Reionization in N-body dArk mattEr simUlationS), which couples a state-of-the-art N-body simulation with the semi-analytical galaxy evolution DELPHI and the semi-numerical reionization scheme CIFOG. ASTRAEUS includes all the key processes of galaxy formation and evolution (including accretion, mergers, supernova and radiative feedback) and follows the time and spatial evolution of the ionized regions in the intergalactic medium (IGM). Importantly, it explores different radiative feedback models that cover the physically plausible parameter space, ranging from a weak and delayed to a strong and immediate reduction of gas mass available for star formation. From our simulation suite that covers the different radiative feedback prescriptions and ionization topologies, we find that radiative feedback continuously reduces star formation in galaxies with $M_h<10^{9.5}M_{\odot}$ upon local reionization; larger mass halos are unaffected even for the strongest and immediate radiative feedback cases during reionization. For this reason, the ionization topologies of different radiative feedback scenarios differ only on scales smaller than $1-2$Mpc, and significant deviations are only found when physical parameters (e.g. the escape fraction of ionizing photons) are altered based on galactic properties. Finally, we find observables (the ultra-violet luminosity function, stellar mass function, reionization histories and ionization topologies) are hardly affected by the choice of the used stellar population synthesis models that either model single stars or binaries.


INTRODUCTION
The Epoch of Reionization (EoR) represents the last major phase transition of hydrogen in the history of the Universe. Its beginning is marked by the appearance of the first stars and galaxies, whose Lyman continuum photons (with energy E > 13.6eV) gradually ionize the neutral hydrogen (H I ) in the intergalactic medium (IGM). The growing ionized bubbles around galaxies merge and expand until the IGM is completely ionized by z 6 (e.g. Fan et al. 2006;Becker et al. 2015). A rising number of high-redshift galaxy observations are providing us with increasing hints on the properties and numbers of star-formation driven ionizing sources a.k.hutter@rug.nl (e.g. Smit et al. 2014;Bouwens et al. 2015;Smit et al. 2018;Ouchi et al. 2018;De Barros et al. 2019;Maseda et al. 2020). These galaxy data-sets are complemented by (upper limits on) the 21cm emission from H I in the IGM during reionization obtained by experiments such as LOFAR 1 (Patil et al. 2017;Mertens et al. 2020), MWA 2 (Li et al. 2019; Barry et al. 2019) and PAPER 3 : (Kolopanis et al. 2019). Over the next decade, this 21cm data will be supplemented by that from state-of-the-art radio interferometers, such as the Square Kilometre Array (SKA; Carilli & Rawlings 2004) and the Hydrogen Epoch of Reionization Array (HERA;DeBoer et al. 2017), which are designed to measure the temporal and spatial evolution of the ionized regions, i.e. the reionization topology (e.g. Greig 2019;Seiler et al. 2019;Elbers & van de Weygaert 2019;Hutter et al. 2017Hutter et al. , 2020. Despite this progress, the reionization topology, the properties of the ionizing sources and the impact of reionization on the evolution of galaxy properties through radiative feedback effects remain key outstanding questions in the field of physical cosmology (for a review see e.g. Dayal & Ferrara 2018).
As the IGM becomes ionized, the associated ultra-violet background (UVB) photo-heats the gas in halos and the IGM to about ∼ 10 4 K. The higher temperature and rising pressure of the gas in a halo causes a fraction of the gas to photo-evaporate into the IGM (Barkana & Loeb 1999;Shapiro et al. 2004) and raises the Jeans mass for galaxy formation (reducing the amount of gas being accreted; Couchman & Rees 1986;Efstathiou 1992;Hoeft et al. 2006). Both mechanisms lead to a reduction of gas mass and the associated star formation rate, particularly in low-mass halos. However, modelling the impact of reionization feedback on galaxy formation remains challenging due to its complex dependence on halo mass and redshift, the patchiness and strength of the UVB and the redshift at which an assembling halo is first irradiated by the UVB (e.g. Gnedin 2000;Sobacchi & Mesinger 2013a).
Early works have studied the effects of radiative (photoheating) feedback on galaxies in cosmological hydrodynamical simulations by quantifying the loss of baryons in low-mass halos in the presence of a homogeneous UVB (e.g. Hoeft et al. 2006;Okamoto et al. 2008;Naoz et al. 2013). However, since reionization is a spatially inhomogeneous and temporal extended process, an increasing number of radiation hydrodynamical simulations have studied the impact of an inhomogeneous and evolving UVB on the galaxy population and found a reduction in the star formation rates in lowmass galaxies with halo mass M h 10 9 M (Gnedin 2000;Hasegawa & Semelin 2013;Gnedin & Kaurov 2014;Pawlik et al. 2015;Ocvirk et al. 2016Ocvirk et al. , 2018Katz et al. 2019;Wu et al. 2019). Most importantly, a number of such radiation hydrodynamical simulations show that the star-formationsuppressing effect of radiative increases with time, even after the Universe has been mostly ionized (Gnedin & Kaurov 2014;Ocvirk et al. 2016Ocvirk et al. , 2018Wu et al. 2019), which could be attributed to a decrease in self-shielding and a slower heating of the gas (Wu et al. 2019). The suppression of star formation is also found to be dependent on the environment, i.e. galaxies in over-dense regions that ionize earlier feature higher star formation rates which declines sharply after local reionization for low-mass halos with M h < 10 9 M ( Dawoodbhoy et al. 2018). Highlighting the interplay between galaxy formation and reionization, Wu et al. (2019) have shown that a stronger stellar feedback reduces the star formation within the galaxy and hence the UVB, weakening the strength of radiative feedback. In order to investigate the signatures of radiative feedback on the ionization topology, a number of works have combined N-body simulations with radiative transfer and used different suppression models for the ionizing emissivities of low-mass halos (e.g. Iliev et al. 2007Iliev et al. , 2012; Dixon et al. 2016). However, since these simulations do not contain a galaxy evolution model, the gas mass in halos below the local Jeans mass of the photo-heated IGM is instantaneously suppressed in ionized regions. Dif-ferent suppression models mostly affect the timing of reionization as compared to the ionization topology (Dixon et al. 2016).
In this paper, our aim is to quantify the effects of radiative feedback, both, on the underlying galaxy population as well the ionization topology during the EoR to answer questions including: When and which galaxies are most affected by radiative feedback? Is the patchiness of reionization imprinted in galaxy observables? How does radiative feedback impact high-redshift observables (including the UV luminosity function, stellar mass function and the redshift evolution of the star formation rate density and stellar mass density) and the 21cm signal from the neutral regions in the IGM? This naturally requires coupling galaxy formation and reionization using large volume simulations with a high-resolution to be able to study the ionization histories of galaxies based, both, on their masses as well as their location in the cosmic web. For example, in an inside-out reionization scenario, low-mass galaxies, can either be located in high-density regions that get ionized quite early on (therefore being strongly affected by UVB feedback) or in low-density regions that are ionized later (resulting in weak to no UVB feedback).
For this reason, we have built the astraeus (seminumerical rAdiative tranSfer coupling of galaxy formaTion and Reionization in N-body dArk mattEr simUlationS) framework that self-consistently couples a state-of-the-art N-body simulation (very small multi-dark; vsmd) with a semi-analytic model of galaxy formation (delphi; Dayal et al. 2014Dayal et al. , 2015Dayal et al. , 2017) and a semi-numerical reionization scheme (cifog; Hutter 2018). While similar approaches have been followed for meraxes (Mutch et al. 2016) and rsage (Seiler et al. 2019), these works have only focused on exploring the suppression of gas mass and star formation in low-mass halos based on 1D radiation hydrodynamical simulations (Sobacchi & Mesinger 2013a). In contrast, in this paper we explore different radiative feedback scenarios that range from a minimum one with gas loss via the characteristic mass approach outlined in Gnedin (2000) to a maximum one where the amount of gas is instantaneously reduced in halos with masses below the local Jeans mass of the ionized region; although similar in spirit to the work of Iliev et al. (2012) and Dixon et al. (2016), our model is an advancement on these works given it uses a much more sophisticated model for galaxy formation and the associated ionizing emissivity. Besides its key strength of supporting multiple radiative feedback models, astraeus comprises (1) a large volume and high-resolution N-body simulation that allows us to simultaneously explore the large-scale reionization topology whilst resolving sources down to the atomic cooling mass at z ∼ 6, (2) a galaxy formation model that uses only three free parameters with feedback being linked to the underlying halo potential, and (3) supports multiple models for the ionizing escape fraction fesc that enable us to cover the physically plausible range of reionization scenarios.
The paper is structured as follows. In Section 2 we describe the underlying N-body simulation and the theoretical galaxy and reionization model, as well as our different models of radiative feedback. In Section 3 we compare our results to observational constraints, such as the luminosity and stellar mass functions, and the Thomson optical depth for reionization. We then use our different models for radia-tive feedback to investigate how the strength and timing of the suppression of star formation in a galaxy depends on its gravitational potential and local reionization in Section 4, how radiative feedback affects the ionization topology and thus the power spectrum of the 21cm signal in Section 5, and whether assuming a different stellar population synthesis models affects any observables in Section 6. We conclude in Section 7.

THE THEORETICAL MODEL
In this Section, we describe our self-consistent, seminumerical model that couples high-redshift galaxy formation and reionization, astraeus 4 . Using the evolving DM density distribution from a high-resolution N-body simulation (Section 2.1), astraeus couples an enhanced version of the semi-analytic galaxy evolution model delphi (Dayal et al. 2014, Section 2.2) to the semi-numerical reionization code cifog 5 (Hutter 2018, Section 2.3). The key novelty of astraeus is that it allows us to explore a wide range of scenarios for the interplay between galaxy formation and reionization using a minimum number of mass-and redshiftindependent free parameters.

N-body simulation
In this work we use the high resolution Very Small Mul-tiDark Planck (vsmdpl) N-body simulation, performed as part of the multidark simulation project 6 . This new simulation, with a box size of 160h −1 Mpc, was run with the same number of particles (3840 3 ) and using the same gadget-2 Tree+PM N-body code (Springel 2005) as in the other Multidark simulations described in Klypin et al. (2016). We also used the same cosmological parameters to set up initial conditions, namely [ΩΛ, Ωm, Ω b , h, ns, σ8] = [0.69, 0.31, 0.048, 0.68, 0.96, 0.83]. The Zeldovich approximation was used to produce the particle positions and velocities at an initial redshift of z = 150. The mass per dark matter particle is 6.2×10 6 h −1 M and the equivalent Plummer's gravitational softening was set to 2h −1 physical kpc at z > 1. A total of 150 different snapshots of the simulation, equally spaced in expansion factor, were stored from z = 25 till z = 0. The Rockstar phase-space halo finder (Behroozi et al. 2013a) was used to identify all halos and subhalos in each of the 150 snapshots, down to a minimum of 20 particles per halo resulting in a minimum resolved halo mass of 1.24 × 10 8 h −1 M . In addition, merger trees from the rockstar halo catalogues were computed using the consistent trees (Behroozi et al. 2013b) method 7 . While the vertical merger trees obtained from consistent 4 astraeus can be built from the source code publicly available under https://github.com/annehutter/astraeus. astraeus includes a new implementation of delphi and uses the cifog library. 5 cifog is publicly available under https://github.com/ annehutter/grid-model. 6 See www.cosmosim.org for further information about the Multidark suite of simulations and access to the simulations database. 7 We caution the reader that although the rockstar merger trees might be incomplete below 10 9 M , we are still complete in terms of the halo mass function down to 10 8.26 M .
trees are well suited to follow the evolution history of a single galaxy, i.e. following the evolution of the progenitors of a galaxy, they do not track the galaxy population on a redshift-step-by-redshift-step basis as required for reionization. In order to use astraeus as a semi-analytic galaxy formation code run on a tree-branch-by-tree-branch basis (i.e. fully vertical) as in sage (Croton et al. 2016) or delphi (Dayal et al. 2014) or on a redshift-by-redshift basis (i.e. fully horizontal) as in meraxes (Mutch et al. 2016), we re-sort the consistent tree outputs as follows. We keep the merger-tree-by-merger-tree order but each merger tree is sorted by redshift (horizontally sorted), and we refer to this as locally-horizontally sorted. We refrain from generating fully horizontal outputs on a galaxy-by-galaxy basis, as such an order would impede the possibility of following the evolution of a single galaxy easily and limit the flexibility of astraeus to be used for non-reionization galaxy studies in the future.

Semi-analytic galaxy modelling
Our semi-analytic galaxy formation model includes all the key baryonic processes of gas accretion, gas and stellar mass being brought in by mergers, star formation and the associated supernovae (SN) feedback and radiative feedback from reionization. At each time step, these are coupled to the merger-and accretion-driven growth of the dark matter halos obtained from the N-body simulations as explained in this section. Throughout this work, we use a Salpeter (Salpeter 1955) initial mass function (IMF) with a slope of α = 2.35 between 0.1 − 100 M .

Gas accretion and mergers
There are two ways in which a galaxy can build up its gas content: through smooth accretion from the IGM and through mergers in the case that a galaxy has progenitors. On the one hand, at the beginning of a time step, a galaxy of halo mass M h (z) that has no progenitors, can, in principle, smoothly accrete an initial gas mass, M i g (z), corresponding to the cosmological baryon-to-dark matter ratio such that M i g (z) = (Ω b /Ωm)M h (z). However, reionization feedback can reduce the initial gas mass by photo-evaporating gas out of the potential. In this case M i g (z) = fg(Ω b /Ωm)M h (z) where fg is the gas fraction not affected by reionization as explain in Sec. 2.3 that follows.
On the other hand, galaxies that have (say Np) progenitors can also gain gas through mergers. In this case, the merged gas mass can be expressed as where Mg,p(z + ∆z) is the final gas mass of the previous time step brought in by the merging progenitors of halo mass M h,p (z + ∆z). The accreted gas mass in this case is given by where we have made the reasonable assumption that accretion of halo mass from the IGM drags in a cosmological fraction of gas mass. Accounting for the impact of reionization feedback, the initial gas mass can be expressed as

Star formation and stellar mass assembly
We assume that at a given time step this initial gas mass, M i g , can form stars with an effective efficiency (f eff ) which is the minimum between that required to eject the rest of the gas from the halo potential (f ej ) and quench star formation and an upper limit (f ∼ 1 − 3%) such that f eff = min f , f ej ; details of the calculation of f ej follow in Section 2.2.3. The newly formed stellar mass at any time step can then be expressed as Physically, the effective efficiency can be thought of as f eff = fs/ts i.e. a fraction (fs) of the gas mass that can form stars over a timescale ts. Given that f eff is linked to the underlying halo potential, our model results in low-mass galaxies (M h < ∼ 10 9.3 M at z = 5) being star formation efficiency limited with f eff = f ej , while larger mass halos form stars with a constant efficiency f eff = f (see also Dayal et al. 2014).
In addition, stellar mass can also be brought in by merging progenitors (M ,p) such that resulting in a total stellar mass Star formation in galaxies has two key physical effects: firstly, at the end of their life, high-mass stars explode as Type II supernovae (SNII) which can eject gas mass (say, M ej g ) from the galaxy. Secondly, star formation provides H I ionizing photons; the fraction of these photons that can escape into the IGM (fesc) contribute to reionizing and heating the IGM as detailed in Sec. 2.3.

Supernova feedback
The explosion of high-mass stars as SNII injects thermal and kinetic energy into the interstellar medium (ISM) that can heat and eject gas from the galaxy, respectively. In our model, we only consider the latter effect that can eject gas out of the galactic environment. We assume each SNII to produce an energy equal to E51 = 10 51 erg of which a fraction (fw) couples to the gas and drives the winds. In this work, we implement a "delayed SN feedback" scheme (see also Mutch et al. 2016;Seiler et al. 2019) that accounts for the mass-dependent (MSN) lifetimes (tSN) of stars before they explode as SNII. We use the MSN − tSN relation found by Padovani & Matteucci (1993) In this case, stars of MSN = 8 (100) M explode as SNII 28.6 (3.23) Myr after star formation starts. astraeus, as many other semi-analytic galaxy evolution models, is based on merger trees that are discrete in time. Hence, the length of the time steps is key when deriving the total SNII energy. In case a time step exceeds 30 Myrs, all high mass stars formed explode as SNII within that time step and SN feedback is "instantaneous". However, as the time steps become increasingly shorter than 30 Myrs, SNII feedback spans over multiple time steps, i.e., is "delayed". The instantaneous SNII feedback scenario can therefore be regarded as a special case of the more general delayed SNII feedback. We now describe our formalism for delayed SNII feedback. At any given redshift, the SNII energy that can couple to gas can be expressed as Here, for a given halo at redshift z, the first term on the right hand side represents the SNII explosions at z from all the stars that formed in the progenitors of that halo, and the second term accounts for SNII from the newly formed stars at z. Further, M new ,p (zj) and νj are the newly formed stellar mass and the fraction of stars that explode as SN in time step j, between zj and zj−1, respectively and Nj is the number of simulation snapshots until and including z. Finally, νz is the fraction of the newly formed stars in the current time-step, M new (z), that explode as SNII. Using the assumed Salpeter IMF, the fraction of stars formed in the time interval [t(zj−1), t(zj)] that explode as supernovae at step z can be calculated as Here MSN,j is the mass of stars that would explode as SN after tSN = t(z) − t(zj) according to equation 7 with t tj, while M star,low = 0.1 M and M star,high = 100 M . In order to derive the amount of gas being ejected from the galaxy, we equate the SNII energy (see equation 8) to the energy required to eject a gas mass (equal to M ej g ) from a galaxy: where ve is the ejection velocity that is related to the rotational velocity of the halo as vc = ve/ √ 2. This implies that the fraction of gas that needs to be converted into stars to eject the rest of the gas from the galaxy corresponds to where at maximum all the gas that is left after star formation can be ejected. In case of instantaneous SN feedback equations 10 and with νz = 0.0077 M −1 for the assumed Salpeter IMF.
In the remainder of this paper we use the delayed SN feedback scheme. As noted, the snapshots of the N-body simulation scale as the logarithm of the scale-factor, resulting in increasingly longer time-steps with decreasing redshift z. Hence, while at z 9 the delayed SN feedback scheme differs significantly from the instantaneous one, these schemes become increasingly similar with decreasing z until there is effectively no difference at z < ∼ 6.

Resulting output of UV and H I ionizing photons
We calculate the spectrum of each galaxy, ξ(ν, t), by convolving its star formation history with the starburst spectrum, ξSP(ν, t), obtained from two stellar population synthesis models (SPS): starburst99 (Leitherer et al. 1999) and bpass that accounts for binaries (Eldridge et al. 2017).
where M new (z ) (or M new (tj)) is the newly formed stellar mass at the time step z (time step tj). This newly formed stellar mass is assumed to form from star formation uniformly distributed over the entire time step. The factor f lin accounts for this and is calculated as where tj+1 and tj are the beginning and end times of the time step with t tj+1.
The intrinsic spectrum of a stellar population sensitively depends on its age (t) and metallicity (Z). In the interest of simplicity, in this paper, we assume all stellar populations to have a stellar metallicity of Z = 0.05 Z and defer a full metallicity calculation to a future paper. In this paper, we focus on two key spectral quantities: (1) the number of H I ionizing photons (λ < 912Å in the rest-frame) that are required to understand the reionization of the IGM, and (2) the UV luminosity (rest-frame 1250 − 1500Å) to validate our model against observed Lyman Break Galaxy (LBG) data.
The intrinsic UV luminosity, Lν [erg s −1 Hz −1 M −1 ], is quite similar in both the starburst99 and bpass models and from the point in time when the stellar population was formed it evolves with time as (18) The total UV luminosity or ionizing photon output over any star formation history can be derived by using ξSP = Lν and ξSP =Q in equation 14, respectively.

The reionization model
Most of the ionizing photons produced by a source (calculated from its star formation history as explained above) are absorbed within the interstellar medium with only a fraction (fesc) escaping and ionizing the IGM. This escaping rate of ionizing photons (the ionizing emissivity) can be expressed aṡ withQ(z) being given by combining equations 14 and 17 (18) for starburst99 (bpass). For the majority of reionization scenarios that we consider in this paper, the ionizing escape fraction fesc is assumed to be constant for all galaxies at all redshifts. However, we also explore a scenario where fesc is coupled to the fraction of gas that is ejected from the galaxy into the IGM. This is supported by a number of simulations that find that SN explosions create under-densities through which ionizing photons can escape (Wise & Cen 2009;Wise et al. 2014;Kimm & Cen 2014;Kimm et al. 2017;Xu et al. 2016), implying that the ionizing escape fraction fesc increases as a larger fraction of gas is pushed into outflows. In this case we model the ionizing escape fraction fesc as, where f 0 esc is a free parameter that can be tuned to adjust the timing of reionization 8 . This ansatz results in a very high escape fractions for low-mass ( < ∼ 10 9.5 M ) galaxies where f eff = f ej . As the gravitational potential of the galaxy deepens, the SN explosions of the stars that are forming can not eject all gas from the galaxy (f eff f ej ) and fesc drops down to a few percent for M h ∼ 10 11 M halos and assuming f 0 esc = 1. These escaping ionizing photons both reionize and heat the IGM. The amount of heating, naturally, critically depends on the energy of the ionizing photons (i.e. the hardness of the source spectrum). For star-forming galaxies the ionized IGM can be heated up to ∼ 10 4 K (e.g. Schaye et al. 2000), which has two important effects: firstly, as the gas residing in halos heats up, the higher pressure causes a fraction 8 For instantaneous SN feedback, equation 20 can be expressed c 0000 RAS, MNRAS 000, 000-000 6 Hutter et al.
of it to photo-evaporate into the IGM, reducing the amount available for star formation (Barkana & Loeb 1999;Shapiro et al. 2004). Secondly, since a higher IGM temperature corresponds to a higher Jeans mass, the minimum mass for galaxy formation increases, thereby reducing the amount of gas being accreted by the galaxy (Couchman & Rees 1986;Efstathiou 1992;Hoeft et al. 2006). These mechanisms lead to a reduction of gas mass, particularly in low-mass halos where the gravitational potential is not deep enough to compensate for the increased pressure of the heated gas. While the rise in gas temperature occurs quasi instantaneously, the gas pressure adjusts over the dynamical time scale of the galaxy (Gnedin 2000), which leads to a time delay between the time of reionization and the onset of the gas (and star formation) suppressing effect of radiative feedback. However, modelling these radiative feedback processes remains challenging due to their complex dependence on the halo mass and redshift, the patchiness and strength of the UVB and the redshift at which an assembling halo is first irradiated by the UV background (Gnedin 2000;Sobacchi & Mesinger 2013b).
In this work, we explore a wide range of physically plausible radiative feedback models in order to study their impact on both the galaxy population and progress of reionization whilst ensuring agreement with all available (galaxy and reionization) data sets. It is essential to account for the patchiness of reionization: for example, the impact of radiative feedback might be more severe on galaxies forming in an over-dense region reionized early-on as compared to those forming in an under-dense region reionized later. As noted before, in order to simulate the galaxy formationreionization interplay, we couple galaxy formation (simulated through the delphi semi-analytic galaxy evolution model) with a semi-numerical reionization code (cifog; 9 Hutter 2018) in a self-consistent manner.
cifog is a MPI-parallelised, semi-numerical reionization code that computes the time-and spatial-evolution of ionized regions in the IGM. Here we provide a brief overview and refer the interested reader to Hutter (2018) for details. Essentially, cifog follows the approach outlined in Furlanetto et al. (2004) where a spherical region is considered to be ionized if the cumulative number of ionizing photons (Nion) emitted exceeds the cumulative number of absorption events (N abs ), and neutral otherwise. Starting with large radii and decreasing the size of the region by reducing the sphere radius R, the central cell of the spherical region is considered ionized if Here,Ṅion(z) is the ionizing emissivity of a galaxy i at redshift z located within the sphere of radius R and N gal (R) is the number of galaxies within that sphere; R indicates that the quantity is averaged over a sphere with radius R. Further, nH,0, V cell andṄrec(z) are the hydrogen density of the cell at z = 0, the comoving volume of the cell and the re-9 https://github.com/annehutter/grid-model combination rate at z, respectively. Applying the ionization criterion in large enough regions ensures that the radiation from neighbouring sources is accounted for. cifog derives the residual H I fraction and recombination rate of each cell from the local gas density and photoionization rate. The code supports two models for the spatially-dependent photoionization rate: one that is based on the mean free path given by the size of the ionized regions (mean free path approach), and one that is based on the flux of the ionizing sources (flux based approach). In this work we utilise the flux based approach 10 .
In this work, we use the density fields that have been obtained by mapping the DM particles on to a 2048 3 grid using the cloud-in-cell (CIC) algorithm. cifog then runs on 512 3 grids that have been obtained by reducing the 2048 3 grids. Furthermore, we note that the time and spatial evolution of the larger ionized regions computed with cifog are hardly affected by the resolution of the grid. However, ionized regions smaller than a grid cell are not spatially resolved, which leads to the cell being reionized when the volume of the cell is ionized. A better resolution would resolve those ionized regions and provide more accurate times of reionization (zreion), particularly around low-mass galaxies.
delphi and cifog are coupled in a self-consistent manner using the following approach at each time step: (i) the ionizing emissivity of each galaxy at zj is computed from its star formation history.
(ii) the ionizing emissivities of the galaxies are fed into cifog. Accounting for the location of the galaxies in the simulated large-scale structure, cifog computes the time evolution of the ionized regions in the IGM on a 512 3 grid from zj to zj+1.
(iii) In the subsequent time step (zj+1), we identify each galaxy whose cell was reionized in previous time steps z > zj+1. For galaxies lying in reionized regions, we track their redshift of reionization zreion and the incident photoionization rate at zreion, and calculate the fraction of gas mass they can retain after radiative feedback (fg) using the prescriptions detailed in Sec. 2.4 that follows; galaxies in neutral regions are naturally unaffected by reionization feedback.

Resulting characteristic masses of supernova and radiative feedback processes
As described above, both supernova and radiative feedback affect the gas content of galaxies, with the feedback efficiency generally decreasing as the gravitational potential of the host halo increases. This results in a "characteristic" mass that can be associated with each feedback processthis is defined as the mass at which the halo can still hold on to some and half of its original gas mass for supernova and radiative feedback, respectively. We now discuss the characteristic masses for SN feedback and the different radiative feedback models that have been implemented in astraeus and are summarised in Table 1. . From left to right the masses are shown for our radiative feedback models Minimum, Heating with T 0 = 2 × 10 4 K and Mc = M F , Photoionization with Γ HI = 10 −12.3 s −1 , and Heating with T 0 = 4×10 4 K and Mc = 8M F . Red, blue and green solid lines correspond to the characteristic masses due radiative feedback when the region has been reionized at z = 8, 11 and 14, respectively. In the Minimum model the filtering mass is independent of reionization and is shown by the black solid line. In the third panel, coloured, dotted lines show the characteristic masses for a Heating model with T 0 = 4 × 10 4 K and Mc = M F . Gray, dash-dotted lines show the characteristic masses for SN feedback, while grey dashed and dash-dotted lines the Jeans mass at mean density and virial over-density for the temperature indicated. The grey shaded area marks the halo masses that are not resolved in our simulation.

Characteristic mass for supernova feedback
Given that SN feedback can eject gas from every star forming galaxy, we define its characteristic mass (M SN c ) as the minimum halo mass which can hold onto a non-zero value of its gas mass after star formation. In the case of instantaneous SN feedback, this can be calculated, by equating the supernova energy coupling to gas to the binding energy of the gas left over after star formation, as The redshift evolution of M SN c is shown as a reference in all panels by the dash-dotted lines in Fig. 1. As seen, SN feedback can eject all of the gas mass in halos less massive than 10 8.3 M at z ∼ 16. As the matter density of the Universe decreases with cosmic time, the gravitational potential corresponding to a given halo mass becomes shallower, leading to an increase in the supernova characteristic mass to M SN c ∼ 10 9 M by z ∼ 5.

Atomic cooling mass (Minimum feedback model)
This is the weakest of our radiative feedback models. Here, ionized IGM gas is assumed to be heated to T = 10 4 K via photo-heating. Only halos massive enough to have virial temperatures exceeding 10 4 K can maintain all of their gas; lower mass halos are assumed to be completely gas-free. The gas fraction left after radiative feedback (fg) is obtained by comparing the halo mass to the (cooling) mass within the virial radius at the critical over-density for collapse, ∆cρc 18π 2 ρc (e.g. Barkana et al. 2001 using Tvir = 10 4 K. Then fg is calculated to be From the first panel in Fig. 1, we see that the characteristic mass for SN feedback always exceeds that for the atomic cooling mass. As a result, our Minimum radiative feedback model shows no impact of radiative feedback on either galaxy formation or reionization.

Temperature-dependent characteristic masses (Heating models)
Gnedin & Hui (1998) introduced a filtering scale, kF , below which baryon fluctuations are suppressed as a consequence of reionization heating. This filtering scale can be linked to a filtering mass (Gnedin 2000;Okamoto et al. 2008;Naoz et al. 2013) where ρ is the average total mass density of the Universe and a the scale factor. Applying the filtering scale for baryons to the continuity equation that describes the linear evolution of perturbations in the dark matter-baryon fluid, it has been shown that whilst the filtering scale is related to the Jeans scale as a function of time, at a given time those two scales are unrelated (Gnedin & Hui 1998); determining kF at a given time therefore requires knowing the evolution of the Jeans scale up to that time. This is due to the fact that the response of the gas density distribution to an immediate temperature change occurs on the dynamical timescale. The filtering mass can be calculated as (at z > 2) where MJ 0 is the Jeans mass at T = 10 4 K and z = 0, and T4(a) is the evolving baryonic temperature in units of 10 4 K.
The Jeans mass at redshift z depends on the Jeans scale kJ = a c −1 s (4πGρ) 1/2 and the linear-theory sound speed cs = 5kBT [3µmp] −1 , and can be calculated as . Furthermore, we model the redshift evolution of the baryonic temperature as if areion a.
These three terms correspond to the epoch after recombination (arec = 1/1100) where gas is still coupled to the cosmic background radiation by Compton heating, the epoch after decoupling (a dec = 1/251) when gas cools adiabatically, and the epoch of reionization and subsequent cooling (Hui & Gnedin 1997), respectively. The IGM temperature, T0, is a free parameter. Although a number of hydrodynamical simulations of galaxy populations with a homogeneous UV background (Gnedin 2000;Okamoto et al. 2008;Naoz et al. 2013) have shown that MF can be related to the characteristic mass Mc (the halo mass that on average retains 50% of its gas mass), the exact relation remains debated: while Gnedin (2000) obtain Mc 8MF , other works yield Mc MF (Naoz et al. 2013). From the characteristic mass, the gas fraction maintained by the galaxy can be found following the fitting formula provided in Gnedin (2000) such that In this paper, we explore three scenarios for such a Heating model: the weakest (Weak Heating) and the strongest (Strong Heating) have been chosen to bracket the physically plausible range of radiative feedback for this model: (i) Weak Heating: Here we assume the reionized IGM is heated to T0 = 2 × 10 4 K and the characteristic mass for radiative feedback is equal to the filtering mass, i.e., Mc = MF . From the second panel in Fig. 1, we can see that only galaxies reionized very early-on (i.e. at z 14) are affected more by radiative feedback as compared to SN feedback. Galaxies reionized later on are only affected if their halo masses are less than ∼ 10 8−9 M .
(ii) Early Heating: In this model we use an IGM temperature of T0 = 4 × 10 4 K and Mc = MF , resulting in similar characteristic masses as the Photoionization model described in the next Section. However, this model is designed to explore the extent to which the impact of radiative feedback can be enhanced for low-mass galaxies (M h 10 9.5 M ) by reionizing them earlier. In order to maximise this effect, whilst remaining in agreement with the Planck optical depth measurements, we assume the ionizing escape fraction fesc to scale with the ejected gas fraction, resulting in a decreasing fesc with halo mass. For identical fesc, f and fw, such a Heating model produces very similar results to the Photoionization model discussed below (c.f. dotted to solid coloured lines in third panel in Fig. 1).
(iii) Strong Heating: As in the Early Heating model, we assume an IGM temperature of T0 = 4 × 10 4 K. However, in order to increase the impact of radiative feedback, we assume the radiative feedback characteristic mass to be 8 times the filtering mass, i.e., Mc = 8MF . From the last panel in Fig.  1, we see that even galaxies reionized later (e.g. z 8) exceed the characteristic mass for SN feedback. Indeed, at z 6, even galaxies with halo masses up to ∼ 10 10 M show suppression of star formation in this model. From the point in time, when a galaxy becomes reionized, radiative feedback dominates over SN feedback for a time period corresponding to ∆z 1 − 1.5.

Photoionization rate dependent characteristic mass (Photoionization model)
Using 1D radiation hydrodynamical simulations and assuming an IGM temperature of T = 0 (10 4 ) K in neutral (ionized) regions, Sobacchi & Mesinger (2013b) have derived the following ansatz for the critical mass using the filtering mass approach proposed in Gnedin (2000)'s: This is motivated by the fact that, inserting their temperature relation (equation 30) into equation 26, we can see that only MJ,0 (or kJ ) is dependent on the temperature T0, which again can be expressed in terms of the ionizing background J21. Quantitatively, the critical mass is found to be (Sobacchi & Mesinger 2013b) with best-fit values of M0 = 2.8 × 10 9 M , a = 0.17, b = −2.1, c = 2, d = 2.5, and J21 = (ΓHI/10 −12 )s −1 where ΓHI is the photoionization rate at the location of the galaxy at the time its environment was reionization. In this case, the the fraction of gas that is retained by a halo can be expressed as We note that the photoionization rate ΓHI is a proxy for the IGM temperature. We find that a photoionization rate  for starburst99 or bpass, respectively) are our model free parameters that are tuned to simultaneously reproduce all high-redshift galaxy and reionization data sets. These model parameters have similar and even identical values, since our radiative feedback models affect only low-mass and faint galaxies where observational constraints are sparse. Extreme models that alter the ionizing emissivities of galaxies, either through suppression of star formation (Jeans Mass) or an fesc depending on the fraction of gas ejected from the galaxy (Early Heating) show higher fesc values.
of ΓHI = 10 −12.3 s −1 corresponds to a temperature of T0 4 × 10 4 K in the Heating model. 11 In this Photoionization model, as shown (by the solid coloured lines) in the third panel in Fig. 1, only galaxies reionized early-on, i.e. at z > ∼ 10, will experience sufficient radiative feedback. Galaxies reionized later are only affected by radiative feedback when they have a smaller gravitational potential, i.e. are less massive than ∼ 10 9 M at z 6. The dotted coloured lines in the third panel in Fig. 1 show results for the corresponding Heating model with an IGM temperature of T0 = 4 × 10 4 K and Mc = MF .

Jeans mass (Maximum feedback model)
The strongest of our radiative feedback models is the Jeans mass model. Here, the gas in the IGM is assumed to be heated to T0 = 4 × 10 4 K via photoheating upon ionization. However, the rise in temperature is assumed to translate immediately into a lower gas density and hence a higher Jeans mass at the virial over-density. The fraction of gas fg that is maintained by a galaxy in an ionized region is given by where MJ is determined by equation 28. This Jeans mass at each redshift is shown (as the dotted grey line) in the fourth panel in Fig. 1. As soon as the cell hosting a galaxy is reionized, galaxies with halo mass M h < ∼ MJ are immediately affected by radiative feedback. 11 The relation between Γ HI and T 0 has been derived from analysing the Γ HI values in over-dense cells at the time of reionization as well as comparing observables such as the UV luminosity and stellar mass functions.

BASELINING THE MODEL AGAINST OBSERVED DATA-SETS
We tune the three free mass-and redshift-independent parameters of our framework (f , fw and fesc) for each radiative feedback model by simultaneously matching to a number of galaxy observables (the UV luminosity functions at z = 5 − 10, the stellar mass functions at z = 5 − 10, the redshift evolution of the stellar mass and star formation rate) and reionization data-sets (constraints on the ionization history inferred using quasars, Lyman-α emitters, Gamma Ray Bursts and the integrated electron scattering optical depth). The best fit values of the free parameters for each radiative feedback model are listed in Table 1. We note that our free parameters should be thought of as the "observed" values, since we calibrate our model to observations without correcting for effects such as dust attenuation.

Redshift evolution of the Ultra-violet luminosity function
For all galaxies in our simulation, we calculate the UV luminosities at 1500Å at z = 5 − 10 from their entire star formation histories (SFH) by inserting equation 16 for ξSP into equation 14. In Fig. 2, we show the UV luminosity functions (UV LFs) for all our radiative feedback models. We start by noting that while our model results are in broad agreement with the observed UV LF at z ∼ 5 − 10, they slightly over-predict the number density of bright galaxies at z < ∼ 6. This is probably due to the fact that our model does not account for the increasing dust attenuation expected for massive galaxies with cosmic time.
In our model the bright end of the UV LF is determined by the threshold star formation efficiency f * , while the faint end of the UV LF (MUV −15) is shaped by a combination of supernova (fw) and radiative feedback (Mc). While the evolution of the bright end is driven by a genuine luminosity evolution as galaxies grow in (halo, gas and stellar) mass through mergers and gas accretion, the evolution of the faint  (Atek et al. 2015(Atek et al. , 2018Bouwens et al. 2015Bouwens et al. , 2016Bouwens et al. , 2017Bowler et al. 2014Bowler et al. , 2015Calvi et al. 2016;Castellano et al. 2010a,b;Finkelstein et al. 2015;Ishigaki et al. 2018;Livermore et al. 2017;McLeod et al. 2015McLeod et al. , 2016McLure et al. 2009McLure et al. , 2013Oesch et al. 2013Oesch et al. , 2014Ouchi et al. 2009;Schenker et al. 2013;Schmidt et al. 2014;Tilvi et al. 2013;van der Burg et al. 2010;Willott et al. 2013;Zheng et al. 2012). end is more complex due to the stochastic star formation in low-mass halos induced by SN feedback and by radiative feedback. Indeed, at the faint end, the UV LF evolution involves a combination of positive and negative luminosity evolution (as low-mass galaxies brighten and fade) and a positive and negative density evolution (as new low-mass galaxies form are consumed by merging) as pointed out in previous works (e.g Dayal et al. 2013).
In order to assess the role of SN feedback in shaping the faint-end of the UV LF, we start by discussing the key features of the UV LFs in our Minimum radiative feedback model (see black solid lines in Fig. 2). In this model, the characteristic mass for radiative feedback only exceeds the halo resolution mass (1.2 × 10 8 h −1 M ) at z < ∼ 5.8 (see Fig.  1) and always remains below the characteristic mass for SN feedback. This results in SN feedback dominating over radiative feedback at all redshifts in this model. As cosmic time proceeds and the density contrast in the Universe increases, the turn over or peak of the UV LF broadens and shifts to fainter UV luminosities. In order to explain these trends, we start by examining the characteristic UV luminosity of galaxies in newly-formed halos at the resolution limit of the simulation (M h 1.2 × 10 8 h −1 M ) indicated by the grey dotted lines in Fig. 2. Initially, these galaxies are gas-rich and have a burst of star formation with the exact UV luminosity depending on the SN feedback efficiency (fw) and redshift; e.g. newly formed halos have a UV magnitude corresponding to MUV −14 at z 10 which increases to MUV −12.5 by z 5). The complete loss of (SN-driven) gas mass, which can not be compensated by accretion in a subsequent time step, results in almost no new star formation after the initial burst. This results in a continual decrease in the UV luminosity of such low-mass galaxies. Indeed, as such galaxies age, the UV luminosity only drops. As a result, fainter UV luminosities can be reached in gas-poor low-mass galaxies as they age, explaining the broadening of the faint end turn over towards fainter UV luminosities with decreasing redshift.
Adding radiative feedback, we find the star formation in low-mass galaxies, that are located in ionized regions, to be increasingly suppressed as reionization proceeds. Hence, the faint end of the UV LF becomes flatter and pushes the turn over at the faint end to lower UV luminosities from z = 10 to z = 5 (see Photoionization and all Heating models). These trends become stronger towards lower redshifts due to two reasons: firstly, as reionization proceeds (and more regions become ionized), a larger fraction of low-mass halos is affected by radiative feedback. Secondly, the effect of radiative feedback increases as the time when the region became ionized lies further in the past (except the Jeans Mass model); the ionized and heated gas in the galaxy has had more time to adjust to its lower "equilibrium" density. The higher is the temperature that the IGM is heated up upon ionization (i.e. the higher the ratio between filtering and Jeans mass at virial over-density), the lower is this "equilibrium" density and hence the stronger is the effect of the radiative feedback. Indeed, in Fig. 2 we see that the turn over at the faint end moves to fainter UV luminosities and that its slope flattens as radiative feedback becomes stronger, i.e. going from the Minimum and Weak Heating models to the Photoionization model to the Strong Heating model (and Jeans Mass model, however we caution the reader that the Jeans Mass model assumes an instantaneous drop in gas density and is discussed in the following).
From Fig. 2, we also find that a stronger feedback goes in hand with a higher UV luminosity below which radiative feedback suppresses star formation; in the following we refer to this characteristic suppressed UV luminosity as MUV,s. At each redshift, the value of MUV,s corresponds to the UV luminosity of a halo whose mass equals the characteristic mass of the respective radiative feedback model for assuming a reionization redshift of zreion 15, i.e. when the first progenitors of these MUV,s galaxies reionized their environ-ment. With decreasing redshift the radiative feedback characteristic mass increases and correspondingly MUV,s shifts to brighter UV luminosities. For example, MUV,s shifts from −13.5 at z = 7 to −14.5 at z = 5 for the Photoionization model, and from −17.5 at z = 7 to −18 at z = 5 for the Strong Heating model. In the Heating models with Mc = k × MF , the characteristic mass approaches k-times the Jeans mass MJ (z) at 1 + z < (1 + zreion)/3.2 where zreion is the reionization redshift. This convergence towards the Jeans mass at virial over-density is also noticeable when comparing the Strong Heating model with an effective temperature of T0 = 16 × 10 4 K to the Jeans Mass model with T0 = 4 × 10 4 K. While MUV,s for the Jeans Mass model corresponds to the UV luminosity of a halo with Jeans mass at all times, MUV,s for the Strong Heating model approaches the Jeans mass at virial over-density only at later times when the gas density has had enough time to adjust to the change in temperature upon reionization.

Stellar mass functions
We now discuss the stellar mass functions (SMFs) at z = 5 − 10 obtained for the different feedback models used in this work, as shown in Fig. 3. We find that the overall normalisation of the SMF increases with decreasing red-shift, as new galaxies form and existing ones grow. As expected, the SMF at the massive end (M h 10 9.5 M or M 10 7 M ) increases with decreasing redshift as these galaxies assemble mass through star formation (at the maximum threshold efficiency f * ) and mergers. In contrast, the low-mass end (M 10 7 M ) shows a decrease with redshift: this is due to a combination of SN and radiative feedback driven decrease in the star formation efficiency, as well as low-mass galaxies moving into higher mass bins with time. When accounting only for supernova feedback (i.e. the Minimum model at z 5.8), the change in the lowmass slope is moderate. Analogous to the UV LF, adding radiative feedback enhances the flattening of the low-mass slope with decreasing redshift as a larger number of lowmass halos are feedback suppressed. The flattening increases as radiative feedback becomes stronger (i.e. the characteristic mass Mc increases), going from the Photoionization to the Strong Heating model), i.e. as low-mass galaxies become increasingly gas-poor and less efficient in forming stars. For the strongest radiative feedback cases, the low-mass slope can become so shallow that there are fewer galaxies with M f eff (Ω b /Ωm) M h 10 5.3 M at low redshifts (z 5 − 6) than at high redshifts (z 7).
Further, similar to the UV LF, the maximum stellar mass M ,s suppressed by radiative feedback, i.e. the mass at which the low-mass slope deviates from that of the Minimum model, increases with the radiative feedback strength. For example, M ,s 10 6 M for the Photoionization and ∼ 10 8 M for the Strong Heating models at z 7, respectively. Analogous to MUV,s, M ,s traces the highest possible radiative feedback characteristic mass at a given redshift, which naturally increases with increasing radiative feedback strength. Again the Jeans Mass model constitutes an exception: here M ,s hardly evolves with cosmic time and corresponds consistently to the stellar mass in a halo of Jeans mass MJ (z). This difference becomes particularly obvious at z = 5 when the radiative feedback characteristic mass Mc of the Strong Heating model exceeds that of the Jeans Mass model. This results in the SMFs (for MUV −13 and MUV −15 in Fig. 3) turning over at higher masses whilst showing a weaker suppression of stellar mass at the low-mass end (M 10 7 M ). Our model results of the SMFs at z = 5−10 are in agreement within the uncertainties of the observations. However, from Fig. 3 we can see that our model SMFs tend to overpredict the low-mass and underpredict the high-mass end of the observed SMF, despite our UV LFs at z = 5 − 10 being in good agreement with the observations (although slightly overpredicting the bright end as noted in Sec. 3.1). There may be multiple reasons for these trends: firstly, from the modelling side, we do not account for dust attenuation when computing the UV luminosities -this can underestimate the star formation efficiency f and the resulting stellar masses. Secondly, the observational data points are subject to modelling uncertainties when deriving the stellar masses from broadband fluxes via spectral energy distribution (SED) fitting: on the one hand, the uncertainty of emission lines contributing to the broadband flux in the spectra of highredshift galaxies could alter the derived observational SMF (Song et al. 2016) but also assumptions on the assumed star formation histories, dust contents and metallicities. On the other hand, the SEDs used to fit the broadband fluxes and  Table 1 for the different radiative feedback models studied in this work: Minimum derive stellar masses assume a fixed initial mass function (IMF). Recently, Chruslinska et al. (2020) have shown that an environment-dependent IMF leads to higher (lower) inferred star formation rates for low-mass (high-mass) galaxies with a cross-over point between such an IMF and a universal one lying at a star formation rate of ∼ 1 M yr −1 . In our simulations this star formation rate corresponds to stellar masses around M ∼ 10 9 M , which also marks the point in the SMFs below (above) which the observations fall below (exceed) our model results. Using an environment dependent IMF would therefore steepen the slope of the M − MUV relation (see e.g. Song et al. 2016) and hence the slope of the theoretical SMFs bringing them into better agreement with data.

Electron scattering optical depth and neutral fraction history
From the cifog ionization fields and the corresponding density fields, we derive the global volume-averaged ( χHI ) and mass-averaged ( χHI (m) ) reionization histories as Here χHI,i and ρi are the neutral hydrogen fraction and density in cell i, while N cell and ρ are the number of grid cells and the mean density of our simulation box. Further, the Thomson integrated electron scattering optical depth is cal- culated as with σT = 6.65 × 10 −25 cm −2 and H(z) being the Thomson cross section and the Hubble parameter at redshift z, respectively. The electron number density ne is determined by the mass-averaged ionization fraction χHI (m) (z) and the hydrogen and helium number densities, nH(z) and nHe(z). We assume that the fraction of singly ionized helium equals the fraction of ionized hydrogen, and that helium is fully ionized at z < 3. As noted before, for each radiative feedback model the fesc value (see Table. 1) has been adjusted to reproduce the optical depth τ (z dec ) for reionization (see Fig. 4). The resulting reionization histories are in agreement with existing constraints from quasars, Lyman-α emitters and GRBs (see Fig. 5).
In the beginning of reionization, radiative feedback has nearly no impact on the number of ionizing photons emitted and hence the evolution of χHI (except for the Jeans Mass model). However, as reionization proceeds, radiative feedback increasingly suppresses star formation in a rising number of low-mass galaxies as more and more ionized regions form. For a fixed optical depth (c.f. Fig. 4), we see from Fig. 5 that the reionization history becomes more extended, as the strength of the radiative feedback increases (from the Minimum, Weak Heating and Photoionization to the Strong Heating models). This is because a stronger radiative feedback model causes suppression of star formation in increasingly massive galaxies, leading to a larger reduction in the number of ionizing photons emitted. As a result, reionization slows down and the Universe becomes fully reionized later (c.f. Photoionization and Strong Heating models in Fig. 5).
Given this trend, it may seem surprising that the Jeans mass model, representing one of our strongest radiative feedback models, shows the same redshift evolution of χHI as the Minimum or Weak Feedback model. The reason for this behaviour originates from the fact that the radiative feedback strength does not increase relative to the SN feedback strength over time, i.e. the ratio of radiative and SN feedback characteristic masses remains constant (see Fig. 1). Hence, the ratio between the ionizing emissivity and halo mass is only lowered by a constant factor, and the lower production of ionized photons in low-mass halos can be compensated by an overall higher escape fraction fesc.
Moving from a constant fesc scenario to one where fesc decreases with halo mass, reionization is driven more by galaxies in low-mass halos, resulting in ionized regions being more centred around the smallest over-density peaks in the density field (c.f. higher mass-to-volume averaged neutral hydrogen fraction in top panel of Fig. 5). Since our simulations are tuned to reproduce the same optical depth value which depends on the mass averaged ionization fraction (and thus on the correlation strength between the IGM density and the time of reionization zreion), we find the model where fesc scales with the ejected gas fraction to accelerate towards the end of reionization, resulting in an earlier ionization of the IGM.

THE IMPACT OF RADIATIVE FEEDBACK ON EARLY GALAXY POPULATIONS
In Section 2.4 we have seen that the degree by which star formation in a galaxy is suppressed by radiative feedback depends sensitively on its individual reionization history and the redshift evolution of the characteristic mass for radiative feedback (Mc). Here we explore how the star formation histories of a representative sample of galaxies depend on their location in the cosmic web, as a function of their gravitational potentials and reionization histories.

Global star formation rate density for different galaxy masses
We start by discussing the global star formation rate density (SFRD) for galaxies with different halo masses, shown in Fig. 6. As expected from hierarchical galaxy formation, low-mass halos appear earlier and are always more abundant than high-mass halos throughout cosmic time. While at high redshifts (z 12 − 13) low-mass halos (with M h = 10 8−9 M ) provide the majority (about 90% at z = 15) of the star formation rate density in all our models, more massive halos (M h 10 9 M ) start to dominate the SFRD as time proceeds (z 11 − 12). The reason for this discrepancy between halo abundance and SFRD are radiative and SN feedback processes that contribute to the suppression of star formation in low-mass halos.
Indeed, from Fig. 6 we can see that for the lowest-mass galaxies (M h = 10 8−9 M , solid lines), the SFRD rises with time at z > 10 before turning over at z 9 − 10. The initial rise of the SFRD is due to the increasing number of galaxies that emerge with time. However, as time proceeds, SN feedback suppresses star formation in increasingly more massive galaxies: the location of the peak in the SFRD, z 9 − 10, marks the time when at least half of these low-mass galaxies are fully affected by SN feedback, i.e. their entire gas mass is ejected (see Fig. 1) in addition to these low-mass galaxies inheriting a small, or even no, gas content from their progenitors. While the drop in the SFRD at z 9 is entirely due to SN feedback for the Minimum model, models including radiative feedback show an additional drop in the SFRD. In case of the Heating and Photoionization models this drop increases with time, since radiative feedback causes an additional suppression of star formation in a rising number of galaxies as an increasing fraction of the IGM becomes ionized. In contrast, the Jeans Mass model shows an overall lower SFRD with an increasing logarithmic difference with decreasing redshift compared to the other models. This increasing logarithmic difference in the SFRD with decreasing redshift is caused by the increasing number of low-mass galaxies (below the Jeans mass) in ionized regions that are affected by radiative feedback as soon as their environment is reionized.
For more massive galaxies, M h = 10 9−10 M (long dashed lines) and M h = 10 10−11 M (dash-dotted lines), the turn over in the SFRD at lower redshifts, z 7 − 8 and (4 − 5), is also caused by the increasing impact of SN feedback in decreasing the gas mass brought in by their merging progenitors with cosmic time. However, in contrast to the lowest-mass halos, only a fraction of the gas in the galaxies is ejected, shifting the turn-over to increasingly lower redshifts with increasing halo mass. Furthermore, for M h = 10 9−10 M , we see that only the strongest radiative feedback models (Strong Heating and Jeans mass) cause a noticeable decrease in the SFRD. Only in those models, the radiative feedback characteristic mass exceeds 10 9 M significantly for galaxies whose environment was reionized at z 7−8 (when the majority of the volume becomes ionized). While theoretically, the star formation in M h = 10 10−11 M halos reionized very early on could also be suppressed by radiative feedback in the Strong Heating model, this is not the case in our reionization scenarios where the IGM was predominantly reionized at z 9 (see also the reionization history in Fig. 5). For high-mass galaxies, the Jeans mass model shows an overall lower SFRD, albeit with a constant logarithmic offset. While the lower SFRD is again due to the immediate effect of radiative feedback upon reionization, the constant logarithmic offset (and not an increasing one as for low-mass sources) is caused by the fact that the fraction of 10 9−10 M halos affected by radiative feedback remains constant over time. Again, here the lower gas content in 10 9−10 M halos is caused by their progenitors being gas-poorer for stronger radiative feedback.
We note that similar trends of the global SFRDs for various halo mass bins have also been found in the CoDaI simulation (c.f. Ocvirk et al. 2016;Dawoodbhoy et al. 2018, and see Section A2.2).

Impact of the reionization redshift on the SFR
In this section, we start by discussing the (average) redshiftdependent stellar mass (M ) and halo mass (M h ) assembly histories for galaxies at z = 5 as a function of their reionization redshift, as shown in Fig. 7. In the same figure, we also show the number of galaxies (N gal ) occupying this z − zreion plane, over which the stellar mass and halo mass assembly histories have been averaged.
We find N gal to reflect the hierarchical structure formation scenario, where a massive galaxy at a given time has formed and (in an inside-out reionization scenario) reionized its environment earlier than a less massive galaxy. For M h = 10 10−10.5 M halos at z = 5, N gal remains almost constant at z < zreion (see fourth panel in the first row). In contrast, as we go to the least massive galaxies in our simulation, N gal starts to rise towards smaller z values on the x-axis (e.g. for M h = 10 8−8.5 M and zreion = 12, there are about N gal 10 2 galaxies at z = 11, while N gal 10 5 at z = 5). Whether N gal remains constant at z < zreion also provides an indication of whether the galaxies in the respective halo mass bin have been able to ionize their grid cell in the simulation alone over the course of their life. From the figure we see that in our simulations the threshold for having emitted enough ionized photons to ionize the respec- Figure 7. As a function of redshift, we show the number of galaxies (first row), the averaged halo mass histories (second row), and averaged stellar mass histories (third row) binned by the reionization redshift of the galaxy for the Photoionization model. We show results for different final halo masses at z = 5, as marked above each column. For a given halo mass M h and reionization redshift z reion bin we average over the halo or stellar mass summed over all its progenitors at redshift z. The black solid line marks z = z reion . tive grid cell lies around M h 10 9 M , i.e. when star formation is not severely suppressed by feedback. From the same figure, we also see the build-up of stellar and the halo mass of galaxies. In our simulations a galaxy with M h 10 10.5 M (with M 10 8 M ) at z = 5 forms early on at z 17 (i.e. zreion 17) with a mass of M h 10 8 M (M 10 5.5 M ) and continuously accumulates mass reaching M h 10 9.5 M (M 10 7 M ) at z 11. Furthermore, we note that in our model, reionization proceeds in an insideout fashion, i.e. the ionization fronts move from over-to under-dense regions. Consequently, zreion is correlated to the underlying density field and thus the halo mass. Halo masses close to the upper limit of the halo mass bin are more likely to be found at high zreion values, while those close to the lower limit of the halo mass bin have lower zreion values. This effect can be seen in the second row in Fig. 7 where we show the average halo mass history of z = 5 galaxies in the given zreion and halo mass bins, e.g. for M h = 10 10−10.5 M the final halo mass at z = 5 is ∼ 10 10 M for zreion = 6 and ∼ 10 10.5 M for zreion = 15. In the following we will refer to this effect as the positive zreion − M h correlation effect. Also, in more massive halos with M h 10 9 M , the stellar mass follows the growth of the halo mass. However, in less massive halos (M h 10 9 M ), the stellar mass does not follow the growth of the halo mass but drops as star formation is increasingly suppressed by SN and radiative feedback (c.f. first panels in Fig. 7).
We now discuss the (average) redshift evolution of the star formation rate histories (SFH) for galaxies in a given halo mass bin at z = 5 as a function of the redshift zreion when the surrounding region of the galaxy became ionized in Fig. 8. We remark that for z > zreion the shown star formation rate (SFR) is averaged over increasingly fewer galaxies with increasing redshift (see last row in Fig. 7) and not affected by radiative feedback. Given our interest in studying the impact of radiative feedback on the SFH, we limit our discussion to z zreion (i.e. galaxies above the black line).
Since the Minimum model effectively corresponds to the case of SN feedback only, the SFR at redshift z is basically independent of zreion (when accounting for the zreion − M h correlation effect) and only depends on the mass of the galaxy and its redshift z (first row in Fig. 8). In agreement with the SN filtering mass (see dash-dotted grey line in Fig.  1), we can see that SN feedback suppresses star formation only in halos with M h 10 9.5−10 M with the suppression for a given M h increasing with decreasing redshift. For more massive halos at z = 5, the SN energy is not large enough to push out most of the gas, and the SFR continues to rise with decreasing redshift as galaxies becomes more massive. For galaxies in halos with M h = 10 9.5−10 , we see that the SFR peaks again around z 7 − 8 (see also Fig. 6): due to the increasing characteristic mass for SN feedback with decreasing redshift, these halos become increasingly SFR suppressed by SN feedback with time.
Including radiative feedback changes the SFHs since they become dependent on the reionization redshift zreion. For low-mass galaxies, the SFR at redshift z decreases as the galaxy is located in a region that has been reionized earlier (higher zreion value) due to the higher radiative feedback characteristic mass Mc (see discussion in Section 2.4). For example, in the case of the Photoionization model, M h = 10 8−8.5 M halos at z = 7 have a SFR ∼ 10 −3 M yr −1 for zreion = 7, while the SFR drops to a value as low as ∼ 10 −5 M yr −1 for zreion = 15. Furthermore, the strength and time of radiative feedback depends strongly on the gravitational potential of the galaxy. The less massive a galaxy is, the shallower is its potential, the earlier and the more its star formation is suppressed by radiative feedback: from low to high halo masses (left to right in Fig. 8) the relative drop of the SFR due to radiative feedback (not SN feedback) since its maximum decreases: while the SFR drops by roughly 2-4 (1-4) orders of magnitude for M h = 10 8−9.5 M (M h = 10 8−8.5 M ) halos, there is nearly no drop for M h = 10 9.5−10 M (M h = 10 9−9.5 M ) in the Strong Heating (Photoionization) model. As the radiative feedback model becomes more efficient, the drop in the SFR becomes not only more pronounced but extends also to galaxies with higher halo masses. While the maximum halo mass being affected by radiative feedback at z 5 is M h 10 9 M for the Photoionization model, it increases to M h 10 10 M for the Strong Heating model, respectively.
We also note that the SFHs for low-mass (10 8−9 M ) and mid-mass (10 9−10 M ) halos differ from the point of reionization: while the SFR drops continuously for low-mass halos (see M h 10 9 M for (Photoionization and) Heating models), it rises first and drops then for mid-mass halos (see 10 9 M M h 10 10 M for Heating models). These findings are consistent with those from the radiation hydrodynamic simulation CoDaI (Dawoodbhoy et al. 2018). Besides radiative feedback, galaxies in low-mass galaxies are also subject to SN feedback, which causes the SFR to decline from their start. In contrast, galaxies in mid-mass halos are only marginally affected by SN feedback but subject to the timedelayed radiative feedback in the Photoionization and Heating models: the SFR keeps rising as long as the halo mass of the galaxy exceeds the radiative feedback characteristic mass Mc, and starts dropping continuously as soon as Mc has surpassed the galaxy's halo mass. Since Mc is a function of the reionization redshift zreion of a chosen galaxy, the resulting peak in the SFR also depends on zreion. As the strength of radiative feedback of a model increases, its redshift shifts closer to zreion, since Mc surpasses the considered halo mass shorter after the reionization of the galaxy. During reionization, we find that Mc never approaches the halos with M h 10 10 M , and consequently SFHs are not affected by radiative feedback.
Finally, we comment on the results from the Jeans Mass model. This model differs from the Photoionization and Heating models, since the degree of star formation suppression does not depend on when a galaxy's environment became reionized (zreion) but on its redshift and the galaxy's halo mass (similar to SN feedback). Consequently, although this model results in a strong suppression of star formation in low mass halos, the SFHs effectively mimic a higher SN characteristic mass.

Impact of the reionization redshift on the gas fraction
In this Section we discuss how the relation between the halo mass and star formation rate and halo mass and the initial gas-fraction M i g /M h evolve with redshift and the reionization feedback model used. As we have detailed in Sec. 2.2, the halo mass is determined by mergers and accretion, while the initial gas mass (i.e. gas mass present after mergers and accretion but before star formation and SN feedback) and hence the star formation rate depends on the star formation and assembly history of a galaxy.
While in real galaxies, processes such as gas accretion, gas ejection and star formation overlap at times, the complex interplay of the physical processes in galaxies impedes a direct all-encompassing solution but forces us to execute them sequentially. Hence, in our model, gas accretion (including mergers) and evaporation happen at the beginning of a time step, while gas ejection and star formation is executed at the end. In this Section, we comment on the initial gas fraction M i g /M h . In Fig. 9 we show the median SFR (top panels) and gas fractions (bottom panels) as a function of the halo mass M h for the Minimum and Strong Heating models (left and right panels, respectively). We see that the gas fraction and SFR show the same increasing trend with halo mass M h (c.f. bottom and top panel in Fig. 9) as expected from our model (see Section 2.2.2). However, while the SFR continuously rises with increasing M h , the gas fraction increases with M h before saturating to a value that is around 60% of the cosmological baryon fraction Ω b /Ωm at M h 10 10.5−11.5 for z = 10 − 6. While the SFR is a direct tracer of the gas mass in the galaxy, the increasing gas fraction with rising halo mass reflects the deepening of the gravitational potential and the associated decrease in the efficiency with which SN explosions can eject gas from galaxies. The saturation of the gas fraction to a value lower than the cosmological baryon fraction towards and for massive galaxies stems from the gas loss of their lower mass progenitors. When radiative feedback is added, we find the gas fraction to decrease, with the decrease being strongest for the lowest halo masses (c.f. Strong Heating to Minimum model in Fig. 9). Galaxies in shallower potentials have shorter dynamical timescales and hence they are the quickest to adjust their gas density to the heating by reionization and to increase their Jeans mass or radiative feedback characteristic mass. Similarly, this buildup of the galaxy's actual Jeans mass (or increasing radiative feedback characteristic mass) leads to a stronger drop in the gas fractions, relative to the Minimum model, as redshift decreases.
We note that the median gas fraction in low mass halos increases when going from halos with M h 10 8.5 M to those with M h 10 8 M . This is an artefact of our assumption that each newly formed halo starts with an initial gas fraction that equals the cosmological baryon-matter-ratio. However, by the end of the time step SN feedback will push out most/all gas, resulting in the gas fractions in subsequent time steps to include all key physical processes of gas accretion, evaporation, ejection and loss due to star formation. If the "final" gas fraction Mg/M h were plotted, they would drop further for low-to mid-mass halos (i.e. where SN feedback ejects a noticeable fraction of gas from the galaxies) and even to zero for M h 10 9 M .
Finally, we find our gas fractions to be lower than those obtained from radiation hydrodynamical simulations (e.g. Hasegawa & Semelin 2013;Okamoto et al. 2008;Gnedin & Kaurov 2014;Pawlik et al. 2015;Wu et al. 2019), since our semi-analytical galaxy evolution model accounts only for the kinetic but not temperature effects of SN feedback, i.e. all the gas that does not form stars can be ejected out of the potential rather than just being heated as in radiation hydrodynamical simulations. Ignoring gas ejection from SN feedback and accounting only for gas evaporation and reduced accretion from radiative feedback (fgΩ b /Ωm versus M h ), our model yields results comparable to those of radiation hydrodynamical simulations. A more detailed discussion can be found in the Appendix A2.3.

THE IMPACT OF RADIATIVE FEEDBACK ON THE 21CM POWER SPECTRUM
Current and forthcoming radio interferometers will provide measurements of the 21cm signal that arises from the spin flip transition in neutral hydrogen H I . One of the key quantities to be measured is the 21cm power spectrum that provides a tracer of the size distribution of the ionized and neutral regions. Here we discuss whether radiative feedback has an impact on the sizes of the ionized regions and the associated redshift evolution of the 21cm power spectra. For this purpose, we show the 21cm power spectra at redshifts z = 9.0, 8.0, 7.0, 6.7 and 6.5 (Fig. 10) corresponding roughly to χHII ∼ 0.1, 0.25, 0.5, 0.75 and 0.9, respectively, as well as the ionization maps of our simulation slices at z = 9.2 and z = 7.6 ( Fig. 11). First, we consider all reionization scenarios where a constant ionizing escape fraction has been assumed (i.e. all except Early Heating). For these scenarios, the ionization fields and 21cm power spectra are very similar, however we note a few small differences. Firstly, while the 21cm power spectra are effectively identical during the first two thirds of reionization, the increasing suppression of star formation (and hence the emission of ionizing photons) in low-mass galaxies extends reionization as the radiative feedback strength is increased from the Minimum to the Strong Heating model as discussed in Sec. 3.3. With the amplitude of the 21cm power spectrum tracing the mean neutral hydrogen fraction, we find the 21cm power spectrum amplitude to reflect the reionization histories, with the Strong Heating model having a marginally higher amplitude than the weaker radiative feedback models and the Jeans Mass model. Secondly, in terms of the reionization topology, a stronger radiative feedback translates into less ionizing photons being emitted from low-mass galaxies, resulting in smaller ionized regions around these sources. Focusing on the smallest ionized regions in Fig. 11, one can see that ionized regions of cell size are more ionized and abundant in the Minimum than in the Jeans Mass model. In case of the Strong Heating model, a difference to the Minimum model is barely seen, as radiative feedback has not become effective at these redshifts.
Across our reionization scenarios, we note that the ionization fields and 21cm power spectra of the Early Heating model differ strongly from all other models during the first half of reionization: due to the higher (lower) abundance of small-sized (large-sized) ionized regions respectively, the large scale 21cm power is reduced. Given that the ionizing escape fraction in this scenario scales with the ejected gas fraction, and hence decreases with halo mass, this change is expected and is in agreement with previous works (e.g. Seiler et al. 2019).
If radiative feedback was even stronger than our Strong Heating model, i.e. star formation in low-mass sources was more suppressed and star formation in even higher-mass galaxies became affected, it would have the two effects. Firstly, this would result in more-small sized ionized regions, since low-mass galaxies forming in neutral regions are not affected by radiative feedback (and delayed SN feedback) in the first ∼ 20 Myrs of their life. Secondly, we would find larger-sized ionized regions to be smaller and less connected, since a stronger radiative feedback reduces star formationand hence the number of ionizing photons produced -not only in more massive galaxies, but also in low-mass galaxies located in the vicinity to the most massive galaxies that contribute to the connectivity of the ionized regions.
In summary, even the strongest radiative feedback hardly affects the galaxies that determine the ionization topology on > 1 − 2 Mpc scales, which is in agreement with previous findings (e.g. Dixon et al. 2016). Different functional forms of the ionizing escape fraction dominate over the intrinsic ionizing emissivity distribution of the underlying galaxy population and can lead to significantly different ionization topologies. We will explore this finding in more detail in future works.

IMPACT OF STELLAR POPULATION SYNTHESIS MODELS
In all the previous sections we have presented the results using the single stellar population synthesis (SPS) model starburst99 (Leitherer et al. 1999). Over the past decade, it has been increasingly recognised that the effects of stellar multiplicity can not be neglected when modelling the evolution and observable properties of young stellar populations, since about 70% of massive stars are part of binaries (Sana et al. 2012) and exchanging mass with their companions affects their structure and evolution. For this reason, we also have run our simulations where the ionizing emissivity of a starburst is derived with the bpass SPS model 12 (Eldridge et al. 2017). While the initial ionizing emissivity at 4Myrs is similar in both, starburst99 and bpass, models, at later times the ionizing emissivity decreases slower in time in the bpass model (see equations 17 and 18). Thus, we find the intrinsic ionizing luminosities of all galaxies to be about a factor 10 higher. However, to fit the Planck optical depth, this increased ionizing emissivity from galaxies is compensated by a lower escape fraction of these ionizing photons, f BPASS esc 0.1 × f S99 esc . This re-normalisation causes the effective ionizing emissivities to be about equal, and reionization histories as well as the sizes, shapes and distribution of the ionized regions agree well with each other. This strong agreement for a constant relation between f BPASS esc and f S99 esc arises because the number of ionizing photons at each time is dominated by the youngest stellar populations, given that the SFR remains constant within the current time step (10−30 Myr) and drops strongly in the subsequent time step. Since we derive the distribution of the ionized regions from the cumulative number of ionizing photons, we do not expect that smaller time steps would result in significantly differing ionization topologies.
It is interesting to note that Rosdahl et al. (2018) found higher ionizing escape fractions when using binary stellar populations in their radiation hydrodynamical simulation, while our simulations require lower escape fractions in order to agree with observations. A reason for this disagreement might be that we model the global gas properties of galaxies, while the simulations of Rosdahl et al. (2018) resolve the overall gas density distribution and follow the radiation emitted by the stellar populations in the galaxy.
While the inclusion of binaries may not significantly affect the ionization topology, their inclusion would lead to an increase in certain emission line strengths of galaxies for a constant ionizing escape fraction scenario. We aim to explore this effect in a future work.

CONCLUSIONS
We have developed the semi-analytical/semi-numerical model astraeus that self-consistently derives the evolution of galaxies and the reionization of the IGM based on the merger trees and density fields of a DM-only N-body simulation. astraeus models gas accretion, star formation, SN feedback, the time and spatial evolution of the ionized regions, accounting for recombinations, H I fractions and photoionization rates within ionized regions, and radiative feedback. astraeus aims at studying the galaxy-reionization interplay in the first billion years.
In this paper, we have focused on the impact of radiative feedback on the interlinked processes of galaxy evolution and reionization. We have considered 6 radiative feedback and reionization models that cover the physically plausible parameter space. Our models (cf. Table 1) include scenarios (1) where the IGM is heated between 2 × 10 4 K and 4 × 10 4 K, (2) where the gas fraction a galaxy can maintain upon ionization is described by the filtering mass according to Gnedin (2000); Naoz et al. (2013) or the Jeans mass, and (3) where the ionizing escape fraction ranges from being constant to scaling with the fraction of gas ejected from the galaxy (i.e. increasing with decreasing halo mass). Comparing the results from these models, we find: (i) During the Epoch of Reionization, radiative feedback affects only galaxies with halo masses less than M h ∼ 10 9.5−10 M corresponding to stellar masses less than M ∼ 10 7.5−8 M ( Fig. 3 and Section 3.2) and UV luminosities lower than MUV −16 ( Fig. 2 and Section 3.1). Hence, increasing the strength of radiative feedback, both, the faint end slope of the UV LFs as well as the low-mass end of the SMFs, flatten increasingly as a larger fraction of the IGM becomes ionized (ii) Radiative feedback causes lower gas fractions and suppressed star formation in low-mass galaxies. The strength of radiative feedback increases both with decreasing halo mass and the longer the IGM gas around the galaxy has been ionized: while the temperature of the gas in the galaxy is instantaneously increased when the region becomes ionized, its density adapts to the raised temperature only on the dynamical time scale of the galaxy. The gradual decrease of gas density results in an increasing Jeans mass, steadily lowering the amount of gas the galaxy can maintain. Moderate radiative feedback (e.g. Sobacchi & Mesinger 2013b) leads to a full suppression of star formation in galaxies with M h = 10 8−9 M , while stronger radiative feedback can reduce the gas content even in galaxies with M h = 10 9−10 M ( Fig. 8 and Section 4).
(iii) The radiative feedback strength and hence the degree of star formation suppression depends on the time when the environment of a galaxy becomes ionized; the earlier the IGM around a galaxy is reionized, the more star formation is affected by radiative feedback (Fig. 8 and Section 4). Hence, the patchiness of reionization is imprinted in the star formation histories of low-mass galaxies. However, due to the delayed onset of radiative feedback (Fig. 1) and the rapidness of reionization (Fig. 5), for the Photoionization and Heating radiative feedback models the differences in observables such as the UV LFs are minor. This is discussed in more detail in an accompanying paper (Ucci et al. 2020).
(iv) Radiative feedback does not affect the ionization topology on scales larger than 1 − 2 comoving Mpc, as it affects only the star formation in low-mass galaxies during reionization. Consequently, we find the power spectrum of the 21cm signal to be hardly altered on such scales for even the strongest radiative feedback models ( Fig. 10 and Section 5). Galactic properties such as the ionizing escape fraction fesc, however, are far more crucial in determining the ionization topology (Fig. 11).
(v) Changing the stellar population synthesis model, i.e. going from stellar populations modelling single stars to those including binaries, has no visible effect on the ionization topology after the escape fraction is tuned to yield reionization histories in agreement with observations.
(vi) Our results on the impact of radiative feedback on galactic properties obtained with astraeus are in agreement with the findings in radiation hydrodynamical simulations (e.g. Ocvirk et al. 2016;Wu et al. 2019), such as the galaxies affected by radiative feedback, the degree of star formation suppression in low-mass halos as well as its dependence on the time when a galaxy's environment was reionized (c.f. Section 4 and Appendix A). In particular, the results for the Photoionization or a Medium or Early Heating (T0 = 4 × 10 4 K, Mc = MF ) model agree best with the results of radiation hydrodynamical simulations. However, even for those models, we find astraeus to have lower baryon fractions in low-mass galaxies compared to radiation hydrodynamical simulations, since astraeus accounts only for the kinetic but not temperature effects of SN feedback and not for both as numerical simulations do.
In the following we list some caveats. Firstly, our framework assumes the same low metallicity for all stellar populations and does not account for the evolution of the metallicity in gas and stars, the evolution of the dust content in galaxies, and the associated changes in the stellar populations. This assumption probably begins to break down towards lower redshifts z 6, and we will present a selfconsistent treatment of the metallicity evolution within the astraeus framework in a follow-up paper (Ucci et al., in prep.). Secondly, the redshift of reionization of each grid cell that determines the onset of radiative feedback depends on the grid resolution. In our simulations, a galaxy needs to ionize at least a volume of ∼ 0.05 Mpc −3 before its cell is considered reionized. A finer resolution would lead to more accurate results for fainter sources that would then be able to ionize their smaller grid cell earlier.
Finally, we note that astraeus reproduces the key observations at z 5 using only three mass-and redshiftfree parameters and requires much less computation time than comparable radiation hydrodynamical simulation. The underlying code is publicly available 13 , and written in a modular fashion that allows to incorporate new physics easily. This and especially its short computing time make astraeus well equipped to pursue quick and efficient parameter studies, exploring the impact of varying galaxy properties on the ionization topology and their visibility with current and forthcoming telescopes. 13 https://github.com/annehutter/astraeus ACKNOWLEDGMENTS All the authors acknowledge support from the European Research Council's starting grant ERC StG-717001 ("DEL-PHI"). PD acknowledges support from the NWO grant 016.VIDI.189.162 ("ODIN") and the European Commission's and University of Groningen's CO-FUND Rosalind Franklin program. GY acknowledges financial support from MICIU/FEDER under project grant PGC2018-094975-C21. We thank Peter Behroozi for creating and providing the rockstar merger trees of the VSMPL simulation, and AH thanks Manodeep Sinha for discussions about the benefits of local horizontal merger trees. The authors wish to thank V. Springel for allowing us to use the L-Gadget2 code to run the different Multidark simulation boxes, including the VSMDPL used in this work. The VSMDPL simulation has been performed at LRZ Munich within the project pr87yi. The CosmoSim database (www.cosmosim.org) provides access to the simulation and the Rockstar data. The database is a service by the Leibniz Institute for Astrophysics Potsdam (AIP).
sources and needs to boost their star formation efficiency to produce the sufficient stellar mass. Wu et al. (2019) showdespite a DM particle mass of mDM = 1.4 × 10 6 M -the lowest SFRD at higher redshifts; in contrast to the other mentioned works, their density threshold for star formation is given by an absolute physical density and not a relative one to the mean, increasingly reducing star formation towards higher redshifts.

A2.2 Star formation rate per halo
Comparing the SFR as a function of halo mass, we find our results to be in agreement with the results from radiation hydrodynamical simulations presented in Pawlik et al. (2015), Wu et al. (2019) and Ocvirk et al. (2018) ;Hasegawa & Semelin (2013) show SFRs that are about 0.5-1 dex higher, possibly due to their SN feedback implementation and an IMF that results in a lower SN rate per unit stellar mass formed. In particular, after most of the Universe has been ionized, all simulations find a drop in the SFR below 10 9−9.5 M , which is in agreement with our Photoionization model. Ocvirk et al. (2018) find a slightly stronger suppression in low mass halos, corresponding to a radiative feedback model that lies between our Photoionization and Strong Heating models.
In Wu et al. (2019), we see the same effect that we find in our SFR histories: switching the UVB on at z ∼ 10.7 causes an early (z > 6) suppression of star formation in low mass halos. However, when deriving the radiation field selfconsistently, the suppression becomes only visible towards lower redshifts. As the Universe becomes more ionized, more galaxies are subject to radiative feedback, and the more time has been passed since a galaxy's environment has been ionized and adapted its gas density to the new temperature. Our SFR histories also echo the results in Dawoodbhoy et al. (2018), which are based on the CoDaI simulation (Ocvirk et al. 2016). In CodaI, the Universe reionizes just at z 4.6 and consequently lower SFRs for halos with M h < 10 9 M are found compared to CoDaII (Ocvirk et al. 2018) and our simulations. However, they find the same trends of the SFR with halo mass and reionization redshift: SFR in low mass halos (M h = 10 8−9 M ) is suppressed immediately upon reionization, SFR in mid mass halos (M h = 10 9−10 M ) increases upon reionization but drops later in time, and SFR in massive halos (M h > 10 10 M ) is not affected by radiative or SN feedback and continues to increase with time. The mass ranges for these three regimes shift to higher masses as the strength of the radiative feedback is increased, allowing us to pinpoint their radiative feedback strength to lie between our Photoionization and Strong Heating models.
distribution of the SN energy to neighbouring particles results in distributing the energy to too much mass, leading to a smaller increase in temperature, a shorter cooling time, and a weaker stellar feedback causing increased star formation (overcooling problem). Hasegawa & Semelin (2013) use this description, which possibly contributes to their higher SFRD.

A2.3 Baryon or gas fractions:
It is important to note that our stellar feedback model represents an upper limit of SN feedback on gas in galaxies: firstly, our model does not account for the cooling of gas within galaxies, and secondly, the gas energised by SN is fully ejected from the halo. This is different to the SN feedback prescriptions in hydrodynamical simulations where gas heated by SN can cool (thermal SN feedback) or gas accelerated by SN energy does not reach the edge of the halo (kinetic SN feedback). Hence, in our model, more gas leaves the galaxy than in hydrodynamical simulations, and consequently we find lower gas or baryon fractions in galaxies that are located in lower mass halos. It is interesting to notice that if we consider the gas fraction fg in Equation 3, which is purely due to radiative feedback and given by the respective filtering mass, we find the Photoionization or a Early Heating (T0 = 2 × 10 4 K, Mc = MF ) model to agree best with the findings at z 6 − 7 in Hasegawa & Semelin (2013); Okamoto et al. (2008); Gnedin & Kaurov (2014) 16 ; Pawlik et al. (2015) and the UVB model in Wu et al. (2019). This finding indicates that, indeed, SN feedback in the aforementioned hydrodynamical simulations is not strong enough to eject the gas from the halo that has been heated or accelerated by SN explosions (but is able to suppress star formation by either heating the gas or decreasing the gas density in star formation sites through winds).

APPENDIX B: STAR FORMATION RATE DENSITY & STELLAR MASS DENSITY
For all our models, we show the star formation rate densities (SFRDs) and stellar mass functions for different UV luminosity selection criteria, i.e. MUV < −13, −15, −17, −18 in Fig. B1 and B2, respectively. Applying the same selection criterion as the observations (i.e. MUV < −17), we find all our models to be in good agreement with the observational data points for both the SFRDs and SMDs, tracking the overall increase in star formation and stellar mass with cosmic time as new galaxies form and existing ones continuously grow. We can see that the SFRDs and SMDs of our models agree well with each other, since firstly the suppression of star formation by radiative feedback is mostly found in lowmass and UV faint galaxies where star formation is already limited by SN feedback, and secondly the reduction of star formation by radiative feedback is time delayed. Only the Jeans Mass model deviates when faint galaxies are included: radiative feedback suppresses star formation instantaneously when the environment of the galaxy becomes ionized. In addition, it affects only galaxies in low-mass halos, and hence reduces the SFRD and SMDs only for MUV > −13 and −15.