11institutetext: Scuola Normale Superiore, Piazza dei Cavalieri 7, 56125 Pisa,
PI, Italy

# The importance of stochasticity in determining galaxy emissivities and UV
LFs during cosmic dawn and reionization

Ivan Nikolić  Andrei Mesinger  James E. Davies  David Prelogović

The stochastic nature of star formation and photon propagation in high-
redshift galaxies can result in sizable galaxy-to-galaxy scatter in their
properties. Ignoring this scatter by assuming mean quantities can bias
estimates of their emissivity and corresponding observables. We construct a
flexible, semi-empirical model, sampling scatter around the following mean
relations: (i) the conditional halo mass function (CHMF); (ii) the stellar-to-
halo mass relation (SHMR); (iii) galaxy star formation main sequence (SFMS);
(iv) fundamental metallicity relation (FMR); (v) conditional intrinsic
luminosity; and (vi) photon escape fraction. In our fiducial model, ignoring
scatter in these galaxy properties overestimates the duration of the EoR,
delaying its completion by $\Delta z\sim$ 1–2. We quantify the relative
importance of each of the above sources of scatter in determining the
ionizing, soft-band X-ray and Lyman Werner (LW) emissivities as a function of
scale and redshift. We find that scatter around the SFMS is important for all
bands, especially at the highest redshifts where the emissivity is dominated
by the faintest, most ”bursty” galaxies. Ignoring this scatter would
underestimate the mean emissivity and its standard deviation computed over 5
cMpc regions by factors of up to $\sim$2–10 at $5\lesssim z\lesssim 15$.
Scatter around the X-ray luminosity to star formation rate and metallicity
relation is important for determining X-ray emissivity, accounting for roughly
half of its mean and standard deviation. The importance of scatter in the
ionizing escape fraction depends on its functional form, while scatter around
the SHMR contributes at the level of $\sim$10–20%. Other sources of scatter
have a negligible contribution to the emissivities. Although scatter does
flatten the UV luminosity functions, shifting the bright end by 1–2
magnitudes, the level of scatter in our fiducial model is insufficient to
fully explain recent estimates from JWST photometry (consistent with previous
studies). We conclude that models of the EoR should account for the burstiness
of star formation, while models for the cosmic 21cm signal should additionally
account for scatter in intrinsic X-ray production.

###### Key Words.:

Galaxies: high-redshift – intergalactic medium – Cosmology: diffuse radiation
– dark ages, reionization, first stars – X-rays: diffuse background

##  1 Introduction

The Universe underwent dramatic changes during the first billion years.
Following cosmic recombination, the Universe was cold, dark and fairly empty.
During the Cosmic Dawn (CD) when the first galaxies formed, their ultraviolet
(UV) and X-ray radiation spread out, heating and ionizing the intergalactic
medium (IGM). This culminated in the final major phase change of our Universe:
the epoch of reionization (EoR; see for example reviews in Zaroubi 2013,
Mesinger 2016 and Dayal & Ferrara 2018).

Understanding how the first galaxies heated and ionized the Universe requires
modeling their UV and X-ray emission, and constraining these models with data
(e.g. Qin et al., 2021; Abdurashidova et al., 2022). The emission of any
single galaxy is highly variable, depending on the time evolution of star
formation, feedback, and geometry of interstellar absorption (e.g. Tacchella
et al., 2016; Barrow et al., 2017; Lovell et al., 2021; Pallottini & Ferrara,
2023). These processes are not known from first principles, and are extremely
challenging to simulate for a single galaxy, let alone for a cosmological
sample of galaxies.

Luckily, the relevant cosmic radiation fields are sourced by the combined
contribution from many galaxies, which allows us to take advantage of the
Central Limit Theorem and use only average scaling relations to connect galaxy
properties to their host dark matter halos (whose abundances and evolution are
reasonably well-known). This is the general approach taken by many analytic,
semi-numerical and numerical models of the EoR/CD (e.g. Haiman et al., 2000;
Ciardi et al., 2003; Furlanetto et al., 2004; Mesinger et al., 2011; Holzbauer
& Furlanetto, 2012; Fragos et al., 2013; Ross et al., 2017; Mirocha et al.,
2021; Schaeffer et al., 2023).

However, it is not clear when is it safe to ignore galaxy to galaxy scatter
(i.e. stochasticity). Stochasticity can be important even when estimating
globally-averaged quantities such as the mean EoR history. Assuming
population-averaged quantities (e.g. ionizing escape fraction, stellar to halo
mass relation, etc.) can give biased results for correlated distributions
(e.g. the average of a product is not the same as the product of the averages;
c.f. Appendix B for simple examples). Moreover, some measurements (e.g. 21-cm
interferometry, Ly$\alpha$ forest, kinetic Synaev-Zel’dovich signal, etc.) are
sensitive to the spatial fluctuations in the galaxy emissivity, on some range
of spatial scales. As that scale is reduced, there are fewer galaxies over
which to average, and stochasticity becomes more important (e.g. Davies &
Furlanetto, 2016). The importance of stochasticity also increases at high
redshifts, where sources are rarer and more biased. It has been evoked to
explain controversial claims at $z>10$, such as a rapid redshift evolution of
the global 21cm signal during the CD (e.g. Kaurov et al., 2018), and an
overabundance of massive galaxy candidates from JWST photometry at $z>10$
(e.g. Mirocha & Furlanetto, 2023; Mason et al., 2023; Shen et al., 2023).

Here we construct a model of galaxy emissivity in the bands that are relevant
for interpreting current and upcoming observations of the EoR and the CD: (i)
ionizing UV (which drives the EoR and determines the residual HI fraction in
the ionized IGM); (ii) soft X-ray (which heats and partially ionizes the IGM
during the CD); (iii) Lyman Werner (which determines when H2 cooling stops
being efficient in the first galaxies). We compute the distribution of these
multi-frequency emissivities as a function of scale and redshift. Our model
samples the largest expected sources of stochasticity, including: the
abundance of dark-matter halos, stellar-to-halo mass relation, galaxy main
sequence, fundamental metallicity relation, luminosity and escape fraction
scalings. We quantify the relative importance of each term to the total
emissivity in each of the considered bands. We also quantify the importance of
these stochastic terms for simple estimates of the EoR history, as well as the
high redshift UV luminosity functions (UV LFs). Our results can be used to
improve estimates of cosmic radiation fields and guide models of the EoR/CD by
highlighting the most important sources of scatter.

The structure of this paper is as follows. In Section 2 we introduce our model
for calculating galaxy emissivities. In Section 3 we present the resulting UV,
X-ray, LW emissivity distributions, quantifying the relative importance of
each source of stochasticity. In Section 4 we show two analytic estimates of
the EoR history, quantifying the relative impact of ignoring galaxy to galaxy
scatter. In Section 5 we show the UV LFs implied by our fiducial model,
comparing them to observational estimates from photometric candidates.
Finally, we conclude in Section 6. All quantities are presented in comoving
units unless stated otherwise. Throughout this work, we assume standard
$\Lambda$CDM cosmological parameters
($\Omega_{\textrm{m}},\Omega_{\textrm{b}},\Omega_{\Lambda},h,\sigma_{8},n_{\textrm{s}}=0.310,0.049,0.689,0.677,0.81,0.963$),
consistent with the latest estimates from Planck Collaboration et al. (2020).

##  2 Computing emissivities at high redshifts

If galaxy properties could be written as deterministic functions of the mass
of their host halos and/or redshift, we could write the emissivity (e.g. erg
s-1 cMpc-3) in some spectral band, $i$, at a redshift $z$ as:

$\varepsilon_{i}(z)=\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h},z)}{\textrm{d}M_{h}}L_{i}(M_{h},z)\
f_{\rm esc,i}(M_{h},z)\leavevmode\nobreak\ .$ (1)

Here $\frac{\textrm{d}n}{\textrm{d}M_{\rm h}}$ is the number density of halos
per unit mass (i.e. the halo mass function; HMF), $L_{i}$ is the intrinsic
luminosity of a galaxy hosted in a halo of mass $M_{h}$ at redshift $z$, and
$f_{{\rm esc},i}$ is the fraction of photons that escape the galaxy to make it
into the IGM.

However, we know that the above relations are not deterministic functions of
halo mass and redshift. Complex physics of galaxy evolution and radiative
transfer induces a spread around relations linking different galaxy
properties. Nevertheless, there are empirically well-established relations
that characterize some of the main correlations of galaxy properties.
Therefore, a more general form for the mean of the emissivity would
marginalize over these relations. Specifically, we can write the emissivity in
a spectral band, $i$, at redshift $z$, averaged over comoving volumes
$(4/3)\pi R_{\rm nl}^{3}$, as:111Throughout we use the subscript ”nl” to
indicate non-linear (Eularian) quantities and the subscript ”0” to indicate
Lagrangian quantities linearly evolved to $z=0$ (following convention). We
recall that all length scales are in comoving units, unless otherwise
specified.

$\displaystyle\bar{\varepsilon}_{i}(R_{\rm nl},z)=$
$\displaystyle\int\textrm{d}M_{h}\int\textrm{d}\delta_{0}\frac{\textrm{d}n(M_{h},z\leavevmode\nobreak\
|\leavevmode\nobreak\
R_{0},\delta_{0})}{\textrm{d}M_{h}}p_{z}(\delta_{0}\leavevmode\nobreak\
|\leavevmode\nobreak\ R_{\rm nl})$ (HMF) $\displaystyle\times$
$\displaystyle\int\textrm{d}M_{\ast}\ p(M_{\ast}\leavevmode\nobreak\
|\leavevmode\nobreak\ M_{\rm h})$ (SHMR) $\displaystyle\times$
$\displaystyle\int\textrm{dSFR}\ p_{z}(\textrm{SFR}\leavevmode\nobreak\
|\leavevmode\nobreak\ M_{\ast})$ (SFMS) $\displaystyle\times$
$\displaystyle\int\textrm{d}Z\ p_{z}(Z\leavevmode\nobreak\
|\leavevmode\nobreak\ \textrm{SFR},M_{\ast})\ $ (FMR) $\displaystyle\times$
$\displaystyle\int\ \textrm{d}L_{i}\ L_{i}\ p(L_{i}\leavevmode\nobreak\
|\leavevmode\nobreak\ \textrm{SFR},Z)$ (L) $\displaystyle\times$
$\displaystyle\int\textrm{d}f_{\rm esc,i}\ f_{\rm esc,i}\ p(f_{\rm esc,i})$
(EF)

