The Star Formation Efficiency of Galaxies During the Epoch of Reionization

Constraints on the Star Formation Efficiency of Galaxies During the Epoch of Reionization


Reionization is thought to have occurred in the redshift range of , which is now being probed by both deep galaxy surveys and CMB observations. Using halo abundance matching over the redshift range and assuming smooth, continuous gas accretion, we develop a model for the star formation efficiency of dark matter halos at that matches the measured galaxy luminosity functions at these redshifts. We find that peaks at at halo masses  M, in qualitative agreement with its behavior at lower redshifts. We then investigate the cosmic star formation histories and the corresponding models of reionization for a range of extrapolations to small halo masses. We use a variety of observations to further constrain the characteristics of the galaxy populations, including the escape fraction of UV photons. Our approach provides an empirically-calibrated, physically-motivated model for the properties of star-forming galaxies sourcing the epoch of reionization. In the case where star formation in low-mass halos is maximally efficient, an average escape fraction can reproduce the optical depth reported by Planck, whereas inefficient star formation in these halos requires either about twice as many UV photons to escape, or an escape fraction that increases towards higher redshifts. Our models also predict how future observations with JWST can improve our understanding of these galaxy populations.

early universe, reionization–galaxies: evolution–galaxies: high-redshift

1 Introduction

The epoch of reionization, occurring approximately 500 million to 1 billion years after the Big Bang (), was the last major phase transition in the history of the universe, during which the neutral intergalactic medium (IGM) gradually transformed into the highly ionized state that we observe today. At the same time, the first stars and galaxies were forming from primordial gas clouds in the Universe. Although there is a clear connection between these two events (Barkana & Loeb 2001; Bromm & Yoshida 2011; Loeb & Furlanetto 2013; Robertson et al. 2013), their detailed relation is unknown, thanks to uncertainties in the properties of the galaxy populations. Fortunately, understanding the epoch of reionization itself will also shed light on the formation and evolution of early galaxies.

A number of diverse observational probes have helped us develop a preliminary picture of reionization over the past decade. Direct imaging of high-redshift star-forming galaxies has measured the abundance, luminosity distribution, and emission lines of galaxies out to , from which valuable information about the budget of ionizing photons might be drawn (e.g., Ellis et al. 2013; Oesch et al. 2015; Bouwens et al. 2015; Robertson et al. 2015; Atek et al. 2015). Attempts have also been made to measure the ionization state of the IGM directly. Reionization is relatively well-constrained to be completed at by analysis of the Gunn-Peterson effect (Gunn & Peterson, 1965) in the spectra of high-redshift QSOs (e.g., Fan et al. 2006; Willott et al. 2007; Bolton et al. 2011; McGreer et al. 2015) and gamma-ray bursts (GRBs, e.g., Chornock et al. 2013). Preliminary measurements of the IGM neutral fraction at –8 have been made using the Ly emission fraction of star-forming galaxies (e.g., Stark et al. 2010; Treu et al. 2013; Pentericci et al. 2014; Tilvi et al. 2014; Faisst et al. 2014; Schenker et al. 2014). On the other hand, the integrated Thomson scattering optical depth, , due to the scattering of CMB photons off free electrons after reionization, places an important integral constraint on the extended reionization history. The recently reported value by the Planck Collaboration et al. (2015) is significantly lower than previously measured by the Wilkinson Microwave Anisotropy Probe (WMAP) satellite (Hinshaw, 2013).4 It is also necessary for these early ionizing sources to map smoothly onto their post-reionization counterparts at , whose contribution to the global ionizing background has been characterized through measurements of the Ly forest (Kuhlen & Faucher-Giguère, 2012; Becker & Bolton, 2013).

Meanwhile, theorists have been trying to to develop plausible models for reionization that synthesize the results of these observations. For many years, the major challenge was reconciling the large optical depth observed by WMAP with the ionizing emissivity provided by galaxies at redshift (Robertson et al., 2013). Both extrapolations of the estimated post-reionization ionizing background at (Bolton & Haehnelt, 2007; Calverley et al., 2011) and measurements of the star formation rate evolution of Lyman break galaxies (LBGs, e.g., Bouwens et al. 2012; Schenker et al. 2013; Oesch et al. 2015) suggested a paucity of ionizing photons necessary for reproducing . Resolving this dilemma required adjusting one or more of the many unknown parameters that map galaxy observations to reionization: the escape fraction of UV photons, , the minimum halo mass to host a galaxy, , or other, more subtle parameters (Bromm & Yoshida, 2011; Kuhlen & Faucher-Giguère, 2012). For example, a common solution was to increase the overall ionizing emissivity by assuming an evolving escape fraction, , of ionizing photons (e.g., Alvarez et al. 2012; Kuhlen & Faucher-Giguère 2012), which is motivated by both numerical simulations of star formation in the high-redshift universe (Wise & Cen, 2009; Ferrara & Loeb, 2013) and the wide range of values estimated from galaxy observations at lower redshifts (Siana et al., 2010; Nestor et al., 2013; Mostardi et al., 2015). While the tension between the optical depth and ionizing emissivity has been largely resolved by Planck (Robertson et al., 2015; Bouwens et al., 2015), there remains a great deal of freedom in setting these important parameters of both galaxy formation and reionization.

The primary goal of this paper is to investigate how current observations inform models of high-redshift galaxies and cosmic reionization. We will construct a simple analytic model of galaxy formation whose parameters can be constrained by existing observations but that also allows a theoretically-motivated extrapolation to earlier epochs. The halo abundance matching approach (Vale & Ostriker, 2004), which empirically associates galaxies and dark matter halos by matching their number densities, has been widely used to model the star formation activity in galaxies across cosmic time (Conroy & Wechsler, 2009; Yang et al., 2012; Behroozi & Silk, 2015). Taking advantage of its simplicity for analytic models, several recent studies have demonstrated its utility for empirically predicting the luminosity functions of high-redshift galaxies responsible for reionization (Mashian et al., 2016; Mason et al., 2015; Visbal et al., 2015).

In this paper, we will first give a more thorough analysis of the star formation efficiency implied by abundance matching, assuming continuous star formation in high-redshift galaxies. In contrast to most treatments of reionization, our approach allows the star formation efficiency to evolve with both halo mass and redshift and thus effectively alters the overall ionizing efficiency. Whereas many analytic models simply treat as an arbitrary constant, our approach carefully constrains with abundance matching, providing a more accurate way of calibrating models of reionization to observations. Because it is calibrated to bright galaxies at –8, our approach will provide a “baseline” model for as-yet-undiscovered galaxy populations during the “Cosmic Dawn”: while the extrapolation (to both higher redshifts and fainter luminosities) is by no means a firm theoretical prediction, it at least provides a clear understanding of the role of known galaxy populations in inferences about the reionization process. Here we explore these implications in some detail through comparisons of our extrapolations to several other observables, and we comment on how additional populations or effects (such as feedback or Population III stars) may affect the results.

This paper is organized as follows. In Section 2, we describe our simple analytic model for continuous mode star formation in high-redshift galaxies based on the dust-corrected, rest-frame UV luminosity functions and the mass assembly history of dark matter halos. We also provide a brief discussion of the halo abundance matching (HAM) technique. We compare our modeled star formation histories to observations of high-redshift galaxies and present our predictions for JWST in Section 3. Then, in Section 4, we show how to use our star formation models to calculate simple reionization histories. Section 5 is dedicated to a joint analysis of galaxy populations and observational constraints, in particular the optical depth measured by Planck. In Section 6, we discuss a few possible variations around our baseline model, including physical processes such as photo-suppression of low-mass galaxies and the formation of very massive Pop III stars. Finally, we briefly conclude in Section 7.

Throughout this paper, we assume a flat, cosmology with , , , , , and , consistent with the most recent measurement from Planck (Planck Collaboration et al., 2015). We also assume a Salpeter initial mass function (IMF) from 0.1 to 100 and a stellar metallicity of for calculations of star formation and reionization. All the magnitudes cited are expressed in terms of (Oke, 1974) and the base-10 logarithm is assumed unless stated otherwise.

2 Modeling the Cosmic Star Formation History

We use a simple but physically intuitive approach to model the average star formation rate (SFR) inside a galaxy at redshift . When the timescale of star formation is much less than the dynamical time on which galaxies grow (i.e. stars form “instantaneously”), the star formation rate can be approximated as a balance between the gas fueling rate, , and the rate, , at which galactic outflows deplete the gas (Muñoz, 2012). We introduce a free parameter, , which measures the star formation efficiency, namely, the fraction of accreted baryons that forms stars:5


where is the total halo accretion rate.

