The Star Formation Efficiency during Reionization as Inferred from the Hubble Frontier Fields

A recent ultraviolet luminosity function (UVLF) analysis in the Hubble Frontier Fields, behind foreground lensing clusters, has helped solidify estimates of the faint-end of the z ∼ 5–9 UVLF at up to 5 mag fainter than in the field. These measurements provide valuable information regarding the role of low-luminosity galaxies in reionizing the universe and can help in calibrating expectations for JWST observations. We fit a semiempirical model to the lensed and previous UVLF data from Hubble. This fit constrains the average star formation efficiency (SFE) during reionization, with the lensed UVLF measurements probing halo mass scales as small as M ∼ 2 × 109 M ⊙. The implied trend of SFE with halo mass is broadly consistent with an extrapolation from previous inferences at M ≳ 1010 M ⊙, although the joint data prefer a shallower SFE. This preference, however, is partly subject to systematic uncertainties in the lensed measurements. Near z ∼ 6, we find that the SFE peaks at ∼20% between ∼1011 and 1012 M ⊙. Our best-fit model is consistent with the Planck 2020 determinations of the electron scattering optical depth, and most current reionization history measurements, provided the escape fraction of ionizing photons is f esc ∼ 10%–20%. The joint UVLF accounts for nearly 80% of the ionizing photon budget at z ∼ 8. Finally, we show that recent JWST UVLF estimates at z ≳ 11 require strong departures from the redshift evolution suggested by the Hubble data.


INTRODUCTION
Measurements of the galaxy luminosity function provide a fundamental input for models of galaxy formation and evolution, with Hubble Space Telescope (HST) observations probing the ultraviolet luminosity function (UVLF) out to z ∼ 10 and ongoing JWST studies reaching still earlier phases in our cosmic history.In current models of structure formation, dark matter halos collapse under gravity and galaxies form as gas falls into these halos, cools, and fragments to form stars (White & Rees 1978).This process is thought to be regulated by feedback from supernova explosions, stellar radiation, and via the energy injected into the gas from active galactic nuclei (AGN) outflows and radiation.While the abundance of dark matter halos as a function of mass and redshift is relatively well understood from N-body simulations and analytic theory (Sheth & Tormen 2002), galaxy formation and the feedback processes which regulate it, are challenging to model.
The UVLF measurements provide important empirical guidance for efforts to model galaxy formation: af-ter anchoring to the well understood abundance of dark matter halos (the "halo mass function"), the UVLF can be used to determine correlations between the UV luminosity of a galaxy and the mass of the dark matter halo which hosts it.This, in turn, is closely connected to the star formation efficiency (SFE), which quantifies the fraction of the baryons accreting onto a halo that are converted into stars.In so-called "semi-empirical models", these connections between UV luminosity, star formation rate, and host halo mass are extracted from UVLF observations and can be used to cross-check or constrain galaxy formation models (e.g.Vale & Ostriker 2004;Behroozi & Silk 2015;Sun & Furlanetto 2016).
Among other topics, the inferred star formation efficiencies have important implications for our understanding of the Epoch of Reionization (EoR).The EoR is the time period during which the first luminous sources form, photo-ionize surrounding neutral hydrogen, and gradually fill the universe with ionized gas (Loeb & Furlanetto 2013).Important goals for reionization studies include determining the volume-filling factor of ionized gas in the intergalactic medium (IGM) as a function of redshift and characterizing the properties of the ionizing sources.Joint measurements of the reionization history and the galaxy UVLF can address whether we have an accurate census of the sources of reionization, or whether missing populations of faint galaxies, accreting black holes, or other more exotic possibilities are required.
Measurements of the UVLF from HST combined with estimates of the electron scattering optical depth to cosmic microwave background (CMB) photons, τ e , from WMAP (Hinshaw et al. 2013) and Planck (Planck Collaboration et al. 2018) have started to address this question.Additional information comes from measurements of the average ionization fraction as a function of redshift (McGreer et al. 2015;Davies et al. 2018;Mason et al. 2018Mason et al. , 2019;;Umeda et al. 2023), and inferences of the emissivity of ionizing photons from the Lymanalpha (Ly-α) forest (Becker et al. 2021;Gaikwad et al. 2023).A number of previous studies have found that recent HST galaxy UVLF measurements are in accord with Planck τ e estimates provided that a significant fraction -around 10 − 20% -of ionizing photons escape their host galaxies (e.g.Robertson et al. 2015;Mashian et al. 2016;Sun & Furlanetto 2016;Yung et al. 2020).Since the escape fraction is hard to measure, especially in reionization-era galaxies, it remains unclear whether this requisite escape fraction is reasonable.There are also still sizable uncertainties related to how one extrapolates the UVLF down the faint-end and in redshift, as well as in the spectral shape of the UV emission from galaxies.
The intensity of UV radiation and its redshift evolution also have important consequences for 21 cm surveys.Specifically, UV photons redshifting into Lymanseries resonances are thought to couple the 21 cm spin temperature to the gas temperature, allowing the 21 cm signal to be observable in absorption against the CMB during the Cosmic Dawn (before e.g.X-rays heat the gas temperature above the CMB temperature leading to an emission signal; Furlanetto et al. 2006.)The long-anticipated JWST is starting to revolutionize our understanding of early galaxy populations and further address these questions in some detail.Indeed, after only months of operation, early JWST detections of z ∼ 11 − 17 photometric candidate galaxies already suggest a surprisingly large abundance of UVbright galaxies in the early universe (e.g.Donnan et al. 2022;Naidu et al. 2022;Bouwens et al. 2023b;Harikane et al. 2023b).Additionally, the stellar mass estimates for some high redshift JWST candidate galaxies appear quite a bit higher than expected (Labbé et al. 2023;Boylan-Kolchin 2023).On the other hand, the JWST analyses are still in their early stages -for one, the total number of candidates identified and the volume surveyed are still small and so statistical uncertainties remain large.
In this paper, we focus our attention on interpreting HST UVLF measurements including a recent analysis of the UVLF in the Hubble Frontier Fields (HFF) from Bouwens et al. (2022) (hereafter B22).The HFF program (Coe et al. 2015;Lotz et al. 2017) took deep multi-band images around six galaxy clusters (and six parallel fields): the foreground galaxy clusters lens distant background galaxies, which allows measurements of the reionization-era UVLF at unprecedentedly small luminosities.As we review in §2, there are a number of systematic concerns with lensed UVLF measurements: however, the community has made progress in accounting for such worries and B22 presents a comprehensive analysis with a large number of reassuring cross-checks.The B22 measurements probe absolute UV magnitudes as faint as M UV = −12.75,nearly five magnitudes fainter than probed in UVLF measurements towards typical regions (i.e., without the magnification boost from a foreground cluster).The B22 HFF analysis can also be combined with an earlier analysis from Bouwens et al. (2021) (hereafter B21) which determined the UVLF in typical regions.Naturally, the B21 measurements span more cosmic volume, contain a larger sample of galaxies, and hence provide better estimates of the bright end of the UVLF than does B22.Taken together, the combined B21+B22 data provide estimates of the UVLF across a large dynamic range in luminosity and over a decent fraction of the EoR (out to z ∼ 10).However, the combined B21+B22 results have yet to be compared with semi-empirical UVLF models.
The goal of the present work is hence to fit semiempirical models to the joint B21 and B22 UVLF measurements.Of particular interest here is the increased leverage in luminosity from the faint-end UVLF measurements of B22, which allow determinations of the SFE in small mass halos.We then explore the resulting implications for the reionization history of the universe, the electron scattering optical depth, measurements of the ionizing emissivity from the Ly-α forest, and for the timing of the 21 cm absorption signal.We will also compare our HST-calibrated models to early results from JWST.
Specifically, the outline of this paper is as follows.In §2 we briefly summarize the UVLF data from B21+B22.In §3 we describe our semi-empirical modeling approach and our procedure for fitting models to the UVLF measurements.We present the results of our fits in §4, while §5 explores their implications, including for: the SFE versus halo mass ( §5.1), the reionization history ( §5.2), our census of the sources of reionization ( §5.3), the ion-izing emissivity ( §5.4), the electron scattering optical depth ( §5.5), and the timing of the 21 cm absorption signal ( §5.6).In §6 we compare our models with early JWST results.Finally, §7 considers extensions of our fiducial model, while §8 presents our conclusions and mentions possible future research directions.Throughout, we assume a LCDM model described by the following parameters: Ω m = 0.31, Ω Λ = 0.69, Ω b = 0.049, h = 0.677, σ 8 = 0.818, n s = 0.96, in agreement with Planck 2018 results (Planck Collaboration et al. 2018).

ULTRAVIOLET LUMINOSITY FUNCTION DATA
In this work, we primarily compare models with the state-of-the-art pre-JWST UVLF measurements from the HST, using the analyses behind foreground lensing clusters from the HFF in B22 and estimates in the "field" (i.e. in typical regions sometimes referred to as "blank-fields") from B21.Here we briefly summarize the properties of these data sets, and describe potential systematic concerns with the UVLF estimates behind lensing clusters.
The HFF program obtained deep multi-band images around six galaxy clusters with well-characterized mass distributions, along with observations in six parallel fields (Coe et al. 2015;Lotz et al. 2017).This program was designed to probe the faint end of the UVLF during the EoR, exploiting gravitational lensing magnification by foreground galaxy clusters to access distant and low luminosity background galaxies.Although this technique is powerful, fully robust UVLF measurements behind lensing clusters require mitigating or accounting for a number of potential sources of systematic error.These include the need: to accurately subtract light contamination from the lensing cluster itself, to properly account for the size distribution of background galaxies (with the source size impacting the detection efficiency behind a lensing cluster), and to quantify the effect of uncertainties in the lensing models (see e.g.B22 and references therein).A great deal of progress, however, in addressing these concerns has been made through a number of HFF analyses over the years (e.g.Atek et al. 2015Atek et al. , 2018;;Bouwens et al. 2017bBouwens et al. , 2022;;Bhatawdekar et al. 2019;Castellano et al. 2016;Ishigaki et al. 2018;Livermore et al. 2017;Yue et al. 2018).
Here we make use of the most recent and comprehensive lensed analysis in the HFF from B22, and the field based measurements in B21.We refer the reader to B22 and B21 for detailed comparisons with other related works.Regarding the systematic concerns mentioned above, a few key points are as follows.First, recent studies have found that faint-end source sizes are smaller than previously expected (Bouwens et al. 2017a) and so minor completeness corrections are adequate in estimating the faint-end UVLF.Second, B22 use median magnification maps compiled from a suite of lensing models for each HFF cluster (e.g.Livermore et al. 2017;Bouwens et al. 2017b).Furthermore, B22 adopt a forward modeling approach which can account for at least much of the uncertainty in the lensing models (Bouwens et al. 2017b;Atek et al. 2018).Specifically, the authors construct mock UVLF data with one lensing model as input and analyze the maps using alternative magnification models.This procedure allows them to assess the UVLF error budget owing to uncertainties in the lensing models. 1 Finally, an encouraging result is that B22 demonstrates consistency between the faint-end slopes of the UVLF estimates behind lensing clusters and in the field.
Our following analyses make use of only measurements in redshift bins centered at or above z = 5, i.e. we consider galaxy populations during or slightly after reionization.The sample of B21 UVLF measurements used originate from 5,405 galaxies in redshift bins of width ∆z = 1 with bin centers spanning z = 5 − 10, where the data centered on z = 10 comes from Oesch et al. (2018).These measurements typically cover absolute magnitudes of M UV = −24 to −17.The B22 lensed UVLF sample used has bin centers at z = 5−9 of widths ∆z = 1 and is comprised of 525 galaxies with absolute magnitudes of about M UV = −19 to −13.
In order to account for both cosmic variance and systematic uncertainties in the lensed UVLF results, B22 recommend incorporating an additional ∼ 20% normalization uncertainty on their lensed UVLF estimates. 2In most of what follows, we neglect this uncertainty, but return to address its possible impact in §7.3.
1 Note, however, that we use the binned UVLF estimates from that study, which do not account directly for the uncertainties in the magnification maps.The magnification uncertainties are only included in the error budget for the alternate parametric fit in that work.However, we discuss the impact of uncertainties in the lensed UVLF results further in §7.3. 2 Their recommended uncertainty increases to ∼ 22% in the z = 5 bin where fewer HFF clusters are included in the analysis.(2016).As mentioned in the Introduction, these approaches are motivated by the notion that galaxy formation is driven by the gravitational collapse of dark matter halos, with gas subsequently falling into dark matter potential wells, cooling, condensing into the halo center, and fragmenting to form stars (White & Rees 1978).This process is, in turn, regulated by poorly understood feedback effects involving supernova explosions, active galactic nuclei (AGN), and photoionization, among other processes (see e.g.Silk 2011 for a brief review).Given this context, our semi-empirical modeling involves several key inputs and assumptions.Specifically, we: i) adopt a fitting formula for the halo mass function (Sheth & Tormen 2002), as derived from N-body simulations (and inspired by analytic calculations), ii) assume that the star formation rate (SFR) in a halo of mass M is proportional to the average baryonic accretion rate onto such halos, with the accretion rate determined from numerical simulations of cosmological structure formation (McBride et al. 2009;Springel et al. 2005), and iii) assume a conversion factor relating SFR and UV luminosity.As made more explicit below, the average SFE in a halo of mass M , f (M ), relates the baryonic accretion rate and the SFR, which then maps to the UV luminosity via step iii). 3Further supposing a functional form for f (M ) (or UV luminosity) and allowing for stochasticity in this relationship, we can then model the UVLF.That is, our model anchors to the well understood dark matter halo mass function and relatively well characterized accretion rates onto such halos, while it parameterizes our ignorance regarding how UV luminous galaxies populate their host halos.
As further motivated in the Introduction, this semiempirical approach calibrates the physically interesting f (M, z) relation from UVLF measurements, and this then allows predictions beyond the range in redshift, luminosity, and halo mass in which the model is calibrated.This is valuable for understanding the effects of feedback on galaxy formation, which are thought to shape f (M, z), and also for determining the properties of the sources of reionization and related topics.This section lays out the main ingredients of this methodology in more detail, including the halo mass function and accretion rate models ( §3.1), the connection between SFR, SFE, and UV luminosity ( §3.2), the conditional luminosity function approach to map between UV luminosity and halo mass ( §3.3), and alternative models 3 Sometimes the star formation efficiency is defined instead by the ratio of the stellar mass in a halo to the total baryonic mass of the halo.We denote this alternative quantity by f (M ) = M /(M Ω b /Ωm) as in e.g.Furlanetto et al. (2017).See also §5.1.
( §3.4).In addition, it is important to account for the attenuation of each galaxy's intrinsic UV luminosity by dust grains ( §3.5).The final step in our methodology is to compare models with the UVLF data: here we use the Monte-Carlo Markov Chain (MCMC) technique to efficiently sample the posterior distribution across the multi-dimensional parameter space of interest ( §3.6).

