A PHYSICAL MODEL FOR THE EVOLVING ULTRAVIOLET LUMINOSITY FUNCTION OF HIGH REDSHIFT GALAXIES AND THEIR CONTRIBUTION TO THE COSMIC REIONIZATION

Our approach is in line with the scenario worked out by Granato et al. (2004) and further elaborated by Lapi et al. (2006, 2011), Mao et al. (2007), and Cai et al. (2013). It exploits the outcomes of many intensive N-body simulations and semi-analytic studies (e.g., Zhao et al. 2003; Wang et al. 2011; Lapi & Cavaliere 2011) according to which pre-galactic halos undergo an early phase of fast collapse, including a few major merger events during which the central regions rapidly reach a dynamical quasi-equilibrium; a subsequent slower accretion phase follows, which mainly affects the halo outskirts. We take the transition between the two phases as the galaxy formation time.

The model envisages that during the fast collapse phase the baryons falling into the newly created potential wells are shock-heated to the virial temperature. This assumption may be at odds with analytic estimates and SPH simulations showing that the amount of shock-heated gas is relatively small for halo masses below ∼10¹² M_☉ (Kereš et al. 2005; Dekel & Birnboim 2006), that is, in the mass range of interest here (see below). However, the analytic estimates assume spherical collapse and the SPH simulations consider a smooth initial configuration. In both cases, the results do not necessarily apply to the halo build up during the fast collapse phase, which is dominated by a few major mergers. A crucial test of the star formation history implied by the model is provided by the data discussed in the following.

An alternative picture for the formation of galaxies is the "cold stream driven" scenario, according to which massive galaxies (halo masses above ∼10¹² M_☉) formed from steady, narrow, cold gas streams, fed by dark matter filaments from the cosmic web, that penetrate the shock-heated media of massive dark matter haloes (Dekel et al. 2009). However, this mechanism may require an implausibly high star formation efficiency (Silk & Mamon 2012; Lapi et al. 2011) and a too high bias factor of galaxies at z ≃ 2 (Davé et al. 2010; Xia et al. 2012). In addition, observational evidence of cold flows has been elusive. The recent adaptive optics–assisted SINFONI observations of a z = 2.3285 star-forming galaxy obtained at the Very Large Telescope (Bouché et al. 2013) showed evidence of cold material at about one third of the virial size of the halo. The gas appears to be co-rotating with the galaxy and may eventually inflow toward the center to feed star formation. However, this gas does not necessarily come from the intergalactic space; it may well be the result of cooling of the gas inside the halo, as implied by our scenario. In any case, as we will see in the following, massive galaxies are only of marginal importance in the present context.

In our framework, the galaxy LF is directly linked to the halo formation rate, d³N(M_vir, τ)/dM_vir dτ dV, as a function of halo mass, M_vir, and cosmic time, τ (see Figure 1). The halo formation rate is approximated by the positive term of the derivative with respect to τ of the cosmological halo mass function (Sasaki 1994). Lapi et al. (2013), using an excursion set approach, showed that this is a good approximation. For the halo mass function, we used the analytic approximation by Sheth & Tormen (1999).

Zoom In Zoom Out Reset image size

In converting the halo formation rate into the UV LF, we also need to take into account the survival time, t_destr, of the halos that are subject to merging into larger halos. A short halo survival time would obviously decrease the number density of halos of given mass that are on at a given time. At the high redshifts of interest here, we expect that a large fraction of halos undergo major mergers, thus losing their individuality; therefore the survival time could impact our estimate of the LF.

The survival time is difficult to define unambiguously. We adopt the value obtained from the negative term of the time derivative, (∂_tN)₋, of the Press & Schechter (1974) or of the Sheth & Tormen (1999) mass functions, that is, t_destr ≡ N/(∂_tN)₋ = t_Hubble(z), independent of halo mass. In our calculations of the UV LF (Equation (9)) we cut the integration over time at t_destr. As we will see in the following in the mass and redshift ranges of interest here the duration of the UV bright phase of galaxy evolution is generally substantially shorter than t_Hubble(z) and therefore also shorter than t_destr. Therefore, the results are only weakly affected by the cut and a more refined treatment of this tricky issue is not warranted.

The history of star formation and accretion into the central black hole (BH) were computed by numerically solving the set of equations laid down in the Appendix of Cai et al. (2013). These equations describe the evolution of the gas along the following lines. Initially, a galactic halo of mass M_vir contains a cosmological gas mass fraction f_b = M_gas/M_vir = 0.17, distributed according to a Navarro et al. (1997) density profile, with a moderate clumping factor and primordial composition. The gas is heated to the virial temperature at the virialization redshift, z_vir, then cools and flows toward the central regions. The cooling is fast, especially in the denser central regions where it triggers a burst of star formation. We adopt an Chabrier (2003) initial mass function (IMF). The gas cooling rate is computed by adopting the cooling function given by Sutherland & Dopita (1993). The evolution of the metal abundance of the gas is followed using the classical equations and stellar nucleosynthesis prescriptions, as reported, for instance, in Romano et al. (2002).

An appreciable amount of cooled gas loses its angular momentum via the drag by the stellar radiation field, settles down into a reservoir around the central supermassive BH, and eventually accretes onto it by viscous dissipation, powering the nuclear activity. The supernova (SN) explosions and the nuclear activity feed energy back to the gaseous baryons, and regulate the ongoing star formation and BH growth.

The feedback from the active nucleus has a key role in quenching the star formation in very massive halos, but is only of marginal importance for the relatively low halo masses of interest here. This implies that the results presented in this paper are almost insensitive to the choice for the parameters governing the BH growth and feedback. The really important model parameters are those that regulate the star formation rate (SFR), and thus the chemical evolution. These are the clumping factor C, which affects the cooling rate; the ratio between the gas inflow timescale, t_cond, and the star formation timescale, t_⋆, s = t_cond/t_⋆; and the SN feedback efficiency, _SN.

For the first two of the latter parameters, clumping factor C and ratio s, —as well as for all the parameters controlling the growth of and the feedback from the active nucleus that, albeit almost irrelevant, are nevertheless included in our calculations—we adopt the values determined by Lapi et al. (2006) and also used by Lapi et al. (2011) and Cai et al. (2013) (i.e., C = 7 and s = 5). A discussion on the plausible ranges for all the model parameters can be found in Cai et al. (2013).

