# Computationally Efficient Reionization in a Large Hydrodynamic Galaxy
Formation Simulation

James E. Davies1, Simeon Bird2, Simon Mutch3,4, Yueying Ni5,6, Yu Feng7,
Rupert Croft6, Tiziana Di Matteo6, and J. Stuart B. Wyithe3,4,
1Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
2Department of Physics & Astronomy, University of California, Riverside, 900
University Ave., Riverside, CA 92521, USA
3School of Physics, The University of Melbourne, Parkville, Victoria 3010,
Australia
4ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)
5Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA
6McWilliams Center for Cosmology, Carnegie Mellon University, Pittsburgh PA,
15213
7Berkeley Center for Cosmological Physics, University of California at
Berkeley, Berkeley, CA, 94720, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Accuracy in the topology and statistics of a simulated Epoch of Reionization
(EoR) are vital to draw connections between observations and physical
processes. While full radiative transfer models produce the most accurate
reionization models, they are highly computationally expensive, and are
infeasible for the largest cosmological simulations. Instead, large
simulations often include EoR models that are pre-computed via the initial
density field, or post-processed where feedback effects are ignored. We
introduce Astrid-ES, a resimulation of the Astrid epoch of reionisation
$20>z>5.5$ which includes an on-the-fly excursion-set reionization algorithm.
Astrid-ES produces more accurate reionization histories without significantly
impacting the computational time. This model directly utilises the star
particles produced in the simulation to calculate the EoR history and includes
a UV background which heats the gas particles after their reionization. We
contrast the reionization topology and statistics in Astrid-ES with the
previously employed parametric reionisation model, finding that in Astrid-ES,
ionised regions are more correlated with galaxies, and the 21cm power-spectrum
shows an increase in large scale power. We calculate the relation between the
size of HII regions and the UV luminosity of the brightest galaxy within them.
Prior to the overlap phase, we find a power-law fit of
$\mathrm{log}(R)=-0.314M_{\mathrm{UV}}-2.550\mathrm{log}(1+z)+7.408$ with a
standard deviation $\sigma_{R}<0.15\mathrm{dex}$ across all mass bins. We also
examine the properties of halos throughout reionization, finding that while
the properties of halos in the simulation are correlated with the redshift of
reionisation, they are not greatly affected by reionisation itself.

###### keywords:

dark ages, reionization, first stars – intergalactic medium – early universe

††pubyear: 2023††pagerange: Computationally Efficient Reionization in a Large
Hydrodynamic Galaxy Formation Simulation–Computationally Efficient
Reionization in a Large Hydrodynamic Galaxy Formation Simulation

##  1 Introduction

Constraining the epoch of reionization requires theoretical modelling which
spans a vast range of physical scales, from the growth of large scale
structure and formation of early galaxies (see e.g. Dayal & Ferrara (2018)),
to the interactions between hydrogen atoms and the Lyman Continuum, Lyman-
Werner, and X-ray photon backgrounds from these early sources, as well as the
CMB (see e.g Furlanetto et al. (2006)). As such there is always a trade-off
between the level of detail in our physical models and their computational
cost. More expensive models involve detailed prescriptions for the generation
and propagation of ionising photons (Kannan et al., 2021; Ocvirk et al., 2020;
Lewis et al., 2022), allowing comparison with a wider range of observational
data. On the other hand, more computationally cheap modelling frameworks use
more approximate methods to simulate the connection between early ionising
sources and the topology of reionisation (Mesinger et al., 2011; Mutch et al.,
2016; Molaro et al., 2019; Choudhury & Paranjape, 2018; Hutter et al., 2021).
As a result these models can be run on larger volumes, at higher resolutions,
and over a wider range of model parameters with the same computational
resources.

The most detailed models of reionization are radiative transfer models, which
explicitly model photons within a simulation and follows their propogation
through the IGM, accounting for absorption along the way (Mihalas & Mihalas,
1984). They take into account the complex geometry and distributions of both
ionising sources and absorbers. However the high dimensionality of the
equations which discretise over position, time, angle and frequency, as well
as the modelling of several chemical processes in the intergalactic medium
make this approach highly computationally expensive.

Figure 1: Galaxy stellar mass functions for Astrid-ES (solid lines), compared
with observational data from Stefanon et al. (2017); Grazian et al. (2015);
González et al. (2011); Duncan et al. (2014), and Song et al. (2016). We also
show the stellar mass function from Astrid (Bird et al., 2021) as a grey
dashed line, which is consistent with Astrid-ES, as the reionisation timing
does not have a significant effect on the total stellar mass in most galaxies
(see section 4.3). As in Bird et al. (2021), we find that galaxy stellar
masses are $\sim 1\sigma$ below the bulk of observations at redshifts $7\leq
z\leq 8$.

The largest hydrodynamic simulations such as Bluetides (Feng et al., 2016) and
Astrid (Bird et al., 2021) are too large to utilise radiative transfer models
whilst remaining computationally feasible, although recent works have coupled
radiative transfer with hydrodynamic simulations up to a cube of side length
$<70h^{-1}\mathrm{cMpc}$ (Kannan et al., 2021; Ocvirk et al., 2020). Instead,
these large simulations use parametric models which pre-determine reionization
history according to the initial conditions of the simulation (Battaglia et
al., 2013; Oñorbe et al., 2019). While these methods were shown to reproduce
certain large-scale characteristics of a radiative transfer simulation, the
topology of reionization in these models is disconnected from the ionising
sources in the simulation volume. There are also efficient post-processed
radiative transfer models which calculate an inhomogeneous EoR history from
the sources in a simulation, either ignoring inhomogeneous photoheating
feedback onto the sources (Keating et al., 2018) or requiring a second
hydrodynamic simulation to fully capture the effects of reionization on
galaxies (Puchwein et al., 2023).

The excursion-set reionization algorithm has been used extensively to produce
realistic ionisation topologies in large volumes whilst being very
computationally efficient, and occupies a middle ground between radiative
transfer models and analytic models. The algorithm can be applied to an
evolved density field in semi-numerical models (Mesinger & Furlanetto, 2007;
Mesinger et al., 2011; Choudhury & Paranjape, 2018; Molaro et al., 2019; Park
et al., 2019; Oñorbe et al., 2019) or discrete ionising sources in semi-
analytic galaxy evolution models (Mutch et al., 2016; Hutter, 2018; Hutter et
al., 2021). The radiation field is approximated by constructing grids of gas
density and ionising emissivity, filtering them on a series of scales, and
flagging cells as ionised in regions where the number of ionising photons
generated by the stars within the filter scale is greater than the number of
neutral hydrogen atoms. The large scale structure of reionization produced by
similar models has been shown to compare favorably to radiative transfer
models (Hutter, 2018; Hassan et al., 2022).

In this paper, we have adapted the 21cmFAST (Mesinger et al., 2011) excursion-
set algorithm to perform on-the-fly reionisation calculations within MP-
Gadget. We take the distribution of gas and stellar particles as input, and
direct its output of the UV background to the existing cooling model, which
evolves the ionisation and thermal state of each particle. This algorithm
allows us to study more realistic histories of reionization which are directly
dependent on the sources of ionising photons within the simulation and feeds
back onto those sources as the gas is heated. The excursion-set algorithm is
computationally cheap compared to the hydrodynamics, such that the overall
run-time of any simulation will be negligibly affected.

