SCORCH. I. THE GALAXY–HALO CONNECTION IN THE FIRST BILLION YEARS

Cosmic reionization is a frontier topic in cosmology with plenty of scientific richness for theoretical and observational explorations. The epoch of reionization (EoR) is uniquely marked by the emergence of the first luminous sources: stars, galaxies, and quasars in the first billion years. It is possibly the only time that luminous sources directly and dramatically alter the state of the entire universe, converting the cold and neutral intergalactic medium (IGM) into a warm and highly ionized one. Studying the EoR will reveal how the first generation of stars, galaxies, and black holes formed and evolved. It can also provide complementary constraints on cosmological parameters similar to studies of the cosmic microwave background (CMB). See Loeb & Furlanetto (2013) for a recent review.

Current observations suggest that reionization was already in significant progress by z ∼ 9 and must have ended by z ∼ 6. The photoionization and photoheating of hydrogen leaves many observable imprints. Neutral hydrogen can be detected through 21 cm radiation from radio observations (e.g., Bowman & Rogers 2010; Parsons et al. 2014), in absorption as the Lyα forest (e.g., Fan et al. 2002, 2006), and also through Lyα scattering of high-redshift galaxies (e.g., Kashikawa et al. 2006; Ouchi et al. 2009). Free electrons will scatter CMB photons and produce patchy Thomson scattering (e.g., Dvorkin & Smith 2009; Su et al. 2011; Natarajan et al. 2013), kinetic Sunyaev–Zel'dovich temperature anisotropy on arcminute scales (e.g., Sunyaev & Zeldovich 1970; Zahn et al. 2012; Sievers et al. 2013), and polarization anisotropy at low multipoles (e.g., Kogut et al. 2003; Hinshaw et al. 2013; Planck Collaboration et al. 2015). Heating of the IGM affects the thermal broadening of the Lyα forest (e.g., Theuns et al. 2002; Cen et al. 2009; Lidz & Malloy 2014). Ongoing and upcoming observations will provide multi-wavelength constraints on the uncertain timing and duration of the EoR.

Current theory suggest that large-scale, overdense regions near radiation sources are generally reionized earlier than large-scale, underdense regions far from sources. Three main approaches are used to model the EoR (e.g., Trac & Gnedin 2011). Cosmological simulations combining N-body, hydro, and radiative transfer (RT) algorithms are the most accurate, but also most expensive approach to solve the coupled evolution of the dark matter, baryons, and radiation (e.g., Gnedin 2000; Trac et al. 2008; Gnedin 2014; Norman et al. 2015). RT post-processed on saved N-body or hydro simulated data can be more cost effective while capturing important features of nonlinear structure formation (e.g., Ciardi et al. 2003; Iliev et al. 2006; McQuinn et al. 2007; Finlator et al. 2009; Aubert & Teyssier 2010). Semi-analytical/numerical methods provide an approximate and efficient approach to explore the vast parameter space of reionization (e.g., Furlanetto et al. 2004; Zahn et al. 2007; Santos et al. 2010; Mesinger et al. 2011; Battaglia et al. 2013b). Recently, there is renewed emphasis on using radiation-hydrodynamic simulations to model the complex astrophysics of sources and sinks and using semi-numerical methods to make theoretical predictions and mock observations. These two approaches in tandem provide our best option in making forward progress in synergy with observations.

SCORCH (Simulations and Constructions of the Reionization of Cosmic Hydrogen) is a new project to study the EoR and provide useful theoretical tools and predictions to facilitate more accurate and efficient comparison between observations and theory. To build a robust theoretical framework, we need to further investigate how the distribution and properties of radiation sources and sinks affect the complex reionization process. The RadHydro code (Trac & Pen 2004; Trac & Cen 2007) will be used to produce N-body + hydro + RT simulations with subgrid physics modeled using the latest observations and simulations. Mock observations will be constructed by mapping higher-resolution, smaller-volume radiation-hydrodynamic simulations onto lower-resolution, larger-volume N-body simulations (e.g., Battaglia et al. 2013a, 2013b; Natarajan et al. 2013; La Plante et al. 2014).

Currently, there are ≈1500 galaxy candidates observed in the redshift range 6 ≲ z ≲ 10 from Hubble Space Telescope observations (e.g., Grogin et al. 2011; Koekemoer et al. 2011; Windhorst et al. 2011; Ellis et al. 2013; Illingworth et al. 2013; Schmidt et al. 2014). The galaxy luminosity function is well fit by a Schechter function and the redshift dependence of the parameters have been quantified (e.g., Bouwens et al. 2014; Finkelstein et al. 2014; Oesch et al. 2014; Bouwens et al. 2015a). However, it is unclear what the physical origin for the evolution is and how is it connected to the growth of structure, particularly that of dark matter halos. The abundance of dark matter halos are accurately quantified using N-body simulations (e.g., Jenkins et al. 2001; Tinker et al. 2008), but previous work related to the EoR (e.g., Lukić et al. 2007; Reed et al. 2007) have not extensively studied atomic cooling halos ($T\gtrsim {10}^{4}$ K, $M\gtrsim {10}^{8}\ {M}_{\odot }$) capable of hosting high-redshift galaxies. Similarly, the growth of dark matter halos have also been measured using N-body simulations (e.g., Wechsler et al. 2002; Fakhouri & Ma 2008; Correa et al. 2015b), but generally for higher mass and at lower redshift.

In Paper I, we quantify the connection between observed high-redshift galaxies and simulated dark matter halos in order to better understand the abundance and evolution of the primary source of ionizing radiation. Section 2 describes the N-body code and halo finder used to simulate and locate dark matter halos. Sections 3 and 4 describe how the abundance of dark matter halos as a function of mass M, accretion rate $\dot{M},$ and redshift z are quantified. In Sections 5 and 6, we use the abundance matching technique to estimate the luminosity–mass relation and the luminosity–mass–accretion-rate relation and infer the star formation efficiency. In Sections 7 and 8, we make predictions for the high-redshift galaxy luminosity function and forecast galaxy counts for the upcoming James Webb Space Telescope (JWST). Section 9 discusses the implications of our results on cosmic reionization and Section 10 concludes with a summary of major results. We adopt the concordance cosmological parameters: ${{\rm{\Omega }}}_{{\rm{m}}}=0.27,$ ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73,$ ${{\rm{\Omega }}}_{{\rm{b}}}=0.045,$ h = 0.7, σ₈ = 0.8, and ${n}_{{\rm{s}}}=0.96.$

N-body simulations are run using a new particle–particle–particle-mesh code (P³M; e.g., Hockney & Eastwood 1981) containing updated algorithms. The gravitational potential is decomposed into short and long-range components using a Gaussian kernel following Bagla (2002). The long-range potential is computed using a particle-mesh algorithm where Poisson's equation is efficiently solved using Fast Fourier Transforms. Particles are assigned to a mesh using the triangular-shaped-cells scheme to reduce mesh anisotropy. The short-range force is computed for particle–particle interactions using direct summation, which is only moderately expensive at high redshifts when the large-scale structure is not highly clustered. Adaptive time stepping is performed using the kick-drift-kick scheme from Springel (2005).

Dark matter halos are found using a new hybrid finder that has two major components. First, a friends-of-friends (FoF; e.g., Jenkins et al. 2001) algorithm is used to locate the peaks of halo candidates. A linking length b = 0.08 times the mean interparticle separation is chosen to pick out only the inner regions and prevent over-merging (e.g., Cohn & White 2008). For each peak, the particle with the minimum value of the density-weighted gravitational potential $\rho \phi$ is chosen as the center.

Second, a spherical overdensity (SO; e.g., Lacey & Cole 1994) algorithm is used to measure halo properties. A halo with mass ${M}_{{\rm{\Delta }}}$ within radius ${R}_{{\rm{\Delta }}}$ has an average density that is factor of ${{\rm{\Delta }}}_{{\rm{halo}}}$ times the average matter density ${\bar{\rho }}_{{\rm{m}}}$ such that

$$\begin{eqnarray}&&{M}_{{\rm{\Delta }}}=\displaystyle \frac{4}{3}\pi ({{\rm{\Delta }}}_{{\rm{halo}}}{\bar{\rho }}_{{\rm{m}}}){R}_{{\rm{\Delta }}}^{3}.\end{eqnarray} \tag{ 1 }$$

Throughout the paper, M and R refer to M₂₀₀ and R₂₀₀, respectively. Halos are allowed to overlap, but if two candidates have centers within the larger halo's radius, one candidate is chosen as the primary halo and the other is considered a subhalo. The primary halo is the one that satisfies more of the following criteria: larger masses and more negative potentials within R₂₀₀ and R₅₀₀. Only halos with at least 400 particles ($1/\sqrt{400}=0.05$) are counted to avoid resolution effects and ensure completeness near the low-mass end of each simulation.