As a result, the star formation history of galaxies can be directly related to the mass assembly history of dark matter halos. The evolution of has been studied using catalogs of halos from cosmological simulations (see Section 2.1), allowing us to describe the SFR once can be characterized. For high-redshift galaxies dominated by massive, young stars, the rest-frame UV 1500  luminosity is a good indicator of the SFR if any significant fluctuations in the star formation rate occur on a timescale longer than the evolution timescale of the massive stars, which is about 100 Myr (Kennicutt, 1998; Madau & Dickinson, 2014), and provided that dust attenuation is properly accounted for. Specifically, the dust-corrected SFR of a galaxy is proportional to its rest-frame UV continuum () luminosity by


where we assume a fiducial constant evaluated for continuous mode star formation with a Salpeter IMF.6 The factor introduces an overall uncertainty in the normalization of , due to its dependence on the unknown IMF and metallicity of these stars (or a systematic difference if the conversion is mass or redshift dependent). Note that we make the necessary conversions when comparing to other models.

2.1 Properties of Dark Matter Halos

The halo mass function measured from cosmological N-body simulations is usually expressed in the form


where is the average matter density and is the density variance smoothed over the scale . is a function determined by the particular fit (or often through analytic arguments, as in the original Press & Schechter 1974). Here we adopt the Tinker et al. (2008) mass function, corrected for the high-redshift universe by Behroozi et al. (2013b) and implemented by the online calculator HMFcalc7 (Murray et al., 2013). We note that using different functional forms does change our parameterization of the star formation efficiency slightly (see Fig. 2), but given the close agreement with high-redshift observations shown in Behroozi et al. (2013b) and the much greater uncertainties arising from the luminosity function, we do not regard the mass function as a systematic worry in our final results.

Many empirical and analytical models have been developed to characterize the growth history of dark matter halos (McBride et al., 2009; Behroozi & Silk, 2015; Correa et al., 2015). Noting the general similarity among their predicted mass accretion histories and accretion rates, we adopt the two-parameter model given by McBride et al. (2009), which fits the halo’s baryonic mass accretion rate as


where McBride et al. (2009) and Goerdt et al. (2015) found at high redshifts. The -dependence can be understood from the redshift scaling of the dark matter halo mass function (Dekel et al., 2013). The mass dependence has only been tested in simulations over a limited range of halo masses (typical of galaxies), and the power-law assumption is probably not a good assumption for very small or very large halos. However, over the limited range of halo masses relevant to galaxy formation at (where the mass function is very steep), we have found a power-law approximation to be reasonably accurate.

We note that our continuous accretion model averages over galaxy mergers, which can affect the inferred star formation efficiency (Tissera, 2000; Drory & Alvarez, 2008; Conroy & Wechsler, 2009; Pawlik et al., 2011; Vulcani et al., 2016). Goerdt et al. (2015) have shown that, over a broad range of halo masses, a clear majority of accretion occurs through the continuous mode, at least at moderate redshifts, while Behroozi & Silk (2015) have found similar results at higher redshifts. According to equation (4), the mean mass accretion rate increases much more rapidly than the Hubble expansion rate toward high redshifts, so halos grow very rapidly which likely decreases the stochasticity due to mergers. Below we will allow for some scatter in the halo mass-luminosity relationship, which can partly be due to the effects of mergers. Fortunately, its presence does not affect our qualitative results (and has only a small quantitative effect), though a model relying entirely on mergers shifts the overall abundance matching relation substantially (Visbal et al., 2015). Nevertheless, because the merger rate increases slowly with halo mass (Behroozi & Silk, 2015), mergers may affect our inferences about the brightest galaxies, and we intend to examine their effects in more detail in future work.

To determine the minimum halo mass in which galaxies can form, , we take the criterion given by Okamoto et al. (2008), which takes account of galaxy mass loss by incorporating a spatially uniform and time dependent UV background in their cosmological hydrodynamical simulations. In their model, is evaluated using the equilibrium temperature of the gas at the edge of the halo as (Noh & McQuinn, 2014)


where is the halo’s virial density, and .

2.2 Halo abundance matching

Given that the UV luminosity traces the star formation rate of galaxies (and hence accretion onto dark matter halos, according to our model), we wish to map the UV luminosity function onto the mass function of dark matter halos. We use the halo abundance matching (HAM) technique (Vale & Ostriker, 2004) to assign a unique halo mass to each UV luminosity by solving the equation


for , where and are the luminosity function of galaxies (or more precisely the intrinsic luminosity function derived from the observed one after taking the scatter in the relation into account) and the mass function of dark matter halos respectively.8 We take the latest Schechter parameterizations (see Table 1) of the observed UV luminosity functions at redshift to from Bouwens et al. (2015); Atek et al. (2015) have recently gone deeper at using strongly lensed clusters, and we will comment on the implications of their results as well. It is also worth noting that Bowler et al. (2015) find that a Schechter function might underestimate the bright end of the luminosity function at and thus might not be the best functional form at high redshifts. However, other studies find it an acceptable fit, and the bright end of the luminosity function does not dominate the integral quantities most relevant to reionization, so we simply adopt the Schechter form in this paper. Throughout these calculations, we use the observed detection limits from the UDF12 program () to determine the minimum luminosity in our abundance matching.

Research suggests that dust attenuation is non-trivial in galaxies at (Smit et al., 2012; Bouwens et al., 2014; Finkelstein et al., 2015), so we apply the method introduced by Smit et al. (2012) to perform dust corrections on the observed luminosity functions using the well-established relation between the UV continuum slope and (Bouwens et al., 2014) and the linear fit of infrared excess versus established by Meurer et al. (1999). We note, however, that the latter relation was calibrated by starburst galaxies in local universe and therefore its application to galaxies at redshift is uncertain.

( )

These are determined at a rest-frame wavelength 1600.

Table 1: Schechter parameters derived for the rest-frame UV luminosity functions at redshift , , and
Figure 1: Effects of the dust correction and stochasticity in the relation on the derived mass-magnitude relation from halo abundance matching. The solid curves show our full model, while the dashed and dot-dashed curves ignore the scatter and dust correction, respectively. A limiting magnitude of -17.7 mag is assumed in light of current detection limit.

The traditional abundance matching technique assumes that the galaxy luminosity or stellar mass is a monotonic function of halo mass, and vice versa. This ignores the possibility of intrinsic scatter in the galaxy-dark matter halo relation as well as complications from satellite galaxies (Behroozi et al., 2010; Moster et al., 2010). The latter appears to be a small effect (Wang et al., 2006; Moster et al., 2010), and at high redshifts it is difficult to observationally distinguish satellite and central galaxies thanks to the very small sizes of the halos. To incorporate the former effect, we adopt the deconvolution method described in Behroozi et al. (2010), which assumes the distribution to be a lognormal in


with being independent of the halo mass. The observed luminosity function described by the Schechter parameters can therefore be related to an unknown intrinsic luminosity function by solving the convolution


We note that the logarithmic spread in galaxy luminosity at a given halo mass, , is determined to be a constant both analytically (McBride et al., 2009; Muñoz, 2012) and observationally (More et al., 2009), so we take .

Fig. 1 shows the halo mass versus the galaxy’s observed absolute magnitude given by the abundance matching approach, as well as how the relationship depends on the dust correction and . Here the solid curves show our full model. The dash-dotted curves assume no dust correction. In particular, the dust correction is larger for larger halo masses (i.e. brighter galaxies), resulting in a net flattening of the observed luminosity function. Abundance matching without the luminosity scatter is represented by the dashed curves. As the scatter flattens the bright end of the intrinsic luminosity function, abundance matching to the observed luminosity function directly slightly underestimates the at a given luminosity.

Figure 2: Star formation efficiency as a function of halo mass and redshift. Top: Fits in the observed redshift range (–8). Data points with error bars show the uncertainties from 1000 Monte Carlo fits to the observed Schechter functions. Solid curves show the -dependent fits to the data, while the thick dashed curve shows the -independent fit to the same data. The vertical arrow indicates the detection limit of the UDF12 survey at (; this increases slightly toward lower redshifts). Below this limit, we show three potential extrapolations of the abundance matching relation (see text). Middle: Solid cures show the evolution of in our -dependent fit, assuming the extrapolation of Model III. Bottom: Variations in model results, shown at . The curves vary the input assumptions (see legend and text). The red shaded region shows inferred from Atek et al. (2015) and .

2.3 Evolution of

Fig. 2 shows the star formation efficiency, , as a function of halo mass and redshift, that results from combining equations (1) to (5). Here the uncertainties in are determined from the quoted errors on the three Schechter parameters and via 1000 Monte Carlo realizations. We note that the error bars provided represent the uncertainty at a specific halo mass estimated from Monte Carlo realizations and do not consider the correlation between data points resulting from abundance matching. We also ignore correlations between the measured Schechter parameters for simplicity, finding that they can change the inferred shape marginally but remain within the uncertainty envelope illustrated here. The impact of choosing different mass functions and halo accretion models is insignificant compared with the observationally-driven uncertainties in . However, including the dust correction raises by a factor of and slightly increases the peak halo mass, since the former is related to the amplitude of whereas the latter is determined by the slope (see Fig. 1).

