1. Introduction
Despite recent advances in our understanding of the abundances and properties of high-redshift galaxies (e.g. Eisenstein et al. Reference Eisenstein2023; Adams et al. Reference Adams2024; Finkelstein et al. Reference Finkelstein2024; Donnan et al. Reference Donnan2024; Harikane et al. Reference Harikane2024), their ionising properties remain highly uncertain. Recent efforts have focused on directly measuring
$\xi_{\textrm{ion}}$
using JWST (Simmonds et al. Reference Simmonds2023; Atek et al. Reference Atek2024; Pahl et al. Reference Pahl2024; Simmonds et al. Reference Simmonds2024), and studying
$f_{\textrm{esc}}$
in low-redshift analogs of reionisation-era galaxies. One goal of the latter is to discover correlations with other observables that can be used to infer
$f_{\textrm{esc}}$
at higher redshifts (Pahl et al. Reference Pahl, Shapley, Steidel, Chen and Reddy2021, Chisholm et al. Reference Chisholm2022, Flury et al. Reference Flury2022, Jaskot et al. Reference Jaskot2024b, Reference Jaskot2024a). However, a complementary approach is to measure the collective ionising output of galaxies independently by leveraging constraints on the ionisation state of the IGM, and then infer the required ionising properties of galaxies (e.g. Becker & Bolton Reference Becker and Bolton2013; Becker et al. Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021; Gaikwad et al. Reference Gaikwad2023; Bosman & Davies Reference Bosman and Davies2024).
The net ionising photon emissivity of the galaxy population is given by
$\begin{equation} \dot{N}_{\textrm{ion}} = \rho_{\textrm{UV}} \langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}\end{equation}$
where
$\rho_{\textrm{UV}}$
is the integrated UV luminosity density. The quantity
$\langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
is the UV-luminosity (
$L_{\textrm{UV}}$
)-weighted product of two quantities that together quantify the ionising properties of galaxies. The first,
$\xi_{\textrm{ion}}$
, is the ‘intrinsic ionising efficiency’, or the rate of HI-ionising photon production per unit UV luminosity. The second,
$f_{\textrm{esc}}$
, is the HI ionising escape fraction, which is the fraction of ionising photons produced within a galaxy that escape into the IGM. The galaxy UV luminosity function (UVLF) has been measured using JWST up to
$z = 14$
(e.g. Pérez-González et al. Reference Pérez-González2023; Adams et al. Reference Adams2023; Donnan et al. Reference Donnan2024; Harikane et al. Reference Harikane2024), facilitating direct measurements of
$\rho_{\textrm{UV}}$
up to that redshift. As such, an independent measurement of
$\dot{N}_{\textrm{ion}}$
using constraints on the IGM ionisation state would enable a measurement of
$\langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
. One such approach is to combine information about the photo-ionisation rate in the ionised IGM inferred from the Ly
$\alpha$
forest of high-redshift quasars (Wyithe & Bolton Reference Wyithe and Bolton2011; Calverley et al. Reference Calverley, Becker, Haehnelt and Bolton2011; Becker & Bolton Reference Becker and Bolton2013; D’Aloisio et al. Reference D’Aloisio, McQuinn, Davies and Furlanetto2018; Bosman et al. Reference Bosman2022) with measurements of the mean free path (MFP) to ionising photons (Prochaska,Worseck, & O’Meara Reference Prochaska, Worseck and O’Meara2009;Worseck et al. Reference Worseck2014; Becker et al. Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021; Zhu et al. Reference Zhu2023). The former contains information about the number of ionising photons in the IGM, and the latter about how quickly those photons are being absorbed. They can be combined to indirectly measure
$\dot{N}_{\textrm{ion}}$
at
$z \leq 6$
, where measurements of both quantities are available.
The exercise described above has been carried out by a number of previous works in the context of both HI and He II reionisation (e.g. Bolton & Haehnelt Reference Bolton and Haehnelt2007; Kuhlen & Faucher-Giguère Reference Kuhlen and Faucher-Giguère2012; Becker & Bolton Reference Becker and Bolton2013; Khaire et al. Reference Khaire2017).Footnote
a
Most recently, Bosman & Davies (Reference Bosman and Davies2024) (hereafter B24) used measurements of
$\rho_{\textrm{UV}}$
from Bouwens et al.(Reference Bouwens2021) to measure
$\langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
given an updated measurement of
$\dot{N}_{\textrm{ion}}$
at
$z = 4 - 6$
by Gaikwad et al.(Reference Gaikwad2023). Assuming galaxies produce ionising photons down to a limiting magnitude of
$M_{\textrm{UV}} = -11$
, they found
$\log \langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}} = 24.28^{+0.21}_{-0.20}$
erg
$^{-1}$
Hz at
$z = 5$
and
$\log \langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}} = 24.66^{+0.18}_{-0.47}$
erg
$^{-1}$
Hz at
$z = 6$
. They also report an upper limit at
$z = 4$
of
$\log \langle f_{\textrm{esc}} \xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
$\lt 24.11$
erg
$^{-1}$
Hz. Notably, these measurements are up to a factor of
$\sim 5$
lower than estimated in Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) (for
$M_{\textrm{UV}} \lt -11$
, Figure 1 of B24) using direct measurements of
$\xi_{\textrm{ion}}$
from JWST Simmonds et al.(Reference Simmonds2023) and inferences on
$f_{\textrm{esc}}$
derived from low-redshift measurements Chisholm et al.(Reference Chisholm2022), suggesting a downward revision on the photon budget in the late stages of reionisation.
Collection of measurements of
$\Gamma_{\textrm{HI}}$
(top) and
$\lambda_{912}^{\textrm{mfp}}$
(bottom) used to measure
$\dot{N}_{\textrm{ion}}$
in this work. The red dashed curves shows the maximum-likelihood fits to each set of measurements. The thin black curves show random draws from the posteriors of the model parameters. The cyan dot-dashed curve in the bottom panel shows the ‘ionised phase’ MFP estimated using Equation (7). Measurements of
$\Gamma_{\textrm{HI}}$
are from Becker & Bolton (Reference Becker and Bolton2013), Bosman et al. (Reference Bosman2022), Gaikwad et al. (Reference Gaikwad2023). Measurements of
$\lambda_{912}^{\textrm{mfp}}$
are from Fumagalli et al. (Reference Fumagalli, O’Meara, Xavier Prochaska and Worseck2013), O’Meara et al. (Reference O’Meara, Xavier Prochaska, Worseck, Chen and Madau2013), Worseck et al. (Reference Worseck2014), Zhu et al. (Reference Zhu2023), Gaikwad et al. (Reference Gaikwad2023), Gao et al. (Reference Gao2024).
The measurements of
$\dot{N}_{\textrm{ion}}$
by B24, and most other works at these redshifts, neglect the effects of radiative transfer (RT) and reionisation. They rely on the so-called ‘local source approximation’ (LSA), which assumes that the absorption rate of ionising photons in the IGM is equal to the emission rate. This is a good approximation during the bulk of reionisation, when the MFP is very short owing to much of the IGM being significantly neutral. However, at
$5 \lt z \lt 6$
, when reionisation is likely in its ending stages (Kulkarni et al. Reference Kulkarni, Keating, Haehnelt, Bosman, Puchwein, Chardin and Aubert2019; Keating et al. Reference Keating, Kulkarni, Haehnelt, Chardin and Aubert2020; Nasir & D’Aloisio Reference Nasir and D’Aloisio2020), both the ionising background and the MFP evolve very rapidly. At this point, photons may experience a significant delay between emission and absorption. In this regime, an explicit treatment of RT effects may be required to accurately recover
$\dot{N}_{\textrm{ion}}$
.
Another effect is that of possible ‘islands’ of neutral hydrogen, which may still be present in the IGM down to redshifts as low as
$z = 5$
. Typically, measurements of
$\dot{N}_{\textrm{ion}}$
assume that the IGM is fully ionised. However, if neutral islands are present at
$z \gt 5$
, these likely affect measurements of the MFP (Roth et al. Reference Roth, D’Aloisio, Cain, Wilson, Zhu and Becker2024; Satyavolu et al. Reference Satyavolu, Kulkarni, Keating and Haehnelt2024; Chen, Fan, & Avestruz Reference Chen, Fan and Avestruz2024), and possibly the ionising background. This is of particular concern at
$z = 6$
, when reionisation was likely still ongoing and the global neutral fraction may be as high as 20% (Zhu et al. Reference Zhu2024b; Spina et al. Reference Spina, Bosman, Davies, Gaikwad and Zhu2024). As such, accounting for the presence of neutral islands at
$z = 6$
may be necessary.
Measurements of
$\dot{N}_{\textrm{ion}}$
at these redshifts provide valuable insight into the evolution of galaxy properties and the tail end of the reionisation process. Recently, JWST has begun to measure
$\xi_{\textrm{ion}}$
for a large sample of galaxies up to
$z \approx 9$
(e.g. Simmonds et al. Reference Simmonds2023; Zhu et al. Reference Zhu2024a). For the first time, these measurements statistically characterise the dependence of
$\xi_{\textrm{ion}}$
on both redshift and
$M_{\textrm{UV}}$
at redshifts extending into the reionisation epoch. Combining these new results with forest-based measurements of
$\dot{N}_{\textrm{ion}}$
will allow us to measure the average escape fraction of the galaxy population up to
$z = 6$
more robustly than has been previously possible. Measuring the ionising output of galaxies at
$z \gt 5$
will also begin to constrain the number of ionising photons required to complete reionisation – the ionising photon budget. The photon budget is sensitive to the ionising opacity of the ionised IGM, a major source of uncertainty in reionisation studies (So et al. Reference So, Norman, Reynolds and Wise2014; Becker et al. Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021; Cain et al. Reference Cain, D’Aloisio, Gangolli and Becker2021; Davies et al. Reference Davies, Bosman, Furlanetto, Becker and D’Aloisio2021).
In this work, we present updated measurements of
$\dot{N}_{\textrm{ion}}$
at
$2.5 \leq z \leq 6$
, taking into account RT effects, and present new formalism to account for the presence of neutral islands in measurements of
$\dot{N}_{\textrm{ion}}$
. We interpret our findings in the context of high-redshift galaxy observations, and infer estimates of the escaping ionising efficiency and the population-averaged escape fraction. This work is organised as follows. In Section 2, we present the formalism used to measure
$\dot{N}_{\textrm{ion}}$
. We describe the observations used in our analysis and our modelling uncertainties in Section 3, present our main results in Section 4, and conclude in Section 5. Throughout, we assume the following cosmological parameters:
$\Omega_m = 0.305$
,
$\Omega_{\Lambda} = 1 - \Omega_m$
,
$\Omega_b = 0.048$
,
$h = 0.68$
,
$n_s = 0.9667$
and
$\sigma_8 = 0.82$
, consistent with Planck Collaboration et al. Reference Collaboration2020 results. All distances are co-moving unless otherwise specified.
2. Formalism
2.1. Radiative transfer equation
In the LSA, the IGM photo-ionisation rate
$\Gamma_{\textrm{HI}}$
can be expressed as
$\begin{equation} \Gamma_{\textrm{HI}} = (1+z)^3\int d\nu \dot{N}_{\textrm{ion}}^{\nu} \lambda_{\nu} \sigma_{\textrm{HI}}^{\nu}\end{equation}$
where
$\dot{N}_{\textrm{ion}}^{\nu}$
is the co-moving ionising emissivity output by sources per unit frequency
$\nu$
,
$\lambda_{\nu}$
is the MFP to ionising photons, and
$\sigma_{\textrm{HI}}^{\nu}$
is the HI-ionising cross-section. Equation (2) supposes that (1) photons are absorbed within a short time after being emitted and (2) that the IGM is highly ionised. A more general expression that does not rely on the first assumption is
$\begin{equation} \Gamma_{\textrm{HI}} = (1+z)^3\int d\nu \dot{N}_{\textrm{abs}}^{\nu} \lambda_{\nu} \sigma_{\textrm{HI}}^{\nu}\end{equation}$
where we have replaced the emission rate
$\dot{N}_{\textrm{ion}}^{\nu}$
with the absorption rate,
$\dot{N}_{\textrm{abs}}^{\nu}$
, no longer assuming these to be equal. This distinction is important whenever the timescale for an ionising photon to travel one mean free path,
$\lambda_{\nu}/c$
, is significant compared to the timescales over which IGM properties evolve. The absorption rate at time t is given by
$\begin{equation} \dot{N}_{\textrm{abs}}^{\nu}(t) = \int_{0}^{t} dt' \dot{N}_{\textrm{ion}}^{\nu}(t',\nu') G(t,t',\nu,\nu')\end{equation}$
where photons emitted at time
$t' \leq t$
(
$z' \geq z$
) with frequency
$\nu' = \frac{1+z'}{1+z}\nu$
red-shifting to observed frequency
$\nu$
at time t, and
$\dot{N}_{\textrm{ion}}^{\nu}(t',\nu')$
is the emissivity per unit co-moving frequency (
$\nu$
) evaluated at
$\nu'$
. Here, G is given by
$\begin{equation} G(t,t',\nu,\nu') = c \kappa_{\nu}(t) \exp\left[-\int_{t'}^{t} dt'' c \kappa_{\nu''}(t'')\right]\end{equation}$
where
$\kappa_{\nu} \equiv \lambda_{\nu}^{-1}$
is the absorption coefficient, c is the speed of light, and
$\nu'' = \frac{1+z''}{1+z}\nu$
. G is the Green’s function for the cosmological RT equation (e.g. in Haardt & Madau Reference Haardt and Madau1996, Reference Haardt and Madau2012) with the usual proper emissivity and angle-averaged intensity replaced by the co-moving ionising photon emissivity and absorption rate, respectively. It quantifies how the ionising background at time t is built up by photons emitted at earlier times,
$t' \lt t$
. In the limit that
$1/(c \kappa_{\nu}) \rightarrow 0$
(the short MFP limit),
$G(t,t',\nu,\nu')$
becomes the Dirac delta function
$\delta_{\textrm{D}}(t-t',\nu-\nu')$
, in which case we recover the LSA. We show in Appendix A that our formulation is equivalent to the usual way of expressing the solution to the cosmological RT equation. The advantage to our formulation is that Equations (3)–(5) show directly the relationships between
$\lambda$
,
$\Gamma_{\textrm{HI}}$
, and
$\dot{N}_{\textrm{ion}}$
.
2.2. Accounting for neutral islands
Equation (3) is only accurate if all photons are absorbed by the ionised IGM and contribute to its total ionisation rate. Photons absorbed by any neutral islands will not contribute to
$\Gamma_{\textrm{HI}}$
in the ionised IGM, which is what the Ly
$\alpha$
forest probes. As such, it is unclear whether the standard approach to measuring
$\dot{N}_{\textrm{ion}}$
will give the right answer in an IGM containing neutral islands. A modified version of Equation (3) accounting for this is
$\begin{equation} \Gamma_{\textrm{HI}} = \int d\nu \dot{N}_{\textrm{abs,ionised}}^{\nu} \lambda_{\nu,\textrm {ionised}} \sigma_{\textrm{HI}}^{\nu}\end{equation}$
where
$\dot{N}_{\textrm{abs,ionised}}^{\nu}$
is the absorption rate in ionised gas, and
$\lambda_{\nu}^{\textrm{ionised}}$
is the MFP in ionised gas. These quantities are defined such that the contribution of neutral islands to both is removed. We can estimate
$\lambda_{\nu,\textrm{ionised}}$
by
$\begin{equation} \lambda_{\nu,\textrm{ionised}}^{-1} = \lambda_{\nu}^{-1} - \lambda_{\nu,\textrm{neutral}}^{-1}\end{equation}$
where
$\lambda_{\nu,\textrm{neutral}}^{-1}$
is the contribution to the absorption coefficient from neutral islands, i.e. not counting the opacity of the ionised gas. Similarly, we can estimate
$\dot{N}_{\textrm{abs,ionised}}^{\nu}$
using
$\begin{equation} \dot{N}_{\textrm{abs,ionised}}^{\nu} = \dot{N}_{\textrm{abs}}^{\nu} - \dot{N}_{\textrm{ abs,neutral}}^{\nu}\end{equation}$
where the absorption rate by neutral islands only is given by
$\begin{equation} \int d\nu \dot{N}_{\textrm{abs,neutral}}^{\nu} = -\frac{dx_{\textrm{HI}}^{\textrm{m}}}{dt} (1+\chi) n_{\textrm{H}}\end{equation}$
Here,
$x_{\textrm{HI}}^{\textrm{m}}$
is the mass-weighted HI fraction in the IGM,
$n_{\textrm{H}}$
is the cosmic mean hydrogen density, and the factor of
$1+\chi \equiv 1.082$
accounts for single ionisation of He. Equation (9) simply counts the net rate at which the fully neutral IGM is being ionised, which is equivalent to the absorption rate by neutral islands only. Since the global reionisation history is unknown, in our analysis we must assume it (and the MFP to neutral islands,
$\lambda_{\nu,\textrm{neutral}}$
) from a simulation, which we discuss further in the next section. We explore the effect of accounting for neutral islands in this way in Appendix B.
In sum, our RT formalism takes into account several effects that are missing in the LSA. First, it accounts for red-shifting of ionising photons past the Lyman Limit, which becomes important at
$z \lesssim 4$
(Becker & Bolton Reference Becker and Bolton2013). Second, it accounts for the rapid buildup of the ionising background immediately following reionisation, which necessitates
$\dot{N}_{\textrm{ion}}^{\nu} \gt \dot{N}_{\textrm{abs}}^{\nu}$
. Lastly, including neutral islands accounts for the possibility that reionisation may be ongoing, such that not all absorptions occur in the highly ionised IGM.
3. Method
3.1. Measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda$
We use measurements of
$\Gamma_{\textrm{HI}}$
and the Lyman Limit MFP,
$\lambda_{912}^{\textrm{mfp}}$
, to measure
$\dot{N}_{\textrm{ion}}$
at
$2.5 \lt z \lt 6$
. We show these measurements in Figure 1 – see the legend and caption for references.Footnote
b
In practice, solving Equations (3)–(5) for
$\dot{N}_{\textrm{ion}}$
requires a smooth functional form for both
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
(see Section 3.2). We fit the collection of measurements to smooth functions of redshift using an MCMC approach, accounting for the reported
$1\sigma$
error bars on the measurements using a standard Gaussian likelihood.Footnote
c
We use a 5th-order polynomial in redshift to fit the
$\log_{10}\Gamma_{\textrm{HI}}$
and a double power-law to fit the
$\lambda_{912}^{\textrm{mfp}}$
measurements (see Appendix C for details). The maximum likelihood fit is shown by the red dashed curve in each panel of Figure 1. The thin black lines are random draws from the recovered posteriors on the best-fit parameters, which approximately capture the error in the best-fits to the measurements (see Section 3.3). The cyan dot-dashed curve in the bottom panel shows the best fit with the neutral island opacity subtracted (Equation 7, and see below).
A key source of uncertainty is the choice of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements used at
$z \geq 5$
. For both
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
, the measurements in Figure 1 represent two different approaches to measuring these quantities. The standard way to measure
$\Gamma_{\textrm{HI}}$
is to match the mean transmission of the Ly
$\alpha$
forest in hydrodynamical simulations with that measured in the spectra of high-redshift quasars. Typically, these simulations assume that reionisation is over and that the UV background (UVB) is homogeneous. Significant sources of uncertainty in such measurements include the thermal history of the IGM, which is typically marginalised over (Becker & Bolton Reference Becker and Bolton2013), and numerical convergence in both spatial resolution and box size (Doughty et al. Reference Doughty, Hennawi, Davies, Lukic and Oñorbe2023). The measurements of Becker & Bolton (Reference Becker and Bolton2013), Bosman et al. (Reference Bosman, Fan, Jiang, Reed, Matsuoka, Becker and Haehnelt2018),Footnote
d
and Becker et al. (Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021) in the top panel were done in this way. An important caveat is that the Bosman et al. (Reference Bosman2022) measurements were done assuming an early reionisation model with a relatively low IGM temperature, and thus may be biased high (see below).
Direct measurements of
$\lambda_{912}^{\textrm{mfp}}$
use the Lyman Continuum (LyC) spectrum of high-redshift quasars. The observed MFP at a given redshift is defined as the distance travelled by a photon emitted at the redshift of the quasar that reaches a LyC opacity of unity when it redshifts to
$912\, \mathring{A} $
. This quantity is estimated by stacking and fitting LyC spectra of bright quasar spectra (Prochaska, Worseck, & O’Meara Reference Prochaska, Worseck and O’Meara2009). All the measurements in Figure 1 use this approach except those of Gaikwad et al. (Reference Gaikwad2023, see below). One important caveat is that measurements using this method at
$z \gt 5$
are challenging and must account for the effect of the proximity zone of the quasar on the LyC spectrum. This may introduce additional model-dependence and potential systematic uncertainty in the measurement (Becker et al. Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021).
Recently, Gaikwad et al. (Reference Gaikwad2023, hereafter G23) proposed a new method of jointly measuring
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
based on the fluctuations in the Ly
$\alpha$
forest on large spatial scales (see also Davies et al. Reference Davies2024b).Their approach is somewhat complementary to the standard approaches, since it uses different information available in the quasar spectra. It relies on semi-numerical simulations that assume a relationship between
$\lambda_{912}^{\textrm{mfp}}$
and large-scale fluctuations in both
$\Gamma_{\textrm{HI}}$
and the IGM density, which makes it (to some degree) model-dependent. Crucially, however, it takes into account spatial fluctuations in the UVB and the possible presence of neutral islands near the end of reionisation (
$z \gt 5$
). Their measurements are the blue points in both panels of Figure 1. Note that all the most recent
$z \gt 5$
measurements of both
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
use quasar spectra from the XQR-30 data set (D’Odorico et al. Reference D’Odorico2023).
To gauge the differences between
$\dot{N}_{\textrm{ion}}$
implied by these two types of measurements at
$z \gt 5$
, we also fit two sub-sets of the data. The first includes only the G23 measurements (for both quantities) at
$z \geq 5$
and the second excludes only the G23 results. Both sub-sets include the
$\Gamma_{\textrm{HI}}$
measurements of Becker & Bolton (Reference Becker and Bolton2013) at
$z \lt 5$
and all direct
$\lambda_{912}^{\textrm{mfp}}$
measurements at
$z \lt 5$
. These fits are shown in Appendix C and are used in the analysis below. The parameters for the maximum-likelihood fits to our full data set and these two reduced data sets are also in Appendix C. We note that differences in
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
between the two sub-sets of measurements contribute at similar levels to differences in measured quantities at
$z \gt 5$
.
3.2. Estimation of
$\dot{N}_{\textrm{ion}}$
To measure
$\dot{N}_{\textrm{ion}}$
from
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
we need to the frequency (
$\nu$
) dependence of the spectrum of the emitted ionising radiation, and the frequency dependence of
$\lambda_{\nu}$
. We model the former as a power law in
$\nu$
, such that
$\begin{equation} \dot{N}_{\textrm{ion}}^{\nu} \propto \frac{\nu^{-\alpha}}{h_p \nu} \propto \nu^{-(\alpha+1)}\end{equation}$
where
$h_p$
is Planck’s constant, for photon energies between 1 and 4 Ryd. The cutoff at
$E \gt 4$
Ryd is appropriate for galaxy spectra, and a reasonable approximation for quasars prior to the onset of He II reionisation. At
$z \lesssim 4$
, when He II reionisation is likely underway, emission at energies greater than 4 Ryd from quasars contributes to the ionising background in most of the IGM. We note, however, that this effect on our
$\dot{N}_{\textrm{ion}}$
measurements is
$\lt 15\%$
over most of our parameter space, well below the total error budget (see next section).Footnote
e
Assuming the HI column density distribution has the form of a power law with slope
$\beta_{\textrm{N}}$
, we can write (McQuinn, Oh, & Faucher-Giguère 2011),
$\begin{equation} \lambda_{\nu} \propto (\sigma_{\textrm{HI}}^{\nu})^{-(\beta_{\textrm{N}}-1)} \propto \nu^{2.75(\beta_{\textrm{N}} - 1)}\end{equation}$
where we have approximated
$\sigma_{\textrm{HI}} \propto \nu^{-2.75}$
in the frequency range of interest. The limit that
$\beta_{\textrm{N}} = 1$
(no frequency dependence) is that of the IGM opacity being dominated by highly opaque (
$\tau \gt\gt 1$
) gas. In the opposite limit (
$\beta_{\textrm{N}} = 2$
), it is dominated by the diffuse, highly ionised IGM. Smaller
$\alpha$
corresponds to harder (more energetic) ionising spectra, which results in lower
$\Gamma_{\textrm{HI}}$
at fixed
$\dot{N}_{\textrm{ion}}$
because of the frequency dependence of
$\sigma_{\textrm{HI}}$
. Thus, assuming a smaller
$\alpha$
requires a larger
$\dot{N}_{\textrm{ion}}$
at fixed
$\Gamma_{\textrm{HI}}$
. A smaller
$\beta_{\textrm{N}}$
decreases the frequency-averaged MFP at fixed
$\lambda_{912}^{\textrm{mfp}}$
, which also increases
$\dot{N}_{\textrm{ion}}$
at fixed
$\Gamma_{\textrm{HI}}$
. Following B24, we assume fiducial values of
$\alpha = 2$
and
$\beta_{\textrm{N}} = 1.3$
. We will vary these parameters in the ranges
$1 \leq \alpha \leq 3$
and
$1 \leq \beta_{\textrm{N}} \leq 2$
. Note that our choice of
$\alpha = 2$
is motivated by models of metal-poor (Pop II) galaxy SEDs (Bressan et al. Reference Bressan, Marigo, Girardi, Salasnich, Cero, Rubele and Nanni2012; Choi, Conroy, & Byler Reference Choi, Conroy and Byler2017), but is also reasonable for quasars (Lusso et al. Reference Lusso, Worseck, Hennawi, Prochaska, Vignali, Stern and O’Meara2015), which probably dominate the ionising output of the source population at
$z \lt 4$
(see, however, Madau et al. Reference Madau, Giallongo, Grazian and Haardt2024 for higher z). Our choice of
$\beta_{\textrm{N}}$
is commonly assumed in the literature (e.g. Gaikwad et al. Reference Gaikwad2023) and is motivated by the best-fit to the HI column density distribution Becker & Bolton (Reference Becker and Bolton2013).
The integration in Equation (4) runs over all
$t' \lt t$
, so formally it should start from
$z = \infty$
(
$t = 0$
). In practice, it is not possible to do this, since we have measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
only up to
$z = 6$
. To approximate Equation (4), we extrapolate our maximum-likelihood fits for
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
to
$z = 6.5$
, and assume that at that redshift the LSA is valid. This lets us estimate
$\dot{N}_{\textrm{ion}}(z = 6.5)$
using Equation (2). We then set the lower limit of the integration in Equation (4) to
$t(z = 6.5)$
and adjust the value of
$\dot{N}_{\textrm{ion}}$
at each successively lower redshift until Equation (3) returns the measured
$\Gamma_{\textrm{HI}}$
. In this way, we find the
$\dot{N}_{\textrm{ion}}(z)$
down to
$z = 2.5$
that satisfies Equations (3)–(5) for our assumed
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
. Our results at
$2.5 \leq z \leq 6$
are only sensitive at the few-percent level (or less) to how
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
are extrapolatedFootnote
f
to
$z = 6.5$
.
Since reionisation may be ongoing at
$z \gt 5$
, we also use Equations (6)–(9) to measure
$\dot{N}_{\textrm{ion}}$
assuming a late reionisation history ending at
$z \approx 5$
. To do this, we need (1) a reionisation history, to evaluate Equation (9), (2) a neutral island MFP to evaluate Equation (7) and (3) a functional form for the spectrum of ionising radiation absorbed by neutral islands (to extract
$\dot{N}_{\textrm{abs,neutral}}^{\nu}$
from the integral in Equation 9). To satisfy (1) and (2), we use the late start/late end model from Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024), which is a ray-tracing RT simulation of the EoR run with the FlexRT code of Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024). The
$\dot{N}_{\textrm{ion}}$
in the simulation is calibrated to produce agreement with the mean Ly
$\alpha$
forest transmission measurements of Bosman et al. (Reference Bosman2022) at
$5 \leq z \leq 6$
. It has a volume-averaged neutral fraction of 30% at
$z = 6$
, making it a fairly extreme case of late reionisation (for discussion, see Zhu et al. Reference Zhu2024b).As such, our estimate of the effect of neutral islands on measurements of
$\dot{N}_{\textrm{ion}}$
is likely an over-estimate, to be viewed as a rough upper limit. However, by comparing models with and without neutral islands, we will see that they have only a mild effect on our results. We evaluate the MFP to neutral islands in the simulation by setting the opacity in ionised gas to 0 and then estimating
$\lambda$
using the definition described in Appendix C of Chardin et al. (Reference Chardin, Haehnelt, Aubert and Puchwein2015) (Equation 5 Cain et al. Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024). For (3), we assume (for simplicity) that the the spectrum of
$\dot{N}_{\textrm{abs,neutral}}^{\nu}$
is the sameFootnote
g
as that of
$\dot{N}_{\textrm{abs}}^{\nu}$
, allowing us to straightforwardly evaluate Equation (8).
3.3. Errors on
$\dot{N}_{\textrm{ion}}$
Errors on
$\dot{N}_{\textrm{ion}}$
arise from uncertainty in
$\Gamma_{\textrm{HI}}$
,
$\lambda_{912}^{\textrm{mfp}}$
,
$\alpha$
, and
$\beta_{\textrm{N}}$
. To estimate uncertainties on
$\dot{N}_{\textrm{ion}}$
, we draw
$1\,000$
random sets of parameters from the posteriors on the MCMC fits to
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
and calculate
$\dot{N}_{\textrm{ion}}$
for each combination. For each draw, we also randomly draw values of
$\alpha$
and
$\beta_{\textrm{N}}$
from the uniform distributions
$\alpha \in [1,3]$
and
$\beta_{\textrm{N}} \in [1,2]$
. At each redshift, we treat the 13–87% (2.5–97.5%) range of the resulting
$\dot{N}_{\textrm{ion}}$
values as the
$1 \sigma$
(
$2 \sigma$
) spread around our maximum likelihood fiducial measurement. Note that our
$1\sigma$
and
$2\sigma$
ranges are not, in general, symmetric around our maximum likelihood result.Footnote
h
We do this for our fiducial fit using all data points, and our reduced data sets including and excluding the G23 measurements.
To quantify how important each parameter is to determining the errors in
$\dot{N}_{\textrm{ion}}$
, we re-run this analysis allowing only one parameter at a time to vary, holding all the others fixed to their maximum likelihood fits or fiducial values. In Table 1, we report the two-sided
$1\sigma$
(
$2\sigma$
in parentheses) logarithmic errors on
$\dot{N}_{\textrm{ion}}$
arising from each parameter individually at
$z = 2.5$
, 4, 5, and 6. The bottom row reports the total errors on
$\dot{N}_{\textrm{ion}}$
. At
$z = 2.5$
, uncertainties in
$\Gamma_{\textrm{HI}}$
dominate the error budget. At
$z = 4$
and 5,
$\Gamma_{\textrm{HI}}$
and
$\beta_{\textrm{N}}$
dominate, and contribute at roughly equal levels. At
$z = 6$
, uncertainties in
$\lambda_{912}^{\textrm{mfp}}$
become important at a level similar to
$\Gamma_{\textrm{HI}}$
and
$\beta_{\textrm{N}}$
. Uncertainties from
$\alpha$
are sub-dominate at all redshifts. The full redshift dependence of these one-parameter errors is shown in Appendix D.
Estimate of the error budget for our fiducial
$\dot{N}_{\textrm{ion}}$
measurement at several redshifts. The bottom row reports the total logarithmic errors on the measurements, and the rows above give an estimate of the contribution from each uncertain quantity in the analysis. We report
$\pm 1\sigma$
errors and
$2 \sigma$
errors in parentheses.
We mention here a couple of caveats that likely render our uncertainties on
$\dot{N}_{\textrm{ion}}$
too small. The first is that they rely on parametric fits to measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
. These do not take into account the full covariance between uncertainties on different measurements, and the scatter in the posteriors is likely too small on account of the relatively small number of parameters (
$4-5$
) used in the fits. Another caveat is that the two different types of measurements discussed in Section 3.1 are systematically offset from each other at
$z \gt 5$
, as can be clearly seen in Figure 1. Using both in the MCMC fits will thus produce overly tight posteriors at
$z \gt 5$
. In what follows, we will show results using each set of measurements individually, allowing us to quantify better the real uncertainties in
$\dot{N}_{\textrm{ion}}$
at these redshifts. We discuss these caveats in more detail in Appendix D.
4. Results
4.1. Measurements of
$\dot{N}_{\textrm{ion}}$
We show our fiducial measurements of
$\dot{N}_{\textrm{ion}}$
at
$2.5 \lt z \lt 6$
in the top panel of Figure 2. The solid black curve shows our max-likelihood measurement, and the shaded region indicates the
$1\sigma$
uncertainties. Since our measurement and error bars at each redshift are a result of a cumulative integral over higher redshifts, we show them as a continuous curve rather than as discrete points at particular redshifts. We show, for comparison, the
$\dot{N}_{\textrm{ion}}$
measurements from G23 at
$z = 5$
and 6 (blue points) and those by Becker & Bolton (Reference Becker and Bolton2013) at
$2 \lt z \lt 5$
(green points). Notably, our measurement is
$\approx 0.8\sigma$
higher than that of G23 at
$z = 5$
and
$\approx 2 \sigma$
higher at
$z = 6$
. At
$z \lt 5$
, our measurement is close to the central values of the Becker & Bolton (Reference Becker and Bolton2013) measurements.
Measurements of
$\dot{N}_{\textrm{ion}}$
compared to previous measurements. Top: our fiducial measurement (black solid curve) at
$2.5 \lt z \lt 6$
. The dark (light) shaded region indicates the approximate
$1\sigma$
(
$2\sigma$
) uncertainties. The red dashed curve shows estimates of
$\dot{N}_{\textrm{ion}}$
using the LSA, and the black dotted curve includes neglects red-shifting of ionising photons past the Lyman Limit, but still accounts for the finite travel time of photons. The cyan dot-dashed curve shows our revised measurement accounting for the presence of neutral islands at
$z \gt 5$
. Bottom: the same measurement (including only
$1\sigma$
uncertainties) using our reduced data sets with only the G23 points at
$z \geq 5$
(blue dashed curve) and the excluding only the G23 points (red dot-dashed curve). The differences between these are small at
$z \lt 5$
, but become significant at
$z \gt 5$
. In the former case,
$\dot{N}_{\textrm{ion}}$
remains nearly flat up to
$z = 6$
, while in the latter,
$\dot{N}_{\textrm{ion}}$
grows by a factor of
$\approx 4$
between
$z = 5$
and 6. In our fiducial measurement, this increase is a factor of
$\approx 2$
.
In our fiducial measurement,
$\dot{N}_{\textrm{ion}}$
falls by a factor of
$\approx 2$
between
$z = 6$
and 5. This is because the measured
$\lambda_{912}^{\textrm{mfp}}$
grows with time more quickly over that redshift range than
$\Gamma_{\textrm{HI}}$
(since
$\dot{N}_{\textrm{ion}} \sim \Gamma_{\textrm{HI}}/\lambda_{912}^{\textrm{mfp}}$
). We find roughly constant
$\dot{N}_{\textrm{ion}}$
between
$z = 5$
and 4, with a small dip around
$z = 3-3.5$
and a sharp increase towards
$z = 2.5$
. This last feature is driven by the slight upturn in the Becker & Bolton (Reference Becker and Bolton2013)
$\Gamma_{\textrm{HI}}$
measurements at
$z \lt 3.5$
.Footnote
i
Our
$\pm 1\sigma$
range is a factor of
$\approx 2$
at most redshifts, but increases considerably at
$z \gt 5.5$
due to rising uncertainty in
$\lambda_{912}^{\textrm{mfp}}$
(see Table 1).
The red-dashed curve in Figure 2 shows the effect of using the LSA to compute
$\dot{N}_{\textrm{ion}}$
instead of the full formalism accounting for RT. At
$z \approx 6$
, this agrees well with the full measurement, thanks to the short MFP. However, at
$z \lt 6$
, they start to diverge, and by
$z = 5$
the LSA gives a result
$\approx 20\%$
below the full calculation. This effect explains over half the difference between our measurement and that of G23, and the rest arises from differences in the assumed
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
. At
$z \lt 5$
, the difference between the LSA and the full measurement increases, reaching a factor of
$\approx 3$
by
$z = 2.5$
. Previous works (including Becker & Bolton Reference Becker and Bolton2013) have noted that at
$z \lt 4$
, the LSA is expected to fail because it ignores the red-shifting of ionising photons past the Lyman Limit before they can ionise a neutral atom. The black dotted curve isolates this effect by explicitly neglecting the red-shifting of photons while still accounting for the distance they travel through the IGM before being absorbed.Footnote
j
This result differs from the full measurement by at most
$\approx 25\%$
, and explains less than half of the difference between the black solid and red dashed curves most of the time. This indicates that the steady buildup of ionising photons in the IGM due to the lengthening MFP contributes more to the failure of the LSA than does redshifting, even at low redshifts.
The cyan dot-dashed curve includes the effect of neutral islands at
$z \gt 5$
, when reionisation is ongoing in our assumed model (vertical shaded region). We see that taking islands into account changes the result by at most 15%. We can understand why by returning to Equation (6). In that equation,
$\dot{N}_{\textrm{abs,ionised}}^{\nu}$
is smaller than the total absorption rate,
$\dot{N}_{\textrm{ion}}^{\nu}$
, since some ionising photons are consumed by the neutral islands (Equation 9). However, the MFP in ionised gas is also larger than that in the IGM at large (Equation 7). It turns out that the product
$\dot{N}_{\textrm{abs,ionised}}^{\nu} \lambda_{\nu,\textrm {ionised}}$
is similar to
$\dot{N}_{\textrm{abs}}^{\nu} \lambda_{\nu}$
, meaning we recover the nearly the same total
$\dot{N}_{\textrm{ion}}$
that we would if we had ignored islands. We investigate this result more carefully in Appendix B, and find it to be true even when
$\dot{N}_{\textrm{abs,neutral}} \gt \dot{N}_{\textrm{abs,ionised}}$
(that is, when islands dominate the global absorption rate).
The bottom panel of Figure 2 shows how our results change when only the G23 measurements (Figure 1) are used at
$z \geq 5$
(blue dashed) and when they are excluded (red dot-dashed). While the results are similar at
$z \lt 5$
, they diverge considerably at
$z \gt 5$
. The maximum likelihood result using only the G23 measurements declines slightly with redshift at
$z \gt 5$
. However, without G23, we find a factor of
$\approx 4$
increase in
$\dot{N}_{\textrm{ion}}$
from
$z = 5$
to 6 in our maximum likelihood result, and scenarios with flat or decreasing
$\dot{N}_{\textrm{ion}}$
are clearly disfavoured. The spread between these two sets of results better captures the uncertainties in
$\dot{N}_{\textrm{ion}}$
at
$z \gt 5$
, which are artificially tightened when combining all measurements (see Section 3.3).
The differences between the blue and red curves arise from the fact that the
$\Gamma_{\textrm{HI}}$
measured by G23 decreases more quickly, and
$\lambda_{912}^{\textrm{mfp}}$
less quickly, with redshift than in the other subset of measurements. Notably, their measurement of
$\lambda_{912}^{\textrm{mfp}}$
at
$z = 6$
is a factor of 2 higher than the direct measurements from Zhu et al. (Reference Zhu2023). Their
$z = 6$
$\Gamma_{\textrm{HI}}$
measurement is also a factor of
$\approx 2$
below that of Becker et al. (Reference Becker, D’Aloisio, Christenson, Zhu, Worseck and Bolton2021) at
$z = 6$
. In the G23-only subset,
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
decrease with redshift at the roughly the same rate, keeping
$\dot{N}_{\textrm{ion}}$
approximately constant. In the other subset,
$\Gamma_{\textrm{HI}}$
decreases much more slowly than
$\lambda_{912}^{\textrm{mfp}}$
, causing
$\dot{N}_{\textrm{ion}}$
to grow rapidly. The difference shrinks to
$\approx 20\%$
at
$z = 5$
and disappears by
$z = 4$
.
4.2. Ionising photon budget
Several recent works have suggested, based on Ly
$\alpha$
forest observations, that reionisation ended at
$z = 5 - 5.5$
. If this is the case, we can use our measurements to estimate the ionising photon output per H atom,
$N_{\gamma/H}$
, at the tail end of reionisation. Assuming reionisation ended at
$z = 5.3$
, as suggested by Bosman et al. (Reference Bosman2022), we define the
$5.3 \lt z \lt 6$
photon budget to be the number of ionising photons produced per H atom between
$z = 5.3$
and 6,
$N_{\gamma/H}^{z = 5.3-6}$
. We measure
$N_{\gamma/H}^{z=5.3-6} = 1.10_{-0.39}^{+0.21}$
for our fiducial measurement at
$1\sigma$
confidence, and
$0.61_{-0.22}^{+0.13}$
(
$1.95_{-0.98}^{+0.93}$
) using our sub-sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
with (without) the G23 results, respectively.
If the universe were 20% neutral at
$z = 6$
, as suggested by recent models, the minimum budget required to complete reionisation by
$z = 5.3$
in the absence of recombinations is
$N_{\gamma/H}^{z = 5.3-6} = 0.216$
(which accounts for single ionisation of helium). We can consider the tail-end of reionisation to be ‘absorption-dominated’ if the actual budget is at least twice this value (see Davies et al. Reference Davies, Bosman, Furlanetto, Becker and D’Aloisio2021), which is true even when using only G23 measurements at
$z \gt 5$
. In our most extreme scenario, the tail of reionisation is absorption-dominated by a factor of
$\approx 9$
. This suggests that absorption by star-forming galaxies and/or small-scale intergalactic structure in the ionised IGM may dominate the reionisation budget, at least at reionisation’s tail end (Davies et al. Reference Davies, Bosman, Furlanetto, Becker and D’Aloisio2021; Cain et al. Reference Cain, D’Aloisio, Gangolli and Becker2021, Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024), perhaps requiring an increased photon budget to finish reionisation (Muñoz et al. Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024; Davies et al. Reference Davies, Bosman and Furlanetto2024a).
If the trend suggested by our fiducial measurement – that
$\dot{N}_{\textrm{ion}}$
increases with redshift – holds true to higher redshifts, then a
$z \sim 5.3$
end to reionisation would likely require the entire process to be absorption-dominated. Indeed, this would be necessary in a scenario like that proposed by Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024), in which galaxies produce many more ionising photons that needed to re-ionise the universe by this time. This would be consistent with the recent measurements of the IGM clumping factor by Davies et al. (Reference Davies, Bosman and Furlanetto2024a) – they find
$C \sim 12$
at
$z \leq 5$
and an upward trend towards
$z = 6$
. However, our measurement using the G23 data allows for scenarios in which
$\dot{N}_{\textrm{ion}}$
decreases at
$z \gt 6$
, in which case only the end stages may be dominated by recombinations.
4.3. Comparison to simulations
There have been several recent attempts to reproduce Ly
$\alpha$
forest and MFP measurements at
$z \leq 6$
using numerical simulations of reionisation. These include semi-numerical approaches (e.g. Qin et al. Reference Qin2024), post-processing RT (e.g. Cain et al. Reference Cain, D’Aloisio, Gangolli and Becker2021), and RT coupled to hydrodynamics and galaxy formation (e.g. Kannan et al. Reference Kannan, Garaldi, Smith, Pakmor, Springel, Vogelsberger and Hernquist2022). In some cases,
$\dot{N}_{\gamma}(z)$
is calibrated or fitted to reproduce these observations (e.g. Kulkarni et al. Reference Kulkarni, Keating, Haehnelt, Bosman, Puchwein, Chardin and Aubert2019; Cain et al. Reference Cain, D’Aloisio, Gangolli and Becker2021; Qin et al. Reference Qin2024), and in others it is predicted from an assumed galaxy model (e.g. Ocvirk et al. Reference Ocvirk, Lewis, Gillet, Chardin, Aubert, Deparis and Thélie2021; Lewis et al. Reference Lewis2022; Garaldi et al. Reference Garaldi, Kannan, Smith, Springel, Pakmor, Vogelsberger and Hernquist2022). The
$\dot{N}_{\textrm{ion}}$
in models that reproduce quasar observations sets a theoretical expectation that we can compare to our measurements.
We compare our fiducial measurement to
$\dot{N}_{\textrm{ion}}$
from numerical simulations in the top panel of Figure 3. The black solid curve and shaded region is our measurement and
$1\sigma$
range, and the faded curves are results from several recent numerical simulations in the literature (referenced in the caption). All these agree reasonably well with the mean transmission of the Ly
$\alpha$
forest and its large-scale fluctuations at
$z \lt 6$
. However, our fiducial measurement is a factor of
$\sim 2$
above the simulations at
$z \geq 5$
, with most well below the
$1\sigma$
range at all redshifts. At face value, this hints at a possible tension between simulations of re-ionisation’s end and direct
$\dot{N}_{\textrm{ion}}$
measurements.
Comparison of our fiducial measurement of
$\dot{N}_{\textrm{ion}}$
to simulations that agree with the properties of the Ly
$\alpha$
forest at
$z \lt 6$
. Top: our fiducial measurement, with
$1\sigma$
uncertainties, compared with simulation results from Kulkarni et al. (Reference Kulkarni, Keating, Haehnelt, Bosman, Puchwein, Chardin and Aubert2019), Keating et al. (Reference Keating, Kulkarni, Haehnelt, Chardin and Aubert2020), Ocvirk et al. (Reference Ocvirk, Lewis, Gillet, Chardin, Aubert, Deparis and Thélie2021), Yeh et al. (Reference Yeh2023), Gaikwad et al. (Reference Gaikwad2023), Asthana et al. (Reference Asthana, Kulkarni, Haehnelt, Bolton, Keating and Simmonds2024a), Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024), Qin et al. (Reference Qin2024). Our measurement is a factor of
$\sim 2$
above most simulation results, suggesting a possible tension between measurements and simulations. Bottom: mock measurement of
$\dot{N}_{\textrm{ion}}$
using the
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
from the late start/late end model of Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024) (green dot-dashed) run with the FlexRT code compared to the simulation result (gray dashed). We assume
$\beta_{\textrm{N}} = 1.9$
and
$\alpha = 1.5$
, consistent with the IGM and source properties in the simulation. The agreement between these validates our formalism. The red dotted curve shows the same calculation assuming
$\beta_{\textrm{N}} = 1.3$
, which lies above the simulation by a factor of
$\sim 1.5$
, potentially explaining some of the difference between the simulation and our measurement. We also show, for reference, our fiducial measurement (black solid, same as in the top panel), which assumes
$\beta_{\textrm{N}} = 1.3$
and the measured
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
.
The bottom panel of Figure 3 explores the origin of this apparent disagreement. The gray-dashed curve shows (as in the top panel)
$\dot{N}_{\textrm{ion}}$
from the late start/late end model of Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024). We then take
$\lambda_{912}^{\textrm{mfp}}$
and
$\Gamma_{\textrm{HI}}$
from the simulation and repeat the procedure described in Section 2 to get a mock ‘measurement’ of
$\dot{N}_{\gamma}$
, which is shown as the green dot-dashed curve. In this calculation, we assume
$\beta_{\textrm{N}} = 1.9$
, consistent with the typical value of
$\beta_{\textrm{N}}$
seen for the ionising opacity model used in FlexRT (see Appendix B and Figure B2 of Cain et al. Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024 for details). We also use
$\alpha = 1.5$
, the value used in the simulation. Our mock measurement agrees well with
$\dot{N}_{\textrm{ion}}$
from the simulation, validating the formalism used in this work. The red dotted curve is the same calculation, but assuming
$\beta_{\textrm{N}} = 1.3$
, corresponding to a larger contribution to the IGM opacity from high column-density absorbers. The inferred
$\dot{N}_{\textrm{ion}}$
is a factor of
$\sim 1.5$
higher than assuming
$\beta_{\textrm{N}} = 1.9$
, which is a large fraction of the difference between simulations and our fiducial measurement. For reference, we also show the fiducial measurement from the top panel as the black solid curve. The difference between this and the red dotted curve arises from (
$10-20\%$
-level) differences in
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
between the simulation and the measurement.
This comparison indicates that uncertainty in
$\beta_{\textrm{N}}$
may explain the difference between our measurement and
$\dot{N}_{\textrm{ion}}$
in simulations that reproduce the Ly
$\alpha$
forest (see also Asthana et al. Reference Asthana, Kulkarni, Haehnelt, Bolton, Keating and Simmonds2024b). If
$\beta_{\textrm{N}}$
is close to 2, as it is in FlexRT, the IGM opacity is likely dominated by low column density, highly ionised absorbers (McQuinn, Oh, & Faucher-Giguère Reference McQuinn, Peng Oh and Faucher-Giguère2011). A value closer to 1 would indicate a large contribution from high column, self-shielding absorbers, and would demand a higher ionising output from galaxies. The true column density distribution at these redshifts is poorly understood. It likely is not well-described by a single power law, evolves in a complicated way during reionisation (Nasir et al. Reference Nasir, Cain, D’Aloisio, Gangolli and McQuinn2021), depends on the dynamics of small-scale structures that are challenging to resolve in reionisation simulations (Park et al. Reference Park, Shapiro and Choi2016; D’Aloisio 2020; Chan et al. Reference Chan, Bentez-Llambay, Theuns, Frenk and Bower2024; Gnedin Reference Gnedin2024). The considerable uncertainty in measurements caused by
$\beta_{\textrm{N}}$
motivates further studies of the HI column density distribution.
4.4. Measurements of
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
We can translate our
$\dot{N}_{\textrm{ion}}$
measurements into constraints on
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
using Equation (1). Following B24, we compute
$\rho_{\textrm{UV}}(z)$
at
$2.5 \lt z \lt 6$
using the measured UVLFs from Bouwens et al. (Reference Bouwens2021), and we use two limiting UV magnitudes,
$M_{\textrm{UV}} = -17$
and
$-11$
. In the top left panel of Figure 4, we show our constraints on
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
using both cutoffs (black solid and magenta-dotted curves, respectively). Following B24, we assume that at
$z \gt 4$
, it is reasonable to treat
$\dot{N}_{\textrm{ion}}$
as dominated by the galaxy population. At
$z \lt 4$
, the ionising output of quasars likely begins to contribute significantly (Kulkarni,Worseck, & Hennawi Reference Kulkarni, Worseck and Hennawi2019; Finkelstein & Bagley Reference Finkelstein and Bagley2022; Smith et al. Reference Smith2024) and may even dominate the ionising budget (Boutsia et al. Reference Boutsia2021). As such, our estimates of
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
must be interpreted as upper limits. At
$z \gt 4$
, the shaded regions show the
$1\sigma$
uncertainties, including errors from the measurements of
$\dot{N}_{\textrm{ion}}$
and
$\rho_{\textrm{UV}}$
. For consistency with B24, we estimate the latter from the reported errors on the amplitude of the UVLF in Bouwens et al. (Reference Bouwens2021). The black and magenta points show measurements from B24 assuming the same limiting values of
$M_{\textrm{UV}}$
. At
$z \lt 4$
, we show only shaded regions denoting upper limits, as annotated in the figure.
Constraints on
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
at
$2.5 \lt z \lt 6$
. Top Left: our fiducial constraints for
$M_{\textrm{UV}} \lt -17$
(black solid curve) and
$M_{\textrm{UV}} \lt -11$
(dotted magenta curve). The shaded regions at
$z \gt 4$
indicate
$1\sigma$
uncertainties, which include errors from both
$\dot{N}_{\textrm{ion}}$
and
$\rho_{\textrm{UV}}$
measurements. At
$z \lt 4$
, we treat our constraints as strict upper limits, since AGN likely dominate the ionising output of the source population at those redshifts. As such, we show shaded regions extending down to 0 at
$z \lt 4$
. The black and magenta points are constraints from B24. Our measurements are close to those of B24 at
$z = 6$
, nearly
$1\sigma$
higher at
$z = 5$
, and at
$z = 4$
our upper limit is slightly above theirs. Top Right: the red curves show that for
$M_{\textrm{UV}} \lt -16$
, our measurement roughly agrees with the model used in Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024), which uses
$\xi_{\textrm{ion}}$
measurements from Simmonds et al. (Reference Simmonds2023) and the
$f_{\textrm{esc}}-\beta_{\textrm{UV}}$
relation from Chisholm et al. (Reference Chisholm2022) (thin red dashed curve). The blue curves show the same comparison, but using the updated fit to measurements of
$\xi_{\textrm{ion}}$
from the complete JADES sample in Simmonds et al. (Reference Simmonds2024). In this case, we find a similar level of agreement for
$M_{\textrm{UV}} \lt -13.5$
, which is more consistent with constraints on the faint-end cutoff of the UVLF (see text). Bottom Left: Same as in the top left panel, but measuring
$\dot{N}_{\textrm{ion}}$
including only G23 measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
at
$z \geq 5$
. Here,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
is nearly flat at
$z \gt 5$
for
$M_{\textrm{UV}} \lt -17$
, declines slightly for
$M_{\textrm{UV}} \lt -11$
, and is well below the B24 measurements at
$z = 6$
. Bottom Right: like the bottom left, but excluding only the G23 data. In this case,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
rises more steeply than in our fiducial measurement, and is more than
$1 \sigma$
above the
$z = 6 $
B24 measurements.
At
$z = 6$
, the central values of our measurements are similar to those of B24, but our error bars are smaller. The agreement is coincidental, since our assumed
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
are higher than theirs at this redshift. At
$z = 5$
, our measurement is almost
$1\sigma$
above that of B24, reflecting the difference between our
$\dot{N}_{\textrm{ion}}$
and that of G23 seen in Figure 2. This arises in part from our more accurate treatment of RT effects, and also from a difference in assumed
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements. Our
$1\sigma$
upper limits are slightly higher than those of B24 at
$z = 4$
. At
$z \lt 4$
, our limits evolve little with redshift until
$z = 2.5$
, reflecting the lack of evolution in
$\rho_{\textrm{UV}}$
and
$\dot{N}_{\textrm{ion}}$
. In reality,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
probably continues to decline at
$z \lt 4$
due to the increasing fraction of
$\dot{N}_{\textrm{ion}}$
sourced by AGN. Note that we have not included any correction for the presence of neutral islands in Figure 4, since we have shown that this correction is small (and depends on the uncertain reionisation history).
The top right panel compares our results to those predicted by the empirically-motivated model of Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024). They combined measurements of
$\xi_{\textrm{ion}}$
from Simmonds et al. (Reference Simmonds2023)Footnote k with the
$f_{\textrm{esc}}-\beta_{\textrm{UV}}$
relation calibrated by Chisholm et al. (Reference Chisholm2022) to estimate
$\dot{N}_{\textrm{ion}}$
during reionisation. Both the
$\xi_{\textrm{ion}}$
measurements by Simmonds et al. (Reference Simmonds2023) and
$f_{\textrm{esc}}$
predicted by the
$\beta_{\textrm{UV}}$
-
$f_{\textrm{esc}}$
increase with
$M_{\textrm{UV}}$
(that is, for fainter galaxies). As such, the prediction for
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
of the Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) model increases for fainter
$M_{\textrm{UV}}$
cutoffs, which is the opposite of the dependence of our measurements. Given this, we can ask what cutoff produces the best agreement between their model and our measurement. The red solid and thin dashed curves show that our measurement agrees best with their model for
$M_{\textrm{UV}} \lt -16$
. This agreement echos the conclusion in their work that relatively bright UV cutoffs are required to bring their model into agreement with observations supporting a late end to reionisation.
The recent re-measurements of
$\xi_{\textrm{ion}}$
presented in Simmonds et al. (Reference Simmonds2024) predict significantly lower
$\xi_{\textrm{ion}}$
, and weaker dependence on
$M_{\textrm{UV}}$
, than found in Simmonds et al. (Reference Simmonds2023). Updating the
$\xi_{\textrm{ion}}$
used in the Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) model with their new result and assuming
$M_{\textrm{UV}} \lt -13.5$
gives the thin blue dashed curve, which agrees reasonably well with our measurement assuming the same cutoff (solid blue). This indicates that the updated
$\xi_{\textrm{ion}}$
measurements from Simmonds et al. (Reference Simmonds2024) relieve some of the apparent tension between the Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) model and a late end to reionisation. This cutoff is fainter than currently probed by lensed galaxies (Atek et al. Reference Atek, Richard, Kneib and Schaerer2018), and as such does not violate existing constraints on the faint-end turnover of the UVLF at these redshifts.
At
$z \gt 5$
, our measurements imply steeper redshift evolution of
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
than in the Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) model. This is a result of the rapid decline in the MFP at
$z \gt 5$
, which pushes
$\dot{N}_{\textrm{ion}}$
up quickly approaching
$z = 6$
. We measure a factor of
$\approx 2$
increase in
$\dot{N}_{\textrm{ion}}$
between
$z = 5$
and 6, which translates into a similar increase in
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
, since
$\rho_{\textrm{UV}}$
stays roughly constant. This evolution is qualitatively similar to that required in several prior works that match the evolution of the observed mean transmission of the Ly
$\alpha$
forest (Kulkarni et al. Reference Kulkarni, Keating, Haehnelt, Bosman, Puchwein, Chardin and Aubert2019; Keating et al. Reference Keating, Kulkarni, Haehnelt, Chardin and Aubert2020; Ocvirk et al. Reference Ocvirk, Lewis, Gillet, Chardin, Aubert, Deparis and Thélie2021, see also Cain et al. Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024).
In the bottom panels, we show how our results change when we use different sub-sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements. The bottom left panel is the same as the top left, but using only the G23 measurements, while the bottom right excludes them. In the former case,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
increases by a factor of
$\approx 2$
from
$z = 2.5$
to 5, then plateaus at
$z \gt 5$
for
$M_{\textrm{UV}} \lt -17$
(and declines slightly for
$M_{\textrm{UV}} \lt -11$
). In the latter,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
rises more steeply than in our fiducial measurement, eclipsing the
$z = 6$
measurements from B24 by
$\approx 2\sigma$
. In this scenario,
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
increases by a factor of
$\sim 3-4$
(depending on the
$M_{\textrm{UV}}$
cutoff) between
$z = 5$
and 6. These different sets of measurements have very different implications for the ionising properties of
$z \gt 5$
galaxies, reflect the large uncertainty in our knowledge of galaxy ionising properties at these redshifts. Thus, converging on precise high-z measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda$
is critical for understanding the ionising output of galaxies as reionisation is ending.
4.5. Measurements of
$f_{\textrm{esc}}$
Lastly, we turn to perhaps the most uncertain parameter in reionisation models,
$f_{\textrm{esc}}$
. We can constrain the population averaged
$f_{\textrm{esc}}$
by taking advantage of the redshift and
$M_{\textrm{UV}}$
-dependent fits to measurements of
$\xi_{\textrm{ion}}$
provided in Simmonds et al. (Reference Simmonds2023) and Simmonds et al. (Reference Simmonds2024). We can combine measurements of the UVLF with this result to measure the intrinsic ionising photon production rate averaged over the galaxy population,
$\begin{equation} \dot{N}_{\textrm{ ion }}^{\textrm{intr.}} = \int_{-\infty}^{M_{\textrm{UV}}^{\textrm{cut}}} dM_{\textrm{UV}} \frac{dn}{dM_{\textrm{UV}}} L_{\textrm{UV}} \xi_{\textrm{ion}}(z,M_{\textrm{UV}})\end{equation}$
where
$\frac{dn}{dM_{\textrm{UV}}}$
is the UVLF. Then the population-averaged escape fraction is
$\begin{equation} \langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}} = \frac{\dot{N}_{\textrm{ion}}}{\dot{N}_{\textrm{ion}}^{\textrm{intr.}}}\end{equation}$
The sub-script
$\dot{n}_{\gamma}^{\textrm{intr.}}$
indicates that the average is weighted by the intrinsic ionising output of individual galaxies. A caveat is that the measured ionising efficiency, which we call
$\xi_{\textrm{ion}}^0$
, is related to the true one by
$\xi_{\textrm{ion}}^0 = \xi_{\textrm{ion}}(1-f_{\textrm{esc}})$
, since escaping ionising radiation does not contribute to the recombination emission used to measure
$\xi_{\textrm{ion}}$
. We correct for this approximately by assuming that
$\xi_{\textrm{ion}} = \xi_{\textrm{ion}}^0/(1 - \langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}})$
, which ignores any dependence of
$f_{\textrm{esc}}$
on
$M_{\textrm{UV}}$
. Under this approximation, it is straightforward to show that
$f_{\textrm{esc}} = f_{\textrm{esc}}^0/(1 + f_{\textrm{esc}}^0)$
, where
$f_{\textrm{esc}}^0$
is Equation (13) evaluated assuming
$\xi_{\textrm{ion}} = \xi_{\textrm{ion}}^0$
. Since Equation (12) uses measured galaxy properties to get
$\dot{N}_{\textrm{ion}}^{\textrm{intr.}}$
, the
$\dot{N}_{\textrm{ion}}$
that is applied in Equation 13 should be the contribution from galaxies only, without the AGN contribution. To avoid complications and uncertainties associated with quantifying the AGN contribution, we instead use our measured
$\dot{N}_{\textrm{ion}}$
and interpret our
$z \lt 4$
results as upper limits, as we do in Figure 4.
We show our fiducial results for this quantity in the top left panel of Figure 5 for
$M_{\textrm{UV}} \lt -17$
and
$-11 $
(thick curves). These measurements use the updated
$\xi_{\textrm{ion}}$
estimates from Simmonds et al. (Reference Simmonds2024). In the top-right panel, we show the same results assuming the prior findings of Simmonds et al. (Reference Simmonds2023). The black point shows the
$f_{\textrm{esc}} = 0.085$
measured at
$z \sim 3$
by Pahl et al. (Reference Pahl, Shapley, Steidel, Chen and Reddy2021) for their sample of faint (
$L \lt L_{\ast}$
) galaxies, which are more likely to be analogous to the galaxies that dominated the ionising photon budget during reionisation (Atek et al. Reference Atek2024, although see arguments to the contrary in Naidu et al. Reference Naidu2022; Matthee et al. Reference Matthee2022). We include uncertainties in the power-law intercept reported in the
$\xi_{\textrm{ion}}$
fits in Simmonds et al. (Reference Simmonds2023, Reference Simmonds2024) in our error budgetFootnote
l
on
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
, along with the errors on
$\dot{N}_{\textrm{ion}}$
and
$\rho_{\textrm{UV}}$
already included in Figure 4.
Constraints on the
$\dot{n}_{\gamma}^{\textrm{intr.}}$
-weighted average escape
$f_{\textrm{esc}}$
of the galaxy population at
$2.5 \lt z \lt 6$
. Top left: Fiducial constraints for
$M_{\textrm{UV}} \lt -17$
and
$-11$
. At
$z \gt 4$
we report our constraints as measurements, and as upper limits at
$z \lt 4$
(as in Figure 4). The black point shows
$f_{\textrm{esc}} = 0.085$
measured by Pahl et al. (Reference Pahl, Shapley, Steidel, Chen and Reddy2021) at
$z \sim 3$
for faint (
$L \lt L_{\ast}$
) galaxies. Top right: the same, but using the prior results for
$\xi_{\textrm{ion}}$
from Simmonds et al. (Reference Simmonds2023). Bottom Left: same as the top left, but using only
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements from G23 at
$z \gt 5$
. Bottom Right: the same, but excluding the G23 results. See text for details.
The redshift evolution we infer for
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
is qualitatively similar to that of
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
. At
$z \lt 4$
and
$M_{\textrm{UV}} \lt -17$
, our upper limit is not much higher than the Pahl et al. (Reference Pahl, Shapley, Steidel, Chen and Reddy2021) measurement of
$8.5\%$
at
$z = 3$
, and is actually slightly below this value for
$M_{\textrm{UV}} \lt -11$
. Between
$z = 4$
and 5, our measurements increase by a factor of
$\approx 2$
, and by another factor of
$\approx 2$
from
$z = 5$
to 6. This steep increase occurs because
$\rho_{\textrm{UV}}$
decreases between
$z = 4$
and 6 and
$\xi_{\textrm{ion}}$
remains flat, such that an increasing
$f_{\textrm{esc}}$
is required to explain the increase in
$\dot{N}_{\textrm{ion}}$
over this range. Note that because the AGN contribution to
$\dot{N}_{\textrm{ion}}$
grows between
$z = 4$
and
$2.5$
,
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
likely declines with cosmic time over this range.
The top right panel shows the same measurements, but using the older results for
$\xi_{\textrm{ion}}$
from Simmonds et al. (Reference Simmonds2023). Note that the error bars are much larger because of the higher uncertainties on
$\xi_{\textrm{ion}}$
in that work. We find systematically lower
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
, reflecting the higher
$\xi_{\textrm{ion}}$
values measured in that work. We also see noticeably different redshift evolution, especially for the
$M_{\textrm{UV}} \lt -11$
case. The evolution of
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
is significantly flatter, and does not reach 10% until
$z = 6$
in that case. The dependence on the assumed
$M_{\textrm{UV}}$
cutoff is also much stronger, with a factor of
$\approx 2$
difference between
$M_{\textrm{UV}} \lt -11$
and
$M_{\textrm{UV}} \lt -17$
, compared to a 15–20% difference in the upper left. This arises from the fact that in the Simmonds et al. (Reference Simmonds2023) results, fainter galaxies have much higher
$\xi_{\textrm{ion}}$
than bright galaxies. In the updated measurements,
$\xi_{\textrm{ion}}$
does not depend strongly on
$M_{\textrm{UV}}$
. This finding further highlights the much lower ionising output expected for galaxies based on the new Simmonds et al. (Reference Simmonds2024), which helps relieve the “too many photons” problem described in Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024).
The bottom left is the same as the top left, but using only the G23 measurements at
$z \gt 5$
to measure
$\dot{N}_{\textrm{ion}}$
. Consistent with Figure 4, we find a more modest evolution in
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
than our fiducial result. We find
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}} \lt 15\%$
at
$z = 2.5$
and 4,
$\approx 15\%$
at
$z = 5$
, and
$\approx 19\%$
at
$z = 6$
for
$M_{\textrm{UV}} \lt -17$
. These numbers become 12%, 11%, and 9% for
$M_{\textrm{UV}} \lt -11$
. Within the uncertainties, these results are consistent with no evolution in
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
at
$2.5 \lt z \lt 6$
. In the bottom right, we show the result without G23 measurements included. We find steep evolution in
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
at
$z \gt 4$
, rising to 40–60% by
$z = 6$
(depending on the cutoff
$M_{\textrm{UV}}$
).Footnote
m
Our upper limits are consistent with population-averaged escape fractions of
$\lesssim 10\%$
at
$z \leq 4$
. This agrees with recent efforts to measure
$f_{\textrm{esc}}$
directly at these redshifts (e.g. Smith et al. Reference Smith2018, Reference Smith2020; Pahl et al. Reference Pahl, Shapley, Steidel, Chen and Reddy2021; Kerutt et al. Reference Kerutt2024) and with indirect determinations based on simulations (Finkelstein et al. Reference Finkelstein2019; Yeh et al. Reference Yeh2023; Choustikov et al. Reference Choustikov2024). Indeed, escape fractions at these redshifts may be well below 10% if the ionising background is sustained mainly by the quasar population (Boutsia et al. Reference Boutsia2021). In this case, a self-consistent measurement of galaxy ionising properties at
$z \lt 4$
should also take into a account the luminosity function and ionising properties of the quasar population. Such an endeavour could be aided by using the properties of the HeII Ly
$\alpha$
forest to distinguish between galaxy and quasar ionising output (e.g. McQuinn &Worseck Reference McQuinn and Worseck2014; D’Aloisio et al. Reference D’Aloisio, Upton Sanderbeck, McQuinn, Trac and Shapiro2017; Gaikwad, Davies, & Haehnelt Reference Gaikwad, Davies and Haehnelt2025). We defer such an investigation to future work.
At
$z \gt 4$
, when
$f_{\textrm{esc}}$
becomes impossible to measure directly, we find a factor of
$\sim 1.5-2$
increase between
$z = 4$
and 5 in our fiducial measurement. Indeed, the actual evolution is probably steeper than this, because AGN contribute more to
$\dot{N}_{\textrm{ion}}$
at
$z = 4$
than at
$z = 5$
. This finding suggests modest but significant evolution in galaxy ionising properties, perhaps due to evolution in the properties of the ISM/CGM between
$z = 4$
and 5 (e.g. Kakiichi & Gronke Reference Kakiichi and Gronke2021; Kimm et al. Reference Kimm, Bieri, Geen, Rosdahl, Blaizot, Michel-Dansac and Garel2022). At
$z \gt 5$
, our results are sensitive to the choice of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements used. Using the only indirect measurements of G23 suggests flat evolution in
$f_{\textrm{esc}}$
at
$z \gt 5$
, while ignoring these measurements gives a factor of
$\sim 3$
increase between
$z = 5$
and 6. These findings motivate further efforts to reduce uncertainty on measurements of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
at
$z \gt 5$
, and to understand why different measurement methods give significantly different results.
4.6. Comparison to indirect determinations of
$f_{\textrm{esc}}$
We briefly compare our
$f_{\textrm{esc}}$
results to several indirect observational and theoretical determinations. These include empirically motivated estimates of
$f_{\textrm{esc}}$
based on measurements at lower redshift and estimates from numerical simulations. In Figure 5, we show the comparison to our maximum likelihood measurements of
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
for all three sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements with
$M_{\textrm{UV}} \lt -11$
. The thin dot-dashed green curve shows the empirical model of Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) (for
$M_{\textrm{UV}} \lt -11$
), based on the
$\beta_{\textrm{UV}}$
-
$f_{\textrm{esc}}$
relation from Chisholm et al. (Reference Chisholm2022), and the
$\beta_{\textrm{UV}}-M_{\textrm{UV}}$
relation from Zhao & Furlanetto (Reference Zhao and Furlanetto2024). Unlike Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024), however, we use the updated
$\xi_{\textrm{ion}}$
measurements from Simmonds et al. (Reference Simmonds2024). We also show the population-averaged
$f_{\textrm{esc}}$
from THESAN red dashed, Yeh et al. Reference Yeh2023) and SPHINX blue dotted, Rosdahl et al. Reference Rosdahl2022), and the global
$f_{\textrm{esc}}$
from Finkelstein et al. (Reference Finkelstein2019) for all galaxies (black dot-dashed) and faint (
$M_{\textrm{UV}} \gt -15$
) galaxies (cyan dot-dashed). Note that we only show these down to
$z = 4$
, since our constraints at
$z \lt 4$
are upper limits.
All the simulation results lie below the measurements to varying degrees at
$z \gt 4$
, although their redshift evolution is comparable to our fiducial result. The disagreement is greatest for the measurement without the G23 data, which implies steeper redshift evolution than in the
$f_{\textrm{esc}}-\beta_{\textrm{UV}}$
model, or any of the simulation curves. Note that using a brighter cutoff
$M_{\textrm{UV}}$
increases
$f_{\textrm{esc}}$
, worsening this disagreement. At face-value, this suggests that simulations may be under-estimating galaxy escape fractions. Another possibility is that AGN contribute significantly to
$\dot{N}_{\textrm{ion}}$
at these redshifts (Madau et al. Reference Madau, Giallongo, Grazian and Haardt2024; Smith et al. Reference Smith2024; Dayal et al. Reference Dayal2024, although see Jiang et al. Reference Jiang, Jiang, Sun, Liu and Fu2025), or current measurements of
$\xi_{\textrm{ion}}$
are under-estimates. Both scenarios would reduce our measured
$f_{\textrm{esc}}$
.
The empirical
$f_{\textrm{esc}}-\beta_{\textrm{UV}}$
relation agrees reasonably well with our measurement for
$M_{\textrm{UV}} \lt -11$
, predicting
$f_{\textrm{esc}} \sim 10\%$
at
$z\lesssim 5$
(see also upper right panel of Figure 4). At
$5 \lt z \lt 6$
it agrees best with the fiducial measurement, while the measurements only with (without) data from G23 at
$z \gt 5$
fall slightly below (above) the model. One caveat to this agreement is that the measured
$f_{\textrm{esc}}$
becomes significantly higher for
$M_{\textrm{UV}} \lt -17$
(Figure 5), while the same is not true of the Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) model. As such, our fiducial measurement would lie above the prediction of the
$\beta_{\textrm{UV}}$
-
$f_{\textrm{esc}}$
relation at
$z \gt 5$
if the UVLF cuts off significantly brighter than
$M_{\textrm{UV}} = -11$
.
Recently, Qin et al. (Reference Qin2024) used semi-numerical simulations of reionisation to constrain both
$\dot{N}_{\textrm{ion}}$
(red solid curve in Figure 3) and the dependence of the escape fraction on both redshift and galaxy properties. Their constraints included UVLF measurements, the CMB optical depth, and Ly
$\alpha$
forest data, although they found the latter to be by far the most constraining on galaxy ionising properties. They found that the escape fraction is much higher in lower-mass halos, with a normalisation that increases with decreasing redshift,
$f_{\textrm{esc}} \propto (1+z)^{-1.6}$
(note that they do not plot
$\langle f_{\textrm{esc}}\rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
, so we cannot show it on Figure 6). The latter is the opposite the trend we find - that
$f_{\textrm{esc}}$
grows towards higher redshifts, especially from
$z = 5-6$
(although the strong dependence on halo mass may make the
$\dot{n}_{\textrm{intr.}}$
-weighted average evolve less quickly). This is in part due to the fact that our
$\dot{N}_{\textrm{ion}}$
increases rapidly with redshift across this redshift range while theirs remains flat (Figure 3). It is also unclear whether the
$\xi_{\textrm{ion}}$
assumed in their model matches the parameterisation we use here based on JWST measurements, which could further explain the difference between our results.
Comparison of our
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
measurements to several indirect determinations. The thin green dot-dashed curve shows the model used in Muñoz et al. (Reference Muñoz, Mirocha, Chisholm, Furlanetto and Mason2024) based off the
$f_{\textrm{esc}}$
-
$\beta_{\textrm{UV}}$
relation calibrated by Chisholm et al. (Reference Chisholm2022) at lower redshifts, and using the latest
$\xi_{\textrm{ion}}$
measurements from Simmonds et al. (Reference Simmonds2024). The other curves show simulation results from Finkelstein et al. (Reference Finkelstein2019), Rosdahl et al. (Reference Rosdahl2022), Yeh et al. (Reference Yeh2023). See text for discussion.
5. Conclusions
In this work, we have provided new measurements of the ionising emissivity of the galaxy population at
$2.5 \lt z \lt 6$
. Our measurements take into account the effects of RT without using the local source approximation. We also developed a formalism to account for the possible presence of neutral islands at
$z \gt 5$
, when reionisation may be ongoing. We present measurements of the global ionising emissivity, the average escaping ionising efficiency, and the average escape fraction of the high-redshift galaxy population. Our main conclusions are summarised below:
-
• We measure
$\dot{N}_{\textrm{ion}} = 10.04_{-4.66}^{+2.72}$
,
$7.05_{-2.17}^{+1.41}$
,
$8.70_{-2.55}^{+1.15}$
, and
$17.67_{-8.29}^{+5.62}$
$\times 10^{50}$
s
$^{-1}$
cMpc
$^{-3}$
at
$z = 2.5$
, 4, 5, and 6, respectively, at
$1\sigma$
confidence. Our measurement of
$\dot{N}_{\textrm{ion}}$
at
$z = 6$
is nearly
$2\sigma$
higher than that of G23, largely due to the shorter
$\lambda_{912}^{\textrm{mfp}}$
assumed here. At
$z = 5$
, our measurement is higher than theirs by
$\approx 0.8\sigma$
. The latter is in large part due to our inclusion of RT effects. At
$z \leq 4.5$
, our measurements are consistent with those of Becker & Bolton (Reference Becker and Bolton2013). We find that including the opacity from neutral islands at a level consistent with expectations from reionisation simulations has a
$\lesssim 15$
% effect on these measurements. The change is modest because the effect of islands on
$\dot{N}_{\textrm{abs}}$
and
$\lambda$
approximately cancels in Equation (6). This suggests that the standard approach to measuring
$\dot{N}_{\textrm{ion}}$
may work even while reionisation is ongoing. -
• We measure the
$5.3 \lt z \lt 6$
photon budget, likely covering the last 10–20% of reionisation, to be
$N_{\gamma/H}^{z = 5.3-6} = 1.10_{-0.39}^{+0.21}$
. This is a factor of
$\approx 4$
higher than the budget required to complete the last 20% of reionisation in the absence of absorptions in the ionised IGM. This elevated photon budget is reflected in our measurements of
$f_{\textrm{esc}}$
at
$z \gt 5$
, which are a factor of
$\gtrsim 2$
higher than predicted by simulations. Our findings suggest that at least the tail-end of reionisation may have been absorption-dominated. -
• We find a factor
$\sim 2$
tension between our fiducial measurement of
$\dot{N}_{\textrm{ion}}$
and the emissivity required in reionisation simulations that reproduce the observed transmission of the
$z \lt 6$
Ly
$\alpha$
forest. A significant fraction (roughly half) of that difference may arise from differences in the shape of the HI column density distribution in the simulations compared to that assumed in the measurement. -
• We constrain the collective escaping ionising efficiency of the galaxy population, assuming the UVLF from Bouwens et al. (Reference Bouwens2021), to be
$\log_{10}(\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}/[{\textrm {erg}}^{-1}{\textrm {Hz}}]) \lt 24.61$
,
$\lt 24.49$
,
$= 24.67_{-0.17}^{+0.10}$
and
$= 25.12_{-0.28}^{+0.15}$
at
$z = 2.5$
, 4, 5, and 6, respectively, for
$M_{\textrm{UV}} \lt -17$
and at
$1\sigma$
confidence. These numbers become
$\lt 24.49$
,
$\lt 24.32$
,
$= 24.47_{-0.17}^{+0.10}$
, and
$= 24.75_{-0.28}^{+0.15}$
, respectively, for
$M_{\textrm{UV}} \lt -11$
. Our measurement of
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
is similar to that of B24 at
$z = 6$
, but is almost
$1\sigma$
higher at
$z = 5$
. At
$z = 4$
, our upper limit is slightly higher than theirs. We find that our upper limits on
$\langle f_{\textrm{esc}}\xi_{\textrm{ion}}\rangle_{L_{\textrm{UV}}}$
remain roughly flat at
$2.5 \lt z \lt 4$
, but we find a factor of
$1.5-2$
increase in escaping efficiency between
$z = 4$
and 5, and another factor of 2 between
$z = 5$
and 6. The evolution between
$z = 5$
and 6 owes to the fact that the measured
$\lambda_{912}^{\textrm{mfp}}$
declines at those redshifts more quickly than does
$\Gamma_{\textrm{HI}}$
. -
• We constrain the population-averaged ionising escape fraction
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}}$
, leveraging the redshift and
$M_{\textrm{UV}}$
-dependent fit to measurements of
$\xi_{\textrm{ion}}$
provided by Simmonds et al. (Reference Simmonds2023, Reference Simmonds2024). For
$M_{\textrm{UV}} \lt -17$
, we measure
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}} \lt 0.172$
,
$\lt 0.131$
,
$= 0.178_{-0.058}^{+0.048}$
, and
$= 0.385_{-0.186}^{+0.168}$
at
$z = 2.5$
, 4, 5, and 6, respectively, at
$1\sigma$
confidence. These numbers become
$\langle f_{\textrm{esc}} \rangle_{\dot{n}_{\gamma}^{\textrm{intr.}}} \lt 0.138 $
,
$\lt 0.096$
,
$= 0.126_{-0.041}^{+0.034}$
, and
$= 0.224_{-0.108}^{+0.098}$
for
$M_{\textrm{UV}} \lt -11$
. At
$z \lesssim 4$
, our upper limits are consistent with measurements of
$f_{\textrm{esc}}$
in galaxies at those redshifts suggesting values
$\lesssim 10\%$
. However, we find a factor of
$2-3$
increase between
$z = 4$
and 6, suggesting evolution in the ionising properties of the galaxy population at these redshifts. -
• We repeated our measurements using two sub-sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements. The first excludes all measurements at
$z \gt 5$
except those of G23, and the second set includes all measurements except those of G23. The evolution of
$\Gamma_{\textrm{HI}}$
(
$\lambda_{912}^{\textrm{mfp}}$
) with redshift is steeper (shallower) with redshift in the G23-only case, resulting in much flatter
$\dot{N}_{\textrm{ion}}$
evolution. Conversely, leaving out the G23 measurements results in a steeper evolution in
$\dot{N}_{\textrm{ion}}$
. These differences translate into very different estimates of the redshift evolution of
$f_{\textrm{esc}}$
at
$z \gt 5$
- nearly flat in the former case, and a factor of 4 increase to
$\approx 60\%$
(
$40 \%$
), for
$M_{\textrm{UV}} \lt -17$
(
$-11$
), in the latter.
We find that RT has a modest, but significant effect on measurements of the ionising output of galaxies at
$4 \lt z \lt 6$
. They also highlight uncertainties in measurements of IGM conditions at
$z \gt 5$
, which translate into large uncertainties in the evolution of galaxies ionising properties near reionisation’s end. Perhaps most notably, the implied evolution of
$f_{\textrm{esc}}$
varies considerably at
$z \gt 5$
when different
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
data sets are used to calculate
$\dot{N}_{\textrm{ion}}$
. These findings motivate further efforts to pin down these measurements at
$z \gt 5$
.
Acknowledgments
The authors thank Sarah Bosman, Simeon Bird, and Rolf Jansen for helpful discussions. We also thank George Becker, Frederick Davies, Steven Finkelstein, Brent Smith, Yongda Zhu, and Rogier Windhorst for helpful comments on the draft version of this manuscript.
Data availability statement
The data underlying this article will be shared upon reasonable request to the corresponding author.
Funding statement
CC was supported by the Beus Center for Cosmic Foundations. AD was supported by grants NSF AST-2045600 and JWSTAR02608.001-A. JBM was supported by NSF Grants AST-2307354 and AST-2408637, and by the NSF-Simons AI Institute for Cosmic Origins.
Competing interests
The authors are not aware of any competing interests connected to this work.
Ethical standards
The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
Appendix A. Connection to the standard radiative transfer formalism
In this appendix, we show that our RT formulation (Equations 3–5) is equivalent to the standard formalism presented e.g. in Haardt & Madau (Reference Haardt and Madau1996). The angle-averaged specific intensity at redshift z and frequency
$\nu$
is
$\begin{equation} J_{\nu}(z) = \frac{c}{4 \pi} \int_{z}^{\infty} \frac{dz'}{(1+z')H(z')}\frac{(1+z)^3}{(1+z')^3} \epsilon_{\nu'}(\nu',z') e^{-\tau_{\textrm{eff}}(\nu,z,z')},\end{equation}$
where primed and un-primed variables have the same meaning as in Section 2. Here,
$\epsilon_{\nu'}(\nu',z')$
is the proper energy emissivity of sources per unit
$\nu'$
, and
$\tau_{\textrm{eff}}$
is the effective optical depth between z’ and z, given by (see e.g. Equation 8 of McQuinn Reference McQuinn2016),
$\begin{equation} \tau_{\textrm{eff}}(\nu,z,z') = \int_{x(z')}^{x(z)}\frac{dx}{1+z(x)} \lambda_{\nu''}^{-1} = \int_{t'}^{t} dt'' c \kappa_{\nu''},\end{equation}$
where the first integral is over co-moving distance x.Footnote
n
Note that
$\nu'' = \nu \frac{1+z''}{1+z}$
, analogously to the singly primed variable. In the second equality, we have used the fact that
$c dt = \frac{dx}{1+z(x)}$
and
$\kappa_{\nu} \equiv \lambda_{\nu}^{-1}$
, which recovers the expression in the exponential of Equation (5). Finally, the ionisation rate is given by (Equation 10 of Haardt & Madau Reference Haardt and Madau1996)
$\begin{equation} \Gamma_{\textrm{HI}} = 4 \pi \int d\nu \frac{J_{\nu}}{h_{\textrm{p}} \nu} \sigma_{\textrm{HI}}^{\nu}\end{equation}$
To show that Equations (A1)–(A4) are equivalent to Equations (3)–(5), we first change variables to t in Equation (A1), recognising that
$dt = -\frac{dz}{(1+z)H(z)}$
, and substitute Equation (A2) into Equation (A1), which gives
$\begin{equation} J_{\nu}(z) = \frac{c (1+z)^3}{4 \pi} \int_{0}^{t} dt'\frac{\epsilon_{\nu'}(\nu',t')}{(1+z')^3} e^{-\int_{t'}^{t} dt'' c \kappa_{\nu''}}\end{equation}$
The proper emissivity
$\epsilon_{\nu'}(\nu',t')$
can be re-written in terms of the co-moving photon emissivity per unit observed frequency
$\nu'$
(see Equation 4) as
$\begin{equation} \epsilon_{\nu'}(\nu',t') = (1+z')^3 \frac{1+z}{1+z'} h_{\textrm{p}} \nu' \dot{N}_{\textrm{ion}}^{\nu}(t',\nu')\end{equation}$
where the factor of
$\frac{1+z'}{1+z}$
arises from the fact that
$\dot{N}_{\textrm{ion}}^{\nu}$
is a derivative with respect to co-moving frequency
$\nu$
rather than
$\nu'$
. Putting this into Equation (A4) gives
$\begin{equation} J_{\nu}(z) = \frac{c (1+z)^3}{4 \pi} \int_{0}^{t} dt' \frac{1+z}{1+z'} h_{\textrm{p}} \nu' \dot{N}_{\textrm{ion}}^{\nu}(t',\nu') e^{-\int_{t'}^{t} dt'' c \kappa_{\nu''}}\end{equation}$
Recognising that
$\nu = \frac{1+z}{1+z'}\nu'$
and putting Equation (A6) into Equation (A3) yields
$\begin{equation} \Gamma_{\textrm{HI}} = 4 \pi \int \frac{d\nu}{h_{\textrm{p}} \nu} \sigma_{\textrm{HI}}^{\nu} \frac{c (1+z)^3}{4 \pi} \int_{0}^{t} dt' h_{\textrm{p}} \nu \dot{N}_{\textrm{ion}}^{\nu}(t',\nu') e^{-\int_{t'}^{t} dt'' c \kappa_{\nu''}}\end{equation}$
which simplifies to
$\begin{equation} \Gamma_{\textrm{HI}} = c(1+z)^3 \int d\nu \sigma_{\textrm{HI}}^{\nu} \int_{0}^{t} dt' \dot{N}_{\textrm{ion}}^{\nu}(t',\nu') e^{-\int_{t'}^{t} dt'' c \kappa_{\nu''}}\end{equation}$
Finally, we can multiply inside the integral over
$\nu$
by
$\lambda_{\nu}(t) \kappa_{\nu}(t) = 1$
and group terms to obtain
$\begin{equation} \Gamma_{\textrm{HI}} = (1+z)^3 \int d\nu \lambda_{\nu} \sigma_{\textrm{HI}}^{\nu} \int_{0}^{t} dt' \dot{N}_{\textrm{ion}}^{\nu}(t',\nu') \left[c \kappa_{\nu} e^{-\int_{t'}^{t} dt'' c \kappa_{\nu''}}\right]\end{equation}$
The term in brackets is
$G(t,t',\nu,\nu')$
(Equation 5). It is straightforward to see that Equation (A10) is the result of directly substituting Equation (5) into Equation (4), and then putting the result into Equation (3).
Appendix B. Insensitivity of
$\dot{N}_{\textrm{ion}}$
measurements to neutral islands
In Figure B1, we show the ratios of the black dotted, red dashed, and cyan dot-dashed curves in Figure 2 with the maximum likelihood fiducial measurement (black solid curve). We see that during reionisation, accounting for neutral islands (cyan curve) changes the measurement by at most 15% while reionisation is ongoing, and only at the couple-percent level afterwards. This is smaller than the error incurred by neglecting redshifting (
$\approx 20\%$
) or using the LSA (a factor of 2) at
$z = 2.5$
. This is despite the fact that our assumed reionisation scenario is somewhat extreme, such that we likely over-estimate the importance of islands for the measurement at
$5 \lt z \lt 6$
.
Relative effect on our measurement of neglecting redshifting (black dotted curve), using the LSA (red dashed curve) and accounting for neutral islands at
$z \gt 5$
(cyan dot-dashed curve). The effect of accounting for islands is smaller during reionisation than that of the other two effects at
$z \approx 2.5$
, even for the somewhat extreme reionisation scenario assumed here.
We can understand the reason for this result by considering the physical meaning of the MFP. In Appendix C of Cain et al. (Reference Cain, D’Aloisio, Lopez, Gangolli and Roth2024), we showed that the frequency-averaged MFP is well-approximated by
$\begin{equation} \lambda^{-1} = \frac{\langle \Gamma_{\textrm{HI}} n_{\textrm{HI}} \rangle_{\textrm{V}}}{F_{\gamma}}\end{equation}$
where the numerator is the volume-averaged absorption rate in the IGM and
$F_{\gamma}$
is the average incident ionising photon flux. The former is simply
$\dot{N}_{\textrm{abs}}$
, and the latter is given by
$F_{\gamma} = N_{\gamma} c$
, where
$N_{\gamma}$
is the mean number density of ionising photons. So, we can write
$\begin{equation} \lambda^{-1} = \frac{\dot{N}_{\textrm{abs}}}{N_{\gamma} c}\end{equation}$
That is, the IGM absorption coefficient is the absorption rate divided by the ionising flux. It is straightforward to define ‘ionised’ and ‘neutral’ components of Equation (B2) such that
$\begin{equation} \lambda^{-1} = \lambda_{\textrm{ionised}}^{-1} + \lambda_{\textrm{neutral}}^{-1} = \frac{\dot{N}_{\textrm{abs,ionised}}}{N_{\gamma} c} + \frac{\dot{N}_{\textrm{abs,neutral}}}{N_{\gamma} c}\end{equation}$
Indeed, the first equality is just Equation (7) and the second is Equation (8). It follows that treating the IGM absorption as one component (Equation B1) or two components (ionised and neutral, Equation B2) should give the same result, provided
$\lambda$
is appropriately defined. Of course, it is unclear whether the quantity recovered by the observations precisely matches the ‘theoretical’ definition of the MFP [Roth et al.(Reference Roth, D’Aloisio, Cain, Wilson, Zhu and Becker2024)]. Still, this shows the physical reason why we should not expect treating absorption by neutral islands separately to significantly affect our measurement.
In the top panel of Figure B2, we break down the contribution of the different components of the IGM (ionised and neutral) to the total absorption rate,
$\dot{N}_{\textrm{abs}}$
(Equations 4 and 8), in our fiducial measurement. The black solid curve shows
$\dot{N}_{\textrm{abs}}$
without accounting for neutral islands (Equation 4), and the cyan dashed curve shows the same accounting for islands (Equation 8). The red dotted and green dot-dashed curves show contributions to
$\dot{N}_{\textrm{abs}}$
in this case arising from only ionised gas (
$\dot{N}_{\textrm{abs,ionised}}$
), and from only neutral islands (
$\dot{N}_{\textrm{abs,neutral}}$
), respectively. Although the measured absorption by ionised gas when including the effects of islands is smaller than it is without including them, the additional absorption by the islands almost offsets this difference, yielding nearly the same total
$\dot{N}_{\textrm{abs}}$
. Another interesting result is that
$\dot{N}_{\textrm{abs,neutral}}$
is always less than half of
$\dot{N}_{\textrm{abs,ionised}}$
, even at
$z = 6$
when the IGM is 30% neutral in our assumed reionisation model. Our fiducial measurement therefore suggests that the majority of absorption during reionisation’s tail end takes place in the highly ionised IGM - the so-called ‘absorption-dominated’ scenario (Davies et al. Reference Davies, Bosman, Furlanetto, Becker and D’Aloisio2021, Reference Davies, Bosman and Furlanetto2024a).
Contribution to the total absorption rate by the ionised IGM and neutral islands. Top: The black solid curve shows
$\dot{N}_{\textrm{abs}}$
for our fiducial measurement without accounting for neutral islands at
$z \gt 5$
, and the cyan dashed curve shows the same accounting for islands. The red dotted and green dot-dashed curves show the absorption rate by only ionised, and only neutral gas, respectively. Bottom: the same, but for our G23-only measurement. In this case,
$\dot{N}_{\textrm{ion,ionised}}$
and
$\dot{N}_{\textrm{ion,neutral}}$
are comparable at
$z = 6$
, but still add to a value close to the black solid curve. This shows that the standard formalism works well even if absorption by ionised gas does not dominate the total budget.
In the bottom panel, we show the same thing, but for our G23-only measurements. In this case,
$\dot{N}_{\textrm{abs,ionised}} \approx \dot{N}_{\textrm{ion,neutral}}$
at
$z \geq 5.8$
, such that reionisation is not absorption-dominated. However, even in this case, the total is very close to the black solid curve. This suggests that even if the absorption rate is not dominated by ionised gas, the standard formalism remains a good approximation. This is good news for measurements of
$\dot{N}_{\textrm{ion}}$
using this technique, since the reionisation history at
$5 \lt z \lt 6$
is poorly constrained observationally.
Appendix C. Fitting functions for
$\boldsymbol{\Gamma}_{\textbf{HI}}$
and
$\boldsymbol{\lambda}_{\textbf{912}}^{\textbf{mfp}}$
Here, we give the maximum likelihood fits to the sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements used in our analysis. Our fiducial result (including all measurements) is shown in Figure 1. We show fits to the two subsets described in the main text in Figure C1, in the same format as Figure 1. The left column shows our fit using only G23 data at
$z \gt 5$
, and the left panel excludes the G23 at these redshifts. In the first case,
$\Gamma_{\textrm{HI}}$
decreases with redshift more rapidly than in the fiducial fit, and
$\lambda_{912}^{\textrm{mfp}}$
declines less rapidly. The opposite is true in the right column. In this case, the
$z \gt 5$
evolution of
$\lambda_{912}^{\textrm{mfp}}$
is determined entirely by the direct measurements from Zhu et al. (Reference Zhu2023).
Fits to alternative sets of
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements. Left column: fit using only G23 measurements at
$z \gt 5$
, in the same format as Figure 1. Right column: the same, but excluding only the G23 points at
$z \gt 5$
. In the left column, we see faster (slower) redshift evolution in
$\Gamma_{\textrm{HI}}$
(
$\lambda_{912}^{\textrm{mfp}}$
) than in the fiducial case, and the opposite is true in the right column. This has a significant effect on inferred galaxy ionising properties at
$z \gt 5$
, as shown in the main text.
We fit the
$\Gamma_{\textrm{HI}}$
measurements, in units of s
$^{-1}$
, to a fifth-order polynomial of the form
$\begin{equation} \log_{10}(\Gamma_{\textrm{HI}}/[\text{s}^{-1}]) = \sum_{n = 0}^{4} a_n z^n\end{equation}$
where
$(a_0,a_1,a_2,a_3,a_4) = ({-}0.47,-11.50,4.09,-0.62,0.03)$
are the maximum likelihood parameters for our fiducial fit. Using only G23 data at
$z \gt 5$
, we obtain
$(a_0,a_1,a_2,a_3,a_4) =$
$ ({-}6.52,-5.08,1.62,-0.20,0.01)$
, and excluding the G23 measurements,
$(a_0,a_1,a_2,a_3,a_4) = ({-}5.18,-6.60,2.24,-0.32,0.02)$
.
Our fit to
$\lambda_{912}^{\textrm{mfp}}$
measurements, in units of
$h^{-1}$
cMpc, has the functional form
$\begin{equation} \log_{10}\left([\lambda_{912}^{\textrm{mfp}} + 1]/[h^{-1}\text{cMpc}]\right) = \frac{a_1 z^{b_1}}{1 + (z/a_2)^{b_2}}\end{equation}$
where
$(a_1,b_1,a_2,b_2) = (3.95, -0.46, 5.92, 13.90)$
is the maximum likelihood result for our fiducial fit. Using only G23 measurements at
$z \gt 5$
, we find
$(a_1,b_1,a_2,b_2) = (3.98,-0.46,6.08,13.38)$
, and without those measurements, we find
$(a_1,b_1,a_2,b_2) = (3.98,-0.46,5.77,16.77)$
. We emphasise that these fitting functions should not be extrapolated much beyond
$2.5 \lt z \lt 6$
.
Appendix D. Error budget for
$\dot{N}_{\textrm{ion}}$
Here, we take a closer look at the breakdown of our error budget, and compare our errors more closely to those of previous works that measure
$\dot{N}_{\textrm{ion}}$
. Figure D1 shows our fiducial measurement with errors from the top panel of Figure 2, and compares this to the errors in Becker & Bolton (Reference Becker and Bolton2013), which are several times larger than ours. In that work, the authors combined statistical and systematic errors from a number of sources to measure
$\dot{N}_{\textrm{ion}}$
, including uncertainty in
$\alpha$
and
$\beta_{\textrm{N}}$
as we do in this work (see their Table 3). They combined uncertainties from difference sources by linearly adding the logarithmic errors from each component, which assumes maximal correlation between sources of error and is the most conservative assumption. Our analysis assumes independence between different sources of error. If we make the same assumption using the reported errors in Becker & Bolton (Reference Becker and Bolton2013), and add their fractional errors in quadrature, we recover the magenta points, which have
$\pm 1\sigma$
ranges about half as large as their reported errors.
A more careful comparison of our error bars to those of Becker & Bolton (Reference Becker and Bolton2013) at
$z \lt 5$
. We show our fiducial measurement and error bars from the top panel of Figure 2, alongside the reported Becker & Bolton (Reference Becker and Bolton2013) measurements and error bars (green points) and two revisions of their errors. The purple points show what happens if the error from different sources in their analysis are combined in quadrature rather than linearly, which is more consistent with the assumptions made in our work. The red points further adjust the size of some of their systematic errors to better reflect those in our work (see text). With these revisions, their errors are only
$\approx 30\%$
larger than ours – the rest of the difference is likely explained by lack of flexibility in our parametric fits to
$\Gamma_{\textrm{HI}}$
and
$\lambda_{912}^{\textrm{mfp}}$
measurements.
There are a couple additional differences between our analysis and that of Becker & Bolton (Reference Becker and Bolton2013) can explain some of the remaining difference between our error bars and theirs. First, they do not assume a 4 Ryd cutoff in their ionising spectrum, which results in their uncertainties from
$\alpha$
being larger than ours. They also include a systematic uncertainty for the effects of recombination radiation that is not included in our analysis. If we further adjust their errors to reflect these differences, we get the red points, which have uncertainties only
$\approx 30\%$
larger than ours. The remaining difference is likely explained by the fact that we fit a parametric function to a large number of
$\Gamma_{\textrm{HI}}$
measurements, which may be artificially constricting our uncertainties. Another factor is that we do not include the full error covariance matrix for
$\Gamma_{\textrm{HI}}$
measurements from Becker & Bolton (Reference Becker and Bolton2013) (which we note is unavailable in Bosman et al. Reference Bosman2022 and Gaikwad et al. Reference Gaikwad2023). Despite these factors, the relatively good agreement between the red points in Figure D1 and our errors is encouraging.
In Figure D1, we show the full breakdown of contributions to our
$1\sigma$
errors across all redshifts from uncertainty in
$\alpha$
,
$\beta_{\textrm{N}}$
,
$\Gamma_{\textrm{HI}}$
, and
$\lambda_{912}^{\textrm{mfp}}$
(see Table 1). In each panel, we show the fiducial fit and total uncertainty (
$2\sigma$
) compared to the error arising from varying only one parameter at a time. As seen from Table 1,
$\beta_{\textrm{N}}$
and
$\Gamma_{\textrm{HI}}$
dominate the uncertainty budget at most redshifts, and
$\alpha$
is always sub-dominant. Uncertainties from
$\lambda_{912}^{\textrm{mfp}}$
are negligible at
$z \lt 5$
, but become important at
$z \gt 5$
and are comparable to the other major sources of uncertainty at
$z = 6$
.
Breakdown of the contribution to the total error budget for
$\dot{N}_{\textrm{ion}}$
from each source of uncertainty. Each panel shows the full errors (black curve and shaded region,
$2\sigma$
) and the contribution from one parameter at a time. We see that
$\Gamma_{\textrm{HI}}$
and
$\beta_{\textrm{N}}$
are comparable and dominate sources of error at most redshifts, and
$\alpha$
is always sub-dominant. Uncertainties in
$\lambda_{912}^{\textrm{mfp}}$
are only important at
$z \gt 5$
.
Open access