Halo Mass Function and Accretion Rate Models
The first key ingredient in our methodology is a model for the halo mass function.The halo mass function, denoted here by n(M ), describes the abundance of dark matter halos per unit mass interval such that n(M )dM is the number of halos per co-moving volume with mass between M and M + dM .We adopt the Sheth-Tormen halo mass function model (Sheth & Tormen 2002), in which: Here ρm is the mean matter density per co-moving volume and σ = σ(M, z) is the rms density fluctuation, according to linear theory, within a top-hat sphere enclosing a mass M at the cosmic mean co-moving matter density.That is, the radius of the top-hat is R = [3M/(4π ρm )] 1/3 .In the Sheth-Tormen model, the function f ST (σ) is determined uniquely by σ(M, z) relative to the threshold for spherical collapse in linear theory, δ c : (2) Here we consider the rms density fluctuations at redshift z, σ(M, z), and δ c = 1.686 is the spherical collapse threshold (at the redshift of interest). 4We drop the redshift dependence in our notation for brevity.The Sheth-Tormen mass function model is motivated by the excursion set formalism (Bond et al. 1991) and the dynamics of ellipsoidal collapse.The parameter values here, a ST = 0.707, and p ST = 0.3, provide good fits to the halo mass function measured from N-body simulations, while A ST = 0.322 is a normalization constant (Sheth & Tormen 2002).Although further simulation tests show that Sheth-Tomen is accurate at high redshifts and low masses, it may over-predict the abundance of rare massive halos at high redshift (see Lukic et al. 2007 for details).In any case, the differences with more accurate models are likely small compared to the uncertainties in the UV luminosity-halo mass relation and the dust corrections, and so we use the Sheth-Tormen model throughout.
The next input to our model is the baryonic accretion rate onto dark matter halos.Following Sun & Furlanetto (2016), we adopt the fitting formula for the matter accretion rate from McBride et al. (2009).Those authors use the Millenium simulation (Springel et al. 2005) dark matter halo catalogs to characterize the average mass accretion rate onto halos as a function of their mass and redshift.The average accretion rate was found to follow a simple scaling relation with Ṁ ∝ M δ (1 + z) η at high redshift (z ≥ 1) with δ = 1.127 and η = 2.5 (McBride et al. 2009).Although this relation was only calibrated on simulated halos with M 10 12 M and z ≤ 6, we will assume that it also applies at higher redshift and in smaller mass halos.Further supposing that the infalling matter has the universal baryon fraction of Ω b /Ω m , the baryonic accretion rate may be written as (Sun & Furlanetto 2016): with δ = 1.127, η = 2.5, and M is the total (dark matter plus baryons) halo mass.In future work, it might be interesting to test this accretion rate formula in the mass and redshift range most relevant to our modeling (see e.g. Figure 3).See Trac et al. (2015) for work in this direction; their simulation results are relatively similar to the model adopted here.

Connection to the Star Formation Efficiency
We can then connect the baryonic accretion rate with the SFR, provided that some fraction of the accreted baryons are rapidly converted into stars.Specifically, denoting this fraction, the average star formation efficiency, in a halo of mass M at redshift z, by f (M, z), we can write: In practice, we aim to calibrate f (M, z) using the UVLF measurements and so need to consider also the relationship between UV luminosity and SFR.Since most UV photons are produced by young massive stars, the UV luminosity of a galaxy is generally a good tracer of star formation.In detail, however, this connection depends on the stellar initial mass function (IMF), the stellar metallicity, stellar binarity, rotation, and the age of the stellar populations (Madau & Dickinson 2014).
Here, as in Sun & Furlanetto (2016); Madau & Dickinson (2014), we assume a constant conversion factor of: with κ UV = 1.15 × 10 −28 M yr −1 /erg s −1 Hz −1 .Here L(M, z) is the specific luminosity at a rest-frame wavelength of 1500 Å (before attenuation by dust) in units of erg s −1 Hz −1 .5As discussed in Furlanetto et al. (2017), this conversion factor is likely uncertain at the factor of a few level owing to the dependencies mentioned above.
As long as κ UV does not evolve strongly across the redshift range considered here, z ∼ 5 − 10, and/or depend significantly on halo mass, we expect the uncertainties in κ UV to mainly impact the overall normalization of f (M, z) rather than the inferred trends with redshift and halo mass.
Although for brevity of notation we have not made it explicit here, the key quantities discussed in this section: Ṁacc , SFR, f , and L, should be thought of, respectively, as the the average accretion rate, SFR, SFE, and UV luminosity in a halo of mass M at redshift z.We will account for scatter in the UV luminosity-halo mass relation when we connect the observable UVLF and the model halo mass function, as described below.

The Conditional Luminosity Function
In practice, we work backwards to infer f (M, z) from the data: given the UVLF measurements ( §2) and the halo mass function model (Equations 1-2), we determine the mapping between UV luminosity and halo mass, while the average SFE then follows from Equations 3 -5.In order to relate UV luminosity and halo mass we use a conditional luminosity function (CLF) model (Yang et al. 2003;Cooray & Milosavljević 2005;Bouwens et al. 2015;Schive et al. 2016).
The CLF quantifies the link between the halo mass function, n(M, z), and the UV luminosity function, φ(L, z).Here φ(L, z)dL is the number of galaxies per co-moving volume with specific UV luminosity between L and L + dL at redshift z.In the fiducial CLF model adopted here, we suppose that each host dark matter halo (above some limiting mass M min ) houses a single UV luminous galaxy at the center of the halo.
That is, we neglect any contribution from satellite galaxies.We further assume a "duty cycle" of unity so that all halos above M min host UV luminous galaxies at any given time.We explore alternative scenarios in §7.Explicitly, the CLF φ c (L|M, z) links the halo mass function and the luminosity function via: Throughout most of this work we consider M min = 10 8 M , comparable to the atomic cooling mass (Barkana & Loeb 2001).We comment on the impact of this assumption in what follows when relevant.Put differently, the quantity φ c (L|M, z) specifies the conditional probability that a halo of mass M at redshift z hosts a galaxy with specific UV luminosity L.
Following previous work, we suppose that the CLF follows a lognormal distribution, parameterized by the median galaxy luminosity in a halo of mass M (and redshift z), L c (M, z), and the scatter around this relation, σ (Cooray & Milosavljević 2005;Bouwens et al. 2015;Schive et al. 2016).Hence, In our fiducial model, we further follow previous work and adopt a five parameter description of the relationship between the median UV luminosity and halo mass.This description is intended to capture a range of possibilities for the feedback-regulated dependence of the SFE, and hence UV luminosity, on halo mass, and the redshift evolution in this relation.Specifically, as in Schive et al. (2016), we take: with L 0 and M 1 controlling the normalization of the relationship, p and p − q the power-law scalings with mass for M << M 1 and M >> M 1 respectively, and r the redshift dependence.
Our fiducial model adopts a mass independent value of the scatter (in the natural logarithm of the UV luminosity), σ = 0.37, equivalent to 0.16 dex of scatter (Schive et al. 2016).Note that the shape of the UV luminosity-halo mass relationship is assumed to be redshift independent in Equation 8.This is a simplifying assumption motivated in part by the limited redshift range probed in our analysis.
This functional form has been found to provide a good fit to previous UVLF measurements (Bouwens et al. 2015;Schive et al. 2016).As shown in more detail in what follows, this parameterization describes cases where the SFE peaks in the general vicinity of M ∼ M 1 and declines towards both smaller and larger masses.
Physically, this behavior may reflect the combined influence of supernova/photoionization feedback at low masses and AGN feedback at large mass scales.A key motivation of the current work is to test whether this description still applies at the small luminosities probed by the B22 HFF results, where feedback effects may be especially strong.We will also explore alternative contrasting models in the low luminosity regime ( §3.4, §7).
Finally, for completeness note that we often work in units of absolute magnitude, M UV (Oke 1974), such that M UV = 51.6 − 2.5 log 10 L/ erg s −1 Hz −1 . (9) Now φ(M UV , z) denotes the abundance of galaxies per unit UV magnitude, as a function of UV magnitude, M UV , whereas φ(L, z) gives the abundance per unit luminosity.The mapping between these luminosity functions is:

Alternative Models Considered
In addition to our fiducial CLF model of Equations 6-8, we consider several potential model variants, many of which imply enhanced SFE relative to our fiducial model at low mass scales.These cases will be discussed further in §7.However in the proceeding sections we will often consider a "flattening model" as a simple illustrative alternative to compare with our fiducial results.In the flattening model, the SFE is assumed to be a constant function of halo mass between M min and a mass scale set by an additional free parameter M 2 (see also Sun & Furlanetto 2016, their "Model I").In §7, we consider further possibilities including models where the SFE is strongly truncated below some mass scale, the scatter in UV luminosity is enhanced at small masses, or where the duty cycle of star formation activity departs from unity.

