Dark-ages reionization and galaxy formation simulation — XXI. Constraining the evolution of the ionizing escape fraction

The fraction of ionizing photons that escape their host galaxies to ionize hydrogen in the inter-galactic medium (IGM) is a critical parameter in analyses of the reionization era. In this paper we use the Meraxes semi-analytic galaxy formation model to infer the mean ionizing photon escape fraction and its dependence on galaxy properties through joint modelling of the observed high redshift galaxy population and existing constraints on the reionization history. Using a Bayesian framework, and under the assumption that escape fraction is primarily related to halo mass, we find that the joint constraints of the UV luminosity function, CMB optical depth, and the Ly 𝛼 forest require an escape fraction of ( 18 ± 5 ) % for galaxies within haloes of 𝑀 ≲ 10 9 M ⊙ and ( 5 ± 2 ) % for more massive haloes. In terms of galaxy properties, this transition in escape fraction occurs at stellar masses of 𝑀 ★ ∼ 10 7 M ⊙ , nearly independent of redshift. As a function of redshift, reionization is dominated by the smaller 𝑀 ★ ≲ 10 7 M ⊙ galaxies with high escape fractions at 𝑧 ≳ 6 and by the larger 𝑀 ★ ≳ 10 7 M ⊙ galaxies with lower escape fractions at 𝑧 ≲ 6. Galaxies with star formation rates of 10 − 2 . 5 M ⊙ yr − 1 to 10 − 1 . 5 M ⊙ yr − 1 provide the dominant source of ionizing photons throughout reionization. Our results are consistent with recent direct measurements of a ∼ 5% escape fraction from massive galaxies at the end of reionization and support the picture of low mass galaxies being the dominant sources of ionizing photons during reionization.


INTRODUCTION
The Epoch of Reionization, during which the progenitors of modern day galaxies formed and reionized hydrogen throughout the intergalactic medium (IGM), is the least explored cosmic epoch.However, over the past two decades there has been significant progress in exploring the stellar growth history of these galaxies (e.g.Bouwens et al. 2017;Finkelstein et al. 2015), with a recent step-change in understanding that has arisen through the James-Webb Space Telescope (e.g.Naidu et al. 2022;Atek et al. 2022).At the same time, observations of the evolution of the ionization state of the inter-galactic medium has been narrowed towards a scenario where reionization completes at, or just below,  ∼ 6 (e.g.Bosman et al. 2022;Qin et al. 2021).This constraint has been achieved through a combination of the integral measure of the reionization history obtained from CMB measurements (e.g.Planck Collaboration et al. 2016b), and measurements of absorption signals against sources emitting Ly- photons, including Ly- forest optical depth, dark gaps and quasar Ly- emission line damping wings (e.g.Becker & Bolton 2013;McGreer et al. 2015;Bañados et al. 2018;Greig et al. 2022).
Star forming galaxies are thought to provide the dominant source of ★ E-mail: smutch@unimelb.edu.au(SJM) ionizing photons to reionize the Universe (e.g.Srbinovsky & Wyithe 2007;Bouwens et al. 2015;Qin et al. 2017).However, we are unable to directly relate the observed star forming galaxy population to the evolution in the ionization state of the IGM due to the unknown fraction of ionizing photons which escape from the ISM of galaxies.
Given the evolution of galaxy density with redshift, and the expected dependence of star-formation rate with host-halo mass, we expect that there is likely a dependence of escape fraction on these quantities (e.g.Gnedin et al. 2008;Wise & Cen 2009).For example, in a scenario where the dominant photon sources are low-mass haloes, reionization starts at higher redshift (e.g.Finkelstein et al. 2019).However, if galaxies in more massive haloes are the more dominant sources, then reionization will be delayed and conclude more rapidly (e.g.Mesinger et al. 2016;Naidu et al. 2020).Moreover, since this escape of radiation is likely through channels of low column density HI and dust, the escape fraction is expected to be direction dependent for individual galaxies and to have a large scatter with galaxy properties (e.g.Cen & Kimm 2015;Yeh et al. 2023).
Given the importance of the escape fraction to our understanding of reionization, there have been a number of studies to i) measure the escape fraction directly in individual galaxies (e.g.Pahl et al. 2022;Begley et al. 2022;Bassett et al. 2022), ii) to directly simulate the ISM of star forming galaxies in a cosmological context (e.g.Yeh et al. 2023;Wise & Cen 2009), and iii) to infer the average escape fraction properties by comparing observed star-formation rates with the reionization history (e.g.Wyithe et al. 2010;Finkelstein et al. 2019;Yung et al. 2020).
In this paper we follow the third approach, using a semi-analytic model coupled with semi-numerical calculations of reionization, which offers a computationally efficient route to calculating the connection between galaxy formation and reionization (e.g.Mutch et al. 2016;Seiler et al. 2019;Visbal et al. 2020;Hutter et al. 2021).We utilize the Meraxes semi-analytic model (Mutch et al. 2016) to provide a realistic description of the galaxy population at  ∼ 4−8, which can be extended to both redshifts higher than observed and to luminosities fainter than observed based on physical descriptions of galaxy formation physics.Meraxes couples this galaxy evolution to the evolution of the ionization state of the IGM with an assumed escape fraction.Our approach is to provide a flexible relationship between galaxy properties and the escape fraction, and then to constrain the parameters describing this relationship via a Bayesian framework.Efficient implementation of the Meraxes semi-analytic model makes this parameter space exploration possible, even though the model includes on the fly calculation of reionization and is directly run on a dark-matter halo population close to the hydrogen cooling limit within a large N-body simulation.
The paper outline is as follows.In Section 2 we outline the relevant features of Meraxes.We then describe our implementation of the statistical sampling of the parameter space for Meraxes in Section 3 before constraining our galaxy formation model parameters in Section 4. In Section 5 we describe the reionization observations we employ before presenting our escape fraction model and the resulting constraints in Section 6.We finally provide our discussion and conclusions in Sections 7 and 8, respectively.Throughout this paper, we adopt a standard ΛCDM cosmology with Ω m = 0.308, Ω Λ = 0.692 and  0 = 67.8km/s Mpc −1 .