Fig. 2 shows that the star formation efficiency evolves strongly with halo mass, peaking at a characteristic halo mass of and decreasing at both smaller and larger masses. At lower redshifts, similar relationships have also been identified (Conroy & Wechsler 2009; Moster et al. 2010; Kravtsov et al. 2014; Behroozi & Silk 2015). At these times, AGN feedback and feedback from massive stars, such as supernovae and stellar winds, are considered to be the dominant physical processes that suppress the star formation at the high- and low-mass ends, respectively (Mo et al. 2010; Kravtsov et al. 2014). At high redshifts, these feedback mechanisms are likely to be less efficient, because (at a fixed halo mass) galaxies are more compact at those times. However, their effects on galaxies are highly uncertain.

Additionally, there might be a non-trivial dependence on redshift, although it is difficult to quantify the dependence given the large errors, especially for smaller masses. The halo mass at which peaks evolves with redshift as for , which is broadly consistent with that found by Behroozi & Silk (2015). This empirically-derived trend might indicate evolution in the feedback mechanisms that control . However, there is no obvious reason for to decrease toward higher redshifts, as the two most important factors in decreasing the efficiency in massive halos at low redshifts (AGNs and heating at the virial shock) are probably less important at high redshifts.

Figure 3: Cumulative emissivity as a function of UV magnitude. The thick, opaque curves show the emissivities integrated down to the corresponding magnitude (-axis), given a minimum halo mass of and our fiducial -dependent model (Model II, see later text). The thin, transparent curves represent the same emissivities but assuming does not evolve with redshift.

More generally, we can examine the overall dependence of on both mass and redshift. Fig. 2 shows our fiducial, six-parameter multiple regression fit


to the abundance matching results, illustrated at –8 in the top panel. The best-fit parameters are and . We intentionally avoid higher-order terms in in order to prevent strong redshift evolution from appearing as an artifact of our fits. Fig. 2 shows the evolution in this model from –16. All the curves intersect near , where the two redshift-dependent terms cancel with each other, leaving the expression independent of redshift. This is nothing but an artifact of the functional form, which models the redshift evolution to the linear order.

For comparison, the thick dashed curve in the top and middle panels of Fig. 2 shows the best-fit result from a regression of the same order in , assuming no redshift dependence (i.e. setting , to 0). The key effects of our fitted form’s redshift dependence are to decrease and to steepen the decline toward higher masses as increases. The bright end of the luminosity function is therefore relatively uncertain at high redshift, but integrated quantities like the global emissivity, which depend primarily on faint galaxies, do not evolve as rapidly. As shown in Fig. 3, the cumulative emissivity is dominated by galaxies well below the characteristic magnitude () at all redshifts, whether or not we allow redshift evolution. However, the -dependent model predicts a much gentler decline in the total emissivity toward higher redshifts, thanks to the decrease of with . For completeness, we also evaluate the importance of the assumed linear redshift dependence (at a fixed ) in equation (9). We find that the uncertainties at induced by deviations from our assumed dependence are comparable in amplitude to the envelope on the curve shown in Fig. 2.

As star formation might be sustained in halos of mass as low as , the cumulative star formation history is sensitive to the extrapolation of with mass as well. In contrast with similar works which treat as a redshift-independent parameter (e.g. Mashian et al. 2016; Mason et al. 2015), we fit as a function of both halo mass and redshift. We then explore three possible scenarios for extrapolating toward the dominant low masses: (i) Model I assumes a constant below the mass cutoff ;9 (ii) Model II assumes a power-law extrapolation of slope 0.5 below the mass cutoff, and (iii) Model III fits a power law to the low-mass end according to the average (i.e. -independent) slope over .10 Note that in Model III, we use the average slope because the evolution of the cutoff mass and the low-mass end slope will bias the result of -dependent, power-law extrapolation. The top panel of Fig. 2 compares these prescriptions at .

Because the behavior at low luminosities is so important, we next consider recent observations from Atek et al. (2015) that extend the measured luminosity function to  mag using strong lensing from clusters. We extend our abundance matching procedure at (the only redshift for which Atek et al. 2015 have substantial data) to this limiting magnitude using the Schechter parameters determined by that group. The resulting is projected on top of our proposed extrapolations in the bottom panel of Fig. 2, which indicates that might scale as or slightly shallower in low-mass halos (Model II). Here we only use the error in the faint-end slope quoted by Atek et al. (2015) in the match.

Figure 4: Comparison of observed UV luminosity functions at and those predicted by equation (9). Here we do not include star formation missing from UV observations because of dust obscuration. 1- confidence intervals are shown by the shaded regions.

Fig. 4 compares our parametric model to the observed UV luminosity functions at . Note that when generating these model luminosity functions, we do not include star formation missing from UV observations because of dust obscuration, so these figures do not show the total star formation rate density in the Universe. The corresponding best-fit coefficients of equation (9) are and . Below, when presenting predictions for the luminosity function we follow the same procedure, but when presenting estimates of the global star formation history and reionization history we always account for dust. As expected, our model agrees very well with the available data, with the only significant difference a modest underprediction of the bright end of the luminosity function at , at least according to the Bouwens et al. (2015) measurements.

2.4 Comparison with previous work

Abundance matching is a widely-used technique at both low and high redshifts, so here we compare our results to similar studies. At lower redshifts, the declining star formation efficiency to lower-mass halos is typically attributed to stellar feedback, either through winds generated by supernovae or by radiation pressure (see Loeb & Furlanetto 2013 and references therein). The detailed effects of such feedback processes are not well-understood, but crudely they should produce a scaling , with , depending on whether the wind momentum or energy is the dominant factor. At , Behroozi et al. (2013a) find that, on average, is proportional to below the characteristic halo mass at which peaks. The mass scaling suggested by the Atek et al. (2015) data is also consistent with this range.

Figure 5: Derived evolution of the ratio of stellar mass to halo mass at using the our model of continuous star formation, compared to the SMHM relations determined at (solid triangles, Behroozi & Silk 2015; solid circles, Harikane et al. 2015) and those at (empty triangles, Behroozi & Silk 2015; empty circles, Harikane et al. 2015) and in local universe (empty squares, Behroozi et al. 2013b). Solid lines represent Model III, whereas dashed and two dotted lines represent our -independent fit, Model II, and Model I, respectively (bottom to top). Stellar masses are determined by tracing the continuous gas accretion of halos with .

Much of the work on abundance matching, especially at lower redshift, focuses on the stellar mass to halo mass (SMHM) relation. We can calculate this ratio by integrating our star formation rates from the time halos reach , assuming a continuous star formation history for each halo. We compare to other work in Fig. 5. As anticipated, the general shape of the SMHM ratio follows that of , reaching the peak at approximately at redshift . This is also qualitatively consistent with the steadily decreasing average masses of galaxies, , 11.6, and 11.4 at , 6, and 7, respectively, found by Finkelstein et al. (2015) through direct abundance matching to the observed UV luminosity function (i.e. without a dust correction). Comparison with local measurements suggests that the peak halo mass might not have evolved significantly since , whereas the peak value itself might have increased slightly with cosmic time, although with large uncertainties.

Figure 6: Comparison of predicted and observed luminosity functions at (top panel) and (bottom panel). The 1– uncertainty level is shown for Model III. The vertical dashed and dotted lines show the UDF12 observational limit and the limiting magnitude of JWST, respectively. The latter assumes a  s integration and a signal-to-noise ratio of 10. Note that the luminosity functions shown are predicted by our model of without correcting for dust attenuation, which is expected to be trivial at .

We also compare our derived evolution of the SMHM ratio with that identified by Behroozi & Silk (2015) (see also Behroozi et al. 2013b) and Harikane et al. (2015). At , we find a similar trend of increasing SMHM ratio with time at , although the increase is less rapid than the evolution inferred by Harikane et al. (2015), who used clustering observations with a modeled halo occupation distribution.

