EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 692 (December 2024) A&A, 692 (2024) A142 Full HTML
Open Access
Issue
A&A
Volume 692, December 2024
Article Number A142
Number of page(s) 15
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202451213
Published online 13 December 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 they 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 allowed 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 and 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 (that is, 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 and stellar to halo mass relation) can give biased results for correlated distributions (e.g., the average of a product is not the same as the product of the averages; cf. Appendix B for simple examples). Moreover, some measurements (e.g., 21-cm interferometry, Lyα forest and kinetic Synaev-Zel’dovich signal) 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 21 cm 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; Sun et al. 2023; Kravtsov & Belokurov 2024).

For this work we constructed 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); and (iii) Lyman Werner (which determines when H2 cooling stops being efficient in the first galaxies). We computed the distribution of these multi-frequency emissivities as a function of scale and redshift. Our model sampled the largest expected sources of stochasticity, including the following: the abundance of dark-matter halos, the stellar-to-halo mass relation, the galaxy main sequence, the fundamental metallicity relation, luminosity, and escape fraction scalings. We quantified the relative importance of each term to the total emissivity in each of the considered bands. We also evaluated the importance of these stochastic terms for simple estimates of the EoR history, as well as the high redshift UV luminosity functions (LFs). Our results can be used to improve estimates of cosmic radiation fields and guide models of the EoR and 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, and 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 ΛCDM cosmological parameters (Ωm, Ωb, ΩΛ, h, σ8, ns = 0.310, 0.049, 0.689, 0.677, 0.81, 0.963), consistent with the latest estimates from Planck Collaboration VI (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:

ε i ( z ) = ∫ d M h d n ( M h , z ) d M h L i ( M h , z ) f esc , i ( M h , z ) . $$ \begin{aligned} \varepsilon _i(z) = \int \mathrm{d}M_h \frac{\mathrm{d}n(M_h, z)}{\mathrm{d}M_h} L_i(M_h, z) \ f_{\rm esc, i}(M_h, z) ~ . \end{aligned} $$(1)

Here d n d M h $ \frac{\mathrm{d}n}{\mathrm{d}M_{\mathrm{h}}} $ is the number density of halos per unit mass – that is, the halo mass function (HMF) – Li is the intrinsic luminosity of a galaxy hosted in a halo of mass Mh at redshift z, and fesc, 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)πRnl3, as follows:1

ε ¯ i ( R nl , z ) = ∫ d M h ∫ d δ 0 d n ( M h , z | R 0 , δ 0 ) d M h p z ( δ 0 | R nl ) ( HMF ) × ∫ d M ∗ p ( M ∗ | M h ) ( SHMR ) × ∫ dSFR p z ( SFR | M ∗ ) ( SFMS ) × ∫ d Z p z ( Z | SFR , M ∗ ) ( FMR ) × ∫ d L i L i p ( L i | SFR , Z ) ( L ) × ∫ d f esc , i f esc , i p ( f esc , i ) . ( EF ) $$ \begin{aligned} \bar{\varepsilon }_i(R_{\rm nl}, z) =&\int \mathrm{d}M_h \int \mathrm{d}\delta _0 \frac{\mathrm{d}n (M_h, z ~|~R_0, \delta _0)}{\mathrm{d}M_h} p_z(\delta _0 ~|~ R_{\rm nl}) \qquad (\mathrm{HMF})\nonumber \\&\quad \times \int \mathrm{d}M_{*} \ p(M_{*}~|~M_{\rm h}) \qquad (\mathrm{SHMR}) \\&\quad \times \int \mathrm{dSFR} \ p_z(\mathrm{SFR}~|~M_*) \qquad {(\mathrm{SFMS})} \nonumber \\&\quad \times \int \mathrm{d}Z \ p_z(Z ~ |~ \mathrm{SFR}, M_{*}) \ \qquad {(\mathrm{FMR})} \nonumber \\&\quad \times \int \ \mathrm{d}L_i \ L_i \ p(L_i~|~\mathrm{SFR}, Z) \\ \qquad {(\mathrm{L})}&\quad \times \int \mathrm{d}f_{\rm esc, i} \ f_{\rm esc, i} \ p(f_{\rm esc, i}). \qquad {(\mathrm{EF})} \end{aligned} $$(2)

Here δ0(R0) is the linear matter overdensity of a spherical volume of Lagrangian radius R0 corresponding to the final Eulerian radius Rnl, p(M* | Mh) is the conditional probability of stellar mass M* for a given Mh, pz(SFR | M*) is the conditional probability of a star formation rate (SFR) for a given stellar mass2, pz(Z | SFR, M*) is the conditional probability of a stellar metallicity Z for a given SFR and M*, p(Li | SFR, Z) is the conditional probability of a luminosity Li in a given wavelength band i for a given SFR and Z, and p(fesc, i) is the probability of an escape fraction fesc, i in band i3. Loosely speaking, the running averages of the conditional probabilities in the first four rows are commonly referred to as the halo mass function (HMF)4, stellar-to-halo mass relation (SHMR), star forming main sequence (SFMS) of galaxies, and fundamental metallicity relation (FMR); while (L) and (EF) represent running averages of conditional distributions of the luminosity and escape fraction at wavelength band i, respectively. 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 ⟨ ( ε i − ε ¯ i ) 2 ⟩ $ \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 (2) 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. (2) 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. (2022) 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 Murmu et al. (2024) 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. 2024, Qin et al. in prep.).

Algorithm 1.Computing a single realization of the emissivity in spectral band i of a region of radius Rnl at redshift z: ε ~ i ( R nl , z ) $ \tilde{\varepsilon}_i(R_{\mathrm{nl}}, z) $

1: sample the linear matter overdensity δ ~ 0 ∼ p z ( δ 0 | R nl ) $ \tilde{\delta}_0 \sim p_z(\delta_0~|~R_{\mathrm{nl}}) $

2: obtain a realization of the halo field by sampling the CHMF: { M h j } ~ ∼ d n ( M h , z | δ ~ 0 , R 0 ) / d M h $ \widetilde{\{M_{h}^{\mathrm{j}}\}} \sim \mathrm{d}n (M_h, z ~|~ \tilde{\delta}_0, R_0)/\mathrm{d}M_h $

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

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

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

6:  CONTINUE

7:  end if

8:  sample stellar mass M ~ * j $ \tilde{M}_{\ast}^j $ ∼ p ( M * | M ~ h j ) $ \sim p(M_{\ast}~|~\tilde{M}_{\mathrm{h}}^j) $

9:  sample star formation rate SFR ~ j $ \tilde{{\mathrm{SFR}}}^j $ ∼ p z ( SFR | M ~ * j ) $ \sim p_z(\mathrm{SFR}~|~\tilde{M}_{\ast}^j) $

10:  sample metallicity Z ~ j $ \tilde{Z}^j $ ∼ p z ( Z | SFR ~ j , M ~ * j ) $ \sim p_z(Z~|~\tilde{\mathrm{SFR}}^j,\tilde{M}_{\ast}^j) $