MERAXES
The Meraxes semi-analytic model (Mutch et al. 2016) is designed to explore galaxy formation across cosmic time, with a particular focus on the Epoch of Reionization and the first galaxies.Some of its key features include a delayed supernova feedback scheme which accounts for the finite lifetime of massive stars and the use of "horizontal" merger trees whereby the entire simulation volume is processed one snapshot at a time (necessary for an efficient treatment of reionization).However, the defining feature of Meraxes is its close coupling to a modified version of the semi-numerical reionization code, 21cmFAST (Mesinger et al. 2011).This allows coupled modelling of reionization with the growth of the galaxy population, both temporally and spatially.Taking into account both local and external ionizing radiation, we can explore how subsequent star formation is affected and how this feeds back on to the progression of reionization.
As with many other semi-analytic models, Meraxes includes treatments for star formation, supernova feedback, metal enrichment, cooling, AGN feedback, and galaxy mergers.In this section we provide a broad overview of the model and any physical prescriptions that are directly relevant for this work.A full description of the model is provided in Mutch et al. (2016) and the subsequent papers referenced below.

Tiamat
Semi-analytic galaxy formation models work by post-processing dark matter halo merger trees either extracted from cosmological N-body simulations or created using Extended Press-Schechter (EPS) theory.In order to fully model reionization, where the spatial distribution of galaxies directly influences their evolution, the former is required.
In this work, we utilise the Tiamat N-body simulation, run as part of the Dark-ages Reionization And Galaxy Observables from Numerical Simulations programme (DRAGONS; Poole et al. 2016).Tiamat consists of 2160 3 collisionless particles following a Planck 2015 cosmology (Planck Collaboration et al. 2016a) in a comoving volume of (100 Mpc) 3 , resulting in a particle mass of 3.89 × 10 6 M ⊙ .This allows us to identify the majority of atomic-cooling mass haloes across all redshifts during reionization within a large enough volume to provide a statistically representative realisation of the large scale process.Tiamat has been run down to a redshift of =1.8 with 101 snapshots output between =49 and 5, evenly spaced in time with a cadence of ∼ 11 Myr, and a further 64 snapshots at  < 5, logarithmically spaced in redshift.This fine snapshot cadence at high redshift allows us to better track the progression of global reionization, which proceeds rapidly during its mid phases.It also allows us to more accurately model the stochastic nature of star formation and its impacts on the reionization process at these high-redshifts (Mutch et al. 2016).For more details on the Tiamat simulation, interested readers are referred to Poole et al. (2016) and Angel et al. (2016).Details on the halo merger tree construction can be found in Poole et al. (2017).

Star formation and supernova feedback
The star formation and supernova feedback prescriptions implemented in Meraxes closely follow those of Croton et al. (2006), with updates presented in Guo et al. (2011).
At each time step in the simulation, we calculate the amount of gas required to be accreted by each host dark matter halo in order to maintain the required baryon fraction in the presence of the intergalactic ionization field.This gas is assumed to be shock-heated to the virial temperature of the halo before cooling down to the central regions to form a rotationally supported exponential disc.Assuming angular momentum conservation of the gas as it cools, the scale radius of this disc is given by  disc =  vir / √ 2, where  is the dimensionless spin parameter of the halo (Bullock et al. 2001).
Following Kennicutt (1989), we assume that star formation in the disc only occurs once the cold gas density reaches a critical threshold, Σ crit , which can be linked to the instantaneous properties of the host dark matter halo.The resulting star formation rate,  * , is calculated by assuming that some fraction,  SF , of the cold gas,  cold , above this threshold is then turned into stars per disc dynamical time,  dyn : where the second term on the right-hand side corresponds to the cold gas mass resulting in the critical surface density for a disc with scale radius  disc .
Galaxy mergers can also drive shocks and turbulence, resulting in enhanced star formation episodes.We model the fraction of cold gas converted into stars in such events following Somerville et al. (2001): where  burst is the mass of newly formed stars,  gal and  parent are the total baryonic masses of the infalling and parent galaxies respectively, and  burst and  burst are free parameters.Following Croton et al. (2006), we fix  burst = 0.7 and  burst = 0.57.Energy injection from supernovae is thought to be the dominant feedback mechanism for controlling the growth of stellar mass (and hence the production rate of ionizing photons) at the high redshifts and low halo masses relevant for reionization (see e.g.Qin et al. 2017).We assume that this feedback mechanism is dominated by massive, short-lived (≲ 40 Myr) stars, ending in energetic Type-II supernovae.The total supernova energy injected into the ISM, Δ total , after a time, , from past star formation episodes of mass  * is given by: where  is the metallicity dependent energy injection from stars of age  −  ′ which we calculate directly from the stellar population synthesis code Starburst99 (Leitherer et al. 1999), assuming a Kroupa (2001) IMF. is a free parameter which controls the efficiency with which the injected energy couples to the surrounding ISM.The amount of gas assumed to be directly heated by this energy injection is modelled as where  is a free parameter commonly referred to as the mass loading factor.
The energy injected into the ISM is used to adiabatically heat the gas in the galaxy to the virial temperature of the host dark matter halo where it joins a hot gas reservoir.The mass of gas added to the hot reservoir, Δ hot , is given by If there is more energy available than is required to reheat the full Δ total worth of gas to the virial temperature (i.e.Δ hot > Δ total ) then the remainder of the energy is used to eject the remaining mass from the system entirely.For more details on the supernova feedback prescription, see Qiu et al. (2019).

Reionization
One of the key features of Meraxes is its treatment of reionization, which is modelled using a modified version of the semi-numerical reionization code 21cmFAST (Mesinger et al. 2011).At each time step in the simulation, Cartesian grids of halo mass, stellar mass, star formation rate, and total matter density are constructed.For this work we employ a grid resolution of 256 3 , corresponding to a side length of 0.39 Mpc.These are then filtered using an excursion set formalism to identify ionized cells where the number of ionizing photons is greater than the number of neutral atoms and recombination events.As discussed in Section 1, the key free parameter controlling this process is the escape fraction of ionizing radiation from the source galaxies,  esc .When a patch of the IGM is reionized by UV photons, the temperature of the gas is raised to ∼ 2 × 10 4 K.This acts to raise the local Jeans mass, reducing the amount of material accreted by low mass haloes.We model this process using a baryon fraction modifier: where Δ baryon is the freshly accreted baryon mass,  b = Ω b /Ω m is the universal baryon fraction,  baryon is the mass of baryons already present in the halo, and  mod is the baryon fraction modifier.The value of  mod is calculated following Sobacchi & Mesinger (2013) and depends on the mass of the halo, the strength of the local UV ionizing background (which is again strongly dependent on the ionizing escape fraction), and the time over which the halo has been exposed to this background.For further details of the reionization prescription, see Mutch et al. (2016).