Note that the difference between the SMHM ratios originally measured by Behroozi et al. 2013b and Harikane et al. (2015) is mostly attributed to the distinct cosmologies and input models used (Behroozi & Silk 2015 includes a redshift-dependent stellar mass correction relative to observed masses in order to fit data over a wide range of redshifts simultaneously), each contributing dex of difference. The difference might also be associated with our simple continuous star formation histories (as opposed to the dark matter merger trees used by Behroozi & Silk 2015, though they argue that the merger channel is not a significant contributor at these redshifts). Accurate determination of the SMHM relation also relies on appropriate treatments of different feedback mechanisms, as demonstrated by recent simulations of dwarf galaxies (e.g. Wheeler et al. 2015) which suggest a lower SMHM ratio than Behroozi et al. (2013b) and Behroozi & Silk (2015). Nonetheless, given the qualitative similarity between the two independent methods and the large uncertainties associated with both of them, we do not pursue the explanation further. Moreover, following our modeled evolution of , we find a nearly constant mass-to-light ratio over consistent with that measured by Barone-Nugent et al. (2014) using galaxy clustering.

Recently, abundance matching and similar techniques have also been applied to the galaxy population in order to understand high- galaxy formation, particularly in light of the recent CMB measurements (Mashian et al., 2016; Mason et al., 2015; Visbal et al., 2015). Compared to these works, we take a more general and rigorous approach to a measurement of , allowing both mass and redshift dependence and allowing a more general extrapolation to the critical faint galaxy population. In contrast, Mashian et al. (2016) averaged over a range of redshifts to obtain , while Mason et al. (2015) calibrated two-parameter halo growth model at . Meanwhile, Visbal et al. (2015) derived a redshift-dependent by averaging the star formation efficiency over halos of different masses at a single redshift and then matching the normalization of the global star formation rate density to the halo growth rate. They then consider two models, one in which they fix the shape of and another in which they take a mass-independent average.

In all of these similar studies, the results are qualitatively similar to ours but differ in the details. Mason et al. (2015) find a stellar mass-halo mass relation very similar to our -independent result, though with a slightly lower normalization. Mashian et al. (2016) find a slightly smaller characteristic mass and a weaker decline at the high-mass end. In both of these cases, the models’ predictions for the bright end of the luminosity function at –10 are very similar to ours, but their redshift-independent models predict a significantly steeper decline in the overall star formation rate toward high redshift (see the next section). Visbal et al. (2015) allow a duty cycle in the star formation rate within galaxies (to represent their assumption that mergers dominate star formation) and therefore find a peak star formation efficiency at , well below our value.

Relative to these previous studies, our results demonstrate the importance of allowing for redshift evolution in and in constraining that evolution in the future. Redshift-independent models are likely oversimplified as recent observations of galaxies at have shown evidence for a non-trivial evolution of the SMHM ratio (Harikane et al., 2015; Rodríguez-Puebla et al., 2016). While our measurements show only tentative evidence for redshift evolution, given the limited range of the observations, our best-fit result has significant implications for the star formation history, as we will show in the following sections. Allowing for such evolution is also important because of the rapidly-evolving halo mass function and the expectation that the feedback mechanisms governing star formation at these epochs are themselves redshift-dependent.

Figure 7: The evolution of the star formation rate density (SFRD) and the UV luminosity density with redshift. The red solid curve shows the SFRD evolution calculated by Model III, with the shaded region representing the 68% confidence interval, integrated down to to match observations from Hubble Frontier Field data, shown by the green (dust-uncorrected, McLeod et al. 2015) and grey (dust-corrected, Bouwens et al. 2015a) circles respectively. Cases where the dust correction is not included (thin dashed curve), Model I/II is assumed (yellow/cyan curves, only the deviated parts from Model II), or -evolution is ignored (thick dashed curve) are shown for comparison. As mentioned in the text, the mass cutoff distinguishing the three models is chosen so that all of them are consistent with the 1– interval of the observed at . The pink diamonds represent the SFRD derived from the observed bright end of UV luminosity functions at and 10 (Bouwens et al. 2015b), whereas the red inverted triangles show our predictions. The black, dotted line shows the maximum-likelihood fit to given by Robertson et al. (2015), assuming a limiting luminosity (or ), whereas the cyan, dotted line represents the prediction by Model II at the same limiting magnitude. The upper limits shown illustrate the evolution of the maximum possible SFRD calculated by extrapolating our best-fit Model I (yellow) and Model III (red) down to the minimum halo mass given by Okamoto et al. (2008).

3 Predictions for High-Redshift Galaxy Observations

Because our model for the star formation efficiency uses theoretical inputs (the halo mass function and accretion rates) that can easily be extrapolated to higher redshifts, it is straightforward to test our results against higher redshift data and make predictions for future observations. Fig. 6 shows the predicted luminosity functions at and in comparison with recent measurements (Oesch et al. 2013; Oesch et al. 2014; Bouwens et al. 2015b). To date, data on the faintest galaxies disfavor Model I, but the constraints are not strong. As the primary tool to study the evolution of luminosity function of high-redshift galaxies, JWST will reach a detection limit about 2 magnitudes fainter than the current one (-17.7 mag, excluding lensing) with  s of integration (shown by the vertical dotted lines), thereby providing a much better understanding of how might evolve in low-mass halos.

Another straightforward observational implication of our model is the evolution of star formation rate density, , which we estimate based on the simple formalism given by equation (1) and our three extrapolations.11 Fig. 7 shows , and the equivalent UV luminosity density, predicted by our baseline models assuming and . The solid and dashed curves in red represent ’s under Model III with and without a dust correction (at these bright limiting luminosities, all three of our models are identical except at very high redshifts). The dust correction becomes increasingly less important between , due of the decreasing abundance of dusty massive galaxies at higher redshifts. The star formation rate density at predicted by our dust-corrected model is in good agreement with those recently determined from observations (green circles, McLeod et al. 2015; orange squares, CANDELS/GOODS/HUDF, Oesch et al. 2014; blue triangles, CLASH cluster searches, Zheng et al. 2012, Coe et al. 2013). The empty blue square shows the integrated star formation rate from Atek et al. (2015) at , which is consistent with all of our models.

There has been some controversy in the literature about the shape of the SFRD at high redshifts (e.g., Robertson et al. 2013; Oesch et al. 2015), reflected by the wide range in estimated star formation rates at higher redshifts shown by the data points in Fig. 7. Model III implies a very slow steepening of with redshift, whereas Models I and II predict a gradually flattening . There is no strong evidence for a sudden drop of the dust-corrected and the overall trend is closer to the linear relation suggested by McLeod et al. (2015), although the exact slope depends on whether or not a dust correction is included. However, forced power-law fits to our dust-corrected model over and yield and , respectively, comparable to scenarios favoring rapid evolution at (Oesch et al., 2015; Mashian et al., 2016). This apparent break is more a result of the fitting method than of any feature in our model, which evolves very smoothly in the context of continuous gas accretion.

We note that the rapidly steepening predicted by various models (Mashian et al., 2016; Mason et al., 2015) might be associated with the assumption of -independent . As shown by the black dashed curve, an fixed in time predicts substantially steeper evolution of SFRD with redshift. Interestingly, the apparent rate depends strongly on whether a dust correction is included, as that is larger for more massive halos and decreases quickly in importance at higher redshifts. Moreover, the extrapolation is very sensitive to assumptions about the behavior of low-mass galaxies. As indicated by the contrasting case of Model I (shown by the uppermost solid curve), different ways of extrapolating leave much uncertainty in the evolution at , calling for more detailed investigation by future observations.

We also show, as a contrasting case, the cosmic predicted by the maximum-likelihood fit of Robertson et al. (2015), which extrapolates the observed evolution to a limiting luminosity of (black, dotted curve), and forces smooth evolution from via a double power-law model. We compare with the prediction of Model II assuming the same limit (cyan, dotted curve). While the two approaches are remarkably similar considering the extrapolations involved, our model has significantly steeper redshift evolution than the purely empirical fit. This occurs because the low-mass halos (which become increasingly dominant at higher redshifts) are inefficient at producing stars in Model II. In the context of our approach, the Robertson et al. fit therefore implicitly assumes that low-mass galaxies become even more efficient at star formation toward higher redshift, demonstrating the importance of relating extrapolations to physical quantities in order to understand their systematic uncertainties.

Finally, we present our predictions for future deep sky surveys that will be conducted by the James Webb Space Telescope (JWST). For observations of the UV continuum, the flux density (or magnitude) can be acquired from the intensity of spatially unresolved objects. Knowing the limiting magnitude consequently allows us to estimate the integration time required to observe the majority of light from the faintest galaxies, as well as to calculate the galaxy number counts per unit area surveyed at a given redshift, which can be expressed as (Pawlik et al., 2011)


where is the comoving volume element per unit solid angle and redshift, and is the minimum observable UV luminosity.

