NEW CONSTRAINTS ON COSMIC REIONIZATION FROM THE 2012 HUBBLE ULTRA DEEP FIELD CAMPAIGN

The UV luminosity density is supplied by short-lived, massive stars and therefore reflects the time rate of change of the stellar mass density ρ_⋆(z). In the context of our model, there are two routes for estimating the stellar mass density. First, we can integrate the stellar mass density supplied by the star formation rate inferred from the evolving UV luminosity density as

$$\begin{equation} \rho _{\star }(z) = (1-R) \int _{\infty }^{z} \eta _{\mathrm{sfr}}(z^{\prime }) \rho _{\mathrm{UV}}(z^{\prime }) \frac{dt}{dz^{\prime }} dz^{\prime }, \end{equation} \tag{ 7 }$$

where η_sfr(z) provides the stellar population model-dependent conversion between UV luminosity and star formation rate (in units M_☉ yr⁻¹ erg⁻¹ s Hz), R is the fraction of mass returned from a stellar population to the ISM (28% for a Salpeter 1955 model after ∼10 Gyr for a 0.1–100 M_☉ IMF), and dt/dz gives the rate of change of universal time per unit redshift (in units of yr). While η_sfr(z) is in principle time-dependent, there is no firm evidence yet of its evolution and we adopt a constant value throughout.

While the evolving stellar mass density can be calculated in our model, the observational constraints on ρ_⋆ have involved integrating a composite stellar mass function constructed from the UV luminosity function and a stellar mass to UV luminosity relation (González et al. 2011; Stark et al. 2013, see also Labbe et al. 2012). Using near-IR observations with the Spitzer Space Telescope, stellar masses of UV-selected galaxies are measured as a function of luminosity. Stark et al. (2013) find that the stellar mass–L_UV relation, corrected for nebular emission contamination, is well described by

$$\begin{equation} \log m_{\star }= 1.433 \log L_{\mathrm{UV}}-31.99 + f_{\star }(z), \end{equation} \tag{ 8 }$$

where m_⋆ is the stellar mass in M_☉, L_UV is the UV luminosity density in erg s⁻¹ Hz⁻¹, and f_⋆(z) = [0, −0.03, −0.18, −0.40] for redshift z ≈ [4, 5, 6, 7].