Here $\delta_{0}(R_{0})$ is the linear matter overdensity of a spherical
volume of Lagrangian radius $R_{0}$ corresponding to the final Eulerian radius
$R_{\rm nl}$, $p(M_{\ast}\leavevmode\nobreak\ |\leavevmode\nobreak\ M_{\rm
h})$ is the conditional probability of stellar mass $M_{\ast}$ for a given
$M_{\rm h}$, $p_{z}({\rm SFR}\leavevmode\nobreak\ |\leavevmode\nobreak\
M_{\ast})$ is the conditional probability of a star-formation rate (SFR) for a
given stellar mass222We use ’$z$’ subscripts to indicate probability
distributions that are also functions of redshift (see below for more
details)., $p_{z}({\rm Z}\leavevmode\nobreak\ |\leavevmode\nobreak\ {\rm
SFR},M_{\ast})$ is the conditional probability of a stellar metallicity $Z$
for a given SFR and $M_{\ast}$, $p(L_{i}\leavevmode\nobreak\
|\leavevmode\nobreak\ \textrm{SFR},Z)$ is the conditional probability of a
luminosity $L_{i}$ in a given wavelength band $i$ for a given SFR and $Z$, and
$p(f_{\rm esc,i})$ is the probability of an escape fraction $f_{\rm esc,i}$ in
band $i$. Loosely speaking, the running averages of the conditional
probabilities in the first four rows are commonly referred to as the halo mass
function (HMF)333Strictly speaking, the average HMF is not the mean of the
conditional halo mass function (CHMF), but its limit as
$R_{0}\rightarrow\infty$ and $\delta_{0}\rightarrow 0$., stellar-to-halo mass
relation (SHMR), star forming main sequence (SFMS) of galaxies, and
fundamental metallicity relation (FMR); we label the corresponding rows in the
equation above with these acronyms, and go through each probability
distribution in more detail below. In principle, the PDFs above could be
conditioned on additional galaxy properties, which could further increase the
importance of stochasticity. Note that higher order moments of the emissivity,
such as its variance
$\langle(\varepsilon_{i}-\bar{\varepsilon}_{i})^{2}\rangle$, can be similarly
expressed in terms of the above conditional probability distributions.

For general distributions, Eqs. 1 and LABEL:eq:full_emissivity do not give the
same mean. This means that even interpreting average quantities like the EoR
history could be biased if not accounting for stochasticity. More
fundamentally, the various sources of scatter in Eq. LABEL:eq:full_emissivity
result in spatial fluctuations in the emissivity which can be important for
many EoR/CD observations. To date, the impact of this scatter on EoR/CD
observables has only been explored in a limited fashion. For example, Hassan
et al. (2021) found that scatter in the intrinsic production rate of ionizing
photons predicted by the Simba simulation (Davé et al., 2019) has only a
modest impact on the EoR morphology. On the other hand, Reis et al. (2022)
used a toy model to characterize the effective scatter in star formation
efficiency, finding in some cases a large impact on the EoR and CD morphology;
though see Shekhar Murmu et al. (2023) for the opposite conclusions using a
different astrophysical model. Indeed, a sizable scatter in the ionizing
emissivity is needed to explain the latest Lyman alpha forest data at $z=5$ –
6.3 (Qin et al. 2021; Gaikwad et al. 2023; Davies et al. 2023a, Qin et al. in
prep.).

We could analytically derive the emissivity distribution,
$p(\varepsilon_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\ R_{\rm nl},z)$,
if the conditional distributions in Eqs. HMF-EF followed simple Gaussian
forms. In the more general case, we can solve for
$p(\varepsilon_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\ R_{\rm nl},z)$ by
numerically sampling the above relations. Specifically, to compute a single
realization (denoted below with a ”tilde”) of the emissivity, we perform the
following Monte Carlo (MC) procedure:

Algorithm 1 Algorithm 1: Computing a single realization of the emissivity in
spectral band i of a region of radius $R_{\rm nl}$ at redshift $z$:
$\tilde{\varepsilon}_{i}(R_{\rm nl},z)$

1: sample the linear matter overdensity $\tilde{\delta}_{0}\sim
p_{z}(\delta_{0}\leavevmode\nobreak\ |\leavevmode\nobreak\ R_{\rm nl})$

2: obtain a realization of the halo field by sampling the CHMF:
$\widetilde{\\{M_{h}^{\rm j}\\}}\sim\textrm{d}n(M_{h},z\leavevmode\nobreak\
|\leavevmode\nobreak\ \tilde{\delta}_{0},R_{0})/\textrm{d}M_{h}$

3: for all halos $j$ with mass $\tilde{M}_{h}^{{\rm j}}$ do

4: sample probability that the halo hosts an actively star-forming galaxy,
$p(t_{\rm duty}\leavevmode\nobreak\ |\leavevmode\nobreak\ \tilde{M}_{h}^{\rm
j})$

5: if halo does not host a star-forming galaxy then

6: CONTINUE

7: end if

8: sample stellar mass $\tilde{M}_{\ast}^{j}$ $\sim
p(M_{\ast}\leavevmode\nobreak\ |\leavevmode\nobreak\ \tilde{M}_{\rm h}^{j})$

9: sample star formation rate $\tilde{{\textrm{SFR}}}^{j}$ $\sim
p_{z}(\textrm{SFR}\leavevmode\nobreak\ |\leavevmode\nobreak\
\tilde{M}_{\ast}^{j})$

10: sample metallicity $\tilde{Z}^{j}$ $\sim p_{z}(Z\leavevmode\nobreak\
|\leavevmode\nobreak\ \tilde{\textrm{SFR}}^{j},\tilde{M}_{\ast}^{j})$

11: sample intrinsic luminosity $\tilde{L}_{i}^{j}$ $\sim
p(L_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\
\tilde{\textrm{SFR}}^{j},\tilde{Z}^{j})$

12: sample escape fraction $\tilde{f}_{\rm esc,i}^{j}$ $\sim p(f_{\rm esc,i})$

13: end for

14: $\tilde{\varepsilon_{i}}=\sum_{{\rm halo}=j}\tilde{L}_{i}^{j}$
$\tilde{f}_{\rm esc,i}^{j}$

We describe each step of the above MC procedure in turn below.

###  2.1 Halo Mass Function (HMF)

In this work we wish to compute the distribution of galaxy emissivities of
regions of a given Eulerian scale at a given redshift,
$p(\varepsilon_{i}|R_{\rm nl},z)$. Our model is anchored by the fact that
galaxies are hosted by dark matter halos, whose relative abundances are
described by conditional halo mass functions.

Here we use the hybrid CHMF proposed by Barkana & Loeb (2004), in which the
analytically-tractable Press-Schecher CHMF (Press & Schechter, 1974) is
normalized to have the same mean as the (non conditional) Sheth-Tormen HMF
(ST, Sheth & Tormen, 1999):

$\displaystyle\begin{split}&\frac{\textrm{d}n}{\textrm{d}M_{\rm h}}(M_{\rm
h},z\leavevmode\nobreak\ |\leavevmode\nobreak\
\delta_{0},R_{0})=\frac{\overline{f}_{\rm ST}}{\overline{f}_{\rm
PS}}\sqrt{\frac{2}{\pi}}\frac{\overline{\rho}}{M_{\rm h}}\times\\\
&\times\frac{\delta_{\rm
c}-\delta_{0}}{\sigma^{2}(M_{h})-\sigma^{2}(R_{0})}\left|\frac{\textrm{d}\sigma}{\textrm{d}M_{\rm
h}}\right|\exp{\left(-\frac{(\delta_{\rm
c}-\delta_{0})^{2}}{\sigma^{2}(M_{h})-\sigma^{2}(R_{0})}\right)}\end{split}$
(3)

In the above, $\overline{f}_{\rm ST}$ and $\overline{f}_{\rm PS}$ correspond
to the mean Sheth-Tormen and Press-Schechter collapsed fractions above the
atomic cooling threshold of $T_{\rm vir}\geq 10^{4}$ K, respectively,
$\delta_{\rm c}(z)$ is the critical linear density from the spherical collapse
model, and $\sigma^{2}(M)$ is the mass variance of the Lagrangian (linear)
density field on scales $M=(4/3)\pi R_{0}^{3}$.

In order to sample from Eq. 3, we need to connect Lagrangian and Eulerian
quantities (see also, e.g. Trapp & Furlanetto 2020). In Lagrangian space,
$p(\delta_{0})$ follows a zero-mean Gaussian distribution whose width is
determined by $\sigma^{2}(R_{0})$. We transform this distribution to Eulerian
space using the spherical collapse model (e.g. Mo & White 1996):

$p(\delta_{\rm nl}\leavevmode\nobreak\ |\leavevmode\nobreak\ R_{\rm
nl},z)\textrm{d}\delta_{\rm nl}=\frac{1}{1+\delta_{0}}f_{\rm
R}(\sigma^{2}\leavevmode\nobreak\ |\leavevmode\nobreak\
R_{0},z)\textrm{d}\sigma^{2}\leavevmode\nobreak\ ,$ (4)
$\delta_{0}=-1.35(1+\delta_{\rm nl})^{-2/3}+0.78785(1+\delta_{\rm
nl})^{-0.58661}\\\ -1.12431(1+\delta_{\rm
nl})^{-1/2}+1.68647\leavevmode\nobreak\ .$ (5)

Here $f_{\rm R}$ is the first-crossing distribution from Sheth (1998), and
$\delta_{\rm nl}=\rho/\bar{\rho}-1$ is the Eulearian (non-linear) overdensity.

Figure 1: Example halo mass functions used in this work at three different
redshifts. The dashed black (solid red) curves correspond to the mean (single
realization) of the halo field conditioned on a region of scale $R_{\rm nl}=5$
Mpc having a density equal to the cosmic average. The sample variance scatter
in the red curve is seen to increase towards large masses and high redshifts,
as the target mean values become smaller. In the top panel we also show the
non-conditional HMF (i.e. the limit as $R_{0}\rightarrow\infty$ at
$\delta_{0}=0$).

With the above relations, we generate a Lagrangian overdensity sample,
$\tilde{\delta}_{0}\sim p(\delta_{0}\leavevmode\nobreak\ |\leavevmode\nobreak\
R_{\rm nl},z)$ (step 1 of the MC procedure in the previous subsection). We
then compute a corresponding realization of the halo field according to the
following procedure. We first sample the total number of halos with masses
above some arbitrary minimum value, obtaining $\tilde{N}(>M_{\rm
min}\leavevmode\nobreak\ |\leavevmode\nobreak\ \tilde{\delta}_{0},R_{0})$, by
assuming a Poisson distribution whose mean is given by the integral of Eq. 3
from $M_{\rm min}$ to infinity. We then assign each halo a mass by sampling
the normalized cumulative mass function, using rejection to ensure the total
mass is within $\pm$10% of the target mean.444Although approximate, our
approach has a couple of notable advantages over other simple MC
implementations of stochasticity in which halo numbers are sampled from
independent Poisson distributions in fixed mass bins (e.g. Reis et al., 2022).
Firstly, by sampling a continuous CDF, we avoid binning halo masses and
forcing them to have discrete values. Moreover, the mean total number of halos
is much larger than the mean number in any given mass bin, validating the
assumption of a Poisson distribution. Furthermore, having a mass error
threshold ensures approximate mass conservation in each realization. The
alternative of not correlating halo samples in neighboring mass bins and not
ensuring mass conservation can significantly overestimate the importance of
stochasticity when halos become rare, which can explain why our results are
different from those in Reis et al. (2022).