In section 2, we briefly introduce the main features of the MP-Gadget
codebase. In section 3, we detail the excursion-set algorithm and ionising
background model. In section 4, we compare the reionization history of the
excursion-set with the previous parametric model, and present new results
regarding the correlation between galaxy properties and the EoR.

##  2 MP-Gadget

MP-Gadget uses a pressure-entropy formulation of Smoothed Particle
Hydrodynamics (Read et al., 2010; Hopkins, 2013). This includes a multi-phase
star formation model (Springel & Hernquist, 2003; Vogelsberger et al., 2013),
cooling via radiative processes (Katz et al., 1996), metal cooling
(Vogelsberger et al., 2014), the effects of molecular hydrogen on star
formation (Krumholz & Gnedin, 2011), a type II supernovae feedback model
(Okamoto et al., 2010), helium II reionization (Upton Sanderbeck & Bird,
2020), and a super-massive black hole treatment including dynamical friction
(Di Matteo et al., 2005; Chen et al., 2021). The $250\mathrm{h^{-1}cMpc}$
Astrid simulation has been shown to produce UV luminosity functions, stellar
mass functions and specific star formation rates consistent with observation
in the redshift range $3<z<10$ (Bird et al., 2021), as well as black hole mass
functions in the redshift range $3<z<6$ (Ni et al., 2022). Both Astrid and our
resimulation of reionisation Astrid-ES introduced in this paper have a dark
matter particle mass of $6.75\times 10^{6}h^{-1}M_{\odot}$ and a gas particle
mass of $1.267\times 10^{6}h^{-1}M_{\odot}$. The Planck 2015 cosmology (Planck
Collaboration et al., 2016) is used for each simulation
$\\{\Omega_{\mathrm{m}},\Omega_{\mathrm{b}},\Omega_{\Lambda},h,\sigma_{8},n_{\mathrm{s}}\\}=\\{0.3089,0.0486,0.6911,0.6774,0.816,0.9667\\}$.

###  2.1 Parametric Reionization

The reionization model in MP-Gadget used in the Astrid simulation is a
parametric model outlined in Battaglia et al. (2013). We expect the EoR to
begin in the dense filaments of the cosmic web, as these areas are likely to
contain both larger numbers of galaxies and brighter galaxies overall.
Battaglia et al. (2013) drew spatial correlations between the large-scale
density field and the expected reionization history by analysing the outputs
of radiative transfer simulations. The result of this work is a Fourier-space
filter that can be applied to the density fields of simulations in order to
predict their reionization history without directly performing the expensive
radiative transfer calculation.

Figure 2: Top panel: Global reionization histories showing the volume (solid
line) and mass (dotted line) weighted neutral fraction of Astrid-ES, as well
as the volume-weighted fraction Astrid (Bird et al., 2021) (orange line). We
compare with observations from Dark pixel QSO damping wing, lyman alpha
emitter fraction, and equivalent width observations. (McGreer et al., 2015;
Inoue et al., 2018; Ouchi et al., 2018; Morales et al., 2021; Wold et al.,
2022; Hoag et al., 2019; Mason et al., 2019; Whitler et al., 2020; Greig &
Mesinger, 2017a; Davies et al., 2018; Greig et al., 2019; Wang et al., 2020).
Second panel: Thomson optical depth of each simulation, compared with the
Planck 2018 observation (Planck Collaboration et al., 2020). Third panel:
volume-weighted ionisation rate for all particles (solid line) and ionised
particles (dotted line) in Astrid-ES, as well as the ionisation rate applied
to ionised regions in Astrid Faucher-Giguère (2020) (dashed orange line).
Ionisation rates are compared with observations (Wyithe & Bolton, 2011;
Calverley et al., 2011; D’Aloisio et al., 2018; Becker et al., 2021). Bottom
panel: Temperature at mean density in the excursion set model, compared with
observations from Gaikwad et al. (2020). We do not show the temperature
history of Astrid, since heating from the initial ionisation of a particle is
not included.

.

In these models, the bias between overdensity and the redshift of reionization
$b_{mz}$ is fit to a functional form:

$b_{mz}(k)=\frac{b_{0}}{(k+k/k_{0})^{\alpha}}$ (1)

where $k_{0}$ and $\alpha$ characterises the scaling of the bias with Fourier
mode $k$. These values were found to be $k_{0}=0.185h^{-1}\mathrm{Mpc}$ and
$\alpha=0.564$. $b_{0}$ is the normalisation of this bias on large scales and
is derived from the extended Press-Schecter formalism (Bond et al., 1991)
$b_{0}=\frac{1}{\delta_{\mathrm{crit}}}=0.593$. The resulting bias is applied
alongside a spherical tophat filter of radius 1 Mpc to the initial conditions
of the simulation, which has been evolved to the desired median redshift of
reionization via second order Lagrangian perturbation theory. The resulting
$z_{\mathrm{reion}}$ field is applied in the simulation by exposing particles
within each cell with $z<z_{\mathrm{reion}}$ to a redshift-dependent UV
background calculated by Faucher-Giguère (2020). Since the bias and
reionization fields are pre-computed, this approach effectively adds no
computational cost to the simulation. However, the model does not directly
correlate the properties of the luminous sources in a simulation with its
reionization history, instead relying on correlation with the underlying
density field.

Figure 3: $250\times 250\times 5h^{-1}\mathrm{Mpc}$ neutral hydrogen fraction
slices in the Astrid-ES excursion set (top) and Astrid parametric (bottom)
reionization models. To account for the differing reionization histories,
columns correspond to different stages of reionisation, with a fraction of
ionised cells $Q_{\mathrm{HII}}\approx 0.22,0.50,0.97$ in each model. The
ionisation structure in the Astrid-ES is composed of overlapping regions which
are highly spherical, centred around bright galaxies. In Astrid, the
ionisation structure closely follows the cosmic filaments of the underlying
density structure.

##  3 Excursion-Set Algorithm

The large volume ($250h^{-1}\mathrm{cMpc}$ side-length) and high mass-
resolution ($6.75\times 10^{6}h^{-1}M_{\odot}$ per dark matter particle) of
Astrid means that a full radiative transfer model to estimate the UV
background is far too computationally expensive. The excursion-set algorithm,
however, is able to produce accurate reionization histories using the ionising
sources within the simulation while requiring very little memory and
computation time compared with the rest of the simulation. As a result it can
be performed during the simulation run with very little additional cost.

Our algorithm for calculating the reionization history within the simulation
is based on the version of 21cmFAST (Mesinger et al., 2011) with UV background
model (Sobacchi & Mesinger, 2013b) implemented in the Meraxes semi-analytic
galaxy evolution model (Mutch et al., 2016). First, we construct total mass,
stellar mass, and star formation rate grids from the particle distribution via
a cloud-in-cell algorithm. A series of real-space spherical tophat filters are
then applied to each grid, beginning from a radius of
$R_{\mathrm{max}}=20.34h^{-1}\mathrm{cMpc}$, decreasing in radius by a factor
of $\Delta R=1.1$ down to the size of the grid cell. These smoothed grids
represent the amount of mass, stars and star-formation within a sphere of
radius R centred on the cell. The maximum filter radius $R_{\mathrm{max}}$ was
chosen to be comparable to extrapolation of measurements of the mean free path
of UV photons in an ionised medium at $z<5$ (Songaila & Cowie, 2010; Worseck
et al., 2014), although the model does not include a physical prescription for
the mean free path (see Davies & Furlanetto (2022)). The interpretation of the
maximum radius as a mean free path is inconsistent with recent measurements
showing a rapid decrease at $z=6$ (Becker et al., 2021). However it is unclear
whether the cause of this decrease is driven by the structure of residual
neutral hydrogen, self-shielding regions on small scales, or the epoch of
reionisation being incomplete at $z=6$. Our maximum filter radius
$R_{\mathrm{max}}$ is meant to represent a maximum scale which a single source
can ionise, and is hence more analogous to the $z=5$ mean free path where we
expect a much more ionised medium.