As shown by the vertical dotted curves in Fig. 6, JWST can significantly improve the determination of the faint-end slope and evolution in faint galaxies, which helps better constrain star formation in the most primitive objects. Fig. 8 shows the predicted galaxy number counts per 10 and the fraction of the total luminosity density at a given redshift accessible to JWST for observations of  s and  s (grey and black lines, respectively). From the upper panel, Model I predicts many more observable galaxies for JWST than the other two models, in which drops rapidly with decreasing halo mass. JWST will only detect substantial numbers of galaxies at in deep integrations if Model I is correct, but it is already in tension with observations at –10 (and with constraints on reionization, as we will see below). For an integration time of  s, JWST is expected to observe roughly 10 UV-bright galaxy per 10 beyond in Model II. On the other hand, a typical second integration will allow JWST to observe more than half of the total UV continuum luminosity from star-forming galaxies at for Model II and for Model III. The apparent flattening (and even increase) of the observable fraction with redshift in Model I is due to the increasing fraction of UV luminosity from faint galaxies, as predicted by the maximally extrapolated , rather than an actual increase of galaxy population or total UV luminosity. Therefore, assuming a star formation efficiency as described by Model II or III, JWST is expected to determine the cosmic star formation rate between to better than 50% and thereby provide much stronger constraints on how faint galaxies might have affected reionization.

Figure 8: Top: Predicted galaxy number counts per for JWST assuming an integration time of seconds and requiring a signal-to-noise ratio of 10 for a detection. The three black curves show the different models for extrapolating to small halo masses. The grey curves represent cases which include photoionization suppression (corresponding to the violet curves in Fig. 10, respectively). Bottom: Fractions of the total luminosity at a given redshift that can be observed by JWST (in the presence of photoionization suppression) with an integration time of  s and  s, respetively. The minimum halo mass is evaluated at the Okamoto et al. (2008) limit.

4 Modeling Cosmic Reionization

To this point, we have focused on direct interpretations of galaxy observations. These are limited in that they cannot yet probe the faint galaxy population at , and Fig. 8 shows that the faint end at will likely remain hidden even after JWST begins observations. We next turn to combining our galaxy model with constraints on the reionization history. Specifically, the IGM ionization state is mostly determined by the emissivity of faint galaxies and is closely related to and . Thus it is helpful to compare models of reionization constructed from our inferred star formation history to various existing observational constraints.

Meanwhile, adding these explicit probes of the faint galaxy population allows us to explore sensibly the ways in which new physics (beyond those relevant to lower redshift galaxies) may affect the star formation history. In this section, we shall introduce toy models for two such effects, the suppression of faint galaxies due to photo-heating during reionization and the introduction of Population III stars. Our models will be highly simplified but will allow us to understand the trends imprinted by these physical processes.

The ionization state of intergalactic hydrogen can be measured by the globally-average ionized fraction, , whose net rate of change is a balance between ionization and recombination,


where is the overall ionizing efficiency, a product of the correction factor for single ionization of helium, , the star formation efficiency, , the escape fraction of ionizing photons, , and the mean number of ionizing photons produced per stellar baryon, . In the second term, is the coefficient for case B recombination12 at electron temperature , is the volume-averaged clumping factor of the IGM, and is the number density of hydrogen nuclei in local universe. In our calculations, we take , which is a reasonable approximation during cosmic reionization (Shull et al. 2012; Pawlik et al. 2015) and .

The value of varies with the stellar population and thus depends on the IMF and metallicity. We take for Population II stars, assuming a Salpeter IMF and . In our baseline model, we assume all halos above form stars in a similar way so that


Note that we use our estimate for in all of our calculations rather than imposing a value drawn from theoretical considerations. The exact value of is not known, so we leave it as a free parameter in our models, scaling to the value provided by Okamoto’s criterion. It is especially important in Model I, where the star formation efficiency remains high inside small halos.

4.1 Suppression of faint galaxies via photo-heating

Simulations suggest that star formation in “photo-sensitive” halos less massive than a few times will be suppressed by photo-heating as the IGM becoming reionized (Thoul & Weinberg 1996; Gnedin 2000; Finlator et al. 2011; Noh & McQuinn 2014). We therefore also consider a more sophisticated model, similar to those discussed by Alvarez et al. (2012) and Visbal et al. (2015), in which we set a cut-off mass, , to distinguish halos in which star formation was gradually suppressed as cosmic reionization proceeded. As a simple model, we assume that halos with cannot form stars inside regions that have already been ionized. For simplicity, we ignore the relative bias of these halos and the ionized regions and we do not impose any delay on the suppression effect for the same reason. Consequently, evolves following


4.2 Population III stars in minihalos

Another often-discussed complication to the early star formation history is the presence of Population III stars. Their contribution is nearly completely unconstrained at present (though see Visbal et al. 2015 for weak constraints based on methods similar to our own). We therefore simply take a toy model of these objects in order to understand their implications for the reionization process. In particular, we take for Population III stars (Bromm et al., 2001). This value assumes that Pop III stars are very massive and hence efficient ionizers. We further assume that they only formed in halos where molecular cooling is efficient (), or halos of mass between and , with an arbitrary constant efficiency that is limited by the observed CMB optical depth. In order to model the transition from Pop III to Pop II star formation, we assume that a fraction


of halos in the relevant mass range (and outside of ionized regions, where photoheating will suppress such sources) are able to form stars. We take to be the fiducial value for the time when PopIII stars started to form and vary .

Because these halos are subject to photoheating suppression, when studying Pop III stars we use the following model for the ionizing emissivity:


Note that, over the range , we let be a free parameter and do not adopt any form of extrapolation used for atomic cooling halos.


.47 {subfigure}.47 {subfigure}.47 {subfigure}.47

Figure 9: Observationally constrained parameter space for the baseline model of reionization, assuming maximally (top-left panel), moderately (top-right panel), and minimally (bottom-left panel) extrapolated (Models I, II, and III, respectively). The parameter space for the redshift-independent model is also shown for comparison (bottom-right panel). The plot is color-coded by the predicted Thomson scattering optical depth, with the thick black contours comparing to the value measured by the Planck satellite, . Yellow curves show model predictions for the logarithm of the emissivity of ionizing photons, , evaluated at redshift , which is measured to be approximately by Becker & Bolton (2013). Dot-hatched region represents the allowed parameter space for which as inferred from the spectra of high-redshift QSOs, where is the redshift at which reionization completed (i.e. ). Portions of parameter space preferred by observational constraints on the evolution of ionization state, , are shown by the regions with slash and back-slash hatch, corresponding to and , respectively (Stark et al. 2010; Treu et al. 2013; Faisst et al. 2014; Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014). White crosses indicate reasonable combinations of and later used in Fig. 10.

5 Joint Analysis of Galaxy Populations and the Reionization Observables

Our models of galaxy growth and reionization have three key free parameters: , , and . In Section 2, we constrained the first, at least over a wide mass range, through observations. We will next use ancillary measurements of the integrated Thomson optical depth, global IGM neutrality, and ionizing background at to quantify the interdependencies and degeneracies between these parameters.

For concreteness, we employ the baseline model described by equation (12). For each model of , we limit the 2D parameter space formed by and . We allow to vary between 0.01 and 1. Observations of star-forming galaxies at suggest a small average ionizing escape fraction from a few percent to smaller than 15%, which is likely to increase towards fainter galaxies (or equivalently towards higher redshifts, e.g., Iwata et al. 2009; Siana et al. 2010; Nestor et al. 2013; Mostardi et al. 2015). In recognition of the redshift dependence of , we model its variation by rescaling as given by Okamoto et al. (2008) by a factor of 0.1–10. For the convenience of discussion, the rescaled is always quoted at .

Figure 10: The IGM neutrality (left), emissivity of ionizing photons (center), and Thomson optical depth (right). The lines are the same in each panel, with the styles indicating Models I, II, III and the -independent model and the colors indicating variations in the stellar populations and feedback prescriptions (the case with both photoionization suppression and PopIII stars is shown in violet for Model II only). See the text for details on the model parameters. We compare the predictions to a number of observational constraints (see text) in each panel. Note that including photoionization suppression makes little difference in Model III and the -independent model and therefore the orange curves are only shown for Models I and II.

5.1 The CMB optical depth

We compare our models to three important measurements of the ionization history. Our first constraint is the Thomson scattering optical depth for CMB photons, which provides an integral measure of the reionization history. This is given by


where or measures the degree of helium ionization.13 The CMB optical depth is useful because it is well-measured by recent experiments (Planck Collaboration 2015, though see below). However, it is an integral constraint over both the faint galaxy population and the overall reionization history, so it still allows a great deal of freedom in models. Fig. 9 is color-coded by and quantifies the degeneracy under Model I (top-left), II (top-right), III (bottom-left) and our -independent model (bottom-right). The thick black contours are marked according to the measurement of reported by the Planck Collaboration et al. (2015) in units of (assuming gaussian statistics); dark regions are disfavored by those observations.