UV luminosities
In order to generate UV luminosities for each galaxy at a time , we sum the luminosities from each past star formation burst of age .These luminosities are calculated by constructing a full spectral energy distribution (SED) and applying a tophat filter centred on 1600Å, with a width of 100Å: where ( where ,  GCD and  are free parameters and  disc and  cold are as introduced above.The resulting transmission is calculated using the two phase model of Charlot & Fall (2000): where  ISM and  BC are again free parameters.Here, photons emitted by young stars with ages less than  BC must travel through both a dense birthcloud and the more extended ISM.After this time, the birthcloud is assumed to have been destroyed and the photons are only attenuated by the ISM.For more detail on both the SED construction and dust model see Qiu et al. (2019).

Unique features for parameter space exploration
One of the key benefits of semi-analytic models over more detailed hydrodynamic simulations is their relatively low computational expense.This opens up the potential for statistical exploration of the underlying free parameter space (e.g.Henriques et al. 2013;Qiu et al. 2019).Meraxes has been designed to maximise this potential by being written to be as efficient as possible without loss of flexibility.Some of the relevant unique features of our model include the pre-loading of input data, the utilisation of graphics processing units (GPUs), and the suite of companion tools which have been developed to programmaticaly interact with the model and its output.
In order to provide a converged sampling of the underlying posterior probability distribution we require many tens-of-thousands of model evaluations for moderate dimensionality problems (∼10 free parameters).In this situation, it is imperative that the model being evaluated can be run as quickly and efficiently as possible.For Meraxes (running on our input N-body simulation, Tiamat) at least ∼70% of the time taken for a single model call is devoted to reading in the trees and dark matter density grids from disk, distributing them amongst all available processors, and processing their contents.Therefore, we have modified the SAM to be able to read in all input data once, store this in memory, and then to run multiple times without the need for any expensive input-output (IO) calls.
A further computational bottleneck for Meraxes is its treatment of reionization using a modified version of the 21cmFast (Mesinger et al. 2011) semi-numerical code.This currently requires the model to filter a number of 3D Cartesian grids using an excursion set formalism at every snapshot.Each of these excursion set calculations entails a large number of Fourier transforms and multiple loops over every cell (of which there are 256 3 for this work) in the gridded volume.To reduce the amount of time spent doing this expensive calculation, and thus allow our analysis to be completed in a timely fashion, the reionization portion of Meraxes has been re-written to take advantage of the Graphics Processing Units (GPUs) commonly available on the majority of modern high performance computing systems.This has sped up the reionization calculation by a factor of 8 for typical combinations of grid dimensions and hardware resources.

PARAMETER SPACE SAMPLING
In order to carry out a statistical exploration of the Meraxes parameter space, we require a method to calculate the likelihood of the model predictions for a given set of parameters and to then provide the code with a new set of trial parameters for the next run.To facilitate this, as well as future parameter space exploration and calibration studies with Meraxes, we have developed a new companion code called Mhysa which provides a flexible two-way communication layer between the model (written in C) and Python.Mhysa was designed to efficiently facilitate calling Meraxes many times, with different parameter values, and allows us to interface any Python optimisation or sampling package with the model.It also allows us to write all of our likelihood calculations in Python and to leverage its many available scientific and numerical libraries, vastly reducing development time.In practice, Mhysa deals with initialising the model, updating its parameters, directly accessing in-memory galaxies and haloes, calling arbitrary Python code at the end of each snapshot, doing calculations on an independent processor in the background whilst the model is running, and orchestrating multiple instances of the model running concurrently.
A further key feature of Mhysa which greatly enhances future science investigations with Meraxes, is that it is not only limited to statistical sampling.Often we wish to simply find the 'best" model which matches our chosen observational constraints, without necessarily investigating parameter degeneracies or uncertainties.Mhysa has been written with a flexible interface that can be used to carry out such a global optimisation of the model parameters using any suitable technique, and leveraging the many publicly available packages.
For this work, we make use of the publicly available Python package Ultranest (Buchner 2021) to carry out our parameter space sampling using the nested sampling (Skilling 2004) Monte Carlo algorithm MLFriends (Buchner 2014(Buchner , 2017)).The benefit of nested sampling over more traditional Markov Chain Monte Carlo (MCMC) methods for this work is four-fold: Nested sampling does not require us to have an estimate of the global maximum likelihood parameter values in order to avoid a lengthy burn-in phase; it can more efficiently deal with multi modality or complicated degeneracies in probability distributions; it typically requires fewer likelihood evaluations to probe the underlying posterior probability distribution to a satisfactory degree of accuracy; and it naturally provides us with an estimate of the Bayes evidence that can be used to quantitatively compare the support for different escape fraction parametrisations.
For each of the parameter space sampling runs presented in this work, we use a minimum of 200 live points.This is the number of points which are maintained throughout the sampling process and iteratively updated to sample the posterior.More live points result in a higher resolution sampling of the posterior, at the expense of more likelihood evaluations.The default Ultranest stopping criteria were employed, with the exception of the maximum fraction of the remaining evidence and minimum effective sample size, which we set to 0.1 and 600, respectively.

THE GALAXY UV LUMINOSITY FUNCTION AND
−  UV RELATION The galaxy UV luminosity function is the most direct observational probe of the reionizing galaxy population at  ≥5.It traces the recent star formation at a fixed number density and hence, the total ionizing photon production rate and its spatial distribution (under the assumption that stars provide the bulk of ionizing photons during the EoR, e.g.Qin et al. 2017).Therefore, if we are to carry out a quantitative investigation of allowed escape fraction models, we must first calibrate our model against the observed UVLF.
The free parameters of the model which we allow to vary when fitting the UVLF are listed in Table 1.These are identical to the free parameters used by Qiu et al. (2019) and include four galaxy formation parameters controlling star formation and supernova feedback, as well as four parameters which control dust extinction and how it varies with cold gas mass.These parameters have been shown in previous works to be the most important for controlling the growth of early galaxies in our model (e.g.Mutch et al. 2016) and reproducing the evolution of the UVLF (Qiu et al. 2019).
The assumed prior distributions for each parameter are presented in Appendix A. They are selected to weakly encode our past experience and physical intuition, however, we have ensured that our conclusions are insensitive to our precise choices except where explicitly stated in the text.
For this work, we utilise the observational luminosity functions from Bouwens et al. (2015).Following Gillet et al. (2020) (and as originally suggested by Bouwens et al. 2017) we enforce a minimum uncertainty of 0.08 dex on all data points.Since reionization is known to be largely complete to within  HI ∼ 0.05 by  ∼ 5 (e.g.Bosman et al. 2022;Qin et al. 2021), we only consider redshifts greater than this.The volume of our input N-body simulation (100 Mpc) 3 also limits the brightest galaxies for which Meraxes can provide statistically significant predictions, due to the low number density of these biased objects.Therefore, at each redshift, we only consider luminosities below which observations indicate there should be at least 10 galaxies.
We calculate the likelihood of a model parameter set, , assuming uncorrelated, (log-)Normally distributed, errors: where  ,obs and  ,mod are the observations and model predictions for data point , respectively,  2  =  2 ,obs +  2 ,mod is the corresponding variance, and   is the total number of data points.Note that we use the appropriate value of  ,obs when the model prediction is greater than or less than the observational data point if the quoted observational uncertainty is non-symmetric.Our choice of uncorrelated log-normal uncertainties is simply our best estimate of the true uncertainty distribution in the absence of further information (such as a covariance matrix) and is a commonly made assumption.For the model variance, we assume a Poisson distribution with a mean equal to the galaxy number counts in each bin.
In order to fully constrain all eight free model parameters and break degeneracies in the dust extinction model, Qiu et al. (2019) found it necessary to additionally constrain against the evolution in the UV slope () as a function of magnitude ( UV ).We do similarly here, utilising the biweight mean observations of Bouwens et al. (2014) and again assuming uncorrelated, (log-)Normally distributed errors as per Equation 12.

Constraints on galaxy formation parameters
We now explore the constraining power of the UVLF and  −  UV relation on the galaxy formation parameters of our model.In Mutch et al. (2016), we demonstrated that reionization feedback only plays a minor role in shaping the evolution of the observable  ≳ 5 galaxy population.We therefore ignore our reionization constraints and turn off reionization feedback for this section, simplifying our interpretation and speeding up our calculations.
The model UVLFs are presented in Figure 1.Shaded bands indicate the 68% and 95% confidence intervals of the model results, Our model does a reasonable job of matching the observations at all fitted redshifts, providing confidence that we can correctly model the evolution of the observed faction of the total ionizing photon budget.See §4.1 for more details.
sampled from the underlying posterior distribution.The model is able to replicate the observational data reasonably well, across all redshifts, however, the 68% range of the model  at fixed  UV is noticeably smaller than that of the constraining observations.This stems from a tension in matching all four redshifts simultaneously, especially around −21 ≲  UV ≲ −20, where the observed / varies non-monotonically.When the model is constrained against a single redshift, the freedom afforded by the model and our choice of free parameters means that the magnitude of our uncertainties are similar to that of the constraining observations.When fitting against multiple redshifts simultaneously, the volume of acceptable model parameters reduces such that we are no longer able to reproduce the full variance of a single redshift.In Figure 2, we show the cor- sess a well-defined peak in their 1D marginalised posteriors which is independent of the imposed priors.There are a number of degeneracies between the galaxy formation and dust model parameters.For example, a higher critical density for star formation (Σ  ) leaves more cold gas available for efficient star formation in massive systems and hence requires more dust attenuation there (controlled by  GCD ) to match the bright end of the UVLF.The similarity of the  0 posterior distribution to its prior at all but the smallest values indicates that the model's supernova feedback energy coupling has a minimum viable efficiency, beyond which there is little change to the predicted UVLF.This is due to saturation effects.Once  0 is large enough, all available cold gas (as determined by the mass loading factor,  0 , and available cold gas reservoirs in the galaxy) is ejected and increasing the value of  0 has no effect.

Low redshift predictions
For this work, we choose to constrain our galaxy formation model parameters against the ≥5 observations alone.Whilst Meraxes is capable of running all the way to =0, there are numerous biases and complications associated with simultaneously constraining against high and low -redshift observations.The abundance of available observations at low redshift, coupled with the typically smaller measurement uncertainties, means that our fit would prefer matching them, at the expense of matching potentially inconsistent high-redshift data.
Having said that, any model which completely fails to reproduce relevant low-redshift observations should be deemed invalid.For this reason, we ran our model with the maximum a-posteriori galaxy formation parameters down to  = 1.8 (the lowest redshift available from the Tiamat simulation) to ensure our predictions look sensible.Figure 3 shows the predicted evolution of the cosmic star formation rate density compared against the UV-derived observational results reported by Katsianis et al. (2021).Figure 4 further shows the predicted  = 2 UVLF compared against the observations of Bouwens et al. (2021).In both cases there is a reasonable level of agreement between our best fit model and the observations, despite the fact that the model was not directly calibrated to achieve this.We are therefore satisfied that our best-fit galaxy formation parameters are valid, at least for the purposes of reproducing the UV evolution of the galaxy population.

OBSERVATIONAL CONSTRAINTS ON REIONIZATION
Whilst the evolution of the observed UV luminosity function allows us to constrain the ionizing photon budget as a function of time, it does not provide any direct constraint on the fraction of these photons which escape their source galaxy and participate in reionization (  esc ).For this, we turn to the following three observations.

Thompson scattering optical depth
The integrated optical depth due to Thompson scattering of cosmic microwave background photons by free electrons can be expressed as: where  T = 6.652 × 10 −25 cm 2 is the Thomson scattering crosssection,  m  is the mass-weighted global ionized fraction of species , and ⟨ H ⟩ = 1.88 × 10 −7 (Ω b ℎ 2 /0.022) cm −3 and ⟨ He ⟩ = 0.148 × 10 −7 (Ω b ℎ 2 /0.022) cm −3 are the average comoving density of hydrogen and helium, respectively (Wyithe & Loeb 2003).In order to calculate   from our model we follow the assumptions outlined in Mutch et al. (2016).In particular, we assume that helium and hydrogen are singly ionized at the same rate (i.e. m HeII () =  m HII () and  m HeIII ( > 4) = 0.0) until  = 4, after which helium is assumed to be doubly ionized (Kuhlen & Faucher-Giguère 2012, i.e.  m HeIII ( ≤ 4) = 1.0).We utilise the Planck Collaboration et al. (2016b) measurement of   = 0.058 ± 0.012 as our single data point.As with the ionizing emissivity constraint above, we use use a standard  2 likelihood and assume Normally distributed errors.We further assume that the uncertainty on the model prediction is negligible and that there is no missing or underestimated systematic uncertainties associated with the observational data.

Ionizing emissivity
The instantaneous ionizing emissivity (the number of meta-galactic ionizing photons per unit time, per Hydrogen atom 1 ) at the latter stages of reionization, can be inferred from Lyman-alpha forest observations (e.g.Bolton & Haehnelt 2007;McQuinn et al. 2011; 1 Often commonly expressed as the number of ionizing photons per unit time, per unit volume Becker & Bolton 2013).Following Mutch et al. (2016), we calculate this quantity from Meraxes using: where  * is the total star formation rate in the full simulation volume,   / ★ is the number of ionizing photons per stellar baryon (fixed at 6000 to be consistent with our assumed Kroupa (2001) IMF,   = Ω  /Ω  is the universal baryon fraction, and  He is the Helium mass fraction.
In this work, we constrain against the emissivities calculated by D' Aloisio et al. (2018).These were derived by combining the most recent measurements of the IGM opacity and gas temperature out to  ∼ 6 with a suite of high-resolution hydrodynamic simulations covering a broad range of IGM thermal histories.For more details on both the methodology used to obtain the observed ionizing emissivity values and the treatment of the associated errors, we refer the reader to D' Aloisio et al. (2018).We utilise a standard  2 likelihood assuming Normally distributed and uncorrelated errors (c.f.equation 12).We choose to restrict ourselves to  > 5 in this work.While observations are available at  < 5, and these ionizing emissivity data would place tighter constraints, the contribution from AGN can become non-negligible at  ≲ 4 (e.g.Trebitsch et al. 2021;Faucher-Giguère 2020).Furthermore, at  = 5 the snapshot cadence of the Tiamat N-body simulation changes from constant in time to logarithmic in redshift, introducing complications in any analysis between  ∼ 4.5 and 5.

Global neutral fraction
The volume averaged global neutral fraction can be probed observationally by a number of indirect methods.These include Lyman- source counts (e.g.McGreer et al. 2015), Lyman- and Lymanbreak galaxy clustering (e.g.Mason et al. 2018;Hoag et al. 2019), and QSO damping-wing absorption features (e.g.Greig et al. 2017;Davies et al. 2018;Greig et al. 2019;Wang et al. 2020).Due to the small available sample sizes and model-dependent nature of these probes, the associated uncertainties on individual measurements are typically large.However, a direct measure of the volume averaged global neutral fraction, xHI , can also be obtained from the fraction of pixels with zero flux in high- quasar spectra (Mesinger 2010).The appearance of these 'dark' pixels in the Lyman- and/or  forests indicates the presence of intervening neutral hydrogen, providing an upper limit to the cosmological neutral hydrogen fraction.
In this work, we utilise a number of different neutral fraction constraints assembled from the recent literature, including both model dependent and independent sources (McGreer et al. 2015;Mason et al. 2018;Hoag et al. 2019;Greig et al. 2017;Davies et al. 2018;Greig et al. 2019;Wang et al. 2020).With the exception of the model independent dark pixel constraints, we restrict ourselves to only utilising observational values with error bars, rather than upper or lower limits.This is largely motivated by the difficulty in sensibly assigning uncertainties to such limits.However, limited investigation has further suggested that the inclusion of these limits, with the most conservative assumption of a step function likelihood, leads to very little additional constraining power on the evolution of our model escape fraction in any case.
We use the 1-sigma uncertainties on each data point as presented and, following Greig & Mesinger (2017), we model upper and lower limits as a one sided Gaussian with a likelihood of 1 below (/above) the observed value.The model predictions are obtained by calculating the fraction of cells in our reionization grid with a neutral fraction less than 10 −4 at the central redshift of each data point.

RESULTS
In this section we use our best fit model (MAP) combined with a flexible halo mass dependent model for escape fraction to constrain the relation between escape fraction and galaxy properties, and to compute their contribution to the ionizing photon budget.

Constraints on escape fraction evolution
Recent numerical studies have indicated that there is a likely correlation between halo mass and the escape fraction of ionizing photons (e.g.Paardekooper et al. 2015;Kimm et al. 2017;Yeh et al. 2023) This is thought to be a result of low mass halos having shallower potential wells, allowing supernova to more efficiently clear high density gas in the IGM and provide channels for ionizing photons to escape the galaxy.In this section, we thus investigate the constraining power of the Thompson scattering optical depth, ionizing emissivity, and neutral fraction evolution on  esc as a function of halo mass.
We fix all galaxy evolution parameters to their maximum aposteriori (MAP) values found in Section 4, where we constrained the model against the evolution of the UVLF.We note that this removes any extra scatter in the recovered escape fraction introduced by uncertainties in the evolution of ionizing source population, simplifying our interpretation.However the observed star-formation rate density is constrained by the luminosity function to a ∼ 10% level at  ∼ 6 (see Figure 3), which is more precise than the observed uncertainties in the ionization rate in the IGM.The 10% uncertainty on the star-formation rate density (which is linearly related to the escape fraction with-respect-to the timing of reionization) is also smaller than the fractional constraints on the escape fraction derived below.We therefore find that the effect of simultaneously allowing our galaxy formation parameters to be free, in addition to the reionization parameters, would be minimal.
We choose a logistic functional form for the escape fraction:  /  log 10  vir midway between the two, at a mass of  offset .
We use non-flat priors for these parameters, both to encode our physical intuition, and to improve the efficiency of our sampling.More details can be found in Appendix A, however, we have ensured that the precise choice of priors does not influence results except where explicitly discussed below.
The results of the model after constraining against the three observations are presented in Figures 5 and the 1-D marginalised posterior distributions of the escape fraction parameters are presented in Figure 6.Coloured lines and shaded regions show the results of fitting the escape fraction to match each observation independently, whilst black lines and shaded regions show the results of constraining against all three observational datasets simultaneously.

Thompson scattering optical depth:
The Thompson scattering optical depth is an integrated quantity spanning the full history of the IGM back to the surface of last scattering.As such, it is sensitive to both the timing and effective mid-point of reionization.From panel (a) of Figure 5 we see that the model is able to reproduce the Planck  15).The red, green and blue curves show the constraints from the CMB optical depth, neutral fraction and ionizing photon rate observations respectively.The black curves show the combined constraints.
Collaboration et al. ( 2016a) observations (dashed line and hatched region) when using a galaxy growth history consistent with observations.
The kernel density estimate (KDE) of the 1-D marginalised constraints on the escape fraction parameters are shown by the red lines in Figure 6.From here we can see that the optical depth provides very little constraining power on the escape fraction.It merely demands that the escape fraction value be between ∼ 5−35% across the full range of contributing halo masses present in the model (10 8 ≤  ⊙ ≲ 10 13 ).
Interestingly, our model disfavours   ≲ 0.52, even when not requiring a simultaneous match to the observed ionizing emissivity and neutral fraction evolution (red shaded regions).This suggests that achieving lower optical depths requires an escape fraction model which has a stronger scaling with redshift than is possible with halo mass alone.
More recent determinations of the optical depth from Planck Collaboration et al. ( 2020) suggest   = 0.0544 +0.007 −0.008 .This is lower than the constraint used in this work (see Section 5.1) and would require a more rapid reionization with a larger minimum  max value than we see in Figure 6.The result would be an increasing tension with our constraints on  ion and  HI , both of which prefer lower  max values.

Ionizing emissivity:
The total ionizing emissivity ( ) directly constrains the number of ionizing photons entering the IGM and thus, when coupled with a fixed UVLF evolution, places the strongest constraints on the escape fraction.From the top two panels of Figure 6 we see that ionizing emissivity prefers a mass independent escape fraction (  min ≈  max ≈ 7%).For escape fractions parametrisations which vary with mass, there is a preference for a large positive slope ( / log 10  vir ) centred around  vir ∼ 10 9 − 10 10 M ⊙ .This corresponds to a rapidly increasing escape fraction with decreasing halo mass in our parametrisation, resulting in a constraint at  vir = 10 8 M ⊙ of  esc ∼ 5 − 25%.This is driven by the observed flattening in the ionizing emissivity beyond  ≲ 6, requiring both a rapid increase in the emissivity using early-forming, low mass haloes and a reduced contribution from high-mass haloes to keep the emissivity as flat as possible at low redshift (see blue shaded regions in panel (b) of Figure 5).

Neutral fraction evolution:
The range of allowed neutral fraction evolution models is shown by the green shaded regions in panel (c) of Figure 5.The dark pixel fraction upper limits of McGreer et al. (2015) ensure that reionization is all but complete by ∼5.5.Combined with the higher redshift constraints, particularly that of Hoag et al. (2019), the escape fraction parameters are driven towards producing a late and rapid reionization scenario, requiring a minimum escape fraction of ∼ 12% across all halo masses (68% confidence).The KDEs of the 1-D marginalised parameter distributions in Figure 6 again indicate that these observations are consistent with a mass independent escape fraction.However, there is a weak preference for a positive slope ( / log 10  vir ) and an increase in escape fraction with increasing halo mass.This is contrary to our physical intuition and results from the desire to have as rapid evolution in the neutral fraction as possible, whilst maintaining a late start by keeping the contribution from earlyforming, low mass haloes as small as allowed.

Combined reionization constraints:
The black shaded regions in Figure 5 and black dashed lines Figure 6 show the result of calibrating the escape fraction model using all three reionization constraints together.For these combined constraints we also present the 2-D marginalised distributions for escape fraction model parameters in Figure 7 and the corresponding summary statistics in Table 3. From Figure 5 we can see that the model is able to successfully match all of the observational constraints simultaneously.By including all constraints, we also greatly reduce the allowed range of model results, most notably those of the Thompson scattering optical depth (panel a) and ionizing emissivity (panel b).Figures 6 and 7 show that a negative  / log 10  vir (increasing escape fraction with decreasing halo mass) is strongly favoured.The midpoint of the transition between  min and  max is also restricted to a much narrower range than was the case when constraining against any one observation alone.
In Figure 8 we show the resulting constraint on the escape fraction as a function of halo mass.At  vir ≳ 10 10.5 M ⊙ the escape fraction is ∼ 5% This is driven largely by a low escape fraction in high mass halos being needed to keep the ionizing emissivity as flat as possible at  = 5 − 6.The lack of ionizing photons entering the IGM from high mass haloes reduces the ability of the model to achieve a rapid reionization, producing optical depths ≳ 0.062 and a slower neutral fraction evolution (see panels (a) and (c) of Figure 5).At lower masses, the escape fraction rises to ∼ 15 − 30%, required to start reionization early enough to match all three constraints.
The dotted line in Figure 8 shows the maximum aposteriori escape fraction parameterisation (see Table 3 for the corresponding parameter values).It is tempting to draw conclusions about physical processes colluding to put some particular significance on halo masses of  vir ≈ 10 9.35 M ⊙ , where the sharp transition in  esc occurs.However, the marginalised 68% highest density interval for the  offset parameter is broad (10 8.38 − 10 9.97 M ⊙ ; Table 3).From Figure 7, we can also see that the position of the step is correlated with the maximum and minimum escape fraction parameters.Indeed, this correlation is what drives the appearance of a gradual evolution in 68 and 95% confidence intervals of  esc vs.  vir in Figure 8.
Figure 9 shows the corresponding cumulative fractional contribution of halos to the instantaneous ionizing photon rate from halos below a given mass at different redshifts.As expected the typical halo mass producing ionizing photons increases during reionization.However, as noted above, the best fit  esc ( vir ) model favours an escape fraction that is dominated by low mass galaxies, while reionization also supresses star formation in low mass halos.Together these effects result in a very rapid transition in the halo mass above which 50% of photons are produced from  vir ≳ 10 8.5 M ⊙ at  ∼ 9 to  vir ≳ 10 10.5 M ⊙ at  ∼ 5.

Escape fraction evolution with galaxy properties
Previous studies of escape fraction evolution, both using simulations and observations, have found mean trends with galaxy properties such as stellar mass and colour (e.g.Pahl et al. 2022;Begley et al. 2022;Yeh et al. 2023).In this section, we use the MAP  esc ( vir ) parameter values to investigate the predicted evolution of the escape fraction with stellar mass and star formation rate.The upper panel of Figure 10 shows the mean escape fraction as a function of galaxy stellar mass resulting from the  esc ( vir ) model (Equation 15) with the MAP parameter values listed in Table 3.The shaded region indicates the 95% confidence intervals on the mean, generated from 1,000 bootstrap samples.Since stellar mass correlates positively with halo mass in our model (Mutch et al. 2016), the results are qualitatively similar to the dependance on halo mass.The typical escape fraction decreases from its maximum of ∼ 20% at  ★ ≲10 6 M ⊙ to ∼ 7% at  ★ ≳10 8 M ⊙ .This decrease begins at  ★ ∼ 10 6.5 M ⊙ at  ∼ 15 and shifts to smaller masses as reionization progresses due to a mild evolution of the typical stellar mass hosted by  offset = 10 9.25 M ⊙ haloes, the pivot point of our MAP  esc ( vir ) model.
The lower panel of Figure 10 shows the corresponding cumulative fractional contribution of galaxies below a given stellar mass to the instantaneous ionizing photon rate.Since stellar mass correlates positively with halo mass we see a similar trend as Figure 9 where the stellar mass above which galaxies dominate the ionizing photon budget increases during reionization.The stellar mass above which 50% of photons are produced is predicted to be close to a constant  2022).We highlight that these measurements were not included as constraints on our model.Lower panel: The corresponding cumulative contribution to the instantaneous ionizing photon rate from halos below a given stellar mass.
Figure 11 shows how the mean escape fraction instead varies by star formation rate for our MAP model.We see a similar evolutionary trend as was found with stellar mass.The escape fraction has a more extended evolution with star formation rate as redshift decreases.The lower panel of Figure 11 illustrates that galaxies with star-formation rates between  ★ ∼ 10 −1.5 M ⊙ yr −1 and  ★ ∼ 10 −2.5 M ⊙ yr −1 contribute the dominant sources of ionizing photons throughout reionization.