Dust Correction
The UV luminous galaxies considered here may harbor substantial amounts of dust, especially at the bright end of the luminosity function, and this will attenuate the galaxies' UV luminosities.We will follow previous work (e.g.Smit et al. 2012;Sun & Furlanetto 2016;Sabti et al. 2022) in applying an average "correction" to account for dust extinction: that is, we will estimate the average dust extinction in each UV magnitude bin of B21 and B22 and use this to determine the "intrinsic" UVLF (and revised error bars) with the impact of dust attenuation removed.We then compare models with the intrinsic UVLF throughout.
We briefly review the dust attenuation correction here, and refer the reader to the literature (e.g.Smit et al. 2012;Vogelsberger et al. 2020;Sabti et al. 2022) for further details and discussion.The dust attenuation estimate makes use of the fact that the portion of the UV emission from star-forming galaxies which is absorbed by dust is then re-radiated in the infrared.The ratio of the far-infrared to observed UV flux is hence an indicator of the amount of UV flux absorbed by dust.At least in the more local universe, this flux ratio -termed IRX for infrared excess -(IRX ≡ L IR /L) correlates with the observed UV continuum spectrum of the star-forming galaxies (Meurer et al. 1999).In other words, the UV spectral slope β (with specific flux f λ ∝ λ β ) correlates with the dust attenuation, A UV , where this latter quantity describes the number of magnitudes of dust attenuation.Hence measurements of the spectral slope β may be used to estimate the intrinsic UV emission from the observed UV luminosity (Meurer et al. 1999;Smit et al. 2012).
In order to determine the average attenuation in a given UV luminosity (or M UV ) bin, we require the mean spectral slope, β , at the luminosity in question and the scatter, σ β .Suppose that the observed dust attenuated luminosity is L obs , while the average intrinsic UV luminosity is simply L (note the brevity of notation, both are specific luminosities), then the average attenuation is: If we then assume a linear correlation with A UV = C 0 + C 1 β, the average of Equation 11 can be calculated assuming a Gaussian distribution of spectral slopes β.Specifically, the average attenuation is (Smit et al. 2012;Vogelsberger et al. 2020;Sabti et al. 2022): where β varies with UV magnitude and redshift.Equation 12 applies as long as A UV ≥ 0, otherwise we set it to zero.The dust attenuation thus depends on the parameters C 0 and C 1 which describe the correlation between A UV and β (Meurer et al. 1999), while β and σ β characterize the distribution of spectral slopes.We adopt the values C 0 = 4.54, C 1 = 2.07, σ β = 0.34 (Sabti et al. 2022;Overzier et al. 2011), and the redshift dependent form of β from Bouwens et al. (2014).We neglect any dust correction in the redshift bins with z ≥ 9, where there is too little data to make a robust estimate of dust attenuation.In any case, the dust corrections are unimportant at UV magnitudes fainter than M UV ∼ −17 at the redshifts considered here (see below for further details).
In order to determine the intrinsic UVLF, we need to apply a correction based on Equation 12 to each of: the observed UV magnitudes at the center of the bins used in the UVLF estimates, to the UV magnitude bin-widths, to the UVLF estimates φ(M UV , z), and to the measurement errors, as in Equations 3.4-3.7 of Sabti et al. (2022).
According to the estimate of Equation 12, the dust attenuation correction has a strong impact at the bright end of the luminosity functions, especially near the low redshift end of our analysis, while the correction is unimportant at the faint end of the UVLF.For example, the attenuation is A UV = 2.49 at M UV = −24 and z = 6, while the correction vanishes for UV magnitudes fainter than M UV = −18 at z = 6.At higher redshifts, the dust corrections are less important and we assume they are negligible at z ≥ 9.

Model Fitting and Parameter Estimation
One of our main goals is to obtain confidence intervals in the multi-dimensional parameter space describing the UV luminosity-halo mass relation (Equation 8).This may be characterized by a parameter vector θ = (p, q, r, M UV,0 , log 10 (M 1 /M )).The data vector d in this case describes the dust-corrected UVLF measurements in different magnitude and redshift bins.Typically, we consider the joint fit to the combination of the B21 and B22 measurements across six redshift bins from z = 5 − 10.We assume that the B21 and B22 measurements are entirely independent of each other and simply combine them.We will also consider the fit to B21 alone (i.e.including only data in the field and ignoring the lensed UVLF measurements from the HFFs).In this case, the data vector is modified accordingly.The posterior probability p(θ|d) follows from Bayes Theorem as p(θ|d) ∝ p(d|θ)p(θ), where p(d|θ) describes the likelihood of the data given the model parameters, and p(θ) characterizes the prior probabilities on the parameters.
We adopt uniform priors on each parameter over the following range θ = (p, q, r, M UV,0 , log 10 (M 1 /M )): [10,14]).This prior range was adopted based on previous CLF fits to earlier UVLF data from Bouwens et al. (2015); Schive et al. (2016).These span a conservative range in parameter values and, as we will see in Figure 2, the likelihood functions are well-peaked within the range spanned by the priors.We hence expect the results to be insensitive to the precise choice of priors here.
For the most part, we assume a Gaussian likelihood function such that (minus twice) the natural logarithm of the likelihood function (or equivalently the χ 2 value) may be written as: where the sum over i runs over both different magnitude and redshift bins.The model UVLF, φ model (i), is determined at the center of the corresponding absolute magnitude and redshift bin using Equations 1, 6-10.We assume that the measurement errors, σ(i), in different magnitude and redshift bins are uncorrelated.
We incorporate two small refinements to the treatment of Equation 13 in order to account for the sometimes asymmetric UVLF measurement errors.Specifically, some UVLF points include only upper limits on the UVLF in a given magnitude bin.In this case, we assume a half-Gaussian likelihood function.In other bins, the error bars are asymmetric with an upper error of σ + and lower error of σ − .In this case, we adopt the treatment of Barlow (2004) following their Equations 13-16.
In order to sample the posterior probability distribution and obtain confidence intervals (given the priors, p(θ), and the likelihood function, p(d|θ), Equation 13), we use the MCMC sampler from the emcee package (Foreman-Mackey et al. 2013).In our emcee runs, we start by initializing 96 walkers in a Gaussian ball with a dispersion of σ init = 0.1 around the best-fit parameters found in Appendix A of Bouwens et al. (2015).We run for over 10,000 iterations in our fiducial model and have cross-checked that the results are stable after doubling the number of iterations.Additionally, as the highest auto-correlation time is τ ∼ 70 removing the first ∼ 1000 samples should remove the influence of the initialization steps.Some of the alternative scenarios in §7 require further iterations, and we have verified convergence in those cases as well.We use the scipy.optimizemodule to find the global best fit parameters.

Conditional Luminosity Fits
Figure 1 compares the (dust-corrected) measurements from B21 and B22 with our models best fit to either the joint B21+B22 data set or to the B21 measurements alone.This comparison spans a redshift range of z ∼ 5− 10 and covers around 12 magnitudes in UV luminosity in many of the redshift bins.
Across the range of redshifts and luminosities probed by the union of the B21+B22 points, the baseline CLF model described by Equations 6-8 provides a fairly good description.Quantitatively, the best fit model to the joint data set has a χ 2 = 177.With five model parameters and 111 data points, the reduced chi-squared is χ 2 ν = 1.67.The best fit model parameters are (p = 1.69, q = 1.58, r = 1.12,M UV,0 = −23.79,M 1 = 10 11.93 M ).
Table 1 provides best fit parameter values and their corresponding marginalized 68% confidence intervals along with metrics describing goodness-of-fit for this and our other model fits.Note that our assumption of a redshift invariant shape for the median UV luminosity-halo mass relationship appears to provide a successful description of the current data across z ∼ 5−10.We will discuss the interpretations of the preferred parameter values further below in §4.3.Interestingly, the combined best fit to B21+B22 is fairly similar to a simple extrapolation from the B21 measurements alone.This is far from a foregone conclusion since the measurements in the HFFs extend down to much smaller luminosities.These low luminosity galaxies likely reside in small mass dark matter halos where the potential wells are shallow and susceptible to feedback: apparently, these effects do not significantly suppress the faint end of the UVLF over the range of luminosities probed (see also §7).In fact, as can be inferred from comparing the blue solid (joint fit) and orange dashed curves (B21 alone) in Figure 1, the combined fit has a slightly steeper faint end slope than the match to B21 alone (i.e. the UVLF appears enhanced rather than suppressed).Note that in the alternative flattening case (green dotted curve), the UVLF is even steeper at the faint end.
For quantitative comparison, the best fit to the B21 points alone yields the parameters (p = 1.84, q = 1.34, r = 0.87, M UV,0 = −23.01,M 1 = 10 11.66 M ), with χ 2 B21 = 75.2and, with five parameters and 50 data points, χ 2 ν,B21 = 1.67.If these same parameters are compared with the joint B21+B22 data set, then the goodness-of-fit values become χ 2 = 198 and χ 2 ν = 1.86.That is, the best fit parameters do shift and the old values are comparatively disfavored but not hugely so.Specifically, approximating the parameter errors from both cases as independent and adding them in quadrature, the shift in the best fit value of p is just under 2 − σ and is less in the other parameters.In particular, including the B22 data leads to a smaller p, which as we will discuss, corresponds to a slightly steeper faint-end UVLF slope.The parameter constraints do, however, tighten significantly after including the B22 lensed data with the error bars on p, M UV,0 , and log 10 (M 1 /M ) all shrinking by about a factor of 2 − 3. The increased lever arm in luminosity provided by the lensed UVLF data strongly sharpens the constraints on the faint end parameter p, which also helps tighten the confidence intervals on other model parameters by breaking degeneracies.
M UV [mag] 10 7 10 5 10 3 10 1 24 22 20 18 16 14 24 22 20 18 16 14 24 22 20 18 16 14 Dust-corrected UVLF measurements and 1 − σ error bars (or upper limits) from z = 5 − 10 compared with best-fits from three of our models.Measurements from B21 are in black and B22 in red.The solid blue curve is our CLF fit (Equations 6-10) to the joint B21+22 data set, the dashed orange curve is a CLF fit to only the B21 measurements, and the green dotted curve is a modified CLF fit to B21+22 with a flat SFE at low mass (see §3.4, §7).Note that although the lensed UVLF measurements probe up to five magnitudes fainter than in the field, the best fit to the combined B21+B22 data is similar to that from B21 alone.However, the preferred flattening model is steeper at the faint-end.