If at any filter radius, a cell in the filtered stellar grid $m_{*}$ and gas
mass grid $M_{\mathrm{tot}}$ satisfies

$\zeta\frac{m_{*}}{M_{\mathrm{tot}}}>1,$ (2)

where

$\zeta=\frac{N_{\gamma}}{f_{b}(1-\frac{3}{4}Y_{\rm{He}})}$ (3)

$m_{*}=f_{10}\sum_{i}^{\mathrm{cell}}m_{i}\left(\frac{M_{\mathrm{halo},\it{i}}}{10^{10}M_{\odot}}\right)^{f_{\alpha}}$
(4)

Then the cell is considered ionised. The escape fraction normalisation
$f_{10}$ and halo mass scaling $f_{\alpha}$ are used to weight each stellar
particles’ mass $m_{i}$ and based on the total mass of the halo in which each
star particle resides $M_{\mathrm{halo},\it{i}}$ and represent the fraction of
ionising photons which escape the host halo into the IGM. In equation 3,
$N_{\gamma}=4000$, $f_{b}=\frac{\Omega_{b}}{\Omega_{m}}\approx 0.165$, and
$Y_{\mathrm{He}}=0.24$ such that $\zeta\approx 29564$ prior to the escape
fraction weighting. We set $f_{10}=0.25$ and $f_{\alpha}=0.5$ so that the
global reionization history generated from the stellar population is
consistent with electron optical depth measurements from Planck Collaboration
et al. (2020), and completes shortly after $z=6$, in agreement with Lyman
alpha forest measurements of the tail of reionization (Qin et al., 2021;
Keating et al., 2020; Nasir & D’Aloisio, 2020; Bosman et al., 2022). The halo
mass scaling of the escape fraction can be driven by different feedback
processes in galaxies which clear out columns of gas (Seiler et al., 2019),
allowing ionising photons to escape. While simulations more commonly predict
that on average, the escape fraction is smaller in larger halos (Ferrara &
Loeb, 2013; Paardekooper et al., 2015; Xu et al., 2016), $f_{\alpha}$ is not
well constrained by current observational data (Park et al., 2019; Mutch et
al., 2023). In the case of Astrid, the stellar mass and UV luminosity
functions sit $\sim 1\sigma$ below the mean of the observations at redshifts
$6<z<8$ (see Bird et al. (2021)), particularly at low stellar masses. This
stellar mass distribution implies a relatively high escape fraction which
scales positively with halo mass, in order to match observations of global
neutral fraction.

Figure 4: Ionised bubble radius versus the UV magnitude of the brightest
galaxies within them in Astrid-ES (left) and Astrid (right). Each point
represents one galaxy which is brightest in its HII region. Points with error
bars represent the binned mean and $95^{\mathrm{th}}$ percentile range in bins
of width 1. We present the snapshots in the left column of figure 3 where the
universe is $\approx 22\%$ ionised. Bubble radius is more closely correlated
with UV magnitude in the excursion-set model, and there are many bright
galaxies outside of ionised regions ($R_{\mathrm{bubble}}=0$) as shown in the
bottom panels, which show the fraction $f_{\mathrm{ntl}}$ of FOF groups
(including those which aren’t brightest in their region) within each bin that
exists in a neutral region.

Once the reionization state of the universe is calculated for a particular
snapshot, we find the specific ionising intensity $J_{21}$ for each fully
ionised cell, on the largest filter scale which satisfies the condition
specified in equation 2.

$J_{21}=\frac{(1+z)^{2}}{4\pi}\lambda_{\mathrm{mfp}}h\alpha_{\mathrm{UV}}\bar{\epsilon}\times
10^{-21}\;\mathrm{erg}\,\mathrm{cm}^{-2}\mathrm{Hz}^{-1}\mathrm{s}^{-1}\mathrm{sr}^{-1},$
(5)

where $h$ is the Planck constant, $\lambda_{\mathrm{mfp}}$ is the mean-free
path of ionising photons (here assumed to be the bubble radius and filter
scale, $r$), and the mean emissivity in the sphere $\bar{\epsilon}$ is given
by

$\bar{\epsilon}=\frac{N_{\gamma}}{\frac{4}{3}\pi r^{3}m_{p}}\frac{\delta
m_{*}}{\delta t}$ (6)

$\frac{\delta m_{*}}{\delta t}$ refers to the star formation rate within the
cell, smoothed on the same scale as the stellar mass grid and weighted by the
same escape fraction. Since the star formation rates of gas particles are very
bursty, we smooth the star formation rate over a fraction of the Hubble time

$\frac{\delta m_{*}}{\delta t}=\frac{m_{*}}{t_{*}H^{-1}(t)}$ (7)

Where we set $t_{*}=0.4$, and $H^{-1}(t)$ is the Hubble time. This smoothing
is necessary in order to function with the cooling model in MP-Gadget which
assumes ionisation equilibrium.

Local ionisation and heating rates can be calculated from $J_{21}$ assuming a
spectral slope $\alpha_{\mathrm{UV}}$ above the HI ionisation threshold
$\nu_{\mathrm{HI}}$ such that
$J_{\nu}=J_{21}\left(\frac{\nu}{\nu_{\mathrm{HI}}}\right)^{-\alpha_{\mathrm{UV}}}$,
by integrating over the source spectra and species cross section
$\sigma_{i}(\nu)$ above the ionisation threshold $\nu_{i}$ for each species
$i=\mathrm{\\{HI,HeI,HeII\\}}$.

$\Gamma_{i}=\int_{\nu_{i}}^{\infty}\frac{4\pi
J_{\nu}}{h_{\mathrm{p}}\nu}\sigma_{i}(\nu)d\nu\;\;\mathrm{s}^{-1},$ (8)
$g_{i}=\int_{\nu_{i}}^{\infty}\frac{4\pi
J_{\nu}}{\nu}(\nu-\nu_{i})\sigma_{i}(\nu)d\nu\;\;\mathrm{erg\,s}^{-1},$ (9)

where we use $\alpha_{\mathrm{UV}}=2$ in the following simulations. The
heating and ionisation rates are input to the existing cooling code (following
Katz et al. (1996) equations 23-32) which uses each $\Gamma_{i}$ to solve for
an equilibrium ionisation state at each timestep and tracks changes in
internal particle energy via the ionisation heating rates $g_{i}$ as well as
recombinations, collisions, Bremsstrahlung, and inverse Compton cooling. The
internal energy of each particle then feeds back onto processes such as gas
cooling and star formation. The ionising background is calculated from the
distribution of particles every $10$ Myr and assumed to be constant between
updates.