However, we have made a different choice for _SN, which is motivated by the fact that the earlier papers dealt with relatively high halo masses (M_vir ≳ 10^11.2 M_☉) whereas, at the high redshifts of interest here, the masses of interest are substantially lower. As pointed out by Shankar et al. (2006), the data indicate that the star formation efficiency in low-mass halos must be lower than implied by the SN feedback efficiency gauged for higher halo masses. In other words, as extensively discussed in the literature, external (UV background) or internal (SN explosions, radiation from massive low-metallicity stars, and, possibly, stellar winds) heating processes can reduce the SFR in low-mass halos (e.g., Pawlik & Schaye 2009; Hambrick et al. 2011; Finlator et al. 2011a; Krumholz & Dekel 2012; Pawlik et al. 2013; Wyithe & Loeb 2013; Sobacchi & Mesinger 2013) and completely suppress it below a critical halo mass, M_crit.

We take this into account by increasing the SN feedback efficiency with decreasing halo mass and introducing a low-mass cutoff (i.e., considering only galaxies with M_vir ⩾ M_crit). We set

$$\begin{equation} \epsilon _{\rm SN} = 0.05 \lbrace 3 + {\rm erf}[-\log (M_{\rm vir}/10^{11}\,M_\odot)/0.5]\rbrace / 2, \end{equation} \tag{ 1 }$$

which converges, for M_vir > 10¹¹ M_☉, to the constant value, _SN = 0.05, used by Cai et al. (2013).

The critical halo mass can be estimated by equating the gas thermal speed in virialized halos to the escape velocity to find M_crit ≃ 2 × 10⁸(T/2 × 10⁴ K)^3/2 M_☉ (Hambrick et al. 2011), where 2 × 10⁴ K corresponds to the peak of the hydrogen cooling curve. If the balance between heating and cooling processes results in higher gas temperatures, M_crit can increase even by a large factor. An accurate modeling of the physical processes involved is difficult and numerical simulations yield estimates of M_crit/M_☉ in the range from 4 × 10⁸ to ≳ 10⁹ (Hambrick et al. 2011; Krumholz & Dekel 2012; Pawlik et al. 2013; Sobacchi & Mesinger 2013), with a trend toward lower values at z > 6 when the IGM was cooler (e.g., Dijkstra et al. 2004; Kravtsov et al. 2004; Okamoto et al. 2008) and the intensity of the UV background is expected to decrease.

If, on one side, heating properties can decrease the SFR in small galactic halos, the production of UV photons per unit stellar mass is expected to substantially increase with redshift, at least at z ≳ 3. Padoan et al. (1997) argued that the characteristic stellar mass, corresponding to the peak in the IMF per unit logarithmic mass interval, scales as the square of the temperature, T_mc, of the giant molecular clouds within which stars form. But the typical local value of T_mc (∼10 K) is lower than the cosmic microwave background (CMB) temperature at z ≳ 3 implying a rapid increase with z of the characteristic stellar mass, and thus of the fraction of massive stars producing UV photons. On the other hand, a higher fraction of massive stars also implies a more efficient gas heating, which leads to a stronger quenching of the SFR. This complex set of phenomena tends to counterbalance each other. Therefore, in the following we adopt, above M_crit, the SFR given by the model and the UV (ionizing and non-ionizing) photon production rate appropriate to the Chabrier (2003) IMF, allowing for the possibility of a correction factor to be determined by comparison with the data.

The model outlined above allows us to compute the evolution with galactic age of the SFR; the mass in stars; the metal abundance of the gas as a function of the halo mass, M_vir; and the halo virialization redshift, z_vir. Examples illustrating four values of the halo mass at fixed z_vir = 6 (left panels) and for four values of z_vir and two values of the halo mass (right panels) are shown in Figure 2. A few points are worth noting. (1) The active star formation phase in the most massive halos (M_vir ≳ 10¹² M_☉) is abruptly terminated by the active galactic nucleus (AGN) feedback, while in less massive halos, where the feedback is dominated by stellar processes, the SFR declines more gently; (2) at fixed halo mass the SFR is initially higher at higher redshifts (mainly because the higher densities imply faster cooling of the gas) and declines earlier; (3) the higher star formation efficiency at higher z implies that the stellar to halo mass ratio and the metallicity are higher for galaxies that formed at higher redshifts. The different behavior of these quantities for different halo masses and different redshifts determines the relative contributions of the various masses to the LFs and their evolution with cosmic time, as discussed in the following.

Zoom In Zoom Out Reset image size

For comparisons with the data we further need to take into account the dust extinction, which is proportional to the dust column density. In turn, the dust column density is proportional to the gas column density, N_gas, times the metallicity, Z_g, to some power reflecting the fraction of metals locked into dust grains. We thus expect that at our reference UV wavelength, 1350 ${\rm \mathring{A}}$, the dust extinction is $A_{1350}\propto \dot{M}_{\mathrm{\star }}^{\alpha } Z_{\rm g}^\beta$, where $\dot M^{\alpha}_{\star}$ is the star formation rate. Mao et al. (2007) found that a relationship of this kind ($A_{1350}\approx 0.35\,(\dot{M}_{\star }/M_{\odot }\,\mathrm{yr}^{-1})^{0.45}\,(Z_{\rm g}/Z_{\odot })^{0.8}$) does indeed provide a good fit for the luminosity–reddening relation found by Shapley et al. (2001). We have adopted a somewhat different relation

$$\begin{equation} A_{1350}\approx 0.75\,\left(\frac{\dot{M}_{\star }}{M_{\odot }\,\mathrm{yr}^{-1}}\right)^{0.25}\, \left(\frac{Z_{\rm g}}{Z_{\odot }}\right)^{0.3}, \end{equation} \tag{ 2 }$$