Posterior Distributions and Parameter Degeneracies
Further, Figure 2 displays the posterior distributions for our five-parameter fits to the B21+B22 data from emcee, showcasing the marginalized 68% confidence intervals on each parameter and illustrating the covariance between different parameters.Each of the model parameters are well constrained by the combined data, with fractional errors of ∼ 10% on q and r, while the other parameters are determined to better than ∼ 5%.In all cases, the posterior distributions are narrow and lie within the range of the input priors ( §3.6).The larger relative errors on q and r arise because these parameters control the bright-end luminosity halo mass relation and the redshift evolution, respectively.The galaxies at the bright end, controlling q, are rare objects and so estimates of their abundance are subject to large Poisson uncertainties.Our current comparison with the UVLF data spans a relatively small range of redshifts, z ∼ 5 − 10, and so this gives a limited handle on the redshift evolution parameter, which is likely related to the larger errors on r.
Turning to consider parameter degeneracies, Figure 2 indicates a very large negative covariance (correlation coefficient = −0.97) between the parameters log 10 (M 1 /M ) and M UV,0 .Each of these parameters is also moderately correlated with p and q: the absolute value of the correlation coefficients between log 10 (M 1 /M ) and M UV,0 with p/q range between ∼ 0.4 − 0.8 (see the correlation matrix in Figure 2).The correlations between p and q themselves are also similar in strength, whereas r is relatively uncorrelated with the four other parameters.
These degeneracies are straightforward to understand: in the simpler case of a pure power-law L c (M ) relationship, there would be a perfect degeneracy between log 10 (M 1 /M ) and M UV,0 .If L c increased as a powerlaw in M , one could compensate for an increase in L 0 (corresponding to a decline in M UV,0 ) by raising M 1 .This degeneracy, while imperfect, still applies (as is evi-p = 1.69 +0.03 0.03 1 .

.
2 .q q = 1.57+0.19  dent from Figure 2) in the slightly more complex case of Equation 8 which departs from power-law type behavior only in the intermediate mass range near M ∼ M 1 .Similar effects lead to degeneracies of p and q with each of M UV,0 and M 1 .The degeneracy between p and q occurs because the bright-end of L c (M ) behaves as a power-law in M with exponent p − q.One can compensate for an increase in p by raising q by the same amount to keep the p − q power-law fixed in this regime.1 and will be used in the figures throughout: the orange dashed curve is our model fit to the B21 measurements alone, while the blue solid curve also includes the B22 lensed UVLF measurements, and the green dotted curve is the alternative flattening scenario for the joint fit to B21+B22.Here we do not extrapolate beyond the brightest and dimmest measurements in the relevant data set.The corresponding star formation rates (from Equation 5) are shown along the y-axis to the right.The B22 measurements probe down to halo mass scales as small as M ∼ 2 × 10 9 M in our fiducial scenario (with still smaller masses probed in the flattening case), about an order of magnitude smaller than possible with the B21 observations alone.
In order to understand the implications of our bestfit parameters, it is helpful to examine the average UV luminosity-halo mass relationships in our models.This is shown in Figure 3 at redshift z = 6 for the joint B21+B22 data fit (solid blue curve), the fit to B21 alone (dashed orange curve), and the flattening fit to B21+22 (dotted green curve).The functional form of L c (M, z) in our model preserves the shape of this relation across redshift (Equation 8) and the redshift evolution is fairly mild.For each curve we show only the range of UV luminosity/halo mass scales directly probed by the data at z = 6 and so the B21 curve spans roughly M UV ∼ −24 to −17, while the fits to B21+B22 extend the reach in UV magnitudes out to M UV ∼ −13.The corresponding average SFRs, assuming Equation 5, span roughly from SFR ∼ 0.01 − 100M yr −1 (as indicated by the yaxis markings on the right-hand side of the figure).Notably, the average host halo mass in the best fit model at the faint end of the B22 measurements is around M ∼ 2×10 9 M (and still smaller in the alternative flattening scenario shown by the green dotted line).This is about an order of magnitude smaller than probed in B21 (according to our best fit models).The best fit fiducial parameters to B21+B22 put the brightest galaxies (from B21 measurements, near M UV ∼ −24) in slightly more massive halos than the fit to B21 alone, and so the blue solid curve extends to slightly larger halo masses than its orange counterpart.As previously mentioned in §3.5, we should bear in mind that the bright-end behavior is sensitive to the uncertain dust correction: the intrinsic z ∼ 6 UV luminosity is brighter by ∼ 2−2.5 magnitudes (i.e.almost an order of magnitude in luminosity) after the dust attenuation adjustment there.
The curves also reveal the double power-law form of the mean UV luminosity-halo mass relationship (in the fiducial, non-flattening scenarios), following L ∝ M p at M << M 1 and L ∝ M p−q at M >> M 16 .The best fit parameters give a transition mass around M 1 ∼ 4 − 9 × 10 11 M , with the precise value depending on which data set we fit to and the model considered (see Table 1).
Let us further comment on the behavior of the UVLF under variations in each of the model parameters.For this purpose, it is helpful to consider some approximate scaling relations derived from the simplified case of a pure power-law halo mass function n(M ) ∝ M −(γ+1) , with γ > 0. For example, in small mass halos (such that M << M 1 ), L ∝ M p , and it then follows from Equation 7 that the UVLF scales as φ(L) ∝ L −(1+γ/p) (see also Schive et al. 2016).That is, the power-law index of the halo mass function on the relevant scales (1+γ) and the low mass end of the L c (M ) relation (p) together determine the faint end slope of the UVLF, α, with φ(L) ∝ L −α .Furthermore, a smaller value of p corresponds to a steeper faint-end slope.For example, using the Sheth-Tormen halo mass function at z = 6 we find γ = 1.14 (1.27) at M = 2 × 10 9 M (M = 10 10 M ), yielding α = 1.67 (1.75) for our best fit parameter value of p = 1.69.
A similar argument can be used to understand the connection between the quantity p − q and the brightend behavior of the UVLF.For large halo masses (where δ c >> σ(M, z), see Equation 2), the halo mass function behaves as a decaying exponential while L ∝ M p−q .
Model p q r MUV,0 log 10  1. Best-fit model parameters with 68% confidence intervals, χ 2 values, the χ 2 per degree of freedom (χ 2 ν ), and ∆AIC values for our model fits to the B21+22 joint UVLF measurement set (see §4.3).The first two rows use our fiducial CLF model, but the parameters in the second row are what are preferred by a fit to the B21 measurements alone.The remaining five models are modifications to the CLF described in §7 that allow for alternative behavior in halos of mass below M2, a new sixth parameter.In the shallow-slope, bursty, and high scatter models the added seventh parameter is s, f duty , or σ, respectively, as described in §7 of the text.
Since p > q, a smaller difference between p and q corresponds to a sharper fall off at the bright-end of the UVLF.
Next, the median luminosity assumed for galaxies residing in halos of mass M is directly proportional to the parameter L 0 (Equation 8).Thus, increasing L 0 boosts the model UVLF in every magnitude bin.In contrast, for the other normalization parameter, M 1 , the reverse is true: L c (M ) is a decreasing function of M 1 for all M and so increasing M 1 lowers the model UVLF across all luminosities.Additionally, although the transition between exponential and power-law behavior in the UVLF is controlled primarily by the behavior of n(M ), M 1 serves as a scale parameter and sets the division between the faint-end and bright-end behaviors in L c (M ).The average galaxy luminosity increases only moderately with host halo mass beyond that point and, as we will see, this implies that the SFE peaks around M ∼ M 1 .
Finally, at a given L, increasing r puts higher redshift galaxies in progressively smaller mass halos.Since the halo mass function drops off less rapidly with increasing redshift on these smaller mass scales, this implies that the UVLF decreases less strongly with z when r is large.
Note that the best fit parameter values and goodnessof-fit numbers for the alternative flattening, shallowslope, bursty (i.e.non-unity duty cycle), high scatter, and truncation models (see §7 for details) are also included in Table 1.To assess whether any improvement in χ 2 is significant given the additional free parameters adopted in these models we also calculate the Akaike information criterion (AIC) differences (Banks & Joyner 2017).We defer our discussion of these results to §7.3.

IMPLICATIONS
We now turn to explore the implications of the CLF fits for our understanding of the SFE versus halo mass and its redshift evolution ( §5.1).This is, in turn, valuable for determining the properties of the sources of reionization, the expected reionization history of the universe ( §5.2), and related topics ( §5.3-5.5).

Star Formation Efficiency
In our model, bounds on the average UV luminosityhalo mass relationship may be directly translated into constraints on the average SFE as a function of halo mass and redshift.Specifically, rearranging Equations 3-5 and inserting some characteristic numbers we find: (14) Here M 1 is the break mass scale from Equation 8, while p and q are the power-law indices which characterize the halo mass dependence on either side of the break.Recall also that r describes the redshift evolution in L c (M, z).Furthermore, δ and η come from the scalings of the baryonic accretion rate with mass and redshift (Equation 3).Finally, e g(σ,MUV,0) is a normalization factor which depends on the scatter assumed in the CLF, σ 7 , and M UV,0 as: (15) In our best fit model at z = 6 the SFE peaks at an efficiency of f = 21% close to a mass scale of ∼ M 1 , as implied by the normalization of Equation 14. 8  The mean of a lognormal is related to its median by a factor of e σ 2 /2 8 Note, however, that f (M, z) does not peak exactly at mass M = M 1 , but the scaling of Equation 14 is nevertheless illustrative.The peak mass is, in general, a function of p and q but will be near M 1 for models similar to our B21+22 best fit.(solid) and z = 8 (dashed) for our three models.As in the previous figure, we show only the range of halo masses probed by each data set.In the fiducial scenario, the inferred f follows the form of Equation 14where at z = 6 it peaks at f = 21% and a halo mass of M = 6 × 10 11 M .The shaded bands represent the effects of uncertainties in model parameters on our z = 6, B21+22 fit at the 1-, 2-, and 3-σ confidence levels.These constraints, while illustrative, only show the allowed range for the UV luminosity-halo mass relationship model of Equation 8and do not capture the uncertainties in the functional form assumed here.The flattening scenario, hence, may lie outside of this band yet is also compatible with the current measurements.The red dotted line gives f (M, z = 6) (see text) for comparison.
Figure 4 shows the average SFE fits in our fiducial models, as well as in the contrasting flattening case.In our fiducial model at z = 6, the combined B21+22 data constrain the star formation across about four orders of magnitude in halo mass.The fits indicate a well defined peak in SFE near a halo mass of M ∼ 6 × 10 11 M , with f (M ) = 0.21 at the peak.The SFE declines towards small halo mass as f (M ) ∝ M p−δ ∝ M 0.56 and drops more steeply as f (M ) ∝ M p−q−δ ∝ M −1.02 at high mass.The grey bands in Figure 4 illustrate the uncertainties in the SFE for our fiducial CLF model at z = 6.Within this class of models, the trends with halo mass are well determined, although we will see that alternative scenarios are still viable.
For comparison we also compute f (M ) = M /(M Ω b /Ω m ), the stellar mass-baryonic mass ratio, as often considered in the literature and sometimes also referred to as the "star formation efficiency" although this is of course distinct from our definition here (see also Furlanetto et al. 2017).This calculation uses Equation 3 to determine the mass history of each halo at earlier times, and we combine this with our knowledge of how star formation rate/efficiency evolves with halo mass and redshift (Equations 4 & 14) to compute the corresponding stellar mass.The resulting f (M ) at z = 6 is shown as the red dotted curve in Figure 4.It is a bit smaller in normalization than f , at least below the peak mass.This is because the stellar mass reflects the cumulative history of the star formation in each halo and the SFE is lower at smaller masses and at earlier times.Nevertheless, the low mass power-law index for f (M ) and f (M ) are nearly identical.
Our results regarding the mass dependence of f and f are suggestive of the feedback-regulated star formation picture, as often invoked at lower redshift (e.g.Behroozi & Silk 2015, as well as in previous studies at high redshift, e.g.Sun & Furlanetto 2016;Furlanetto et al. 2017).In this case, the declining f (M, z) at small masses may reflect the combined impact of supernova and photoionization feedback.In fact, our best fit scaling is very similar to that found in the FIRE zoomin hydrodynamic simulations of galaxy formation (Ma et al. 2018).These authors measure the stellar to total baryonic mass as a function of dark matter halo mass, averaged over their simulation outputs at z ∼ 5 − 12 finding f ∝ M 0.58 .In the FIRE simulations, supernova feedback regulates star formation and leads to a powerlaw index which lies within the 1 − σ allowed range from our fiducial CLF results.9However, the average normalization of f in the FIRE galaxies is a factor of roughly ∼ 2 smaller than in our fits at z ∼ 6.This difference could arise because the FIRE results are averaged over a broad redshift range, but their SFE appears to evolve little over the relevant redshifts.Alternatively, our conversion between UV luminosity and SFR might be inaccurate (Equation 5), or this could point to deficiencies in the sub-grid feedback models adopted in FIRE.In any case, the match between our inferred scaling and the predictions from galaxy formation simulations is suggestive and consistent with the notion that supernova feedback regulates star formation in small mass halos.
At large mass scales, M M 1 , our z = 6 results also show evidence for a decline in SFE.The z = 8 results, however, only reach slightly beyond M 1 and do not show strong evidence for a decline.This may relate to the fact that the UVLF measurements at z ≥ 8 do not yet sample the bright end of the luminosity function.Note that at z = 6 the dust correction (Equation 12) at the bright end is about a factor of ten in luminosity, while a stronger correction would be needed to remove the evidence for a declining SFE towards high mass.At lower redshifts, where a similar decline in SFE at high mass is found (e.g.Behroozi & Silk 2015), it is common to ascribe the drop in SFE at high masses to AGN feedback.An interesting question is whether AGN can be responsible for this effect even at z ∼ 6 given the sharp decline in the AGN luminosity function at high redshifts (Kulkarni et al. 2019).Although the FIRE simulations do not show evidence for a high mass decline in f (M, z) (Ma et al. 2018), these simulations are incomplete at high masses and do not include AGN feedback.
In our best fit fiducial model the normalization of the SFE declines towards higher redshifts.As mentioned earlier, the redshift dependence follows f (M, z) ∝ (1 + z) r−η (Equation 14), reflecting evolution in both the UV luminosity-halo mass relation and the matter accretion rate.In our best fit model r = 1.12, but η = 2.5 (McBride et al. 2009) (see also Equation 3), and so even though the UV luminosity at fixed halo mass increases with redshift this is outpaced by the growing mass accretion rate.The physical origin of our derived SFE evolution is far from obvious and is worthy of future study.Our fits differ from, for example, the FIRE simulations where the SFE appears fairly constant with redshift (Ma et al. 2018).One possibility is that our model may partly mis-attribute evolution in the relationship between UV luminosity and SFR (Equation 5) to evolution in the SFE.
Also, a potential limitation of the assumed functional form in our fiducial model (Equation 14) is that it assumes a redshift invariant double power-law mass dependence.Note that it is only at the low redshift end of our sample that current data really probe the full shape of f (M, z) in the model.This is because probing the turnover mass scale requires sampling the bright end of the UVLF.Therefore, provided that the faint-end slope parameter p also evolves little across the redshift range probed, our description should be a good one.Although the redshift invariant form provides an adequate fit to current data, it will be interesting to consider more flexible models in the future.
In order to explore some possible departures from the double power-law form, we consider various alternatives such as the flattening case in more detail in §7.As a first illustration here, however, the green curves in Figure 4 show best fit flattening model results.Although the z ∼ 6 flattening model lies outside the grey band, it still provides a good fit to the data.That is, the grey band describes only the (±3 − σ) uncertainty within the context of the fiducial double power-law model and does not necessarily disfavor other possibilities.