As expected, given a specific optical depth , and are positively related, and the relationship depends on the relative abundance of faint galaxies encoded by the extrapolation of . Under Model I, in which is maximally extrapolated and therefore faint galaxies are abundant, must grow rapidly with increasing to maintain a constant , whereas under Models II and III is much less sensitive to due to their rapidly diminishing populations of faint galaxies. This suggests that, even though and enter our models as two independent parameters, observations place strong constraints on their joint properties in the context of a particular star formation model. Inferences about the escape fraction also rely on the assumed extrapolation of . In Model I, as small as 0.1 is large enough to reproduce the observed with a moderate comparable to our fiducial model. But in Models II and III, in which declines rapidly with decreasing halo mass, is required to reach just the lower bound of Planck’s at the same (Robertson et al., 2013; Robertson et al., 2015).

If we assume no redshift dependence and take our best fit to the aggregate data from –8, it is very difficult to match the observations without appealing to , which would in practice require changing the stellar populations to be much more efficient ionizers. This is because our -independent fit maintains at  M, far above the characteristic mass at very high redshifts. The shift of toward lower masses in our fiducial model significantly increases the contribution of high- galaxies to the total optical depth.

We emphasize that these conclusions are all in the context of our simple model. For example, the high values of we require (in comparison to measurements at lower redshifts) can be ameliorated by introducing a new population of ionizing sources at very high redshifts, which cannot be modeled properly by our extrapolations.

Fig. 10 illustrates the difficulty of matching the optical depth in our fiducial model. The right-hand panel shows the integrated optical depth, assuming and 0.2 for our three -dependent models and for our -independent model (marked by the white crosses in Fig. 9). Only Model I can easily reach the best-fit value from Planck, though all are within the lower bound. We note that while the dataset of Planck suggests a lower bound (1) of the CMB optical depth as low as (Planck Collaboration et al., 2015), the dataset itself indicates a higher optical depth (i.e., an earlier reionization). The difference between these likely roughly captures the level of systematic uncertainty in CMB measurements (Planck Collaboration et al., 2015). As shown by Fig. 10, the higher optical depth can be reached by assuming more efficient star formation in low-mass halos (close to our Model I) or allowing a higher escape fraction in our fiducial model. While the former is permitted by the existing observational constraints, we will see that the latter poses some difficulties.

5.2 Other Reionization Constraints

We next ask whether other existing observational constraints on further limit the allowed parameter space. These measures of the instantaneous neutral fraction provide different information than , because they help us distinguish the progress of different models across cosmic time. This also implies that they require modeling of the entire ionization history, which is provided by our galaxy evolution model.

The hatched regions stretching from the lower left corner in the panels of Fig. 9 show parameter space permitted by various indirect measures of the neutral fraction at particular redshifts. We show the constraints on themselves as symbols in the left panel of Fig. 10. First, the dot-hatched region in Fig. 9 shows models where reionization completes , as inferred from the spectra of high-redshift QSOs (Fan et al., 2006; Bolton et al., 2011; McGreer et al., 2015); these are also shown by the solid diamonds and open squares in Fig. 10. The regions filled by slash and back-slash hatch mark model where the global ionized fraction are and at and , respectively, in agreement with the measurements of declining Lyman alpha emission fraction of high-redshift galaxies (Stark et al. 2010; Treu et al. 2013; Faisst et al. 2014; Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014; see also the solid pentagons in Fig. 10). The limits on the IGM neutrality inferred from the measurements of quasar near zones (open hexagon, Mortlock et al. 2011) and the GRB damping wings (open circle, Chornock et al. 2013) are also shown in the left panel of Fig. 10. We note that we regard all these measurements as model-dependent thanks to the difficulty in evaluating the systematic modeling uncertainties from which constraints on the neutral fraction are drawn.

The parameter space favored by these measurements roughly traces the to interval of , suggesting a broad agreement between the ionizing efficiency required by the observed history of reionization and that actually supplied by star-forming galaxies. However, there is some tension between the baseline model and these observables, as shown by the left panel of Fig. 10. In particular, Model I is in marginal disagreement with some of these measures, because it relies on faint galaxies that evolve slowly throughout this epoch. The others, with star formation centered on more massive halos, evolve more rapidly and satisfy all the existing constraints.

5.3 The Ionizing Emissivity

Finally, we compare our models to the Lyman-limit photon emissivity , which can be measured from the Lyman- forest at (Kuhlen & Faucher-Giguère 2012, KF12; Becker & Bolton 2013, BB13). In our model, this equals


The total emissivity is useful as a constraint on the integrated population of galaxies at the tail end of reionization, though of course it is an imperfect measure because of potential redshift evolution. The yellow lines in Fig. 9 show the inferred emissivity at in our models. These should be compared to the measurement from BB13 at , .

While it is possible to match both and constraints on without undue difficulty, further matching the ionizing emissivity at , right after the completion of reionization, is challenging, particularly in Models II and III. This tension can also be seen in the middle panel of Fig. 10, which shows the evolution of the ionizing emissivity. Here we show the BB13 measurement as well as that from KF12. The difference results from BB13’s choice of a lower IGM temperature and more accurate treatments of the optical depth and cosmological radiative transfer effects. Therefore, we adopt BB13 as our fiducial limit on the ionizing emissivity at . Models II and III, which rely on massive galaxies to reionize the Universe, have a rapid increase in the ionizing emissivity that only marginally agrees with the observed values at , if they are calibrated to complete reionization at . This demonstrates a tension between the emissivity measurements and : increasing improves agreement with the latter but reduces agreement with the former. This is one motivation to explore additional physics in the following subsections. For comparison, we also present the emissivity derived with the SFRD given by Robertson et al. (2015), which provides a slightly better match to the observations while being largely consistent with our fiducial model.

Note, however, that the low required by Model I avoids over-predicting the post-reionization ionizing background. In particular, Model I fits nicely with recent simulation results that suggest a low time-averaged with negligible mass and redshift dependence (Ma et al., 2015). If Model I or a very gradual extrapolation of is true, then , together with other observational constraints, suggests that might not exceed 0.1. In that case, the minimum halo mass could not differ from Okamoto’s criterion by more than  dex.

6 Implications of Additional Source Physics

Next, we discuss variations from our baseline model. The effects we will explore here are highly uncertain, from the perspectives of both modeling and observations. Thus we simply attempt to provide a qualitative understanding of their relative importance by contrasting them with our baseline model.

6.1 Photoionization Suppression

We begin by considering the effects of the suppression of galaxy formation from photoheating during reionization, as described by equation (13) in Section 4.1. In the following discussion, we assume the critical halo mass below which star formation is “photo-suppressible” to be , consistent with the mass scale for efficient photoionization feedback (Gnedin 2000; Finlator et al. 2011), while all other parameters take the same values as in the fiducial model, unless otherwise stated.

To see how photoionization heating might influence the reionization history, we compare our revised models to the baseline model given a specific combinations of and , 0.2, 0.2 for Model I, II, III, respectively (see the white crosses in Fig. 9), so that the predicted histories of reionization are in general agreement with the observational constraints, in particular the optical depth and the redshift of completion.

The orange curves in Fig. 10 compare Models I and II with photo-heating suppression to the fiducial model discussed previously (shown in black). Including photoionization suppression clearly diminishes the contribution to the total ionizing emissivity from faint galaxies, resulting in a delayed reionization. This effect is particularly significant for Model I (dot-dashed curves), which includes a large population of photo-suppressible galaxies hosted by halos below , whereas in other models the effect of photoionization is almost indiscernible. In these models, the low-mass galaxies already have only a very small contribution to the ionizing emissivity, so their suppression has very little effect on the overall emissivity. Thus the tension between observed and predicted post-reionization emissivities under Model III is not resolved by simply including the feedback from photoionization heating.

This is one of the most important conclusions from our analysis: photoheating can only be a significant feedback mechanism if the shape of differs qualitatively from that at low redshifts, with much more efficient star formation in faint galaxies. While we cannot rule out such a scenario, it is another indication – along with the high required values – that galaxies may undergo substantial evolution at .

Figure 11: Optical depth to electron scattering as a function of Pop III star formation efficiency of molecular cooling halos, in the presence of photoionization suppression and a cessation of Pop III star formation at . Black, blue, and cyan curves (with decreasing thickness) take , respectively.

6.2 Population III Stars

Next, we present a crude analysis of the contribution from massive Pop III stars based on the simple model described in equation (14) and Section 4. We refer readers to Visbal et al. (2015) for a more careful analysis of the implications for Pop III stars in light of the recent measurements from Planck. Here, we only consider the case where , consistent with both models and numerical simulations of Pop III star formation regulated by a variety of feedback mechanisms, which suggest Pop III stars could have formed as early as and a small amount of them might still form until the end of reionization (Scannapieco et al., 2003; Furlanetto & Loeb, 2005; Trenti et al., 2009; Johnson et al., 2013).