We show some example CHMFs in Fig.1. The dashed black curves correspond to the
target mean CHMFs in regions of Eularian scale $R_{\rm nl}$, at mean density,
at redshifts $z=5,10$ and $15$ (top to bottom panels). The red solid curves
show a single realization, computed according to the above procedure. The
impact of stochasticity in the red curves is very evident as the mean number
decreases, i.e. towards high redshifts and high masses. In the top panel we
also show the mean (non conditional) HMF (i.e. the limit as
$R_{0}\rightarrow\infty$ at $\delta_{0}=0$).

Fig.1 also highlights that there are effectively two sources of scatter when
determining the halo abundances in a given volume: (i) the scatter in the mean
value of the CHMF, driven by its dependence on the underlying matter
overdensity (i.e. the difference between the blue and black curves); and (ii)
the scatter due to discrete sampling around the target mean CHMF (i.e. the
difference between the black and red curves). The former determines
cosmological signals like 21cm since it is correlated to the underlying matter
field. The latter on the other hand is effectively a sample noise term. Both
sources of scatter are naturally accounted for in $N$-body simulations,
although periodic boundary conditions mean that (i) is underestimated due to
limited box sizes (e.g. Barkana & Loeb 2004). On the other hand, analytic and
semi-numerical models of inhomogeneous radiation fields account for (i), but
often assume (ii) is negligible in order to reduce computation costs. Below we
confirm the validity of this approximation.

###  2.2 Stellar-to-halo mass relation (SHMR)

Figure 2: Uppermost panel: Our (redshift independent) stellar to halo mass
relation (solid curve) and 2$\sigma$ scatter (shaded region). Lower panels:
Galaxy star-forming main-sequence (solid curves) and 2$\sigma$ scatter (shaded
regions) at $z=5$, 10, 15 (top to bottom).) Coloured symbols represent
galaxies from cosmological simulations, circles for FirstLight (Ceverino et
al., 2018), stars for (Lovell et al., 2021), crosses for ASTRID (Bird et al.,
2022; Davies et al., 2023b) and triangles for SERRA (Pallottini et al., 2022).
For ASTRID we randomly select galaxies in fixed mass bins, to avoid over
crowding the plot, and for SERRA we use their $z=6$ and $z=12$ snapshots for
$z=5$ and 10, respectively.

Both observations and theory established a strong relation between the stellar
and halo masses of galaxies (e.g. Harikane et al., 2016; Ceverino et al.,
2018; Stefanon et al., 2021; Lovell et al., 2021; Kannan et al., 2022;
Pallottini et al., 2022; Di Cesare et al., 2023). Here we assume a log-normal
conditional probability of a galaxy having a stellar mass, $M_{\ast}$, given a
host halo mass, $M_{h}$: $p(\log M_{\ast}|\log M_{h})=\mathcal{N}(\log
M_{\ast}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mu_{M_{\ast}}(M_{h}),\sigma_{M_{\ast}})$. We assume a mass-independent
$\sigma_{M_{\ast}}$ of 0.25 dex (e.g. Ceverino et al., 2018; Lovell et al.,
2021; Pallottini et al., 2022) and a mean given by the following double power
law SHMR:

$\displaystyle\begin{split}\mu_{M_{\ast}}(M_{h})&=\\\ -1.412+&\log{M_{\rm
h}}-\log{\left[\left(\frac{M_{\rm h}}{2.6\times
10^{11}M_{\odot}}\right)^{-0.5}+\left(\frac{M_{\rm h}}{2.6\times
10^{11}M_{\odot}}\right)^{0.6}\right]}\end{split}$ (6)

A standard physical interpretation of the double power law form is that the
low-mass scaling is determined by stellar feedback while the high-mass scaling
is determined by AGN feedback (e.g. Wechsler & Tinker 2018; Behroozi et al.
2019 and references therein). In Eq. (6) the normalization and the low-mass
power-law index correspond to the maximum a posteriori (MAP) values inferred
from a combination of CMB, QSO and high-redshift UV LF observations in Nikolić
et al. (2023), while the high-mass power-law index is taken from the bright-
end UV LF empirical fits in Mirocha et al. (2017). Our results mostly depend
on the former, as the steepness of the HMF at high redshifts means that early
radiation fields are dominated by the faint (low mass) galaxies (e.g. Bouwens
et al., 2015, 2023; Gillet et al., 2020, see also below).

Gas accreting from the IGM onto halos is gravitationally heated, and can also
be photo-heated by the ionizing UV background (UVB) during the EoR. In order
to condense onto the galaxy and form stars, this gas needs to cool. Cooling
can be inefficient in halos with small virial temperatures, with an
exponentially decreasing fraction of halos capable of sustaining star
formation (e.g. Sobacchi & Mesinger, 2013; Xu et al., 2016). Here we account
for this effect by assuming only a fraction $f_{\textrm{duty}}(M_{\rm
h})=\exp[-M_{\textrm{turn}}/M_{\rm h}]$ of halos host star-forming galaxies,
taking $M_{\rm turn}=5\times 10^{8}M_{\odot}$ based on the inference result in
Nikolić et al. (2023). Specifically, for each halo we sample a random variable
uniformly between 0 and 1, and only populate the halo with a star forming
galaxy if the value of the random variable is less than $f_{\rm duty}$.

In the top panel of Fig. 2 we show our mean SHMR and 2$\sigma$ scatter (solid
black line and gray shaded region, respectively). The mean SHMR is a power law
over most of the mass range shown. At the high (low) mass end we see a
flattening due to our parametrization of AGN feedback (inefficient accretion),
as discussed above. For comparison, we also show galaxies from several
hydrodynamic simulations: FirstLight (Ceverino et al., 2018), ASTRID (Bird et
al., 2022; Davies et al., 2023b), and SERRA (Pallottini et al., 2022). The
simulated galaxies are colored according to their redshift, with red for
$z=5$, green for $z=10$ and blue for $z=15$.

We see significant differences in the (mean) SHMR between different
simulations. The cosmological zoom-in SERRA simulations imply a mean SHMR that
is roughly two orders of magnitude higher at the smallest halo masses compared
with the ASTRID simulations. FLARES and FirstLight are somewhere in between
these two extremes, as is our fiducial model. We remind our reader that our
mean relation was inferred from data, as discussed in Nikolić et al. (2023),
and not based on these simulations.

Conversely, our choice of 0.25 dex scatter around the mean relation is roughly
motivated by the galaxy-to-galaxy scatter found in any given hydrodynamic
simulation. Interestingly, despite the fact that different simulations predict
different means, the scatter around the mean is roughly comparable.
Furthermore, we see that the simulations do not show strong evidence of a
redshift evolution of the SHMR, justifying our fiducial model (see also, e.g.
Mutch et al., 2016; Harikane et al., 2016; Tacchella et al., 2016; Ma et al.,
2018; Yung et al., 2019).

###  2.3 Galaxy star formation main sequence (SFMS)

The star-formation rates of galaxies, SFRs, are known to be strongly
correlated with their stellar mass content. The mean of this SFR – $M_{\ast}$
relation is loosely referred to as the galaxy star formation main sequence
(SFMS); galaxies with SFRs significantly above (below) the SFMS are referred
to as bursty (quenched). The SFMS is well established observationally at low
redshifts and (comparably) large masses (e.g. Brinchmann et al., 2004; Santini
et al., 2017; Curtis-Lake et al., 2021; Popesso et al., 2022). The observed
mean relation at small masses follows a power law, whose index is fairly
constant but whose normalization decreases with redshift. This decrease with
redshift is naturally reproduced if one assumes that the star formation time-
scale is related to the free-fall time at the mean viral density of host
halos, $t_{\rm ff}$, which during matter domination scales as the Hubble time:
$t_{\rm ff}\propto H^{-1}(z)$ (e.g. see Park et al., 2019, and references
therein).

At a given $z$, we again assume a log-normal conditional probability
$p_{z}(\log{\rm SFR}\leavevmode\nobreak\ |\leavevmode\nobreak\ \log
M_{\ast})=\mathcal{N}[\log{\rm SFR}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mu_{\rm SFR}(M_{\ast},z),\sigma_{{\rm SFR}}(M_{\ast})]$. For the mean SFMS we
use the model of Park et al. (2019), with the normalization set by the MAP
values in Nikolić et al. (2023):

$\mu_{\rm SFR}(M_{\ast},z)=\log{M_{\ast}}-\log{[0.43H^{-1}(z)]}$ (7)

We assume a mass-dependent scatter that increases towards smaller masses, as
these galaxies are expected to be more bursty555Throughout this work we use
”burstiness” to indicate a wide scatter around the mean SFMS. We do not
investigate what such distributions imply for the star formation histories of
individual galaxies.:

$\sigma_{\rm SFR}(M_{\ast})=\cases{-}0.12\log{M_{\ast}}+1.35&\text{if $\log
M_{\ast}<10$}\\\ 0.19\text{otherwise}{}$ (8)

The normalization and scaling of the scatter was fit to the hydrodynamic
simulations of Ceverino et al. (2018).

We plot our assumed SFMS and 2$\sigma$ scatter in the bottom three panels of
Figure 2. Panels correspond to $z=5$, 10, 15 (top to bottom), with different
symbols indicating values taken from hydrodynamic simulations: FirstLight
(Ceverino et al., 2018), FLARES (Lovell et al., 2021), ASTRID (Bird et al.,
2022), and SERRA (Pallottini et al., 2022). The figure illustrates that our
fiducial model is in general agreement with results from these hydrodynamic
simulations.666Detailed comparisons to other works would require standardizing
definitions. For example, here we define the SFR as an instantaneous quantity,
while elsewhere it could be averaged over $\sim 100$ Myr to allow for a more
direct comparison to photometric observations. Here we are just interested in
confirming that our fiducial choices are reasonable.

###  2.4 Fundamental Metallicity Relation (FMR)

Figure 3: The stellar metalicities, stellar masses and star formation rates of
galaxies from a single realization of a $R_{\rm nl}=5$ cMpc volume at mean
density at $z=6$. Each point corresponds to a single galaxy, with the color
denoting its metallicity: $\tilde{Z}^{j}$ $\sim p_{z}(Z\leavevmode\nobreak\
|\leavevmode\nobreak\ \tilde{\textrm{SFR}}^{j},\tilde{M}_{\ast}^{j})$ (see
text for details). Stars denote the observationally-estimated means in three
bins for $z>6$ galaxies in Curti et al. (2023) (converted from gas
metalicities, see text for details), as well as the metallicity estimate of
GN-z11 from Bunker et al. (2023).