DISCUSSION
Before concluding we compare our population based findings with other recent theoretical results and observational studies.
Simulations typically show values of escape fraction that vary between a few and a few 10s of percent (e.g.Paardekooper et al. 2015;Rosdahl et al. 2022), which is in agreement with our results and with prior empirical comparisons of galaxy luminosity functions and Lyalpha forest constraints (e.g.Bolton & Haehnelt 2007;Wyithe et al. 2010).Hydrodynamic radiative transfer simulations further indicate that escape fractions decrease with halo mass for halos with  vir ≲ 10 9 M ⊙ (e.g.Kimm et al. 2017;Lewis et al. 2020).Most recently Yeh et al. (2023) have computed the escape fractions of reionizing galaxies within the Thesan project.They find that low-mass galaxies having  ★ ≲ 10 7 M ⊙ are the main drivers of reionization above  ≳ 7, while higher-mass galaxies with  ★ ≳ 10 8 M ⊙ dominate the escaped ionizing photon budget at lower redshifts.When calculating the dependence on mean escape fraction as a function of halo mass, Yeh et al. ( 2023) report a value of a few percent in larger halos with a sharp rise towards halo masses below  vir ∼ 10 9 M ⊙ (see also We highlight that these measurements were not included as constraints on our model.Lower panel: The corresponding cumulative contribution to the instantaneous ionizing photon rate from halos below a given stellar mass.Kostyuk et al. 2023).Inference studies by Park et al. (2019) also consider similar observational constraints to those used in this work (the UV luminosity function, CMB optical depth and dark pixels).
Although their UV ionizing escape fraction has a simpler, power-law relation with the host halos mass, they also find  esc around 5-30 per cent for low-mass haloes, with no significant preference on the power-law index.All of these results are in qualitative agreement with our findings, especially when the different simulation techniques and mode of analysis are considered.
In addition to theoretical results we can compare our results in light of observational estimates of the escape fraction.Using stacked spectra of Lyman break galaxies at  ∼ 3 Pahl et al. (2021) (following Steidel et al. 2018) find an average value for the escape fraction of 0.06 ± 0.01.Subsequently Pahl et al. (2022) find an anti-correlation between the escape fraction and stellar mass.In Figures 10 and 11 we reproduce the measurements of escape fraction from Pahl et al. (2022) for comparison with our model constraints (we stress that these measurements were not included as constraints on our model).The masses of these galaxies are at the upper end of those thought to be responsible for reionization ( ★ ≳ 10 8 M ⊙ ).However, our results overlap at the lowest stellar masses and star formation rates present in Pahl et al. (2021), as well as with many of the simulation studies referenced above.Whilst the correlations at high masses and star formation rates appear different, we note that we have constrained our escape fraction to follow a redshift-independent, halo mass dependent form.We have very few galaxies in these high stellar mass / star formation rate regimes due to our limited simulation volume and therefore they provide almost no contribution to our fits.
Finally, we note that, whilst our N-body simulation, Tiamat, provides a good balance between mass resolution and simulated volume, it does not fully resolve all haloes down to the redshift-dependent atomic cooling mass threshold.However, we miss less than 3 per cent of stellar mass at  ≈ 6 and less than 35 per cent at  ≲ 12, therefore we do not expect this to significantly alter our conclusions (see Section A1 of Mutch et al. 2016).