The halo finder is run "on the fly" every 20 million cosmic years. This time interval is comparable to relevant timescales such as the halo dynamical time and the mass accretion timescale at high redshifts. In constructing the halo merger trees, every particle is tracked rather than just the most-bound particle in identified halos. We use both forward and backward linking to carry out a two-step process in progenitor-descendent matching. First, for a given progenitor of mass M₂ at a given time t₂, its descendent of mass M₁ from a previous time t₁ is the one that contributes the largest weight, which is defined to be the sum of the gravitational binding energy from particles in the progenitor that once belonged to the descendent. Second, the selected descendent may contribute particles to other halos, but the weighted contribution must be less than what it gave to the main progenitor. The mass accretion rate is calculated as

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{M}_{2}-{M}_{1}}{{t}_{2}-{t}_{1}}.\end{eqnarray} \tag{ 2 }$$

Since it is more difficult to robustly measure accretion rate than mass, only halos with at least 2500 particles ($1/\sqrt{2500}\;=\;0.02$) are used.

Table 1 lists the N-body simulations used to quantify the abundance of halos as a function of mass, accretion rate, and redshift. Box sizes in the range $10\leqslant L/({\rm{Mpc}}\;{h}^{-1})\leqslant 400$ are chosen to focus on atomic cooling halos and two realizations of each box size are run to reduce sample variance. The nine largest box-size simulations are necessary to reach convergence at the lowest M of interest, while the two smallest box-size simulations are for convergence at the lowest $\dot{M}.$ Each simulation contains 2048³ dark matter particles and has a particle mass resolution $8.72\times {10}^{6}{(L/100)}^{3}\ {M}_{\odot }{h}^{-1}.$ The gravitational softening length is set to 1/16th of the mean interparticle spacing. Initial conditions are generated using second-order Lagrangian perturbation theory (2LPT; Scoccimarro 1998) with a linear transfer function from CAMB (Lewis & Bridle 2002). The starting redshift ${z}_{{\rm{init}}}$ is chosen such that the perturbations, ${\delta }_{{\rm{rms}}}\sim 0.01,$ are still linear.

Table 1.  N-body Simulation Parameters

Name	L	mp		${z}_{{\rm{init}}}$	${M}_{{\rm{halo,min}}}$	${M}_{{\rm{acc,min}}}$
	(${\rm{Mpc}}\;{h}^{-1}$)	(${M}_{\odot }{h}^{-1}$)	(${\rm{kpc}}\;{h}^{-1}$)		(${M}_{\odot }{h}^{-1}$)	(${M}_{\odot }{h}^{-1}$)
NB_L10_N2048	10	$8.72\times {10}^{3}$	0.31	350	$3.49\times {10}^{6}\$	$2.18\times {10}^{7}\$
NB_L15_N2048	15	$2.94\times {10}^{4}$	0.46	350	$1.18\times {10}^{7}\$	$7.36\times {10}^{7}\$
NB_L25_N2048	25	$1.36\times {10}^{5}$	0.76	300	$5.45\times {10}^{7}\$	$3.41\times {10}^{8}\$
NB_L35_N2048	35	$3.74\times {10}^{5}$	1.07	300	$1.50\times {10}^{8}\$	$9.35\times {10}^{8}\$
NB_L50_N2048	50	$1.09\times {10}^{6}$	1.53	300	$4.36\times {10}^{8}\$	$2.73\times {10}^{9}\$
NB_L75_N2048	75	$3.68\times {10}^{6}$	2.29	300	$1.47\times {10}^{9}\$	$9.20\times {10}^{9}\$
NB_L100_N2048	100	$8.72\times {10}^{6}$	3.05	250	$3.49\times {10}^{9}\$	$2.18\times {10}^{10}$
NB_L150_N2048	150	$2.94\times {10}^{7}$	4.58	250	$1.18\times {10}^{10}$	$7.36\times {10}^{10}$
NB_L200_N2048	200	$6.98\times {10}^{7}$	6.10	200	$2.79\times {10}^{10}$	$1.74\times {10}^{11}$
NB_L300_N2048	300	$2.36\times {10}^{8}$	9.16	200	$9.42\times {10}^{10}$	$5.89\times {10}^{11}$
NB_L400_N2048	400	$5.58\times {10}^{8}$	12.2	150	$2.23\times {10}^{11}$	$1.40\times {10}^{12}$

Note. Two realizations of each P³M N-body simulation containing 2048³ dark matter particles are run using 2LPT initial conditions. Dark matter halos are found with a hybrid algorithm and have an average density ${\bar{\rho }}_{{\rm{halo}}}=200{\bar{\rho }}_{{\rm{m}}}.$ For calculating the halo mass and accretion rate functions, the minimum masses correspond to 400 and 2500 particles, respectively.

Download table as:  ASCIITypeset image

The halo mass function in differential form is defined as the comoving number density $n(\gt M,z)$ of halos of mass M per unit mass dM. In extended Press–Schechter theory (EPS; e.g., Bond et al. 1991; Lacey & Cole 1993), it is generally expressed as

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dM}}=f(\sigma )\displaystyle \frac{{\rho }_{0}}{M}\displaystyle \frac{d\mathrm{ln}{\sigma }^{-1}}{{dM}},\end{eqnarray} \tag{ 3 }$$

where ${\rho }_{0}$ is the comoving average cosmic density, σ is the rms fluctuation of the smoothed linear density field, and $f(\sigma )$ is the σ-weighted distribution of random-walk barrier-crossings. For a continuous density field, the variance of the density fluctuations is calculated as

$$\begin{eqnarray}&&{\sigma }^{2}(M,z)=\int P(k,z)| W(k,M){| }^{2}{k}^{2}{dk},\end{eqnarray} \tag{ 4 }$$

where $P(k,z)$ is the linear power spectrum extrapolated to redshift z and $W(k,M)$ is the Fourier transform of the spherical tophat window function of radius $R={[M/(4/3\pi {\rho }_{0})]}^{1/3}.$ Following Warren et al. (2006) and Tinker et al. (2008), the barrier-crossing distribution function is parametrized as

$$\begin{eqnarray}&&f(\sigma )=A\left[1+{\left(\displaystyle \frac{\sigma }{b}\right)}^{-a}\right]{e}^{-c/{\sigma }^{2}},\end{eqnarray} \tag{ 5 }$$

where A sets the overall amplitude, a and b set the slope and amplitude of the low-mass power law, and c sets the exponential decrease scale.

Missing large-scale power due to finite box sizes are corrected using a similar procedure to Reed et al. (2007). An effective variance is calculated as

$$\begin{eqnarray}&&{\sigma }_{{\rm{eff}}}^{2}(M,z)=\displaystyle \sum _{{\boldsymbol{k}}}| \delta ({\boldsymbol{k}},z){| }^{2}| W(k,M){| }^{2},\end{eqnarray} \tag{ 6 }$$

where $\delta ({\boldsymbol{k}},z)$ is the Fourier transform of the linear overdensity field extrapolated to redshift z. An effective mass ${M}_{{\rm{eff}}}$ is calculated by finding the mass in Equation (4) which would yield a variance equal to ${\sigma }_{{\rm{eff}}}^{2}.$ The effective quantities ${M}_{{\rm{eff}}}$ and ${\sigma }_{{\rm{eff}}}$ are used instead of their normal counterparts to calculate the halo mass function.

We combine all of the simulations and bin the halos to calculate the discrete halo mass function,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}N}{{\rm{\Delta }}V{\rm{\Delta }}M}=\displaystyle \sum _{i}^{N}\displaystyle \frac{{w}_{i}}{{V}_{i}{\rm{\Delta }}M}.\end{eqnarray} \tag{ 7 }$$

Each binned halo is given a weight $w=M/{M}_{{\rm{eff}}}$ such that total mass in any given mass bin is conserved. Since simulations with different box sizes have different minimum masses, the effective volumes V differ for the various mass bins.

Figure 1 shows the differential mass functions and barrier-crossing distribution functions at $z=$ 6, 8, and 10. The mass functions have a generic shape with a low-mass power-law and a high-mass decaying tail. The function $f(\sigma )$ shows no redshift dependence after correcting for finite box-size effects. Since our box sizes are reasonably large for the mass and redshift ranges considered, the corrections are typically small, only $\sim 10\%$ at z = 10 and $\sim 1\%$ at z = 6 for moderate masses.

Zoom In Zoom Out Reset image size

