Rapid Reionization by the Oligarchs: The Case for Massive, UV-Bright, Star-Forming Galaxies with High Escape Fractions

TL;DR

This paper presents an empirical model showing that massive, UV-bright, star-forming galaxies with high escape fractions drove rapid cosmic reionization, emphasizing the role of a small subset of galaxies in this process.

Contribution

It introduces a model linking escape fraction to star-formation surface density and highlights the dominant role of bright, massive galaxies in reionization.

Findings

Massive, UV-bright galaxies with high escape fractions are key to reionization.

Less than 5% of galaxies account for over 80% of reionization.

Reionization occurs rapidly within 300 million years, driven by bright galaxies.

Abstract

The protagonists of cosmic reionization remain elusive. Faint star-forming galaxies are leading candidates because they are numerous and may have significant ionizing photon escape fractions ($f_{esc}$). Here we update this picture via an empirical model that successfully predicts latest observations (e.g., the drop in star-formation density at z>8). We generate an ionizing spectrum for each galaxy in our model and constrain $f_{esc}$ using latest measurements of the reionization timeline (e.g., Ly$\alpha$ damping of quasars and galaxies at z>7). Assuming a constant $f_{esc}$, we find $M_{UV}$<-13.5 galaxies need $f_{esc}=0.21^{+0.06}_{-0.04}$ to complete reionization. The inferred IGM neutral fraction is [0.9, 0.5, 0.1] at z=[8.2, 6.8, 6.2]$\pm$0.2, i.e., the bulk of reionization transpires in 300 Myrs. Inspired by the emergent sample of Lyman Continuum (LyC) leakers that overwhelmingly displays higher-than-average star-formation surface density ($\Sigma$), we propose a model relating $f_{esc}$ to $\Sigma$ and find $f_{esc}\propto\Sigma^{0.4\pm0.1}$. Since $\Sigma$ falls by ~2.5 dex between z=8 and z=0, our model explains the humble upper limits on $f_{esc}$ at lower redshifts and its required evolution to ~0.2 at z>6. Within this model, strikingly, <5% of galaxies with $M_{UV}$<-18 (the `oligarchs') account for >80% of the reionization budget. In fact, faint sources ($M_{UV}$>-16) must be relegated to a limited role to ensure high neutral fractions at z=7-8. Shallow faint-end slopes of the UV luminosity function ($\alpha$>-2) and/or $f_{esc}$ distributions skewed toward bright galaxies produce the required late and rapid reionization. We predict LyC leakers like COLA1 (z=6.6, $f_{esc}$~30%, $M_{UV}$=-21.5) become increasingly common towards z~6 and that the drivers of reionization do not lie hidden across the faint-end of the luminosity function, but are already known to us. (abridged)

Tables

Table 1. Table 1: Which galaxies reionized the universe? Mean properties weighted by contribution to reionization ( n ˙ ion subscript ˙ 𝑛 ion \dot{n}_{\mathrm{ion}} ) for Model II (Model I).

Parameter	$z=6$	$z=7$	$z=8$	$z=9$	$z=10$
	$f_{\mathrm{esc}}=1.62\times {\Sigma }_{\mathrm{SFR}}^{0.42}$ ($f_{\mathrm{esc}}$=0.21)
$f_{\mathrm{esc}}$	0.24 (0.21)	0.26 (0.21)	0.24 (0.21)	0.21 (0.21)	0.19 (0.21)
$M_{\mathrm{UV}}\text{<span><sup>\$a\$</sup><span><span><sup>\$a\$</sup><span>footnotemark: </span></span></span></span>}$	-19.7 (-18.9)	-19.5 (-18.6)	-19.2 (-18.2)	-18.7 (-17.8)	-18.3 (-17.3)
$\log \left(M_{\star }/M_{\odot }\right)$	9.5 (9.2)	9.2 (8.9)	8.9 (8.6)	8.6 (8.4)	8.3 (8.0)
$\log \left(M_{\mathrm{halo}}/M_{\odot }\right)$	11.1 (10.9)	10.9 (10.7)	10.7 (10.6)	10.5 (10.4)	10.4 (10.2)
$\mathrm{SFR}/M_{\odot }{\mathrm{yr}}^{-1}$	15.3 (9.3)	12.1 (6.5)	7.7 (4.0)	3.5 (2.1)	2.0 (1.0)
${\Sigma }_{\mathrm{SFR}}/M_{\odot }{\mathrm{yr}}^{-1}{\mathrm{kpc}}^{-2}$	17.7 (8.2)	23.2 (9.7)	17.6 (8.1)	11.6 (5.5)	9.3 (4.3)
$R_{\mathrm{gal}}/\mathrm{kpc}$	0.59 (0.62)	0.46 (0.48)	0.38 (0.39)	0.33 (0.34)	0.27 (0.27)

Table 2. Table 2: A Rapidly Reionizing Universe

Redshift	${\overline{x}}_{\mathrm{HI}}$ (Model I)	${\overline{x}}_{\mathrm{HI}}$ (Model II)
$6$	${0.00}_{-0.00}^{+0.05}$	${0.01}_{-0.01}^{+0.08}$
$7$	${0.58}_{-0.12}^{+0.08}$	${0.64}_{-0.04}^{+0.03}$
$8$	${0.87}_{-0.04}^{+0.03}$	${0.89}_{-0.01}^{+0.01}$
$9$	${0.95}_{-0.01}^{+0.01}$	${0.97}_{-0.01}^{+0.01}$
$10$	${0.99}_{-0.01}^{+0.01}$	${0.99}_{-0.01}^{+0.01}$

Table 3. Table 3: Comparison of Model UVLFs against literature data.

Redshift	${\chi }^2$ ()	${\chi }^2$ ($M_{\mathrm{UV}}>-16$)	Source
$6$	$0.7$	$0.5$	Bouwens et al. (2017)
	$9.9$	$0.1$	Atek et al. (2018)
	$1.5$	$1.7$	Livermore et al. (2017)
$7$	$1.4$	$-$	Bouwens et al. (2015a)
$8$	$1.3$	$-$	Bouwens et al. (2015a)
$10$	$0.9$	$-$	Oesch et al. (2018)

Equations

$${\xi }_{\mathrm{ion}}=\frac{N(H^0)}{L_{1500}}[{\mathrm{s}}^{-1}/erg{\mathrm{s}}^{-1}\mathrm{H}{\mathrm{z}}^{-1}].$$

$${\dot{n}}_{\mathrm{ion}}(z)=f_{\mathrm{esc}}{\xi }_{\mathrm{ion}}{\rho }_{{\mathrm{UV}}_{\mathrm{dust}\mathrm{corr}}}[{\mathrm{s}}^{-1}\mathrm{Mp}{\mathrm{c}}^{-3}],$$

$${\dot{n}}_{\mathrm{ion}}(z)=\sum _{M_{\mathrm{UV}}<-13.5}\frac{f_{\mathrm{esc}}N(H^0)}{10^6}[{\mathrm{s}}^{-1}\mathrm{Mp}{\mathrm{c}}^{-3}].$$

$${\dot{Q}}_{\mathrm{HII}}=\frac{{\dot{n}}_{\mathrm{ion}}}{⟨n_{\mathrm{H}}⟩}-\frac{{\bar{Q}}_{\mathrm{HII}}}{t_{\mathrm{rec}}},$$

$$1/t_{\mathrm{rec}}=C_{\mathrm{HII}}{\alpha }_{\mathrm{B}}\left(1+\left(1-X_{\mathrm{P}}\right)/4X_{\mathrm{P}}\right)⟨n_{\mathrm{H}}⟩{\left(1+z\right)}^3,$$

$$\tau (z)=c⟨n_{\mathrm{H}}⟩{\sigma }_T{\int }_0^z{f}_{\mathrm{e}}{\bar{Q}}_{\mathrm{HII}}H(z{)}^{-1}\left(1+z^{\prime 2}\right)d{z}^{\prime },$$

$${\Sigma }_{\mathrm{SFR}}=\frac{\mathrm{SFR}/2}{\pi R_{\mathrm{gal}}^2}.$$

$$f_{\mathrm{esc}}=min\left(1,1.6_{-0.3}^{+0.3}\times {\left(\frac{{\Sigma }_{\mathrm{SFR}}}{{\Sigma }_{\mathrm{SFR},\max }}\right)}^{0.4_{-0.1}^{+0.1}}\right).$$

$$f_{\mathrm{esc}}\approx 1-f_{\mathrm{cov}}$$

Full text

galaxies: high-redshift — galaxies: evolution — dark ages, reionization, first stars

††software: IPython (Pérez & Granger, 2007), matplotlib (Hunter, 2007), numpy (Oliphant, 2015), scipy (Jones et al., 2001–), jupyter (Kluyver et al., 2016), seaborn (Waskom et al., 2018), FSPS (Conroy et al., 2009, 2010), dynesty (Speagle, 2019), python-FSPS (Foreman-Mackey et al., 2014) ††thanks: Hubble Fellow

1 Introduction

The Epoch of Reionization (EoR) marks the last great phase transition of the universe, during which vast islands of neutral Hydrogen were ionized by the first sources of light (Loeb & Barkana, 2001). The protagonists, topology, and timeline of the EoR are intertwined with our understanding of the early universe and its newly born stellar populations (for a recent review, see Dayal & Ferrara, 2018). Due to the rapidly fading quasar emissivity at $z>3$ (e.g., Kulkarni et al., 2019), star-forming galaxies are favored to dominate reionization (e.g., Bouwens et al., 2015a). Bright star-forming galaxies have not shown much promise of being effective ionizing sources. Until very recently, these galaxies were measured to have negligible ionizing photon escape fractions (e.g. Steidel et al., 2018). This, combined with their observed rarity has meant a reservoir of ultra-faint sources far below current detection limits (modulo highly lensed fields) is widely invoked to drive reionization (e.g., Livermore et al., 2017).

Modeling reionization by star-forming galaxies is typically cast as a tale of three quantities: $\rho_{\mathrm{SFR}}$, $\xi_{\mathrm{ion}}$, and f_{\mathrm{esc}}\(e.g., Madau et al., 1999; Robertson et al., 2015; Bouwens et al., 2015a). The cosmic star-formation rate density, $\rho_{\mathrm{SFR}}$, provides a measure of star-formation in the early universe. It has now been tracked out to $z\sim 10$, with latest studies showing an accelerating drop beyond $z>8$ (Oesch et al., 2018; Ishigaki et al., 2018), consistent with the predictions of simple models that link star-formation rates (SFR) to dark matter accretion (e.g., Tacchella et al., 2013, 2018; Mason et al., 2015; Mashian et al., 2016). The ionizing photon production efficiency, $\xi_{\mathrm{ion}}$, is a conversion factor for how many Hydrogen ionizing photons emerge from each episode of star-formation. Tight constraints on $\xi_{\mathrm{ion}}$ can be placed using H$\alpha$ measurements and some assumptions about $f_{\mathrm{esc}}$, and now exist from direct spectroscopy (Nakajima et al., 2016; Matthee et al., 2017; Shivaei et al., 2018; Tang et al., 2019) and IRAC-excess inferences of H$\alpha$ out to $z\sim 5$ (Bouwens et al., 2016a; Lam et al., 2019).

The escape fraction, $f_{\mathrm{esc}}$, is the fraction of ionizing photons generated in a galaxy that evade photoelectric absorption and dust in the Interstellar Medium (ISM) to escape into the neutral Intergalactic Medium (IGM) and ionize it. While $\rho_{\mathrm{SFR}}$ and $\xi_{\mathrm{ion}}$ will be measured ever more precisely with, e.g., the James Webb Space Telescope (JWST), ionizing radiation and thus f_{\mathrm{esc}}\will never be directly measured in the EoR due to the opacity of the intervening neutral IGM (e.g., Fan et al., 2001; Inoue et al., 2014; McGreer et al., 2015). To make things more challenging, it is also extremely difficult to get a handle on f_{\mathrm{esc}}\through simulations since it depends sensitively on resolving the multi-phase ISM, and the treatment of small-scale processes associated with galaxy formation, which, at present, are modelled in an approximate way (e.g., Wise et al., 2014; Ma et al., 2015; Trebitsch et al., 2017).

However, there is a path forward: since f_{\mathrm{esc}}\is by far the single largest uncertainty in modeling reionization, we can employ $\rho_{\mathrm{SFR}}$ and $\xi_{\mathrm{ion}}$ from state-of-the-art measurements and constrain f_{\mathrm{esc}}\against latest measurements of the timeline of reionization. The fraction of dark pixels in the Ly$\alpha$ and Ly$\beta$ forests provides a model-independent limit on the end of reionization (McGreer et al., 2015). The electron scattering optical depth ($\tau$) reported by Planck Collaboration et al. (2018), much lower and more precise than previous measures of $\tau$ (e.g., Hinshaw et al., 2013), is an integrated probe of the density of ionizing photons, as CMB photons scatter off of electrons knocked out of Hydrogen atoms. The first quasars and large Ly$\alpha$ surveys at $z\gtrsim 7$ allow detailed inferences of the neutral fraction as a function of redshift from the magnitude of Ly$\alpha$ damping (e.g., Bañados et al., 2018; Mason et al., 2018).

Complementing these data on the global history of reionization are clues about f_{\mathrm{esc}}\on a galaxy by galaxy level. For the first time we have a robust sample of individual star-forming galaxies securely detected in Lyman Continuum (LyC) spanning $z\sim 0.3-4$ (“LyC leakers”) (e.g., Naidu et al., 2017; Vanzella et al., 2018; Rivera-Thorsen et al., 2019). The campaigns targeting Green Peas at $z\sim 0.3$ with $HST$/COS have proven remarkably efficient (boasting a $100\%$ success rate for LyC detection; Izotov et al., 2016b, 2018b). Concurrently, a sample at high-$z$ is emerging from deep rest-frame UV spectroscopy (e.g., Steidel et al., 2018) and imaging with HST/UVIS (e.g., Fletcher et al., 2019). Taken together, these leakers provide hints about the galaxy properties that favor LyC leakage during the EoR. For instance, the overwhelming majority of LyC leakers are compact (e.g., Izotov et al., 2018a) and show multi-peaked Ly$\alpha$ (e.g., Verhamme et al., 2017). These insights can be incorporated into models of f_{\mathrm{esc}}\that improve upon previous analyses that assumed a single number across the entire galaxy population.

Meanwhile, the bulk of observational constraints on the average LyC f_{\mathrm{esc}}\have relied on stacking shallow non-detections for individual galaxies to place stringent upper-limits of f_{\mathrm{esc}}$$<10\% out to $z\sim 4$ (e.g., Siana et al., 2010; Rutkowski et al., 2017; Grazian et al., 2017; Japelj et al., 2017; Naidu et al., 2018). Taken at face value, these studies effectively rule out an average f_{\mathrm{esc}}$$>10\% for M_{\mathrm{UV}}$$\lesssim-20 sources and put the focus on fainter galaxies, for which no LyC constraints exist yet, as the drivers of reionization. However, if we consider the anisotropy, stochasticity, and evolution with $z$ of f_{\mathrm{esc}}\that recent simulations have brought to light (e.g., Paardekooper et al., 2015; Trebitsch et al., 2017; Rosdahl et al., 2018) along with a higher CGM+IGM opacity, the limits from these studies are far less stringent and can be relaxed by factors of $2-5\times$ (e.g., Fletcher et al., 2019; Steidel et al., 2018). And indeed, latest studies that emphasize deep spectra and photometry for individual sources find f_{\mathrm{esc}}$$\sim 10\% in stacks of normal $\log(M_{\star}/M_{\odot})\sim 8.5-10$ galaxies at $z\sim 2.5-4$ (Oesch et al. in prep., 2019; Steidel et al., 2018; Marchi et al., 2018). The emerging observational picture is that average f_{\mathrm{esc}}\of $\sim$10$\%$ are possible in relatively bright galaxies (M_{\mathrm{UV}}$$<-20).

In parallel, early hydrodynamical simulations produced a mixture of results: with f_{\mathrm{esc}}\correlating with halo mass (e.g., Gnedin et al., 2008; Wise & Cen, 2009), anti-correlating with halo mass (e.g., Yajima et al., 2011; Paardekooper et al., 2015; Kimm et al., 2017), and with typical time-averaged values far smaller than f_{\mathrm{esc}}$$\sim 20\% (e.g., Ma et al., 2015). However, contemporary hydrodynamical simulations through a combination of feedback, binaries, turbulence, and careful modeling of the multi-phase ISM are able to produce average f_{\mathrm{esc}}\in the $10-20\%$ range in $\gtrsim 10^{8}~{}M_{\odot}$ galaxies (e.g., Ma et al., 2016; Rosdahl et al., 2018). Ma et al. (2015) and Ma et al. (2016) are particularly illustrative of this shift, where these authors first wrote about the difficulty of producing f_{\mathrm{esc}}\$$>5\% due to high absorption in the birth clouds of massive stars but then subsequently found binary models of stellar evolution that destroyed these clouds while retaining highly ionizing sources until late times could achieve f_{\mathrm{esc}}$$\sim 20\%.

The driving impulse of this work is to unite the developments outlined above self-consistently under the same umbrella to see what story they tell about f_{\mathrm{esc}}\and thus reionization. Our umbrella of choice is the empirical galaxy formation model by Tacchella et al. (2018) that incorporates recent developments (e.g., cutting-edge stellar population synthesis models) and accurately predicts latest observations (e.g., the sharp drop in $\rho_{\mathrm{SFR}}$) (§2.1). Leaving f_{\mathrm{esc}}\as a free parameter in the equations of reionization (§2.2), we fit for it against recently derived constraints on reionization that we describe in §2.3. Two models of f_{\mathrm{esc}}\– one constant across all galaxies during reionization, another dependent on star-formation surface density – are justified, set up, and fit in §3 and §4. The implications of the resulting reionization histories – their rapid pace, the concentration of the reionization budget among “oligarch” galaxies, the path forward for observational studies – are discussed in §5. We address open questions and caveats in §6. A summary of our findings and an outlook to the future is presented in §7.

We use f_{\mathrm{esc}}\to denote both the singular and plural “escape fraction” and “escape fractions”. The volume-averaged IGM neutral fraction and ionized fraction are denoted by $\bar{x}_{\mathrm{HI}}$ and $\bar{Q}_{\mathrm{HII}}=1-\bar{x}_{\mathrm{HI}}$ respectively. For cosmological parameters, we adopt the following from Planck Collaboration et al. (2018): $h=0.6772$, $X_{\mathrm{p}}=0.75328$, $\Omega_{\mathrm{b}}=0.02241/h^{2}$, $\rho_{\mathrm{c}}=1.8787\times 10^{-29}h^{2}$. All magnitudes are in the AB system (Oke & Gunn, 1983).

2 Methods

2.1 An Empirical Model for Galaxy Evolution at $z\geq 4$

The foundation of this work is the empirical model introduced in Tacchella et al. (2018). Here we briefly summarize it, and then describe in detail the quantities relevant to reionization.

2.1.1 Model Description

The Tacchella et al. (2018) model is built on top of a $\mathrm{10^{6}~{}Mpc^{3}}$, high-resolution, N-body, dark-matter simulation, color (Sawala et al., 2016; Hellwing et al., 2016). It makes the assumption that the star-formation rate (SFR) of a halo depends on the growth rate of a halo and a star-formation efficiency that is independent of redshift (see also Tacchella et al., 2013; Mason et al., 2015). The halo merger trees self-consistently give rise to star-formation histories for each galaxy from which spectral energy distributions (SEDs) are computed using the Flexible Stellar Population Synthesis code (FSPS, Conroy et al., 2009, 2010; Foreman-Mackey et al., 2014) and the MESA Isochrones and Stellar Tracks which incorporate the effects of rotation (MIST, Dotter, 2016; Choi et al., 2016, 2017).

The star-formation efficiency of the model is tuned to match the $z=4$ Ultraviolet Luminosity Function (UVLF) and then predicts UVLFs out to $z\sim 10$ consistent with the observed data (see Figure 1 in this work and Figure 3 in Tacchella et al. 2018). The faint-end slope of the UVLF ($\alpha_{\mathrm{UV}}$) in our model steepens with redshift and the same trend is seen in observations (e.g., Finkelstein et al., 2015; Bouwens et al., 2015b). For our best-fit Schechter function at $z=$[4, 6, 8, 10, 12] we find $\alpha_{\mathrm{UV}}=$[$-1.63\pm 0.02$, $-1.72\pm 0.03$, $-1.69\pm 0.04$, $-1.84\pm 0.06$, $-1.99\pm 0.24$], consistent with recent empirical and semi-analytical models that find $\alpha_{\mathrm{UV}}=-1.5$ to $-2.0$ at $z=5-10$ (e.g., Yung et al., 2019; Endsley et al., 2019; Behroozi et al., 2019). Our best-fit faint-end slope values are somewhat shallower than what is reported in the recent literature from observations (Bouwens et al., 2017; Livermore et al., 2017; Atek et al., 2018). However, we note the Schechter function parameters are highly degenerate and that our slopes are fit over a different range of M_{\mathrm{UV}}\(we truncate at M_{\mathrm{UV}}\$$<-13.5) compared to observations. Figure 1 shows the actual LF from our model is in decent agreement with the latest observations at $z\sim 6$, down to the faintest limits that can currently be probed in the Hubble Frontier Fields. We provide UVLFs for our model out to $z=12$ and $\chi^{2}$ comparisons against literature data in Appendix A. Note that we compute the ionizing budget via SEDs of individual sources, and not by integrating under the UVLF. We discuss $\alpha_{\mathrm{UV}}$ in detail in §5.3.

The resolution of the dark-matter simulation limits our model to $\mathrm{SFR\gtrsim 0.02\ M_{\odot}/yr}$, corresponding to $\mathrm{M_{UV}\lesssim-13.5}$, roughly where estimates of the faint end of the UVLF begin to diverge due to magnification uncertainties of the lensing models (e.g., Bouwens et al., 2017). In all calculations we integrate down to this limit. Our fiducial model adopts a Salpeter (1955) IMF, a constant metallicity of $\mathrm{Z_{\star}/Z_{\odot}=0.02}$, and the UV continuum dust prescription of Bouwens et al. (2014) (their Table 3). Swapping the MIST models for those that explicitly include binaries (BPASS, Stanway & Eldridge, 2018), or a Chabrier (2003) IMF, or a model with evolving metallicity make no appreciable difference to our results. Dust, the faint-end slope of the UVLF, and the effect of changing the M_{\mathrm{UV}}\cutoff are discussed in §6.

2.1.2 The Ionizing Photon Production Efficiency ($\xi_{\mathrm{ion}}$)

$\xi_{\mathrm{ion}}$ provides a measure of the gross LyC photons produced at a given time in a source. It is typically cast in terms of the rate of ionizing photons ($\mathrm{N(H^{0})}$) per unit UV luminosity, usually measured at $1500\textup{\AA}$ ($\mathrm{L_{1500}}$):

[TABLE]

We compute $\xi_{\mathrm{ion}}$ for each galaxy in the empirical model directly from its SED by integrating the flux produced below the Lyman limit to obtain $N(H^{0})$ and then normalizing by the SED-flux at $\mathrm{1500\textup{\AA}}$. The MIST isochrones that our SEDs are synthesized from include the effects of rotation that boost the ionizing flux production of massive stars akin to, but not exactly like, the effect of binaries (Choi et al., 2017). The harder UV spectra produced by these rotating models (or binaries) are the only kind that self-consistently explain the strong nebular fluxes, high ionization lines, and extreme line widths that are now known to be ubiquitous at $z\sim 2.5-3.5$ (e.g., Steidel et al., 2016; Holden et al., 2016; Reddy et al., 2018b) and also commonly seen at $z>6$ (e.g., Roberts-Borsani et al., 2016; Stark et al., 2015; Mainali et al., 2017).

In Figure 2 we show the $\log_{10}(\xi_{\mathrm{ion}})$ distribution for the galaxies in our model. We predict the median $\xi_{\mathrm{ion}}$ between $z=4-10$ rises by $\sim 40\%$ ($\sim 0.15$ dex) as galaxies get younger at higher redshift and their ionizing spectra become harder. We compare our predictions with Bouwens et al. (2016a), who report a mean $\log_{10}(\xi_{\mathrm{ion}})=25.34^{+0.02}_{-0.02}\ (25.54^{0.12}_{-0.12})$ for a sample spanning $z=3.8-5.0\ (5.0-5.5)$ using an SMC attenuation curve, which Reddy et al. (2018a) show to be the appropriate curve for sub-solar metallicity populations expected at $z>4$. At $z=4$ ($z=5$) we have a median $\mathrm{log_{10}(\xi_{\mathrm{ion}})=25.37^{+0.06}_{-0.06}\ (25.40^{+0.05}_{-0.06}}$) that agrees with their measurements within error-bars. At fixed redshift we find $\xi_{\mathrm{ion}}$ does not vary with mass or $M_{\mathrm{UV}}$. Lam et al. (2019) observe a similar invariance with brightness at $z\sim 4-5$ for M_{\mathrm{UV}}$$<-17.5 galaxies. These trends hold even when using BPASS templates and with an evolving metallicity. While not directly dealing with the redshifts that are the focus of this work, we note that the $\xi_{\mathrm{ion}}$ values reported for $z\sim 2-3$ galaxies are in broad agreement with our model assuming a linear extrapolation to lower redshifts (Nakajima et al., 2016; Matthee et al., 2017; Shivaei et al., 2018).

Typically, reionization studies set $\xi_{\mathrm{ion}}$ to some fixed, redshift-invariant value that lies in the locus our galaxies span in Figure 2. This is a reasonable assumption as evidenced by the narrow spread ($\sim 0.1$ dex at each redshift across $z=4-10$) and gradual evolution during reionization ($z=6-10$). What this means is that when considering the total integrated ionizing output at a particular redshift, there is not much of a distinction between our approach and assuming some reasonable fixed value (note how close the dotted and dashed purple lines are in Figure 3). However, the advantage of our model is that it captures the diversity in $\xi_{\mathrm{ion}}$ on a galaxy by galaxy basis so that we are able to link f_{\mathrm{esc}}\to individual galaxy properties and through the product f_{\mathrm{esc}}$$\times N(H^{0}) probe how much each galaxy contributes to reionization (see also §4, Figure 9, and Table 1).

2.1.3 Cosmic Star-Formation Rate Density ($\mathrm{\rho_{SFR}}$)

At $z=4-8$ the Tacchella et al. (2018) model is in excellent agreement with the consensus $\mathrm{\rho_{SFR}}$ (e.g., Bouwens et al., 2015b; Finkelstein et al., 2015; Madau & Dickinson, 2014). At $z>8$, where various measurements diverge, the model predicts a drop in $\rho_{\rm SFR}$ consistent with the latest HST analyses from Ishigaki et al. (2018) and Oesch et al. (2018). The sharp drop in $\rho_{\mathrm{SFR}}$ in our model comes as the bulk of halos at $z>8$ begin to fall below the halo-mass corresponding to maximal star-formation efficiency ($\mathrm{M_{h}\sim 10^{11}-10^{12}~{}M_{\odot}}$).

The difference between earlier smooth power-law fits for $\mathrm{\rho_{SFR}}$ at $z>4$, which use steeper $\alpha_{\mathrm{UV}}\lesssim-2$ (e.g., $\mathrm{\rho_{SFR}\propto}(1+z)^{-4.2}$; Bouwens et al. 2015b; Finkelstein et al. 2015; Robertson et al. 2015) and our model, which predicts an accelerating decline in $\mathrm{\rho_{SFR}}$, is as large as an order of magnitude at $z=10$ and three orders at $z=14$. This difference is directly reflected in the dearth of LyC photons available for reionization at early times. Earlier works (e.g., Robertson et al., 2015) had already shown that reionization likely proceeds without significant contribution from $z>10$ sources. The dearth of sources in our model (and thus LyC) at $z\sim 8-10$ combined with other data, as we shall see in §5.1, even further compresses the timeline of reionization and pushes it to later times.

2.2 Equations of Reionization

We closely follow the widespread approach that models reionization as an interplay between ionization and recombination (e.g., Madau et al., 1999; Robertson et al., 2013). Here we outline the relevant equations.

We start with the quantity directly inherited from our empirical model: the co-moving production rate of hydrogen-ionizing photons ($\dot{n}_{ion}$), i.e., the gross number of LyC photons escaping into the IGM per unit time per unit volume. $\dot{n}_{ion}$ is usually computed as follows:

[TABLE]

where $\rho_{\mathrm{UV}_{\mathrm{dust\ corr}}}$ is the dust-corrected UV luminosity.

This UV-anchored formulation, where $f_{\mathrm{esc}}$ is the escape fraction of ionizing photons relative to the dust-corrected observed UV, is apt for working with observations, where $\rho_{\mathrm{UV}}$ is the measured quantity around which all else is based. However, in our model which is built on a $\mathrm{10^{6}~{}Mpc^{3}}$ simulation box, we simply sum the LyC photons, $N(H^{0})$, produced by every galaxy from their SEDs directly, reducing Equation 2 to:

[TABLE]

Here $f_{\rm esc}$ is the escape fraction of ionizing photons relative to the total ionizing photons produced in the galaxy.

The IGM ionized fraction, $\bar{Q}_{\mathrm{HII}}$, is evolved as per the following differential equation where the first term represents ionization, and the second recombination,

[TABLE]

where $\langle n_{\mathrm{H}}\rangle=X_{\mathrm{p}}\Omega_{\mathrm{b}}\rho_{\mathrm{c}}$ is the co-moving density of Hydrogen, which depends on the primordial mass-fraction of Hydrogen ($X_{\mathrm{p}}$), the fractional baryon density $\Omega_{\mathrm{b}}$ and the critical density $\rho_{\mathrm{c}}$. $t_{\mathrm{rec}}$, the recombination time of ionized Hydrogen in the IGM, is given by

[TABLE]

where $C_{\mathrm{HII}}=\langle n_{\mathrm{H}}^{2}\rangle/\langle n_{\mathrm{H}}\rangle^{2}$ is the clumping factor that models the inhomogeneity of the IGM which we set to 3, and $\alpha_{\mathrm{B}}=2.6\times 10^{-13}\left(T/10^{4}\mathrm{K}\right)^{0.76}\mathrm{cm^{3}s^{-1}}$ is the Case B recombination coefficient at electron temperature, $T$ that we set to $10^{4}\mathrm{K}$ (Shull et al., 2012; Robertson et al., 2013; Pawlik et al., 2015; Robertson et al., 2015; Sun & Furlanetto, 2016).

The Thomson optical depth, $\tau$, is calculated as

[TABLE]

where $\tau$ is the Thomson optical depth, $c$ is the speed of light, $\sigma_{T}$ is the Thomson scattering cross-section, $f_{e}$ is the number of free electrons for every Hydrogen nucleus in the ionized IGM that we set to $\left(1+(1-X_{\mathrm{P}})/4X_{P}\right)$ (Kuhlen & Faucher-Giguère, 2012), and $H(z)$ is the Hubble parameter.

We note there are caveats: for example, Case B recombination may not be an appropriate description towards the end of reionization when local absorption in dense clumps becomes important (Furlanetto & Oh, 2005), or $C_{\mathrm{HII}}$ is bound to evolve as the universe grows more ionized (Pawlik et al., 2015) though its effect on reionization inference is limited (Mason et al., 2019a; Bouwens et al., 2015a). While we could test the effect of assuming different values for every individual parameter (on top of varying the IMF, metallicity, underlying SED models, and the truncation $M_{\mathrm{UV}}$), our guiding philosophy is to hew to the canonical assumptions from the literature so all our divergent conclusions are clearly attributable to the new data we constrain against and the models we introduce in this work.

The only free parameter in these equations is the escape fraction, $f_{\mathrm{esc}}$. Our strategy to constrain $f_{\mathrm{esc}}$, and thus constrain reionization histories is to solve Equations 3-6 assuming a model for f_{\mathrm{esc}}\and then fit the model parameters against the data described in the following section.

2.3 Observational Constraints on Reionization

Here we enumerate state-of-the-art measurements on the timeline of reionization that we use to constrain our f_{\mathrm{esc}}\models. We briefly describe each measurement and specify how it is included in our inference.

Thomson Optical Depth: $\tau=0.0540\pm 0.0074$ from Planck Collaboration et al. (2018), which is a more precise, downward revision of $\tau=0.066\pm 0.012$ (Planck Collaboration et al., 2016) and a far cry from the WMAP $\tau=0.088\pm 0.014$ (Hinshaw et al., 2013). $\tau$ bears the imprint of free electrons on photons from the last-scattering surface of the CMB and provides a global, integrated, model-independent constraint. The lower value of $\tau$ from Planck Collaboration et al. (2018) allows for the sharp drop in $\rho_{\mathrm{SFR}}$ at $z>8$ (Figure 3) which was disfavored by earlier measurements (e.g., Robertson et al., 2013). 2. 2.

Ly$\alpha$, Ly$\beta$ dark fraction: $\bar{x}_{\mathrm{HI}}\leq 0.06\pm 0.05$ at $z=5.9$, as per the model-independent “dark fraction” in Ly$\alpha$ and Ly$\beta$ forests of quasar spectra (Mesinger, 2010; McGreer et al., 2015). The completely dark pixels in the forests are either due to neutral H I in the IGM and/or astrophysical interlopers and hence provide an assumption-free upper limit to the global neutral fraction. This constraint allows us to impose the “end of reionization” in a more self-consistent fashion than abrupt truncation at some redshift (often fixed to $z=6$ (e.g., Planck Collaboration et al., 2016, 2018)). We adopt it as uniform for $\bar{x}_{\mathrm{HI}}<0.06$ and as a half-Gaussian peaked at $0.06$ with $\sigma=0.05$ elsewhere (Greig & Mesinger, 2017, §3.1). 3. 3.

$z\sim 7-8$ Ly$\alpha$ Equivalent Width (EW) Distributions: $\bar{x}_{\mathrm{HI}}=0.59^{+0.11}_{-0.15}$ at $z\sim 7$, $\bar{x}_{\mathrm{HI}}=0.88^{+0.05}_{-0.10}$ at $z\sim 7.5$, and $\bar{x}_{\mathrm{HI}}>0.76$ at $z\sim 8$ (Mason et al., 2018; Hoag et al., 2019; Mason et al., 2019b). The Ly$\alpha$ line – in particular, its EW and its velocity offset from the systemic redshift – bears the imprint of the neutral IGM that Mason et al. (2018) infer using empirical fits for the ISM, and state-of-the-art IGM and Ly$\alpha$ radiative transfer simulations (Mesinger et al., 2016). While the evolution in the fraction of Ly$\alpha$ emitters in Lyman-break galaxies (e.g., Mesinger et al., 2015) encodes the evolution of $\bar{x}_{\mathrm{HI}}(z)$, Mason et al. (2018) show that a more competitive constraint may be derived by utilizing EW distributions. We use their full posterior PDF on $\bar{x}_{\mathrm{HI}}$ (their Figure 11) and adopt scatter in the redshift at which they report $\bar{x}_{\mathrm{HI}}$ as per the selection function for their sample – centred on $z=6.9$ and $\sigma_{z}=0.5$ (Grazian et al., 2012; Pentericci et al., 2014), which is consistent with the scatter for similarly selected z-band dropouts (e.g., Bouwens et al., 2015b). Mason et al. (2019b) and Hoag et al. (2019) apply the same technique at higher redshifts and we adopt their measurements in similar fashion. 4. 4.

$z>7$ Quasi-stellar Objects (QSOs): $\bar{x}_{\mathrm{HI}}=0.48^{+0.26}_{-0.26}$ at $z=7.09$ and $\bar{x}_{\mathrm{HI}}=0.60^{+0.20}_{-0.23}$ at $z=7.54$ (Davies et al., 2018a) from the IGM Ly$\alpha$ damping wing signature (Miralda-Escudé, 1998) in the quasars ULAS J1120+0641 (Mortlock et al., 2011) and ULAS J1342+0928 (Bañados et al., 2018). This constraint arises from a similar approach as Mason et al. (2018) in that, detailed empirical models (Davies et al., 2018b), IGM simulations (Mesinger et al., 2011) and radiative transfer (Davies et al., 2016) inform the inference of $\bar{x}_{\mathrm{HI}}$ from quasar spectra. We adopt their conservative PDF for $\bar{x}_{\mathrm{HI}}$ (their Figure 11). 5. 5.

$z=6-7$ Ly$\alpha$ Emitter (LAE) Fraction: The drop in number-density of LAEs between $z\sim 6$ and $z\sim 7$ may be interpreted as the universe undergoing drastic evolution in neutrality between these redshifts (Mesinger et al., 2015) but may also be due to survey incompleteness at faint magnitudes (Oyarzún et al., 2017). Greig & Mesinger (2017) conservatively quantify this as implying $\bar{x}_{\mathrm{HI}}(z=7)-\bar{x}_{\mathrm{HI}}(z=6)\geq 0.4$ and we adopt their weak half-Gaussian constraint peaked at $\bar{x}_{\mathrm{HI}}=1$ with $\sigma=0.6$. Note that the dark fraction constraint at $z=5.9$ along with the $z\sim 7$ Ly$\alpha$ EWs already effectively reproduce this sharp change in neutrality. 6. 6.

$z\sim 6.6$ Ly$\alpha$ Emitter Clustering: The observed LAE clustering function at $z\sim 6.6$ (Ouchi et al., 2010) when interpreted in the context of the clustering in detailed reionization simulations suggests $\bar{x}_{\mathrm{HI}}<0.5$ (Sobacchi & Mesinger, 2014) which we implement as a half-Gaussian peaked at zero with $\sigma=0.5$ (Greig & Mesinger, 2017).

We exclude some measurements: Gamma Ray Burst (GRB) damping spectra, while probing some of the highest redshifts (Chornock et al., 2013; Totani et al., 2006; Tanvir et al., 2009), preferentially arise out of low-mass halos (Savaglio et al., 2009). Given the extreme scatter along the lines of sight to such halos across a patchy IGM, bounds on $\bar{x}_{\mathrm{HI}}$ are fated to be weak (McQuinn et al., 2008). Ly$\alpha$-based measurements that do not vary the intrinsic line-width (e.g., Ouchi et al., 2010; Inoue et al., 2018) are likely optimistic, given the vast diversity of ISM conditions in galaxies evident in line-shapes already seen in $z\sim 2-3$ samples (e.g., Trainor et al., 2015, 2016).

In what follows we set up and justify two models for f_{\mathrm{esc}}\that we fit against all the constraints described in this Section.

3 Fitting for $f_{\mathrm{esc}}$ Model I: Constant $f_{\mathrm{esc}}$ During Reionization

Here we assume the $f_{\mathrm{esc}}$ of all galaxies during reionization to be a constant number and denote this as “Model I”. Effectively, we fit for a single normalization factor, $f_{\mathrm{esc}}$, that sets the scale of the emissivity (solid curve in Figure 3). This is the common approach adopted in several reionization studies (e.g. Robertson et al., 2015; Ishigaki et al., 2018). Model I ignores the diversity of galaxies and the highly likely dependence of f_{\mathrm{esc}}\on various galaxy properties. However, this simple model provides a useful benchmark for the “average” escape fraction that observational stacking studies compute. Further, intrinsic galaxy properties (e.g., sizes, average star-formation rates) evolve modestly between $z=6-10$ where the bulk of reionization is expected to occur, hence assuming a constant average is justified.

We assume a uniform prior between 0 and 1 on f_{\mathrm{esc}}\and depict the resulting posteriors projected in various spaces in Figure 4. Combining all constraints we find $f_{\mathrm{esc}}$=$0.21^{+0.06}_{-0.04}$. Simply requiring reionization to be mostly complete by $z=5.9$ via the dark fraction rules out the f_{\mathrm{esc}}$$\lesssim 15\% parameter space (upper left panel of Figure 4). f_{\mathrm{esc}}$$\lesssim 15\% is also disfavored by the Planck $\tau$. Note that as the dark fraction and $\tau$ are model-independent constraints not much can be invoked to allow for f_{\mathrm{esc}}$$\lesssim 15\% (we discuss $M_{\mathrm{UV}}$ truncation and the faint-end slope of the UVLF in §6.2 and §5.3). The most constraining measurements on f_{\mathrm{esc}}\prove to be from the damping wing analysis of quasars and Ly$\alpha$ EW distributions which both require significant neutral fractions at later times ($\bar{x}_{\mathrm{HI}}\sim 0.5$ at $z\sim 7$).

In this constant f_{\mathrm{esc}}\model we make no claims about the f_{\mathrm{esc}}\at $z<6$ – our result is situated in the reionizing universe. f_{\mathrm{esc}}$$=0.21 is larger than the negligible f_{\mathrm{esc}}\measured in deep stacks at $z\sim 0-1$ (e.g., Siana et al., 2010; Rutkowski et al., 2016), where the IGM does not impede observations, and the recently established f_{\mathrm{esc}}$$\sim 10\% at $z\sim 2.5-4$ (Oesch et al. in prep., 2019; Marchi et al., 2017; Steidel et al., 2018; Fletcher et al., 2019).

To self-consistently bridge these findings of f_{\mathrm{esc}}$$\sim 0\% at $z\sim 0$, f_{\mathrm{esc}}$$\sim 10\% at $z\sim 2$, and f_{\mathrm{esc}}$$\sim 20\% at $z>6$ we introduce Model II below, which accounts for an evolving f_{\mathrm{esc}}\while also considering the diversity in properties of individual galaxies.

4 Fitting for $f_{\mathrm{esc}}$ Model II: $f_{\mathrm{esc}}$ as a function of $\Sigma_{\mathrm{SFR}}$

Here we propose a model where f_{\mathrm{esc}}\for each galaxy is solely dependent on its star-formation rate surface density $\Sigma_{\rm SFR}$, $f_{\mathrm{esc}}$= $a\times\Sigma_{\mathrm{SFR}}^{b}$ (where $a$ and $b$ are free parameters which we fit). We justify why this is an apt formulation, specify how it is implemented in our empirical model, and discuss the constraints it yields.

4.1 Motivation: Why $\Sigma_{\mathrm{SFR}}$?

Almost all the individual observed LyC leakers to date spanning $z\sim 0-6.6$ show $\Sigma_{\mathrm{SFR}}$ higher than the average $\Sigma_{\mathrm{SFR}}$ expected at their redshifts. We demonstrate this in Figure 5 where we have compiled all galaxies for which convincing LyC leakage is reported, and that have sizes and UV SFRs available. These include the HST/COS sample at $z\lesssim 0.3$ (Borthakur et al., 2014; Heckman et al., 2011; Izotov et al., 2016a, b, 2018a, 2018b), the HST/F336W, HST/F275W, and ground-based UV spectrograph sources at $z\sim 2-4$ (de Barros et al., 2016; Vanzella et al., 2016; Shapley et al., 2016; Bian et al., 2017; Naidu et al., 2017; Vanzella et al., 2018, Naidu et al. in prep.), and sources that show strong indirect hints of LyC escape at $z\sim 4-6.6$ (low covering fractions: Jones et al. 2013; Leethochawalit et al. 2016, tightly-spaced, double-peaked Ly$\alpha$ resembling the local Green-Pea sample: Matthee et al. 2018). The SFRs are all calculated from the UV because the average $\Sigma_{\mathrm{SFR}}$ vs. $z$ relation we calibrate our model against is derived from the UV (Shibuya et al., 2015). A caveat is that the striking abundance of sources populating the top-left corner of Figure 5 were selected to be Green-Pea-like (e.g., Izotov et al., 2016b) for further follow-up, i.e., with very high $\Sigma_{\mathrm{SFR}}$, but it is nonetheless remarkable that the selection is so successful given the long history of LyC non-detections. While these individual LyC sources provide useful clues about the properties favoring LyC escape, they may be extreme outliers given their rarity. However, Marchi et al. (2018) find that even among normal star-forming galaxies at $z\sim 4$, the UV compact sources (which hence have higher $\Sigma_{\mathrm{SFR}}$) are likelier to be leaking LyC.

Independently, recent state-of-the-art hydrodynamical simulations have put forth the scenario of spatially concentrated star-formation, turbulence, and feedback carving out channels in the ISM through which LyC photons can stream out of the galaxy (e.g., Ma et al., 2016; Sharma et al., 2016; Safarzadeh & Scannapieco, 2016; Trebitsch et al., 2017; Katz et al., 2018; Rosdahl et al., 2018; Kakiichi & Gronke, 2019; Kimm et al., 2019). For instance, Ma et al. (2016) describe supernovae clearing ionized channels in the ISM around stellar birth-clouds. However, the massive stars exploding as supernovae are precisely the ones producing most of the ionizing photons. Invoking effects like binarity and rotation allow a significant population of UV luminous stars to survive longer and pump LyC through the newly cleared ISM (Choi et al., 2017). The ionized channels visible in high-resolution Ly$\alpha$ spectra of LyC leakers (Vanzella et al., 2019) support this picture.

4.2 $\Sigma_{\mathrm{SFR}}$ in our Empirical Model

We use the usual definition of \Sigma_{\mathrm{SFR}}\(e.g., Shibuya et al., 2019):

[TABLE]

The SFR in our model is a function of the halo accretion rate and a redshift-invariant efficiency that converts the halo accretion rate into an SFR. To calculate effective radii ($R_{\mathrm{gal}}$) for the galaxies in our model we assume the angular momenta of the galaxies are a fixed fraction of their DM halo (Mo et al., 1998). In particular, we relate $R_{\mathrm{gal}}=\lambda R_{\mathrm{halo}}$ where $\lambda$ is the spin parameter of the halo. We set $\lambda=0.031$ to reproduce the observed $\Sigma_{\mathrm{SFR}}$-$z$ relation from Shibuya et al. (2015) (see Figure 5). Remarkably, we are able to match the exact evolution of \Sigma_{\mathrm{SFR}}\with redshift via this single parameter that only sets the normalization. This is in line with Shibuya et al. (2015)’s finding that the ratio of galaxy size and halo size does not evolve significantly with redshift. We add log-normal scatter while assigning sizes, $\sigma_{\mathrm{log\lambda}}=0.22$ dex consistent with Kravtsov (2013); Somerville et al. (2018); Jiang et al. (2018).

4.3 Fitting for f_{\mathrm{esc}}$$\propto \Sigma_{\mathrm{SFR}}\

We assume a simple power law dependence f_{\mathrm{esc}}$$=a\times\left(\Sigma_{\mathrm{SFR}}/\Sigma_{\mathrm{SFR,max}}\right)^{b}. $\Sigma_{\mathrm{SFR,max}}=1000\ \mathrm{M_{\odot}yr^{-1}kpc^{-2}}$ is close to the value for the maximum \Sigma_{\mathrm{SFR}}\that can be sustained without radiation pressure instabilities (Thompson et al., 2005; Heiderman et al., 2010; Hopkins et al., 2010). The scatter in $\lambda$ produces a maximum \Sigma_{\mathrm{SFR}}\of typically $\sim 220\ \mathrm{M_{\odot}yr^{-1}kpc^{-2}}$ in our model. We fit for the coefficients $a$ and $b$ by summing the ionizing photon contributions of each individual galaxy as detailed in §2.2 and perform Bayesian inference against the reionization constraints from §2.3 using the dynesty nested sampling package (Speagle, 2019).

We add one additional constraint: the f_{\mathrm{esc}}\of the Steidel et al. (2018) sample, f_{\mathrm{esc}}$$=0.09\pm 0.01 for a stack of $\sim 120$ $M_{\mathrm{UV}}<-19.5$ LBGs that we assume follow the average $\Sigma_{\mathrm{SFR}}$-$z$ relation at $z\sim 3$. While it is impossible to robustly constrain f_{\mathrm{esc}}\for individual sources at high-$z$ due to the stochasticity of the intervening IGM, Steidel et al. (2018) stack in a narrow redshift bin across multiple lines of sight and correct for the mean IGM at that redshift. Further, the extremely deep spectra in their sample ($\sim 10$-hour exposures on a 10m telescope) show weak ISM lines that can be used to fine-tune models to match the covering fraction, correct for attenuation, and produce a robust estimate of $f_{\mathrm{esc}}$. The key assumption here is that the relationship between \Sigma_{\mathrm{SFR}}\and f_{\mathrm{esc}}\at $z\sim 3$ holds at higher-$z$ – we argue that since f_{\mathrm{esc}}\largely depends on the covering fraction of neutral gas at $z>3$, and not dust, this is a justifiable assumption (see §6.1). We do not include any of the individual LyC leakers depicted in Figure 5 in our fits because estimates of f_{\mathrm{esc}}\for any individual source are highly uncertain due to the transmission along a single IGM line of sight being unmeasurable.

We find $a=1.6^{+0.3}_{-0.3}$ and $b=0.4^{+0.1}_{-0.1}$ by deploying a uniform prior over 0 to 5 for $a$ and $-5$ to 5 for $b$. We assume a uniform prior such that 0\leq$$f_{\mathrm{esc}}$$\leq 1, so our best-fit relation effectively is:

[TABLE]

The tight posteriors are driven by the Steidel et al. (2018) constraint that directly links f_{\mathrm{esc}}\to $\Sigma_{\mathrm{SFR}}$, while the constraints from §2.3 are useful in deciding the positive sign of the dependence (Figure 6). We emphasize that in fitting for this power-law we allow for negative powers (i.e., an $f_{\mathrm{esc}}$-\Sigma_{\mathrm{SFR}}\anti-correlation) that are rejected by the evidence since they fail to conclude reionization by $z\sim 6$.

Model I fits for a very similar evolution of the ionizing photon budget, $\dot{n}_{\mathrm{ion}}(z)$, compared to our more physically motivated Model II (top panels of Figure 7). Which is to say, the evolution of $\dot{n}_{\mathrm{ion}}(z)$ and average f_{\mathrm{esc}}\of $\sim 20\%$ in both the models is similar during reionization. However, the way the similar $\dot{n}_{\mathrm{ion}}(z)$ is distributed among galaxies differs radically between the two models in that instead of a constant f_{\mathrm{esc}}$$=0.2 across all galaxies, a minority of galaxies that are more massive and UV bright tend to have high \Sigma_{\mathrm{SFR}}\and thus high escape fractions (bottom panels of Figure 7). The proportion of these high \Sigma_{\mathrm{SFR}}\galaxies as well as the mean \Sigma_{\mathrm{SFR}}\grows with redshift driven primarily by their increasing compactness and naturally explains the evolution of f_{\mathrm{esc}}\from $\sim 0\%$ at $z\sim 0$ to $\sim 10\%$ at $z\sim 2.5-4$ to $\sim 20\%$ at $z>6$ (Figure 8).

5 Discussion

5.1 Rapid Reionization at $z=6-8$

Both our models produce virtually identical reionization histories (top right panel of Figure 7). The mid-point of reionization ($z_{50}$) is $6.83^{+0.24}_{-0.20}$ while the duration of reionization ($z_{90}-z_{10}$) is tightly constrained to be $\Delta_{z}=3.76^{+0.05}_{-0.04}$ (bottom right panel, Figure 4). The universe goes from $90\%$ neutral at $z=8.22^{+0.25}_{-0.22}$ to $10\%$ neutral at $z=6.25^{+0.26}_{-0.22}$, in $\sim 300$ Myrs (see Table 2). This pace is faster than estimated by earlier studies (e.g., Robertson et al., 2015; Planck Collaboration et al., 2018) and is driven by the sharp drop in the ionizing emissivity at early times (Figure 2). The high neutral fraction measurements at late times ($\gtrsim 50\%$ at $z\sim 7$) from Ly$\alpha$ damping combined with the dark fraction requirement for the end of reionization by $z\sim 5.9$ favor this rapid pace. A corollary of this timeline is that efforts to understand the sources of reionization do not have to probe the highest redshifts since more than half of the process occurs between $z=6-7$.

5.2 The ‘Oligarchs’ that Reionized the Universe

In both our models, especially in Model II, we find the reionization budget ($\dot{n}_{\mathrm{ion}}$) is concentrated in a ultra-minority of galaxies with high \Sigma_{\mathrm{SFR}}\at M_{\mathrm{UV}}$$<-18, $\log(M_{\star}/M_{\odot})>8$ (see Figure 9 and Table 1). In Model II less than $5\%$ of galaxies constitute $\gtrsim 80\%$ of the reionization budget. Adopting a popular income inequality measure from Macroeconomics, the Gini coefficient (Gini, 1912), we find the \dot{n}_{\mathrm{ion}}\distributions for Model II at $z=$[6, 7, 8] have Gini coefficients of [0.93, 0.92, 0.90]. For reference, a distribution comprised of equal numbers has a Gini coefficient $\sim 0$ and a distribution where all the density is held by a single member has a coefficient of $\sim 0.99$. Drawing again from the language of wealth concentration, we christen this ultra-minority of galaxies that dominate reionization the “oligarchs”.

In Figure 8 we show the occurrence of the oligarchs grows with redshift as \Sigma_{\mathrm{SFR}}\increases and galaxies become more compact and star-forming. Consequently, the mean f_{\mathrm{esc}}\also grows with redshift but note that it never exceeds an average of $\sim 25\%$ even in the highest mass bin ($\log(M_{\star}/M_{\odot})=9-10$). At $z\sim 4$ we predict a mean f_{\mathrm{esc}}\of $\sim 10\%$ with $\sim 10\%$ of sources displaying f_{\mathrm{esc}}$$>20\% for galaxies at $\log(M_{\star}/M_{\odot})=9-10$. This is a faithful representation of the current observational situation at $z\sim 2-4$ where stacks produce average f_{\mathrm{esc}}\of $\sim 10\%$ (Marchi et al., 2017; Steidel et al., 2018; Oesch et al. in prep., 2019) while a small fraction of sources show f_{\mathrm{esc}}\$$>20\%, even reaching f_{\mathrm{esc}}$$>50\% (e.g., Naidu et al., 2017; de Barros et al., 2016; Vanzella et al., 2018) for galaxies in a similar mass range. In fact, in the Steidel et al. (2018) sample at $z\sim 3$ (which we approximately compare with our predictions at $z\sim 4$) 10/124 sources ($\sim 8\%$) show significant LyC leakage while the stacked mean is $\sim 10\%$ – the agreement in the fraction of sources with high f_{\mathrm{esc}}\is noteworthy since we fit our model against the mean f_{\mathrm{esc}}\and the fraction of f_{\mathrm{esc}}$$>20\% galaxies is a genuine prediction (top-right panel of Figure 8).

In Model I this distribution of the ionizing budget is a direct reflection of the shape of the UVLF arising from our model since $\xi_{\mathrm{ion}}$ does not vary with M_{\mathrm{UV}}\(Figure 2) and all galaxies have the same $f_{\mathrm{esc}}$. Steeper \alpha_{\mathrm{UV}}\keeping all else same, will lead to a lower average $f_{\mathrm{esc}}$, larger luminosity densities at early times, a less oligarchic distribution, and possibly tension with current constraints favoring late and rapid reionization (see §5.3). However, in Model II, $M_{\mathrm{UV}}>-17$ galaxies are limited to very low f_{\mathrm{esc}}\and they constitute a negligible portion of the reionization budget. Truncating as high as $M_{\mathrm{UV}}=-17$ has no effect on the model parameters shown in Figure 6, i.e., the ionizing emissivity between $M_{\mathrm{UV}}=-13.5$ to $-17$ is severely down-weighted by assigning a low $f_{\mathrm{esc}}$. Even with steeper $\alpha_{\mathrm{UV}}$, we expect the oligarch scenario to hold since Model II has the flexibility to ensure late reionization as required by the constraints by setting f_{\mathrm{esc}}\$$\sim 0 for the numerous faint galaxies with low $\Sigma_{\mathrm{SFR}}$. We discuss \alpha_{\mathrm{UV}}\further in §5.3.

5.3 “Democratic” Reionization by Faint Galaxies and the Faint-End Slope of the UVLF in a Rapidly Reionizing Universe

Reionization by oligarchs stands in sharp contrast to “democratic” reionization that is dominated by copious faint sources that lie at M_{\mathrm{UV}}$$>-18 and might potentially have high escape fractions (e.g., Oesch et al., 2009; Bouwens et al., 2011; Wise et al., 2014; Atek et al., 2015; Anderson et al., 2017; Livermore et al., 2017; Finkelstein et al., 2019). Faint galaxies emerged as the candidate-leaders of reionization because the steep slopes (\alpha_{\mathrm{UV}}\$$\leq-2 at $z>6$) of the UVLF measured after the installation of $HST$/WFC3 implied they dominated the luminosity density (e.g., Bouwens et al., 2012). The $\tau$ measurements from WMAP-9 ($0.089\pm 0.014$, $z_{50}=10.5\pm 1.1$) and Planck Collaboration et al. (2016) ($0.066\pm 0.013$, $z_{50}=8.8\pm 1.3$) required significant reionization at $z>8$ and hence large contributions towards the ionizing emissivity from faint galaxies (e.g., Robertson et al., 2013, 2015; Bouwens et al., 2015a). Concurrently, the very low f_{\mathrm{esc}}\reported for bright star-forming galaxies out to $z\sim 4$ (see §1) and the sharply dropping AGN luminosity function (Kulkarni et al., 2019) further shifted the spotlight onto faint star-forming galaxies.

However, the recent constraints on neutral fractions detailed in §2.3 and the latest Planck $\tau$ favor late, rapid reionization between $z=6-8$ (we calculate $z_{50}=6.83^{+0.24}_{-0.20}$ for Model I) i.e., high emissivity from faint galaxies at $z>8$ is no longer required. This, and the high average f_{\mathrm{esc}}\measured even for more massive, M_{\mathrm{UV}}$$<-18 galaxies allow for reionization by the oligarchs. At $z>8$, \dot{n}_{\mathrm{ion}}\must be low enough for the universe to remain significantly neutral ($\gtrsim 90\%$), and between $z=8-6$ it must rise sharply to complete reionization. Since $\xi_{\mathrm{ion}}$ evolves modestly with redshift and across M_{\mathrm{UV}}\(see Figure 2), \dot{n}_{\mathrm{ion}}\effectively depends on $\rho_{\mathrm{SFR}}$ ($\alpha_{\mathrm{UV}}$, M_{\mathrm{UV}}\truncation) and $f_{\mathrm{esc}}$.

Latest studies report $\alpha_{\mathrm{UV}}\lesssim-2$ at $z\geq 6$, albeit with significant uncertainties, that grows steeper with redshift at a rate $d\alpha/dz\sim-0.1$ (e.g., Finkelstein et al., 2015; Livermore et al., 2017; Bouwens et al., 2017; Atek et al., 2018; Ishigaki et al., 2018; Oesch et al., 2018). We compare our reionization histories with models that assume these steep slopes and model $\rho_{\mathrm{SFR}}$ based on Schechter parameters extrapolated from $z<10$ fits in Figure 10. Assuming \alpha_{\mathrm{UV}}$$<-2 and setting f_{\mathrm{esc}}\preferentially higher in the faintest galaxies requires integration down to M_{\mathrm{UV}}$$=-10 and reionizes large volumes of the $z>8$ universe, reaching $\bar{x}_{\mathrm{HI}}\sim 40\%$ at $z=8$ (Finkelstein et al., 2019), in tension with the damping wing measurements ($\bar{x}_{\mathrm{HI}}\sim 90\%$). Using a constant f_{\mathrm{esc}}\across all galaxies with \alpha_{\mathrm{UV}}$$\leq-2 like in Ishigaki et al. (2018) and Robertson et al. (2015) requires integrating to -$$M_{\mathrm{UV}}$$=11-13 and still makes for too low of a neutral fraction ($\sim 60-70\%$) at $z=8$. Simply lowering the constant f_{\mathrm{esc}}\in these models would delay reionization but then it would not conclude by $z\sim 6$ – raising the f_{\mathrm{esc}}\while lowering the M_{\mathrm{UV}}\(trunc.) and/or shallower \alpha_{\mathrm{UV}}\are needed. We illustrate this in the right panel of Figure 10, where we assume Schechter parameters from Finkelstein et al. (2019), $\xi_{\mathrm{ion}}$ from this work, and a constant f_{\mathrm{esc}}\to evaluate how likely various combinations of f_{\mathrm{esc}}\and M_{\mathrm{UV}}\truncation are (as per constraints from §2.3). A truncation M_{\mathrm{UV}}\of $\leq-15$, implying a limited role for fainter galaxies is favored by the constraints. The general feature of \dot{n}_{\mathrm{ion}}\dominated by faint galaxies in the models discussed above is that \dot{n}_{\mathrm{ion}}\is already high at $z=10$, resulting in smooth and early reionization (Figures 3 and 10).

On the other hand, the required late and rapid reionization is naturally produced by shallower ($\alpha_{\mathrm{UV}}\geq-2$) faint-end slopes (Model I, Model II), distributions of f_{\mathrm{esc}}\highly skewed toward brighter galaxies (Model II), and/or a sharp drop in the $z>8$ $\rho_{\mathrm{SFR}}$ in models linking star-formation to dark matter accretion (e.g. Model I, Model II, Mason et al., 2015; Mashian et al., 2016). Truncating at $M_{\mathrm{UV}}=-16\ (-17)$ in Model I (Model II) produces no change in the reionization histories and model parameters, i.e., the ionizing emissivity requires no contributions from $M_{\mathrm{UV}}>-16\ (-17)$ galaxies (Figures 9 and 11). $M_{\mathrm{UV}}>-16$ galaxies are rare in the early universe and their appearance causes \dot{n}_{\mathrm{ion}}\to rise steeply by more than a dex between $z=6-10$ (Figure 3). Thus, while the faint-end slope of the UVLF may indeed be extremely steep, galaxies at $M_{\mathrm{UV}}>-16$ must play only a minimal role in order to achieve rapid and late reionization. Model II explains this as very low f_{\mathrm{esc}}\occurring in these abundant albeit low \Sigma_{\mathrm{SFR}}\galaxies.

5.4 Observing the Oligarchs in Action: Promising Hints and Future Prospects

The luminous Ly$\alpha$ emitter, COLA1 at $z=6.6$ (f_{\mathrm{esc}}$$\sim 30\%, M_{\mathrm{UV}}$$=-21.5, $\Sigma_{\rm SFR}=100\ \mathrm{M}_{\odot}/\mathrm{yr}/\mathrm{kpc}^{2}$) is a poster-child oligarch (Hu et al., 2016; Matthee et al., 2018). It displays double-peaked Ly$\alpha$ with a low-peak separation reminiscent of local LyC leakers (Verhamme et al., 2017). More statistically, Songaila et al. (2018); Hu et al. (2016); Matthee et al. (2018) find luminous Ly$\alpha$ emitters at $z\sim 7$ have line profiles that are broader (while not being Active Galactic Nuclei (AGN)) and more complex than their lower luminosity counterparts with two sources in a sample of seven showing blue wings despite a highly neutral IGM (e.g., Mason et al., 2018). We speculate these galaxies are oligarchs with high escape fractions that are able to carve out their own ionized bubbles perhaps allowing for their whole line profiles (including blue wings and peaks) to escape unattenuated by neutral gas. The lower luminosity, low f_{\mathrm{esc}}\sources have narrow, less complex Ly$\alpha$ profiles that are perhaps truncated by the neutral gas surrounding them. High-resolution ($R>4500$) Ly$\alpha$ surveys with well-defined selection and completeness functions at $z\sim 0-6$ will help test if these complex Ly$\alpha$ profiles that have been linked to ionized channels and thus LyC f_{\mathrm{esc}}\(Vanzella et al., 2019; Rivera-Thorsen et al., 2019; Herenz et al., 2017) grow more common with redshift and with galaxy properties like $\Sigma_{\mathrm{SFR}}$. Since we do not expect f_{\mathrm{esc}}\to evolve appreciably between $z\sim 6-8$ as \Sigma_{\mathrm{SFR}}\flattens (Figure 5), such a survey can be limited to $z<6$ where the IGM transmission is higher and Ly$\alpha$ is easily observable.

Another intriguing observation is that of an overdensity of 17 HST dropouts at $z\sim 7$. In an extremely long integration (22.5 hrs on VLT/VIMOS), Castellano et al. (2018) find Ly$\alpha$ emission arising only from three UV-bright galaxies among the dropouts while all their faint galaxies are undetected in Ly$\alpha$ despite Ly$\alpha$ EWs generally anti-correlating with brightness. We speculate the bright oligarchs with high f_{\mathrm{esc}}\have reionized their immediate surroundings rendering them transparent to Ly$\alpha$ while the fainter sources lie just outside these ionized bubbles. With JWST’s planned censuses at high-$z$ more such ionized overdensities at $z>6$ will come into view and deep follow-up spectroscopy that reveals features of LyC f_{\mathrm{esc}}\(e.g., multi-peaked Ly$\alpha$) will help test if they are indeed powered by oligarchs.

Our proposed scenario also has strong implications for the topology of reionization, with the distribution of ionized bubble sizes and the patchiness resulting from our model likely lying somewhere intermediate between AGN-driven reionization (e.g., Kulkarni et al., 2017) and reionization by widely distributed, numerous faint sources (e.g., Geil et al., 2016). Upcoming 21cm surveys will provide a strong test of this prediction (e.g., Hutter et al., 2019a, b; Seiler et al., 2019). Our empirical model also tracks the spatial distribution of galaxies and this information can be coupled with models for f_{\mathrm{esc}}\to produce more quantitative predictions. We defer this to future work.

5.5 Related Work: the $f_{\mathrm{esc}}$-\Sigma_{\mathrm{SFR}}\Connection

Heckman et al. (2001) explicitly link f_{\mathrm{esc}}\to a critical value of \Sigma_{\mathrm{SFR}}\$$\sim\mathrm{0.1M_{\odot}/yr/kpc^{2}} above which they observe the occurrence of strong winds becomes common in star-burst galaxies. They hypothesize that these winds are responsible for LyC leakage. Sharma et al. (2016) adopt this idea in the EAGLE simulations (Schaye et al., 2015; Crain et al., 2015) setting f_{\mathrm{esc}}$$=0.2 for star-forming regions above the critical \Sigma_{\mathrm{SFR}}\and averaging these in each galaxy to find an emissivity consistent with Bouwens et al. (2015a). The f_{\mathrm{esc}}$$=0.2 hard upper-limit was motivated by the lack of LyC detections, prior to the recent LyC renaissance detailed in §1.

We do not assume a threshold or an upper bound emboldened by recent discoveries (see Figure 5) and empirically constrain an $f_{\mathrm{esc}}$-\Sigma_{\mathrm{SFR}}\dependence. We use a prior that allows for no relation ($b=0$) or an anti-correlation ($b<0$) and fit against the latest reionization constraints. Note that the Sharma et al. (2016) prescription, even though it invokes \Sigma_{\mathrm{SFR}}\, ends up effectively similar to our Model I that fits $f_{\mathrm{esc}}$=$0.21^{+0.06}_{-0.04}$, since essentially all galaxies at $z\sim 6-8$ at $M_{\mathrm{UV}}<-13.5$ have \Sigma_{\mathrm{SFR}}$$>\mathrm{0.1M_{\odot}yr^{-1}kpc^{-2}} (see bottom left panel in Figure 7) and have sizes on the order of $1$ kpc (their beam size). Furthermore, Sharma et al. (2016) find EAGLE galaxies at $M_{\mathrm{UV}}<-18$ produce $\sim 50\%$ of the reionization budget. Our budget is far more oligarchic ($>80$%), and our reionization far more rapid ($z\sim 6-8$; see §5.1).

6 Caveats and Open Questions

6.1 On Dust and f_{\mathrm{esc}}\at $z>6$

In our Model I, f_{\mathrm{esc}}\implicitly folds in the role of dust and all other processes that may curtail LyC leakage. This is not the case in Model II in which f_{\mathrm{esc}}\is solely a function of $\Sigma_{\mathrm{SFR}}$. Further, we adapt the Steidel et al. (2018) measurement to constrain Model II assuming that the relationship between \Sigma_{\mathrm{SFR}}\and f_{\mathrm{esc}}\at $z\sim 3$ carries over to higher-$z$, unmodulated by dust. We justify these assumptions here.

At $z\sim 3$ already, deep stacks of typical LBGs show that it is not dust, but photoelectric absorption in the ISM that dominates the attenuation of LyC photons to the extent that

[TABLE]

where $f_{\mathrm{cov}}$ is the H I covering fraction (Reddy et al., 2016a, b; Steidel et al., 2018). Moving into the EoR, this approximation is likely even better since the dust attenuation at $z>6$ appears to be lower (Bouwens et al., 2016b; Fudamoto et al., 2017) though significant uncertainties persist (e.g., Casey et al., 2018). Our attenuation prescription bears this out – for instance, on average, an $M_{\mathrm{UV}}=$[$-22$,$-20$,$-18$] galaxy at $z=4$ has $A_{\mathrm{UV}}$=[1.51, 1.05, 0.6] and a $z=8$ galaxy has $A_{\mathrm{UV}}$=[1.40, 0.6, 0.0]. The key physical picture motivating Model II, that of spatially concentrated star-formation, winds, and feedback carving out ionized gas channels is intimately linked with $f_{\mathrm{cov}}$ and thus needs no extra dust parameter, since $f_{\mathrm{cov}}$ essentially determines $f_{\mathrm{esc}}$. However, we note that not explicitly modeling dust in Model II prevents a simple extrapolation of f_{\mathrm{esc}}\using our fit power-law to $z\sim 0$ where attenuation of LyC by dust is highly significant though the qualitative picture and general trend stands.

6.2 Effect of $M_{UV}$ Truncation

We have limited all our calculations to galaxies with $\mathrm{SFR}>0.02\ \mathrm{M_{\odot}}\mathrm{yr^{-1}}$ which corresponds approximately to $M_{\mathrm{UV}}$(observed)$<-13.5$. This limitation arises from the resolution of the dark matter simulations our model is built on as well as the significant uncertainties around the UVLF fainter than this magnitude (Livermore et al., 2017; Bouwens et al., 2017; Atek et al., 2018). What effect does this truncation have on our results?

In Model I extending to fainter magnitudes adds to the ionizing emissivity and should lower the average f_{\mathrm{esc}}\we report. However, since our model has \alpha_{\mathrm{UV}}$$>-2 during reionization we expect this lowering to become negligible at M_{\mathrm{UV}}$$>-13.5, since the differential change to the bulk \dot{n}_{\mathrm{ion}}\asymptotes to zero. We explore this by shifting the limiting magnitude brighter (see Figure 11). We find our f_{\mathrm{esc}}\solution is essentially converged at M_{\mathrm{UV}}$$<-16 since moving down to M_{\mathrm{UV}}$$<-13.5 produces no appreciable change. Moving further down to M_{\mathrm{UV}}$$<-11 should make an even smaller difference especially if the UVLF turns over around these magnitudes due to inefficient star-formation and photo-evaporation in low-mass halos (e.g., Gnedin, 2016).

In Model II the majority of the low-luminosity galaxies have extremely low f_{\mathrm{esc}}\(see bottom panels of Figure 8) to go along with their lower LyC output so the exclusion of M_{\mathrm{UV}}$$>-13.5 galaxies or the faint-end slope have negligible impact on model parameters. We have verified that even truncating as high as M_{\mathrm{UV}}$$<-17 produces very similar reionization histories to those reported here.

6.3 Model Dependent Constraints

The Ly$\alpha$ damping measurements for $z>7$ quasars and galaxies prove to be the most constraining for our Model I. However, these are model dependent constraints in that their reported $\bar{x}_{\mathrm{HI}}$ is a product of several assumptions, e.g., about how Ly$\alpha$ is processed by the ISM at high-$z$ or about the intrinsic spectrum of reionization epoch quasars. These assumptions while reasonable are yet to be tested (e.g., see §2 in Mason et al., 2018). In Model I the model-independent $\tau$ and dark fraction by themselves are unable to zero in on an f_{\mathrm{esc}}\solution, but they rule out f_{\mathrm{esc}}$$\lesssim 15\% and so favor rapid reionization histories. In Model II the Ly$\alpha$ damping measurements are unable to collapse the posterior much further beyond $\tau$ and the dark fraction combined, and the Steidel et al. (2018) measurement proves crucial (see Figure 6). This measurement depends on several assumptions e.g., stellar population model predictions at $<912\mathrm{\AA}$, the IGM+CGM transmission functions, and details of the “hole” and “screen” models for f_{\mathrm{esc}}\developed in their work. These are all sources of systematic uncertainty on the reported absolute $f_{\mathrm{esc}}$. We have tested that the sign of the power ($>0$) recovered for the $f_{\mathrm{esc}}$-\Sigma_{\mathrm{SFR}}\dependence is not sensitive to the exact scale of their reported f_{\mathrm{esc}}\as long as f_{\mathrm{esc}}$$>0, and the dark fraction measurement that ensures the timely conclusion of reionization is used.

6.4 Ionizing Emissivity at $z<6$ and the Role of AGN

In this work we have focused on reionization driven by galaxies and fit for parameters that show them satisfying all constraints outlined in §2.3 without invoking AGN. This is supported by latest determinations of AGN luminosity functions that limit their contribution to $<5\%$ at $z=6$ (e.g., Kulkarni et al., 2019, though see Giallongo et al. 2015; Boutsia et al. 2018). In a companion work we deploy a similar framework, but assume nothing about the underlying ionizing population and fit a non-parametric model to recover \dot{n}_{\mathrm{ion}}$$(z) (Mason et al., 2019b). A sharp drop in \dot{n}_{\mathrm{ion}}$$(z) at $z=6-8$ fully consistent with the galaxy-only models presented in this work is recovered.

A related issue is whether our models overrun the constraints on ionizing emissivity at lower-$z$ plotted in the top-left panel of our Figure 7 when combined with the AGN emissivity at lower redshifts (Becker & Bolton, 2013; Kuhlen & Faucher-Giguère, 2012). Note that in Model I we make no claims about f_{\mathrm{esc}}\at $z<6$ and fit a $\sim 20\%$ f_{\mathrm{esc}}\during the short window when reionization transpires. In Model II we have an evolving f_{\mathrm{esc}}\that falls from $\sim 20\%$ at $z>6$ to $\sim 10\%$ ($\sim 4\%$) at $z=3$ ($z=2$). This causes \dot{n}_{\mathrm{ion}}\to flatten and turn-over at lower-$z$ so that AGN can dominate the emissivity – we begin to see this in the top-left panel of Figure 7 (dotted purple curve).

6.5 Cosmic Variance and Completeness

Due to the finite volume of the N-body simulation on which our empirical model is built ($10^{6}$ Mpc3), we miss some of the most massive halos (see Figure 20 and Appendix B in Tacchella et al. 2018). At $z>6$ these halos also tend to be the most star-forming and UV bright. Comparing to an analytical halo mass function (Sheth et al., 2001) and applying a completeness correction produces a $\lesssim 0.3$ dex difference at the brightest end of the UVLF at $z=10$. The correction is smaller at lower redshifts, where the bulk of reionization occurs and hence the magnitude of the effect is likely small. For Model I, the mean f_{\mathrm{esc}}\estimated would slightly shrink due to the missing luminosity density. For Model II, extrapolating from the trends shown in Figure 8 for the proportion of oligarchs as a function of galaxy mass at $z>6$, including these massive halos would make the reionization budget even more oligarchic. An update to the empirical model using a larger box that also self-consistently includes AGN is currently under preparation. The larger box will also allow us to address cosmic variance by resampling smaller volumes and computing the resulting scatter introduced in $f_{\mathrm{esc}}$.

7 Summary

Using an empirical model that accurately predicts observations (e.g., the sharp drop in $\rho_{\rm SFR}$ at $z>8$, UVLFs at $z>4$, $\xi_{\mathrm{ion}}$ at $z\sim 4-5$, the $z$-evolution of $\Sigma_{\mathrm{SFR}}$) and leveraging recent measurements of the timeline of reionization (e.g., $\bar{x}_{\mathrm{HI}}$ at $z\gtrsim 7$, the Planck $\tau$) we constrain $f_{\mathrm{esc}}$, the most uncertain parameter in reionization calculations. Deploying two models – one assuming a constant f_{\mathrm{esc}}\across all galaxies during reionization (Model I), another linking f_{\mathrm{esc}}\to \Sigma_{\mathrm{SFR}}\(Model II) – we find the following:

• •

In both our models, M_{\mathrm{UV}}$$<$$-13.5 star-forming galaxies need an average f_{\mathrm{esc}}$$\sim20$\%$ at $z>6$ to conclude reionization, a factor of only 2 higher than the recently measured f_{\mathrm{esc}}$$\sim 10\% at $z\sim 2.5-4$. [Figures 4, 7]

• •

Our Model II explains this evolution in f_{\mathrm{esc}}\by appealing to \Sigma_{\mathrm{SFR}}\that decreases by $\sim$2.5 dex between $z=8$ (f_{\mathrm{esc}}$$\sim 0.2) and $z=0$ (f_{\mathrm{esc}}$$\sim 0) by fitting f_{\mathrm{esc}}$$\propto($$\Sigma_{\mathrm{SFR}}$$)^{0.4}. This $f_{\mathrm{esc}}$-\Sigma_{\mathrm{SFR}}\connection is inspired by the newly emerging sample of LyC leakers that show higher \Sigma_{\mathrm{SFR}}\than the average galaxy population at their redshifts. Latest hydrodynamical simulations qualitatively support the idea of spatially concentrated star-formation blowing channels in the ISM through which LyC produced by long-lived, rotating, binary stars escapes. [Figures 5, 6]

• •

The universe goes from $90\%$ neutral to $10\%$ neutral in a short span of $\sim 300$ Myrs between $z\sim 6-8$, and favored by Ly$\alpha$ damping measurements requiring a $\gtrsim 50\%$ neutral universe at $z\sim 7$. This conclusion stands even considering only model-independent constraints ($\tau$, dark fraction) that rule out f_{\mathrm{esc}}$$\lesssim 15\%. [Figure 4, Table 2]

• •

The bulk of the reionization budget ($\sim$50$\%$ in Model I, $\sim$80$\%$ in Model II) is concentrated among a small number ($<$5$\%$) of galaxies (the “oligarchs”). This is due to the faint-end slopes of the UVLF (\alpha_{\mathrm{UV}}$$>-2) in our model and the distribution of f_{\mathrm{esc}}\skewed toward high $\Sigma_{\mathrm{SFR}}$, massive galaxies. The oligarchs are compact ($R_{\rm gal}\sim 0.5$ kpc), have higher \Sigma_{\mathrm{SFR}}\than average ($\sim$10$-$20 $\mathrm{M_{\odot}yr^{-1}kpc^{-2}}$), are relatively massive (\log(M_{\star}/M_{\odot})$$>8) and are UV bright (M_{\mathrm{UV}}$$<$$-18). The fraction of these oligarchs grows with redshift, while keeping the average f_{\mathrm{esc}}\to $\sim$20$\%$. Extrapolating to $z\sim 3-4$ we match the current situation where a small fraction of galaxies ($\lesssim$10$\%$) display f_{\mathrm{esc}}$$>0.2, some even exceeding $50\%$, while the average f_{\mathrm{esc}}\stays at $\sim$10$\%$. [Figures 8, 9, Table 1]

• •

Faint galaxies are disfavored to drive reionization. When faint galaxies with steep $\alpha_{\mathrm{UV}}<-2$ dominate the emissivity, they ionize large volumes of the universe at $z=7-8$, in tension with Ly$\alpha$ damping constraints that require a $60-90\%$ neutral universe at these redshifts. Shallower faint-end slopes ($\alpha_{\mathrm{UV}}>-2$) and/or f_{\mathrm{esc}}\distributions skewed toward massive galaxies like in our models ensure high neutral fractions at late times while also completing reionization by $z=6$. Concurrently, the motivation for excluding galaxies at $M_{\mathrm{UV}}<-18$ with high f_{\mathrm{esc}}\as the protagonists of reionization has grown weaker as the observational picture has shifted to these galaxies being able to produce f_{\mathrm{esc}}\$$\sim 10\% at $z=2.5-4$. [Figure 5, 10]

Our predictions are eminently testable since the oligarchs are bright, currently observable galaxies. Deep Ly$\alpha$ surveys at high-resolution ($R>4500$) spanning $z\sim 0-6$ should show a growing incidence of galaxies with multi-peaked Ly$\alpha$ at higher-$z$. These peaks represent ionized channels for LyC escape, as seen in the $z\lesssim 4$ LyC leakers. Upcoming 21 cm experiments should infer a bubble size distribution with a high Gini coefficient as the first ionization fronts form predominantly around the oligarchs.

We are grateful to the referee, Brant Robertson, for a thorough report that greatly improved the clarity of this work. It is a pleasure to acknowledge illuminating discussions with Steve Finkelstein, Jorryt Matthee, Daniel Schaerer, David Sobral, Irene Shivaei, Ben Johnson, and Chuck Steidel that enriched this work. We thank Daniel Eisenstein, Avi Loeb, and Lars Hernquist for participating in a robust discussion of these results during RPN’s research exam. We thank Steve Finkelstein for sharing the posteriors from Finkelstein et al. (2019) featured in Figure 10. We are grateful to Eros Vanzella for sharing their data on Ion3 and for useful comments on Naidu et al. (2018) that sparked parts of this work. RPN thanks Michelle L. Peters for a hilarious proofreading session that highlighted the ludicrous verbiage we deploy in the astrophysics community. RPN thanks the organizers and attendees of IAU Symposium 352 “Uncovering early galaxy evolution in the ALMA and JWST era” held in Viana do Castelo for providing much-needed affirmation that he is not alone in the often testing endeavour that is working on $f_{\mathrm{esc}}$. RPN gratefully acknowledges an Ashford Fellowship and Peirce Fellowship granted by Harvard University. ST is supported by the Smithsonian Astrophysical Observatory through the CfA Fellowship. CAM acknowledges support by NASA Headquarters through the NASA Hubble Fellowship grant HST-HF2-51413.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. SB acknowledges Harvard University’s ITC Fellowship. PAO acknowledges support from the Swiss National Science Foundation through the SNSF Professorship grant 157567 “Galaxy Build-up at Cosmic Dawn”. The Cosmic Dawn Center is funded by the Danish National Research Foundation. CC acknowledges funding from the Packard foundation.

Appendix A UV Luminosity Functions

In addition to Figure 1, where we compare against literature data at $z=6$ (the redshift at which the faint-end of the UVLF has been measured), we provide UVLFs predicted by our empirical model for $z=7-12$. We also quantify agreement between our model and observational data in Table 3. In general, our predicted UVLFs show good agreement with the data, both at the bright and faint ends. The only tension is at $z=6$ with Atek et al. (2018) at $M_{\mathrm{UV}}<-16$ where they estimate several more galaxies than predicted by our model – we find better agreement with Bouwens et al. (2015b) and Livermore et al. (2017) in the same luminosity regime, while also noting that the excess of brighter galaxies would further support this paper’s point of view about the reionization budget.