When a region is first ionised, it is heated to a high temperature $>10^{4}K$
depending on the background ionising spectrum and the speed of the ionisation
front (D’Aloisio et al., 2019). Since our timesteps between excursion set
calculations are too long to reliably calculate ionisation front speed, we set
the post-reionisation temperature of each particle to
$T_{\mathrm{reion}}=15000K$, corresponding to a scenario where ionisation
fronts are slowed significantly by unresolved absorbers. $T_{\mathrm{reion}}$,
the star formation rate timescale $t_{*}$, and the spectral slope of sources
$\alpha_{\mathrm{UV}}$ are chosen in order to be simultaneously consistent
with high-redshift IGM temperature measurements from Gaikwad et al. (2020) and
ionisation rate measurements (Wyithe & Bolton, 2011; Calverley et al., 2011;
D’Aloisio et al., 2018; Becker et al., 2021).

##  4 Results

We present reionization realizations from Astrid with the existing parametric
reionization model on a $248^{3}$ grid, and from a re-run of reionisation from
$z=20$ to $z=5.5$ labelled Astrid-ES, which includes the excursion-set model
performed on a $500^{3}$ grid and described in section 3. Both simulations
have a dark matter particle mass of $6.75\times 10^{6}h^{-1}M_{\odot}$ and a
gas particle mass of $1.267\times 10^{6}h^{-1}M_{\odot}$, with initial
conditions containing $2\times 5500^{3}$ particles within a cube of side
length $250h^{-1}\mathrm{Mpc}$. We show the galaxy stellar mass functions of
Astrid-ES at redshifts $z=8$, $z=7$ and $z=6$ in Figure 1, which are in
agreement with Astrid (Bird et al., 2021).

Figure 2 shows the global reionization history of Astrid and Astrid-ES,
compared with the Planck 2018 electron optical depth measurements (Planck
Collaboration et al., 2020) and observational constraints on $X_{\mathrm{HI}}$
from dark pixels (McGreer et al., 2015), Lyman alpha emitter fraction (Inoue
et al., 2018; Ouchi et al., 2018; Morales et al., 2021; Wold et al., 2022),
lyman alpha equivalent width (Hoag et al., 2019; Mason et al., 2019; Whitler
et al., 2020) and Quasar damping wings (Greig & Mesinger, 2017a; Bañados et
al., 2018; Davies et al., 2018; Greig et al., 2019; Wang et al., 2020). The
parametric model reionises earlier, with a midpoint of $z\approx 7.5$, chosen
to match the Planck Collaboration et al. (2020) electron optical depth
measurements. We also show that the volume-averaged ionisation rate
$\Gamma_{12}$ and the temperature at mean density $T_{0}$ of the excursion set
model are consistent with observations (Wyithe & Bolton, 2011; Calverley et
al., 2011; D’Aloisio et al., 2018; Becker et al., 2021; Gaikwad et al., 2020).
The late reionisation of the excursion-set model is consistent with inference
on large-scale Lyman alpha optical depth fluctuations (Qin et al., 2021;
Choudhury et al., 2021; Bosman et al., 2022) which suggest that the final
overlap stages of reionisation $x_{\mathrm{HI}}\lesssim 5\%$ should occur at
late times $z<5.6$.

###  4.1 Comparison of models

We first present a comparison between the reionization topologies of the
parametric reionization model in Astrid and the excursion-set algorithm in
Astrid-ES. Figure 3 shows the reionization topology of the simulation between
redshifts $z=8$ and $z=5.5$. When comparing these models, we select snapshots
from each simulation that have a fraction of ionised cells (those satisfying
the condition in equation 2) that are approximately equal in order to present
each simulation at a similar stage of cosmic reionization.

As expected, the parametric model proceeds outward from the dense cosmic
filaments, regardless of where stars form in the simulation. In the excursion-
set algorithm, spherical ionised regions form around areas with high stellar
mass. The highly spherical nature of the ionised regions in the excursion set
arises from our choice of escape fraction which scales positively with halo
mass, and is consistent with other reionization models driven by larger
sources (e.g. McQuinn et al. (2007); Geil et al. (2016)). To show the
connection between star formation and reionization quantitatively, we examine
the cross-correlation between the density, stellar mass, and galaxy fields in
Figure 5. The Astrid-ES ionisation field has a stronger negative correlation
with the stellar mass, reaching $\sim-0.75$ at scales of
$1h^{-1}\mathrm{Mpc}$, compared to $\sim-0.6$ in Astrid. While the cross-
correlation between all Friends-of-Friends groups and the ionization field is
stronger in Astrid, this includes many low-mass structures which better trace
the density field than the stellar mass field. Placing a UV Magnitude cut on
the objects at $M_{UV}<-18$ returns an even stronger correlation in both
simulations due to galaxy bias, however in Astrid-ES the increase is sharper,
due to both the direct dependence of reionization history on stellar mass as
well as the positive escape fraction weighting with halo mass, which increases
the effect of the largest structures on reionization. The UV magnitude is
calculated per Friends-of-Friends group from its star formation rate (Stringer
et al., 2011):

$\mathrm{M_{UV}}=-2.5\mathrm{log}_{10}(\psi)-18.45,$ (10)

where $\psi$ is the star formation rate in solar masses per year. We can also
see the connection between galaxies and the ionised structure by comparing the
size of the ionised regions to the UV magnitude of the brightest galaxies
within them. We measure the relationship between galaxy UV magnitude and HII
region size in Figure 4, where the radius of each ionised region is calculated
by taking the average of the distance in $5000$ random directions from a
galaxy to an ionisation boundary $X_{\mathrm{HI}}>0.5$, and eliminating
galaxies that are not the brightest (lowest $\mathrm{M_{UV}}$) within their
radius.

Figure 5: The cross-correlation functions between the $x_{\mathrm{HI}}$ field
and other fields in each simulation at a neutral fraction of
$Q_{\mathrm{HI}}\approx 0.22$, using a $0.5\mathrm{Mpc}$ resolution. We show
the cross-correlation of neutral fraction with density (orange), stellar mass
(blue), All FOF groups (green) and FOF groups with a UV magnitude below -18
(red). Cross-correlations in Astrid-ES are shown with solid lines and those in
Astrid are showin with dased lines. On all scales, the Astrid-ES ionization
field is more strongly correlated with the stellar mass and brighter galaxies,
whereas in Astrid there is has a slightly stronger correlation with the
density field.

Astrid-ES has a much tighter correlation between ionised region size and the
UV magnitude of their central galaxies. Astrid also has many bright galaxies
outside of HII regions ($R=0$), including $\sim 10\%$ of galaxies at
$\mathrm{M_{UV}}=-20$, up to $\sim 20\%$ of galaxies at $\mathrm{M_{UV}}=-17$.
This is an effect of the parametric model not directly depending on galaxy
properties, such that relatively bright galaxies outside of dense regions will
not ionise their surroundings. However, it should be noted that our choice of
escape fraction which scales positively with halo mass will also contribute to
the tight correlation between $\mathrm{M_{UV}}$ and HII region size, since
brighter galaxies will contribute a greater proportion of the ionising photons
within each HII region.

One of the key observables during reionization is the 21cm power spectrum.
This is calculated from the 21cm brightness temperature, which depends on the
neutral fraction $x_{\mathrm{HI}}$ and overdensity $\delta$ grids. Here we use
(Furlanetto et al., 2006):

$\delta
T_{\mathrm{b,21}}=27\mathrm{K}\left(\frac{\Omega_{b}h^{2}}{0.023}\right)\left(\frac{0.15}{\Omega_{m}h^{2}}\frac{1+z}{10}\right)^{1/2}\left(1-\frac{T_{\mathrm{CMB}}}{T_{\mathrm{s}}}\right)(1+\delta)x_{\mathrm{HI}}$
(11)