CONCLUSIONS
We have used the Meraxes semi-analytic model to infer the mean ionizing photon escape fraction and its dependence on galaxy properties through joint modelling of the observed high redshift UV galaxy luminosity function and the available constraints on the reionization history with the Meraxes semi-analytic model.Meraxes couples galaxy formation to a semi-numerical description of reionization, allowing the exploration of reionization and the growth of the galaxy population both temporally and spatially.
In this paper we present the application of nested sampling to Meraxes, allowing statistical exploration of the underlying free parameter space, including identification of the best fit model as well as statistical confidence intervals.Using a flexible model that relates the escape fraction to halo mass, we use the joint constraints of the UV luminosity function, CMB optical depth, and the Ly forest to evaluate the allowed range of escape fraction values and halo mass dependence.Our results show that available constraints require an escape fraction of (18±5)% for galaxies within halos of  ≲ 10 9 M ⊙ and (5 ± 2)% for larger mass halos.
Because Meraxes provides the star-formation rates and stellar masses for each dark matter halo at each simulation snapshot, we can use these results to estimate the dependence of escape fraction on stellar mass and redshift.We find that there is a transition from an escape fraction of ∼ 18% to a smaller value of ∼ 5% which occurs at stellar masses of  ★ ∼ 10 7 M ⊙ .We find that this transition value is quite insensitive to redshift.Early in reionization it corresponds to star formation rates of  ★ ∼ 10 −1 M ⊙ yr −1 , and decreases towards  ★ ∼ 10 −2 M ⊙ yr −1 by the end of reionization at  ∼ 5.
When considering the contribution of galaxies to the ionizing photon budget as a function of redshift, we find that reionization is dominated by the smaller  ★ ≲ 10 7 M ⊙ galaxies with high escape fractions at  ≳ 6, and by the larger galaxies  ★ ≳ 10 7 M ⊙ with lower escape fractions at  ≲ 6.These stellar masses correspond to galaxies with star formation rates of between  ★ ∼ 10 −1.5 M ⊙ yr −1 and  ★ ∼ 10 −2.5 M ⊙ yr −1 , which represent the dominant sources throughout reionization.In sum, these results agree with recent direct measurements of a ∼ 5% escape fraction from massive galaxies at the end of reionization, and support results from hydrodynamic modelling and direct observation indicating that the escape of ionizing photons is dominated by low mass galaxies during reionization.).These are compared with the corresponding prior distributions (blue lines; see Tables 3 & 1).A similarity between posterior distributions and priors would indicate a lack of constraining power.