The stellar mass density will then involve an integral over the product of the UV luminosity function and the stellar mass m_⋆(L_UV). However, the significant scatter in the m_⋆(L_UV) relation (σ ≈ 0.5 in log m_⋆; see González et al. 2011; Stark et al. 2013) must be taken into account. The Gaussian scatter p[M' − M(log m_⋆)] in luminosity contributing at a given log m_⋆ can be incorporated into the stellar mass function dn_⋆/dlog m_⋆ with a convolution over the UV luminosity function. We write the stellar mass function as

$$\begin{equation} \frac{dn_{\star }}{d\log m_{\star }} = \frac{d M}{d\log m_{\star }} \int _{-\infty }^{\infty } \Phi (M^{\prime }) p[M^{\prime }-M(\log m_{\star })] d M^{\prime }. \end{equation} \tag{ 9 }$$

For vanishing scatter, Equation (9) would give simply dn_⋆/dm_⋆ = Φ[M(m_⋆)] × dM/dm_⋆. The stellar mass density ρ_⋆ can be computed by integrating this mass function as

$$\begin{equation} \rho _{\star }(<m_{\star }, z) = \int _{-\infty }^{\log m_{\star }} \frac{dn_{\star }}{d\log m_{\star }^{\prime }} m_{\star }^{\prime } d\log m_{\star }^{\prime }. \end{equation} \tag{ 10 }$$

A primary feature of the stellar mass function is that the stellar mass–UV luminosity relation and scatter flattens it relative to the UV luminosity function. Correspondingly, the stellar mass density converges faster with decreasing stellar mass or luminosity than does the UV luminosity density (Equation (5)). The stellar mass density then serves as an additional, integral constraint on $\dot{n}_{\mathrm{ion}}$.

Once the evolving production rate of ionizing photons $\dot{n}_{\mathrm{ion}}$ is determined, the reionization history Q_H ii(z) of the universe can be calculated by integrating Equation (1). An important integral constraint on the reionization history is the electron scattering optical depth τ inferred from observations of the CMB. The optical depth can be calculated from the reionization history as a function of redshift z as

$$\begin{equation} \tau (z) = \int _{0}^{z} c \langle n_{\mathrm{H}} \rangle \sigma _{\mathrm{T}}f_{\mathrm{e}} Q_{\mathrm{H\,{\scriptsize{II}}}}(z^{\prime }) H(z^{\prime }) (1+z^{\prime })^{2} dz^{\prime }, \end{equation} \tag{ 11 }$$

where c is the speed of light, σ_T is the Thomson cross section, and H(z) is the redshift-dependent Hubble parameter. The number f_e of free electrons per hydrogen nucleus in the ionized IGM depends on the ionization state of helium. Following Kuhlen & Faucher-Giguère (2012) and other earlier works, we assume that helium is doubly ionized (f_e = 1 + Y_p/2X_p) at z ⩽ 4 and singly ionized (f_e = 1 + Y_p/4X_p) at higher redshifts. To utilize the observational constraints on τ as a constraint on the reionization history, we employ the posterior probability distribution p(τ) determined from the Monte Carlo Markov Chains (MCMC) used in the nine-year WMAP results⁹ as a marginalized likelihood for our derived τ values. This method is described in more detail in Section 5.

We wish to use the evolving UV luminosity density ρ_UV to constrain the availability of ionizing photons $\dot{n}_{\mathrm{ion}}$ as given in Equation (4). Given how significantly ρ_UV increases with the limiting absolute magnitude M_UV, we must evaluate the full model presented in Section 2 for chosen values of M_UV. We select M_UV = −17, M_UV = −13, and M_UV = −10 to provide a broad range of galaxy luminosities extending to extremely faint magnitudes. The observations probe to M_UV = −17, and this magnitude therefore serves as a natural limit. We consider limits as faint as M_UV = −10 as this magnitude may correspond to the minimum mass dark matter halo able to accrete gas from the photoheated intergalactic medium or the minimum mass dark matter halo able to retain its gas supply in the presence of supernova feedback. These scenarios have so far proven impossible to distinguish empirically (Muñoz & Loeb 2011). To make use of the constraints of the UV luminosity density on $\dot{n}_{\mathrm{ion}}$, we need to adopt a likelihood for use with Bayesian inference on parameterized forms of ρ_UV(z). Given the parameterized form of the Schechter (1976) function, we find that the marginalized posterior distributions of ρ_UV are skewed at z ≳ 6, with tails extending to larger log ρ_UV values than a Gaussian approximation would provide. To capture this skewness, when constraining the cosmic reionization history in Section 5 below, we therefore use the full marginalized posterior distribution of ρ_UV provided by the integrated LF determinations calculated in Section 4.

Figure 3 shows the marginalized posterior distribution of the UV luminosity density ρ_UV for limiting magnitudes of M_UV = −17 (red lines), M_UV = −13 (orange lines), and M_UV = −10 (blue lines) for our Schechter function fits to the z ∼ 5–6 LFs of Bouwens et al. (2007), the z ∼ 7–8 LFs of Schenker et al. (2012a), and the z ∼ 7–8 LFs of McLure et al. (2012). Although Schenker et al. (2012a) and McLure et al. (2012) use the same data sets, their luminosity function determinations are based on different selection techniques. They therefore represent independent determinations of the high-redshift luminosity functions and are treated accordingly. We additionally use constraints from z ∼ 4 (Bouwens et al. 2007) and z ∼ 9 (McLure et al. 2012). The posterior distributions for ρ_UV(z ∼ 9) we calculate are likely underestimated in their width owing to an assumption of fixed M_⋆ and α values, but this assumption does not strongly influence our results. In what follows, these posterior distributions on ρ_UV are used as likelihood functions when fitting a parameterized model to the evolving UV luminosity density.

Zoom In Zoom Out Reset image size

To infer constraints on the reionization process, we must adopt a flexible parameterized model for the evolving UV luminosity density. The model must account for the decline of ρ_UV apparent in Figures 2 and 3, without artificially extending the trend in ρ_UV at z ≲ 6 to redshifts z ≳ 8. The rapid decline in ρ_UV between z ∼ 4 and z ∼ 5 suggests a trend of dlog ρ_UV/dlog z ∼ −3, but the higher redshift values of ρ_UV flatten away from this trend, especially for faint limiting magnitudes. Beyond z ∼ 10, we expect that some low-level UV luminosity density will be required to reproduce the Thomson optical depth, to varying degrees depending on the chosen limiting magnitude.

With these features in mind, we have tried a variety of parameterized models for ρ_UV. We will present constraints using a three-parameter model given by

$$\begin{equation} \rho _{\mathrm{UV}}(z) = \rho _{\mathrm{UV},z=4}\left(\frac{z}{4}\right)^{-3} + \rho _{\mathrm{UV},z=7}\left(\frac{z}{7}\right)^{\gamma }. \end{equation} \tag{ 14 }$$

The low-redshift amplitude of this model is anchored by the UV luminosity density at z ∼ 4, ρ_{UV, z = 4} with units of erg s⁻¹ Hz⁻¹ Mpc⁻³. The high-redshift evolution is determined by the normalization ρ_{UV, z = 7} at z ∼ 7, provided in units of erg s⁻¹ Hz⁻¹ Mpc⁻³, and the power-law slope γ. Using this model, we compute the evolving ρ_UV(z) and evaluate the parameter likelihoods as described immediately above.

We have examined other models, including general broken (double) power laws and low-redshift power laws with high-redshift constants or fixed slope power laws. Single power-law models tend to be dominated by the low-redshift decline in the ρ_UV and have difficulty reproducing the Thomson optical depth. Generic broken power-law models have a degeneracy between the low- and high-redshift evolution, even for fixed redshifts about which the power laws are defined, and are disfavored based on their Bayesian information relative to less complicated models that can also reproduce the ρ_UV constraints. The model in Equation (14) is therefore a good compromise between sufficient generality and resulting parameter degeneracies.

However, when using this model, we use a prior to limit γ < 0 to prevent increasing ρ_UV beyond z ∼ 9. The best current constraints on the abundance of z > 8.5 galaxies is from Ellis et al. (2013), where we demonstrated that the highest-redshift galaxies (with M_UV ≲ −19) continue the smooth decline in abundance found at slightly later times. Correspondingly, the constraint on γ amounts to a limit on introducing a new population of galaxies arising at z ≳ 9. The usefulness of this potentiality for the reionization of the universe has been noted elsewhere (e.g., Cen 2003; Alvarez et al. 2012). While imposing no constraint besides γ < 0, we typically find that nearly constant, low-level luminosity densities at high redshift are nonetheless favored (see below). Results for models that feature low-redshift power-law declines followed by low-level constant ρ_UV at high redshift will produce similar quantitative results. All parameterized models we have examined that feature a declining or constant ρ_UV and can reproduce the Thomson optical depth constraint produce similar results, and we conclude that our choice of the exact form of Equation (14) is not critical.

Figure 4 shows constraints on the reionization process calculated using the model described in Section 2. We perform our Bayesian inference modeling assuming M_UV < −17 (close to the UDF12 limit at z ∼ 8, the maximum likelihood model is shown as a dashed line in all panels), M_UV < −10 (an extremely faint limit, with the maximum likelihood model shown as a dotted line in all panels), and an intermediate limit M_UV < −13 (colored regions show 68% credibility regions, while the white lines indicate the maximum likelihood model). We now will discuss the UV luminosity density, stellar mass density, ionized filling fraction, and electron scattering optical depth results in turn.

Zoom In Zoom Out Reset image size

The upper left hand panel shows parameterized models of the UV luminosity density (as given in Equation (14)), constrained by observations of the UV luminosity density (inferred from the measured luminosity functions, see Figures 2 and 3) integrated down to each limiting magnitude. The figure shows error bars to indicate the 68% credibility width of the posterior distributions of ρ_UV at redshifts z ∼ 4–9 integrated to M_UV < −13, but the likelihood of each model is calculated using the full marginalized ρ_UV posterior distributions appropriate for its limiting magnitude. In each case, the ρ_UV(z) evolution matches well the constraints provided by the luminosity function extrapolations, which owes to the well-defined progression of declining ρ_UV with redshift inferred from the UDF12 and earlier data sets. For reference, the maximum likelihood values for the parameters of Equation (14) for M_UV < −13 are log ρ_{UV, z = 4} = 26.50log erg s⁻¹ Hz⁻¹ Mpc⁻³, log ρ_{UV, z = 7} = 25.82log erg s⁻¹ Hz⁻¹ Mpc⁻³, and γ = −0.003.

What drives the constraints on ρ_UV for these limiting magnitudes? Some tension exists at high redshift, as the models prefer a relatively flat ρ_UV(z) beyond the current reach of the data. In each case, the maximum likelihood models become flat at high redshift and reflect the need for continued star formation at high redshift to sustain a low-level of partial IGM ionization. To better answer this question, we have to calculate the full stellar mass density evolutions, reionization histories, and Thomson optical depths of each model.

The lower left hand panel shows how the models compare to the stellar mass densities extrapolated from the results by Stark et al. (2013), using the method described in Section 2. The error bars in the figure reflect the stellar mass densities extrapolated for contributions from galaxies with M_UV < −13, but the appropriately extrapolated mass densities are used for each model. The error bars reflect bootstrap-calculated uncertainties that account for statistically possible variations in the best-fit stellar mass–UV luminosity relation given by Equation (8), and are σ ≈ 0.3 dex at z ∼ 7. To balance the relative constraint from the stellar mass density evolution with the single optical depth constraint (below), we assign the likelihood contributions of the stellar mass density and Thomson optical depth equal weight. The stellar mass densities calculated from the models evolve near the 1σ upper limits of the extrapolated constraints, again reflecting the need for continued star formation to sustain a low-level of IGM ionization.¹¹

Assuming a ratio of Lyman continuum photon production rate to UV luminosity log ξ_ion = 25.2log erg⁻¹ Hz for individual sources consistent with UV slopes measured from the UDF12 data (Dunlop et al. 2012b), an ionizing photon escape fraction f_esc = 0.2, and an intergalactic medium clumping factor of C_H ii = 3, the reionization history Q_H ii(z) produced by the evolving ρ_UV can be calculated by integrating Equation (1) and is shown in the upper right panel of Figure 4. We find that the currently observed galaxy population at magnitudes brighter than M_UV < −17 (dashed line) can only manage to reionize fully the universe at late times z ∼ 5, with the IGM 50% ionized by z ∼ 6 and <5% ionized at z ∼ 12. Contributions from galaxies with M_UV < −13 reionize the universe just after z ∼ 6, in agreement with a host of additional constraints on the evolving ionized fraction (see Section 6 below) and as suggested by some previous analyses (e.g., Salvaterra et al. 2011). These galaxies can sustain a ∼50% ionized fraction at z ∼ 7.5, continuing to ∼12–15% at redshift z ∼ 12. Integrating further down to M_UV < −10 produces maximum likelihood models that reionize the universe slightly earlier.

The Thomson optical depths resulting from the reionization histories can be determined by evaluating Equation (11). To reproduce the nine-year WMAP τ values, most models display low levels of UV luminosity density that persist to high redshift. Even maintaining a flat ρ_UV at the maximum level allowed by the luminosity density constraints at redshifts z ≲ 9, as the maximum likelihood models do, the optical depth just can be reproduced. When considering only the currently observed population M_UV < −17, the UV luminosity density would need to increase toward higher z > 10 redshifts for the universe to be fully reionized by redshift z ∼ 6 and the Thomson optical depth to be reproduced.

We remind the reader that while the exact results for, e.g., the Q_H ii evolution of course depends on the choices for the escape fraction f_esc or the ionizing photon production rate ξ_ion. If f_esc or ξ_ion are lowered, then the evolution of Q_H ii is shifted toward lower redshift. For instance, with all other assumptions fixed, we find that complete reionization is shifted to z ∼ 5 for f_esc = 0.1 and z ∼ 4.75 for log ξ_ion = 25.0log erg⁻¹ Hz. Complete reionization can also be made earlier by choosing a smaller C_H ii or delayed by choosing a larger C_H ii.

Given the above results, we conclude that under our assumptions (e.g., f_esc = 0.2 and C_H ii ≈ 3) galaxies currently observed down to the limiting magnitudes of deep high-redshift surveys (e.g., UDF12) do not reionize the universe alone, and to simultaneously reproduce the ρ_UV constraints, the τ constraint, and to reionize the universe by z ∼ 6 requires yet fainter populations. This conclusion is a ramification of the UDF12 UV spectral slope constraints by Dunlop et al. (2012b), which eliminate the possibility that increased Lyman continuum emission by metal-poor populations could have produced reionization at z > 6 (e.g., Robertson et al. 2010). However, too much additional star formation beyond the M_UV < −13 models shown in Figure 4 will begin to exceed the stellar mass density constraints, depending on f_esc. It is therefore interesting to know whether the M_UV < −13 models that satisfy the ρ_UV, ρ_⋆, and τ constraints also satisfy other external constraints on the reionization process, and we now turn to such an analysis.

The best known Lyα forest constraints come from Fan et al. (2006b), who measured the effective optical depth evolution along lines of sight taken from Sloan Digital Sky Survey (SDSS; including both Lyα and higher-order transitions, where available). Those authors then made assumptions about the IGM temperature and the distribution of density inhomogeneities throughout the universe to infer the evolution of the neutral fraction (these assumptions are necessary for comparing the higher-order transitions to Lyα as well, because those transitions sample different parts of the IGM). The corresponding limits on the neutral fraction according to their model are shown by the filled triangles in Figure 5.

The original data here are the transmission measurements in the different transitions. Transforming those into constraints on the ionized fraction requires a model that: (1) predicts the temperature evolution of the IGM (Hui & Haiman 2003; Trac et al. 2008; Furlanetto & Oh 2009); (2) describes the distribution of gas densities in the IGM (Miralda-Escudé et al. 2000; Bolton & Becker 2009); (3) accounts for spatial structure in the averaged transmission measurements (Lidz et al. 2006); and (4) predicts the topology of ionized and neutral regions at the tail end of reionization (Furlanetto et al. 2004b; Choudhury et al. 2009). These are all difficult, and the Fan et al. (2006b) constraints—based upon a simple semi-analytic model for the IGM structure—can be evaded in a number of ways. In particular, their model explicitly ignores the possibility that Q_H ii < 1, working only with the residual neutral gas inside a mostly ionized medium.

Unsurprisingly, our models do not match these Lyα forest measurements very well. The crude approach of Equation (1) fails to address effects from the detailed density structure of the IGM on the evolution of the mean ionization fraction near the end of the reionization process. Our model is therefore not sufficient to model adequately the tail end of reionization, and only more comprehensive models that account for both radiative transfer effects and the details of IGM structure will realistically model these data at z < 6.

More useful to us is a (nearly) model-independent upper limit on the neutral fraction provided by simply counting the dark pixels in the spectra (Mesinger 2010). McGreer et al. (2011) present several different sets of constraints, depending on how one defines "dark" and whether one uses a small number of very deep spectra (with a clearer meaning to the dark pixels but more cosmic variance) versus a larger set of shallower spectra. Their strongest constraints are roughly 1 − Q_H ii < 0.2 at z = 5.5 and 1 − Q_H ii < 0.5 at z = 6, shown by the open triangles in Figure 5. Interestingly, this model-independent approach permits rather late reionization.

Another use of the Lyα forest is to measure the ionizing background and thereby constrain the emissivity of galaxies: the ionization rate Γ∝λ, where λ is the Lyman continuum mean free path. Bolton & Haehnelt (2007b) attempted such a measurement at z > 5. The method is difficult as it involves interpretation of spectra near the saturation limit of the forest. Nonetheless, the Bolton & Haehnelt (2007b) analysis shows that the ionizing background falls by about a factor of 10 from z ∼ 3 to z ∼ 6 (though see McQuinn et al. 2011). The resulting comoving emissivity is roughly the same as at z ∼ 2, corresponding to 1.5–3 photons per hydrogen atom over the age of the Universe (see also Haardt & Madau 2012).

Additional detailed constraints were derived by Faucher-Giguère et al. (2008), who used quasar spectra at 2 ≲ z ≲ 4 to infer a photoionization rate for intergalactic hydrogen. Combining these measurements with the previous work by Bolton & Haehnelt (2007b) and estimates of the intergalactic mean free path of Lyman continuum photons (Prochaska et al. 2009; Songaila & Cowie 2010), Kuhlen & Faucher-Giguère (2012) calculated constraints on the comoving ionizing photon production rate. These inferred values of $\log \dot{n}_{\mathrm{ion}}< 51 \log \mathrm{s}^{-1} \mathrm{Mpc}^{-3}$ at z ∼ 2–6 are lower than those produced by naively extrapolating typical models of the evolving high-redshift UV luminosity density that satisfy reionization constraints. To satisfy both the low-redshift IGM emissivity constraints and the reionization era constraints then available, Kuhlen & Faucher-Giguère (2012) posited an evolving escape fraction

$$\begin{equation} f_{\mathrm{esc}}(z) = f_{0} \times [(1+z)/5]^{\kappa }. \end{equation} \tag{ 15 }$$

If the power-law slope κ is large enough, then the IGM emissivity constraints at z < 6 can be satisfied with a low f₀ while still reionizing the universe by z ∼ 6 and matching previous WMAP constraints on the electron scattering optical depth. Other previous analyses have emphasized similar needs for an evolving escape fraction (e.g., Ferrara & Loeb 2013; Mitra et al. 2013).

To study the effects of including the IGM emissivity constraints by permitting an evolving f_esc, we repeat our calculations in Section 5, additionally allowing an escape fraction given by Equation (15) with varying f₀ and κ (and maximum escape fraction f_max = 1). With a limiting magnitude M_UV < −13, and utilizing the updated constraints from UDF12 and WMAP, we find maximum likelihood values f₀ = 0.054 and κ = 2.4 entirely consistent with the results of Kuhlen & Faucher-Giguère (2012, see their Figure 7). When allowing for this evolving f_esc (with f_max = 1), we sensibly find that the need for low-level star formation to satisfy reionization constraints is reduced (the maximum likelihood model has a high-redshift luminosity density evolution of ρ_UV∝z^−1.2) and the agreement with the stellar mass density constraints is improved.

However, the impact of the evolving f_esc(z) depends on the maximum allowed f_esc. We allow only f_esc to vary (not the product f_escξ_ion) owing to our UV spectral slope constraints, and the maximum f_esc = 1 (and, correspondingly, the maximum ionizing photon production rate per galaxy) is reached by z ∼ 15. In this model with maximum f_max = 1, the escape fraction at z ∼ 7–8 (where the ionized volume filling factor Q_H ii is changing rapidly) is f_esc ∼ 0.17–0.22, very similar to our fiducial constant f_esc = 0.2 adopted in Section 5. The rate of escaping ionizing photons per galaxy in this evolving f_esc model with maximum f_max = 1 trades off against the high-redshift UV luminosity density evolution. The need for a low-level of star formation out to high redshift is somewhat reduced as the increasing escape fraction with redshift compensates to maintain a partial IGM ionization fraction and recover the observed WMAP τ. Without placing additional constraints on the end of reionization, this evolving f_esc model with f_max = 1 produces a broader redshift range of z ∼ 4.5–6.5 for the completion of the reionization process than the constant f_esc model described in Section 5.

If we instead adopt an evolving f_esc model with f_max = 0.2, then we recover ρ_UV, ρ_⋆, and Thomson optical depth evolutions that are similar to the constant f_esc model from Section 5, but can also satisfy the additional low-redshift IGM emissivity constraints. Compared with the comoving ionizing photon production rates $\dot{n} = [3.2, 3.5, 4.3]\times 10^{50}\ \mathrm{s}^{-1}\ \mathrm{Mpc}^{-3}$ at redshifts z = [4.0, 4.2, 5.0] inferred by Kuhlen & Faucher-Giguère (2012) from measurements by Faucher-Giguère et al. (2008), Prochaska et al. (2009), and Songaila & Cowie (2010), this maximum likelihood model with evolving f_esc recovers similar values ($\dot{n} = [3.3, 3.5, 4.7]\times 10^{50}\ \mathrm{s}^{-1} \ \mathrm{Mpc}^{-3}$ at z = [4, 4.2, 5]). In contrast, the constant f_esc model in Section 5 calculated without consideration of these IGM emissivity constraints would produce $\dot{n}\sim 8\hbox{--}13\times 10^{50}\ \mathrm{s}^{-1}\ \mathrm{Mpc}^{-3}$ at z ∼ 4–5. The reduction of the escape fraction at z ∼ 4–6 compared with the baseline model with constant f_esc allows full ionization of the IGM to occur slightly later in the evolving f_esc model (at redshifts as low as z ≈ 5.3, although in the maximum likelihood model reionization completes within Δz ≈ 0.1 of the redshift suggested by the maximum likelihood constant f_esc model of Section 5).

With a similar range of assumptions to Kuhlen & Faucher-Giguère (2012) for a time-dependence in the escape fraction, the model presented in Section 5 can therefore also satisfy the low-redshift IGM emissivity constraints without assuming an escape fraction of f_esc ∼ 1 at high redshifts. While the evolving covering fraction in galaxy spectra (Jones et al. 2012) provides additional observational support for a possible evolution in the escape fraction at z ≲ 5, such empirical support for evolving f_esc does not yet exist at z ≳ 5.

Another use of the Lyα line is to measure precisely the shape of the red damping wing of the line: the IGM absorption is so optically thick that the shape of this red absorption wing depends upon the mean neutral fraction in the IGM (Miralda-Escudé 1998), though the interpretation depends upon the morphology of the reionization process, leading to large intrinsic scatter and biases (McQuinn et al. 2008; Mesinger & Furlanetto 2008a). Such an experiment requires a deep spectrum of a bright source in order to identify the damping wing.

One possibility is a gamma-ray burst, which has an intrinsic power-law spectrum, making it relatively easy to map the shape of a damping wing. The disadvantage of these sources is that (at lower redshifts) they almost always have damped Lyα (DLA) absorption from the host galaxy, which must be disentangled from any IGM signal (Chen et al. 2007).

To date, the best example of such a source is GRB050904 at z = 6.3, which received rapid follow-up and produced a high signal-to-noise spectrum (Totani et al. 2006). McQuinn et al. (2008) studied this spectrum in light of patchy reionization models. Because it has intrinsic DLA absorption, the constraints are relatively weak: they disfavor a fully neutral IGM but allow Q_H ii ∼ 0.5. We show this measurement with the open circle in Figure 5.

Quasars provide a second possible set of sources. These are much easier to find but suffer from complicated intrinsic spectra near the Lyα line. Schroeder et al. (2013) have modeled the spectra of three SDSS quasars in the range z = 6.24–6.42 and found that all three are best fit if a damping wing is present (see also Mesinger & Haiman 2004). Although the spectra themselves cannot distinguish an IGM damping wing from a DLA, they argue that the latter should be sufficiently rare that the IGM must have Q_H ii ≲ 0.1 at 95% confidence. We show this point as the open square in Figure 5. This conclusion is predicated on accurate modeling of the morphology of reionization around quasars and the distribution of strong IGM absorbers at the end of reionization.

Any ionizing source that turns on inside a mostly neutral medium will carve out an H ii region whose extent depends upon the total fluence of ionizing photons from the source. If one can measure (or guess) this fluence, the extent of the ionized bubble then offers a constraint on the original neutral fraction of the medium. This is, of course, a difficult proposition, as the extremely large Lyα IGM optical depth implies that the spectrum may go dark even if the region is still highly ionized: in this case, you find only a lower limit to the size of the near zone (Bolton & Haehnelt 2007a).

Carilli et al. (2010) examined the trends of near-zone sizes in SDSS quasars from z = 5.8–6.4. The sample shows a clear trend of decreasing size with increasing redshift (after compensation for varying luminosities) by about a factor of two over that redshift range. Under the assumption that near zones correspond to ionized bubbles in a (partially) neutral medium, the volume V ∝ (1 − Q_H ii)⁻¹, or, in terms of the radius, $R_{\rm NZ}^{-3} \propto (1-Q_{\mathrm{H\,{\scriptsize{II}}}})$. In that case, a twofold decrease in size corresponds to an order of magnitude increase in the neutral fraction. However, if the Fan et al. (2006b) Lyα forest measurements are correct, then the neutral fraction at z ∼ 5.8 is so small that the zones are very unlikely to be in this regime. In that case, the trend in sizes cannot be directly transformed into a constraint on the filling factor of ionized gas. We therefore do not show this constraint in Figure 5, as its interpretation is unclear.

The recently discovered z = 7.1 quasar has a very small near-zone (Mortlock et al. 2011). Bolton et al. (2011) used a numerical simulation to analyze it in detail. They concluded the spectrum was consistent both with a small ionized bubble inside a region with Q_H ii ≲ 0.1 and with a highly ionized medium (1 − Q_H ii) ∼ 10⁻⁴–10⁻³ if a DLA is relatively close to the quasar. They argue that the latter is relatively unlikely (occurring ∼5% of the time in their simulations). Further support to the IGM hypothesis is lent by recent observations showing no apparent metals in the absorber, which would be unprecedented for a DLA (Simcoe et al. 2012). A final possibility is that the quasar is simply very young (≲ 10⁶ yr) and has not had time to carve out a large ionized region. We show this constraint, assuming that the absorption comes from the IGM, with the filled square in Figure 5.

Recent small-scale temperature measurements have begun to constrain the contribution of patchy reionization to the kinetic Sunyaev–Zel'dovich (kSZ) effect, which is generated by CMB photons scattering off coherent large-scale velocities in the IGM. These scatterings typically cancel (because any redshift gained from scattering off of gas falling into a potential well is canceled by scattering off gas on the other side), but during reionization modulation by the ionization field can prevent such cancellation (Gruzinov & Hu 1998; Knox et al. 1998). There is also a contribution from nonlinear evolution, see Ostriker & Vishniac (1986). Both the Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT) have placed upper limits on this signal (Dunkley et al. 2011; Reichardt et al. 2012).

Because the patchiness and inhomogeneity of the process induce this signal, the kSZ signal grows so long as this patchy contribution persists. As a result, it essentially constrains the duration of reionization. Using SPT data, Zahn et al. (2012) claim a limit of Δz < 7.2 (where Δz is the redshift difference between Q_H ii = 0.2 and Q_H ii = 0.99) at 95% confidence. The primary limiting factor here is a potential correlation between the thermal Sunyaev–Zel'dovich effect and the cosmic infrared background, which pollutes the kSZ signal. If that possibility is ignored, then the limit on the duration falls to Δz < 4.4. Mesinger et al. (2012) found similar limits from ACT data. They showed that the limit without a correlation is difficult to reconcile with the usual models of reionization. Intensive efforts are now underway to measure the correlation.

We show two examples of the slowest possible evolution allowed by these constraints (allowing for a correlation) with the straight dashed lines in Figure 5. Note that the timing of these curves is entirely arbitrary (as is the shape: there is absolutely no reason to expect a linear correlation between redshift and neutral fraction). They simply offer a rough guide to the maximum duration over which substantial patchiness can persist.

The final class of probes to be considered here relies on Lyα emission lines from galaxies (including some of those cataloged in the HUDF09 and UDF12 campaigns). As these line photons propagate from their source galaxy to the observer, they pass through the IGM and can suffer absorption if the medium has a substantial neutral fraction. This absorption is generally due to the red damping wing, as the sources most likely lie inside of ionized bubbles, so the line photons have typically redshifted out of resonance by the time they reach neutral gas.

One set of constraints come from direct (narrowband) surveys for Lyα emitting galaxies. As the IGM becomes more neutral, these emission lines should suffer more and more extinction, and so we should see a drop in the number density of objects selected in this manner (Santos 2004; Furlanetto et al. 2004a, 2006; McQuinn et al. 2007; Mesinger & Furlanetto 2008b). Such surveys have a long history (e.g., Hu et al. 2002; Malhotra & Rhoads 2004; Santos et al. 2004; Kashikawa et al. 2006). To date, the most comprehensive surveys have come from the Subaru telescope. Ouchi et al. (2010) examined z = 6.6 Lyα emitters and found evidence for only a slight decline in the Lyα transmission. They estimate that Q_H ii ≳ 0.6 at that time. Ota et al. (2008) examined the z = 7 window. Very few sources were detected, which may indicate a rapid decline in the population. They estimate that Q_H II ≈ 0.32–0.64. We show these two constraints with the open pentagons in Figure 5.

The primary difficulty with measurements of the evolution of the Lyα emitter number density is that the overall galaxy population is also evolving, and it can be difficult to determine if evolution in the emitter number counts is due to changes in the IGM properties or in the galaxy population. An alternate approach is therefore to select a galaxy sample independently of the Lyα line (or at least as independently as possible) and determine how the fraction of these galaxies with strong Lyα lines evolves (Stark et al. 2010). Such a sample is available through broadband Lyman break searches with HST and Subaru. Several groups have performed these searches to z ∼ 8 (Fontana et al. 2010; Pentericci et al. 2011; Schenker et al. 2012b; Ono et al. 2012b) and found that although the fraction of Lyα emitters increases slightly for samples of similar UV luminosity from z ∼ 3–6, there is a marked decline beyond z ∼ 6.5 (see also Treu et al. 2012). This decline could still be due to evolutionary processes within the population (e.g., dust content, see Dayal & Ferrara 2012), but both the rapidity of the possible evolution and its reversal from trends now well established at lower redshift make this unlikely. Assuming that this decline can be attributed to the increasing neutrality of the IGM, it requires Q_H ii ≲ 0.5 (McQuinn et al. 2007; Mesinger & Furlanetto 2008b; Dijkstra et al. 2011). We show this constraint with the filled pentagon in Figure 5.

There is one more signature of reionization in the Lyα lines of galaxies. A partially neutral IGM does not extinguish these lines uniformly: galaxies inside of very large ionized bubbles suffer little absorption, while isolated galaxies disappear even when Q_H ii is large. This manifests as a change in the apparent clustering of the galaxies, which is attractive because such a strong change is difficult to mimic with baryonic processes within and around galaxies (Furlanetto et al. 2006; McQuinn et al. 2007; Mesinger & Furlanetto 2008b). Clustering is a much more difficult measurement than the number density, requiring a large number of sources. It is not yet possible with the Lyα line sources at z ∼ 7, but at z ∼ 6.6 there is no evidence for an anomalous increase in clustering (McQuinn et al. 2007; Ouchi et al. 2010), indicating Q_H ii ≳ 0.5 at that time. We show this constraint with the filled diamond in Figure 5.

Figure 5 also shows the reionization history in our preferred model, with the associated confidence intervals, that extrapolates the observed luminosity function to M_UV = −13. Amazingly, this straightforward model obeys all the constraints we have listed in this section, with the exception of the Fan et al. (2006b) Lyα forest measurements. However, we remind the reader that our crude reionization model fails at very small neutral fractions because it does not properly account for the high gas clumping in dense systems near galaxies, so we do not regard this apparent disagreement as worrisome to any degree.

It is worth considering the behavior of this model in some detail. A brief examination of the data points reveals that the most interesting limits come from the Lyα lines inside of galaxies (the filled pentagon, from spectroscopic follow-up of Lyman break galaxies, and the open pentagons, from direct narrowband abundance at z ∼ 6.6 and 7). These require relatively high neutral fractions at z ∼ 7 in order for the IGM to affect the observed abundance significantly; the model shown here just barely satisfies the constraints. In a model in which we integrate the luminosity function to fainter magnitudes (e.g., M_UV < −10, dotted line), the ionized fraction is somewhat higher at this time and may prevent the lines from being significantly extinguished, while a model with the minimum luminosity closer to the observed limit (e.g., M_UV < −17, dashed line) reionizes the universe so late that the emitter population at z ∼ 6.6 should have been measurably reduced.

Clearly, improvements in the measured abundance of strong line emitters at this epoch (perhaps with new multi-object infrared spectrographs, like LUCI on LBT or MOSFIRE on Keck, or with new wide-field survey cameras, like HyperSuprimeCam) will be very important for distinguishing viable reionization histories. Equally important will be improvements in the modeling of the decline in the Lyα line emitters, as a number of factors both outside (the morphology of reionization, resonant absorption near the source galaxies, and absorption from dense IGM structures) and inside galaxies (dust absorption, emission geometry, and winds) all affect the detailed interpretation of the raw measurements (e.g., Santos 2004; Dijkstra et al. 2011; Bolton & Haehnelt 2013).

In order to have a relatively high neutral fraction at z ∼ 7, a plausible reionization history cannot complete the process by z ∼ 6.4 (to which the most distant Lyα forest spectrum currently extends), unless either (1) the Lyman continuum luminosity of galaxies, per unit 1500 Å luminosity, evolves rapidly over that interval, (2) the escape fraction f_esc increases rapidly toward lower redshifts, or (3) the abundance of galaxies below M_UV ∼ −17 evolves rapidly. Otherwise, the UDF12 luminosity functions have sufficient precision to fix the shape of the Q_H ii(z) curve over this interval (indeed, to z ∼ 8), and as shown in Figure 5, the slope is relatively shallow. Indeed, in our best-fit model the universe still has 1.0 − Q_H ii ∼ 0.1 at z ∼ 6. This suggests intensifying searches for the last neutral regions in the IGM at this (relatively) accessible epoch.

Although reionization ends rather late in this best-fit model, it still satisfies the WMAP optical depth constraint (at 1σ) thanks to a long tail of low-level star formation out to very high redshifts, as suggested by Ellis et al. (2013), although the agreement is much easier if the WMAP value turns out to be somewhat high. This implies that the kSZ signal should have a reasonably large amount of power. Formally, our best-fit model has Δz ∼ 5, according to the definition of Zahn et al. (2012), which is within the range that can be constrained if the thermal Sunyaev–Zel'dovich contamination can be sorted out. We also note that the WMAP results (Hinshaw et al. 2012) indicate that when combining the broadest array of available data sets, the best-fit electron scattering optical depth lowers by ∼5% to τ ≈ 0.08.

In summary, we have shown that the UDF12 measurements allow a model of star formation during the cosmic dawn that satisfies all available constraints on the galaxy populations and IGM ionization history, with reasonable extrapolation to fainter systems and no assumptions about high-redshift evolution in the parameters of star formation or UV photon production and escape. Of course, this is not to say such evolution cannot occur (see Kuhlen & Faucher-Giguère 2012, Alvarez et al. 2012, and our Section 6.2 for examples in which it does), but it does not appear to be essential.