Figure 6: Top Panel: Dimensional power spectrum $\overline{\delta
T}_{\mathrm{b,21}}(z)^{2}\Delta_{21}^{2}$ in Astrid-ES (solid lines) and
Astrid (dotted lines). Snapshots for each simulation excursion-set model are
chosen at matching stages of reionisation, as in figure 3, when the universe
is $\approx 22\%$, $50\%$, and $95\%$ ionised. Bottom Panel: Ratio of the
power spectra in the excursion set and parametric models at the matched
snapshots.

We assumed the spin temperature of the gas is saturated
$T_{\mathrm{s}}>>T_{\mathrm{CMB}}$, which is valid in the scenario where X-ray
heating from early quasars occurs prior to the bulk of reionization
(Furlanetto et al., 2006). This formula also ignores redshift space
distortions due to the peculiar velocities of the gas.

The spherically averaged 21cm power spectrum is then calculated from the
Fourier transform of the brightness temperature grid on a scale $k$ (e.g (Lidz
et al., 2008)):

$\Delta_{21}^{2}(k,z)=\frac{k^{3}}{2\pi
V}\left<|\delta_{21}(k,z)|^{2}_{k}\right>,$ (12)

where

$\delta_{21}(k,z)=\frac{\delta T_{\mathrm{b,21}}(k,z)}{\overline{\delta
T}_{\mathrm{b,21}}(z)}-1$ (13)

and $\overline{\delta T}_{\mathrm{b,21}}(z)$ is the average brightness
temperature in the simulation at redshift $z$, and $V$ is the simulation
volume. Due to the difficulty in measuring the 21cm signal from the EoR, only
upper limits have been placed on the 21cm power spectrum so far (Trott et al.,
2020; Mertens et al., 2020; Patil et al., 2017). However, the 21cm power
spectrum contains a lot of information on the scale of ionised and neutral
regions during reionization, and will be one of the most important observables
used to constrain our models in the near future (Greig et al., 2021a; Greig et
al., 2021b).

We show the resulting power spectra in Figure 6. There are two main
differences in the power-spectra between the excursion-set and parametric
reionisation models. Firstly, the excursion set contains fewer and larger
ionised regions at a fixed $X_{\mathrm{HI}}$, which increases large scale
power, especially toward the end of reionisation where the last remaining
neutral regions are much larger. This is caused by a combination of the escape
fraction scaling and the clustering of stellar particles compared to the
density field. Secondly, the more spherical ionised regions change the shape
of the power-spectrum, showing a slight increase in power on scales
$0.1\lesssim k\lesssim 0.5h\mathrm{Mpc}^{-1}$, corresponding to the size of
the ionised regions $5\lesssim R\lesssim 30h^{-1}\mathrm{Mpc}$. The more
asymmetric shape of ionised regions in the parametric model spreads this power
over a larger range of $k$.

Figure 7: $250\times 250\times 5h^{-1}\mathrm{Mpc}$ slices ofAstrid-ES between
redshifts $z=8$ and $z=5.5$. Left: Matter overdensity, Centre: Neutral
hydrogen fraction. Right: Gas temperature. Spherical ionised regions form and
expand from the dense filaments of the cosmic web, heating the gas to
$T=1.5\times 10^{4}\mathrm{K}$. Figure 8: HII region radius versus the UV
magnitude of the brightest galaxy within, at redshifts $z=8.5$, $z=7$, and
$z=6.5$ in Astrid-ES. Points with error bars represent the mean and
$95^{\mathrm{th}}$ percentile range in radius, within bins of width 1. The red
dashed line shows a power law fit to this relation
$\mathrm{log}(R)=-0.314M_{\mathrm{UV}}-2.550\mathrm{log}(1+z)+7.408$ using the
snapshots at $z=8.0$, $z=7.5$, $z=7.0$, and $z=6.5$.

###  4.2 Excursion-Set Results

In this section we present more detailed results for Astrid-ES. Figure 7 shows
the reionization history of a $250\times 250\times 5h^{-1}\mathrm{Mpc}$ slice
of the simulation, showing the neutral fraction, matter overdensity, and gas
temperature at a resolution of $0.5h^{-1}\mathrm{Mpc}$. The first ionised
regions grow in the overdense regions of the simulation where most of the
stars form. Reionization proceeds outward from these regions as new bubbles
form, expand and overlap to cover the entire volume shortly after redshift
$z=6$. As the ionisation front passes, gas is heated to $15000\mathrm{K}$ (see
section 3). While we stop the simulation at the last stages of reionisation,
we can see some of the earliest ionised regions begin to cool toward a
temperature closely correlated with their density (Upton Sanderbeck et al.,
2016).

We show the correlation between HII region size and UV magnitude at redshifts
$z=8$, $z=7$, and $z=6.5$ for Astrid-ES in Figure 8. We find that a power-law
fit of $\mathrm{log}(R)=-0.314M_{\mathrm{UV}}-2.550\mathrm{log}(1+z)+7.408$
well describes the data pre-overlap, with a standard deviation of
$\sigma_{R}<0.15\mathrm{dex}$ across all mass bins and redshifts. Since
brightest galaxies within a HII region are selected by excluding galaxies
which have a neighbor within $\mathrm{R}$ with a lower $M_{\mathrm{UV}}$, the
fit becomes less accurate during overlap, where HII regions are less
spherical, and dimmer galaxies within a larger HII structure are included in
the sample.

We present the dimensional 21cm power spectrum $\overline{\delta
T}_{\mathrm{b,21}}(z)^{2}\Delta_{21}^{2}$ throughout reionization in our
simulation in Figure 9. While the universe is entirely neutral, the shape of
the 21cm power spectrum resembles the matter power spectrum. Once the first
ionised regions start to form, power increases on large scales $k\lesssim
1h\mathrm{Mpc}^{-1}$, corresponding to the size of the bubbles. Once the
universe is mostly reionised, the power on all scales decreases as the 21cm
signal from atomic hydrogen disappears. The qualitative shape of the power
spectrum is similar to many reionization models (Lidz et al., 2008), but the
increase in power on large scales, which results from the distribution in HII
region sizes, will take a different form depending on factors such as escape
fraction dependencies and feedback processes. This means that predictions for
the 21cm power spectrum can be useful in ruling out certain reionization
models (Geil et al., 2016; Greig et al., 2021a; Greig et al., 2021b).

Figure 9: The dimensional power spectrum $\overline{\delta
T}_{\mathrm{b,21}}(z)^{2}\Delta_{21}^{2}$ in Astrid-ES between redshifts $z=8$
and $z=6$. The power spectrum shape initially follows the density field, power
on large scales then increases as ionised regions form and grow. Once
reionization nears completion, power on all scales decreases as the 21cm
brightness temperature approaches zero.

###  4.3 Reionization Feedback

The photo-ionisation of intergalactic hydrogen is accompanied by rapid
heating, as the excess energy of ionising photons is converted to kinetic
energy in the IGM (Miralda-Escudé & Rees, 1994; Hui & Gnedin, 1997). This
heating will affect star-formation in lower mass galaxies, as gas from the IGM
is less likely to cool toward the centre of a halo (Couchman & Rees, 1986;
Sobacchi & Mesinger, 2013a). Other works have found that while reionization
feedback does not greatly affect the progress of reionization compared to
other feedback mechanisms, it significantly decreases star formation in small
halos over time (Sobacchi & Mesinger, 2013b; Mutch et al., 2016; Hutter et
al., 2021).