To set , we consider in Fig. 11 the total contribution of this population to , assuming a fixed escape fraction for Pop III stars since effects like supernova blowout and weaker dust attenuation could allow more photons to escape these fragile halos. When , all models reach or exceed Planck’s upper limit of . We use this calculation to choose as a plausible star formation efficiency in the minihalo population. It is also consistent with the results of Visbal et al. (2015), which suggest the formation efficiency of Pop III stars is tightly constrained by Planck data.

Now taking and , Fig. 10 contrasts our fiducial model with one in which Pop III stars are added following our simple model with the same as more massive halos (magenta curves, only shown for Model II). The formation of Pop III stars effectively enhances the integrated ionizing efficiency, so that the reionization history develops a long tail toward higher redshifts. In our simple model, it does so in such a way that the additional contribution is largely independent of the emissivity at the end of reionization, which provides one route toward simultaneously matching the ionizing emissivity at and . However, the optical depth measurements still require the efficiency of star formation in these halos to be relatively low.

Case [] Sup PopIII
1 0.5 0.01 No No
2 0.5 0.1 Yes
3 0.3 0.05 No No
4 0.3 0.1 Yes

This is evaluated at .

Table 2: Parameter choices for the sample cases shown in Fig. 12, which include a step-wise evolution of .
Figure 12: The upper and lower left panels are identical to those in Fig. 10, except using the variations of Model II listed in Table 2. The lower right panel shows the time evolution of the escape fraction, , averaged over halo mass. The thin dashed curves are taken from Fig. 10, assuming both photoionization suppression and the formation of PopIII stars. Allowing to evolve with halo mass (i.e. cosmic time) effectively reduces the post-reionization UV background.

6.3 Evolution in

As mentioned earlier, we find that the preferred parameter spaces for Models II and III only marginally agree with the UV ionizing background measured at . This happens because the two cases are dominated by star formation in relatively massive halos, whose abundance is increasing very rapidly (whether or not photo-heating suppression is included). The post-reionization ionizing background is always dominated by these massive halos and is sensitive to the escape fraction assigned to those massive halos. Thus making the escape fraction either a function of halo mass or of redshift can potentially reconcile the tension.

To this end we consider a simple method that was used to match to the higher measured by WMAP (see Alvarez et al. 2012), which assigns different to halos according their masses. Effectively, this allows to evolve with time as well, because of evolution in the halo mass function. Fig. 12 shows a few sample cases based on Model II (and specified in Table 2) that assume a stepwise evolution of with halo mass. For the sake of simplicity, we assume the critical mass for evolution to be the same as the critical mass scale, , for photo-suppressible halos and assign two distinct escape fractions, and , to halos with mass below (the escape fraction of molecular halos is fixed at 0.5) and above , respectively. All of these choices reduce the emissivity at so as to be consistent with observations while remaining consistent with all or nearly all of our other constraints. They also increase the optical depth to electron scattering so that the models agree better with the best-fit result from Planck.

The curves in the lower right panel of Fig. 12 show the time evolution the mean escape fraction, which is averaged by the total number of ionizing photons produced by star formation. By assumption, there will be no star formation in photo-sensitive halos after reionization ends. Consequently, the low assigned to photo-insensitive halos yields a moderate emissivity at , independent of the of photo-sensitive halos. On the contrary, the much higher quickly elevates the mean escape fraction to the level required to reionize the IGM by ( for Model II as indicated by the horizontal dotted line, see also Fig. 9) as we trace back from the end of reionization, when photo-sensitive halos begin to contribute.

7 Summary

We have presented a model of the star formation efficiency of galaxies during reionization developed from the abundance matching technique and the assumption that galactic growth is driven by continuous mass accretion of dark matter halos. A variety of existing observations allow us to constrain the parameters of our models of cosmic star formation and reionization, shedding light onto the ionizing efficiency of known galaxy populations and their influence on reionization.

In contrast to many previous works, which simply treat the star formation efficiency as a constant, we are able to more accurately model its evolution with halo mass and redshift by applying the halo abundance matching technique to the dust-corrected, UV luminosity functions at –8. We find that evolves strongly with halo mass at , in qualitative agreement with its evolution at lower redshifts, and peaks at a characteristic halo mass . We model the evolution of in halos more massive than with a third order polynomial in as given by equation (9) and investigate different extrapolations to lower masses. The mass to light ratio derived from our models is in good agreement with observations. Because the luminosity function in this period is so steep, the behavior of low-mass halos is particularly important. Recent data from the Hubble Frontier Field (Atek et al., 2015) suggests a rapid decline in the star formation efficiency below the peak, with , albeit with large uncertainties. We therefore allow for a range of extrapolations to small halo masses.

By empirically extrapolating our model of to both smaller masses and earlier times, we predict the UV luminosity function and the star formation rate density at –10 in reasonable agreement with current observations. Our fiducial power-law extrapolations of in low-mass halos suggests a roughly linear or even slightly flattening in the log–log space between , with little evidence of rapid evolution at . We note, however, that the evolution pattern can be complicated by the dust correction, the assumed extrapolation of , and the underlying model of galactic growth. We also provide forecasts for future JWST observations based on our star formation histories. With a  s observation, JWST is expected to measure the faint-end slope of UV luminosity function to . For a power-law extrapolation of with slope 0.5, we predict that JWST may detect 10 UV-bright galaxy per 10 beyond and more than half of the total UV luminosity emitted at .

The qualitative behavior of , with a peak at near  M, is similar to analogous studies at lower redshifts. There, the decline toward higher masses is attributed to a combination of gravitational shock heating and feedback from active galactic nuclei (AGNs), but these effects have not been extensively studied at high redshifts. Any explanation of this trend at high masses will have to account for both the clear (though gentle) redshift dependence of the peak, which moves to smaller masses at higher redshifts, and the decreasing star formation efficiency at fixed mass as redshift increases.

At small halo masses, the star-formation efficiency is usually assumed to be set by stellar feedback, either through supernovae or winds and radiation pressure. Those processes certainly operate at these high redshifts, and the shape inferred from the Atek et al. (2015) data is consistent with expectations from simple models. This picture is also consistent with the increasing star formation efficiency in faint halos toward higher redshifts, as (at a given mass) the binding energy of halos increases with redshift.

We emphasize that our model is intentionally simple, ignoring a number of complications that must be addressed in the future. We associate star formation with halo accretion and ignore stochasticity in that relation (including mergers and a duty cycle for star formation). We also assume that the stellar populations remain constant within halos (including the IMF and metallicity). Evolution in any of these parameters could explain some or all of the trends we see: for example, if massive galaxies have significantly higher metallicities or less top-heavy IMFs than small galaxies, our model would erroneously assign them a smaller star formation efficiency. But the trend would have to be quite strong to account for all of the effect we see.

Modeling as a function of both halo mass and redshift, we find a slower evolution of the cosmic star formation rate density compared to that predicted by models that fix the star formation efficiency as a function of redshift (see also Mashian et al. 2016; Mason et al. 2015). Consequently, it is easier to explain the observed electron scattering optical depth using the star formation histories predicted by our redshift-dependent models, especially if one desires a small escape fraction comparable to that observed at lower redshifts. Our baseline models show that with a single average of and a minimum halo mass of , known galaxy populations are able to ionize the IGM by and reproduce the optical depth recently observed by the Planck satellite. However, the ionizing background at might be overestimated, especially when star formation in massive halos dominates the emissivity (Models II and III). We consequently explore the possibilities of photoionization heating, massive Pop III stars, and evolution in , in addition to our baseline model. We find that photoheating can only be a significant feedback mechanism when the shape of differs qualitatively from that at low redshifts, allowing faint galaxies to contribute a substantial fraction of ionizing photons. Our crude analysis of Pop III stars demonstrates that their presence effectively enhances the ionizing efficiency while having little impact on the post-reionization emissivity, which causes a tail of ionization to high redshifts and improves the agreement with the electron scattering optical depth. Meanwhile, the recent measurement from Planck places an upper limit on the efficiency at which those earliest stars could form (see also Visbal et al. 2015). Finally, the tension between our predicted emissivity of ionizing photons and that estimated from observations might be resolved by allowing the escape fraction to evolve with halo mass or redshift (see also e.g., Alvarez et al. 2012). Assuming that the escape fraction is typically higher in faint galaxies, we suggest that bright galaxies with an absolute escape fraction less than 10% can reproduce the observed reionization history without over-predicting the emissivity or relying on efficient star formation in low-mass halos.

