PROBING THE ULTRAVIOLET LUMINOSITY FUNCTION OF THE EARLIEST GALAXIES WITH THE RENAISSANCE SIMULATIONS

The most fundamental characteristics of the earliest galaxies are challenging to determine directly. These galaxies lie at the edge of observability, or beyond, for even the largest current ground- and space-based telescopes, and inferences based on local stellar populations (from, e.g., the Milky Way's stellar halo) are uncertain at best. In particular, it is difficult to determine the characteristics of the luminosity function of high-redshift galaxies, which directly influences how reionization proceeds in the early universe. A recent example of this are efforts to use the Hubble Ultra Deep Field and its extensions (Beckwith et al. 2006; Koekemoer et al. 2013; Finkelstein et al. 2014; Bouwens et al. 2015), and, separately, lensed observations of high-redshift galaxies from the Hubble Frontier Fields (Atek et al. 2015; Ishigaki et al. 2015), to estimate the luminosity function for galaxies at $z\;\simeq \;$ 6–10. These observations succeeded in directly measuring the stellar luminosity function for galaxies with luminosities within a few magnitudes of L^*, but no deeper.

The work of Bouwens et al. (2015) and Atek et al. (2015) provides only weak constraints on the faint-end slope of the luminosity function, which has profound theoretical implications (e.g., Robertson et al. 2010; Fontanot et al. 2014). If the slope of the luminosity function is very steep, this implies a vast number of dim galaxies, unobservable by current instruments, that can produce more than enough ionizing photons to complete and sustain reionization with no additional sources (such as high-redshift black hole populations) needed. If, on the other hand, the slope of the high-redshift luminosity function is shallow, there will be fewer galaxies than one might expect, with too few ionizing photons and thus the need for additional sources of ionizing radiation.

A second question related to the shape of the high-redshift luminosity function is that of its ending point: how bright are the smallest high-redshift galaxies? Simulations of Population III stars and the transition to metal-enriched star formation suggest that the smallest halos to form stars have masses of around 10⁶–10⁷ ${M}_{\odot }$ (O'Shea & Norman 2007; Johnson et al. 2008; Smith et al. 2009; Wise et al. 2012b; Crosby et al. 2013), with correspondingly few stars and low luminosities. These simulations' predictive capabilities, however, suffer from challenges relating to small-number and small-volume statistics.

An additional issue relating to the study of high-redshift galaxies relates to cosmic variance. The James Webb Space Telescope (JWST) will certainly be able to see substantial numbers of galaxies at z ≳ 10, but has a small field of view (FOV; Gardner et al. 2006). As a result, interpretation of any JWST survey must by necessity take into account cosmic variance (e.g., Trenti & Stiavelli 2008; Moster et al. 2011). While this can be done by using many small, widely spaced fields instead of one larger field (Sibthorpe et al. 2013), this may be undesirable, and thus a deeper understanding based on theory may be necessary to correctly interpret high-z survey results.

In this paper, we address these pressing issues regarding the high-redshift galaxy luminosity function using the Renaissance Simulations—a trio of physics-rich simulations of high-redshift galaxy formation that resolve several hundred galaxies apiece, and which explore a range of cosmic environments at $z\geqslant 8$. Using these simulations (described in Section 2), we show results (in Section 3) that predict the high-redshift galaxy luminosity function will flatten at magnitudes that can be probed by the next generation of telescopes. We discuss the implications of this work in Section 4.

The simulations were carried out using Enzo ⁶ (Bryan et al. 2014), an open-source adaptive mesh refinement code that has been extensively used for simulating cosmological structures, and in particular high-redshift structure formation (e.g., Abel et al. 2002; O'Shea & Norman 2007; Turk et al. 2009; Wise et al. 2012a, 2012b, 2014; Xu et al. 2013). Notably, H₂-photodissociating and ionizing radiation from stellar populations is followed using the Moray radiation transport solver (Wise & Abel 2011). The properties of hydrogen and helium are calculated using a nine-species primordial non-equilibrium chemistry and cooling network (Abel et al. 1997), supplemented by metal-dependent cooling tables (Smith et al. 2009). Prescriptions for Population III and metal-enriched star formation and feedback are employed, using the same density and metallicity criteria as Wise et al. (2014) but with a Population III characteristic mass of $40{\text{}}{M}_{\odot }$.

We simulated a region of the universe 28.4 Mpch⁻¹ on a side using the WMAP7 best-fit cosmology. Initial conditions were generated at z = 99 using MUSIC (Hahn & Abel 2011), and a low-resolution simulation was run to z = 6 to find regions suitable for re-simulation. The volume was then smoothed on a physical scale of 5 comoving Mpc, and regions of high ($\langle \delta \rangle \equiv \langle \rho \rangle /({{\rm{\Omega }}}_{M}{\rho }_{C})-1\simeq 0.68$), average ($\langle \delta \rangle \simeq 0.09$), and low ($\langle \delta \rangle \simeq -0.26$) mean density were chosen for re-simulation. These subvolumes, designated the "Rare peak," "Normal," and "Void" regions, with comoving volumes of 133.6, 220.5, and 220.5 Mpc³, were resimulated with an effective initial resolution of 4096³ grid cells and particles in the region of interest, giving a dark matter mass resolution of $2.9\times {10}^{4}$ ${M}_{\odot }$. We allowed further refinement based on baryon or dark matter overdensity for up to 12 total levels of refinement (i.e., a maximum comoving resolution of 19 pc). For more simulation details, see Xu et al. (2013, 2014), Chen et al. (2014). These simulations were evolved to $z=(15,12.5,\mathrm{and}\;8)$ for the Rare peak, Normal, and Void simulations. In these simulations, the halo mass function is well-resolved down to $\simeq 2\times {10}^{6}$ ${M}_{\odot }$ (≃70 particles/halo), and the simulations contained $(822,758,458)$ galaxies having at least 1000 particles $({M}_{\mathrm{vir}}\simeq 2.9\times {10}^{7}$ ${M}_{\odot }$). The simulations were analyzed using the yt analysis tool (Turk et al. 2011).

Figure 1 shows the ultraviolet (UV) luminosity function for the galaxies within the refined region of our simulations. This luminosity function is calculated by extracting the mass, formation time, and metallicity of all star particles in each simulated galaxy at the final timestep of each simulation, and using the stellar population synthesis models of Bruzual & Charlot (2003) to create a spectrum for that galaxy (with each star particle representing a star cluster with a Chabrier mass function). The UV luminosity in a window of ${\rm{\Delta }}\lambda =100$ Å centered on λ = 1600 Å is extracted from the resulting spectra. In each panel, the redshift-dependent fit to the Schechter function from Bouwens et al. (2015) is shown at the same redshift as the simulation data, with the shaded region representing the 1–σ uncertainty in their fit parameters. The figure shows the luminosity function derived for all three simulations at z = 15, the Normal and Void simulations at z = 12, and the redshift evolution of the Normal simulation from z = 18 to z = 12 (212–375 Myr). In each panel, we have calculated the magnitude of the faintest galaxy that could in principle be detected at that redshift in an un-lensed Hubble Frontier Field, the Hubble Ultra-deep Field, and the planned JWST ultradeep field, as well as the limiting magnitude for a JWST observation assuming a field with a magnification factor of $\mu =10$. For all observational estimates we assume a uniform K-correction of −2 mag. In all simulations the UV luminosity function matches observations reasonably well at the bright end, but there is a flattening of the luminosity function in dim galaxies–galaxies that currently cannot be observed, but which in the future may be probed by JWST. We also see that the UV luminosity function has a termination point, with no galaxies dimmer than ${M}_{\mathrm{UV}}\simeq -2$. This behavior is captured reasonably well by a simple toy model, shown in the bottom panel and described in Section 4.

Zoom In Zoom Out Reset image size

Figure 2 displays the relationship between a given galaxy's stellar mass and the halo's virial mass or UV luminosity, combined at their final redshift. In the left panel, a 2D histogram of stellar mass versus halo virial mass, lines of constant stellar mass fraction (defined as ${f}_{*}\;\equiv$ M_*/M_vir), ranging from ${f}_{*}={10}^{-1}$ to 10⁻⁴, are plotted. The purple solid line and accompanying equation is a fit to the median stellar fraction as a function of virial mass,

$$\begin{eqnarray}&&{f}_{*}=1.26\times {10}^{-3}{\left(\displaystyle \frac{{\text{}}{M}_{\mathrm{vir}}}{{10}^{8}\;{\text{}}{M}_{\odot }}\right)}^{0.74},\end{eqnarray} \tag{ 1 }$$

which is made by fitting a line to the distribution of log stellar mass and log virial mass for halos with ${10}^{7}\lt$ M_vir/M${}_{\odot }\lt {10}^{8.5}$. At virial masses larger than ${10}^{8.5}$ ${M}_{\odot }$, there is some deviation from this relationship, with the stellar fraction in the most massive halos being lower than this power-law fit would predict. The red dashed line indicates the stellar fraction as a function of virial mass predicted by Behroozi et al. (2013) but extrapolated to z = 15—a fit that agrees with the simulated galaxies to a remarkable degree, including the flattening of the stellar fraction at high mass. The right panel shows a 2D histogram of halo virial mass versus absolute UV magnitude for the same galaxies, and indicates a positive (albeit noisy) correlation between halo mass and star formation rate (SFR). The quantities in both panels are insensitive to redshift during the epochs considered by this work.

Zoom In Zoom Out Reset image size

Figure 3 shows the fraction of halos containing metal-enriched stellar populations as a function of halo virial mass for all three simulations at their final redshift. All halos with ${M}_{\mathrm{vir}}\gtrsim 5\times {10}^{7}$ ${M}_{\odot }$ contain stars, and all halos with ${M}_{\mathrm{vir}}\gtrsim 2\times {10}^{8}$ ${M}_{\odot }$ have formed stars in the last 20 Myr (and thus are emitting significant UV light). No halos with virial masses below $\simeq 7\times {10}^{6}$ ${M}_{\odot }$ contain stars. This figure also contains a fitting function for the fraction of halos that have active stars (modeled on the form used by Okamoto et al. 2008):

$$\begin{eqnarray}&&{f}_{\mathrm{occ}}(M)={\left[1+\left({2}^{\alpha /3}-1\right){\left(\displaystyle \frac{M}{{M}_{c}}\right)}^{-\alpha }\right]}^{-3/\alpha }\end{eqnarray} \tag{ 2 }$$

with ${M}_{c}=6.0\times {10}^{7}$ ${M}_{\odot }$ and $\alpha =1.5$.

Zoom In Zoom Out Reset image size

Figure 4 shows projections of matter overdensity and electron fraction in the refined simulation subvolumes at their final redshifts. In each panel of the top row, the black bar indicates the physical size subtended by a FOV $0\buildrel{\,\prime}\over{.} 5$ on a side at that redshift. The projection area is slightly smaller than the $2\buildrel{\,\prime}\over{.} 2\times 2\buildrel{\,\prime}\over{.} 2$ JWST NIRCam FOV. In both rows, all galaxies with ${M}_{1600}\lt -10$ are shown, with glyphs showing galaxies which could be seen in the 2012 Hubble Ultra Deep Field (red triangles; Koekemoer et al. 2013), the un-lensed (blue squares) and lensed (green circles; μ = 10) JWST ultradeep field (Gardner et al. 2006), and dimmer galaxies (black dots). The spatial distribution of galaxies in each field demonstrates cosmic variance. The rare peak field shows many clustered galaxies above the lensed JWST limit, whereas the normal field has very few galaxies outside of a clustered region. The void region is bounded on the right by an overdense sheet containing several galaxies above the JWST unlensed sensitivity limit. In the void, there is little collapsed structure that can sustain efficient star formation, leading to a lack of galaxies with ${M}_{1600}\lt -10$. The projected electron fraction in the lower row depicts the H ii regions with radii up to ∼1 comoving Mpc associated with these galaxies. Scattered around the volume are relic H ii regions created by Population III stars that form in substantially smaller halos (${M}_{\mathrm{vir}}\sim {10}^{6}{\text{}}{M}_{\odot }$) than the galaxies shown in the figure, marginally adding to the overall ionized fraction. The overall ionized morphology is consistent with an inside-out scenario, previously noted by several groups (e.g., Gnedin 2000; Trac et al. 2008; Finlator et al. 2009; So et al. 2014) during the initial phases of reionization. This is most apparent in the void simulation, because the ionization front from the overdense sheet has not yet propagated into that region and there are no void galaxies to reionize their surroundings.

Zoom In Zoom Out Reset image size

The key results presented in this paper are as follows.

• 1.  
  The UV luminosity function of our simulated galaxies is consistent with observations of high-z galaxy populations at the bright end of the luminosity function (${M}_{1600}\lesssim -17$), but at lower luminosities is essentially flat rather than rising steeply as inferred by a Schechter function fit to observational data, and has a clearly-defined lower limit in UV luminosity.
• 2.  
  This behavior of the UV luminosity function is due to two factors: (i) the strong dependence of the SFR on halo virial mass in our simulated galaxy population, with lower-mass halos having systematically lower SFRs and thus lower UV luminosities; and; and (ii) the fact that halos with virial masses below $\simeq 2\times {10}^{8}$ ${M}_{\odot }$ do not universally contain stars, with the fraction of halos containing stars falling with decreasing virial mass and reaching zero at $\simeq 7\times {10}^{6}$ ${M}_{\odot }$.
• 3.  
  The brightest of our simulated galaxies may be visible to current and future ultra-deep space-based surveys at $z\sim 12$ (and likely at lower redshifts as well), particularly if lensed regions are chosen for observation.

The primary result of this paper—the flattening of the mass function—is significant for our understanding of the reionization of the universe. Observations of high-redshift galaxies provide poor constraints on the low-luminosity end of the galaxy luminosity function, and thus make it challenging to accurately account for the full budget of ionizing photons during that epoch. Our work suggests that there are far fewer faint galaxies at high redshift than one would infer from fitting a Schechter function to the high-luminosity end, as has been done by, e.g., Bouwens et al. (2015). Taken at face value, a smaller number of galaxies would result in an under-production of photons relative to what is needed to reionize the universe. However, the Planck 2015 results (Planck Collaboration et al. 2015), particularly the Thomson scattering optical depth $\tau =0.066\pm 0.012$ corresponding to an instantaneous reionization redshift ${z}_{r}={8.8}_{-1.1}^{+1.2}$, relieve the previous tension between reionization constraints from polarization measurements of the CMB and galaxy observations (Robertson et al. 2015).

There are many complicating factors when computing a reionization history from a galaxy luminosity function—in particular, the UV photon escape fraction and its relationship to the halo virial mass (see, e.g., Fontanot et al. 2014; Wise et al. 2014)—and we defer a detailed analysis to future work (H. Xu et al. 2015, in preparation). Recent work has shown that low-mass galaxies (${M}_{\mathrm{vir}}\lesssim {10}^{8}{\text{}}{M}_{\odot }$; ${M}_{\star }\lesssim {10}^{5}{\text{}}{M}_{\odot }$) prior to reionization have typical escape fractions in the range 5%–50% with the upper range occurring when supernova blastwaves and ionization fronts create low-density channels that allow ionizing radiation to freely escape into the neutral intergalactic medium (Paardekooper et al. 2013, 2015; Wise et al. 2014). In more massive halos, the simulations of Kimm & Cen (2014) showed that the time-averaged UV escape fraction for a particular galaxy is only ∼10% (although Ma et al. 2015 suggest that this fraction may be considerably lower over halos with a mass range of 10⁹–10¹¹ ${M}_{\odot }$, with a time-averaged value of around $\simeq 5\%$). This suggests that low-luminosity galaxies may provide a non-negligible amount of the ionizing photon budget, given their higher escape fractions and number densities.

The flattening of the high-redshift galaxy luminosity function at the faint end is a robust result (see Appendix B) in our simulations (and has been suggested by other, albeit less physics-rich, simulations—see Jaacks et al. 2013). What is less robust, however, is what can be inferred from the overall normalization of the luminosity function from our calculations. Given JWST's small FOV, cosmic variance will present a significant challenge (e.g., Trenti & Stiavelli 2008). We have sampled regions of varying mean density, which gives some sense of cosmic variance. However, a precise discussion of this subject is difficult because of the small cosmological volumes sampled, as well as the somewhat arbitrary nature of our choice of resimulated regions.

We note that the strongly evolving stellar mass fraction as a function of halo mass has an intriguing similarity with observations of Local Group dwarf galaxies (e.g., McConnachie 2012), where there is a strong positive correlation between a given galaxy's mass-to-light ratio and its bolometric luminosity. The smallest local dwarf galaxies are strongly dark matter-dominated, having similar stellar fractions to our lower-mass galaxies, and the trends are similar as well. While it is hard to make precise comparisons between Local Group observations and high-redshift galaxy simulations, the similarity in trends suggests that there may be universal, halo mass-dependent behaviors in galaxy formation.

Figure 3 indicates something intriguing—namely, that there is some physical effect that quenches star formation in low mass halos in the early universe. The cause is not reionization—Figure 4 demonstrates that the universe is patchily reionized at this epoch—but rather H₂-photodissociating Lyman–Werner radiation (e.g., O'Shea & Norman 2008). The universe is transparent to this radiation, and thus star formation anywhere in the volume can suppress molecular hydrogen elsewhere, with the level of suppression growing with the comoving SFR density. The consequence of this suppression is that low-mass halos of primordial composition, which cool primarily by H₂ line emission, must grow to be more massive and hotter before the gas in the halo core can cool efficiently, which suppresses Population III star formation (and, thus, metal-enriched star formation). Once halos reach this threshold mass (which depends on the strengh of the Lyman–Werner background; see O'Shea & Norman 2008), stars can then form. We note that at later times, when the universe is transparent to hydrogen-ionizing photons, a similar but more pronounced effect must exist due to this photon population.

Figures 2 and 3 strongly suggest that the flattening of the galaxy luminosity function is due to a combination of decreased SFR at low halo virial masses and the lack of stars in halos below $\sim {10}^{8}$ ${M}_{\odot }$. A simple toy model (see Appendix A) using the analytic halo mass function convolved with a constant specific star formation rate (sSFR) and the halo stellar occupation fraction from Figure 3 is shown as the magenta line in the bottom panel of Figure 1, and matches the properties of the luminosity function of our simulated galaxies quite well, accurately capturing the deviation from a Schechter function, the overall flattening at intermediate luminosity, and the termination at the low-luminosity end.

The dependence of stellar fraction on halo mass is well-known observationally (e.g., Behroozi et al. 2013), and is also seen in simulations of more massive galaxies at lower redshifts (e.g., Hopkins et al. 2014). We speculate that this behavior is due to a combination of cooling and stellar feedback (as has been discussed in a different context by Voit et al. 2014). Gas in low-mass halos cools inefficiently due to the low virial temperature, and stellar feedack can efficiently remove metal-enriched gas from these low-mass halos (see, e.g., Chen et al. 2014), which further decreases the efficiency with which gas can cool and form stars. The suppression of star formation in halos below $\sim {10}^{8}$ ${M}_{\odot }$ (and complete absence below $\sim {10}^{7}$ ${M}_{\odot }$) is likely due to inefficient cooling, destruction of H₂ via Lyman–Werner radiation from neighboring halos, and the resulting Jeans screening (e.g., Gnedin & Hui 1998; Wise & Abel 2008; Wise et al. 2014). This suggests that there is an effective lower limit on the masses of high-redshift galaxies, and thus a definitive lower end to the UV luminosity function. We will explore this in future work.

This research is part of the Blue Waters project using NSF PRAC OCI-0832662. This research was supported by NSF and NASA grants PHY-0941373, AST-1109243, AST-1211626, AST-1333360, NNX12AC98G, HST-AR-13261.01 A, HST-AR-13895.001. The authors thank Devin Silvia and Eric Bell for useful discussions.

The simple luminosity function model shown by the magenta line in the bottom panel of Figure 1 was created using the following procedure.

• 1.  
  Calculate the comoving number density of cosmological halos as a function of halo virial mass m, $\frac{{dN}(m)}{{dV}}$, at z = 12 for the WMAP7 best-fit cosmology, using the fitting function of Warren et al. (2006).
• 2.  
  For each halo mass bin, calculate the recent star formation history of the halo assuming a constant sSFR of 10⁻⁸ years⁻¹ over the last 100 Myr. This results in a constant SFR, under the assumption of constant halo mass.
• 3.  
  Given this star formation history and the assumption of a simple Salpeter IMF with parameters identical to those in the z = 0 Milky Way (consistent with our simulation treatment of metal-enriched stars; see Section 2), calculate the IMF-integrated bolometric luminosity of the stars that have not yet undergone supernovae at the epoch under consideration. Given that this luminosity is dominated by the most massive and hottest stars (since luminosity scales approximately as stellar mass to the fourth power), we make the assumption that the UV luminosity of the stellar population in each halo is roughly its bolometric luminosity. This results in a comoving number density of halos as a function of UV magnitude, $\frac{{dN}}{{dVdM}}$, where M represents absolute UV magnitude.
• 4.  
  Finally, we multiply the resulting quantity by the fraction of halos that have experienced recent star formation as a function of halo mass that is shown in Figure 3, using the crude approximation shown by the gray line. This fit captures the essence of the figure—namely, universal recent star formation above M${}_{\mathrm{vir}}\geqslant 2\times {10}^{8}$ ${M}_{\odot }$, no recent star formation below $\sim {10}^{7}$ ${M}_{\odot }$, and a steady decrease in star formation with decreasing halo mass between those endpoints.

This toy model captures several of the essential features of the luminosity function shown by our simulations: it has a significant slope at the high-luminosity end, flattens out at ${M}_{1600}\gt -12$, and goes entirely to zero at ${M}_{1600}\simeq -2$. Furthermore, the properties of this toy model are robust to several of the simplifying assumptions made here. The most critical of these are as follows.

Specific star formation rate: the sSFR evolves as both a function of redshift and halo mass. This evolution is most pronounced at low redshift ($z\lt 2$)—at higher redshifts, such as those considered in this work, the sSFR evolves very slowly with redshift and is relatively insensitive to halo mass (see, e.g., Appendix F of Behroozi et al. 2013), making the approximation of a constant sSFR sensible. We note that at the redshifts and halo masses of interest to the work presented here this quantity is extremely hard to constrain observationally; however, extrapolating the observed trends to $z\gt 8$ suggests little evolution. Deviation from a constant sSFR will modify the predicted UV luminosity function in proportion to the deviation.

Halo formation history: cosmological halos grow over time, and do so particularly quickly at the redshifts considered in this work. Assuming that the sSFR is directly proportional to halo mass, this implies that the total (as opposed to specific) SFR in a given halo is increasing rapidly as well. The period of 100 Myr over which the SFR is assumed to be constant is a significant fraction of the Hubble time at high redshift ($\simeq 250$ Myr at z = 16, $\simeq 380$ Myr at z = 12), suggesting that this assumption may be suspect. However, the UV luminosity is dominated by massive stars having lifetimes of millions of years, and thus the production of UV photons is dominated almost entirely by the very recent star formation history of the halo. As a result, the UV luminosity is quite robust to the exact details of the halo's formation history, with deviations from constant SFR contributing mildly to variation in magnitude.

The simulated behavior of the UV luminosity function at the faint-end can be suspectible to limited numerical resolution. To check for any spurious effects, we compare our results to our previous work (Wise et al. 2012a, 2012b, 2014) that focused on 32 first-generation galaxies in a small comoving volume of (1 Mpc)³ with a DM particle mass resolution of 1840 ${M}_{\odot }$, 15 times smaller than the dark matter particles in the simulations presented in this paper. The major shortcoming of our previous work was the limited galaxy sample size, which was rectified in the Renaissance Simulations. These two simulation suites share the same subgrid physics models—chemical reaction networks, metal-free and metal-enriched star formation, and radiative and supernova feedback. Although the Renaissance Simulations cannot capture Population III star formation in the smallest minihalos ($\sim 2\times {10}^{5}{\text{}}{M}_{\odot }$), properties of metal-enriched star formation, in particular the stellar mass–halo mass relation, agree very well between the two simulation suites (see Figure 4 in Chen et al. 2014). There we showed that this relationship in both simulations were within 1σ of each other in halos with ${M}_{\mathrm{vir}}\geqslant {10}^{7}{\text{}}{M}_{\odot }$, which are the smallest halos considered in our analysis. Furthermore, Wise et al. (2014) found that the UV luminosity function was flat at magnitudes ${M}_{1600}\gtrsim -12$, consistent with our findings in the Renaissance Simulations. Thus, we are confident that the flattening of the luminosity function is not a result of limited numerical resolution but is caused by the suppression of star formation from radiative and supernova feedback effects in halos with masses ${M}_{\mathrm{vir}}\lesssim 2\times {10}^{8}{\text{}}{M}_{\odot }$.