Ionization History
Our inferences of the SFE as a function of halo mass and redshift can be used to predict the reionization his-tory of the universe, as will be discussed below.These predictions can be compared with current and future measurements.
Although numerous studies have investigated related questions (e.g.Robertson et al. 2015;Mashian et al. 2016;Sun & Furlanetto 2016;Yung et al. 2020), we consider the recent UVLF data from B22 which probe the role of individually dim galaxies in reionizing the universe.
The evolution of the average ionization fraction in the IGM, x i (z) , is controlled by the difference between the rate of production of ionizing photons per hydrogen atom, ṅion , and the rate of recombinations, here characterized by a recombination time, t rec , as: (Shapiro & Giroux 1987;Madau et al. 1999): In this context, note that ṅion receives contributions from all of the ionizing sources, including any sources of radiation that are too faint to detect individually, even in the deep lensed UVLF samples of B22.
In what follows we neglect the effect of alternative sources of ionizing photons, such as AGN, but do extrapolate our models to cover star formation in all halos above a minimum mass M min .We assume that ṅion is proportional to the star formation rate density, i.e., to the average star formation rate per unit co-moving volume.This quantity, denoted here by ρ SFR , is set by the average SFR-halo mass relationship in our models as: As usual, the SFR follows directly from the SFE, as determined by our CLF fits, along with the matter accretion rate (see Equations 3-4).
In order to predict the ionizing photon production rate per hydrogen atom (reaching the IGM), we need to further specify the typical ionizing spectrum of each galaxy, and the escape fraction of ionizing photons.Here we parameterize the spectrum through N γ , the number of ionizing photons produced per stellar baryon, which then connects ṅion and ρ SFR . 10As in Sun & Furlanetto (2016) we adopt N γ = 4, 000.The escape fraction of 10 Note, however, that for this calculation it is unnecessary to go through the intermediate step of modeling the star formation rate density.Instead, we could have directly predicted the ionizing photon production rate from the UV luminosity density.We nevertheless phrase our results in terms of ρ SFR , since this is a relatively familiar and intuitive quantity.Put differently, here ṅion only depends on the product of κ UV and Nγ rather than on each of these individually.
ionizing photons, f esc , accounts for the fact that many of the ionizing photons produced are absorbed within galaxies and only a fraction of such photons escape to ionize atoms within the IGM.Unfortunately, the escape fraction during reionization is empirically and theoretically uncertain but we adopt f esc ∼ 0.1 − 0.2 as our baseline range for comparison.We also assume that N γ and f esc are independent of halo mass and redshift.
Inserting some characteristic numbers, the ionizing photon production rate per hydrogen atom is: Here A He = 4 4−3Yp , with Y p the primordial mass fraction of helium, is a correction factor due to helium (e.g.(Sun & Furlanetto 2016)), Ω b is the baryon density parameter, and ρ SFR (z = 6) = 3.7 × 10 −2 M Mpc −3 yr −1 in our best fit model (to the B21+B22 data).Since the age of the universe at z ∼ 6 is comparable to 1 Gyr, the characteristic rate at z = 6 above corresponds to several photons per hydrogen atom over the age of the universe.
Recombinations have a sub-dominant effect on the average ionization history: the second term on the righthand side of Equation 16 is typically 20-30% as large as the ionization term during most of reionization.More specifically, the average time between recombinations is: Here T e is the temperature of the IGM (near the cosmic mean density), the scaling arises from the atomic physics of hydrogen recombination, and C HII is the clumping factor of ionized hydrogen in the IGM, defined as C HII ≡ Figure 5 shows the resulting ionization history implied by the best fit parameters in each of our three models, under the usual color code and with solid and dashed curves corresponding to f esc = 0.2 and f esc = 0.1 respectively.These models are compared to current ionization fraction measurements/bounds from: dark segments in the Ly-α and Ly-β forests (McGreer et al. 2015), Ly-α damping-wing observations towards high redshift quasars (Davies et al. 2018), observations of Lyα emission lines towards Lyman-break galaxies (Mason et al. 2018(Mason et al. , 2019)), and recent damping-wing measurements towards a stacked sample of UV luminous galaxies from JWST (Umeda et al. 2023).We find better agreement in the B21 and B21+B22 fits with f esc = 0.2, while f esc = 0.1 is slightly preferred in the flattening case.Assuming these values for the escape fraction, in broad agreement with previous work (e.g.Sun & Furlanetto 2016, Robertson et al. 2015), we find consistency with current reionization history measurements, although the measurement errors at still too large to provide a strong test.In detail, our fiducial best fit to the lensed UVLF prefers a slightly earlier completion to reionization (assuming a fixed escape fraction) than in the fit to B21 alone.This is driven by the preference for a slightly steeper faint end UVLF after including the B22 data ( §4.1).In the flattening case, a still slightly earlier completion to reionization is achieved.This results because of the more abundant faint galaxy populations in this scenario (see Figure 1 and §7.3).Hence, in addition to the large uncertainties in the escape fraction, even after including the lensed UVLF data, the importance of low luminosity galaxies remains somewhat unclear.

The Census of Ionizing Photons
We can further address the extent to which the HFF UVLF measurements, and those in the field, capture the galaxies responsible for reionizing the universe.Do In the fiducial CLF models, the peak contribution at z = 6 arises from halos with mass slightly larger than M ∼ 10 11 M and the peak mass scale shifts towards lower mass at increasing redshift.In the flattening model, lower mass halos also make an important (at z = 6) or even dominant contribution (e.g. in the z = 8 example).
Right: The fraction of the ionizing photon rate production accounted for by current UVLF measurements as a function of redshift (assuming Mmin = 10 8 M ).In our fiducial model, the lensed UVLF measurements account for ≥ 80% of the ionizing photons at z 8.The field only (B21) observations are less complete, while the contributions from low luminosity galaxies below the reach of current surveys are more important in the flattening model.The untouched region beyond z 10 is being probed by JWST.
the abundant but individually dim sources largely drive reionization or are the more luminous yet rarer sources most important?As already alluded to, this answer depends somewhat on the model assumed.
To approach this question, in the left panel of Figure 6 we consider the relative contributions to ṅion from different (logarithmic) mass bins.This characterizes the relative importance of galaxies in reionizing the universe as a function of their host halo mass (see also Yung et al. 2020).Since the results will, in general, depend on the precise SFR(M ) relation assumed.
The figure shows that, in our fiducial CLF models, the peak in d ṅion /d ln M is at a mass slightly larger than M ∼ 10 11 M at z ∼ 6, while it moves to a mass scale slightly smaller than M ∼ 10 11 M by z = 8.Although smaller mass halos are more abundant, their SFE is sufficiently small so such halos play a sub-dominant role.The B21+B22 curves extend to lower halo mass than B21 alone, showing how the additional reach towards low luminosity in the HFFs helps to enumerate the contribution from galaxies in smaller mass halos.
However, as discussed further in §7.3 the flattening case is still allowed by the present data.In this case, the z ∼ 8 d ṅion /d ln M curve (dashed green in the lefthand panel of Figure 6) actually peaks at a halo mass beneath the reach of even the lensed UVLF measurements.Consequently, although ṅion is insensitive to the precise M min in our fiducial CLF model since low mass halos have increasingly small star formation efficiencies, the results in the flattening case become dependent on M min (see below for numbers).
The right-hand panel of Figure 6 compares the fraction of the total ionizing photon production rate, at each redshift, from halos of masses accounted for by the current UVLF measurements (using Equations 8 and 20).Remarkably, in our baseline CLF model the inclusion of the lensed UVLF measurements allows about 90%/90%/80%/70% of the ionizing photon budget to be accounted for at z = 6/7/8/9, i.e. this implies only a small fraction comes from still fainter galaxies.This is improved from the corresponding numbers from B21 alone (i.e.without the lensed measurements) of 80%/60%/40%/30% at z = 6/7/8/9.However, even the lensed UVLF estimates are less complete in the flattening case and the precise fraction of the ionizing photon budget accounted for in this case depends on M min .Specifically, the fractions become 80%/70%/60%/50% at z = 6/7/8/9 for M min = 10 8 M and 70%/50%/40%/30% at z = 6/7/8/9 for M min = 10 7 M .Moreover, it is also possible that the escape fraction increases towards small halo mass, as might N γ ; in such cases, the galaxies in small mass halos would play a more important role than under the assumptions of Figure 6.It will likely require improved measurements of the reionization history, along with refined UVLF measurements, to fully assess the role of faint galaxies in reionizing the universe.Estimates of the escape fraction of ionizing photons, and its dependence on galaxy properties, will also be important.
On the other hand, we can expect improvements at z 9 from the JWST.Measurements of the lensed UVLF and in the field from the JWST may extend the current census of Figure 6 out to z ∼ 12 − 15 or so at a comparable level of completeness, with the precise limits dependent on the still-uncertain z 10 UVLF and the volume probed in JWST lensing fields (see §6).

Ionizing Emissivity & The Escape Fraction
It is also instructive to model the ionizing emissivity, i.e. the rate of production of ionizing photons per unit co-moving volume as a function of redshift.At z 6 this quantity can be inferred from observations of the Ly-α forest.More specifically, the mean transmission through the Ly-α forest can be estimated and this, in turn, implies constraints on the average hydrogen photoionization rate, Γ HI .This step generally requires comparisons with numerical simulations of the Ly-α forest.One further requires knowledge of the mean-free path to ionizing photons which can also be estimated from Ly-α forest data.
The ionizing emissivity inferred from this procedure at z ∼ 5, 6 by Becker et al. (2021) is shown in Figure 7 (see also that work for a discussion of the uncertainties in these measurements).The total rate of ionizing photon production in our model per unit volume may be computed as: where nb is the average co-moving abundance of baryons.For present purposes, the most interesting aspect of comparing this model with the measurements is that it can be used to determine the uncertain product of f esc and N γ (Kuhlen & Faucher-Giguère 2012;Sun & Furlanetto 2016).This then provides a consistency test regarding the reionization history model of the previous two sections, which required f esc ∼ 0.1 − 0.2 (assuming N γ = 4, 000). 11 Figure 7 shows the ionizing emissivity and its redshift evolution in our fiducial CLF fit (to B21+B22) for each of f esc = 0.2 (solid blue), 0.1 (dashed purple), and 0.05 (dotted cyan).Encouragingly, the f esc = 0.2 model lies within the 2−σ error range from the Becker et al. (2021) measurements at each of z = 5 and z = 6.However, the recent Gaikwad et al. (2023) points prefer a smaller value of f esc ∼ 0.05 − 0.1.
When combined with the reionization history inferences of the previous section and the electron scattering optical depth results of the next section, this might hint that the escape fraction is evolving and is smaller by the redshifts probed with the Ly-α forest than during most of reionization.Alternatively, it could indicate an inadequacy in our modeling or systematic errors in one or more of the measurements.
Overall, the errors on the emissivity measurements are still relatively large and so we conclude that our model for the ionizing sources -which matches the UVLF data from B21+B22 -is still in general agreement with the z ∼ 5 − 6 ionizing emissivity inferred from the Ly-α forest in addition to the reionization history measuremtents ( §5.2) when we adopt a constant f esc ∼ 0.1 − 0.2 and N γ = 4, 000.This is similar to the conclusion drawn by Sun & Furlanetto (2016) from earlier UVLF measurements in the field.A further test is provided by measurements of the optical depth of CMB photons to Thomson scattering, τ e .This provides an integral constraint on the ionization history and so bounds the cumulative effects from all of the ionizing sources during the EoR.It complements the UVLF measurements, which are not directly sensitive to the ionization state of the IGM, and current measurements of the reionization history/ionizing emissivity (see Figure 5-6) which are confined mostly to z ≤ 8.The statistical errors on τ e are also much smaller than the current reionization history/ionizing emissivity measurement errors.
To name one example of where the τ e constraints are especially powerful, note that a sufficiently large optical depth could indicate an additional population of high redshift z 10 ionizing sources; the B22 UVLF and current reionization history measurements would be entirely blind to such sources.Instead, as is widely appreciated, the Planck 2018 τ e measurements are consistent with a relatively rapid and late completion to reionization.Exotic early source populations are not required, and are in fact bounded, by the τ e measurements.Indeed, these measurements are entirely consistent with our baseline UVLF fits and f esc ∼ 0.1 − 0.2 model.
More quantitatively, the average electron-scattering optical depth out to a redshift z is given by: (22) Here σ T is the Thomson scattering cross section.The quantity N He (z) accounts for helium, which we as-sume to be singly-ionized along with hydrogen and then doubly-ionized below z = 3 so that: N He (z > 3) = Yp 4(1−Yp) ∼ 0.08 and N He (z < 3) = Yp 2(1−Yp) ∼ 0.16.The observable total optical depth τ e follows from taking the large z limit of Equation 22, while the cumulative contribution out to redshift z, τ e (< z), is helpful for understanding the impact of sources at varying redshifts.
Figure 8 shows τ e (< z) for our B21+22, B21, and flattening models for escape fractions of f esc = 0.1 and 0.2, along with Planck 2018's best fit τ e and the 1 − σ, 2 − σ confidence intervals (from their TT,TE,EE+lowE+lensing+BAO results, Planck Collaboration et al. 2018).The B21+22 (B21) model gives τ e = 0.062 (τ e = 0.059) for f esc = 0.2, each of which lies within the 1 − σ allowed region from Planck 2018, while the flattening model predicts τ e = 0.065 for f esc = 0.1, just outside of the 1 − σ range.Hence all three cases are consistent with the current τ e measurements.In the flattening model, low luminosity sources beyond the reach of even B22 play an important role and this leads to a more extended reionization history (see Figure 5) and a larger τ e , although this scenario is still entirely consistent with Planck 2018 observations for f esc = 0.1.The flattening case with f esc = 0.2 lies just a little bit beyond the 2 − σ confidence region from Planck 18.