Figure 10: Top: Mean gas mass (left), stellar mass (middle) and star formation
rate (right) in Star formation rate distributions in Astrid-ES between
redshifts $5<z<9$. The mean SFR for halos within large-scale neutral regions
are shown as dotted lines, and for those within ionised within ionised regions
as solid lines. The distributions are calculated on halo mass bins within
$30\%$ of $2\times 10^{9}$, $5\times 10^{9}$, $1\times 10^{10}$, and $5\times
10^{10}$ solar masses. Bottom: Ratio of each halo property in neutral regions
to ionised regions, at the same masses and redshifts.

In Astrid, Bird et al. (2021) found up to a $50\%$ reduction in star formation
in within $2\times 10^{9}M_{\odot}$ halos which reside in large-scale ionised
regions compared to halos in neutral regions, toward the end of reionization.
On the other hand, large halos had similar SFR regardless of whether or not
they are in an ionised region. Since the UV background in Astrid-ES is
included on-the-fly and results directly from the sources in the simulation,
we repeat and expand on this analysis to determine if these correlations
persist in a more self-consistent reionisation model. The UV background given
by the excursion set (equations 8 \- 9) is applied to the gas particles in an
ionised cell, as well as the sharp temperature increase of
$T_{\mathrm{reion}}=15000K$ at the time of reionisation. This heating can
reduce the amount of cold gas available for star formation, particularly in
less massive halos.

We examine the same results in Astrid-ES in Figure 10, additionally including
gas mass and stellar mass in halos. The correlation between star formation
rates and ionisation state is similar to Astrid in the lowest masses, but is
weaker toward higher masses, and reversed in the highest mass bin with $\sim
5\times 10^{10}M_{\odot}$ halos having mean star formation rates up to $75\%$
lower in neutral regions. This is due to the direct connection between the
star formation in these groups and reionization structure in the excursion-set
model, and is likely exaggerated by our increasing escape fraction toward
larger halos. In order for a halo of this size to exist in a neutral region,
it must have a relatively low star formation, otherwise it would have ionised
its surroundings. The median stellar mass follows a similar pattern, being
slightly higher in lower mass halos, and much lower in high-mass halos within
neutral regions. The median gas mass shows the opposite trend at low masses,
being lower in neutral regions at all masses.

The differences between these properties in ionised and neutral regions of the
universe can be driven by a number of effects. The heated gas left behind
after reionization may not cool as quickly onto halos or form stars. This
means that existing halos will form stars more slowly, and newly formed halos
will contain fewer stars than those forming in neutral regions. This drives an
anti-correlation between star formation and ionisation for smaller halos with
shallower potential wells, as well as at later times when halos have spent
more time on average surrounded by a heated IGM. Highly star forming galaxies
are also more likely to ionise their surroundings, and therefore exist within
ionised regions, causing a strong correlation between star formation and
ionisation at higher masses while reionisation is ongoing.

Figure 11: Mean halo gas mass (Top), stellar mass fraction (Middle) and
specific star formation rate (Bottom) in halos of mass
$M_{\mathrm{halo}}\approx 2\times 10^{9}h^{-1}M_{\odot}$ which ionise around
$z=10$ (blue line), $z=8$ (orange line), and $z=6$ (green line). The
reionisation timing of each group is marked with a vertical dotted line of the
corresponding colour. While there are correlations between halo properties and
reionisation timing at fixed halo mass, there is no sudden change or
acceleration in each property at the redshift of reionisation, meaning that
reionisation feedback likely plays a subdominant role in the evolution of
these properties.

In order to examine these effects more closely, we focus on the lowest mass
bin $M_{\mathrm{halo}}\approx 2.0\times 10^{9}$ solar masses, which we expect
to experience the largest effect from reionisation feedback. We plot the
median gas mass fraction, stellar mass fraction, and specific star formation
rate for groups of halos which ionised at redshifts $z\approx 10$, $z\approx
8$, and $z\approx 6$ in figure 11. The $z_{\mathrm{reion}}\approx 10$ halos
contain up to $\sim 5\%$ more gas, $\sim 5\%$ less stellar mass, and $\sim
10\%$ less star formation. Suggesting that low mass halos which ionise earlier
form stars from their gas reservoirs more slowly. However, these correlations
persist both before and after the IGM surrounding the halos ionise, and do not
show a sudden change or acceleration post-ionisation. This suggests that the
less efficient star formation does not result from reionisation feedback, and
that the heating of the IGM on large scales by reionisation does not greatly
affect the cooling of gas onto halos within the timeframe of our simulation.
Instead, it is likely that the differences come from environmental effects,
where the regions that form the first ionising sources at higher masses form
stars more slowly in lower mass halos. Since we expect reionization feedback
to have a cumulative effect on galaxies, running our simulation to lower
redshift may reveal a larger effect as halos are exposed to the UV background
for longer. It is also worth mentioning that the reionisation temperature
$T_{\mathrm{reion}}=15000K$ chosen to match observations of IGM temperature at
mean density is low compared to many models (D’Aloisio et al., 2018), a higher
reionisation temperature may also show a stronger effect on small halos.
Comparing similar halos over a longer period of time, or matching halos across
simulations with different reionisation histories (such as Astrid and Astrid-
ES) would provide a more thorough study of feedback and environmental effects,
however we leave this analysis to future work.

##  5 Conclusion

In this paper we have introduced the addition of an excursion-set reionization
model based on 21cmFAST into the hydrodynamic code MP-Gadget. This model is
efficient enough to run on large simulations such as Astrid without
significantly impacting the computational cost. We have performed a re-run of
the epoch of reionization $5.5<z<20$ in Astrid labeled Astrid-ES. The global
history, ionisation rates and mean temperatures in Astrid-ES are consistent
with high-redshift quasar spectra and CMB measurements. We find that the
excursion-set model correlates more strongly with the ionising sources in the
simulation when compared against against the parametric model applied in the
Astrid production run, in which the redshift of reionisation is calculated
based on the initial mass density grid. With our choice of escape fraction,
reionization in the excursion-set model contains fewer, larger ionised
regions, which is reflected in the 21cm power spectrum. The model also
includes the effect of reionization feedback, where an ionising background
dependent on the galaxy star formation rate is applied to each particle. While
there are significant correlations between halo properties and reionisation
timing, they are likely driven by environmental effects, rather than by
reionisation feedback, since they exist before and after the reionisation of
the halo.

Future large simulations run with MP-Gadget will have the option to enable the
excursion-set model, allowing for more accurate reionization topology and for
further studies of the connection between galaxy properties and reionization.
With a larger simulation, analyses of 21cm cross correlations (Kubota et al.,
2018; Vrbanec et al., 2020; Cox et al., 2022) or predictions for 21cm images
from the SKA (Ghara & Choudhury, 2020; Davies et al., 2021) can be performed
using a more accurate reionization model, producing more accurate predictions
for the sensitivity of future observations to the parameters of reionization.

The source-driven reionization model will also allow the addition of further
physics to the model, so that a wider range of observations can be predicted
and calibrated against. The size of simulations such as Astrid, makes them
ideal for studies of large scale X-ray heating during cosmic dawn (Greig &
Mesinger, 2017b). A future work will include the X-ray heating of the IGM from
galaxies, allowing us to make predictions for the spin temperature, $T_{s}$,
include the effects on the EoR of early photoheating, and extend the range of
predictions to include 21cm measurements toward higher redshifts, such as the
data obtained by global 21cm experiments (Bowman et al., 2018; Singh et al.,
2022).