Finally, we must stress that the models of presented in this paper demonstrate that the known galaxy population can self-consistently reproduce the histories of both galactic growth and cosmic reionization using a simple prescription for the halo assembly history. In the future, improved observations and numerical simulations with careful treatments of feedback will improve our understanding of the star formation efficiency trends with halo mass and redshift.


The authors would like to thank the anonymous referee, as well as Richard Ellis, Phil Hopkins, Yu Lu, Jordan Mirocha, and Brant Robertson for extremely helpful comments on the manuscript. GS also thanks Jamie Bock for his kind support and Jianfei Shen for discussions of the statistical inference. This research was completed as part of the University of California Cosmic Dawn Initiative. We acknowledge support from the University of California Office of the President Multicampus Research Programs and Initiatives through award MR-15-328388. SRF was partially supported by NASA grant NNX15AK80G, administered through the ATP program, and by a Simons Fellowship in Theoretical Physics. SRF also thanks the Observatories of the Carnegie Institute of Washington for hospitality while much of this work was completed.


  1. footnotemark:
  2. pagerange: LABEL:000LABEL:000
  3. pubyear: 0000
  4. The Planck measurement actually relies on a measurement of CMB lensing, as the estimate from the primary CMB anisotropies alone is somewhat higher (Planck Collaboration et al., 2015). Recently, Addison et al. (2015) have questioned the self-consistency of the Planck parameter extraction. We will therefore comment on the implications of a higher optical depth as well.
  5. Note that is defined in an “instantaneous” sense, proportional to the current accretion rate and depending on the current halo mass and redshift. As a halo grows, will therefore change, and our results require averaging over halo growth histories in order to calculate a halo’s net star formation efficiency.
  6. The adopted is a compromise value assuming an evolving , smaller than the value from Kennicutt (1998) which assumes solar metallicity (see Madau & Dickinson 2014 for details).
  8. is more commonly expressed in the magnitude form , where , , and are the normalization factor, the characteristic magnitude and the faint-end slope, respectively.
  9. The value is chosen here so that our extrapolations yield reasonable star formation rate densities at in conformity with observations, see Section 3 and Fig. 7 for details.
  10. This mass dependence is considerably steeper than expected from the simplest theoretical models of supernova wind feedback (e.g., Davé et al. 2011; Dayal et al. 2014). We emphasize again that the slope is not well-constrained by observations.
  11. Note that the evolution of is given by equation (9) only at . Extrapolation is required below this mass threshold.
  12. In the early universe when ionizations were distributed fairly uniformly throughout the IGM, ionizing photons regenerated by recombinations to the ground state would not be re-absorbed locally in dense clumps but ionize other IGM atoms. Thus case B recombination is reasonable, though toward the end of the process recombinations in dense systems become important (Furlanetto & Oh, 2005). However, this is largely degenerate with uncertainty in the evolution of .
  13. Here we assume helium is completely reionized at , in agreement with current observations (Fechner et al. 2006; Shull et al. 2010; Syphers & Shull 2014).


  1. Addison G. E., Huang Y., Watts D. J., Bennett C. L., Halpern M., Hinshaw G., Weiland J. L., 2015, ArXiv e-prints
  2. Alvarez M. A., Finlator K., Trenti M., 2012, ApJL, 759, L38
  3. Atek H., Richard J., Jauzac M., Kneib J.-P., Natarajan P., Limousin M., Schaerer D., Jullo E., Ebeling H., Egami E., Clement B., 2015, ApJ, 814, 69
  4. Barkana R., Loeb A., 2001, Phys. Rep., 349, 125
  5. Barone-Nugent R. L., Trenti M., Wyithe J. S. B., Bouwens R. J., Oesch P. A., Illingworth G. D., Carollo C. M., Su J., Stiavelli M., Labbe I., van Dokkum P. G., 2014, ApJ, 793, 17
  6. Becker G. D., Bolton J. S., 2013, MNRAS, 436, 1023
  7. Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  8. Behroozi P. S., Silk J., 2015, ApJ, 799, 32
  9. Behroozi P. S., Wechsler R. H., Conroy C., 2013a, ApJL, 762, L31
  10. Behroozi P. S., Wechsler R. H., Conroy C., 2013b, ApJ, 770, 57
  11. Bolton J. S., Haehnelt M. G., 2007, MNRAS, 382, 325
  12. Bolton J. S., Haehnelt M. G., Warren S. J., Hewett P. C., Mortlock D. J., Venemans B. P., McMahon R. G., Simpson C., 2011, MNRAS, 416, L70
  13. Bouwens R. J., Illingworth G. D., Oesch P. A., Caruana J., Holwerda B., Smit R., Wilkins S., 2015, ApJ, 811, 140
  14. Bouwens R. J., Illingworth G. D., Oesch P. A., Labbé I., van Dokkum P. G., Trenti M., Franx M., Smit R., Gonzalez V., Magee D., 2014, ApJ, 793, 115
  15. Bouwens R. J., Illingworth G. D., Oesch P. A., Trenti M., Labbé I., Bradley L., Carollo M., van Dokkum P. G., Gonzalez V., Holwerda B., Franx M., Spitler L., Smit R., Magee D., 2015, ApJ, 803, 34
  16. Bouwens R. J., Illingworth G. D., Oesch P. A., Trenti M., Labbé I., Franx M., Stiavelli M., Carollo C. M., van Dokkum P., Magee D., 2012, ApJL, 752, L5
  17. Bouwens R. J., Oesch P. A., Labbe I., Illingworth G. D., Fazio G. G., Coe D., Holwerda B., Smit R., Stefanon M., van Dokkum P. G., Trenti M., Ashby M. L. N., Huang J.-S., Spitler L., Straatman C., Bradley L., Magee D., 2015, ArXiv e-prints
  18. Bowler R. A. A., Dunlop J. S., McLure R. J., McCracken H. J., Milvang-Jensen B., Furusawa H., Taniguchi Y., Le Fèvre O., Fynbo J. P. U., Jarvis M. J., Häußler B., 2015, MNRAS, 452, 1817
  19. Bromm V., Kudritzki R. P., Loeb A., 2001, ApJ, 552, 464
  20. Bromm V., Yoshida N., 2011, ARA&A, 49, 373
  21. Calverley A. P., Becker G. D., Haehnelt M. G., Bolton J. S., 2011, MNRAS, 412, 2543
  22. Chornock R., Berger E., Fox D. B., Lunnan R., Drout M. R., Fong W.-f., Laskar T., Roth K. C., 2013, ApJ, 774, 26
  23. Coe et al. 2013, ApJ, 762, 32
  24. Conroy C., Wechsler R. H., 2009, ApJ, 696, 620
  25. Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015, MNRAS, 450, 1514
  26. Davé R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 415, 11
  27. Dayal P., Ferrara A., Dunlop J. S., Pacucci F., 2014, MNRAS, 445, 2545
  28. Dekel A., Zolotov A., Tweed D., Cacciato M., Ceverino D., Primack J. R., 2013, MNRAS, 435, 999
  29. Drory N., Alvarez M., 2008, ApJ, 680, 41
  30. Ellis R. S., McLure R. J., Dunlop J. S., Robertson B. E., Ono Y., Schenker M. A., Koekemoer A., Bowler R. A. A., Ouchi M., Rogers A. B., Curtis-Lake E., Schneider E., Charlot S., Stark D. P., Furlanetto S. R., Cirasuolo M., 2013, ApJL, 763, L7
  31. Faisst A. L., Capak P., Carollo C. M., Scarlata C., Scoville N., 2014, ApJ, 788, 87
  32. Fan X., Strauss M. A., Becker R. H., White R. L., Gunn J. E., Knapp G. R., Richards G. T., Schneider D. P., Brinkmann J., Fukugita M., 2006, AJ, 132, 117
  33. Fechner C., Reimers D., Kriss G. A., Baade R., Blair W. P., Giroux M. L., Green R. F., Moos H. W., Morton D. C., Scott J. E., Shull J. M., Simcoe R., Songaila A., Zheng W., 2006, A&A, 455, 91
  34. Ferrara A., Loeb A., 2013, MNRAS, 431, 2826
  35. Finkelstein S. L., Song M., Behroozi P., Somerville R. S., Papovich C., Milosavljević M., Dekel A., Narayanan D., Ashby M. L. N., Cooray A., Fazio G. G., Ferguson H. C., Koekemoer A. M., Salmon B., Willner S. P., 2015, ApJ, 814, 95
  36. Finkelstein et al. 2015, ApJ, 810, 71
  37. Finlator K., Davé R., Özel F., 2011, ApJ, 743, 169
  38. Furlanetto S. R., Loeb A., 2005, ApJ, 634, 1
  39. Furlanetto S. R., Oh S. P., 2005, MNRAS, 363, 1031