The middle panels of Figure 1 show a comparison with the best-fit results from Tinker et al. (2008), which are calibrated for the mass range ${10}^{11}\lt M/({M}_{\odot }{h}^{-1})\lt {10}^{15}$ and redshift range $0\leqslant z\leqslant 2.5.$ As suggested in their paper and through personal communication, we use a maximum redshift z = 2.5 in their equations for the redshift evolution of the fitting parameters to obtain $A=0.156,$ $a=1.36,$ $b=2.54,$ and c = 1.19. The agreement is remarkably good overall considering the extrapolation in both mass and redshift. It is better at lower masses and at lower redshifts, but differences of $\approx 20\%$ are seen at the high-mass end where bright high-redshift galaxies are expected to reside.

Since $f(\sigma )$ has a universal form for the scales and redshifts considered, we choose to fit the z = 6 results only because of the larger range and higher signal-to-noise. Furthermore, the parameters a and b are kept unchanged because the simulated halo catalogs do not sample the high-σ (low-mass) power law portion of the function. The best-fit barrier-crossing distribution function is

$$\begin{eqnarray}&&f(\sigma )=0.150\left[1+{\left(\displaystyle \frac{\sigma }{2.54}\right)}^{-1.36}\right]{e}^{-1.14/{\sigma }^{2}}.\end{eqnarray} \tag{ 8 }$$

The bottom panels of Figure 1 show that the residuals for the new fit are $\lesssim 5\%$ across the mass and redshift ranges of interest. The larger deficits at the highest mass and particularly at higher redshifts arise because we cannot perfectly correct for the missing large-scale power in finite simulation boxes. In the exponential tail of the halo mass function, a small mass change can lead to a relatively large change in the number density. While we have started the simulations at high enough initial redshifts such that the second-order displacement corrections are themselves small, starting at even high redshifts may reduce the mass deficits.

The halo mass function can be calculated for any cosmological model given the linear power spectrum. For our set of cosmological parameters, the tophat variance can be fit to $\lesssim 2\%$ accuracy for ${10}^{7}\lt M/({M}_{\odot }{h}^{-1})\lt {10}^{13}$ with

$$\begin{eqnarray}&&{\sigma }^{-1}(z)={D}^{-1}(z)\left[0.11+0.077{\left(\displaystyle \frac{M}{{10}^{8}\ {M}_{\odot }{h}^{-1}}\right)}^{0.18}\right],\end{eqnarray} \tag{ 9 }$$

where D(z) is the linear growth factor. Note that a double power law plus a constant term can accurately fit the entire mass range spanning mini-halos ($M\sim {10}^{5}\ {M}_{\odot }{h}^{-1}$) to massive clusters ($M\sim {10}^{15}\ {M}_{\odot }{h}^{-1}$).

The halo accretion rate function in differential form is defined as the comoving number density $n(\gt \dot{M},z)$ of halos with accretion rate $\dot{M}$ per unit accretion rate $d\dot{M}.$ In principle, its functional form in terms of accretion rate and redshift can be motivated with EPS theory. Alternatively, we take a novel approach and express it as

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{d\dot{M}}=\displaystyle \frac{{dn}}{{dM}}\displaystyle \frac{{dM}}{d\dot{M}},\end{eqnarray} \tag{ 10 }$$

where dn/dM is the differential mass function (Equation (3)). The advantage of using this functional form is that the redshift dependence of the halo mass function is already well known. To use Equation (10), we need to understand how mass and accretion rate are related and be able to calculate the derivative ${dM}/d\dot{M}.$

Figure 2 shows how the mass accretion rate depends on mass at z = 6, 8, and 10. The average accretion rates for the mass range ${10}^{8}\lt M\lt {10}^{13}$ and redshift range $6\leqslant z\leqslant 10$ are well fit by

$$\begin{eqnarray}&&\langle \dot{M}\rangle =0.21\;{M}_{\odot }\;{{\rm{yr}}}^{-1}{\left(\displaystyle \frac{M}{{10}^{8}\;{M}_{\odot }}\right)}^{1.06}{\left(\displaystyle \frac{1+z}{7}\right)}^{2.5}.\end{eqnarray} \tag{ 11 }$$

The residuals are typically $\lesssim 5\%$ across the mass and redshift ranges of interest. The larger deficits at the highest masses and particularly at higher redshifts are related to those seen in the halo mass function in Figure 1 since the accretion rates are also affected by the missing large-scale power in finite simulation boxes.

Zoom In Zoom Out Reset image size

Our mass power-law slope of 1.06 is slightly shallower than the value of 1.1 from Fakhouri et al. (2010) and slightly steeper than the linear mass dependence from Correa et al. (2015a, 2015b). We find similar results for $M\gtrsim {10}^{10}\ {M}_{\odot }$ where our mass range overlaps with theirs, but our accretion rates become relatively larger (smaller) at lower mass because of the shallower (steeper) slope. We all find the same redshift power-law slope of 2.5 at high z. Fakhouri et al. (2010) measure accretion rates from N-body simulations, but they use a FoF halo finder instead. Since the ratio of M₂₀₀ and ${M}_{{\rm{FoF}}}$ is neither constant with mass nor redshift, it is not surprising to find differences between ${\dot{M}}_{200}$ and ${\dot{M}}_{{\rm{FoF}}}.$ Their fit is done at low redshift $z\sim 0$ and it is not clear how well their best fit compares to their own data for $6\leqslant z\leqslant 10.$ Correa et al. (2015a, 2015b) use EPS theory and N-body simulations to derive and measure the mass accretion history, respectively. They find a linear dependence on mass from their analytical work and assume this same mass dependence when fitting the numerical results. Thus, it is not clear if the latter would have had a slightly higher mass power-law slope. Nonetheless, there is overall good agreement between all of the results given the differences in details.

The distribution of mass accretion rate at any given mass is positively skewed. The best-fit relations for the variance ${\mu }_{2}\equiv \langle {(\dot{M}-\langle \dot{M}\rangle )}^{2}\rangle$ and skewness ${\mu }_{3}\equiv \langle {(\dot{M}-\langle \dot{M}\rangle )}^{3}\rangle$ are

$$\begin{eqnarray}&&{\mu }_{2}^{1/2}=0.28\;{M}_{\odot }\;{{\rm{yr}}}^{-1}\;{\left(\displaystyle \frac{M}{{10}^{8}\;{M}_{\odot }}\right)}^{1.0}{\left(\displaystyle \frac{1+z}{7}\right)}^{1.1},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\mu }_{3}^{1/3}=0.42\;{M}_{\odot }\;{{\rm{yr}}}^{-1}\;{\left(\displaystyle \frac{M}{{10}^{8}\;{M}_{\odot }}\right)}^{1.0}{\left(\displaystyle \frac{1+z}{7}\right)}^{0.6}.\end{eqnarray} \tag{ 13 }$$

Equations (11)–(13) can be used with, for example, an Edgeworth expansion to model the probability distribution function of accretion rates at a given mass and redshift. Mergers are responsible for the positive tail of the distribution. Thus, the similarity in mass dependence and difference in redshift evolution between $\langle \dot{M}\rangle ,$ ${\mu }_{2}^{1/2},$ and ${\mu }_{3}^{1/3}$ reflect those between the average accretion rate and the average merger rate. See Fakhouri et al. (2010) for recent work on halo merger rates.

Figure 3 shows the abundance of dark matter halos as a function of mass accretion rate at z = 6, 8, and 10. The halo accretion rate function and the halo mass function are similar in shape, with each having a low-end power law and a high-end exponential decline. We can model the halo accretion rate function using Equation (10) and by combining the halo mass function with the best-fit mediating mass,

$$\begin{eqnarray}&&M=4.2\times {10}^{8}\;{M}_{\odot }{\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{0.91}{\left(\displaystyle \frac{1+z}{7}\right)}^{-2.4}.\end{eqnarray} \tag{ 14 }$$

The bottom panel of Figure 3 shows that the residuals for the best fit are typically $\lesssim 10\%$ for the mass and redshift ranges of interest. The larger deficits at the highest $\dot{M}$ and particularly at higher redshifts are related to those seen in the halo mass function in Figure 1 since the accretion rates are also affected by the missing large-scale power in finite simulation boxes. Note that Equation (14) is not simply obtained from Equation (11) because the distribution of accretion rate at a given mass is skewed. However, they both have the convenient property of separable power-law dependences on M, $\dot{M},$ and z.

Zoom In Zoom Out Reset image size

The halo accretion rate function does not have a simple redshift dependence according to Equation (10). For a given accretion rate, both the mediating mass and its derivative increase at lower redshifts. This is a consequence of accretion slowing down as seen in Figure 2. In contrast, the halo mass function at a higher mass decreases in amplitude because more massive halos are rarer in hierarchical structure formation. We expect the redshift dependence in Equations (11)–(14) to hold for redshifts $z\gtrsim 6$ relevant to the EoR. However, our single redshift power law scaling is not appropriate at lower redshifts. McBride et al. (2009) have shown that a more complex scaling is required to parametrize the entire mass accretion history. More work is required to examine the validity and accuracy of extrapolating our results to lower redshifts.

Abundance matching is performed by equating the number density of galaxies to the number density of halos in differential form,

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{{\rm{gal}}}}{{{dL}}_{{\rm{UV}}}}=\displaystyle \frac{{{dn}}_{{\rm{halo}}}}{{dX}}\displaystyle \frac{{dX}}{{{dL}}_{{\rm{UV}}}},\end{eqnarray} \tag{ 15 }$$