UV Luminosity Density & The Detectability of the 21 cm Signal
An extrapolation of our model fits to higher redshifts also has implications for the 21 cm signal from Cosmic Dawn.In particular, the 21 cm signal is expected to be observable in absorption against the CMB at early times (Furlanetto et al. 2006).This is the case provided the gas kinetic temperature is cooler than the CMB temperature at the redshifts of interest, and as long as the excitation temperature (aka the "spin temperature") of the 21 cm line is well coupled to the gas temperature.The expectation is that UV photons from the first luminous sources will redshift into Lyman-series resonances, mix the hyperfine states, and couple the spin temperature to the gas temperature (Wouthuysen 1952;Field 1958;Pritchard & Loeb 2012).This process is referred to as the "Wouthuysen-Field" (WF) effect: the upshot here is that a sufficiently intense UV radiation field is required to keep the 21 cm spin temperature from equilibrating with the CMB temperature and to yield a measurable 21 cm absorption signal.Therefore, we can extrapolate the UV luminosity density in our model, calibrated to the Hubble UVLF measurements, predict the WF coupling strength, and hence the redshift onset of the 21 cm absorption signal.This calculation can be compared to the possible 21 cm absorption signal identified by the EDGES experiment (Bowman et al. 2018) (although see Singh et al. 2022 for a recent non-detection, inconsistent with EDGES at ∼ 2 − σ).It also has implications for the design of future global 21 cm surveys, as well as measurements of 21 cm fluctuations.Here we make use of the study of Madau (2018) which gives the range in UV luminosity densities required to achieve WF coupling, as a function of the redshift at which the coupling occurs.We can then compare this threshold with the extrapolated UV luminosities predicted in our various models (see also e.g Mirocha & Furlanetto 2019;Bera et al. 2022;Meiksin 2023;Hassan et al. 2023).Specifically, the UV luminosity density may be computed in our models as: while the threshold luminosity density to achieve WF coupling is approximately (Madau 2018): (24) This is the UV luminosity density corresponding to a WF coupling coefficient of x α = 1.The quantity g relates to a sum over different Lyman-series resonances and accounts for redshift evolution between the emission of UV photons and scattering in one of the resonances (see Madau 2018 for details).This factor is normalized to a fiducial value of g = 0.06.As in Madau (2018), we suppose that the onset of the 21 cm absorption signal corresponds to a coupling constant of x α = 1 − 3 (i.e.spanning between the value in Equation 24 and three times that number).We can then test at which redshifts our models (from Equation 23) cross these thresholds.In comparison to the earlier calculations of Madau (2018), our analysis incorporates the more recent lensed UVLF results of B22 and also adopts a physicallymotivated model to extrapolate from the UVLF measurements to higher redshifts (while Madau 2018 assumes that the UV luminosity density itself evolves as a power-law in redshift).
Figure 9 shows the UV luminosity density / star formation rate density in our models extrapolated to higher redshifts, as compared to the threshold densities required to achieve WF coupling.We show results for each of the CLF fits to B21+B22 (blue), B21 alone (orange), and in the flattening model fit to B21+B22 (green).For contrast, we also show (points labeled "observed" in legend) ρ UV implied by the models without extrapolating to fainter luminosities than currently observed.]) Figure 9. Implications for the onset of WF coupling between the 21 cm spin temperature and the gas temperature.
The curves give the UV luminosity density in our models, following the usual color code.The lines include extrapolations down to arbitrarily small luminosities (assuming Mmin = 10 8 M ), while the points only consider luminosities measured in the respective data sets and redshift bins.
The purple shaded band indicates the minimum range in UV luminosity densities required for UV radiation to couple the 21 cm spin temperature to the gas temperature Madau (2018).The red labeled region indicates the approximate redshift range for the onset of WF coupling implied by the EDGES measurement (Bowman et al. 2018).In the flattening model one can achieve coupling at the EDGES implied redshifts, but the onset of the 21 cm absorption signal is expected later (near z ∼ 13) in our other models.
In our fiducial CLF models, including faint sources and extrapolating trends in redshift, we find that coupling is achieved at z ∼ 12 − 13.In these scenarios, we expect the absorption signal to be observable only at lower redshifts than implied by EDGES.Here, the sweet spot for global 21 cm measurements appears to be around a frequency of ν = 1420 MHz/(1 + z) ∼ 100 − 110 MHz, which unfortunately lies in the FM radio band.However, one should keep in mind that the 21 cm absorption signal also depends on heating from e.g.X-rays, and so the WF coupling coefficient provides an incomplete descriptor of the signal.Interestingly, in the flattening model one expects an onset redshift close to that of the EDGES signal, and so models with rather efficient star formation in small mass halos could achieve the timing of the EDGES absorption feature (see also Mirocha & Furlanetto 2019).Even in this case, however, the depth and the shape of the EDGES feature are surprising (e.g.Mirocha & Furlanetto 2019).
The discrete points, without extrapolation towards low halo mass, are relevant for the interpretation of current 21 cm fluctuation measurements.For example, the HERA collaboration recently used upper limits on the 21 cm fluctuation power spectrum at z = 7.9 and z = 10.4 to place a lower bound on the amount of X-ray heating at these redshifts (Abdurashidova et al. 2023).That is, a cold IGM might lead to larger 21 cm fluctua-tions than observed.This argument requires that there are significant amounts of neutral gas at these redshifts (as strongly suggested by e.g. the Planck τ e measurements, §5.5), and that the 21 cm spin temperature is coupled to the gas temperature.The discrete points in the figure show that currently observable galaxies should easily yield WF coupling at z ∼ 7.9 and that coupling at z ∼ 10.4 requires only a modest extrapolation in redshift and luminosity beyond those of current UVLF measurements.This further validates the bounds in Abdurashidova et al. (2023), although we caution that the average WF coupling considered here provides only a loose figure of merit.

COMPARISON WITH EARLY JWST RESULTS
The JWST has already started to make revolutionary new measurements of high redshift galaxy populations (Donnan et al. 2022;Finkelstein et al. 2023;Bouwens et al. 2023a;Harikane et al. 2023b).Interestingly, these early results suggest large populations of UV luminous galaxies even at z ∼ 11 − 17.In addition, some galaxies have estimated stellar masses as large as M 10 10 − 10 11 M at redshifts as high as z ∼ 9 (Labbé et al. 2023).Most of these detections are photometric candidate galaxies, selected via the Lyman-break technique, and so their redshift estimates await spectroscopic confirmation.Moreover, it will be important to further assess the purity of the z 10 JWST photometric samples: since the abundance of potential lower redshift interloping galaxies likely outnumbers that of the z 10 galaxies by a large factor, even interlopers with rare/unusual spectra could be an important source of contamination (Furlanetto & Mirocha 2022).For example, a z ∼ 16 candidate galaxy (Donnan et al. 2022) turned out to be a dusty z ∼ 5 star-forming galaxy with strong H-α and [OIII] emission lines (Arrabal Haro et al. 2023), as suggested earlier based on a secondary redshift solution and nearby z ∼ 5 neighboring galaxies (Naidu et al. 2022) (see also Zavala et al. 2023).In addition, updated instrumental calibrations have led to downward revisions in the photometric redshift estimates for at least some candidate galaxies from early JWST images (Adams et al. 2023).Also, the current UVLF estimates are sensitive to the precise likelihood threshold adopted for accepting candidate galaxies (Bouwens et al. 2023a).
Finally, the total number of JWST candidate galaxies and the volume probed are still relatively small.Hence, while the early JWST results may suggest more z 10 star formation than previously thought, future efforts are required to achieve fully robust UVLF measurements in this era.
Nevertheless, we can place the early JWST UVLF estimates in context by comparing them with our models calibrated to the latest HST data from B21+B22.Here, we find relatively good agreement between the HST UVLF estimates and the new JWST measurements at z 10.However, the JWST UVLF estimates at higher redshifts land significantly above a simple extrapolation, to z 10, of our HST-calibrated best fitting model.One possibility is that the SFE stays constant or starts to increase beyond the redshifts wellprobed by HST, even though our best fit model has the SFE decreasing across the redshift range we fit to (from z = 5 − 10).To explore this, we suppose that the redshift evolution parameter transitions abruptly from the value r (see Equation 8) to an alternate value, r at z > 8 only.That is, we suppose the median UV luminosity scales as, for a given increased parameter r , L c,r (M, z) = L c (M, z) × 1+z 9 r −r at z > 8.Note that our fiducial best fit model has r = 1.12 (Table 1).
Specifically, Figure 10 compares the HST points and the new JWST UVLF measurements from Donnan et al. (2022); Harikane et al. (2023b); Bouwens et al. (2023a); Finkelstein et al. (2023) for models with r = r (our fiducial model), r = 2.5, and r = 5 at z = 7 − 17.As alluded to above, the new JWST measurements are consistent with the fiducial model and the HST points at z 10, but this model significantly underpredicts the data at higher redshifts.In the alternative r = 2.5 model, the SFE becomes redshift independent at z > 8 (recall that f (M, z) ∝ (1 + z) r −2.5 for our model, Equation 14), yet this model still falls quite a bit below the JWST data at z 13.In order to better reconcile with the JWST case, a rather extreme scenario with r = 5 is required.In this case, the implied star formation efficiencies become quite high: for example, in galaxies of absolute magnitude M UV = −20, the r = 5 scenario implies f = 13% at z = 16, compared to f = 3% in our fiducial scenario at the same redshift and mass scale.Additionally, if the shape of f (M, z) is redshift invariant as in our fiducial model, the SFE at M UV = −24 for r = 5 is f = 44% at z = 16.Further, the r = 5 case implies an abrupt reversal in the redshift evolution of the SFE to match both this behavior and the trend suggested by the HST data.
Although the early JWST results at z 10 require a rapidly increasing SFE towards high redshift, note that such scenarios do not generally overproduce the electron scattering optical depth constraints from Planck (Planck Collaboration et al. 2018).Specifically, the r = 5 scenario is consistent at the 1.4 − σ level, assuming an escape fraction of f esc = 0.2.The τ e boost is only small as star-forming halos are still relatively rare at these red-Another potential explanation considered is the presence of a "top-heavy" IMF at these high redshifts, where the formation of low-mass stars is suppressed by a higher gas temperature (Inayoshi et al. 2022;Mason et al. 2023;Harikane et al. 2023b;Finkelstein et al. 2023).The higher CMB temperature and the less efficient cooling of the lower metallicity gas at early times might both contribute to this effect (Chon et al. 2022).In this case, the conversion factor between UV luminosity and SFR, κ UV (Equation 5), would be smaller in the early universe and galaxies would be more luminous at similar star formation efficiencies.Other possibilities include an incomplete understanding of dust attenuation (Mason et al. 2023;Ferrara et al. 2022;Finkelstein et al. 2023), larger variance in the star formation rates in small halos (Mirocha & Furlanetto 2023;Mason et al. 2023;Yajima et al. 2022), a greater prevalence of AGN than expected (Harikane et al. 2023a), or more exotic options including enhanced structure formation from primordial black holes or axion dark matter mini-clusters (Liu & Bromm 2022;Hütsi et al. 2023).A closely related problem is that the stellar mass estimates for some of the JWST galaxies exceed expectations informed by structure formation in LCDM with Planck Collaboration et al. ( 2018) parameters (Labbé et al. 2023;Boylan-Kolchin 2023;Steinhardt et al. 2022).Among other possibilities, this could be related to systematic errors in the stellar mass estimates, or it might provide another indication of efficient early star formation and/or a topheavy IMF.
Further work is required to solidify the early JWST estimates, especially: obtaining additional spectroscopic confirmations of galaxy candidates at z 10 (Harikane et al. 2023a;Curtis-Lake et al. 2023;Bunker et al. 2023;Fujimoto et al. 2023), tests of sample purity, and improved statistics.Our HST-calibrated models will continue to provide useful baseline model expectations, including for upcoming lensed measurements with JWST which should access the faint-end of the UVLF at z 10.