The large box size and black hole formation model within Astrid also makes it
ideal for the study of the contribution of quasars to reionisation. The
excursion set model can be extended to include quasars (e.g: Qin et al.
(2017)). While most works find that quasars have a subdominant contribution to
hydrogen reionisation (Trebitsch et al., 2023), their inclusion into the model
would allow us to study the effect of rare, bright sources on the topology of
reionisation (Chardin et al., 2017; Eide et al., 2020), and early quasar
proximity zones (Chen & Gnedin, 2021). This would also provide a smooth
coupling to the existing helium reionisation model within MP-GADGET, which
will be useful to study the post-reionisation UV background.

## Acknowledgements

JD acknowledges support from the Ministry of Universities and Research (MUR)
through the PRO3 project “Data Science methods for Multi-Messenger
Astrophysics & Multi-Survey Cosmology”. This research was supported by the
Australian Research Council Centre of Excellence for All Sky Astrophysics in 3
Dimensions (ASTRO 3D), through project number CE170100013. This work utilised
the OzSTAR national facility at Swinburne University of Technology. OzSTAR is
funded by Swinburne University of Technology and the National Collaborative
Research Infrastructure Strategy (NCRIS). SB acknowledges funding support from
NASA-80NSSC21K1840. The authors acknowledge the Frontera computing project at
the Texas Advanced Computing Center (TACC) for providing HPC and storage
resources that have contributed to the research results reported within this
paper. Frontera is made possible by National Science Foundation award
OAC-1818253. URL: http://www.tacc.utexas.edu.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author. E-mail<EMAIL_ADDRESS>
## References

  * Bañados et al. (2018) Bañados E., et al., 2018, ApJ, 856, L25
  * Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, The Astrophysical Journal, 776, 81
  * Becker et al. (2021) Becker G. D., D’Aloisio A., Christenson H. M., Zhu Y., Worseck G., Bolton J. S., 2021, MNRAS, 508, 1853
  * Bird et al. (2021) Bird S., Ni Y., Di Matteo T., Croft R., Feng Y., Chen N., 2021, arXiv e-prints, p. arXiv:2111.01160
  * Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Bowman et al. (2018) Bowman J. D., Rogers A. E. E., Monsalve R. A., Mozdzen T. J., Mahesh N., 2018, Nature, 555, 67
  * Calverley et al. (2011) Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543
  * Chardin et al. (2017) Chardin J., Puchwein E., Haehnelt M. G., 2017, Monthly Notices of the Royal Astronomical Society, 465, 3429
  * Chen & Gnedin (2021) Chen H., Gnedin N. Y., 2021, ApJ, 911, 60
  * Chen et al. (2021) Chen N., Ni Y., Tremmel M., Di Matteo T., Bird S., DeGraf C., Feng Y., 2021, arXiv e-prints, p. arXiv:2104.00021
  * Choudhury & Paranjape (2018) Choudhury T. R., Paranjape A., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3821
  * Choudhury et al. (2021) Choudhury T. R., Paranjape A., Bosman S. E. I., 2021, MNRAS, 501, 5782
  * Couchman & Rees (1986) Couchman H. M. P., Rees M. J., 1986, MNRAS, 221, 53
  * Cox et al. (2022) Cox T. A., Jacobs D. C., Murray S. G., 2022, MNRAS, 512, 792
  * D’Aloisio et al. (2018) D’Aloisio A., McQuinn M., Davies F. B., Furlanetto S. R., 2018, Monthly Notices of the Royal Astronomical Society, 473, 560
  * D’Aloisio et al. (2019) D’Aloisio A., McQuinn M., Maupin O., Davies F. B., Trac H., Fuller S., Upton Sanderbeck P. R., 2019, Astrophysical Journal, 874, 154
  * Davies & Furlanetto (2022) Davies F. B., Furlanetto S. R., 2022, MNRAS, 514, 1302
  * Davies et al. (2018) Davies F. B., et al., 2018, ApJ, 864, 142
  * Davies et al. (2021) Davies J. E., Croft R. A. C., Di-Matteo T., Greig B., Feng Y., Wyithe J. S. B., 2021, MNRAS, 501, 146
  * Dayal & Ferrara (2018) Dayal P., Ferrara A., 2018, Phys. Rep., 780, 1
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  * Duncan et al. (2014) Duncan K., et al., 2014, MNRAS, 444, 2960
  * Eide et al. (2020) Eide M. B., Ciardi B., Graziani L., Busch P., Feng Y., Di Matteo T., 2020, MNRAS, 498, 6083
  * Faucher-Giguère (2020) Faucher-Giguère C.-A., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1614
  * Feng et al. (2016) Feng Y., Di-Matteo T., Croft R. A., Bird S., Battaglia N., Wilkins S., 2016, Monthly Notices of the Royal Astronomical Society, 455, 2778
  * Ferrara & Loeb (2013) Ferrara A., Loeb A., 2013, MNRAS, 431, 2826
  * Furlanetto et al. (2006) Furlanetto S. R., Oh S. P., Briggs F. H., 2006, Physics Reports, 433, 181
  * Gaikwad et al. (2020) Gaikwad P., et al., 2020, MNRAS, 494, 5091
  * Geil et al. (2016) Geil P. M., Mutch S. J., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, Monthly Notices of the Royal Astronomical Society, 462, 804
  * Ghara & Choudhury (2020) Ghara R., Choudhury T. R., 2020, Monthly Notices of the Royal Astronomical Society, 496, 739
  * González et al. (2011) González V., Labbé I., Bouwens R. J., Illingworth G., Franx M., Kriek M., 2011, ApJ, 735, L34
  * Grazian et al. (2015) Grazian A., et al., 2015, A&A, 575, A96
  * Greig & Mesinger (2017a) Greig B., Mesinger A., 2017a, Monthly Notices of the Royal Astronomical Society, 465, 4838
  * Greig & Mesinger (2017b) Greig B., Mesinger A., 2017b, Monthly Notices of the Royal Astronomical Society, 472, 2651
  * Greig et al. (2019) Greig B., Mesinger A., Bañados E., 2019, Monthly Notices of the Royal Astronomical Society, 484, 5094
  * Greig et al. (2021a) Greig B., Trott C. M., Barry N., Mutch S. J., Pindor B., Webster R. L., Wyithe J. S. B., 2021a, Monthly Notices of the Royal Astronomical Society, 500, 5322
  * Greig et al. (2021b) Greig B., et al., 2021b, Monthly Notices of the Royal Astronomical Society, 501, 1
  * Hassan et al. (2022) Hassan S., Davé R., McQuinn M., Somerville R. S., Keating L. C., Anglés-Alcázar D., Villaescusa-Navarro F., Spergel D. N., 2022, ApJ, 931, 62
  * Hoag et al. (2019) Hoag A., et al., 2019, ApJ, 878, 12
  * Hopkins (2013) Hopkins P. F., 2013, Monthly Notices of the Royal Astronomical Society, 428, 2840
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, Monthly Notices of the Royal Astronomical Society, 292, 27
  * Hutter (2018) Hutter A., 2018, MNRAS, 477, 1549
  * Hutter et al. (2021) Hutter A., Dayal P., Yepes G., Gottlöber S., Legrand L., Ucci G., 2021, MNRAS, 503, 3698
  * Inoue et al. (2018) Inoue A. K., et al., 2018, PASJ, 70, 55
  * Kannan et al. (2021) Kannan R., Garaldi E., Smith A., Pakmor R., Springel V., Vogelsberger M., Hernquist L., 2021, arXiv e-prints, p. arXiv:2110.00584
  * Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, Astrophysical Journal, Supplement, 105, 19
  * Keating et al. (2018) Keating L. C., Puchwein E., Haehnelt M. G., 2018, MNRAS, 477, 5501
  * Keating et al. (2020) Keating L. C., Weinberger L. H., Kulkarni G., Haehnelt M. G., Chardin J., Aubert D., 2020, MNRAS, 491, 1736
  * Krumholz & Gnedin (2011) Krumholz M. R., Gnedin N. Y., 2011, Astrophysical Journal, 729, 36
  * Kubota et al. (2018) Kubota K., Yoshiura S., Takahashi K., Hasegawa K., Yajima H., Ouchi M., Pindor B., Webster R. L., 2018, MNRAS, 479, 2754
  * Lewis et al. (2022) Lewis J. S. W., et al., 2022, MNRAS, 516, 3389
  * Lidz et al. (2008) Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Hernquist L., 2008, Astrophysical Journal, 680, 962
  * Mason et al. (2019) Mason C. A., et al., 2019, MNRAS, 485, 3947
  * McGreer et al. (2015) McGreer I. D., Mesinger A., D’Odorico V., 2015, MNRAS, 447, 499
  * McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
  * Mertens et al. (2020) Mertens F. G., et al., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1662
  * Mesinger & Furlanetto (2007) Mesinger A., Furlanetto S., 2007, Astrophysical Journal, 669, 663
  * Mesinger et al. (2011) Mesinger A., Furlanetto S., Cen R., 2011, Monthly Notices of the Royal Astronomical Society, 411, 955
  * Mihalas & Mihalas (1984) Mihalas D., Mihalas B. W., 1984, Foundations of radiation hydrodynamics
  * Miralda-Escudé & Rees (1994) Miralda-Escudé J., Rees M. J., 1994, Monthly Notices of the Royal Astronomical Society, 266, 343
  * Molaro et al. (2019) Molaro M., Davé R., Hassan S., Santos M. G., Finlator K., 2019, Monthly Notices of the Royal Astronomical Society, 489, 5594
  * Morales et al. (2021) Morales A. M., Mason C. A., Bruton S., Gronke M., Haardt F., Scarlata C., 2021, ApJ, 919, 120
  * Mutch et al. (2016) Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, Monthly Notices of the Royal Astronomical Society, 462, 250
  * Mutch et al. (2023) Mutch S. J., Greig B., Qin Y., Poole G. B., Wyithe J. S. B., 2023, arXiv e-prints, p. arXiv:2303.07378
  * Nasir & D’Aloisio (2020) Nasir F., D’Aloisio A., 2020, MNRAS, 494, 3080
  * Ni et al. (2022) Ni Y., et al., 2022, MNRAS, 513, 670
  * Oñorbe et al. (2019) Oñorbe J., Davies F. B., Lukić Z., Hennawi J. F., Sorini D., 2019, Monthly Notices of the Royal Astronomical Society, 486, 4075
  * Ocvirk et al. (2020) Ocvirk P., et al., 2020, Monthly Notices of the Royal Astronomical Society, 496, 4087
  * Okamoto et al. (2010) Okamoto T., Frenk C. S., Jenkins A., Theuns T., 2010, Monthly Notices of the Royal Astronomical Society, 406, 208
  * Ouchi et al. (2018) Ouchi M., et al., 2018, PASJ, 70, S13
  * Paardekooper et al. (2015) Paardekooper J.-P., Khochfar S., Dalla Vecchia C., 2015, MNRAS, 451, 2544
  * Park et al. (2019) Park J., Mesinger A., Greig B., Gillet N., 2019, MNRAS, 484, 933
  * Patil et al. (2017) Patil A. H., et al., 2017, Astrophysical Journal, 838, 65
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, Astronomy & Astrophyscis, 641, A6
  * Puchwein et al. (2023) Puchwein E., et al., 2023, MNRAS, 519, 6162
  * Qin et al. (2017) Qin Y., et al., 2017, Monthly Notices of the Royal Astronomical Society, 472, 2009
  * Qin et al. (2021) Qin Y., Mesinger A., Bosman S. E. I., Viel M., 2021, MNRAS, 506, 2390
  * Read et al. (2010) Read J. I., Hayfield T., Agertz O., 2010, Monthly Notices of the Royal Astronomical Society, 405, 1513
  * Seiler et al. (2019) Seiler J., Hutter A., Sinha M., Croton D., 2019, MNRAS, 487, 5739
  * Singh et al. (2022) Singh S., et al., 2022, Nature Astronomy,
  * Sobacchi & Mesinger (2013a) Sobacchi E., Mesinger A., 2013a, MNRAS, 432, L51
  * Sobacchi & Mesinger (2013b) Sobacchi E., Mesinger A., 2013b, Monthly Notices of the Royal Astronomical Society, 432, 3340
  * Song et al. (2016) Song M., et al., 2016, ApJ, 825, 5
  * Songaila & Cowie (2010) Songaila A., Cowie L. L., 2010, Astrophysical Journal, 721, 1448
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, Monthly Notices of the Royal Astronomical Society, 339, 289
  * Stefanon et al. (2017) Stefanon M., Bouwens R. J., Labbé I., Muzzin A., Marchesini D., Oesch P., Gonzalez V., 2017, ApJ, 843, 36
  * Stringer et al. (2011) Stringer M., Cole S., Frenk C. S., Stark D. P., 2011, Monthly Notices of the Royal Astronomical Society, 414, 1927
  * Trebitsch et al. (2023) Trebitsch M., Hutter A., Dayal P., Gottlöber S., Legrand L., Yepes G., 2023, MNRAS, 518, 3576
  * Trott et al. (2020) Trott C. M., et al., 2020, Monthly Notices of the Royal Astronomical Society, 493, 4711
  * Upton Sanderbeck & Bird (2020) Upton Sanderbeck P., Bird S., 2020, Monthly Notices of the Royal Astronomical Society, 496, 4372
  * Upton Sanderbeck et al. (2016) Upton Sanderbeck P. R., D’Aloisio A., McQuinn M. J., 2016, Monthly Notices of the Royal Astronomical Society, 460, 1885
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, Monthly Notices of the Royal Astronomical Society, 436, 3031
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, Monthly Notices of the Royal Astronomical Society, 444, 1518
  * Vrbanec et al. (2020) Vrbanec D., Ciardi B., Jelić V., Jensen H., Iliev I. T., Mellema G., Zaroubi S., 2020, MNRAS, 492, 4952
  * Wang et al. (2020) Wang F., et al., 2020, ApJ, 896, 23
  * Whitler et al. (2020) Whitler L. R., Mason C. A., Ren K., Dijkstra M., Mesinger A., Pentericci L., Trenti M., Treu T., 2020, MNRAS, 495, 3602
  * Wold et al. (2022) Wold I. G. B., et al., 2022, ApJ, 927, 36
  * Worseck et al. (2014) Worseck G., et al., 2014, MNRAS, 445, 1745
  * Wyithe & Bolton (2011) Wyithe J. S. B., Bolton J. S., 2011, MNRAS, 412, 1926
  * Xu et al. (2016) Xu H., Wise J. H., Norman M. L., Ahn K., O’Shea B. W., 2016, ApJ, 833, 84