Figure 1 .
Figure1.The best-fit galaxy UVLF (shaded regions) obtained by fitting against the observed evolution of the UVLF alone (error bars).Dark and light shading indicates the 68 and 95% confidence intervals, respectively.Our model does a reasonable job of matching the observations at all fitted redshifts, providing confidence that we can correctly model the evolution of the observed faction of the total ionizing photon budget.See §4.1 for more details.

Figure 2 .Figure 3 .Figure 4 .
Figure 2. The best-fit galaxy - UV relation (shaded regions) obtained by fitting against the observed evolution of the UVLF and - UV relation alone (error bars).Dark and light shaded regions indicate the 68 and 95% confidence intervals, respectively.See §4.1 for more details.

Figure 5 .
Figure 5.The reionization constraints used in this work (tan colour), along with the mass dependent escape fraction model constrained against each observation independently (red, blue and green shading) and simultaneously (black shading).Dark and light shading indicate the 68 and 95% confidence intervals, respectively.The Meraxes galaxy formation parameters are fixed to their maximum a posteriori values as presented in table 2. See §6.1 for discussion of these results.Panel (a): The Thompson scattering optical depth ( §5.1).Panel (b): The evolution of the total ionizing emissivity ( §5.2).The inset panel shows a zoom-in of =5-6.Panel (c): The evolution of the volume weighted neutral fraction ( §5.3).