or in cumulative form,

$$\begin{eqnarray}&&{n}_{{\rm{gal}}}(\gt {L}_{{\rm{UV}}},z)={n}_{{\rm{halo}}}(\gt X,z).\end{eqnarray} \tag{ 16 }$$

If X is taken to be the halo mass M, then we infer the luminosity–mass relation ${L}_{{\rm{UV}}}(M,z).$ This procedure does not account for scatter in the mass-to-light ratio and the relation should be interpreted as an average luminosity for a given mass. If X is taken to be the mass accretion rate $\dot{M},$ then we infer the luminosity–accretion rate relation ${L}_{{\rm{UV}}}(\dot{M},z).$ This procedure accounts for the episodic nature of star formation and scatter in the mass-to-light ratio. Furthermore, it is more logical since the star formation rate should be more highly correlated with the halo growth rate rather than halo mass.

The galaxy luminosity function is generally parametrized with a Schechter function in luminosity form,

$$\begin{eqnarray}&&\phi ({L}_{{\rm{UV}}})={\phi }_{\star }{\left(\displaystyle \frac{{L}_{{\rm{UV}}}}{{L}_{\star }}\right)}^{\alpha }\mathrm{exp}\left(-\displaystyle \frac{{L}_{{\rm{UV}}}}{{L}_{\star }}\right),\end{eqnarray} \tag{ 17 }$$

or in magnitude form,

$$\begin{eqnarray}\phi ({M}_{{\rm{UV}}}) & = & (0.4\mathrm{ln}10){\phi }_{\star }{\left[{10}^{0.4({M}_{\star }-{M}_{{\rm{UV}}})}\right]}^{\alpha +1}\\ & & \times \mathrm{exp}\left[{10}^{0.4({M}_{\star }-{M}_{{\rm{UV}}})}\right],\end{eqnarray} \tag{ 18 }$$

where ${\phi }_{\star }$ is an overall normalization, ${L}_{\star }$ is a characteristic luminosity, ${M}_{\star }$ is a characteristic magnitude, and α is the slope of the faint-end power law. The conversion between UV magnitude and luminosity is given by the standard AB relation,

$$\begin{eqnarray}&&{M}_{{\rm{UV}}}=-2.5\mathrm{log}\left(\displaystyle \frac{{L}_{{\rm{UV}}}}{4.487\times {10}^{20}\ \mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Hz}}^{-1}}\right).\end{eqnarray} \tag{ 19 }$$

The cumulative galaxy luminosity function is obtained by integrating Equation (17) analytically to get the incomplete gamma function ${\phi }_{\star }{L}_{\star }{\rm{\Gamma }}(1+\alpha ,{L}_{{\rm{UV}}}/{L}_{\star }).$ The cumulative halo mass function and mass accretion rate function do not have simple analytical forms, but are easily obtained by numerical integration.

Bouwens et al. (2015a) have made the most recent measurements of the high-redshift UV luminosity functions. For three redshift samples with $\langle z\rangle \approx 5.9,$ 7.9, and 10.4, there are 867, 217, and 6 galaxy candidates observed, respectively. For reference, their best-fit Schechter parameters are:

$$\begin{eqnarray*}\langle z\rangle \approx 5.9,{\phi }_{\star } & = & 5.0\times {10}^{-4},\ {M}_{\star }=-20.94,\ \alpha =-1.87\\ \langle z\rangle \approx 7.9,{\phi }_{\star } & = & 2.1\times {10}^{-4},\ {M}_{\star }=-20.63,\ \alpha =-2.02\\ \langle z\rangle \approx 10.4,\ {\phi }_{\star } & = & 8.0\times {10}^{-5},{M}_{\star }=-20.92,\alpha =-2.27\\ z\to 10.4,\ {\phi }_{\star } & = & 3.0\times {10}^{-4},\ {M}_{\star }=-20.92,\ \alpha =-2.27\end{eqnarray*}$$

where ${\phi }_{\star }$ has units of comoving Mpc⁻³. Note that the last entry above is not a best fit to the galaxy counts, but an extrapolation of the luminosity functions from lower redshifts. They allowed the parameters $\mathrm{log}{\phi }_{\star },$ ${M}_{\star },$ and α to vary linearly with z and fitted the galaxy counts in the redshift range $4\lesssim z\lesssim 8$ in order to estimate the results for higher redshifts. Using the derived redshift dependence, they were able to fix ${M}_{\star }$ and α and then fit for ${\phi }_{\star }$ at $\langle z\rangle \approx 10.4.$ We note that ${M}_{\star }=-20.92$ at $z\approx 10$ is rather bright considering that the best-fit characteristic magnitude changed from −20.94 to −20.63 for z ≈ 6–8. Furthermore, the extrapolation from a different redshift range $5\lesssim z\lesssim 8$ yields ${M}_{\star }\approx -20.22$ (Bouwens et al. 2015b), which is more in line with the best-fit values for $z\approx 6$–8.

To estimate upper and lower bounds on a given galaxy luminosity function, we vary the three Schechter parameters using their uncertainties, add the weighted Schechter functions to the ${M}_{{\rm{UV}}}$-ϕ grid, and calculate the 68% confidence level. For the weight, we use a trivariate Gaussian likelihood that assumes uncorrelated errors. However, there is degeneracy between the parameters and their errors are correlated. To prevent overestimating the confidence region, we reduce the errors in the Schechter parameters to 2/3 of their original values. This procedure allows us to match the $1\sigma$ errors in the binned measurements of the galaxy luminosity function. A more rigorous statistical analysis would require using the full likelihood for the galaxy luminosity function fit, which is not available.

Figure 4 shows how the UV luminosity and magnitude depend on mass and accretion rate at $z\approx 6,$ 8, and 10 (i.e., 5.9, 7.9, and 10.4). The shaded curves are illustrative of the $1\sigma$ (68%) uncertainty in the current galaxy luminosity function. Our error analysis is able to capture the trend of increasing uncertainty at higher redshifts. Neither the luminosity–mass relation nor the luminosity–accretion-rate relation is represented by a single power law, but appears to be composed of several power law portions.

Zoom In Zoom Out Reset image size

To qualitatively understand the shape of the luminosity–mass and luminosity–accretion-rate relations, we also plot their effective power-law slopes in the bottom panels of Figure 4. The cumulative galaxy and halo number densities can be written in effective power form as ${n}_{{\rm{gal}}}(L)\propto {L}^{1+{\alpha }_{{\rm{eff}}}},$ ${n}_{{\rm{halo}}}(M)\propto {M}^{1+{\beta }_{{\rm{eff}}}},$ and ${n}_{{\rm{halo}}}(\dot{M})\propto {\dot{M}}^{1+{\gamma }_{{\rm{eff}}}}.$ The luminosity–mass and luminosity–accretion-rate relations can be written as ${L}_{{\rm{UV}}}\propto {M}^{{\delta }_{{\rm{eff}}}}$ and ${L}_{{\rm{UV}}}\propto {\dot{M}}^{{\epsilon }_{{\rm{eff}}}}.$ From Equation (15), the slopes are related through

$$\begin{eqnarray}&&{\delta }_{{\rm{eff}}}=\displaystyle \frac{d\mathrm{ln}{L}_{{\rm{UV}}}}{d\mathrm{ln}M}=\displaystyle \frac{d\mathrm{ln}{n}_{{\rm{halo}}}}{d\mathrm{ln}M},//\displaystyle \frac{d\mathrm{ln}{n}_{{\rm{gal}}}}{d\mathrm{ln}{L}_{{\rm{UV}}}}=\displaystyle \frac{1+{\beta }_{{\rm{eff}}}}{1+{\alpha }_{{\rm{eff}}}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\epsilon }_{{\rm{eff}}}=\displaystyle \frac{d\mathrm{ln}{L}_{{\rm{UV}}}}{d\mathrm{ln}\dot{M}}=\displaystyle \frac{d\mathrm{ln}{n}_{{\rm{halo}}}}{d\mathrm{ln}\dot{M}},//\displaystyle \frac{d\mathrm{ln}{n}_{{\rm{gal}}}}{d\mathrm{ln}{L}_{{\rm{UV}}}}=\displaystyle \frac{1+{\gamma }_{{\rm{eff}}}}{1+{\alpha }_{{\rm{eff}}}}.\end{eqnarray} \tag{ 21 }$$