which we have found provides a better fit for the UV LFs of LBGs (Reddy & Steidel 2009; Bouwens et al. 2007; McLure et al. 2013, see below) while still being fully consistent with the Shapley et al. (2001) data (see Figure 3), as well as with observational estimates of the escape fractions ($f^{\rm esc}_\lambda =\exp (-A_\lambda /1.08)$, see below) of UV and Lyα photons (Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As illustrated by the bottom left panel of Figure 2, the UV extinction is always low for the least massive galaxies (M_vir ≲ 10¹¹ M_☉), but increases rapidly with increasing M_vir. This implies a strong correlation between dust attenuation and halo/stellar mass, thus providing a physical explanation for the correlation reported by Heinis et al. (2014, see their Figure 3). The bottom right panel of the figure shows that at fixed halo/stellar mass the model implies a mild increase of the attenuation with increasing z, paralleling the small increase in the gas metallicity due to the higher star formation efficiency.

The observed absolute magnitude in the AB system of a galaxy at our reference UV wavelength, λ = 1350 Å, is related to its SFR and to dust extinction A_λ as

$$\begin{equation} M_{1350}^{\rm obs}=51.59 - 2.5\, \log {L_{1350}^{\rm int}\over \mathrm{erg\,s}^{-1}\,\mathrm{Hz}^{-1}} + A_{1350}, \end{equation} \tag{ 3 }$$

with

$$\begin{equation} {L_{1350}^{\rm int}} = k_{\rm UV} {\dot{M}_\star }, \end{equation} \tag{ 4 }$$

where $L_{1350}^{\rm int}$ is the intrinsic monochromatic luminosity at the frequency corresponding to 1350 Å. On account of the uncertainties on the IMF at high z, mentioned above, we treat the normalization factor $k_{\rm UV}\, ({\rm erg}\, {\rm s}^{-1}\, {\rm Hz}^{-1}\, M_\odot ^{-1}\, {\rm yr})$ as an adjustable parameter. Using the evolutionary stellar models of Fagotto et al. (1994a, 1994b), we find that k_UV varies with galactic age and gas metallicity as illustrated by Figure 5. The main contributions to the UV LF come from galactic ages ≳ 10⁷ yr, when log (k_UV) ≳ 27.8, and is somewhat higher for lower metallicity. We adopt log (k_UV) = 28 as our reference value. We note that the relationship between the absolute AB magnitude at the effective frequency ν, M_ν, and the luminosity is $\log [(\nu \,L_\nu)/L_\odot ]=2.3995-\log (\lambda /1350\,{\rm \mathring{A}})-0.4\,M_{\nu }$, where L_☉ = 3.845 × 10³³ erg s⁻¹.

Zoom In Zoom Out Reset image size

We adopt the standard "dust screen" model for dust extinction, so that the observed luminosity is related to the intrinsic one by

$$\begin{equation} L^{\rm obs}(\lambda) = L^{\rm int}(\lambda) 10^{-0.4A_\lambda }, \end{equation} \tag{ 5 }$$

with the Calzetti et al. (2000) reddening curve holding for 0.12 μm ⩽ λ ⩽ 0.63 μm,

$$\begin{eqnarray} A_\lambda /E_s(B-V) &\simeq& 2.659 (-2.156 + 1.509/\lambda \nonumber\\ &&- 0.198/\lambda ^2 + 0.011/\lambda ^3) + R_V, \end{eqnarray} \tag{ 6 }$$

where E_s(B − V) is the color excess of the galaxy stellar continuum and R_V = 4.05. This gives A₁₃₅₀/E_s(B − V) ≃ 11.0 and A₁₂₁₆ ≃ 1.08A₁₃₅₀. Oesch et al. (2013c) found that the dust reddening at z ≃ 4 is better described by a Small Magellanic Cloud (SMC) extinction curve (Pei 1992) rather than by the Calzetti curve. However, for our purposes the only relevant quantity is the A₁₂₁₆/A₁₃₅₀ ratio. Using the Pei (1992) fitting function for the SMC curve, we find (A₁₂₁₆/A₁₃₅₀)_SMC = 1.14, which is very close to the Calzetti ratio (1.08). Thus our results would not appreciably change adopting the SMC law.

Some average properties, which are weighted by the halo formation rate, of galaxies at z_obs = 2, 6, and 12 are shown as a function of the observed luminosity in Figures 6 and 7. The weighted averages have been computed as follows. We sample the formation rate of halos with mass M_vir virializing at t_vir, $\dot{n}(M_{\rm vir}, t_{\rm vir}) \equiv d^3N(M_{\rm vir}, t_{\rm vir})/d\log M_{\rm vir}\,dt\,dV$, in steps of Δlog M_vir = 0.08 and cosmic time interval Δt_vir = 4 Myr, and select those having some property X, like the UV, L₁₃₅₀, luminosity; the Lyα, L₁₂₁₆, luminosity; or the stellar mass (actually we set X = log (quantity)) in a given range, X_i ∈ (X − ΔX/2, X + ΔX/2], at the selected z_obs. The average value, μ_Y(X, ΔX), of some other property Y (e.g., age of stellar populations, t_age, gas metallicity, Z_g, SFR, etc.) is then computed as

$$\begin{eqnarray} &&\mu _Y(X, \Delta X) = \frac{\sum _i Y_i \,\dot{n}_i \, \theta _{\rm H}(X - \Delta X / 2 < X_i \leqslant X + \Delta X / 2)}{ \sum _i \dot{n}_i \, \theta _{\rm H}(X - \Delta X / 2 < X_i \leqslant X + \Delta X / 2)}, \nonumber\\ \end{eqnarray} \tag{ 7 }$$

where θ_H (θ_H(x) = 1 if x is true, θ_H(x) = 0 otherwise) is the Heaviside function. The comoving number density of galaxies satisfying the selection criterion is

$$\begin{eqnarray} n(X, \Delta X) &=& \frac{1}{\Delta X} \sum _i \dot{n}_i \cdot \Delta \log M_{\rm vir} \cdot \Delta t_{\rm vir} \cdot \theta _{\rm H}(X - \Delta X / 2 < X_i \nonumber\\ &\leqslant& X + \Delta X / 2). \end{eqnarray} \tag{ 8 }$$

The data we compare with are specified in the panels and in the captions of the figures. Note that these data were not used to constrain the model parameters. As illustrated by Figure 6, the model compares favorably with observational estimates of various galaxy properties at different redshifts. Stellar population ages are, on average, substantially younger for the brighter galaxies, which are associated with the most massive halos. This is fully consistent with the well established "downsizing" scenario, since it happens because massive galaxies have high SFRs that may reach SFR ∼ 100 M_☉ yr⁻¹ and therefore develop their stellar populations faster. The chemical enrichment and the associated dust obscuration also grow faster in these galaxies that soon become UV-faint. This causes a downward trend of t_age with increasing $L^{\rm obs}_{1350}$, which is consistent with the data (upper left panel of Figure 6).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compare the predictions of our model with observations on the metallicity–luminosity and metallicity–stellar mass relations at different redshifts, we have converted the values of [12 + log (O/H)] given by Maiolino et al. (2008) and Liu et al. (2008) to Z_g in solar units adopting a solar oxygen abundance [12 + log (O/H)]_☉ = 8.69 (Asplund et al. 2009). In Figure 7 the stellar masses estimated by Maiolino et al. (2008) assuming a Salpeter IMF have been divided by a factor of 1.4 to convert them to the Chabrier IMF used in the model. We have also plotted the metallicity estimates by Erb et al. (2006) for a sample of 87 LBGs at z ∼ 2.2, as revised by Maiolino et al. (2008) using an improved calibration. Cullen et al. (2013) found an offset in the relation for z ∼ 2 galaxies compared to Erb et al. (2006). They argue that the difference originates from the use of different metallicity estimators with locally calibrated metallicity relations that may not be appropriate for the different physical conditions of star-forming regions at high redshifts.

The dependence of gas-phase metallicity on stellar mass, with lower M_⋆ galaxies having lower metallicities, is accounted for by the model, at least at z_obs = 2.2 (Figure 7). This trend is also present in the data on z_obs = 3.4 galaxies by Maiolino et al. (2008), which indicate a decrease by 0.6 dex of the amplitude of the Z_g–M_⋆ relation compared to that observed at z_obs = 2.2. Such a strong evolution is at odds with model predictions, according to which the evolution should be quite weak, as shown in the bottom right panel of Figure 7. However, as argued by Cullen et al. (2013), the strong evolution may be an artifact due to the use of locally calibrated metallicity indicators, which do not account for evolution in the physical conditions of star-forming regions.

Dust extinction is differential in luminosity as is clearly seen by comparing the solid black lines (no extinction) in Figure 8 with the dot–dashed blue lines (extinction included). This follows from the faster metal enrichment in the more massive objects (see Figure 2), where the SN and radiative feedbacks are less effective in depressing the SFR.

Figure 2 also shows that the least massive galaxies (M_vir < 10¹¹ M_☉) have low extinction throughout their lifetime. In conjunction with the fast decline of the number density of more massive halos, this results in an increasing dominance of the contribution of low mass galaxies to the UV LF at higher and higher z. Figure 8 indeed shows that the observed portion of the UV LF at z ⩾ 7 is accounted for by galaxies with M_vir < 10^11.2 M_☉, and the effect of dust extinction is negligible over most of the observed luminosity range. This offers the interesting possibility of reconstructing the halo formation rate from the UV LF without the complication of uncertain extinction corrections.

At the faint end the LF has a break corresponding to M_crit, the mass below which the star formation is suppressed by external and internal heating processes (see Section 2). In Figure 8 M_crit = 10^8.5 M_☉, but in the panels comparing the model with the data at various redshifts the thin dotted and solid blue lines show the effect of increasing M_crit to 10^9.8 M_☉ and 10^11.2 M_☉, respectively. This shows that the data only constrain M_crit to be ≲ 10^9.8 M_☉, which is consistent with estimates from simulations as summarized in Section 2.1.

Interestingly, Muñoz & Loeb (2011), using a substantially different model, find that the observed UV LFs at z = 6, 7, and 8 are best fit with a minimum halo mass per galaxy log (M_min/M_☉) = 9.4(+ 0.3, −0.9) in good agreement with our estimate. Whereas these authors simply assume negligible dust extinction at these redshifts, our model provides a quantitative physical account of the validity of this assumption. Raičević et al. (2011), using the Durham GALFORM semi-analytical galaxy formation model, also find that the minimum halo mass that substantially contributes to the production of UV photons at these redshifts is log (M_min/h⁻¹ M_☉) ≃ 9. A lower value, log (M_min/ M_☉) ≃ 8.2, was obtained by Jaacks et al. (2012a) from their cosmological hydrodynamical simulations. Trenti et al. (2010) derive from abundance matching that the galaxies currently detected by HST live in dark matter halos with M_H ≳ 10¹⁰ M_☉, and they predict a weak decrease of the star formation efficiency with decreasing halo mass down to a minimum halo mass for star formation M_H ∼ 10⁸ M_☉, under the assumption that the luminosity function remains a steep power law at the faint end.

The higher feedback efficiency for less massive halos (Equation (1)) makes the slope of the faint end of the UV LF (above the break) flatter than that of the low mass end of the halo formation rate function (as illustrated by Figure 1). However, the former slope reflects to some extent the steepening with increasing redshift of the latter. Just above the low-luminosity break corresponding to M_crit, the slope $-d\log \Phi (L^{\rm obs}_{1350},z)/d\log L^{\rm obs}_{1350}$ (where Φ is per unit $d\log L^{\rm obs}_{1350}$) steadily increases from ≃ 0.5 at z = 2 to ≃ 0.7 at z = 6 and ≃ 1 at z = 10 (see the bottom right panel of Figure 8). The steepening of the faint end of the UV LF with increasing z implies an increasing contribution of low luminosity galaxies to the UV background. The dust extinction, differential with halo mass, further steepens the bright end of the observed UV LF.

As shown by the bottom right panel of Figure 8, the model implies that the evolution of the UV LF from z = 2 to z = 6 is weak. In particular the fast decrease with increasing redshift of the high mass halo formation rate (Figure 1) is not mirrored by the bright end of the LF. This is partly due, again, to the fast metal enrichment of massive galaxies, which translates into a rapid increase of dust extinction and, consequently, in a short duration of their UV bright phase (Figure 2). Thus the contribution to the UV LF of galaxies in halos more massive than ∼10^11.2 M_☉ (thin, solid blue line in Figure 8) decreases rapidly with increasing redshift and galaxies with M_vir ≫ 10^11.2 M_☉ contribute very little. The UV LF at high z is therefore weakly sensitive to the rapid variation of the formation rate of the latter galaxies. The milder decrease of the formation rates of less massive galaxies is largely compensated by the increase of the star formation efficiency due to the faster gas cooling (see Figure 2). At z > 6, however, the decrease of the formation rate even of intermediate mass galaxies prevails and determines a clear negative evolution of the UV LF.

As illustrated by the bottom left panel of Figure 2, the star formation in massive galaxies is dust enshrouded over most of its duration. For example, the figure shows that the extinction at 1350 Å of a galaxy with M_vir = 10¹² M_☉ is already growing at galactic ages of a few times 10⁷ yr, reaching two magnitudes at an age of 4 × 10⁷ yr. The total duration of the star formation phase for these galaxies is ≃ 7 × 10⁸[(1 + z)/3.5]^−1.5 yr (Cai et al. 2013), which implies that they show up primarily at far-IR/(sub-)millimeter wavelengths. Thus the present framework entails a continuity between the high-z galaxy SFR function derived from the UV LF (dominated by low-mass galaxies) and that derived from the infrared (IR; 8–1000 μm) LF (dominated by massive galaxies), investigated by Lapi et al. (2011) and Cai et al. (2013). This continuity is borne out by observational data at z = 2 and 3, as shown by Figure 9. Once proper allowance for the effect of dust attenuation is made, the model accurately matches the IR and the UV LFs and, as expected, UV and IR data cover complementary SFR ranges.

Zoom In Zoom Out Reset image size

The SFR histories inferred from UV and far-IR data are compared in Figure 10 (see also Fardal et al. 2007). This figure shows that the ratio of dust-obscured to unobscured SFR increases with increasing redshift until it reaches a broad maximum at z ∼ 2–3 and decreases afterward. This entails a prediction for the evolution of the IR luminosity density, ρ_IR, beyond z ≃ 3.5, where it cannot yet be determined observationally. At these redshifts, according to the present model, ρ_IR decreases with increasing z faster than the UV luminosity density, ρ_UV. In other words, the extinction correction needed to determine the SFR density from the observed ρ_UV is increasingly small at higher and higher redshifts. This is because the massive halos that can reach high values of dust extinction become increasingly rare at high z.

Zoom In Zoom Out Reset image size

The ionizing photons ($\lambda \leqslant 912\,{\rm \mathring{A}}$) can be absorbed by both dust and H i. About 2/3 of those absorbed by H i are converted into Lyα photons (Osterbrock 1989; Santos 2004), so that the LFs of LAEs provide information on their production rate. The Lyα line luminosity before extinction is then

$$\begin{equation} L^{\mathrm{int}}_{{\mathrm{Ly}}\alpha }= {2\over 3}\,\dot{N}^{\mathrm{int}}_{912 }\,f^{\rm dust}_{912}\, \left(1-f^{{\rm H\,\scriptsize{I}}}_{912} \right) h_{\rm P}\nu _{{\mathrm{Ly}}\alpha }, \end{equation} \tag{ 10 }$$

where $\dot{N}^{\mathrm{int}}_{912 }$ is the production rate of ionizing photons, while $f^{\rm dust}_{912}=\exp (-\tau _{\rm dust,ion})$ and $f^{{\rm H\,\scriptsize{I}}}_{912}=\exp (-\tau _{{\rm H\,\scriptsize{I}}})$ are the fractions that escape absorption by dust and by H i, respectively, and h_P is the Planck constant. We model $\tau _{{\rm H\,\scriptsize{I}}}$ as

$$\begin{equation} \tau _{{\rm H\,\scriptsize{I}}} = \tau ^0_{{\rm H\,\scriptsize{I}}} \left(\frac{\dot{M}_\star }{M_\odot \,\rm yr^{-1}} \right)^{\alpha _{{\rm H\,\scriptsize{I}}}} + \beta _{{\rm H\,\scriptsize{I}}}, \end{equation} \tag{ 11 }$$

where the first term on the right-hand side refers to the contribution from high density star-forming regions, while the second term refers to the contribution of diffuse H i. Since the Calzetti et al. (2000) law does not provide the dust absorption optical depth of ionizing photons, we set

$$\begin{equation} \tau _{\rm dust,ion}=A_{912}/1.08=\gamma \,A_{1350}/1.08 \end{equation} \tag{ 12 }$$

where A₁₃₅₀ is given by Equation (2) and γ is a parameter of the model.

The evolution with the galactic age $\dot{N}^{\mathrm{int}}_{912 }$, for a Chabrier (2003) IMF, a constant SFR of 1 M_☉ yr⁻¹, and three values of the gas metallicity is illustrated by Figure 11. We adopt an effective ratio $k_{\rm ion}\equiv \dot{N}^{\rm int}_{912}/\dot{M}_\star = 4.0 \times 10^{53}\, {\rm photons}\, {\rm s}^{-1}\, M^{-1}_\odot \ {\rm yr}$, appropriate for the typical galactic ages and metallicities of our sources. The figure also shows the evolution of the intrinsic ratio of Lyman-continuum ($\lambda \leqslant 912\,{\rm \mathring{A}}$; $L^{\rm int}_{912} = \dot{N}^{\rm int}_{912} h_{\rm P} \rm \ erg\ s^{-1}\ Hz^{-1}$) to UV luminosity at 1350 Å, $R_{\rm int}\equiv L^{\rm int}_{912}/L^{\rm int}_{1350}$. For our choice of k_ion and k_UV (Equation (4)) we have R_int ≃ 0.265.

Zoom In Zoom Out Reset image size

The interstellar dust attenuates $L^{\mathrm{int}}_{{\mathrm{Ly}}\alpha }$ by a factor $f^{\rm dust}_{{\rm Ly}\alpha }= \exp (-\tau ^{\mathrm{dust}}_{{\mathrm{Ly}}\alpha })$, where $\tau ^{\mathrm{dust}}_{{\mathrm{Ly}}\alpha }$ is the dust optical depth at the Lyα wavelength (1216 Å). The physical processes governing the escape of Lyα photons from galaxies are complex. Dust content, neutral gas kinematics, and the geometry of the neutral gas seem to play the most important roles. For objects with low dust extinction, such as those relevant here, the detailed calculations of the Lyα radiation transfer by Duval et al. (2014, see their Figures 18 and 19) show that the Lyα is more attenuated than the nearby UV continuum by a factor ≃ 2, consistent with the observational result (e.g., Gronwall et al. 2007) that the SFRs derived from the Lyα luminosity are ≃ 3 times lower than those inferred from the rest-frame UV continuum. With reference to the latter result, it should be noted that part of the attenuation of the Lyα luminosity must be ascribed to the IGM (see below) and that the discrepancy between Lyα and UV continuum SFRs may be due, at least in part, to uncertainties in their estimators. A good fit of the Lyα line LFs is obtained setting $f^{\rm dust}_{{\rm Ly}\alpha }\simeq f^{\rm dust}_{1216}/1.6$, where $f^{\rm dust}_{1216} = \exp (-A_{1216}/1.08)$, consistent with the above results.

The fraction $f^{\mathrm{IGM}}_{{\mathrm{Ly}}\alpha }=\exp (-\tau ^{\mathrm{IGM}}_{{\mathrm{Ly}}\alpha })$ of Lyα photons that survive the passage through the IGM was computed following Madau (1995), taking into account only the absorption of the blue wing of the line

$$\begin{equation} f^{\rm IGM}_{\rm Ly\alpha } = 0.5 \lbrace 1 + \exp \, [-0.0036(1+z)^{3.46}]\rbrace. \end{equation} \tag{ 13 }$$

The strong attenuation by dust within the galaxy and by H i in the IGM implies that estimates of the SFR of high-z LBGs and LAEs from the observed Lyα luminosity require careful corrections and are correspondingly endowed with large uncertainties. Conversely, the statistics of LAEs and LBGs at high redshifts are sensitive absorption/extinction probes.

Thus the observed Lyα line luminosity writes

$$\begin{eqnarray} L_{{\mathrm{Ly}}\alpha }^{\mathrm{obs}} &\simeq& 4.36 \times 10^{42}\,\left(\frac{\dot{M}_{\star }}{M_{\odot }\,\mathrm{yr}^{-1}}\right) \nonumber\\ &&\times f^{\rm dust}_{912}\,(1-f^{\mathrm{H\,\scriptsize{I}}}_{912}) \,f^{\mathrm{dust}}_{{\mathrm{Ly}}\alpha }\,f^{\mathrm{IGM}}_{{\mathrm{Ly}}\alpha }\, \mathrm{erg\,s}^{-1}. \end{eqnarray} \tag{ 14 }$$

On the whole, this equation contains four parameters: the three parameters in the definition of $\tau _{{\rm H\,\scriptsize{I}}}$ (Equation (11)), that is, $\tau ^0_{{\rm H\,\scriptsize{I}}}$, $\alpha _{{\rm H\,\scriptsize{I}}}$, and $\beta _{{\rm H\,\scriptsize{I}}}$, plus γ (Equation (12)). We have attempted to determine their values fitting the Lyα line LFs by Blanc et al. (2011) at z ∼ 3 and by Ouchi et al. (2008) at z ∼ 3.8. However, we could not find an unambiguous solution because of the strong degeneracy among the parameters. To break the parameters' degeneracy, we fixed $\alpha _{{\rm H\,\scriptsize{I}}}= 0.25$, analogous to Equation (2), and $\tau ^0_{{\rm H\,\scriptsize{I}}} = 0.3$. Furthermore, we discarded the solutions that imply too high emission rates of ionizing photons and too low electron scattering optical depth (see below). Based on these somewhat loose criteria, we chose $\beta _{{\rm H\,\scriptsize{I}}}\simeq 1.5$ and γ ≃ 0.85.

A check on the validity of our choices is provided by recent Lyman-continuum imaging of galaxies at z ≃ 3. Nestor et al. (2013) measured the average ratios of non-ionizing UV to Lyman-continuum flux density corrected for IGM attenuation, η_obs, for a sample of LBGs and a sample of LAEs, both at z ∼ 3. Such ratios were found to be $\eta _{\rm LBG}=18.0^{+34.8}_{-7.4}$ for LBGs (rest-frame UV absolute magnitudes −22 ⩽ M_UV ⩽ −20) and $\eta _{\rm LAE}=3.7^{+2.5}_{-1.1}$ for LAEs (−20 < M_UV ⩽ −18.3). Mostardi et al. (2013) probed the Lyman-continuum spectral region of 49 LBGs and 70 LAEs spectroscopically confirmed at z ∼ 2.85, as well as 58 LAE photometric candidates in the same redshift range. After correcting for foreground galaxy contamination and H i absorption in the IGM, the average values for their samples are η_LBG = 82 ± 45, η_LAE = 7.6 ± 4.1 for the spectroscopic sample, and η_LAE = 2.6 ± 0.8 for the full LAE sample.

The observed ratios are equal to the intrinsic ones $\eta _{\rm int} \simeq R_{\rm int}^{-1}$ attenuated by dust and H i absorption (the latter only for ionizing photons)

$$\begin{equation} \eta _{\rm obs}=\eta _{\rm int}{f^{\rm dust}_{1350}\over f^{\rm dust}_{912}\cdot f^{\mathrm{H\,\scriptsize{I}}}_{912}}. \end{equation} \tag{ 15 }$$

The intrinsic ratio we have adopted, $\eta _{\rm int} \simeq R_{\rm int}^{-1}\simeq 3.77$, is within the range measured for LAEs in both the Nestor et al. (2013) and the Mostardi et al. (2013) samples, implying that the attenuation by dust and H i is small, as expected in the present framework (see Figure 6). LBGs in both samples have SFRs in the range of 5–250 M_☉ yr⁻¹, with median values around 50 M_☉ yr⁻¹, and gas metallicities Z_g ≃ 0.7 ± 0.3 Z_☉. After Equation (15) the optical depth for the absorption of ionizing photons by H i is $\tau _{{\rm H\,\scriptsize{I}},\rm LBG} = \ln (\eta _{\rm obs,LBG}/\eta _{\rm int}) - \ln (f^{\rm dust}_{1350}/f^{\rm dust}_{912}) = \ln (\eta _{\rm obs,LBG}/\eta _{\rm int}) - (\gamma - 1) A_{1350} / 1.08$, giving $\tau _{{\rm H\,\scriptsize{I}},\rm LBG} \simeq 1.8^{+1.3}_{-0.6}$ for $\eta _{\rm obs,LBG} = 18.0^{+34.8}_{-7.4}$ (Nestor et al. 2013) and $\tau _{{\rm H\,\scriptsize{I}},\rm LBG} \simeq 3.3^{+0.6}_{-0.9}$ for η_{obs, LBG} = 82 ± 45 (Mostardi et al. 2013). With our choice for the parameters, Equation (11) gives $\tau _{{\rm H\,\scriptsize{I}}} = \tau ^0_{{\rm H\,\scriptsize{I}}} \dot{M}_\star ^{\alpha _{{\rm H\,\scriptsize{I}}}} + \beta _{{\rm H\,\scriptsize{I}}} \simeq 2.3^{+0.4}_{-0.4}$, which is consistent with both observational estimates.

As illustrated by Figure 12, the Lyα line LFs yielded by the model compare quite well with observational determinations at several redshifts, from z = 3 to z = 7.7. The bottom right panel shows that the intrinsic evolution of the Lyα line LF is remarkably weak from z = 2 to z = 6, it is even weaker than in the UV case (Figure 8), and similarly to the case of the UV LF, its faint portion is predicted to get steeper with increasing redshift. The observed evolution at high-z is more strongly negative than the intrinsic one due to the increasing attenuation by the IGM. The figure also demonstrates that, at z ⩾ 5.7, observational estimates based on photometric samples (open symbols), only partially confirmed in spectroscopy, tend to be systematically higher than those based on purely spectroscopic samples (filled symbols). Therefore, more extensive spectroscopic confirmation is necessary before firm conclusions on the high-z evolution of the Lyα line LF can be drawn.

Zoom In Zoom Out Reset image size

The average properties, weighted by the halo formation rate, of LAEs at z = 2, 6, and 12 are shown as a function of the observed Lyα luminosity in the right-hand panels of Figure 6. Compared to LBGs (left panels of the same figure), LAEs have somewhat younger ages, implying somewhat lower stellar masses, SFRs, and metallicities at given M_vir. The latter two factors combine to give a substantially lower dust extinction.

The injection rate of ionizing photons into the IGM is

$$\begin{equation} \dot{N}^{\rm esc}_{912}=\dot{N}^{\rm int}_{912} f^{\rm esc}_{912} = k_{\rm ion} {\dot{M}_\star } f^{\rm esc}_{912} \; {\rm photons}\;{\rm s}^{-1}, \end{equation} \tag{ 16 }$$

where $f^{\rm esc}_{912} \equiv f^{\rm dust}_{912} \times f^{{\rm H\,\scriptsize{I}}}_{912}$ is the fraction of ionizing photons emerging at the galaxy boundary. The dependencies on UV magnitude and redshift of $f^{\rm dust}_{912}$, $f^{{\rm H\,\scriptsize{I}}}_{912}$, and $f^{\rm esc}_{912}$, weighted by the halo formation rate, are illustrated in Figure 13.

Zoom In Zoom Out Reset image size

As shown by Equations (2) and (11), the optical depths for absorption by both dust and H i (and thus the corresponding escape fractions) are determined by the intrinsic properties of the galaxies (the SFR and, in the case of dust absorption, the metallicity), mostly controlled by the halo mass. They are thus weakly dependent on redshift. Lower-mass galaxies have lower SFRs and, more importantly, lower metallicities. This translates into substantially lower optical depths (<1) (i.e., substantially higher escape fractions). For brighter galaxies, which have optical depths ⩾1, the escape fractions are exponentially sensitive to the weak redshift dependence of the metallicity at a given halo mass (see Figure 2). Thus the small decrease with z of the metallicity translates, for bright galaxies, into a significant increase of $f^{\rm dust}_{912}$. These luminosity and redshift dependencies do not apply to $f^{{\rm H\,\scriptsize{I}}}_{912}$, which, not being affected by metallicity, is almost redshift independent at all luminosities. In the bottom right panel of Figure 13 the escape fractions of ionizing photons given by the model for two ranges of observed UV luminosity, are compared with observational estimates at several redshifts. The agreement is generally good.

The redshift-dependent emission rate function, $\phi (\log \dot{N}^{\rm esc}_{912}, z)$, can be constructed in the same way as the UV LF (Section 3). The average injection rate of ionizing photons into the IGM per unit comoving volume at redshift z is

$$\begin{equation} \langle \dot{N}^{\rm esc}_{912} \rangle (z) = \int \dot{N}^{\rm esc}_{912} \phi (\log \dot{N}^{\rm esc}_{912}, z) d\log \dot{N}^{\rm esc}_{912}. \end{equation} \tag{ 17 }$$

Figure 14 compares the average injection rate of ionizing photons into the IGM, $\langle \dot{N}^{\rm esc}_{912}\rangle (z)$, as a function of redshift yielded by the model for two choices of the critical halo mass (M_crit = 10^8.5 M_☉ and M_crit = 10¹⁰ M_☉) with observational estimates. The original data refer to different quantities such as the proper hydrogen photoionization rate, $\Gamma _{{\rm H\,\scriptsize{I}}}(z)$; the average specific intensity of UV background, J_ν(z); and the comoving spatially averaged emissivity, _ν(z). The conversion of these quantities into $\langle \dot{N}^{\rm esc}_{912}\rangle (z)$ was done using the formalism laid down by Kuhlen & Faucher-Giguère (2012). The values of $\langle \dot{N}^{\rm esc}_{912}\rangle$ obtained from the model are above the estimates by Kuhlen & Faucher-Giguère (2012), but consistent with (although on the high side of) the more recent results by Becker & Bolton (2013) and Nestor et al. (2013).

Zoom In Zoom Out Reset image size

The reionization of the IGM is described in terms of the evolution of the volume filling factor of H ii regions, $Q_{{\rm H\,\scriptsize{II}}}(t)$, which is ruled by (Madau et al. 1999)

$$\begin{equation} \frac{dQ_{{\rm H\,\scriptsize{II}}}}{dt} \simeq \frac{\langle \dot{N}^{\rm esc}_{912}\rangle }{\bar{n}_{\rm H}} - \frac{Q_{{\rm H\,\scriptsize{II}}}}{\bar{t}_{\rm rec}}, \end{equation} \tag{ 18 }$$

where $\bar{n}_{\rm H} = X \rho _{\rm c} \Omega _{\rm b} / m_{\rm p} \simeq 2.5 \times 10^{-7} X (\Omega _{\rm b} h^2/0.022)$ cm⁻³ is the mean comoving hydrogen number density in terms of the primordial mass fraction of hydrogen X = 0.75, of the present day critical density ρ_c = 1.878h² × 10⁻²⁹ g cm⁻³, and of the proton mass m_p = 1.67 × 10⁻²⁴ g. The mean recombination time is given by (Madau et al. 1999; Kuhlen & Faucher-Giguère 2012)

$$\begin{eqnarray} \bar{t}_{\rm rec} & = & \frac{1}{f_e \bar{n}_{\rm H} (1+z)^3 \alpha _{\rm B}(T_0) C_{{\rm H\,\scriptsize{II}}}} \nonumber \\ &=& {0.51\over f_e} \left (\frac{\Omega _{\rm b}}{0.049} \right)^{-1} \left(\frac{1+z}{7} \right)^{-3} \left(\frac{T_0}{2 \times 10^4\ {\rm K}} \right)^{0.7} \nonumber\\ &&\times \left(\frac{C_{{\rm H\,\scriptsize{II}}}}{6} \right)^{-1}\ {\rm Gyr}, \end{eqnarray} \tag{ 19 }$$

where f_e is the number of free electrons per hydrogen nucleus in the ionized IGM, assumed to have a temperature T₀ = 2 × 10⁴ K, and $C_{{\rm H\,\scriptsize{II}}} \equiv \langle \rho ^2_{{\rm H\,\scriptsize{II}}} \rangle / \langle \rho _{{\rm H\,\scriptsize{II}}} \rangle ^2$ is the clumping factor of the ionized hydrogen. f_e depends on the ionization state of helium; we assumed it to be doubly ionized (f_e = 1 + Y/2X ≃ 1.167) at z < 4 and singly ionized (f_e = 1 + Y/4X ≃ 1.083) at higher redshifts (Robertson et al. 2013).

The clumping factor has been extensively investigated using numerical simulations (see Finlator et al. 2012 for a critical discussion of earlier work). A series of drawbacks have been progressively discovered and corrected. The latest studies accounting for the photo-ionization heating that tends to smooth the diffuse IGM, and for the IGM temperature, which could suppress the recombination rate further, generally find relatively low values of $C_{{\rm H\,\scriptsize{II}}}$ (Pawlik et al. 2009; McQuinn et al. 2011; Shull et al. 2012; Finlator et al. 2012; Kuhlen & Faucher-Giguère 2012; see the upper left panel of Figure 15). Alternatively, the clumping factor can be computed as the second moment of the IGM density distribution, integrating up to a maximum overdensity (Kulkarni et al. 2013). We adopt, as our reference, the model $C_{\rm H\,\scriptsize{II},\rm T_b, x_{{\rm H\,\scriptsize{II}}}>0.95}$ by Finlator et al. (2012), but we will also discuss the effect of different choices.

Zoom In Zoom Out Reset image size

Additional constraints on the reionization history are set by the electron scattering optical depth, τ_es, measured by CMB anisotropy experiments. The optical depth of electron scattering up to redshift z is given by

$$\begin{eqnarray} \tau _{\rm es}(\leqslant z) &=& \int ^z_0 dz^{\prime } \left| \frac{dt}{dz^{\prime }} \right| c \sigma _{\rm T} n_e(z^{\prime })\nonumber\\ &=& \int ^z_0 dz^{\prime } \frac{c (1+z^{\prime })^2}{H_0 \sqrt{\Omega _{\Lambda } + \Omega _{\rm m} (1+z^{\prime })^3}} f_e Q_{{\rm H\,\scriptsize{II}}}(z^{\prime }) \sigma _{\rm T} \bar{n}_{\rm H}, \nonumber\\ \end{eqnarray} \tag{ 20 }$$

where $n_e \simeq f_e Q_{{\rm H\,\scriptsize{II}}} \bar{n}_{\rm H} (1+z)^3$ is the mean electron density, c is the speed of light, and σ_T ≃ 6.65 × 10⁻²⁵ cm² is the Thomson cross section. WMAP 9 yr data alone give τ_es = 0.089 ± 0.014, which slightly decrease to τ_es = 0.081 ± 0.012 when they are combined with external data sets (Hinshaw et al. 2013). The combination of the Planck temperature power spectrum with a WMAP polarization, low-multipole likelihood results in an estimate closely matching the WMAP 9 yr value, $\tau _{\rm es}=0.089^{+0.012}_{-0.014}$ (Planck Collaboration 2013b). However, replacing the WMAP polarized dust template with the far more sensitive Planck/HFI 353 GHz polarization map lowers the best fit value to τ_es = 0.075 ± 0.013 (Planck Collaboration 2013a); this result, however, has to be taken as preliminary.

The evolution of $Q_{{\rm H\,\scriptsize{II}}}(t)$ and the electron scattering optical depth τ_es(⩽ z) yielded by the model adopting the critical halo mass M_crit = 10^8.5 M_☉, the effective escape fraction $f^{\rm esc}_{912}$ laid down before, and the evolutionary law for the IGM clumping factor $C_{{\rm H\,\scriptsize{II}}}$ by Finlator et al. (2012), are shown by the thick, solid black lines in the main panel of Figure 15. Although this model gives a τ_es consistent with the determination by Planck Collaboration (2013b) and less than 2σ below those by Hinshaw et al. (2013) and Planck Collaboration (2013b), it yields a more extended fully ionized phase than indicated by the constraints on $Q_{{\rm H\,\scriptsize{II}}}$ at z ⩾ 6, which were collected by Robertson et al. (2013) who, however, cautioned that they are all subject to substantial systematic or modeling uncertainties. Indications pointing to a rapid drop of the ionization degree above z ≃ 6 include hints of a decrease of the comoving emission rates of ionizing photons (see Figure 14), sizes of quasar near zones, and the Lyα transmission through the IGM (see Robertson et al. 2013 for references).

As we have seen before, the observed UV LFs imply M_crit ≲ 10¹⁰ M_☉. Adopting the latter value and keeping our baseline $f^{\rm esc}_{912}(z)$ and $C_{{\rm H\,\scriptsize{II}}}(z)$ shortens the fully ionized phase that now extends only up to z ∼ 6 (dotted red line in the main figure of Figure 15), lessening the discrepancy with constraints on $Q_{{\rm H\,\scriptsize{II}}}$ at the cost of τ_es, dropping almost ≃ 3 σ below the best fit WMAP value and 2 σ below the best fit value of Planck Collaboration (2013a). This conclusion is unaffected by different choices for $C_{{\rm H\,\scriptsize{II}}}(z)$, as illustrated in the upper left panel of Figure 15, as long as we keep our baseline $f^{\rm esc}_{912}(z)$. This is because the production rate of ionizing photons (first term on the right-hand side of Equation (18)) is always substantially larger than the recombination rate.

The minimum SFR density required to keep the universe fully ionized at the redshift z is (Madau et al. 1999)

$$\begin{equation} \dot{\rho }_\star \simeq 8.2 \times 10^{-4} \left (\frac{C_{{\rm H\,\scriptsize{II}}}}{f_{\rm esc}} \right) \left (\frac{\Omega _{\rm b} h^2}{0.022} \right)^2 \left (\frac{1+z}{7} \right)^3\ M_\odot \ {\rm yr}^{-1}\ {\rm Mpc}^{-3}. \end{equation} \tag{ 21 }$$

This is shown in Figure 10 (gray area) for $3 \lesssim C_{{\rm H\,\scriptsize{II}}}/f_{\rm esc} \lesssim 30$. A similar figure was presented by Finkelstein et al. (2012, their Figure 3). A comparison of their green curve with our blue curve illustrates the difference between the UV luminosity density implied by our model for sources brighter than M_UV = −18 and that from the hydrodynamic simulations of Finlator et al. (2011a).

A sort of tradeoff between the constraints of $Q_{{\rm H\,\scriptsize{II}}}$ and those of τ_es is obtained if the cooling rate of H ii increases rapidly with decreasing z for z ≲ 8. This might be the case if, for example, the clumping factor climbs in this redshift range, as in the C₋₁ L6N256 no-reheating simulation by Pawlik et al. (2009; dot–dashed blue line in the main figure of Figure 15). The drawbacks suffered by these simulations were, however, pointed out by Pawlik et al. (2009) and Finlator et al. (2012).

The tension between the determinations of τ_es(⩽ z) and constraints on $Q_{{\rm H\,\scriptsize{II}}}(z)$ has been repeatedly discussed in recent years (e.g., Kuhlen & Faucher-Giguère 2012; Haardt & Madau 2012; Robertson et al. 2013). Our conclusions are broadly consistent with the earlier ones, as can be seen by all models matching the observed high-z UV LFs. However, our model differs from the others because our LFs and the escape fraction as a function of redshift are physically grounded, whereas the quoted models are based on phenomenological fits to the data. This means that our model is more constrained; for example, the extrapolations of the LFs outside the luminosity and redshift ranges covered by observations come out from our equations, rather than being controlled by adjustable parameters. In other words we explore different parameter spaces. As a result, we differ in important details such as the slope of the faint tail of the UV LF and its redshift dependence, the total number density of high-z UV galaxies, and the redshift dependence of the escape fraction of ionizing photons.

The need for a redshift- or mass dependence of $f^{\rm esc}_{912}$ has also been deduced by other studies (e.g., Alvarez et al. 2012; Mitra et al. 2013). Mitra et al. (2013) combined data-constrained reionization histories and the evolution of the LF of early galaxies to find an empirical indication of a 2.6 times increase of the average escape fraction from z = 6 to z = 8. Alvarez et al. (2012) argued that there are both theoretical and observational indications that $f^{\rm esc}_{912}$ is higher at lower halo masses and proposed that a faint population of galaxies with host halo masses of ∼10^8–9 M_☉ dominated the ionizing photon budget at redshifts z ≳ 9 due to their much higher escape fractions, which were again empirically estimated. Our model provides a physical basis for the increase of the effective $f^{\rm esc}_{912}$ with mass and redshift.

The present data do not allow us to draw any firm conclusion on the reionization history because they may be affected by substantial uncertainties. Those uncertainties on data at z = 6–7 are discussed by Robertson et al. (2013). Those on τ_es are illustrated by the finding (Planck Collaboration 2013a) that different corrections for the contamination by polarized foregrounds may lower the best fit value by about 1 σ.