The galaxy emissivity also depends on the metallicity of the stellar
population. Here we relate the metallicity to the SFR and stellar mass of a
galaxy, taking advantage of the well-studied fundamental metallicity relation
(FMR; Mannucci et al. 2010; Curti et al. 2020). Specifically, we assume a log-
normal conditional probability of a galaxy having a stellar metallicity $Z$,
given its SFR and stellar mass, $p_{z}(\log Z\leavevmode\nobreak\
|\leavevmode\nobreak\ \log{\rm SFR},\log M_{\ast})=\mathcal{N}[\log
Z\leavevmode\nobreak\ |\leavevmode\nobreak\ \mu_{Z}(M_{\ast},{\rm
SFR},z),\sigma_{Z}]$. We assume a constant scatter of $\sigma_{Z}=0.1$ dex,
and a mean given by the following (c.f. Curti et al. 2020):

$\mu_{Z}(M_{\ast},{\rm
SFR},z)=0.296\left(1+\left(\frac{M_{\ast}}{M_{0}}\right)^{-2.1}\right)^{-0.148}10^{\Delta_{z}}\
Z_{\odot}\leavevmode\nobreak\ ,$ (9)

where $M_{0}(\rm SFR)\equiv 10^{10.11}\times({\rm
SFR/M_{\odot}{yr^{-1}})}^{0.56}M_{\odot}$, and $\Delta_{z}=-0.056z+0.064$
accounting for putative redshift evolution (Curti et al., 2023). In the above,
we converted from gas phase to stellar metalicities using
$Z/Z_{\odot}=10^{(12+\log{(\textrm{O}/\textrm{H})}-8.69)}$ with solar
metallicity $Z_{\odot}=0.02$ (Asplund et al., 2004), and adjusting for gas
phase metallicities being higher by a factor of $\approx 2.63$ on average
(Strom et al., 2018).

In Fig. 3 we show galaxies from a single realization of a comoving volume with
radius $R_{\rm nl}=5$ cMpc, at mean density at redshift $6$. Each point
denotes a single galaxy with the color corresponding to its typical stellar
metallicity. Note that the apparent scatter in the metallicity at a fixed
$M_{\ast}$ is considerably larger than the 0.1 dex scatter we set around the
mean FMR at a given $M_{\ast}$ and SFR. This is because the intrinsic scatter
in the combination of SHMR and SFMS (i.e. the width of $p_{z}({\rm
SFR},M_{\ast})$ from the previous sections) dominates over our choice of
scatter in the metallicity given these properties (i.e. the width of
$p_{z}(Z\leavevmode\nobreak\ |\leavevmode\nobreak\ {\rm SFR},M_{\ast})$; see
also e.g. Garcia et al. 2024). In the figure we also show the binned values
for the metallicity of $z>6$ galaxies from Curti et al. (2023) as well as the
metallicity estimate of GN-z11 from Bunker et al. (2023). Although these
observations span a range of redshifts, they are generally consistent with our
samples.

###  2.5 Luminosity scalings

The intrinsic luminosity of a galaxy depends primarily on the SFR and its
history, as well as the metallicity of the stellar population (e.g. Brammer et
al., 2008; Allende Prieto et al., 2018; Stanway & Eldridge, 2018; Lehmer et
al., 2021; Fragos et al., 2023). Here we describe how we compute the intrinsic
luminosities for each of the wavelength bands of interest: X-ray, ionizing and
Lyman Werner.

####  2.5.1 Soft-band X-ray luminosity

Figure 4: Upper panel: scaling of soft-band, high-mass X-ray binary luminosity
with SFR. Red and blue lines with the corresponding shaded regions represent
the mean and $2\sigma$ range for metalicities $Z=-3.6,-2.0$, respectively.
Dash-dotted and dashed lines correspond to galaxy-wide high mass X-ray binary
LF estimates discussed in Lehmer et al. (2021) (mean and enclosing $2\sigma$,
respectively). Green stars correspond to values from local star forming
galaxies, discussed in Brorby et al. (2016). Lower panels: SFR scaling of the
integrated Lyman-Werner (11.2 – 13.6 eV) luminosity (left panel) and specific
ionizing luminosity at the Lyman-limit (right panel). Shaded regions represent
the $2\sigma$ scatter around the mean relation. Both panels assume that
metallicity follows the mean FMR.

Soft-band777Here we define the soft band to be $0.5-2\ \textrm{keV}$. Roughly
speaking, photons with higher energies do not interact with the high-$z$ IGM
(e.g. Oh, 2001; Xu et al., 2014; Madau & Fragos, 2017), while photons with
lower energies get absorbed inside the host galaxies (e.g. Das et al., 2017,
see also Section 2.6). X-rays emerging from the first galaxies are responsible
for heating and partially ionizing the IGM during the cosmic dawn (e.g.
McQuinn, 2016), which can have a dramatic imprint in the cosmic 21cm signal
(e.g. Mesinger et al., 2013; Abdurashidova et al., 2022). It is likely that
the X-ray emissivity of $z>6$ galaxies is dominated by high mass X-ray
binaries (HMXBs; e.g. Furlanetto, 2006; Fragos et al., 2013; Pacucci et al.,
2014; Eide et al., 2018). HMXBs are massive stars accreting onto a compact
companion. The total X-ray output of a galaxy from HMXBs should therefore
scale with the SFR of the galaxy (due to the rapid stellar evolution
timescales of massive stars) and its metallicity (which determines the
efficiency of radiative-driven winds and the resulting mass loss of the
massive companion). Indeed we observe a strong dependence of the X-ray
luminosity on the galaxy’s SFR and metallicity in local galaxies and in stacks
out to $z\sim 2.5$ (e.g. Brorby et al., 2016; Lehmer et al., 2016; Fornasini
et al., 2019; Lehmer et al., 2021).

Here we assume a log-normal conditional probability of a galaxy having an
intrinsic soft-band X-ray luminosity, $L_{\rm X}$ (in units of erg s-1), given
a SFR and metallicity: $p(\log L_{\rm X}\leavevmode\nobreak\
|\leavevmode\nobreak\ \log{\rm SFR},\log Z)=\mathcal{N}(\log L_{\rm
X}\leavevmode\nobreak\ |\leavevmode\nobreak\ \mu_{\rm X}({\rm
SFR},Z),\sigma_{\rm X})$. We assume a constant $\sigma_{\rm X}$ of 0.5 dex and
a mean given by:

$\displaystyle\begin{split}\mu_{\rm X}({\rm
SFR},Z)=&\quad(-0.11\times\log{(Z/Z_{\odot})}+1.30)\ \log{\textrm{SFR}}\\\
&+(-0.31\times\log{(Z/Z_{\odot})}+38.44)\leavevmode\nobreak\ .\end{split}$
(10)

These fiducial choices are based on the galaxy-averaged HMXB luminosity
functions in Lehmer et al. (2021), and are roughly consistent with empirical
fits to local galaxies (e.g. Brorby et al. 2016). Note that we changed the
dependence from the gas-phase metallicity to the stellar one using the
prescription from Section 2.4. We have also converted the hard band X-ray
luminosity in Lehmer et al. (2021) ($0.5$-$8$ keV) to the soft band one by
multiplying the luminosities by a factor of $0.3$, consistent with
observational estimates (e.g. Basu-Zych et al., 2013) and corresponding to an
intrinsic SED with a power-law index of $\Gamma=2.0$ (e.g. Mineo et al.,
2012).

We show this dependence of the X-ray luminosity with SFR, for two different
metallicity values, in the top panel of Figure 4. For comparison, we include
some Chandra observations of local actively star-forming galaxies from Lehmer
et al. (2021) and Lyman-Break analogues from Brorby et al. (2016). We see that
our fiducial model is consistent with current data; however, it is highly
uncertain how these relations scale to the first galaxies whose metallicity
ranges are not sampled by current observations (e.g. Magg et al., 2022; Kaur
et al., 2022).

####  2.5.2 Ionizing and Lyman Werner luminosities

We use the Binary Population and Spectral Synthesis (BPASS) code to compute
intrinsic ionizing and Lyman Wener UV luminosities (Stanway & Eldridge, 2018;
Byrne et al., 2022). BPASS provides a deterministic prediction for the UV
luminosity, $L_{\rm UV}$, as a function of SFR, Z, and SFR history. For the
latter we assume that our sampled SFR is exponentially declining towards
higher redshifts, as implied by Equation 7. Therefore $L_{\rm UV}$ is sampled
assuming a log-normal conditional probability $p(\log L_{\rm
UV}\leavevmode\nobreak\ |\leavevmode\nobreak\ \log{\rm SFR},\log
Z)=\mathcal{N}(\log L_{\rm UV}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mu_{\rm UV,BPASS}({\rm SFR},Z),\sigma_{\rm L})$ where $\mu_{\rm UV,BPASS}$ is
the predicted luminosity from BPASS. We add an additional scatter of
$\sigma_{\rm L}=0.1$ dex around the mean to compensate for unaccounted sources
of stochasticity, e.g. the mean IMF, alpha-element distribution, etc. (Byrne &
Stanway, 2023). However, this level of scatter is negligible compared to the
scatter of the bulk galaxy properties like SFR and stellar mass. We show the
scaling relation of $L_{\rm UV}$ with $\rm SFR$ in the bottom panels of Fig. 4
for the 11.2–13.6 eV Lyman-Werner band (left panel; in units of erg s-1) and
Lyman-limit (right panel; in units of erg s${}^{-1}\AA^{-1}$ evaluated at 13.6
eV).

###  2.6 Escape fractions

Our final step in computing the emissivity is determining what fraction of the
produced photons manage to escape the host galaxy into the IGM. This is
referred to as the escape fraction. We use different prescriptions for the
escape fraction in our three bands of interest. We describe each in turn
below.

Both hydrodynamic simulations (e.g. Cen & Kimm, 2015; Xu et al., 2016; Barrow
et al., 2020; Yeh et al., 2022; Kostyuk et al., 2023) and direct observations
of low redshift galaxies (e.g. Izotov et al., 2016; Grazian et al., 2017;
Steidel et al., 2018; Pahl et al., 2023) show sizable stochasticity in the
ionizing escape fraction, though there is no consensus on what is an
appropriate distribution. Here we take two scenarios. Our fiducial model
assumes a log-normal distribution for the ionizing escape fraction with a
width of 0.3 dex (c.f. Mascia et al. 2023), while we also show a bimodal
distribution in which galaxies have an ionizing escape fraction of either 0 or
1 (resulting in maximum scatter). In both cases we take the inference result
of Nikolić et al. (2023): $\bar{f}_{\rm esc}=0.053$ for the mean. We do not
assume a dependence of the escape fraction on galaxy properties, as such
relations are not yet well established at high redshifts. However, our
framework can easily be extended to include putative dependencies on galaxy
properties (e.g. Mascia et al., 2023), as well as accommodating different
functional distributions (e.g. Kreilgaard et al. 2024).

For the X-ray escape fraction, we adopt the results of Das et al. (2017),
where they computed the X-ray opacities of simulated high-$z$ galaxies,
finding that most photons with energies above 0.5 keV manage to escape.
Following that work, we assume an escape fraction of unity above 0.5 keV and
zero below that value. Similarly, we assume values of unity for the Lyman
Werner escape fraction, given the typical low opacities of such photons
through the host ISM (e.g. Haiman et al., 2000; Wolcott-Green et al., 2011).