At low ${L}_{{\rm{UV}}},$ M, and $\dot{M},$ where the galaxy luminosity function and halo mass and accretion rate functions are approximately power laws (${\alpha }_{{\rm{eff}}}\approx \alpha ,$ ${\beta }_{{\rm{eff}}}\approx -2,$ ${\gamma }_{{\rm{eff}}}\approx -2$), the effective power-law slopes of the luminosity–mass and luminosity–accretion-rate relations are approximately constant ${\delta }_{{\rm{eff}}}\approx {\epsilon }_{{\rm{eff}}}\approx -1/(1+\alpha ).$ At intermediate scales, the halo mass and accretion rate functions deviate from their power law forms before the galaxy luminosity function does, both ${\beta }_{{\rm{eff}}}$ and ${\gamma }_{{\rm{eff}}}$ become more negative, and both ${\delta }_{{\rm{eff}}}$ and ${\epsilon }_{{\rm{eff}}}$ increase. As the galaxy luminosity function deviates from its power law form, ${\alpha }_{{\rm{eff}}}$ becomes more negative, and both ${\delta }_{{\rm{eff}}}$ and ${\epsilon }_{{\rm{eff}}}$ decrease again. At the very bright and massive end, both the galaxy and halo number densities are exponentially decreasing, but such that ${\alpha }_{{\rm{eff}}}\lt {\beta }_{{\rm{eff}}}$ and ${\alpha }_{{\rm{eff}}}\lt {\gamma }_{{\rm{eff}}},$ resulting in ${\delta }_{{\rm{eff}}}$ and ${\epsilon }_{{\rm{eff}}}$ being approximately constant. The effective power-law slopes suggest that both relations can be fit with, for example, a triple power law,

$$\begin{eqnarray}&&{L}_{{\rm{UV}}}={L}_{0}{\left(\displaystyle \frac{X}{{X}_{a}}\right)}^{a}{\left(1+\displaystyle \frac{X}{{X}_{b}}\right)}^{b-a}{\left(1+\displaystyle \frac{X}{{X}_{c}}\right)}^{c-b},\end{eqnarray} \tag{ 22 }$$

where L₀ is an overall amplitude, a, b, and c are three power-law slopes, and ${X}_{a},$ ${X}_{b},$ and ${X}_{c}$ are three characteristic scales.

For the luminosity–mass relation, brighter galaxies are found in more massive halos at any given redshift as expected from the abundance matching procedure. At any given mass, brighter galaxies are found at higher redshifts and this is a consequence of the relative evolution of the galaxy luminosity and halo mass functions. Our results are similar to Kuhlen & Faucher-Giguère (2012), where they performed abundance matching using older fits to the galaxy luminosity function and halo mass function.

For the luminosity–accretion-rate relation, brighter galaxies are found in halos with larger accretion rates at any given redshift by construction. The z ≈ 6–8 results are remarkably very similar for the entire accretion rate of interest. In fact, the luminosity–accretion-rate relation is consistent with no evolution for the redshift range $6\lesssim z\lesssim 10$ when taking into account the current uncertainties in the galaxy luminosity functions. In Section 7, we show that our results at $z\approx 10$ are consistent with the binned measurements of the galaxy luminosity function (e.g., Oesch et al. 2014; Bouwens et al. 2015a).

The luminosity–accretion-rate relation at $z\approx 6$ can be fit with a triple power law,

$$\begin{eqnarray}{L}_{{\rm{UV}}} & = & 1.5\times {10}^{25}\ \mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Hz}}^{-1}{\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{1.15}\\ & & \times {\left(1+\displaystyle \frac{\dot{M}}{3.5\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{0.25}\\ & & \times {\left(1+\displaystyle \frac{\dot{M}}{1000\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{-0.95},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}{M}_{{\rm{UV}}} & = & -11.3-2.88\mathrm{log}\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)\\ & & -0.63\mathrm{log}\left(1+\displaystyle \frac{\dot{M}}{3.5\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)\\ & & +2.38\mathrm{log}\left(1+\displaystyle \frac{\dot{M}}{1000\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right).\end{eqnarray} \tag{ 24 }$$

Because of the degeneracy between the triple power law parameters, we choose to fix the low, middle, and high-$\dot{M}$ slopes at 1.15, 1.40, and 0.45, respectively based on Figure 4. There are still degeneracies between the amplitude and the characteristic accretion rates. The residuals for the fit are $\lesssim 5\%$ across the accretion rate range of interest. If we assume no evolution in the luminosity–accretion-rate relation, then Equations (23) and (24) can be used as a universal EoR template to construct a fiducial model for the evolution of the galaxy luminosity function as discussed in Section 7.

To better understand the high-redshift galaxy–halo connection from abundance matching, we examine the implications of the inferred luminosity–mass and luminosity–accretion rate relations on the star formation rate and efficiency. The star formation rate ${\dot{M}}_{{\rm{s}}}$ can be estimated from the UV luminosity using the standard relation from Madau et al. (1998):

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{s}}}=\left(\displaystyle \frac{{L}_{{\rm{UV}}}}{8\times {10}^{27}\ \mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Hz}}^{-1}}\right){M}_{\odot }\;{{\rm{yr}}}^{-1},\end{eqnarray} \tag{ 25 }$$

assuming a Salpeter initial mass function (IMF) truncated at 0.1 ${M}_{\odot }$ and 125 ${M}_{\odot }.$ The normalization is uncertain and depends on the stellar IMF and formation history, but we adopt this commonly used relation to allow a wider comparison with other work.

We calculate the star formation efficiency in different ways for the two abundance matching results. In the case of the luminosity–accretion-rate relation, the star formation efficiency is defined as

$$\begin{eqnarray}&&{\varepsilon }_{\dot{M}}\equiv \displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{{\dot{M}}_{{\rm{b}}}},\end{eqnarray} \tag{ 26 }$$

where the baryonic mass accretion rate,

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{b}}}=\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{b}}}}{{{\rm{\Omega }}}_{{\rm{m}}}}\dot{M}\propto {(1+z)}^{2.5},\end{eqnarray} \tag{ 27 }$$

is calculated assuming the cosmic baryon fraction. In the case of the luminosity–mass relation, we assume no information about accretion rates and define the efficiency as

$$\begin{eqnarray}&&{\varepsilon }_{M}\equiv \displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{{M}_{{\rm{b}}}/{t}_{{\rm{dyn}}}},\end{eqnarray} \tag{ 28 }$$

where the baryonic mass consumption rate,

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{b}}}}{{t}_{{\rm{dyn}}}}=\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{b}}}}{{{\rm{\Omega }}}_{{\rm{m}}}}M\right){\left(\displaystyle \frac{32G{\bar{\rho }}_{{\rm{halo}}}}{3}\right)}^{1/2}\propto {(1+z)}^{1.5},\end{eqnarray} \tag{ 29 }$$

depends on the redshift-dependent halo dynamical time. Equation (29) can be considered the three-dimensional (3D) analog of the observed Schmidt–Kennicutt relation and is often used in semi-analytical models and cosmological simulations of galaxy formation. The star formation efficiency is expected to range from ∼0.01 to ∼0.1, but keep in mind that the normalization is uncertain.

We impose a star formation limit based on the criterion that galaxies form in dark matter halos where the gas cools efficiently through atomic transitions. For a minimum virial temperature ${T}_{{\rm{min}}}={10}^{4}$ K, the minimum mass for a halo to host a galaxy is defined as

$$\begin{eqnarray}&&{M}_{{\rm{min}}}(z)\equiv {\left(\displaystyle \frac{2{{kT}}_{{\rm{min}}}}{\mu {m}_{H}G}\right)}^{3/2}{\left(\displaystyle \frac{4\pi }{3}{\bar{\rho }}_{{\rm{halo}}}\right)}^{-1/2}\\ &&\approx 1.6\times {10}^{8}\ {M}_{\odot }{\left(\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{halo}}}}{200}\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}{h}^{2}}{0.132}\right)}^{-1/2}{\left(\displaystyle \frac{1+z}{7}\right)}^{-3/2},\end{eqnarray} \tag{ 30 }$$