Figure 6 .
Figure 6.Kernel density estimate plots of the 1-D marginalised constraints on the escape fraction parameters (Equation15).The red, green and blue curves show the constraints from the CMB optical depth, neutral fraction and ionizing photon rate observations respectively.The black curves show the combined constraints.

Figure 7 .
Figure 7.The 1 and 2-D marginalised posterior distributions for our escape fraction model parameters when fitting to all three observational constraints simultaneously.Blue cross-hairs indicate the maximum a-posteriori parameter set.The 2-D contours represent 1 and 2- probability density regions.

Figure 9 .
Figure 9.The cumulative contribution to the instantaneous ionizing photon rate from halos below a given mass at different redshifts resulting from the best fit  esc (  vir ) model constrained to match all reionization datasets simultaneously.

Figure 10 .
Figure10.Upper panel: The evolution of mean escape fraction vs. galaxy stellar mass resulting from the best fit  esc (  vir ) model constrained to match all reionization datasets simultaneously.Results are shown at 5 values of redshift.Shaded regions represent the 95% confidence intervals of the mean, calculated from 1,000 bootstrapped samples.Data points are the observed ∼3 results ofPahl et al. (2022).We highlight that these measurements were not included as constraints on our model.Lower panel: The corresponding cumulative contribution to the instantaneous ionizing photon rate from halos below a given stellar mass.