##  3 Results: emissivities

Here we present our distributions for the emissivity in each of the three
bands in turn. We show the full distributions as a function of redshift,
before quantifying the relative importance of each source of scatter. For the
latter, we compute the mean and standard deviations of the emissivity PDF when
one source of scatter is removed (i.e. using only the corresponding mean
relation with zero scatter), normalized to the values of the full distribution
containing all sources of scatter: $\mu_{\epsilon}$/$\mu_{\epsilon}^{\rm
full}$ and $\sigma_{\epsilon}$/$\sigma_{\epsilon}^{\rm full}$. As mentioned
above, we consider the following sources of scatter:

  1. (i)

spatial dependence of the mean CHMF on the large scale matter density

  2. (ii)

Poisson sample variance in halo number around the target mean CHMF

  3. (iii)

scatter around the SHMR

  4. (iv)

scatter around the SFMS

  5. (v)

scatter around the FMR

  6. (vi)

scatter in the mapping of the intrinsic luminosity to SFR, $M_{\ast}$ and $Z$

  7. (vii)

scatter in the escape fraction

We need to define a comoving volume over which to sum up the contributions of
galaxies, in order to compute the emissivity PDFs. Here we chose a fiducial
scale of $R_{\rm nl}=5$ cMpc. This is roughly comparable to several relevant
scales during the EoR and CD: (i) the typical HII bubble sizes during the
early-middle stages of reionization (e.g. McQuinn et al., 2007; Lin et al.,
2016); (ii) the resolution of 21cm maps achievable after a 1000h observation
with SKA1-low (Koopmans et al., 2015; Prelogović et al., 2022); (iii) the
Lyman limit mean free path at $z\sim 6$ (e.g. Becker et al., 2021); and (iv)
the field of view of JWST (e.g. Treu et al., 2022; Finkelstein et al., 2023;
Bunker et al., 2023). Our emissivity PDFs are generated from $10000$
realizations of such volumes. In Appendix A we vary this scale and demonstrate
that the estimated mean emissivities have converged to within a few percent.

Figure 5: Distribution of Lyman limit emissivities for regions with a radius
of 5 cMpc. Violin plots correspond to the full emissivity PDFs, while the
crosses and horizontal bars demarcate the mean and 99th percentiles,
respectively. The rectangle on the bottom with the matching number represents
the fraction of 5 cMpc regions with zero emissivity. On the left axis we show
the specific emissivity at the Lyman limit, while on the right axis we show
the corresponding number of ionizing photons (¿13.6 eV) per baryon per Gyr.
The blue shaded region at the top demarcates the approximate criteria for a 5
cMpc to ionize: having an emissivity greater than one (dot dashed line) or two
(dashed line) ionizing photons per baryon in the age of the Universe. Assuming
a threshold value of two ionizing photons per stellar baryon (e.g. Bolton &
Haehnelt 2007; Sobacchi & Mesinger 2014).) we see that roughly half of 5 cMpc
regions can self-ionize by $z\sim$ 7, consistent with the latest estimates of
the EoR history.

###  3.1 Ionizing UV emissivity

In Figure 5 we show the distributions of the ionizing emissivity in our
fiducial model, sampling all of the above-mentioned sources of stochasticity.
On the left axis we report the specific emissivity in units of erg s-1 Hz-1
evaluated at the Lyman limit, while on the right axis we show the total number
of ionizing photons above the Lyman limit per baryon per Gyr. Red violins show
the ionizing emissivity PDF, with crosses (horizontal bars) demarcating the
mean (99% C.L.) of the distributions.The fraction of our $R_{\rm nl}$ = 5 cMpc
realizations that have a zero emissivity is denoted at the bottom of each
violin.

Figure 6: The fractional contribution of different sources of scatter (see
legend) to the mean and standard deviation of the ionizing emissivity,
computed over $R_{\rm nl}=5$ Mpc volumes. The top (bottom) panel shows the
mean (standard deviation) when removing one source of stochasticity,
normalized to the fiducial value that includes all scatter.

As galaxies become rarer towards higher redshifts, the mean ionizing
emissivity decreases and the region to region scatter increases. The
emissivity PDF becomes bimodal, with some regions having an emissivity of zero
while those that have a non-zero emissivity show an approximately log-normal
distribution. At $z\gtrsim 18$ the majority of $R_{\rm nl}$ = 5 cMpc volumes
are expected not to have any galaxies that are actively emitting ionizing
photons, in this fiducial scenario with a log-normal $p(f_{\rm esc})$. If
instead we assume a binomial $p(f_{\rm esc})$ distribution, the majority of
$R_{\rm nl}$ = 5 cMpc have zero ionizing emissivity already by $z\gtrsim 15$.

We now quantify the main sources of scatter driving the variance in Fig. 5. As
discussed above, we do this by repeating our emissivity calculation but
omitting one source of stochasticity (i.e. only using the corresponding mean
relation with no scatter). In Figure 6 we plot the corresponding mean (top
panel) and standard deviation888For numerical stability, we calculate the
standard deviation by fitting a log-normal to the non-zero distribution of
emissivities. (bottom panel), normalized to the corresponding values from the
full calculation shown in Fig. 5. Note that since most of our sources of
scatter are log-normal, assuming a mean relation instead of the full
distribution would underestimate the mean of the emissivity shown in the top
panel (see Appendix B).

The most important source of scatter is the escape fraction, if the escape
fraction is binomial. In this scenario, assuming only the mean escape fraction
for all galaxies would underestimate the standard deviation (std) by 60 – 70%
throughout the EoR and CD (gray dash-dotted curves). However, the mean
emissivity is unchanged (since the bimodal distribution has the same median
and mean; see Appendix B). If we instead assume that the escape fraction is
log-normally distributed (red circled curves), not including scatter in this
quantity only underestimates the mean and std of the emissivity by $\sim
10\%$.

Another important source of stochasticity is the burstiness of star formation.
Assuming all galaxies follow the mean SFMS without scatter (blue dashed
curves) would underpredict the mean (std) of the ionizing emissivity by 40%
(15%) at $z\sim 5$. This underprediction in the std rises to $\sim$ 50%
towards $z\sim 20$, as the typical galaxies have smaller stellar masses and
therefore a broader $p({\rm SFR}\leavevmode\nobreak\ |\leavevmode\nobreak\
M_{\ast})$ (c.f. the bottom panels in Fig. 2). In other words, the increased
”bursty” nature of star formation at higher redshifts (at which the emissivity
is dominated by galaxies with smaller masses) drives a correspondingly larger
spatial variance in the ionizing emissivity.

On the other hand, ignoring scatter around the SHMR results in an
underprediction of the mean and std of the emissivity by only 10%. Other
sources of scatter have a negligible impact on the mean and variance of the
ionizing emissivity. In particular, we note that only $\sim 5$% of our
realizations of 5 Mpc regions at $z=20$ contain fewer than 10 actively star
forming galaxies. Therefore, it is not surprising that Poisson scatter in the
halo number is unimportant in determining emissivities. We note that here we
only consider galaxies above the atomic cooling threshold; had we considered
an additional population of molecular-cooling galaxies, Poisson scatter would
have been even less important since their expected mean number density is much
larger.

###  3.2 X-ray emissivity

Figure 7: Like Fig. 5, but for soft band (0.5–2 keV) X-ray emissivity. Figure
8: Like Fig. 6, but for the soft band (0.5–2 keV) X-ray emissivity.

In Figure 7 we show the distributions of soft-band X-ray emissivities for our
fiducial model, accounting for all sources of stochasticity. Red violins again
represent PDFs averaged over comoving volumes with $R_{\rm nl}=$ 5 cMpc. As
expected the means and widths of the distributions decrease towards higher
redshifts.

We see that the X-ray emissivities have broader distributions compared with
the ionizing emissivities in Fig. 5. For example, at $z\sim 10$ the region-to-
region std of X-ray emissivities is $300\%$ of the mean, while for ionizing
emissivities it is only $50\%$ of the mean. This is primarily due to the fact
that the HMXB LFs that source the X-ray emission in our model are fairly
shallow (see Lehmer et al., 2021). Thus the galaxy-averaged X-ray luminosity
is sensitive to sample variance as it can be determined by a small number of
HMXBs. This is evident by comparing the widths of the conditional
$p(L\leavevmode\nobreak\ |\leavevmode\nobreak\ {\rm SFR},M_{\ast},Z)$
distributions for X-rays and ionizing photons in Fig. 4. The additional
stochasticity in the ionizing emissivity due to the ionizing escape fraction
(assuming it is log-normarly distributed) is sub-dominant compared with the
wider X-ray intrisic luminosity distribution.

We isolate the relative importance of each source of scatter to the X-ray
emissivity in Figure 8. As in the previous section, we show
$\mu_{i}$/$\mu_{\rm full}$ and $\sigma_{i}$/$\sigma_{\rm full}$ in the upper
and lower panels, respectively.

The biggest impact on the mean and std comes from the scatter in the
SFR$-M_{\ast}$ relation (blue dashed curves). Ignoring the scatter around the
SFMS results in an underprediction of the mean (standard deviation) of the
X-ray emissivity by 20% at $z=5$ rising to a factor of $80\%$ at $z=20$. As in
the previous section, this is driven by the mass-dependence of scatter in
SFMS. The physical interpretation is the same: increased burstiness of star
formation in small mass galaxies (that dominate at higher redshifts) boosts
the variance of the X-ray emissivity. The scatter around the SFMS is even more
important for X-ray emissivity, compared with ionizing emissivity, due to the
strong dependence of the intrinsic X-ray luminosity on the SFR (see Fig. 4 and
associated discussion).

Another important source of scatter is the $L_{\rm X}-$SFR relation (violet
dash-dotted curves). The relative difference in standard deviations is roughly
$35\%$ and this value is constant with redshift given the assumption that the
scatter around the scaling relation is constant with mass. At $z=5$ it is more
important than the scatter in SFR$-M_{\ast}$ relation. Note that complex
physics of the formation of binary stars could induce an additional redshift
dependence in this scatter, with different IMF’s giving different populations
of binary stars. This would go in the direction of increasing the importance
of modeling $L_{\rm X}-$SFR scatter at earlier times.

Scatter in the SHMR has a $\sim 10$% effect, again without redshift dependence
since we chose a constant width for $p(M_{\ast}\leavevmode\nobreak\
|\leavevmode\nobreak\ M_{h})$. Scatter in the other terms has a negligible
impact on the X-ray emissivity.

###  3.3 Lyman Werner emissivity

Figure 9: Like Fig. 5, but for the Lyman-Werner (11.2-13.6 eV) emissivity.
Figure 10: Like Fig. 6, but for the Lyman Werner (11.2-13.6 eV) emissivity.

Soft UV photons are important during the cosmic dawn as they regulate the
abundances of H2 (which provides an important cooling channel for the first
galaxies) and the excited spin state of HI (which determines the cosmic 21cm
signal). For concreteness, here we evaluate the emissivity in the Lyman-Werner
band ($11.2$-$13.6$eV) noting that our conclusions would remain the same
regardless of the specific soft UV range of interest.