where the mean atomic mass is set as $\mu =0.6$ for a fully ionized intrahalo medium. For our fiducial model assuming a universal EoR luminosity–accretion rate relation, the corresponding limiting magnitude can be fit with the linear relation,

$$\begin{eqnarray}&&{M}_{{\rm{UV,SF}}}(z)=-10.0-0.35(z-6).\end{eqnarray} \tag{ 31 }$$

Thus, we truncate all relevant curves in our figures at the star formation limit.

Figure 5 shows the inferred star formation efficiencies from abundance matching ${L}_{{\rm{UV}}}$ with M and with $\dot{M}.$ The uncertainties in the efficiencies correspond to the those in the luminosity–mass and luminosity–accretion-rate relations shown in Figure 4. The efficiency is not monotonic with mass nor accretion rate at any given redshift, but has a maximum value at a characteristic peak scale near where the galaxy luminosity function transitions from a power law to an exponential decline. This peak occurs at a characteristic mass $\sim 2\times {10}^{11}\ {M}_{\odot }$ and a characteristic accretion rate $\sim 6\times {10}^{2}\ {M}_{\odot }\;{{\rm{yr}}}^{-1}$ at $z\approx 6.$

Zoom In Zoom Out Reset image size

The dependence of the star formation efficiency on mass and accretion rate has a physical explanation. The reduced efficiency at higher M and $\dot{M}$ is consistent with relatively inefficient atomic cooling and cold gas accretion in larger halos. The atomic cooling rate is relatively low at $T\gtrsim {10}^{6}$ K (e.g., Sutherland & Dopita 1993) and this coincides with the observed peak mass and accretion rate. The reduced efficiency at lower M and $\dot{M}$ is consistent with feedback effects from photoheating and supernova. Smaller halos are more severely affected by feedback because of the lower binding energy, which scales as ${M}^{5/3}.$ It is possible that the assumed Schechter form for the galaxy luminosity function, in particular the bright-end exponential decline and the extrapolated faint-end power law, may not be accurate and therefore would bias the inferred star formation efficiency. However, there is support for the nonmonotonic efficiency, which has been observed at lower redshifts (e.g., Guo et al. 2010; Leauthaud et al. 2012; Behroozi et al. 2013b).

We find that the star formation efficiency ${\varepsilon }_{\dot{M}}$ inferred from abundance matching ${L}_{{\rm{UV}}}$ and $\dot{M}$ is more consistent with having no redshift evolution than the efficiency ${\varepsilon }_{M}$ inferred from abundance matching ${L}_{{\rm{UV}}}$ and M. Behroozi et al. (2013a) have also found that there is lack of evolution in the star formation efficiency ${\varepsilon }_{\dot{M}}$ for the redshift range $0\lesssim z\lesssim 8.$ Finkelstein et al. (2015) have shown that the stellar baryon fraction increases at higher redshifts over the range $4\lesssim z\lesssim 8,$ but this is consistent with both of our abundance matching results and does not obviously favor one particular scenario. Our results suggest that the star formation rate evolves more like ${(1+z)}^{2.5}$ rather than the commonly assumed ${(1+z)}^{1.5}$ scaling from the 3D analog of the Schmidt–Kennicutt relation. More precise measurements of the galaxy luminosity function and a more rigorous statistical analysis are required to strengthen this argument.

We note that the comparison of the exponential tails in the galaxy and halo abundances is very sensitive to the value of ${M}_{\star }$ in the Schechter luminosity function. In Figure 5, the star formation efficiency curves are shifted to higher M and $\dot{M}$ because having ${M}_{\star }=-20.92$ at $z\approx 10$ is rather bright considering that the best-fit characteristic magnitude changed from −20.94 to −20.63 for z ≈ 6–8. Recall that the extrapolation of the observed galaxy luminosity functions from a different redshift range $5\lesssim z\lesssim 8$ yields ${M}_{\star }\approx -20.22$ (Bouwens et al. 2015b). With this fainter characteristic magnitude, we expect the peaks in the star formation efficiency ${\varepsilon }_{\dot{M}}$ in Figure 5 (right) to be even more aligned.

The efficiency at $z\approx 6$ can be fit with a triple power law,