ALTERNATE MODELS FOR THE FAINT-END BEHAVIOR AND SYSTEMATIC ERROR TESTS
As discussed in §3 our parameterization of L c (M, z), as a double power-law in halo mass and single powerlaw in redshift is motivated by previous treatments (Bouwens et al. 2015;Schive et al. 2016).These are, in turn, informed by models for the build up of dark matter halos and stellar feedback effects, which are broadly borne out by simulations (Ma et al. 2018).Yet, our relatively limited observational handle on low luminosity galaxies at high redshifts leaves open the possibil-ity of alternative behavior in this regime.Furthermore, bursty star formation may actually lead to a higher SFE in small mass halos than otherwise expected (Furlanetto & Mirocha 2022) 12 , and this will likely be accompanied by an enhanced scatter.In the following subsections we consider cases where the star formation is enhanced or suppressed, or the scatter around the UV luminosityhalo mass relationship is varied, in small mass dark matter halos relative to our fiducial model.We also quantify the impact of systematic uncertainties and cosmic variance on the lensed UVLF measurements.

Enhancing Star Formation at Low Mass
First, we consider four model extensions in which the SFE is enhanced below some critical mass, M 2 .In the first case, the SFE flattens.Second, we allow a free spectral index between M min and M 2 : we find that this prefers a shallower slope while a steeper dependence is disfavored.In the next two variants, the scatter in the UV luminosity-halo mass relation is boosted; we will see that this also increases the SFEs implied by our fits.In the first enhanced scatter model, we simply increase the variance of the CLF (Equation 7), while in the second case we allow a duty cycle to account for periods of bursty star formation.

Flattening and Shallow-Slope Scenarios
In the flattening and shallow-slope models, we suppose that the median UV luminosity-halo mass relationship scales as L c (M, z) → L c (M, z) × (M/M 2 ) δ+s−p between M min ≤ M ≤ M 2 , where the new parameters s and M 2 help quantify the departure from pure power-law behavior at small masses.Here s is defined such that the SFE scales as f (M ) ∝ M s at low masses (see Equation 14); in our fiducial model, the best fit result is f (M ) ∝ M 0.56 .Alongside the five fiducial CLF parameters we can fit for M 2 and fix s = 0 (allow s to vary) in our flattening (shallow-slope) models.Note that our prior on s in the latter case is [-0.5, 1.5], which also allows for the possibility of a steepening f (M ) relation (s > p−δ).A steeper relation might occur if feedback effects give a stronger suppression than in our baseline model (see also §7.2): however, we find that the data disfavor this possibility.
The results of these fits are found in Table 1, with f (M ) ∝ M 0.27 at small mass preferred in the shallowslope case (the 2 − σ confidence interval is 0.14 < s < 0.40 and disfavors a steepening).As detailed in the table, the other parameters vary somewhat to compensate for the changes in the model UVLFs.The largest changes are in p, but note that this parameter now reflects the slope of the UVLF at intermediate luminosities rather than the faint-end behavior.Notably, in both cases, the alternate behavior occurs at < M 2 ∼ 3 × 10 10 M : these are mass scales almost exclusively probed by the lensed UVLF measurements.Interestingly, these scenarios are, in fact, preferred (even accounting for a penalty for the greater model complexity) over our fiducial models with ∆AICs of 12.1 and 0 for the flattening and shallow-slope cases compared to the fiducial value of ∆AIC = 24.3.That is, these cases are relatively favored at 3.5 − σ and 4.9 − σ respectively.This evidence, is however, subject to systematic uncertainties in the lensed UVLF measurements ( §7.3).shows the best fit model with a shallower slope below M2, the purple dot-dashed curve allows a non-unity duty cycle at small masses, and the brown dash-dot-dotted case adopts a higher scatter at low masses.The red dashed shallowslope model is most preferred by the data.Bottom: The effect on the z = 6 UVLF from sharp truncations in star formation in halos of mass M < 5 × 10 9 M (red dashed) and M < 2.6 × 10 10 M (purple dotted, the mass scale below which we found a flattening in SFE is preferred) compared with our fiducial best fit model.The lensed UVLF strongly disfavors these and similar truncation models.
The top panel of Figure 11 shows the SFEs of our best fit enhancement models.Additionally, note that Figure 1 shows the UVLF in the best fit flattening model as an illustrative case.By construction, the flattening model has a larger SFE in small halo masses and this leads to a steeper faint-end slope in the UVLF.Specifically, the z = 6 SFE flattens at ∼ 4%, which is a factor of several larger than the fiducial best fit SFE in the lowest masses probed by B22.In the shallow-slope model at z = 6 the SFE is 2.5% for 2 × 10 9 M and 4% at 10 10 M .For further reference, the faint-end UVLF slopes become α = 2.01 (α = 1.82) for the flattening (shallow-slope) case, as compared to α = 1.67 in the fiducial model at z = 6, M = 2 × 10 9 M .Although the flattening and shallow-slope models boost the rate of ionizing photon production in small mass halos relative to our fiducial model, the changes in the reionization history are relatively small and can be compensated by lowering the uncertain escape fraction in these alternative scenarios (Figures 5,8).However, the UVLF measurements provide a less complete census of the sources of reionization (Figure 6) and as emphasized previously, the onset of WF coupling also occurs earlier in these models (Figure 9) (z ∼ 15 in the shallow-slope case).Although not shown in the figures, the corresponding behavior in the shallow-slope case is always intermediate between those in the fiducial and flattening models.

Increased Stochasticity
In our first model with enhanced scatter in small mass halos, we allow the variance parameter, σ, in φ c (L|M, z) (Equation 7) to take on a new value below M 2 .In this case, our fitting prefers a value of σ = 1.51 (from a prior range of [0.05, 2.5]) below M 2 = 5 × 10 9 M compared to σ = 0.368 in the baseline model, with ∆AIC= 9.4 compared to ∆AIC= 24.3 in the fiducial model.So the data moderately prefer enhanced scatter, although only relatively near the smallest halo mass scales probed.Table 1 shows only minor changes in most other fitting parameters.In this model, the average SFE is also increased.Specifically, Equations 14-15 show that f (M ) ∝ e σ 2 /2 and so the SFE is larger in the best fit enhanced scatter model by nearly a factor of ∼ 3 below M 2 , as compared to the fiducial case, see the top panel of Figure 11.
Next, we adopt a simplified model for a constant, nonunity duty cycle in small mass halos motivated in part by Mirocha (2020);Furlanetto & Mirocha (2022).If 0 ≤ f duty ≤ 1 represents the fraction of time galaxies in those halos are actively forming stars, in the CLF this corresponds to φ(L|M, z) → f duty φ(L|M, z) for M min ≤ M ≤ M 2 .In this model the best fit is f duty = 0.48 below M 2 = 10 11 M .As shown in Table 1, some of the other fitting parameters have moderate shifts in this case -see the top panel of Figure 11 for the impact of this on f (M ).As in the previous enhanced scatter model, the average SFE increases somewhat in small halos.At z = 6 the SFE is 2.3% for 2×10 9 M and 3.6% at 10 10 M .Similar to the other alternatives explored, this bursty star formation case is preferred relative to our fiducial model (∆AIC= 4.5 compared to our fiducial ∆AIC= 24.3).
In reality, a more complicated treatment for modeling the faint end SFE may be required.A combination of the previous extensions, a smoother transition between the efficiency in small and intermediate mass halos, and additional evolution (with mass/redshift) in σ or f duty may be warranted.Ultimately, we find many promising initial candidates for extensions to the CLF model, especially the case of a shallower power-law in SFE for the lowest mass halos, which gives χ 2 ν = 1.32, and the smallest AIC value (comparatively our fiducial model having ∆AIC= 24.3) (see Table 1).However, as discussed further in §7.3 these hints are subject to systematic concerns in the lensed UVLF measurements and so we defer further exploration to future work.

Truncation of Star Formation at Low Mass
The UVLF measurements can also be used to strongly rule out scenarios in which the SFE drops sharply in halos below some critical mass scale, M 2 .As we have discussed, feedback effects from supernovae and photoionization heating may strongly reduce the SFE in low mass halos.In our fiducial model, feedback leads to a gradual reduction in the SFE at small masses, while we adopt M min = 10 8 M as the scale below which it is strongly suppressed, as this is close to the atomic cooling mass.However, the SFE may actually be strongly reduced on mass scales somewhat larger than the atomic cooling mass: for example, photo-heating from reionization may prevent gas from accreting onto halo masses somewhat larger than our M min value (see e.g.Loeb & Furlanetto 2013 and references therein).B22 compares their UVLF measurements with various feedback models in the literature to constrain such cases.Here, rather than comparing with particular feedback scenarios, we simply bound the possibility that the SFE is truncated below M 2 and vary this scale as a free parameter.A related point is that in alternatives to CDM such as warm dark matter and fuzzy dark matter, the halo mass function itself may be suppressed at small masses (e.g.Schive et al. 2016.)The bounds on these scenarios from the new lensed UVLF measurements will be considered in future work.
In the simple truncation scenario considered here we set L c (M, z) → 0 below M 2 .Upon varying M 2 in our truncation model, we find that it is preferred to lie below the mass range probed by the measurements. 13For example, a truncation at M 2 2 × 10 9 M is excluded at 2 − σ and M 2 5 × 10 9 M is excluded at 5 − σ.The bottom panel of Figure 11 shows example truncation models with strong suppression in halos either below M 2 = 5 × 10 9 M or M 2 = 2.6 × 10 10 M , the preferred M 2 in the flattening scenario.The other parameters have been allowed to vary in the these cases.The dramatic effect on the UVLF illustrates why a truncation is excluded in the mass range probed by the B21+22 data.
Additionally, the shallow-slope model of §7.1 allowed for the possibility of a steepening of f (M ) below M 2 which could have indicated a more mild suppression in star formation.Instead, a shallower relation, indicating an SFE enhancement, is preferred.
These findings are similar to those in B22, which compared various UVLF models in the literature that suppress the faint end of the luminosity function.Here, we determine bounds on the mass scale of any suppression which should be agnostic to the details of particular models.Based on the B22 UVLF measurements, viable feedback models must not strongly suppress the SFE above M 2 × 10 9 M (at 2 − σ confidence).