In Fig. 9 we show the distribution of LW emissivities in 5 cMpc regions. We
see that LW emissivities are more uniform (i.e. with narrower PDFs) than both
ionizing or X-ray emissivities from the previous subsections. This is to be
expected, as the latter bands are sensitive to stochasticity in the ionizing
escape fraction and HMXB LFs, neither of which contribute to the LW
emissivity.

In Fig. 10 we show the fractional contribution of different sources of scatter
to the mean and std of LW emissivity. Again, the SFMS (blue dashed lines) is
the most important contributing source to the variance of emissivity, but less
so compared to X-rays (as could be expected from Fig. 4). At $z\sim 20$ the
burstiness of SFR contributes at a $\sim 50$% level to the std of the
distribution, but this drops to $\sim 10$% at $z\sim 5$.

Also important is the scatter in the SHMR (green full curves) which
contributes at a $\sim 10$% level to the mean and std for all redshifts. Other
sources of scatter are negligible.

##  4 Results: EoR history

Figure 11:  The relative contribution of galaxy stochasticity to the EoR
history. In the left panel we show the common approximation of an EoR history
calculated assuming a constant clumping factor (c.f. Eq. 11), while in the
right we show the fraction of 5 cMpc regions that exceed the threshold of two
ionizing photons per baryon per age of the Universe (c.f. Fig. 5). Black
curves correspond to our fiducial model, including all sources of scatter.
Green / blue curves ignore scatter around the SHMR / SFMS, while the orange
curves ignore all sources of scatter. We see that not accounting for the
burstiness of star formation (i.e. scatter around the SFMS) can result in EoR
histories that are delayed by $\Delta z\sim$ 0.5–1. Not accounting for any
scatter and assuming only mean galaxy properties delays the completion of the
EoR history by $\Delta z\sim$ 2.

The ionizing emissivities shown in the previous section can be used to
estimate the redshift evolution of the volume filling factor of ionized
regions, $Q_{\rm HII}(z)$ – the EoR history. Even though it is an average
quantity, computing the EoR history accurately requires accounting for the
spatial and temporal co-evolution of sources and sinks of ionizing photons,
and is therefore best done numerically. However popular analytic
approximations exist and can provide insight into the relative impact of
scatter in galaxy properties.

Here we compute two proxies for the EoR history. The first is the most common
approximation in the literature, obtained by:

$\frac{\textrm{d}Q_{\rm HII}}{\textrm{d}t}=\dot{n}_{\rm ion/b}-\alpha_{\rm A}\
C\ \langle n_{H}\rangle\ Q_{\rm HII}.$ (11)

Here $\dot{n}_{\rm ion/b}$ is the ionizing emissivity per baryon predicted by
our model999We make the standard assumption that helium is singly ionized by
stellar sources together with hydrogen, due to their comparable ionization
thresholds., $\alpha_{\rm A}$ is the case-A recombination coefficient,
$\langle n_{H}\rangle$ is the mean hydrogen density, and $C\equiv\langle
n_{H}^{2}\rangle/\langle n_{H}\rangle^{2}$ is the so-called ”clumping factor”
computed only over the ionized (not self-shielded) gas. By assuming a constant
clumping factor, this equation ignores the correlation between sources and
sinks of ionizing photons.101010In reality, most recombinations will come from
the earliest patches of the IGM to ionize, which are those with the highest
densities of galaxies. As a result, the growth of HII regions surrounding the
highest galaxy densities begins to stall as reionization progresses, with an
increasing fraction of ionizing photons required to balance recombinations.
This process naturally results in a ”soft landing”, with the recombinations
starting to balance ionizations in the late EoR stages, smoothly transitioning
to the post-EoR regime (e.g. Sobacchi & Mesinger 2014). Estimates of the EoR
history obtained with eq. (11) underpredict the duration of the EoR by $\Delta
z\sim$ 1–2, with the error increasing towards the end stages (see, e.g. Figure
6 in Sobacchi & Mesinger 2014). Here we take $C=2$, noting that we are only
interested in the relative impact of galaxy stochasticity on the EoR history.

We show the resulting estimates in the left panel of Fig. 11. The black curve
corresponds to our fiducial model, in which we account for all of the
aforementioned sources of scatter. The green (blue) curve is computed ignoring
scatter around the mean SHMR (SFMS). The orange curve does not account for any
scatter, taking only the mean values for each relation. We see that scatter in
the SHMR only delays the EoR history by $\Delta z\sim$ 0.1. Ignoring scatter
around the SFMS has a bigger impact, delaying the EoR history by $\Delta
z\sim$ 0.5 – 1. Ignoring scatter in all galaxy properties underestimates the
duration of the EoR and delays the end stages by up to $\Delta z\sim$ 2.

Our second proxy for the EoR history is obtained directly from Fig. 5.
Specifically, we compute the fraction of 5 cMpc regions whose emissivities are
larger than two ionizing photons per baryon per age of the Universe at that
redshift, i.e. $\dot{n}_{\rm ion/b}t_{age}>2$ (shown by the dashed line in Fig
5; c.f. Bolton & Haehnelt 2007; Sobacchi & Mesinger 2014). This approximation
of the EoR history assumes that each 5 cMpc region of the Universe is
instantaneously ionized when this criterion is reached, and that each such
region is independent. However, it does correctly compute the spatial
variation in the emissivity, allowing us to account for local recombinations
by increasing the required ionizing photon threshold. Again, we stress that
here we are only interested in the relative impact of galaxy stochasticity on
the EoR history.

We show the resulting estimate in the right panel of Fig. 11, for the same
models as shown in the left panel. The qualitative evolution of this quantity
is different from the one in the left panel. By its definition, taking only
mean values (orange curve) would result in a step function at the end of the
EoR. Importantly however, the relative impact of ignoring scatter in the SFMS
and SHMR is the similar in both panels.

Regardless of the proxy used in estimating the EoR history, neglecting scatter
around the SFMS results in a delayed EoR history by $\Delta z\sim$ 0.5–1.
Neglecting all galaxy to galaxy scatter by assuming mean values (c.f. Eq. 1)
results in an overextended EoR history and delays its completion by $\Delta
z\sim$ 1–2.

Our results suggest that inferring galaxy properties from EoR history data
without accounting for stochasticity could bias recovery towards brighter
galaxies or higher escape fractions. We will investigate this further in
future work, using 3D simulations that can more accurately capture the
evolution of photon sinks and thus better predict the EoR history.

##  5 Results: UV luminosity functions

The framework developed in Section 2 also allows us to compute the
corresponding galaxy UV LFs. Specifically, for each galaxy realization, we
measure its rest frame magnitude using the luminosity from $1450\AA$ to
$1550\AA$ directly from BPASS (see Section 2.5.2). For simplicity, we do not
account for nebular emission, nor dust attenuation (e.g. Ferrara et al.,
2023). We will include these in future work focused on interpreting UV LFs.

Figure 12: High redshift UV luminosity functions. Our fiducial model including
all of the aforementioned sources of scatter is shown with the red solid lines
(mean values) and surrounding shaded regions (68% C.L.s). The solid green
curves correspond to UV LFs calculated using only mean relations without any
scatter. Also shown in each panel are various observational estimates from HST
and JWST (see text for details).

In Fig. 12 we plot the mean UV LF at each redshift (red line) along with the
68% C.L. (red shaded region). In green, we show the UV LFs calculated assuming
only mean relations without any scatter. Also shown are various observational
estimates from HST (Bouwens et al., 2015 (B15); Bouwens et al., 2016 (B16);
Bouwens et al., 2021 (B21); Livermore et al., 2017 (L17); Ishigaki et al.,
2018 (I18); Oesch et al., 2016 (O16); Oesch et al., 2018 (O18); Leethochawalit
et al., 2023 (L22); Kauffmann et al., 2022 (K22)) and JWST (Naidu et al., 2022
(N22); Finkelstein et al., 2023 (F22); Donnan et al., 2023 (D23); Donnan et
al., 2024 (D24); Pérez-González et al., 2023 (P23); Robertson et al., 2023
(R24); Harikane et al., 2024 (H23); McLeod et al., 2024 (M24); Willott et al.,
2024 (W24)).

Comparing the green and red curves, we see that including scatter shifts the
mean to brighter magnitudes and flattens the UV LFs. This is a well-known
effect of upscattering some fraction of the more abundant faint galaxies to
brighter magnitudes. For our fiducial model, the shift is roughly 1-2
magnitudes at $M_{\rm UV}\sim$ -18, consistent with other estimates of the
impact of stochasticity on UV LFs (e.g. Mason et al., 2023; Shen et al., 2023;
Gelli et al., 2024).

Comparing the red curve to the observational estimates, we see that our
fiducial model is consistent with UV LFs at $z\lesssim 10$. However, the mean
underpredicts the recent estimates of UV LFs at $z\gtrsim 11$ from broad-band
JWST photometry. Although the observational data points are mostly within the
68% C.L. of our model, they are systematically higher. Assuming there is no
observational bias and that there are no correlations between the magnitude
bins, this systematic underprediction would imply that our fiducial model is
strongly disfavored by the data at $z\gtrsim 11$. This is qualitatively
consistent with previous conclusions from the literature that larger than
expected levels of scatter would be required to explain JWST results, provided
the $z>10$ photometric estimates are accurate (e.g. Mirocha & Furlanetto 2023;
Mason et al. 2023; Shen et al. 2023; Pallottini & Ferrara 2023; Gelli et al.
2024). Alternately, a correlation between observational estimates in different
magnitude bins (for example through cosmic variance; e.g. Willott et al.,
2024) could alleviate this apparent tension.

The empirical framework we developed here is very flexible, and allows us to
explicitly define the mean and scatter in every fundamental relation that
leads to the 1500 Å UV magnitude. In future work we will use our model
combined with physically-motivated priors to infer these conditional
distributions from JWST UV LFs and other observational data.

##  6 Conclusions

Here we quantify how does the galaxy to galaxy scatter in their properties
impact estimates of their emissivities and related observables. We use a semi-
empirical model that explicitly defines scatter around well-studied mean
relations: (i) the conditional halo mass function (CHMF); (ii) the stellar-to-
halo mass relation (SHMR); (iii) galaxy star formation main sequence (SFMS);
(iv) fundamental metallicity relation (FMR); (v) conditional intrinsic
luminosity; and (vi) photon escape fraction. We compute the corresponding
multi-frequency (ionizing UV, X-rays, LW) emissivities, EoR histories, and UV
LFs, quantifying the relative importance of the above sources of scatter.