Figure 11 .
Figure11.Upper panel: The evolution of mean escape fraction vs. instantaneous star formation rate resulting from the best fit  esc (  vir ) model constrained to match all observational datasets simultaneously.Results are shown at 5 values of redshift.Shaded regions (only visible for =15) represent the 95% confidence intervals of the mean, calculated from 1,000 bootstrapped samples.Data points are the observed ∼3 results ofPahl et al. (2022).We highlight that these measurements were not included as constraints on our model.Lower panel: The corresponding cumulative contribution to the instantaneous ionizing photon rate from halos below a given stellar mass.

Table 1 .Figure A1 .
Figure A1.Kernel density estimates of the marginalised PDFs for each galaxy formation parameter (orange; see Section 4.1) compared with the corresponding prior distributions (blue lines; see Tables2 & 1).A similarity between posterior distributions and priors indicates a lack of constraining power (e.g.log 10 (  0 )).

Figure A2 .
Figure A2.Kernel density estimates of the marginalised PDFs for each  esc model parameter when fitting against the combined reionization constraints (orange; see Section 6.1).These are compared with the corresponding prior distributions (blue lines; see Tables3 & 1).A similarity between posterior distributions and priors would indicate a lack of constraining power.

Figure A3 .
Figure A3.The 1 and 2-D marginalised posterior distributions for our galaxy formation model parameters when constraining against the observed evolution of the UVLF and - UV relation.Blue cross-hairs indicate the maximum a-posteriori parameter set.The 2-D contours represent 1 and 2- probability density regions.
, ) is the metallicity () dependant star formation rate of a burst () of age   ,   (  , ) is the wavelength and metallicity dependent luminosity of a simple stellar population of age   , and   (  ) is the dust transmission function.We calculate   (  , ) us- ′ −

Table 1 .
The free galaxy formation and dust model parameters used in this work when constraining against the observed UVLF and  −  UV relations.See §4 for details.

Table 2 .
The maximum a-posteriori (MAP) for each galaxy formation parameter when constraining against the observed evolution of the UVLF and - UV relation.The 1-D marginalised MAP, mean, and 68 and 95% highest density intervals (HDI) are also presented for each parameter.Plots of the 1-D and 2-D marginalised posterior distributions are presented in FigureA3.

Table 3 .
The maximum a-posteriori (MAP) for each escape fraction model parameter (Equation15) when constraining against the observed Thompson scattering optical depth, ionizing emissivity evolution, and global neutral fraction evolution.The 1-D marginalised mode, mean, and 68% and 95% highest density intervals (HDI) are also presented for each parameter.Plots of the 1 and 2-D marginalised posterior distributions are presented in Figures6 and 7, respectively.The predicted escape fraction evolution when constrained against the Thompson scattering optical depth, ionizing emissivity, and neutral fraction evolution.The galaxy formation parameters are fixed to the maximum a-posteriori (MAP) values found in Section 4. Dark and light shaded regions indicate the 68 and 95% confidence intervals, respectively.The dotted line shows the MAP escape fraction evolution.A general trend of increasing escape fraction with decreasing halo mass can be seen, in good agreement with both theory and other works.