Impact of Systematic Errors on Lensed UVLF Measurements
As emphasized in §2, B22 recommend an additional 20% (22% at z = 5) relative normalization error between the lensed and field UVLF data.This is intended to account for both systematic errors in the lensed measurements as well as cosmic variance in the relatively small fields probed behind foreground lensing clusters.This additional contribution to the error budget has been ignored in our analysis thus far, in part because it is non-trivial to include.In particular, while cosmic variance contributions should be essentially uncorrelated across different redshift bins, systematic errors may be correlated.Hence, a fully rigorous treatment requires knowledge of the covariance between the uncertainties in different redshift bins.We can nevertheless estimate the impact of this additional normalization error by re-running our analysis on two modified data sets, denoted here by B21+22-U and B21+22-D, where the lensed B22 measurements are all shifted up or down 13 log 10 (M 2 /M ) = 8.68± 0.46 0.46 are the best fit and 1 − σ errors (Table 1).However, the likelihood is nearly flat for M 2 10 9 M and so these values alone are not as useful for determining the constraints on a truncation in the SFE.by 20% (22% at z = 5) respectively, keeping the error bars fixed. 14These tests may overestimate the impact of this added normalization uncertainty if the offsets in fact vary with redshift (or luminosity).
First, we adopt our fiducial CLF parameterization and fit to the shifted B21+22-U and B21+22-D data sets.We find that the baseline model performs better at explaining B21+22-D (χ 2 ν = 1.47) but somewhat worse at matching B21+22-U (χ 2 ν = 2.04), compared to the unmodified case (χ 2 ν = 1.67).In these new fits all of the best-fit parameters vary only slightly except p which controls the faint-end behavior.In the up-shifted case the value lowers from p = 1.69 to p = 1.60 (corresponding to f (M ) ∝ M 0.47 at low mass) and in the downshifted case to p = 1.83 (f (M ) ∝ M 0.7 at low mass).The modified low mass trends for f (M ) can be compared with the baseline result of f (M ) ∝ M 0.56 .These modifications lead to only minor changes in our conclusions regarding reionization, shifting the timing of reionization by around ∆z = ±0.2 at fixed escape fraction, and the onset redshift of WF coupling by ∆z = ±0.9.
Next we turn to the effects the added normalization uncertainties have on our alternative faint-end models.If the lensed measurements are ∼ 20% overestimates of their true values, the evidence for an enhancement in star formation at low-mass is much less strong.For instance, although the flattening and shallow-slope models still perform well at explaining the data (χ 2 ν = 1.43 and χ 2 ν = 1.25 respectively), in an AIC sense they are less strongly preferred over the fiducial model fit to B21+22-D.For example, the shallow-slope case is still preferred (although now s = 0.35, so f (M ) ∝ M 0.35 , compared to s = 0.27 earlier), but the preference over the fiducial case drops to ∆AIC= 15.7 (from ∆AIC= 24.3).Yet, even accounting for potential systematic uncertainties in the lensed UVLF, the data seem to prefer some degree of SFE enhancement in low-mass halos relative to our fiducial case.
Alternatively, if the lensed measurements are ∼ 20% under-estimates of their true value this would yield stronger evidence of enhanced SFE at low mass.For instance, a shallowing slope of s = 0.19 would be preferred over the fiducial best fit at ∆AIC= 30.7 with respect to B21+22-U.
These shifts do not impact the conclusions regarding the truncation models.In all cases, the data disfavor strong suppressions in the SFE at low halo mass.Quantitatively, even the down-shifted data set disfavors 14 Some lower bounds are trimmed to forbid φ(m UV ) < 0 at 1−σ.
For measurements that are only 1−σ upper limits the upper limit is shifted instead.
strong suppression for halos of mass M ≥ 4 × 10 9 M at 2 − σ confidence.In summary, even accounting for this additional systematic uncertainty, our analysis indicates that a simple power law extrapolation from f at M ∼ 10 10 − 10 12 M may underestimate the SFE in M 10 10 M halos.The exact form of the SFE in such halos is not, however, wellconstrained by the current data.Further progress will likely require refinements in the lensed UVLF measurements and larger survey volumes.

CONCLUSION
The aim of this work has been to fit semi-empirical models to updated analyses of the UVLF in the HFF from B22, combined with measurements in the field from B21, together spanning z = 5 − 10.We explored the implications for our understanding of the SFE and its trends with halo mass and redshift, and the resulting consequences for our understanding of cosmic reionization, the redshifted 21 cm signal, and early JWST results.
Our main findings can be summarized as follows.First, a simple extrapolation of the best fit results to the B21 data continues to provide a reasonably good match to the lensed UVLF measurements from B22, even though these probe up to five magnitudes fainter in UV luminosity.The data, however, show a mild preference for an enhancement in the SFE at small masses relative to this simple extrapolation.Alternatively, the data may prefer an enhanced scatter in the correlation between UV luminosity and halo mass or a non-unity duty cycle at small mass scales.The combination of an enhanced average SFE and increased scatter at low masses is also plausible.However, the precise modification of the SFE in small mass halos is not yet well constrained after accounting for systematic and cosmic variance uncertainties in the lensed UVLF data.
Next, the lensed UVLF measurements probe the SFE in halos down to M ∼ 2×10 9 M , yet show no indication of a strong suppression on these mass scales from feedback effects.This constrains previously viable feedback models (see also B22).Further, the lensed UVLF measurements should also constrain departures from CDM in which the halo mass function is suppressed on small mass scales (e.g Schive et al. 2016).This will be the subject of upcoming work.
We also consider the implications of our UVLFcalibrated models for the EoR and Cosmic Dawn.Along the lines of previous work (e.g.Robertson et al. 2015;Sun & Furlanetto 2016), we find that we can reproduce Planck 18 τ e and current neutral fraction measurements provided the escape fraction of ionizing photons is f esc = 0.1 − 0.2.This scenario also roughly matches inferences of the ionizing emissivity from the Ly-α forest at z ∼ 5 − 6 (Becker et al. 2021), although the latest results from Gaikwad et al. (2023) prefer a slightly lower escape fraction of f esc = 0.05 − 0.1.The consistency between our models and reionization history measurements suggests that HST measurements have, in fact, identified the sources primarily responsible for reionization.Quantitatively, even in a model with a flat SFE in small mass halos the lensed UVLF measurements account for 50% of the total ionizing photon budget at z ≤ 9, with some caveats discussed in §5.3.In our fiducial scenario, WF coupling (related to the onset of the 21 cm absorption signal) occurs near z ∼ 13, while our flat SFE model gives an onset redshift (z ∼ 18) similar to that suggested by the possible EDGES 21 cm absorption signal (Bowman et al. 2018); see also Mirocha & Furlanetto (2019).
Finally, our HST-calibrated models agree well with early JWST results at z ≤ 10, but current UVLF estimates from JWST lie quite a bit above the HSTcalibrated models at higher redshift.This might be reconciled if the SFE reverses the decreasing trend with redshift suggested by the HST data, and increases at the higher redshifts only probed by JWST (z ≥ 10).However, dramatic evolution would be required.We discuss alternative resolutions further in §6.
In the near future, upcoming JWST lensed UVLF measurements should extend the reach of the lensed es-timates beyond z ≥ 10.This will provide further valuable information regarding the SFE at high redshifts and low halo masses, while early JWST candidates at z ≥ 10 have thus far been confined to the bright-end of the UVLF.In addition, we look forward to spectroscopic confirmations of photometric candidates from JWST, and further assessments of the purity of Lyman-break selected samples at z 10.The UVLF estimates can also be combined with Balmer line measurements from JWST, which provide a more direct handle on the ionizing emissivity, while upcoming line-intensity mapping surveys may probe the collective impact of individually faint galaxies (Bernal & Kovetz 2022).Further, we anticipate improved measurements of the reionization history itself from a wide variety of observational probes.Finally, additional comparisons between semi-empirical models and first-principles hydro-dynamical simulations of galaxy formation will also improve our understanding of the sources of reionization.

Figure 2 .
Figure 2. Posterior distributions in the five-dimensional parameter space of our CLF model fits from emcee.See Equation 8 and the accompanying text for a description of the model parameters, which describe the trend of the median UV luminosity with halo mass.The dashed black lines in the histogram panels indicate the 1-σ bounds on each parameter.The contours show the ∼12, 39, 68, and 86% confidence intervals in each 2D plane.The blue lines and squares show the best fit parameter values from scipy.optimize.The table in the upper right hand corner gives the correlation coefficients between the different parameter constraints.See the text for further discussion regarding the best fit parameter values and degeneracy directions.

Figure 3 .
Figure3.Mean UV luminosity-halo mass relationship implied by the CLF model fits at z = 6.The line and color styles are similar to those in Figure1and will be used in the figures throughout: the orange dashed curve is our model fit to the B21 measurements alone, while the blue solid curve also includes the B22 lensed UVLF measurements, and the green dotted curve is the alternative flattening scenario for the joint fit to B21+B22.Here we do not extrapolate beyond the brightest and dimmest measurements in the relevant data set.The corresponding star formation rates (from Equation5) are shown along the y-axis to the right.The B22 measurements probe down to halo mass scales as small as M ∼ 2 × 10 9 M in our fiducial scenario (with still smaller masses probed in the flattening case), about an order of magnitude smaller than possible with the B21 observations alone.

Figure 4 .
Figure4.Constraints on the SFE versus halo mass at z = 6 (solid) and z = 8 (dashed) for our three models.As in the previous figure, we show only the range of halo masses probed by each data set.In the fiducial scenario, the inferred f follows the form of Equation14where at z = 6 it peaks at f = 21% and a halo mass of M = 6 × 10 11 M .The shaded bands represent the effects of uncertainties in model parameters on our z = 6, B21+22 fit at the 1-, 2-, and 3-σ confidence levels.These constraints, while illustrative, only show the allowed range for the UV luminosity-halo mass relationship model of Equation8and do not capture the uncertainties in the functional form assumed here.The flattening scenario, hence, may lie outside of this band yet is also compatible with the current measurements.The red dotted line gives f (M, z = 6) (see text) for comparison.
adopt redshift independent values of C HII = 3 and T e = 10 4 K (see e.g.Sun & Furlanetto 2016 and references therein for a discussion).

Figure 5 .
Figure5.The reionization history from our best fit models, as compared to current measurements.The curves show the average ionization fraction as a function of redshift, and follow the usual color code of the previous figures.The solid curves assume that the escape fraction of ionizing photons is fesc = 0.2 while the dashed curves adopt fesc = 0.1.Reionization completes slightly earlier (i.e. at higher z) after including the B22 data since the faint-end luminosity function is marginally steeper in these observations, while reionization occurs at still higher redshifts in the flattening model.The data points show recent measurements fromMcGreer et al. (2015);Davies et al. (2018);Mason et al. (2018),Mason et al. (2019) andUmeda et al. (2023).

Figure 6 .
Figure6.Census of ionizing photons from current UVLF measurements.Left: The contribution to the ionizing photon product rate budget per logarithmic interval in halo mass, d ṅ/d ln M vs M .The color codes/line-styles for the curves are the same as in previous figures.The escape fractions are assumed to be independent of halo mass and set to fesc = 0.2, 0.2, 0.1 for B21+B22, B21, and the flattening model, respectively.In the fiducial CLF models, the peak contribution at z = 6 arises from halos with mass slightly larger than M ∼ 10 11 M and the peak mass scale shifts towards lower mass at increasing redshift.In the flattening model, lower mass halos also make an important (at z = 6) or even dominant contribution (e.g. in the z = 8 example).Right: The fraction of the ionizing photon rate production accounted for by current UVLF measurements as a function of redshift (assuming Mmin = 10 8 M ).In our fiducial model, the lensed UVLF measurements account for ≥ 80% of the ionizing photons at z 8.The field only (B21) observations are less complete, while the contributions from low luminosity galaxies below the reach of current surveys are more important in the flattening model.The untouched region beyond z 10 is being probed by JWST.

Figure 7 .
Figure 7.The ionizing emissivity as a function of redshift compared to inferences from Ly-α forest observations at z ≤ 6.The dotted cyan, dashed purple, and solid blue curves show the emissivity predicted in our fiducial B21+B22 model fit for fesc = 0.05, 0.1, and 0.2.The black squares with 1 − σ (black) and 2 − σ (grey) error bars are constraints from Becker et al. (2021) and red circles with 1 − σ (red) and 2 − σ (pink) error bars come from Gaikwad et al. (2023).Points at z = 5.1 and z = 6 that overlap have been slightly displaced in redshift for visual clarity.The Gaikwad et al. (2023) measurements prefer fesc ∼ 0.05 − 0.1 over higher escape fractions.

5. 5 .
Thomson Scattering Optical Depth 11 Strictly speaking, both measurements are sensitive only to the overall product of Nγ and fesc.Although we fix Nγ = 4, 000 and vary fesc around fesc ∼ 0.1 − 0.2, one can think of this step instead as exploring variations in the product of these quantities around Nγ fesc = 400 − 800.

Figure 8 .
Figure 8. Electron scattering optical depths at redshifts < z in our models compared to Planck 2018 τe measurements.Models are labeled with the usual color code, and escape fractions of fesc = 0.1(0.2) are shown as solid (dashed) curves.In the flattening model an escape fraction of fesc 0.15 is disfavored (at 2 − σ confidence), while fesc ∼ 0.1 − 0.2 are consistent with current measurements in the fiducial models.The best fit optical depth inferred from Planck 18 measurements is shown by the black dotted line, while the grey bands give the 1− and 2 − σ measurement errors (Planck Collaboration et al. 2018).

Figure 11 .
Figure11.Effect of alternative models on the SFE and faint-end UVLF.Top: The SFE at redshift z = 6 in our four cases of enhanced faint-end star formation compared with our fiducial best fit model (solid blue).The green dotted curve shows our usual flattening model, the red dashed case shows the best fit model with a shallower slope below M2, the purple dot-dashed curve allows a non-unity duty cycle at small masses, and the brown dash-dot-dotted case adopts a higher scatter at low masses.The red dashed shallowslope model is most preferred by the data.Bottom: The effect on the z = 6 UVLF from sharp truncations in star formation in halos of mass M < 5 × 10 9 M (red dashed) and M < 2.6 × 10 10 M (purple dotted, the mass scale below which we found a flattening in SFE is preferred) compared with our fiducial best fit model.The lensed UVLF strongly disfavors these and similar truncation models.