We find that the burstyness of star formation (i.e. scatter around the mean
SFMS) is important for all emissivities. Because we assume burstyness
increases towards smaller mass galaxies, the scatter around the SFMS becomes
increasingly important at higher redshifts. Neglecting this source of
stochasticity could underpredict the mean and std of emissivities by factors
of up to few-10 during the EoR and CD. Stochasticity in the ionizing escape
fraction can dominate the spatial scatter in the ionizing emissivity if its
distribution is binomial. If instead the escape fraction is log-normally
distributed, its contribution to the total emissivity scatter is only of order
$\sim$ 10%. For the X-ray emissivity, one must account for scatter in the
intrinsic luminosity, which in our fiducial model is driven by high mass X-ray
binary luminosity functions.

We find that neglecting stochasticity overestimates the duration of
reionization, delaying its completion by $\Delta z\sim$ 1–2. Neglecting only
scatter around the mean SFMS results in a delay of the EoR history by $\Delta
z\sim$ 0.5 – 1. This suggests that inferring galaxy properties from EoR
history data without accounting for stochasticity could bias recovery towards
brighter galaxies or higher escape fractions.

We recover the well-known effect of stochasticity flattening the UV LFs. In
our fiducial model, this results in a shift of 1–2 UV magnitudes at $M_{\rm
UV}\sim$ -18. Our UV LFs are consistent with observational data at $z\leq 10$
but consistently under-predict recent estimates at at higher redshifts. This
is qualitatively in line with other studies, and implies that larger scatter
is required in order for it to be the sole explanation for photometric
estimates at $z>10$.

We conclude that models of the EoR and CD should at least account for scatter
around the SFMS. Simulating the X-ray background during these epochs (for
example when computing the 21cm signal) additionally requires accounting for
scatter in the intrinsic X-ray luminosities of galaxies.

The semi-empirical framework we use here is flexible and transparent. It can
easily be extended to accommodate additional observables, different functional
distributions, and/or dependencies on additional galaxy properties.

###### Acknowledgements.

We thank Yuxiang Qin for the help with collecting UV LF observations. We
gratefully acknowledge computational resources of the Center for High
Performance Computing (CHPC) at SNS. A.M. acknowledges support from the
Italian Ministry of Universities and Research (MUR) through the PRIN project
”Optimal inference from radio images of the epoch of reionization” as well as
the PNRR project ”Centro Nazionale di Ricerca in High Performance Computing,
Big Data e Quantum Computing”.

## References

  * Abdurashidova et al. (2022) Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 924, 51
  * Allende Prieto et al. (2018) Allende Prieto, C., Koesterke, L., Hubeny, I., et al. 2018, A&A, 618, A25
  * Asplund et al. (2004) Asplund, M., Grevesse, N., Sauval, A. J., Allende Prieto, C., & Kiselman, D. 2004, A&A, 417, 751
  * Barkana & Loeb (2004) Barkana, R. & Loeb, A. 2004, ApJ, 609, 474
  * Barrow et al. (2020) Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJ, 902, L39
  * Barrow et al. (2017) Barrow, K. S. S., Wise, J. H., Norman, M. L., O’Shea, B. W., & Xu, H. 2017, MNRAS, 469, 4863
  * Basu-Zych et al. (2013) Basu-Zych, A. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2013, ApJ, 774, 152
  * Becker et al. (2021) Becker, G. D., D’Aloisio, A., Christenson, H. M., et al. 2021, MNRAS, 508, 1853
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143
  * Bird et al. (2022) Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703
  * Bolton & Haehnelt (2007) Bolton, J. S. & Haehnelt, M. G. 2007, MNRAS, 382, 325
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34
  * Bouwens et al. (2016) Bouwens, R. J., Oesch, P. A., Labbé, I., et al. 2016, ApJ, 830, 67
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47
  * Bouwens et al. (2023) Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Brorby et al. (2016) Brorby, M., Kaaret, P., Prestwich, A., & Mirabel, I. F. 2016, MNRAS, 457, 4081
  * Bunker et al. (2023) Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, 677, A88
  * Byrne & Stanway (2023) Byrne, C. M. & Stanway, E. R. 2023, MNRAS, 521, 4995
  * Byrne et al. (2022) Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329
  * Cen & Kimm (2015) Cen, R. & Kimm, T. 2015, ApJ, 801, L25
  * Ceverino et al. (2018) Ceverino, D., Klessen, R. S., & Glover, S. C. O. 2018, MNRAS, 480, 4842
  * Ciardi et al. (2003) Ciardi, B., Ferrara, A., & White, S. D. M. 2003, MNRAS, 344, L7
  * Curti et al. (2023) Curti, M., Maiolino, R., Carniani, S., et al. 2023, arXiv e-prints, arXiv:2304.08516
  * Curti et al. (2020) Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944
  * Curtis-Lake et al. (2021) Curtis-Lake, E., Chevallard, J., Charlot, S., & Sandles, L. 2021, MNRAS, 503, 4855
  * Das et al. (2017) Das, A., Mesinger, A., Pallottini, A., Ferrara, A., & Wise, J. H. 2017, MNRAS, 469, 1166
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827
  * Davies et al. (2023a) Davies, F. B., Bosman, S. E. I., Gaikwad, P., et al. 2023a, arXiv e-prints, arXiv:2312.08464
  * Davies & Furlanetto (2016) Davies, F. B. & Furlanetto, S. R. 2016, MNRAS, 460, 1328
  * Davies & Furlanetto (2022) Davies, F. B. & Furlanetto, S. R. 2022, MNRAS, 514, 1302
  * Davies et al. (2023b) Davies, J. E., Bird, S., Mutch, S., et al. 2023b, MNRAS, 525, 2553
  * Dayal & Ferrara (2018) Dayal, P. & Ferrara, A. 2018, Phys. Rep, 780, 1
  * Di Cesare et al. (2023) Di Cesare, C., Graziani, L., Schneider, R., et al. 2023, MNRAS, 519, 4632
  * Dixon et al. (2016) Dixon, K. L., Iliev, I. T., Mellema, G., Ahn, K., & Shapiro, P. R. 2016, MNRAS, 456, 3011
  * Donnan et al. (2023) Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011
  * Donnan et al. (2024) Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, arXiv e-prints, arXiv:2403.03171
  * Eide et al. (2018) Eide, M. B., Graziani, L., Ciardi, B., et al. 2018, MNRAS, 476, 1174
  * Ferrara et al. (2023) Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986
  * Finkelstein et al. (2023) Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13
  * Fornasini et al. (2019) Fornasini, F. M., Kriek, M., Sanders, R. L., et al. 2019, ApJ, 885, 65
  * Fragos et al. (2023) Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45
  * Fragos et al. (2013) Fragos, T., Lehmer, B. D., Naoz, S., Zezas, A., & Basu-Zych, A. 2013, ApJ, 776, L31
  * Furlanetto (2006) Furlanetto, S. R. 2006, MNRAS, 371, 867
  * Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1
  * Gaikwad et al. (2023) Gaikwad, P., Haehnelt, M. G., Davies, F. B., et al. 2023, MNRAS, 525, 4093
  * Garcia et al. (2024) Garcia, A. M., Torrey, P., Ellison, S., et al. 2024, arXiv e-prints, arXiv:2403.08856
  * Gelli et al. (2024) Gelli, V., Mason, C., & Hayward, C. C. 2024, arXiv e-prints, arXiv:2405.13108
  * Gillet et al. (2020) Gillet, N. J. F., Mesinger, A., & Park, J. 2020, MNRAS, 491, 1980
  * Grazian et al. (2017) Grazian, A., Giallongo, E., Paris, D., et al. 2017, A&A, 602, A18
  * Haiman et al. (2000) Haiman, Z., Abel, T., & Rees, M. J. 2000, ApJ, 534, 11
  * Harikane et al. (2024) Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56
  * Harikane et al. (2016) Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123
  * Hassan et al. (2021) Hassan, S., Davé, R., McQuinn, M., et al. 2021, arXiv e-prints, arXiv:2109.03840
  * Holzbauer & Furlanetto (2012) Holzbauer, L. N. & Furlanetto, S. R. 2012, MNRAS, 419, 718
  * Ishigaki et al. (2018) Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73
  * Izotov et al. (2016) Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016, MNRAS, 461, 3683
  * Kannan et al. (2022) Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005
  * Kauffmann et al. (2022) Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65
  * Kaur et al. (2022) Kaur, H. D., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 513, 5097
  * Kaurov et al. (2018) Kaurov, A. A., Venumadhav, T., Dai, L., & Zaldarriaga, M. 2018, ApJ, 864, L15
  * Koopmans et al. (2015) Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 1
  * Kostyuk et al. (2023) Kostyuk, I., Nelson, D., Ciardi, B., Glatzle, M., & Pillepich, A. 2023, MNRAS, 521, 3077
  * Kreilgaard et al. (2024) Kreilgaard, K. C., Mason, C. A., Cullen, F., Begley, R., & McLure, R. J. 2024, arXiv e-prints, arXiv:2405.10364
  * Leethochawalit et al. (2023) Leethochawalit, N., Roberts-Borsani, G., Morishita, T., Trenti, M., & Treu, T. 2023, MNRAS, 524, 5454
  * Lehmer et al. (2016) Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7
  * Lehmer et al. (2021) Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2021, ApJ, 907, 17
  * Lin et al. (2016) Lin, Y., Oh, S. P., Furlanetto, S. R., & Sutter, P. M. 2016, MNRAS, 461, 3361
  * Livermore et al. (2017) Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113
  * Lovell et al. (2021) Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127
  * Ma et al. (2018) Ma, X., Hopkins, P. F., Garrison-Kimmel, S., et al. 2018, MNRAS, 478, 1694
  * Madau & Fragos (2017) Madau, P. & Fragos, T. 2017, ApJ, 840, 39
  * Magg et al. (2022) Magg, M., Reis, I., Fialkov, A., et al. 2022, MNRAS, 514, 4433
  * Mannucci et al. (2010) Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115
  * Mascia et al. (2023) Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155
  * Mason et al. (2023) Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, 521, 497
  * McLeod et al. (2024) McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, 527, 5004
  * McQuinn (2016) McQuinn, M. 2016, ARA&A, 54, 313
  * McQuinn et al. (2007) McQuinn, M., Lidz, A., Zahn, O., et al. 2007, MNRAS, 377, 1043
  * Meriot & Semelin (2024) Meriot, R. & Semelin, B. 2024, A&A, 683, A24
  * Mesinger (2016) Mesinger, A. 2016, Astrophysics and Space Science Library, Vol. 423, Understanding the Epoch of Cosmic Reionization
  * Mesinger et al. (2013) Mesinger, A., Ferrara, A., & Spiegel, D. S. 2013, MNRAS, 431, 621
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955
  * Mineo et al. (2012) Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 419, 2095
  * Mirocha & Furlanetto (2023) Mirocha, J. & Furlanetto, S. R. 2023, MNRAS, 519, 843
  * Mirocha et al. (2017) Mirocha, J., Furlanetto, S. R., & Sun, G. 2017, MNRAS, 464, 1365
  * Mirocha et al. (2021) Mirocha, J., La Plante, P., & Liu, A. 2021, MNRAS, 507, 3872
  * Mo & White (1996) Mo, H. J. & White, S. D. M. 1996, MNRAS, 282, 347
  * Muñoz et al. (2022) Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657
  * Mutch et al. (2016) Mutch, S. J., Geil, P. M., Poole, G. B., et al. 2016, MNRAS, 462, 250
  * Naidu et al. (2022) Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14
  * Nikolić et al. (2023) Nikolić, I., Mesinger, A., Qin, Y., & Gorce, A. 2023, MNRAS, 526, 3170
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105
  * Oesch et al. (2016) Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129
  * Oh (2001) Oh, S. P. 2001, ApJ, 553, 499
  * Pacucci et al. (2014) Pacucci, F., Mesinger, A., Mineo, S., & Ferrara, A. 2014, MNRAS, 443, 678
  * Pahl et al. (2023) Pahl, A. J., Shapley, A., Steidel, C. C., et al. 2023, MNRAS, 521, 3247
  * Pallottini & Ferrara (2023) Pallottini, A. & Ferrara, A. 2023, A&A, 677, L4
  * Pallottini et al. (2022) Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621
  * Park et al. (2019) Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933
  * Pérez-González et al. (2023) Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJ, 951, L1
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Popesso et al. (2022) Popesso, P., Concas, A., Cresci, G., et al. 2022, arXiv e-prints, arXiv:2203.10487
  * Prelogović et al. (2022) Prelogović, D., Mesinger, A., Murray, S., Fiameni, G., & Gillet, N. 2022, MNRAS, 509, 3852
  * Press & Schechter (1974) Press, W. H. & Schechter, P. 1974, ApJ, 187, 425
  * Qin et al. (2021) Qin, Y., Mesinger, A., Bosman, S. E. I., & Viel, M. 2021, MNRAS, 506, 2390
  * Reis et al. (2022) Reis, I., Barkana, R., & Fialkov, A. 2022, MNRAS, 511, 5265
  * Robertson et al. (2023) Robertson, B., Johnson, B. D., Tacchella, S., et al. 2023, arXiv e-prints, arXiv:2312.10033
  * Ross et al. (2017) Ross, H. E., Dixon, K. L., Iliev, I. T., & Mellema, G. 2017, MNRAS, 468, 3785
  * Santini et al. (2017) Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76
  * Schaeffer et al. (2023) Schaeffer, T., Giri, S. K., & Schneider, A. 2023, MNRAS, 526, 2942
  * Shekhar Murmu et al. (2023) Shekhar Murmu, C., Datta, K. K., Majumdar, S., & Greve, T. R. 2023, arXiv e-prints, arXiv:2311.17062
  * Shen et al. (2023) Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254
  * Sheth (1998) Sheth, R. K. 1998, MNRAS, 300, 1057
  * Sheth & Tormen (1999) Sheth, R. K. & Tormen, G. 1999, MNRAS, 308, 119
  * Sobacchi & Mesinger (2013) Sobacchi, E. & Mesinger, A. 2013, MNRAS, 432, L51
  * Sobacchi & Mesinger (2014) Sobacchi, E. & Mesinger, A. 2014, MNRAS, 440, 1662
  * Stanway & Eldridge (2018) Stanway, E. R. & Eldridge, J. J. 2018, MNRAS, 479, 75
  * Stefanon et al. (2021) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29
  * Steidel et al. (2018) Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123
  * Strom et al. (2018) Strom, A. L., Steidel, C. C., Rudie, G. C., Trainor, R. F., & Pettini, M. 2018, ApJ, 868, 117
  * Tacchella et al. (2016) Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790
  * Trapp & Furlanetto (2020) Trapp, A. C. & Furlanetto, S. R. 2020, MNRAS, 499, 2401
  * Treu et al. (2022) Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110
  * Wechsler & Tinker (2018) Wechsler, R. H. & Tinker, J. L. 2018, ARA&A, 56, 435
  * Willott et al. (2024) Willott, C. J., Desprez, G., Asada, Y., et al. 2024, ApJ, 966, 74
  * Wolcott-Green et al. (2011) Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2011, MNRAS, 418, 838
  * Xu et al. (2014) Xu, H., Ahn, K., Wise, J. H., Norman, M. L., & O’Shea, B. W. 2014, ApJ, 791, 110
  * Xu et al. (2016) Xu, H., Wise, J. H., Norman, M. L., Ahn, K., & O’Shea, B. W. 2016, ApJ, 833, 84
  * Yeh et al. (2022) Yeh, J. Y. C., Smith, A., Kannan, R., et al. 2022, arXiv e-prints, arXiv:2205.02238
  * Yung et al. (2019) Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019, MNRAS, 490, 2855
  * Zaroubi (2013) Zaroubi, S. 2013, in Astrophysics and Space Science Library, Vol. 396, The First Galaxies, ed. T. Wiklind, B. Mobasher, & V. Bromm, 45