$$\begin{eqnarray}{\varepsilon }_{\dot{M}} & = & 1.1\times {10}^{-2}\;{\left(\displaystyle \frac{\dot{M}}{{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{0.15}\\ & & \times {\left(1+\displaystyle \frac{\dot{M}}{3.5\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{0.25}\\ & & \times {\left(1+\displaystyle \frac{\dot{M}}{1000\;{M}_{\odot }\;{{\rm{yr}}}^{-1}}\right)}^{-0.95}.\end{eqnarray} \tag{ 32 }$$

Equation (32) corresponds to the best fit for the luminosity–accretion-rate relation in Equation (23). This universal EoR template can be used in semi-analytical calculations and cosmological simulations to model star formation in high-redshift galaxies during the EoR. Recently, Mashian et al. (2015) have also applied the abundance matching technique to calibrate a relation between star formation rate and halo mass. They find a double power law scaling relation for $M\gt {10}^{10}{M}_{\odot },$ which is consistent with our results when the same mass range is considered.

To understand how the abundance of galaxies as a function of luminosity and redshift is connected to the abundance and growth of dark matter halos, it is informative to relate the galaxy luminosity function to the halo accretion rate function and the luminosity–accretion-rate relation. The luminosity function in magnitude form can be calculated as

$$\begin{eqnarray}&&\phi ({M}_{{\rm{UV}}},z)=\displaystyle \frac{{dn}}{d\dot{M}}\displaystyle \frac{d\dot{M}}{{{dM}}_{{\rm{UV}}}},\end{eqnarray} \tag{ 33 }$$

where ${dn}/d\dot{M}$ is obtained using Equation (10). Since the luminosity–accretion-rate relation and the star formation efficiency are consistent with no evolution for the redshift range $6\lesssim z\lesssim 10,$ it is intriguing to predict the galaxy luminosity function using the universal EoR template for ${M}_{{\rm{UV}}}(\dot{M})$ given by Equation (24). If more precise measurements of the luminosity function from upcoming observations clearly deviate from this fiducial model, then it would point to exciting, additional astrophysics in star and galaxy formation.

Figure 6 shows the predicted galaxy luminosity function at $z\approx 6,$ 8, and 10 for our fiducial model. The shaded prediction curves reflect the $1\sigma$ (68%) uncertainty in the luminosity–accretion rate relation. The binned observational measurements at $6\lesssim z\lesssim 10$ from Bouwens et al. (2015a) and at $z\approx 10$ from Oesch et al. (2014) are shown for comparison. Our results match at $z\approx 6$ by construction and are in very good agreement at $z\approx 8.$ At $z\approx 10,$ our model generally has larger amplitude than the few observational data points. This explains the apparent redshift dependence in the luminosity–accretion rate relation (Figure 4) and the star formation efficiency (Figure 5). However, our model is still consistent with the highly uncertain observations. Thus, the luminosity–accretion rate relation (Equation (24)) and the star formation efficiency (Equation (32)) are highly consistent with no evolution and can be used as universal EoR templates for the EoR. Mashian et al. (2015) and Mason et al. (2015) have also predicted similar evolution for the galaxy luminosity function by relating observed star formation rate to halo mass.

Zoom In Zoom Out Reset image size

Table 2 shows the best-fit Schechter parameters for our fiducial model for the redshift range $6\leqslant z\leqslant 15.$ At higher redshifts, the normalization ${\phi }_{\star }$ decreases dramatically, the characteristic magnitude ${M}_{\star }$ is more positive (fainter), and the faint-end slope α is more negative (steeper). While there are degeneracies between the best-fit parameters, all three must vary with redshift to have a good fit. Note that the differences:

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{log}\phi (z)=\mathrm{log}\phi (z)-\mathrm{log}\phi (6),\end{eqnarray}$$

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{{\rm{UV}}}(z)={M}_{{\rm{UV}}}(z)-{M}_{{\rm{UV}}}(6),\end{eqnarray}$$

and

$$\begin{eqnarray}&&{\rm{\Delta }}\alpha (z)=\alpha (z)-\alpha (6)\end{eqnarray}$$

are more general and can be applied to estimate the evolution of the galaxy luminosity function when a different calibration at $z\approx 6$ is used. For better accuracy, we suggest repeating the abundance matching procedure using the new galaxy luminosity function and the halo accretion rate function.

Table 2.  Galaxy Luminosity Function Parameters

z	${\phi }_{\star }$	${M}_{\star }$	α
6	$4.9\times {10}^{-4}$	−20.92	−1.88
7	$3.6\times {10}^{-4}$	−20.76	−1.94
8	$2.4\times {10}^{-4}$	−20.61	−2.01
9	$1.5\times {10}^{-4}$	−20.46	−2.08
10	$8.9\times {10}^{-5}$	−20.32	−2.15
11	$5.0\times {10}^{-5}$	−20.18	−2.22
12	$2.7\times {10}^{-5}$	−20.04	−2.30
13	$1.4\times {10}^{-5}$	−19.91	−2.38
14	$7.0\times {10}^{-6}$	−19.78	−2.45
15	$3.4\times {10}^{-6}$	−19.64	−2.53

Note. The best-fit Schechter parameters for our fiducial model. The normalization ${\phi }_{\star }$ has units of Mpc⁻³.

Download table as:  ASCIITypeset image

We can explain the redshift dependence of the Schechter parameters by considering what happens to the halo accretion rate function (Figure 3) with increasing redshift. The overall amplitude of the halo accretion rate function decreases and the normalization ${\phi }_{\star }$ also decreases accordingly. The characteristic accretion rate, at which the halo abundance undergoes exponential decline, shifts to lower values. Correspondingly, the characteristic luminosity ${L}_{\star }$ decreases and the characteristic magnitude ${M}_{\star }$ becomes more positive (fainter). Halos at a given accretion rate are rarer and the slope of the halo accretion rate function becomes steeper. Consequently, the faint-end slope α becomes more negative.

If more precise measurements of the galaxy luminosity function find different amplitudes at higher redshifts compared to our fiducial model, then it could be due to a number of reasons. It may simply be that the galaxy luminosity function at $z\approx 6$ used to calibrate our model may be overestimated or underestimated given the current measurement errors. Or perhaps it is due to additional astrophysics in star and galaxy formation. Lower metallicity and less metal cooling could lead to a larger Jeans mass for molecular clouds and therefore result in a more top-heavy IMF. Less photoheating and less Jeans smoothing of the IGM at an earlier stage in reionization could result in more efficient gas accretion and star formation. This would affect the formation and abundance of dwarf galaxies and lead to different values for the faint-end slope. More accurate measurements are required to understand the dependence of the star formation efficiency on mass, accretion rate, and redshift and to shed light on the galaxy formation process.

Future observations with the JWST will provide an exciting and important window to the EoR. It will allow more precise measurements of the galaxy luminosity function, particularly at lower luminosities and higher redshifts. According to Windhorst et al. (2014, and through personal communication), JWST has the sensitivity to reach AB apparent magnitude ${m}_{{\rm{AB}}}\approx 31$ with deep observations and even ${m}_{{\rm{AB}}}\approx 32$ with ultra deep observations. They have suggested the strategy of observing a large number of deep fields and a much larger number of medium-deep surveys on gravitational lensing foreground targets.

Figure 7 and Table 3 show the forecast based on our fiducial model for the evolution of the galaxy luminosity function. From z = 6 to z = 15, the limiting absolute magnitude ${M}_{{\rm{AB}}}$ gets brighter by $\approx 2.3$ mag as the luminosity distance increases by a factor of $\approx 2.8.$ The comoving number density of observable galaxies $n(\lt {M}_{{\rm{AB}}})$ drops by about 4–5 orders of magnitude. The number of observable galaxies per square degree per unit redshift ${dN}(\lt {M}_{{\rm{AB}}})/{dz}$ changes similarly since the differential comoving volume dV/dz only changes by a factor of $\approx 2.2$. Mashian et al. (2015) and Mason et al. (2015) have also made similar forecasts for JWST using their predictions for the evolution of the galaxy luminosity function.

Zoom In Zoom Out Reset image size

Table 3.  Forecast of Galaxy Counts for JWST

z	M31	M32	$n({\lt M}_{31})$	$n({\lt M}_{32})$	${dN}({\lt M}_{31})/{dz}$	${dN}({\lt M}_{32})/{dz}$
6	−17.9	−16.9	$5.1\times {10}^{-3}$	$1.3\times {10}^{-2}$	50000	130000
7	−18.3	−17.3	$2.4\times {10}^{-3}$	$6.8\times {10}^{-3}$	21000	60000
8	−18.6	−17.6	$9.7\times {10}^{-4}$	$3.2\times {10}^{-3}$	7800	26000
9	−18.9	−17.9	$3.5\times {10}^{-4}$	$1.4\times {10}^{-3}$	2600	9900
10	−19.2	−18.2	$1.2\times {10}^{-4}$	$5.3\times {10}^{-4}$	760	3500
11	−19.4	−18.4	$3.4\times {10}^{-5}$	$1.9\times {10}^{-4}$	210	1200
12	−19.6	−18.6	$8.7\times {10}^{-6}$	$6.1\times {10}^{-4}$	49	340
13	−19.8	−18.8	$2.1\times {10}^{-6}$	$1.8\times {10}^{-5}$	11	94
14	−20.0	−19.0	$4.2\times {10}^{-7}$	$5.0\times {10}^{-6}$	2	24
15	−20.2	−19.2	$7.8\times {10}^{-8}$	$1.2\times {10}^{-6}$	0.4	6

Note. The absolute magnitude ${M}_{x}$ corresponds to an apparent magnitude limit ${m}_{{\rm{AB}}}=x.$ Cumulative comoving number density $n(\lt {M}_{{\rm{AB}}})$ has units of Mpc⁻³ and cumulative galaxy count per angular area per unit redshift ${dN}(\lt {M}_{{\rm{AB}}})/{dz}$ has units of deg⁻².

Download table as:  ASCIITypeset image

JWST has the sensitivity to observe $\gtrsim 11$ (24) unlensed galaxies per square degree per unit redshift at least down to ${M}_{\star }$ at $z\lesssim 13$ (14) with deep (ultra deep) observations. It will be able to probe some portion of the faint end at lower redshifts, but it is still about 8 (7) magnitudes away from the expected star formation limit. At $z\gtrsim 13$ (14), JWST will mainly probe the exponential tail of the luminosity function, where the rarity of galaxies will make them difficult to find. With apparent magnitude limit ${m}_{{\rm{AB}}}\approx 32,$ it can actually reach ${M}_{\star }$ at $z\approx 15.$ Note that our forecasted counts will increase or decrease depending on how the actual galaxy luminosity function at $z\approx 6$ compares to the one used to calibrate our fiducial model.

Our fiducial model for the evolution of the galaxy luminosity function and the corresponding forecast for JWST are for our assumed cosmological parameters. If a different cosmological model is considered, then the results at $z\approx 6$ would remain the same since this is where the calibration is done, but the higher redshift results would change. For example, consider what happens if ${{\rm{\Omega }}}_{{\rm{m}}}$ or ${\sigma }_{8}$ is increased. The ratio of the halo accretion rate function at redshift $z\gt 6$ relative to that at the pivot z = 6 will increase. As a result, the amplitude of the galaxy luminosity function and forecast for JWST at higher redshifts will also increase compared to our fiducial predictions.

The reionization history is determined by the evolving abundance of escaped ionization photons, which depends on the UV luminosity function, the spectral energy distribution, and the radiation escape fraction of high-redshift galaxies. Only a small fraction of the ionizing photon budget comes from currently observable galaxies with ${M}_{{\rm{UV}}}\lesssim -17$ and at $z\lesssim 10.$ In order to calculate the total budget, we need to extrapolate the luminosity function to fainter magnitudes and higher redshifts. With our fiducial model, we can more confidently integrate the luminosity function because of the well-behaved and physically motivated evolution in the Schechter parameters.

Figure 8 shows our prediction for the cumulative ionizing photon number density,

$$\begin{eqnarray}&&{n}_{\gamma }(\gt z)={\displaystyle \int }_{z}^{{z}_{{\rm{max}}}}{\dot{n}}_{\gamma }(z)\left|\displaystyle \frac{{dt}}{{dz}}\right|{dz},\end{eqnarray} \tag{ 34 }$$

where the photon production rate density is calculated as

$$\begin{eqnarray}&&{\dot{n}}_{\gamma }(z)={\displaystyle \int }_{{M}_{{\rm{bright}}}}^{{M}_{{\rm{faint}}}}\phi ({M}_{{\rm{UV}}},z){\dot{N}}_{\gamma }({M}_{{\rm{UV}}}){{dM}}_{{\rm{UV}}}.\end{eqnarray} \tag{ 35 }$$

A galaxy with Population II stars produces ionizing photons at a rate given by (e.g., Bruzual & Charlot 2003; Schaerer 2003)

$$\begin{eqnarray}&&{\dot{N}}_{\gamma }\approx {10}^{45.9-0.4{M}_{{\rm{UV}}}}\;{{\rm{s}}}^{-1}\approx {10}^{25.2}\;{{\rm{s}}}^{-1}\left(\displaystyle \frac{{L}_{{\rm{UV}}}}{\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Hz}}^{-1}}\right).\end{eqnarray} \tag{ 36 }$$

For the integration limits, we choose ${z}_{{\rm{max}}}=25$ and ${M}_{{\rm{bright}}}={M}_{\star }-5.$ For ${M}_{{\rm{faint}}},$ we choose the star formation limit given by Equation (31) or the commonly assumed limit ${M}_{{\rm{UV}}}=-10$ (e.g., Bouwens et al. 2012; Robertson et al. 2013) for comparison. We also show the contribution probed by JWST with apparent magnitude limit ${m}_{{\rm{UV}}}\approx 31$–32.

Zoom In Zoom Out Reset image size

Our fiducial model for the galaxy luminosity function at $6\lesssim z\lesssim 8$ is in very good agreement with current observations. At $z\approx 10,$ it generally has higher amplitude but shallower faint-end slope than the highly uncertain, best-fit observations. Increasing (reducing) the ionizing photon production rate at higher redshifts would hasten (delay) the start of reionization and lengthen (shorten) its duration. Subsequently, this would affect integrated measurements of electron scattering on the CMB such as patchy Thomson scattering, kinetic Sunyaev–Zel'dovich temperature anisotropy, and polarization anisotropy. However, there is still no tension with current EoR constraints since we poorly understand how the radiation escape fraction varies with halo mass, galaxy luminosity, and redshift.

We impose a physically motivated star formation limit based on the criterion that galaxies form in dark matter halos where the gas cools efficiently through atomic transitions. For our fiducial model, we find ${M}_{{\rm{UV,SF}}}\approx -10$ at z = 6, but it becomes more negative (brighter) with increasing redshift. This also reduces the ionizing photon production rate at higher redshifts compared to calculations using the commonly assumed limit ${M}_{{\rm{UV}}}=-10$ at all redshifts. Careful consideration is necessary when integrating the galaxy luminosity function down to the faintest limit, especially if the faint-end slope is steep. For example, if we extrapolate the fitting formula for the evolution of the Schechter parameters from Bouwens et al. (2015a) and integrate down to the commonly assumed magnitude limit, then the photon counts increase by more than an order of magnitude at high redshifts.

The cumulative photon count per hydrogen atom reaches unity at $z\approx 10$ and therefore, reionization is completed below this redshift since we have yet to account for the radiation escape fraction. Assuming reionization ended at $z\gt 6,$ the luminosity-weighted radiation escape fraction is constrained by $\langle {f}_{{\rm{esc}}}\rangle \gt 0.1.$ The lower limit is actually higher than this because additional photons are required to balance recombinations in the clumpy IGM. Since we observe radiation escape fractions ${f}_{{\rm{esc}}}\sim 0.01$ at $z\sim 3$ (e.g., Shapley et al. 2006), this suggests that the escape fraction should vary with redshift. In addition, the average number of recombinations per hydrogen atom is $\lt 10$ because the escape fraction cannot exceed unity by definition. While these are only basic constraints, they already help to narrow down the parameter space of reionization. In an upcoming paper for the SCORCH project, we will make more robust constraints on the radiation escape fraction and hydrogen clumping factor.

In large-scale reionization simulations with box sizes $\gtrsim 50\;{\rm{Mpc}}\;{h}^{-1},$ radiation sources are often identified with dark matter halos from N-body simulations, with properties modeled using simple scaling relations. The source luminosity L is assumed to be proportional to the halo mass M through the light-to-mass ratio, which can be constant (e.g., Iliev et al. 2006), have mass dependence (McQuinn et al. 2007), or have both mass and redshift dependence (e.g., Trac & Cen 2007). However, in all three cases the luminosity is deterministic, without any scatter to account for the episodic nature of starbursting galaxies. In an upcoming paper for the SCORCH project, we will use a physically motivated and observationally constrained approach for modeling galaxy formation in large-scale reionization simulations. Using the luminosity–accretion rate relation, we can study the effects of episodic star formation on the distribution and morphology of H ii regions.

SCORCH is a new project to study the EoR and provide useful theoretical tools and predictions to facilitate more accurate and efficient comparison between observations and theory. In this first paper, we probe the connection between observed high-redshift galaxies and simulated dark matter halos in order to better understand the distribution and evolution of the primary source of ionizing radiation. A series of 22 high-resolution N-body simulations shown in Table 1 is used to quantify the abundance of dark matter halos as a function of mass M, accretion rate $\dot{M},$ and redshift z. The abundance matching technique is used to connect the distribution of observed high-redshift galaxies to the distribution of simulated dark matter halos. The major results are as follows.

• 1.  
  The halo mass function dn/dM can be calculated using a self-similar barrier-crossing distribution function $f(\sigma )$ given by Equation (8). The new fit is $\approx 20\%$ more accurate at the high-mass end where bright galaxies are expected to reside.
• 2.  
  The distribution of mass accretion rate at any given mass is positively skewed. The average accretion rate $\langle \dot{M}\rangle ,$ variance ${\mu }_{2},$ and skewness ${\mu }_{3}$ as functions of mass and redshift are quantified by Equations (11)–(13).
• 3.  
  The halo accretion rate function ${dn}/d\dot{M}$ is related to the halo mass function through a mediating mass relation, and can be calculated using Equations (10) and (14). Both halo abundance functions are similar in shape, with each having a low-end power law and a high-end exponential decline.
• 4.  
  The luminosity–accretion relation ${L}_{{\rm{UV}}}(\dot{M},z)$ is more consistent with no redshift evolution than the luminosity–mass relation ${L}_{{\rm{UV}}}(M,z)$ for the redshift range $6\lesssim z\lesssim 10.$ The former is quantified with a universal EoR template given by Equations (23) and (24).
• 5.  
  The star formation rate ${\dot{M}}_{{\rm{s}}}$ evolves more like ${(1+z)}^{2.5}$ rather than the commonly assumed ${(1+z)}^{1.5}$ scaling. More precise measurements of the galaxy luminosity function and a more rigorous statistical analysis are required to strengthen this argument.
• 6.  
  The star formation efficiency is not monotonic with mass nor accretion rate, but reaches a maximum value at a characteristic mass $\sim 2\times {10}^{11}\ {M}_{\odot }$ and a characteristic accretion rate $\sim 6\times {10}^{2}\ {M}_{\odot }\;{{\rm{yr}}}^{-1}$ at $z\approx 6.$ A corresponding template for the efficiency as a function of accretion rate is given by Equation (32).
• 7.  
  The faintest magnitude corresponding to the star formation limit is ${M}_{{\rm{UV,SF}}}\approx -10$ at z = 6, but it becomes more negative (brighter) with increasing redshift. The redshift dependence is given by Equation (31).
• 8.  
  The fiducial model for the galaxy luminosity function is constructed from the halo accretion rate function and the luminosity–accretion-rate relation using Equation (33). The evolution of the Schechter parameters are shown in Table 2: ${\phi }_{\star }$ decreases, ${M}_{\star }$ is more positive (fainter), and α is more negative (steeper) at higher redshifts.
• 9.  
  Forecast galaxy counts for JWST are shown in Table 3. JWST has the sensitivity to observe $\gtrsim 11$ (24) unlensed galaxies per square degree per unit redshift at least down to ${M}_{\star }$ at $z\lesssim 13$ (14) with deep (ultra deep) observations. It will also be able to probe some portion of the faint end at lower redshifts.

The numerical fits for the abundance matching results: luminosity–accretion-rate relation, star formation efficiency, Schechter luminosity function parameters, and JWST forecasts are based on the current galaxy luminosity function at $z\approx 6.$ When updated observations are available, we suggest repeating the abundance matching procedure using the new galaxy luminosity function and the halo accretion rate function.

We thank Nick Battaglia, Rychard Bouwens, Camila Correa, Jeremy Tinker, and Rogier Windhorst for helpful discussions. We thank Rick Costa, Roberto Gomez, and Bill Wichser for help with computing. H.T. acknowledges support from NASA grant ATP-NNX14AB57G and NSF grant AST-1312991. R.C. acknowledges support from NASA grant NNX12AF91G. Simulations were run at the NASA Advanced Supercomputing (NAS) Division, Pittsburgh Supercomputing Center (PSC), and the Princeton Institute for Computational Science and Engineering (PICSciE).