11:  sample intrinsic luminosity L ~ i j $ \tilde{L}_i^j $ ∼ p ( L i | SFR ~ j , Z ~ j ) $ \sim p(L_i~|~\tilde{\mathrm{SFR}}^j, \tilde{Z}^j) $

12:  sample escape fraction f ~ esc , i j $ \tilde{f}_{\mathrm{esc,i}}^j $ ∼p(fesc, i)

13: end for

14: ε i ~ = ∑ halo = j L ~ i j $ \tilde{\varepsilon_{i}} = \sum_{\mathrm{halo}=j} \tilde{L}_i^j $ f ~ esc , i j $ \tilde{f}_{\mathrm{esc,i}}^j $

We could analytically derive the emissivity distribution, p(εi | Rnl, z), if the conditional distributions in Eqs. (HMF)–(EF) followed simple Gaussian forms. In the more general case, we can solve for p(εi | Rnl, 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 and describe each step in turnin Algorithm 1.

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(εi|Rnl, 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):

d n d M h ( M h , z | δ 0 , R 0 ) = f ¯ ST f ¯ PS 2 π ρ ¯ M h × δ c − δ 0 σ 2 ( M h ) − σ 2 ( R 0 ) | d σ d M h | exp ( − ( δ c − δ 0 ) 2 σ 2 ( M h ) − σ 2 ( R 0 ) ) $$ \begin{aligned}&\frac{\mathrm{d}n}{\mathrm{d}M_{\rm h}} (M_{\rm h}, z ~|~ \delta _0, R_{0}) = \frac{\overline{f}_{\rm ST}}{\overline{f}_{\rm PS}} \sqrt{\frac{2}{\pi }} \frac{\overline{\rho }}{M_{\rm h}}\nonumber \\&\qquad \times \frac{\delta _{\rm c}- \delta _0}{\sigma ^2(M_h) - \sigma ^2(R_{0})} \left|\frac{\mathrm{d}\sigma }{\mathrm{d}M_{\rm h}}\right| \exp {\left(-\frac{(\delta _{\rm c}- \delta _0)^2}{\sigma ^2(M_h) - \sigma ^2(R_{0})}\right)} \end{aligned} $$(3)

In the above, f ¯ ST $ \overline{f}_{\mathrm{ST}} $ and f ¯ PS $ \overline{f}_{\mathrm{PS}} $ correspond to the mean Sheth-Tormen and Press-Schechter collapsed fractions above the atomic cooling threshold of Tvir ≥ 104 K, respectively, δc(z) is the critical linear density from the spherical collapse model, and σ2(M)

is the mass variance of the Lagrangian (linear) density field on scales M = (4/3)πR03.

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(δ0) follows a zero-mean Gaussian distribution whose width is determined by σ2(R0). We transform this distribution to Eulerian space using the spherical collapse model (e.g., Mo & White 1996):

p ( δ nl | R nl , z ) d δ nl = 1 1 + δ 0 f R ( σ 2 | R 0 , z ) d σ 2 , $$ \begin{aligned} p(\delta _{\rm nl} ~|~ R_{\rm nl}, z) \mathrm{d}\delta _{\rm nl} = \frac{1}{1+\delta _{0}} f_{\rm R} (\sigma ^2 ~|~ R_0, z) \mathrm{d}\sigma ^2 ~, \end{aligned} $$(4)

δ 0 = − 1.35 ( 1 + δ nl ) − 2 / 3 + 0.78785 ( 1 + δ nl ) − 0.58661 − 1.12431 ( 1 + δ nl ) − 1 / 2 + 1.68647 . $$ \begin{aligned} \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 ~. \nonumber \end{aligned} $$(5)

Here fR is the first-crossing distribution from Sheth (1998), and δ nl = ρ / ρ ¯ − 1 $ \delta_{\mathrm{nl}} = \rho/\bar{\rho} - 1 $ is the Eulearian (nonlinear) overdensity.

With the above relations, we generate a Lagrangian overdensity sample, δ ∼ 0 ∼ p ( δ 0 | R nl , z ) $ \tilde{\delta}_0 \sim p(\delta_0~|~R_{\mathrm{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 N ∼ ( > M min | δ ∼ 0 , R 0 ) $ \tilde{N}({ > }M_{\mathrm{min}} ~|~ \tilde{\delta}_0, R_0) $, by assuming a Poisson distribution whose mean is given by the integral of Eq. (3) from Mmin 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 ±10% of the target mean5.

We show some example CHMFs in Fig. 1. The dashed black curves correspond to the target mean CHMFs in regions of Eularian scale Rnl, at mean density for infinite realizations, 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, while the red shaded region corresponds to the 95% C.L. of the halo field conditioned on a region of scale Rnl = 5 Mpc. The impact of stochasticity in the red curves is very evident as the mean number decreases, that is, toward high redshifts and high masses. In the top panel we also show the mean (non conditional) HMF (that is, the limit as R0 → ∞ at δ0 = 0).

thumbnail Fig. 1.

Example halo mass functions used in this work at three different redshifts. The dashed black curves and surrounding red regions correspond to the theoretical mean (Eq. (3)) and 95% C.L. of the halo field conditioned on a region of scale Rnl = 5 Mpc having a density equal to the cosmic average. The solid red curve corresponds to a single realization sampled from these distributions. The sample variance scatter in the red curve is seen to increase toward large masses and high redshifts, as the target mean values become smaller. In the top panel we also show in blue the non-conditional HMF (that is, the limit as R0 → ∞ at δ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 (that is, the difference between the blue and black curves); and (ii) the scatter due to discrete sampling around the target mean CHMF (that is, the difference between the black and red curves). The former determines cosmological signals like 21 cm 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)

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*, given a host halo mass, Mh: p(log M*|log Mh) = 𝒩(log M* | μM*(Mh),σM*). We assume a mass-independent σM* 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:

μ M ∗ ( M h ) = − 1.412 + log M h − log [ ( M h 2.6 × 10 11 M ⊙ ) − 0.5 + ( M h 2.6 × 10 11 M ⊙ ) 0.6 ] . $$ \begin{aligned} \mu _{M_{*}}(M_h)&= -1.412 + \log {M_{\rm h}} \nonumber \\ &\quad - \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{aligned} $$(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 fduty(Mh) = exp[−Mturn/Mh] of halos host star-forming galaxies, taking Mturn = 5 × 108 M⊙ 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 fduty.

In the top panel of Fig. 2 we show our mean SHMR and 2σ 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. 2023), 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.

thumbnail Fig. 2.

Stellar-to-halo mass mass relation and star formation main sequence used in this work. Uppermost panel: Our (redshift independent) stellar to halo mass relation (solid curve) and 2σ scatter (shaded region). Lower panels: Galaxy star-forming main-sequence (solid curves) and 2σ 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 FLARES (Lovell et al. 2021), crosses for ASTRID (Bird et al. 2022; Davies et al. 2023) 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.

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. This scatter is driven primarily by stellar/AGN feedback and mergers. As our fiducial choice of scatter is motivated by the simulations, we are implicitly including these effects; however, our parametric approach can be used to infer the mean and scatter in these relations from data in a simulation-agnostic manner. 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 SFRs of galaxies are known to be strongly correlated with their stellar mass content. The mean of this SFR – M* relation is loosely referred to as the galaxy 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. 2023). 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 timescale is related to the free-fall time at the mean viral density of host halos, tff, which during matter domination scales as the Hubble time: tff ∝ 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 pz(logSFR | log M*) = 𝒩[logSFR | μSFR(M*, z),σSFR(M*)]. 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):

μ SFR ( M ∗ , z ) = log M ∗ − log [ 0.43 H − 1 ( z ) ] . $$ \begin{aligned} \mu _{\rm SFR}(M_{*}, z) = \log {M_{*}} - \log {[0.43 H^{-1}(z)]}. \end{aligned} $$(7)

We assume a mass-dependent scatter that increases toward smaller masses, as these galaxies are expected to be more bursty6:

σ SFR ( M ∗ ) = { − 0.12 log M ∗ + 1.35 if log M ∗ < 10 0.19 . otherwise $$ \begin{aligned} \sigma _{\rm SFR}(M_{*}) = {\left\{ \begin{array}{ll} -0.12 \log {M_{*}} + 1.35&\text{ if}\,\log M_{*}\, < \,10 \\ 0.19.&\text{ otherwise} \end{array}\right.} \end{aligned} $$(8)

The normalization and scaling of the scatter was fit to the hydrodynamic simulations of Ceverino et al. (2018), but the mean (that is, Eq. (7)) is the same as the one in Nikolić et al. (2023). Just like with the SHMR, the scatter is mostly driven by galactic feedback and mergers.

We plot our assumed SFMS and 2σ 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 simulations7.

2.4. Fundamental metallicity relation (FMR)

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, pz(log Z | logSFR, log M*) = 𝒩[log Z | μZ(M*, SFR, z),σZ]. We assume a constant scatter of σZ = 0.1 dex, and a mean given by the following (cf. Curti et al. 2020):

μ Z ( M ∗ , SFR , z ) = 0.296 ( 1 + ( M ∗ M 0 ) − 2.1 ) − 0.148 10 Δ z Z ⊙ , $$ \begin{aligned} \mu _{Z} (M_*, \mathrm{SFR}, z) = 0.296 \left( 1+ \left( \frac{M_{*}}{M_0}\right)^{-2.1} \right)^{-0.148} 10^{\Delta _z} \ Z_{\odot }, \end{aligned} $$(9)

where M0(SFR) ≡ 1010.11 × (SFR/M⊙ yr−1) 0.56 M⊙, and Δz = −0.056z + 0.064 accounting for putative redshift evolution (Curti et al. 2024). In the above, we converted from gas phase to stellar metalicities using Z/Z⊙ = 10(12 + log(O/H) − 8.69) with solar metallicity Z⊙ = 0.02 (Asplund et al. 2004), and adjusting for gas phase metallicities being higher by a factor of ≈2.63 on average (Strom et al. 2018).

In Fig. 3 we show galaxies from a single realization of a comoving volume with radius Rnl = 5 cMpc, at mean density at redshift 6. Each point denotes a single galaxy with the color corresponding to its typical stellar metallicity. It is important to note that the apparent scatter in the metallicity at a fixed M* is considerably larger than the 0.1 dex scatter we set around the mean FMR at a given M* and SFR. This is because the intrinsic scatter in the combination of SHMR and SFMS (that is, the width of pz(SFR, M*) from the previous sections) dominates over our choice of scatter in the metallicity given these properties (that is, the width of pz(Z | SFR, M*); 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. (2024) 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.

thumbnail Fig. 3.

Fundamental metallicity relation used in this work. Each point corresponds to a stellar mass and star formation rate of a galaxy from a single realization of a Rnl = 5 cMpc volume at mean density at z = 6, with the color denoting its metallicity: Z ∼ j $ \tilde{Z}^j $ ∼ p z ( Z | SFR ∼ j , M ∼ * j ) $ {\sim} p_z(Z~|~\tilde{\mathrm{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. (2024) (converted from gas metalicities, see text for details), as well as the metallicity estimate of GN-z11 from Bunker et al. (2023).

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

Soft-band8 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 ∼ 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, LX (in units of erg s−1), given a SFR and metallicity: p(log LX | logSFR, log Z) = 𝒩(log LX | μX(SFR, Z),σX). We assume a constant σX of 0.5 dex and a mean given by:

μ X ( SFR , Z ) = log SFR + 40.5 + log [ ( Z / Z ⊙ 0.05 ) 0.64 + 1 ] . $$ \begin{aligned} \mu _{\rm X}(\mathrm{SFR}, Z) = \log {\mathrm{SFR}}+40.5 + \log {\left[\left(\frac{Z/Z_{\odot }}{0.05}\right)^{0.64}+1\right]}. \end{aligned} $$(10)

These fiducial choices are obtained by assuming a double power-law function for the X-ray luminosity function that results in a flattening at lower metallicities and fits the data at the high SFR/metallicity end (e.g., Fragos et al. 2013; Lehmer et al. 2021; Kaur et al. 2022; Geda et al. 2024). They are roughly consistent with empirical fits to local galaxies, which could however suffer from incompleteness at the lowest SFR bins (e.g., Brorby et al. 2016). We have 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 Γ = 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 redshift-binned stacked observations of star-forming galaxies at z ∼ 2 from Fornasini et al. (2019) and Lyman-Break analogs 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).

thumbnail Fig. 4.

Luminosity scalings used in this work. 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σ range for metalicities Z = −3.6, −2.0, respectively. Green stars correspond to values from local star forming galaxies, discussed in Brorby et al. (2016). Purple diamonds are redshift-binned stacks from Fornasini et al. (2019) for z = 1.5 and 2.3. 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σ scatter around the mean relation. Both panels assume that metallicity follows the mean FMR.

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, LUV, as a function of SFR, Z, and SFR history. For the latter we assume that our sampled SFR is exponentially declining toward higher redshifts, as implied by Equation (7). Therefore LUV is sampled assuming a log-normal conditional probability p(log LUV | logSFR, log Z) = 𝒩(log LUV | μUV, BPASS(SFR, Z),σL) where μUV, BPASS is the predicted luminosity from BPASS. We add an additional scatter of σL = 0.1 dex around the mean to compensate for unaccounted sources of stochasticity, for example, 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 LUV with 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Å−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. 2023; 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): f ¯ esc = 0.053 $ \bar{f}_{\mathrm{esc}} = 0.053 $ for the mean. We do not assume that the mean or scatter of the escape fraction depend on galaxy properties, as such dependencies are not yet well established for unbiased galaxy samples, especially at high redshifts. For example, assuming that scatter depends on mass could accelerate or decelerate reionization, by effectively shifting the population-averaged mean as a function of redshift. 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). We defer such studies to future work.

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 (that is, using only the corresponding mean relation with zero scatter), normalized to the values of the full distribution containing all sources of scatter: μ ϵ / μ ϵ full $ \mu_{\epsilon}/\mu_{\epsilon}^{\mathrm{full}} $ and σ ϵ / σ ϵ full $ \sigma_{\epsilon}/\sigma_{\epsilon}^{\mathrm{full}} $. As mentioned above, we consider the following sources of scatter:

  • (i)

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

  • (ii)

    Poisson sample variance in halo number around the target mean CHMF

  • (iii)

    scatter around the SHMR

  • (iv)

    scatter around the SFMS

  • (v)

    scatter around the FMR

  • (vi)

    scatter in the mapping of the intrinsic luminosity to SFR, M* and Z

  • (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 Rnl = 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 21 cm 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 ∼ 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 10 000 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.

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 Rnl = 5 cMpc realizations that have a zero emissivity is denoted at the bottom of each violin.

thumbnail Fig. 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 ∼ 7, consistent with the latest estimates of the EoR history e.g., Qin et al. in prep.

As galaxies become rarer toward 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. This is shown in Fig. 6 with a rectangle at the bottom of the figure representing the probability that a region of Rnl = 5 cMpc has zero emissivity. At z ≳ 18 the majority of Rnl = 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(fesc). If instead we assume a binomial p(fesc) distribution, the majority of Rnl = 5 cMpc have zero ionizing emissivity already by z ≳ 15.

thumbnail Fig. 6.

Fractional contribution of different sources of scatter (see legend) to the mean and standard deviation of the ionizing emissivity, computed over Rnl = 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.

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 (that is, only using the corresponding mean relation with no scatter). In Figure 6 we plot the corresponding mean (top panel) and standard deviation9 (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 ∼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 ∼ 5. This underprediction in the std rises to ∼50% toward z ∼ 20, as the typical galaxies have smaller stellar masses and therefore a broader p(SFR | M*) (cf. 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 ∼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

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 Rnl = 5 cMpc, while the rectangles represent probabilities that a region of size Rnl = 5 cMpc has zero X-ray emissivity. As expected the means and widths of the distributions decrease toward higher redshifts.

thumbnail Fig. 7.

Same as Fig. 5, but for soft band (0.5–2 keV) X-ray emissivity.

We see that the X-ray emissivities have broader distributions compared with the ionizing emissivities in Fig. 5. For example, at z ∼ 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 | SFR, M*, 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 subdominant compared with the wider X-ray intrinsic 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 μi/μfull and σi/σfull in the upper and lower panels, respectively.

thumbnail Fig. 8.

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

The biggest impact on the mean and std comes from the scatter in the SFR−M* 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 60% 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 LX−SFR relation (violet dash-dotted curves). The relative difference in standard deviations is roughly 30% at z = 20, rising to 45% at z = 5. At z = 5 it is more important than the scatter in SFR−M* 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 LX−SFR scatter at earlier times.

Scatter in the SHMR has a ∼10% effect, again without redshift dependence since we chose a constant width for p(M* | Mh). Scatter in the other terms has a negligible impact on the X-ray emissivity.

3.3. Lyman Werner 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 21 cm signal). For concreteness, here we evaluate the emissivity in the Lyman-Werner band (11.2–13.6eV) 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 (that is, 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.

thumbnail Fig. 9.

Same as Fig. 5, but for the Lyman-Werner (11.2–13.6 eV) 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 ∼ 20 the burstiness of SFR contributes at a ∼50% level to the std of the distribution, but this drops to ∼10% at z ∼ 5.

thumbnail Fig. 10.

Same as Fig. 6, but for the Lyman Werner (11.2–13.6 eV) emissivity.

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

4. Results: EoR history

The ionizing emissivities shown in the previous section can be used to estimate the redshift evolution of the volume filling factor of ionized regions, QHII(z) – the EoR history. Even though it is an average quantity, computing the EoR history accurately requires accounting for the spatial and temporal coevolution 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 (e.g., Madau et al. 1999), obtained by:

d Q HII d t = n ˙ ion / b − α A C ⟨ n H ⟩ Q HII . $$ \begin{aligned} \frac{\mathrm{d}Q_{\rm HII}}{\mathrm{d}t} =\dot{n}_{\rm ion/b} - \alpha _{\rm A} \ C \ \langle n_{H} \rangle \ Q_{\rm HII}. \end{aligned} $$(11)

Here ṅion/b is the ionizing emissivity per baryon predicted by our model10, αA is the case-A recombination coefficient, ⟨nH⟩ is the mean hydrogen density, and C ≡ ⟨nH2⟩/⟨nH⟩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 photons11. Estimates of the EoR history obtained with Eq. (11) underpredict the duration of the EoR by Δz ∼ 1–2, with the error increasing toward 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 Δz ∼ 0.1. Ignoring scatter around the SFMS has a bigger impact, delaying the EoR history by Δz ∼ 0.5–1. Ignoring scatter in all galaxy properties underestimates the duration of the EoR and delays the end stages by up to Δz ∼ 2.

thumbnail Fig. 11.

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 (cf. 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 (cf. 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 (that is, scatter around the SFMS) can result in EoR histories that are delayed by Δz ∼ 0.5–1. Not accounting for any scatter and assuming only mean galaxy properties delays the completion of the EoR by Δz ∼ 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, that is, ṅion/b>tage > 2 (shown by the dashed line in Fig. 5; cf. 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 Δz ∼ 0.5–1. Neglecting all galaxy-to-galaxy scatter by assuming mean values (cf. Eq. (1)) results in an overextended EoR history and delays its completion by Δz ∼ 1–2.

Our results suggest that inferring galaxy properties from EoR history data without accounting for stochasticity could bias recovery toward 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 Å to 1550 Å 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.

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. 2024 (R24); Harikane et al. 2024 (H23); McLeod et al. 2024 (M24); Willott et al. 2024 (W24)).

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

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 MUV ∼ −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 ≲ 10. However, the mean underpredicts the recent estimates of UV LFs at z ≳ 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 ≳ 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.

We also calculate the star formation density evolution from galaxies down to MUV < −17. In order to provide a more like-to-like comparison with observations, we do not use the sampled SFRs directly but instead convert from the UV luminosity via a constant conversion factor: SFR(M⊙ yr−1) = KUVLUV(erg s−1 Hz−1), where KUV is the conversion factor which depends on the IMF and star formation history. We take KUV = 1.15 × 10−28 M⊙ yr−1/erg s−1 Hz−1 (Sun & Furlanetto 2016), consistent with other works. The result is shown in Fig. 13. The solid curve represents our fiducial model which includes all sources of stochasticity, while the dashed curve represents the model where only the mean relations are considered. We see that our model reproduces the data very well at z ≲ 12 (Bouwens et al. 2020), though is somewhat lower than some estimates at higher redshifts. As this statistic is fundamentally an integral over the UV LFs, we reach the same qualitative conclusions. Quantitatively, the discrepancy with the data seems less than for the UV LFs. Due to the steepness of the UV LFs, the integrated SFRD is dominated by the faint end limit used to compute it, while the z > 12 JWST observations are more discrepant on the bright end.

thumbnail Fig. 13.

Star formation rate density from bright galaxies (MUV < −17) as a function of density. Our fiducial model including all of the aforementioned sources of scatter is shown with the solid line and surrounding shaded regions (68% C.L.). The dashed line corresponds to the SFRD calculated using only mean relations without any scatter. Also shown are high-redshift estimates from Bouwens et al. (2020), Harikane et al. (2023), Willott et al. (2024).

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

For this work we quantified how the galaxy-to-galaxy scatter in their properties impacts estimates of their emissivities and related observables. We used a semi-empirical model that explicitly defines scatter around well-studied mean relations: (i) the conditional halo mass function (CHMF); (ii) the SHMR; (iii) th galaxy SFMS; (iv) the FMR; (v) the conditional intrinsic luminosity; and (vi) the photon escape fraction. We computed the corresponding multi-frequency (ionizing UV, X-rays, and LW) emissivities, EoR histories, and UV LFs, quantifying the relative importance of the above sources of scatter.

We find that the burstiness of star formation (that is, scatter around the mean SFMS) is important for all emissivities. Because we assume burstiness increases toward 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 a few tens 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 on the order of ∼10%. For the X-ray emissivity, one must account for scatter in the intrinsic luminosity, which in our fiducial model is driven by HMXB luminosity functions.

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

We recovered 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 MUV ∼ −18. Our UV LFs are consistent with observational data at z ≤ 10 but consistently under-predict recent estimates 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 21 cm signal) additionally requires accounting for scatter in the intrinsic X-ray luminosities of galaxies.

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

Data availability

The code related to the work is publicly available at https://github.com/IvanNikolic21/Stochasticitysampler.

1

Throughout we use the subscript “nl” to indicate nonlinear (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.

2

We use “z” subscripts to indicate probability distributions that are also functions of redshift (see below for more details).

3

Note that our choice of conditional distributions is motivated by well-known relations, but this choice is not unique. Furthermore, we do not assume any direct correlation of the ionizing escape fraction to other properties, as there is currently no consensus on what galaxy property would constitute an appropriate basis to characterize the distribution.

4

Strictly speaking, the average HMF is not the mean of the conditional halo mass function (CHMF), but its limit as R0 → ∞ and δ0 → 0.

5

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

6

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

7

Detailed 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 ∼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.

8

Here we define the soft band to be 0.5 − 2 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).

9

For numerical stability, we calculate the standard deviation by fitting a log-normal to the nonzero distribution of emissivities.

10

We make the standard assumption that helium is singly ionized by stellar sources together with hydrogen, due to their comparable ionization thresholds.

11

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

Acknowledgments

We thank the anonymous referee for their insightful comments. 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

  1. Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 924, 51 [NASA ADS] [CrossRef] [Google Scholar]
  2. Allende Prieto, C., Koesterke, L., Hubeny, I., et al. 2018, A&A, 618, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Asplund, M., Grevesse, N., Sauval, A. J., Allende Prieto, C., & Kiselman, D. 2004, A&A, 417, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Barkana, R., & Loeb, A. 2004, ApJ, 609, 474 [NASA ADS] [CrossRef] [Google Scholar]
  5. Barrow, K. S. S., Wise, J. H., Norman, M. L., O’Shea, B. W., & Xu, H. 2017, MNRAS, 469, 4863 [Google Scholar]
  6. Barrow, K. S. S., Robertson, B. E., Ellis, R. S., et al. 2020, ApJ, 902, L39 [NASA ADS] [CrossRef] [Google Scholar]
  7. Basu-Zych, A. R., Lehmer, B. D., Hornschemeier, A. E., et al. 2013, ApJ, 774, 152 [NASA ADS] [CrossRef] [Google Scholar]
  8. Becker, G. D., D’Aloisio, A., Christenson, H. M., et al. 2021, MNRAS, 508, 1853 [NASA ADS] [CrossRef] [Google Scholar]
  9. Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bird, S., Ni, Y., Di Matteo, T., et al. 2022, MNRAS, 512, 3703 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 382, 325 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34 [Google Scholar]
  13. Bouwens, R. J., Oesch, P. A., Labbé, I., et al. 2016, ApJ, 830, 67 [Google Scholar]
  14. Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47 [NASA ADS] [CrossRef] [Google Scholar]
  16. Bouwens, R. J., Stefanon, M., Brammer, G., et al. 2023, MNRAS, 523, 1036 [NASA ADS] [CrossRef] [Google Scholar]
  17. Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 [Google Scholar]
  18. Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151 [Google Scholar]
  19. Brorby, M., Kaaret, P., Prestwich, A., & Mirabel, I. F. 2016, MNRAS, 457, 4081 [NASA ADS] [CrossRef] [Google Scholar]
  20. Bunker, A. J., Saxena, A., Cameron, A. J., et al. 2023, A&A, 677, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Byrne, C. M., & Stanway, E. R. 2023, MNRAS, 521, 4995 [NASA ADS] [CrossRef] [Google Scholar]
  22. Byrne, C. M., Stanway, E. R., Eldridge, J. J., McSwiney, L., & Townsend, O. T. 2022, MNRAS, 512, 5329 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cen, R., & Kimm, T. 2015, ApJ, 801, L25 [Google Scholar]
  24. Ceverino, D., Klessen, R. S., & Glover, S. C. O. 2018, MNRAS, 480, 4842 [NASA ADS] [Google Scholar]
  25. Ciardi, B., Ferrara, A., & White, S. D. M. 2003, MNRAS, 344, L7 [NASA ADS] [CrossRef] [Google Scholar]
  26. Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944 [Google Scholar]
  27. Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024, A&A, 684, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Curtis-Lake, E., Chevallard, J., Charlot, S., & Sandles, L. 2021, MNRAS, 503, 4855 [NASA ADS] [CrossRef] [Google Scholar]
  29. Das, A., Mesinger, A., Pallottini, A., Ferrara, A., & Wise, J. H. 2017, MNRAS, 469, 1166 [NASA ADS] [CrossRef] [Google Scholar]
  30. Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827 [Google Scholar]
  31. Davies, F. B., & Furlanetto, S. R. 2016, MNRAS, 460, 1328 [NASA ADS] [CrossRef] [Google Scholar]
  32. Davies, F. B., & Furlanetto, S. R. 2022, MNRAS, 514, 1302 [NASA ADS] [CrossRef] [Google Scholar]
  33. Davies, J. E., Bird, S., Mutch, S., et al. 2023, MNRAS, 525, 2553 [NASA ADS] [CrossRef] [Google Scholar]
  34. Davies, F. B., Bosman, S. E. I., Gaikwad, P., et al. 2024, ApJ, 965, 134 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dayal, P., & Ferrara, A. 2018, Phys. Rep., 780, 1 [Google Scholar]
  36. Di Cesare, C., Graziani, L., Schneider, R., et al. 2023, MNRAS, 519, 4632 [Google Scholar]
  37. Dixon, K. L., Iliev, I. T., Mellema, G., Ahn, K., & Shapiro, P. R. 2016, MNRAS, 456, 3011 [NASA ADS] [CrossRef] [Google Scholar]
  38. Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011 [Google Scholar]
  39. Donnan, C. T., McLure, R. J., Dunlop, J. S., et al. 2024, MNRAS, 533, 3222 [NASA ADS] [CrossRef] [Google Scholar]
  40. Eide, M. B., Graziani, L., Ciardi, B., et al. 2018, MNRAS, 476, 1174 [NASA ADS] [CrossRef] [Google Scholar]
  41. Ferrara, A., Pallottini, A., & Dayal, P. 2023, MNRAS, 522, 3986 [NASA ADS] [CrossRef] [Google Scholar]
  42. Finkelstein, S. L., Bagley, M. B., Ferguson, H. C., et al. 2023, ApJ, 946, L13 [NASA ADS] [CrossRef] [Google Scholar]
  43. Fornasini, F. M., Kriek, M., Sanders, R. L., et al. 2019, ApJ, 885, 65 [NASA ADS] [CrossRef] [Google Scholar]
  44. Fragos, T., Lehmer, B. D., Naoz, S., Zezas, A., & Basu-Zych, A. 2013, ApJ, 776, L31 [CrossRef] [Google Scholar]
  45. Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45 [NASA ADS] [CrossRef] [Google Scholar]
  46. Furlanetto, S. R. 2006, MNRAS, 371, 867 [NASA ADS] [CrossRef] [Google Scholar]
  47. Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gaikwad, P., Haehnelt, M. G., Davies, F. B., et al. 2023, MNRAS, 525, 4093 [NASA ADS] [CrossRef] [Google Scholar]
  49. Garcia, A. M., Torrey, P., Ellison, S., et al. 2024, MNRAS, 531, 1398 [NASA ADS] [CrossRef] [Google Scholar]
  50. Geda, R., Goulding, A. D., Lehmer, B. D., Greene, J. E., & Kulkarni, A. 2024, ApJ, 965, 67 [NASA ADS] [CrossRef] [Google Scholar]
  51. Gelli, V., Mason, C., & Hayward, C. C. 2024, ApJ, submitted [arXiv:2405.13108] [Google Scholar]
  52. Gillet, N. J. F., Mesinger, A., & Park, J. 2020, MNRAS, 491, 1980 [NASA ADS] [Google Scholar]
  53. Grazian, A., Giallongo, E., Paris, D., et al. 2017, A&A, 602, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Haiman, Z., Abel, T., & Rees, M. J. 2000, ApJ, 534, 11 [NASA ADS] [CrossRef] [Google Scholar]
  55. Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123 [NASA ADS] [CrossRef] [Google Scholar]
  56. Harikane, Y., Ouchi, M., Oguri, M., et al. 2023, ApJS, 265, 5 [NASA ADS] [CrossRef] [Google Scholar]
  57. Harikane, Y., Nakajima, K., Ouchi, M., et al. 2024, ApJ, 960, 56 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hassan, S., Davé, R., McQuinn, M., et al. 2022, ApJ, 931, 62 [NASA ADS] [CrossRef] [Google Scholar]
  59. Holzbauer, L. N., & Furlanetto, S. R. 2012, MNRAS, 419, 718 [NASA ADS] [CrossRef] [Google Scholar]
  60. Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73 [NASA ADS] [CrossRef] [Google Scholar]
  61. Izotov, Y. I., Schaerer, D., Thuan, T. X., et al. 2016, MNRAS, 461, 3683 [Google Scholar]
  62. Kannan, R., Garaldi, E., Smith, A., et al. 2022, MNRAS, 511, 4005 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kauffmann, O. B., Ilbert, O., Weaver, J. R., et al. 2022, A&A, 667, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Kaur, H. D., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 513, 5097 [NASA ADS] [CrossRef] [Google Scholar]
  65. Kaurov, A. A., Venumadhav, T., Dai, L., & Zaldarriaga, M. 2018, ApJ, 864, L15 [NASA ADS] [CrossRef] [Google Scholar]
  66. Koopmans, L., Pritchard, J., Mellema, G., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 1 [Google Scholar]
  67. Kostyuk, I., Nelson, D., Ciardi, B., Glatzle, M., & Pillepich, A. 2023, MNRAS, 521, 3077 [CrossRef] [Google Scholar]
  68. Kravtsov, A., & Belokurov, V. 2024, ArXiv e-prints [arXiv:2405.04578] [Google Scholar]
  69. Kreilgaard, K. C., Mason, C. A., Cullen, F., Begley, R., & McLure, R. J. 2024, A&A, 692, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Leethochawalit, N., Roberts-Borsani, G., Morishita, T., Trenti, M., & Treu, T. 2023, MNRAS, 524, 5454 [NASA ADS] [CrossRef] [Google Scholar]
  71. Lehmer, B. D., Basu-Zych, A. R., Mineo, S., et al. 2016, ApJ, 825, 7 [Google Scholar]
  72. Lehmer, B. D., Eufrasio, R. T., Basu-Zych, A., et al. 2021, ApJ, 907, 17 [NASA ADS] [CrossRef] [Google Scholar]
  73. Lin, Y., Oh, S. P., Furlanetto, S. R., & Sutter, P. M. 2016, MNRAS, 461, 3361 [NASA ADS] [CrossRef] [Google Scholar]
  74. Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113 [Google Scholar]
  75. Lovell, C. C., Vijayan, A. P., Thomas, P. A., et al. 2021, MNRAS, 500, 2127 [Google Scholar]
  76. Ma, X., Hopkins, P. F., Garrison-Kimmel, S., et al. 2018, MNRAS, 478, 1694 [NASA ADS] [CrossRef] [Google Scholar]
  77. Madau, P., & Fragos, T. 2017, ApJ, 840, 39 [Google Scholar]
  78. Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648 [CrossRef] [Google Scholar]
  79. Magg, M., Reis, I., Fialkov, A., et al. 2022, MNRAS, 514, 4433 [NASA ADS] [CrossRef] [Google Scholar]
  80. Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115 [NASA ADS] [CrossRef] [Google Scholar]
  81. Mascia, S., Pentericci, L., Calabrò, A., et al. 2023, A&A, 672, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Mason, C. A., Trenti, M., & Treu, T. 2023, MNRAS, 521, 497 [NASA ADS] [CrossRef] [Google Scholar]
  83. McLeod, D. J., Donnan, C. T., McLure, R. J., et al. 2024, MNRAS, 527, 5004 [Google Scholar]
  84. McQuinn, M. 2016, ARA&A, 54, 313 [NASA ADS] [CrossRef] [Google Scholar]
  85. McQuinn, M., Lidz, A., Zahn, O., et al. 2007, MNRAS, 377, 1043 [Google Scholar]
  86. Meriot, R., & Semelin, B. 2024, A&A, 683, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Mesinger, A. 2016, Astrophys. Space Sci. Libr., 423 [CrossRef] [Google Scholar]
  88. Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955 [Google Scholar]
  89. Mesinger, A., Ferrara, A., & Spiegel, D. S. 2013, MNRAS, 431, 621 [NASA ADS] [CrossRef] [Google Scholar]
  90. Mineo, S., Gilfanov, M., & Sunyaev, R. 2012, MNRAS, 419, 2095 [Google Scholar]
  91. Mirocha, J., & Furlanetto, S. R. 2023, MNRAS, 519, 843 [Google Scholar]
  92. Mirocha, J., Furlanetto, S. R., & Sun, G. 2017, MNRAS, 464, 1365 [NASA ADS] [CrossRef] [Google Scholar]
  93. Mirocha, J., La Plante, P., & Liu, A. 2021, MNRAS, 507, 3872 [NASA ADS] [CrossRef] [Google Scholar]
  94. Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347 [Google Scholar]
  95. Muñoz, J. B., Qin, Y., Mesinger, A., et al. 2022, MNRAS, 511, 3657 [CrossRef] [Google Scholar]
  96. Murmu, C. S., Datta, K. K., Majumdar, S., & Greve, T. R. 2024, JCAP, 2024, 032 [Google Scholar]
  97. Mutch, S. J., Geil, P. M., Poole, G. B., et al. 2016, MNRAS, 462, 250 [NASA ADS] [CrossRef] [Google Scholar]
  98. Naidu, R. P., Oesch, P. A., van Dokkum, P., et al. 2022, ApJ, 940, L14 [NASA ADS] [CrossRef] [Google Scholar]
  99. Nikolić, I., Mesinger, A., Qin, Y., & Gorce, A. 2023, MNRAS, 526, 3170 [CrossRef] [Google Scholar]
  100. Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129 [NASA ADS] [CrossRef] [Google Scholar]
  101. Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105 [Google Scholar]
  102. Oh, S. P. 2001, ApJ, 553, 499 [NASA ADS] [CrossRef] [Google Scholar]
  103. Pacucci, F., Mesinger, A., Mineo, S., & Ferrara, A. 2014, MNRAS, 443, 678 [CrossRef] [Google Scholar]
  104. Pahl, A. J., Shapley, A., Steidel, C. C., et al. 2023, MNRAS, 521, 3247 [NASA ADS] [CrossRef] [Google Scholar]
  105. Pallottini, A., & Ferrara, A. 2023, A&A, 677, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Pallottini, A., Ferrara, A., Gallerani, S., et al. 2022, MNRAS, 513, 5621 [NASA ADS] [Google Scholar]
  107. Park, J., Mesinger, A., Greig, B., & Gillet, N. 2019, MNRAS, 484, 933 [NASA ADS] [CrossRef] [Google Scholar]
  108. Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJ, 951, L1 [CrossRef] [Google Scholar]
  109. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526 [Google Scholar]
  111. Prelogović, D., Mesinger, A., Murray, S., Fiameni, G., & Gillet, N. 2022, MNRAS, 509, 3852 [Google Scholar]
  112. Press, W. H., & Schechter, P. 1974, ApJ, 187, 425 [Google Scholar]
  113. Qin, Y., Mesinger, A., Bosman, S. E. I., & Viel, M. 2021, MNRAS, 506, 2390 [NASA ADS] [CrossRef] [Google Scholar]
  114. Reis, I., Barkana, R., & Fialkov, A. 2022, MNRAS, 511, 5265 [NASA ADS] [CrossRef] [Google Scholar]
  115. Robertson, B., Johnson, B. D., Tacchella, S., et al. 2024, ApJ, 970, 31 [NASA ADS] [CrossRef] [Google Scholar]
  116. Ross, H. E., Dixon, K. L., Iliev, I. T., & Mellema, G. 2017, MNRAS, 468, 3785 [NASA ADS] [CrossRef] [Google Scholar]
  117. Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76 [NASA ADS] [CrossRef] [Google Scholar]
  118. Schaeffer, T., Giri, S. K., & Schneider, A. 2023, MNRAS, 526, 2942 [NASA ADS] [CrossRef] [Google Scholar]
  119. Shen, X., Vogelsberger, M., Boylan-Kolchin, M., Tacchella, S., & Kannan, R. 2023, MNRAS, 525, 3254 [NASA ADS] [CrossRef] [Google Scholar]
  120. Sheth, R. K. 1998, MNRAS, 300, 1057 [CrossRef] [Google Scholar]
  121. Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119 [Google Scholar]
  122. Sobacchi, E., & Mesinger, A. 2013, MNRAS, 432, L51 [NASA ADS] [CrossRef] [Google Scholar]
  123. Sobacchi, E., & Mesinger, A. 2014, MNRAS, 440, 1662 [NASA ADS] [CrossRef] [Google Scholar]
  124. Stanway, E. R., & Eldridge, J. J. 2018, MNRAS, 479, 75 [NASA ADS] [CrossRef] [Google Scholar]
  125. Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29 [NASA ADS] [CrossRef] [Google Scholar]
  126. Steidel, C. C., Bogosavljević, M., Shapley, A. E., et al. 2018, ApJ, 869, 123 [Google Scholar]
  127. Strom, A. L., Steidel, C. C., Rudie, G. C., Trainor, R. F., & Pettini, M. 2018, ApJ, 868, 117 [NASA ADS] [CrossRef] [Google Scholar]
  128. Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417 [CrossRef] [Google Scholar]
  129. Sun, G., Faucher-Giguère, C.-A., Hayward, C. C., et al. 2023, ApJ, 955, L35 [CrossRef] [Google Scholar]
  130. Tacchella, S., Dekel, A., Carollo, C. M., et al. 2016, MNRAS, 457, 2790 [Google Scholar]
  131. Trapp, A. C., & Furlanetto, S. R. 2020, MNRAS, 499, 2401 [CrossRef] [Google Scholar]
  132. Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110 [NASA ADS] [CrossRef] [Google Scholar]
  133. Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435 [NASA ADS] [CrossRef] [Google Scholar]
  134. Willott, C. J., Desprez, G., Asada, Y., et al. 2024, ApJ, 966, 74 [NASA ADS] [CrossRef] [Google Scholar]
  135. Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2011, MNRAS, 418, 838 [NASA ADS] [CrossRef] [Google Scholar]
  136. Xu, H., Ahn, K., Wise, J. H., Norman, M. L., & O’Shea, B. W. 2014, ApJ, 791, 110 [NASA ADS] [CrossRef] [Google Scholar]
  137. Xu, H., Wise, J. H., Norman, M. L., Ahn, K., & O’Shea, B. W. 2016, ApJ, 833, 84 [NASA ADS] [CrossRef] [Google Scholar]
  138. Yeh, J. Y. C., Smith, A., Kannan, R., et al. 2023, MNRAS, 520, 2757 [NASA ADS] [CrossRef] [Google Scholar]
  139. Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019, MNRAS, 490, 2855 [NASA ADS] [CrossRef] [Google Scholar]
  140. Zaroubi, S. 2013, Astrophys. Space Sci. Libr., 396, 45 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Convergence of the mean emissivity with scale

Here we confirm that our fiducial choice of Rnl = 5 cMpc when computing emissivity distributions converges to the correct mean. In principle, the averaging over the scale-dependent overdensity distribution, pz(δ0|Rnl), in the top row of Eq. 2 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 A.1 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.

thumbnail Fig. A.1.

The mean of the ionizing emissivity computed over regions of varying scales, Rnl = 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 Rnl = 1 cMpc by factors of few – 10.

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 such as 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 Rnl= 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:

ε ¯ = ∫ d M h d n ( M h , z ) d M h L ( M h ) , $$ \begin{aligned} \overline{\varepsilon } = \int \mathrm{d}M_{h} \frac{\mathrm{d}n(M_{h},z)}{\mathrm{d}M_{h}} L (M_{h}), \end{aligned} $$(B.1)

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. 2. 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. 2. For simplicity, in this section we reduce the dependencies of luminosity to only L(Mh) with an appropriate log-normal PDF p(log L|Mh) = 𝒩(μL, σL) without loss of generality. The Eq. 2 becomes:

ε ¯ = ∫ d M h d n ( M h ) d M h ∫ d L L p ( log 10 L | log 10 M h ) = ∫ d M h d n ( M h ) d M h ∫ d L L 1 2 π log e 10 σ L L exp ( − ( log 10 L − μ L ) 2 2 σ L 2 ) . $$ \begin{aligned} \overline{\varepsilon } =&\int \mathrm{d}M_{h} \frac{\mathrm{d}n(M_{h})}{\mathrm{d}M_{h}} \int \mathrm{d}L \ L \ p(\log _{10} L | \log _{10} M_{h}) \nonumber \\ =&\int \mathrm{d}M_{h} \frac{\mathrm{d}n(M_{h})}{\mathrm{d}M_{h}} \int \mathrm{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{aligned} $$(B.2)

The integral on the right can be analytically computed:

ε ¯ = ∫ d M h d n ( M h ) d M h 10 ( μ L + log e 10 σ L 2 2 ) . $$ \begin{aligned} \overline{\varepsilon } = \int \mathrm{d}M_h \frac{\mathrm{d}n(M_{h})}{\mathrm{d}M_{h}} 10^{\left(\mu _L + \frac{\log _{e}{10} \ \sigma _L^2}{2}\right)}. \end{aligned} $$(B.3)

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. B.3. The factor 10μL corresponds to the mean of the L(Mh) scaling relation, that is, what one would obtain if one did not consider scatter around the mean. The second part, 10 log e 10 σ L 2 2 $ 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 (for example, one choice for the escape fraction in Sec. 2.6) then the equation becomes:

ε ¯ = ∫ d M h d n ( M h ) d M h ∫ d L L p ( log 10 L | log 10 M h ) × ∫ d f esc f esc ( n k ) P n ( 1 − P ) k , $$ \begin{aligned} \overline{\varepsilon }&= \int \mathrm{d}M_{h} \frac{\mathrm{d}n(M_{h})}{\mathrm{d}M_{h}} \int \mathrm{d}L \ L \ p(\log _{10} L | \log _{10} M_{h}) \nonumber \\&\qquad \times \int \mathrm{d}f_{\rm esc} f_{\rm esc}{{n}\atopwithdelims (){k}} \mathcal{P} ^n(1-\mathcal{P} )^k, \end{aligned} $$(B.4)

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

ε ¯ = ∫ d M h d n ( M h ) d M h ∫ d L L p ( log 10 L | log 10 M h ) f ¯ esc . $$ \begin{aligned} \overline{\varepsilon } =&\int \mathrm{d}M_{h} \frac{\mathrm{d}n(M_{h})}{\mathrm{d}M_{h}} \int \mathrm{d}L \ L \ p(\log _{10} L | \log _{10} M_{h}) \overline{f}_{\rm esc}. \end{aligned} $$(B.5)

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

All Figures

thumbnail Fig. 1.

Example halo mass functions used in this work at three different redshifts. The dashed black curves and surrounding red regions correspond to the theoretical mean (Eq. (3)) and 95% C.L. of the halo field conditioned on a region of scale Rnl = 5 Mpc having a density equal to the cosmic average. The solid red curve corresponds to a single realization sampled from these distributions. The sample variance scatter in the red curve is seen to increase toward large masses and high redshifts, as the target mean values become smaller. In the top panel we also show in blue the non-conditional HMF (that is, the limit as R0 → ∞ at δ0 = 0).
In the text
thumbnail Fig. 2.

Stellar-to-halo mass mass relation and star formation main sequence used in this work. Uppermost panel: Our (redshift independent) stellar to halo mass relation (solid curve) and 2σ scatter (shaded region). Lower panels: Galaxy star-forming main-sequence (solid curves) and 2σ 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 FLARES (Lovell et al. 2021), crosses for ASTRID (Bird et al. 2022; Davies et al. 2023) 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.
In the text
thumbnail Fig. 3.

Fundamental metallicity relation used in this work. Each point corresponds to a stellar mass and star formation rate of a galaxy from a single realization of a Rnl = 5 cMpc volume at mean density at z = 6, with the color denoting its metallicity: Z ∼ j $ \tilde{Z}^j $ ∼ p z ( Z | SFR ∼ j , M ∼ * j ) $ {\sim} p_z(Z~|~\tilde{\mathrm{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. (2024) (converted from gas metalicities, see text for details), as well as the metallicity estimate of GN-z11 from Bunker et al. (2023).
In the text
thumbnail Fig. 4.

Luminosity scalings used in this work. 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σ range for metalicities Z = −3.6, −2.0, respectively. Green stars correspond to values from local star forming galaxies, discussed in Brorby et al. (2016). Purple diamonds are redshift-binned stacks from Fornasini et al. (2019) for z = 1.5 and 2.3. 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σ scatter around the mean relation. Both panels assume that metallicity follows the mean FMR.
In the text
thumbnail Fig. 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 ∼ 7, consistent with the latest estimates of the EoR history e.g., Qin et al. in prep.
In the text
thumbnail Fig. 6.

Fractional contribution of different sources of scatter (see legend) to the mean and standard deviation of the ionizing emissivity, computed over Rnl = 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.
In the text
thumbnail Fig. 7.

Same as Fig. 5, but for soft band (0.5–2 keV) X-ray emissivity.
In the text
thumbnail Fig. 8.

Like Fig. 6, but for the soft band (0.5–2 keV) X-ray emissivity.
In the text
thumbnail Fig. 9.

Same as Fig. 5, but for the Lyman-Werner (11.2–13.6 eV) emissivity.
In the text
thumbnail Fig. 10.

Same as Fig. 6, but for the Lyman Werner (11.2–13.6 eV) emissivity.
In the text
thumbnail Fig. 11.

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 (cf. 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 (cf. 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 (that is, scatter around the SFMS) can result in EoR histories that are delayed by Δz ∼ 0.5–1. Not accounting for any scatter and assuming only mean galaxy properties delays the completion of the EoR by Δz ∼ 2.
In the text
thumbnail Fig. 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 the text
thumbnail Fig. 13.

Star formation rate density from bright galaxies (MUV < −17) as a function of density. Our fiducial model including all of the aforementioned sources of scatter is shown with the solid line and surrounding shaded regions (68% C.L.). The dashed line corresponds to the SFRD calculated using only mean relations without any scatter. Also shown are high-redshift estimates from Bouwens et al. (2020), Harikane et al. (2023), Willott et al. (2024).
In the text
thumbnail Fig. A.1.

The mean of the ionizing emissivity computed over regions of varying scales, Rnl = 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 Rnl = 1 cMpc by factors of few – 10.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.