##  Appendix A Convergence of the mean emissivity with scale

Figure 13: The mean of the ionizing emissivity computed over regions of
varying scales, $R_{\rm nl}$ = 1, 2, 5, 10 cMpc. Our fiducial choice of 5 cMpc
has converged in the mean to within a few percent. The analytical excursion
set model misses massive halos when conditioned smaller scales, resulting in
an underprediction of the mean for $R_{\rm nl}$ = 1 cMpc by factors of few –
10.

Here we confirm that our fiducial choice of $R_{\rm nl}$ = 5 cMpc when
computing emissivity distributions converges to the correct mean. In
principle, the averaging over the scale-dependent overdensity distribution,
$p_{z}(\delta_{0}|R_{\rm nl})$, in the top row of Eq. LABEL:eq:full_emissivity
should self-consistently ensure that the correct mean is recovered, regardless
of the choice of scale. In practice however, conditional halo mass functions
underpredict the numbers of relatively massive halos whose Lagrangian volumes
are close to the conditioning scale.

In Figure 13 we illustrate how the mean emissivity changes with scale. We use
the ionizing emissivity to illustrate the trend; however, the result is the
same for the other two bands of interest. We chose three additional scales: 1,
2, 10 Mpc. The first two choices roughly span the cell sizes used in semi-
numerical (e.g. Mesinger et al. 2011; Muñoz et al. 2022; Schaeffer et al.
2023), or low-resolution (e.g. Dixon et al. 2016; Meriot & Semelin 2024)
radiative transfer simulations of the EoR/CD, while the latter roughly
corresponds to the HII bubble size late in the EoR.

Comparing the red and black curves, we see that the mean for our fiducial
choice of 5 cMpc has converged to within a few percent. We are thus reassured
that our fiducial choice can be used to predict global quantities like the EoR
history and UV LFs.

As the scale is further reduced, we see that the mean emissivity can be
significantly underestimated when using conditional excursion set formalism.
For example, the mean using $R_{\rm nl}=$ 1 cMpc is underestimated by factors
of few - ten. This serves as a caution against computing halo fields only at
the cell level for low-resolution EoR/CD simulations (e.g. Appendix A in
Davies & Furlanetto 2022; Reis et al. 2022). Instead, $N$-body (e.g. Dixon et
al. 2016; Schaeffer et al. 2023; Meriot & Semelin 2024), excursion-set that
accounts for larger scales (e.g. Furlanetto et al. 2004; Mesinger et al. 2011;
Cen & Kimm 2015), or a mixture of the two (e.g. McQuinn et al. 2007, Davies et
al. in prep) should be used.

##  Appendix B Shift in the mean emissivity for correlated log-normal
distributions

In Section 2 we wrote the emissivity in the form:

$\overline{\varepsilon}=\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h},z)}{\textrm{d}M_{h}}L(M_{h})$
(12)

where we ignore the escape fraction for the moment. As we mentioned in that
section, this formula holds if galaxy properties are deterministic functions
of halo mass. However, we know that is not the case and for that reason we
wrote the general formula for the mean in Eq. LABEL:eq:full_emissivity. In our
case, we assume log-normal distribution for most of the scaling relations
(except for the ones relating to the halo abundances) so we can analytically
integrate over some of the distributions in Eq. LABEL:eq:full_emissivity. For
simplicity, in this section we reduce the dependencies of luminosity to only
$L(M_{h})$ with an appropriate log-normal PDF $p(\log
L|M_{h})=\mathcal{N}(\mu_{L},\sigma_{L})$ without loss of generality. The Eq.
LABEL:eq:full_emissivity becomes:

$\displaystyle\begin{split}\overline{\varepsilon}=&\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h})}{\textrm{d}M_{h}}\int\textrm{d}L\
L\ p(\log_{10}L|\log_{10}M_{h})\\\
=&\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h})}{\textrm{d}M_{h}}\int\textrm{d}L\
L\ \frac{1}{\sqrt{2\pi}\log_{e}{10}\
\sigma_{L}L}\exp\left(-\frac{(\log_{10}L-\mu_{L})^{2}}{2\sigma_{L}^{2}}\right)\end{split}$
(13)

The integral on the right can be analytically computed:

$\overline{\varepsilon}=\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h})}{\textrm{d}M_{h}}10^{\left(\mu_{L}+\frac{\log_{e}{10}\
\sigma_{L}^{2}}{2}\right)}$ (14)

Since the halo mass function in general is not an analytic function of halo
mass and mean of the luminosity scaling depends on the halo mass, this
integral cannot be computed analytically. However we can already gain
intuition about the mean looking at the last term of Eq. 14. The factor
$10^{\mu_{L}}$ corresponds to the mean of the $L(M_{h})$ scaling relation,
i.e. what one would obtain if one did not consider scatter around the mean.
The second part, $10^{\frac{\log_{e}10\ \sigma_{L}^{2}}{2}}$ represents the
shift of the mean when integrating over the whole PDF of the distribution.
This is proportional to the width of the distribution, indicating that the
wider the distribution, the larger the corresponding shift in the mean. This
is a general property of asymmetric distributions like the log-normal and has
important implications for interpreting means of scaling relations. This is
clearly seen in Figures. 6,8 and 10 where removing one source of scatter
reduces the mean proportionally to the width of the distribution.

If we instead add an additional term that is binomially distributed (e.g. one
choice for the escape fraction in Sec. 2.6) then the equation becomes:

$\displaystyle\begin{split}\overline{\varepsilon}=\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h})}{\textrm{d}M_{h}}\int\textrm{d}L\
L\ &p(\log_{10}L|\log_{10}M_{h})\times\\\ &\times\int\textrm{d}f_{\rm
esc}f_{\rm esc}{{n}\choose{k}}\mathcal{P}^{n}(1-\mathcal{P})^{k}\end{split}$
(15)

where we have explicitly written out the binomial distribution. In the above,
$n=1$, $k=0$ and $\mathcal{P}=\overline{f}_{\rm esc}$ so the distribution
trivially becomes:

$\displaystyle\begin{split}\overline{\varepsilon}=&\int\textrm{d}M_{h}\frac{\textrm{d}n(M_{h})}{\textrm{d}M_{h}}\int\textrm{d}L\
L\ p(\log_{10}L|\log_{10}M_{h})\overline{f}_{\rm esc}\end{split}$ (16)

Therefore, the mean does not change if scatter is added following a binomial
distribution. This is clearly seen in Fig. 6.