Steep Faint-end Slopes of Galaxy Mass and Luminosity Functions at and the Implications for Reionisation
We present the results of a numerical study comparing photometric and physical properties of simulated galaxies to the observations taken by the WFC3 instrument aboard the Hubble Space Telescope. Using cosmological hydrodynamical simulations we find good agreement with observations in color-color space at all studied redshifts. We also find good agreement between observations and our Schechter luminosity function fit in the observable range, , provided that a moderate dust extinction effect exists for massive galaxies. However beyond what currently can be observed, simulations predict a very large number of low-mass galaxies and evolving steep faint-end slopes from at to at , with a dependence of . During the same epoch, the normalization increases and the characteristic magnitude becomes moderately brighter with decreasing redshift. We find similar trends for galaxy stellar mass function with evolving low-mass end slope from at to at , with a dependence of . Together with our recent result on the high escape fraction of ionizing photons for low-mass galaxies, our results suggest that the low-mass galaxies are important contributor of ionizing photons for the reionisation of the Universe at .
keywords:method : numerical — galaxies : evolution — galaxies : formation - cosmology : theory
When and how galaxies formed throughout the history of the Universe is one of the most fundamental questions of astronomy and astrophysics. As technology improves, astronomers are able to push the redshift frontier of galaxy observation and attempt to obtain a glimpse of galaxy formation processes at high redshift. In particular, upgrades to the Hubble Space Telescope (HST) via the addition of the Wide Field Camera 3 (WFC3) have expanded our view of the Universe to reach redshifts , and dozens of new high-redshift candidates have been identified (e.g., Yan et al., 2009; Oesch et al., 2009; Ouchi et al., 2009; Bouwens et al., 2010; Wilkins et al., 2011). For example, Bouwens et al. (2010) identified Lyman-Break candidates by color-color selection with AB magnitudes and faint-end luminosity slopes of at and at .
While these new observational results are very exciting, candidate galaxies at are very faint and it is still difficult to obtain high-resolution spectra to estimate their physical properties such as stellar masses and star formation rates (SFR). Therefore it would be useful to obtain predictions from theoretical models on their spectro-photometric properties. In particular, cosmological hydrodynamic simulations based on the standard concordance cold dark matter (CDM) model can directly simulate the dynamical evolution of dark matter and baryons from the initial condition that is consistent with the observations of cosmic microwave background. With the implementation of radiative heating and cooling, one can also follow the thermal properties of intergalactic medium and collapse of gas clouds into galaxies. Much work has been done for comparing simulated galaxies at to observations (e.g., Nagamine et al., 2004, 2005a, 2005b; Night et al., 2006; Finlator et al., 2006; Robertson et al., 2007), but the comparisons of photometric properties of galaxies at between simulations and observations have not been performed adequately yet.
|()||(Dark Matter, Gas)||()||()||()|
One of the important missing information from theoretical predictions is the precise determinations of evolution of galaxy luminosity function (LF) and stellar mass function (GSMF) at . Some of the earlier numerical work have suggested very steep faint-end slope of at (Night et al., 2006; Finlator et al., 2006), but its evolution has not been quantified yet at using cosmological hydrodynamic simulations. This is an important topic, because the number density of early dwarf galaxies has significant implications for the reionisation of the Universe at . It is thought that small dwarf galaxies, that are currently beyond detection, may have played an important role in the reionisation of the Universe (e.g., Muñoz & Loeb, 2010; Trenti et al., 2010; Yajima et al., 2010).
In this paper, we compare the results of our cosmological simulations with the recent WFC3 observations. In particular, we will determine the evolution of faint-end (low-mass end) slope () at more precisely in our simulations. The paper is organized as follows: Section 2 will briefly describe the simulation and setup that was used to obtain our data; Section 3 will outline the methods used in processing the photometric data; Sections 4-6 will detail assembled data from simulations; Section 7 will discuss the implications of these findings to the epoch of reionisation. Finally in Section 9 we will summarise our findings.
We use a modified version of the smoothed particle hydrodynamics (SPH) code GADGET-3 (originally described in Springel, 2005). This modified code includes radiative cooling by H, He, and metals (Choi & Nagamine, 2009), heating by a uniform UV background of a modified Haardt & Madau (1996) spectrum (Katz et al., 1996; Davé et al., 1999), an Eisenstein & Hu (1999) initial power spectrum, star formation via the ”Pressure model” (Schaye & Dalla Vecchia, 2008; Choi & Nagamine, 2010), supernova feedback, sub-resolution model of multiphase ISM (Springel & Hernquist, 2003), and the multicomponent variable velocity (MVV) wind model (Choi & Nagamine, 2011). We also explore the effect of a new UV background prescription by Faucher-Giguère et al. (2009) in Section 8.2. Our current simulations do not include the AGN feedback.
We have shown earlier that the metal line cooling enhances star formation across all redshifts by about % (Choi & Nagamine, 2009), and that the Pressure SF model suppresses star formation at high-redshift due to a higher threshold density for SF (Choi & Nagamine, 2010) with respect to the earlier model by Springel & Hernquist (2003). Choi & Nagamine (2011) also showed that the MVV wind model, which is based on the momentum-driven wind, makes the faint-end slope of GSMF slightly shallower compared to the constant velocity galactic wind model of Springel & Hernquist (2003). One of the points of this paper is to report that, even with these improvements to the supernova feedback models, our simulations continue to predict very steep faint-end slopes.
It is known that the collisional ionization equilibrium (CIE) table of Sutherland & Dopita (1993) overestimates the cooling rates compared to the case when the photoionisation by the UVB is taken into account (Efstathiou, 1992; Wiersma et al., 2009), and Choi & Nagamine (2009) have used this CIE table. We have compared our LF results in the two runs with and without the metal cooling, and found that there is little difference in the faint-end slope of the two simulations. Therefore this issue is not a problem for the present work.
For this paper, simulations are setup with either or particles for both gas and dark matter. Multiple runs were made with comoving box sizes of Mpc, Mpc and Mpc. These runs will be referred to as N400L10, N400L34 and N600L100, respectively. Simulation parameters are summarised in Table 1. The adopted cosmological parameters are based on the most recent WMAP data: , , , , , (Komatsu et al., 2010).
At each time step of the simulation, the gas particles that exceed the SF threshold density ( cm) is allowed to spawn star particles with the probability consistent with the intended SFR. Choi & Nagamine (2010) have shown that the Pressure SF model with the above produces favorable results compared to the observed Kennicutt (1998) law, in particular at the low column density end. Galaxies are then identified by a group-finder algorithm using a simplified variant of the SUBFIND algorithm (Springel et al., 2001). Collections of star and gas particles are grouped based on the baryonic density field. Properties such as stellar mass, SFR, metallicity and position, among others are recorded in galaxy property files. A more detailed description of the group-finder process can be found in Nagamine et al. (2004).
Once we identify the individual galaxies, we calculate the spectral energy distribution (SED) of each galaxy by applying the 2007 version of Bruzual & Charlot (2003) stellar population synthesis model that includes the contribution of TP-AGB phase. We have repeated the calculation using the 2003 version of Bruzual & Charlot (2003), and confirmed that the results hardly change. This is expected because we are focusing our analyses on the rest-frame UV wavelengths range at , and the effect of TP-AGB appears only in the infra-red range. In the present work, nebular emission lines are not included in the SEDs, but our results are not affected because this work is mainly concerned about the comparisons in the rest-frame UV, where the optical emission lines do not come into the filters. We will investigate the impact of nebular emission lines on the estimates of stellar masses of high-redshift galaxies in the subsequent work.
After we calculate the intrinsic SEDs, we then account for the dust extinction effect using the Calzetti (1997) extinction law. We adopt uniform constant values of & 0.30. The above extinction values are selected based on the varied data found in recent publications. Bouwens et al. (2009) claimed a trend of decreasing extinction approaching at based on the UV continuum slope. However Schaerer & de Barros (2010) performed SED fits including the effect of nebular emission lines, and showed that dust attenuation of up to ( for ) can be found for galaxies at , and that might be greater for higher mass galaxies. Our final choice for at each redshift is based on a visual match of the LF to the observed data range, staying within the above mentioned range. In Section 5 we also examine the effect of variable extinction models.
The final step to achieving a spectrum as it would appear to HST, is to factor in the absorption by the intergalactic medium (IGM). To do this we apply the model of Madau (1995). Figure 1a shows the relative flux of a representative galaxy redshifted to & 9. The HST filter functions used by the observers are also shown on the X-axis: , , , , centered at , , , , , respectively, are some of the HST filters strategically positioned to identify the Lyman break of high redshift galaxies. Figure 1b shows the effect of Calzetti (1997) extinction on the simulated galaxy spectrum at for & 0.30.
The observer-frame spectrum is then used to calculate the apparent magnitude of each galaxy for each HST filter, as well as the rest-frame UV magnitude centered at 1700Å. This process gives us a rest-frame magnitude in each filter that can then be assembled into color-color plots and luminosity functions (LFs) for direct comparison to those produced by observers.
4 Color-Color Diagrams
In Figures 2, 3 & 4, we present color-color diagrams at , 7 & 8, respectively. For each redshift we make direct comparisons to the observed data. In each figure, the results from three simulations (N400L10, N400L34, N600L100) are shown separately, and the simulated galaxies are represented by blue, green, and red points for , 0.15, & 0.30, respectively. Galaxy evolution tracks are shown by the corresponding solid blue, red and green lines, which were produced by redshifting the average spectrum of simulated galaxies. By doing this we are able to show how the Lyman-break of the galaxy spectrum is passing through the HST filter set (Figure 1a), transitioning into and out of the drop-out candidate selection area which has been defined by observers. For example in Figure 2, for galaxies are substantially redder than at for , creating a threshold above which candidates are selected.
In Figure 2, observed candidate data points are represented by solid black squares (Bouwens et al., 2006) and open black triangles (Yan & Windhorst, 2004). The grey shaded area indicates the color-color selection criteria defined by Bouwens et al. (2006). Due to the abundance of observed points at this epoch, error bars were omitted for clarity.
Figure 3 shows color-color plots for with observed data points from Oesch et al. (2009, black filled squares) and Yan et al. (2009, open triangles). The grey shaded area for color-color selection is taken from Oesch et al. (2009). Error bars were included only for the Oesch et al. (2009) data for clarity. Again we see good agreement with observed data.
Figure 4 shows the color-color diagram for with observed data points from Bouwens et al. (2010, black filled squares) and Yan et al. (2009, open triangles). Expected position for objects is defined by the grey shaded area (Bouwens et al., 2010). Error bars were included only for the Bouwens et al. (2010) data for clarity. Here too simulated galaxies show good agreement with observations.
In all of Figures 2, 3, & 4, massive galaxies are on the lower left side of the data distribution for the simulated galaxies, because higher mass galaxies have higher SFR in our simulations and hence bluer.
In Figure 5 we examine the tight grouping of galaxies in our color-color plots. Figure 5a shows the effect of adding a scatter with a gaussian distribution centered around . This addition does little to increase the scatter of the galaxies in color-color space. In Figure 5b we demonstrate the effect of uncertainties in photometric redshift estimation by adding a scatter with a gaussian distribution centered on . This addition produces a vertical feature which runs throughout the selection box with some horizontal scatter at its top. Finally in Figure 5c we add a scatter with a gaussian distribution to the calculated UV magnitudes as well as including the previous scatter in and redshift. In this case, we reproduce a similar degree of scatter to the observed data. Our results indicate that the current observational color-color selection region can select galaxies with a wide range of extinction values of , and that the scatter in the data points are primarily dominated by the redshift scatter and photometric errors rather than the extinction variance. The tight grouping of galaxies seen in each of the color-color plots is because the simulated galaxies are extracted at exactly and have no errors included in magnitude calculations.
5 Luminosity Functions
The rest-frame UV LFs of simulated galaxies at are shown in Figure 6. To represent a larger range of magnitudes, we construct composite LFs using three simulations with different box sizes, and connect the results at the points where they overlap. At the brightest end, the data correspond to the last bin in which we have data in the N600L100 run. At the faintest end, the minimum corresponds to the minimum allowed by the mass resolution in the N400L10 run (see Section 6). This also corresponds to the point at which the LF begins to turn over, indicating the resolution limit. For example, at , the magnitude coverage of each run are approximately from to mag (N600L100), to mag (N400L34), and to mag (N400L10). These values vary slightly with redshift. In all panels, three sets of results for at & 7, and at & 9 are shown. See Section 3 regarding these extinction choices.
To explore the impact of a varying extinction, we utilize a model presented in Finlator et al. (2006): , where is the galaxy stellar metallicity and is a Gaussian scatter. This relationship is based on the mean trend between metallicity and found in SDSS at . In Figure 6(a) the result of this relationship applied to our simulated galaxies is shown by the cyan solid line. This relationship gives higher extinction for more massive galaxies, therefore our data is consistent with observations only at the very brightest end of the LF with this extinction model.
The disagreement at fainter magnitudes between the cyan line and the yellow observed region has several potential explanations. The Finlator’s model for is based on the observed relationship at , and thus it may not apply at high redshifts. This model also utilizes metallicity taken from star particles, which is set to the same value as the parent gas particle at the time of formation. Therefore the metallicity of star particles in our simulation does not evolve with redshift. We have tested this model utilizing the gas particle metallicity, and found that this change results in lower extinction values, however these values may be artificially lower since this mean metallicity is obtained by averaging over the gas particles including those at the outskirts of galaxies. Finally there is the possibility that if high redshift galaxies are mostly dust free as suggested by recent observations (Bouwens et al., 2009, 2011), then our simulations are systematically overproducing galaxies by dex. This could be true if molecular gas cooling plays an important role in early star formation (see Section 8.2 for further discussion). However the consistent overproduction of galaxies is not a feature found in our GSMF, as we will discuss further in Section 6 and Figure 10.
At (top left panel of Figure 6), the yellow shade represents a range of observations by Yan & Windhorst (2004); Bouwens et al. (2004); Dickinson et al. (2004); Malhotra et al. (2005); Bunker et al. (2004); Beckwith et al. (2006). At this redshift we show good agreement with observed data.
At , the yellow shaded area represents a range of observations by Richard et al. (2006); Mannucci et al. (2007); Bouwens et al. (2008); Oesch et al. (2009); McLure et al. (2010); Oesch et al. (2010). Here too we see good agreement with observations for the N400L34 and N600L100 run.
At , we only show the results from the N400L10 and N400L34 runs as there are no objects found at this redshift in the N600L100 run due to low resolution. Also there has been no published data for robust observed LF estimates at .
5.1 Schechter Function Fit
We perform fit analysis of the Schechter (1976) function to the composite simulated galaxy LFs. The Schechter function, , as a function of magnitude is given by
where , and are the normalization, characteristic magnitude, and faint-end slope.
We minimize the expression
by systematically adjusting each of the three Schechter parameters, where is our simulated data set, is the calculated value of the above Schechter function, and & 10 for & 9, respectively. For a simple power-law was used, because the N600L100 run did not contain any galaxies at this redshift due to limited resolution. The -sigma error of the parameters are determined by fixing two of the best-fit parameters and varying the third from to (Andrae, 2010).
Figure 7 shows our best-fits to the Schechter function along with the observational fits (long-dashed line) by Bouwens et al. (2010). The dotted line represents the -sigma error range of Bouwens et al. (2010) fit. The yellow shaded areas are the same as in Figure 6. For the Schechter fits, we used the simulation data sets with for & 7, and for & 9, as they agree with the yellow shade best. The best-fit Schechter parameters are given in Table 2.
Both our Schechter fit and the simulation LFs agree well with the yellow observed range at the bright end. However our best-fit Schechter LFs are noticeably steeper than the observed one at the faint end. This is caused by numerous faint galaxies with mag in our simulation, which lie beyond the current detection limit. Our results suggest that future surveys with increased sensitivity (such as the JWST) will detect a large number of faint galaxies. Although at different redshifts, this trend of finding more galaxies below the previous magnitude limit seems to be occurring in recent observations. For example, Drory et al. (2009) found an upturn toward a steep faint-end slope of at lower mass at , when they probed down to fainter limits than before in the COSMOS field. Also, Reddy & Steidel (2009) found a steeper slope of at than before using a larger sample of spectroscopic redshifts and 31000 LBGs.
We summarise the evolution of Schechter parameters of simulated galaxies in Figure 8. The yellow shades show the range of observed data points taken from Bouwens et al. (2010). In the top left panel, our simulation shows a clear evolution of from at to at , and our slopes are in general steeper than the current observations by . At , we used a power-law fit to the incomplete data set (missing bright-end in the N600L100 run), therefore the & 9 data points are connected by a dashed line.
Our simulations suggest a steepening trend of
This steepening of towards higher redshift is in congruence with the result of Bouwens et al. (2011), which finds a similar redshift evolution. As we will show in the next section, this steepening of is driven by the simultaneous steepening of galaxy stellar mass function. In a CDM universe, we expect a larger number of lower mass galaxies at higher redshifts, and these low mass galaxies gradually merge to form more massive galaxies. Therefore the steepening of with increasing redshift is naturally expected in our simulations.
The characteristic magnitude in our simulation does not show much evolution, except that it brightens by about 0.3 mag from to . The simulation is brighter than the observed range by mag. The normalization parameter of simulated galaxies increases significantly by about an order of magnitude from to with a redshift dependence of
The simulation is lower than the observed one by about dex.
One of the possible reasons for our lower and brighter than the observational estimates is that our LF does not have a distinct exponential break at the bright-end, and the simulation LF seems more like a single power-law with a slight curvature at the bright-end. Therefore we are somewhat forcing a Schechter fit, which could push the to brighter side to lower side. Future simulations with AGN feedback may make the bright-end break more pronounced, and the values of above two parameters might approach the observed range.
6 Galaxy Stellar Mass Functions
Galaxy stellar mass functions (GSMFs) in our simulations at are shown in Figure 9. As with the luminosity function, data from three simulations have been combined to show a continuous mass function. The yellow shaded area shows the observed mass functions for & 7 (Gonzalez et al., 2010), which was derived directly from the luminosity function obtained by Bouwens et al. (2010). The lower mass limit of for the N400L10 run was determined by converting the faint-end luminosity cut-off of to via the relationship between and in our simulations (Figure 10). In order to stay within the mass resolution limits of each run, we take the following mass range from each run: (N400L10), (N400L34), (N600L100). This ensures that, based on the star particle mass for each run (Table 1), we will have a sufficient number of particles in each galaxies. For example, the N400L10 run has a star particle mass of , which is half of the gas particle mass. Therefore a galaxy of would contain approximately star particles. Simulated galaxies would also contain gas and dark matter particles, so the actual particle number of simulated galaxies would be significantly larger. Further discussion of the simulation resolution limits can be found in Section 8.1.
The Schechter function as a function of galaxy stellar mass is defined as
where , , and are the normalization, characteristic stellar mass, and low-mass (faint-end) slope, respectively. We perform the same analysis for these three Schechter parameters. The best-fit parameters are given in Table 3. Similarly to the LF, we find that the GSMF steepens from at to at , with an evolutionary trend of
Figure 8 includes the results of GSMF Schechter fits, showing a clear trend that galaxies become more massive and numerous from to as they grow with time. The characteristic number density and the associated increase with decreasing redshift, evolving as
These trends are also summarised in the bottom right panel of Figure 9, clearly showing the evolution of GSMF from to .
Interestingly the slope is steeper for the GSMFs than for the LFs by . This can be explained by examining the relationship in our simulations between and (Figure 10). Through a linear fit to the median data point in each bin, we find
where and at , respectively. Since the faint-end slopes of the Schechter luminosity and mass functions (Eqs. [5.1] and ) are and , respectively, the relationship between the two slopes can be derived from as follows:
Therefore plugging in and solving for results in
which clarifies the relationship between the two slopes. Utilizing this relationship, we are able to recover either or to within sigma from one another. The exact value of must depend on the details of the stellar IMF and star formation histories, and we plan to investigate the evolution of as a function of redshift in more detail in future work.
We see in Figure 9 that our simulation results and the observational result of Gonzalez et al. (2010) do not agree well at both faint and bright end, with Gonzalez’s GSMF being shallower than our simulation result. This can be clearly explained through different vs. relationship, as we show in Figure 10. The Gonzalez’s relationship (shown as blue long-dashed line; ; their Fig. 1) is much steeper than what our simulation suggests, and for a given value of , they infer a higher for a brighter , and a lower for a fainter . We note that they derived this relationship from data set, and the limited data at does not seem to follow this relationship tightly. This difference in the slope of vs. relationship is directly reflected in the difference in the GSMF slope. If we translate into luminosity, then our simulation predicts , whereas Gonzalez et al. (2010) obtained at . We plan to study the evolution of relationship between , and SFR in our simulations in our subsequent work.
In the next section, we discuss the implications of these steep mass/luminosity functions in our simulations on the reionisation of the Universe. This is an important subject that is closely related to our current work, as the source of ionising radiation is considered to be coming from these low mass, star forming galaxies at .
7 Implications for Reionisation & Minimum Halo Mass for Hosting Faint Galaxies
In order to quantify the contribution of the faint, low-mass galaxies to the reionisation of the Universe, we examine the SFR density (SFRD) from galaxies with different stellar masses at . Figure 11 shows the SFRD as a function of redshift broken down into galaxy stellar mass bins of (red), (blue), (black), and the total SFRD (dashed dark green). We also show the observational estimates taken from Bouwens et al. (2010, black squares) and Ouchi et al. (2009, magenta triangles). The fractional contribution from each mass range is listed in Table 4, and its redshift evolution is summarised in Figure 12.
The dotted black lines indicate the SFRD required to maintain reionisation (Madau et al., 1999):
which is dependent on redshift and the ratio of IGM clumping factor () and escape fraction () of ionizing photons from galaxies. Although the exact value for is unknown, a range of values has been estimated from numerical simulations: (Gnedin & Ostriker, 1997), (Iliev et al., 2006), and (Pawlik et al., 2009).
In our recent work, Yajima et al. (2010) estimated the values of as a function of halo mass in our simulations, by performing ray-tracing radiative transfer calculations of ionizing photons from star-forming galaxies. They found that decreases with increasing halo mass. In order to obtain a representative value of , we compute the average , weighted by the halo number density, using Sheth & Tormen (2002) halo mass function, and obtain . Here we use this value to calculate for each assumed value of .
For comparison, the yellow and cyan shaded areas in Figure 11 represent a range of model predictions from Muñoz & Loeb (2010). They developed a semi-analytic model of star-forming galaxies using extended Press-Schechter theory, and examined the effect of minimum halo mass that can host a galaxy on the UV LF. They found that the best-fit UV LF to the current observations can be obtained with a minimum halo mass of . The cyan shade indicates the SFRD derived using this best-fit minimum halo mass and varying . The same is done considering , shown by the yellow shade. By comparison, our simulation is more consistent with the case of , as confirmed by the direct data in our simulation (see the later discussion in Section 8.3 and Figure 16). In our N400L10 run which resolves the lowest halo mass among the three simulations, we find , corresponding to the virial temperature of K. This is expected, because the cooling curve used in our simulation cuts off at around this temperature.
Using the conditional luminosity function method, Trenti et al. (2010) found a similar trend of a steepening faint-end slope from at to at , based on the evolution of dark matter halo mass function. Given these steep faint-end slopes, they argued that galaxies between and contribute % of the total luminosity density at and produce enough ionizing photons to maintain reionisation assuming / (, ). In our simulations, this fraction is even higher ( for the same range as above) due to our steeper . However, if we change the integration limit to corresponding to our resolution limit, we find this fraction to be at . See Table 5 for the fractional contribution by the faint-end galaxies to the total UV luminosity density and stellar mass density at .
Examining the total SFRD found in our simulations (long-dashed line), we find that, depending on our assumptions of the values for , galaxies can maintain reionisation as early as for , and as late as for . In our simulations, low-mass galaxies with (red line) and are the primary contributor to the SFRD at as shown in Table 4, and they are capable of achieving reionisation by . Note that our simulations do not include contributions from Pop III stars.
8.1 Numerical Resolution Issues
In order to assess the resolution limit of the simulation, one may compare the particle mass and smoothing length to the Jeans mass and length. Assuming the SF threshold density ( cm) and a temperature range of K, we find that the Jeans length and mass to be pc and . Compared to these numbers, the proper gravitational softening length of the N400L10 run at is pc, which is shorter than by pc. The gas particle mass of the same run ( with ) is lower than by a factor of .
However, Bate & Burkert (1997) argued for a mass resolution criteria that requires the local Jeans mass to be less than , where is the number of particles in the SPH kernel. For the N400L10 run, was used, resulting in . This number is in-between the range of given above, and it suggests that our N400L10 simulation cannot resolve the gravitational collapse of gas clouds properly at temperatures below 2120 K, similarly to other cosmological SPH simulations of galaxy formation with similar box sizes and particle counts. This is one of the reasons why cosmological simulations need to adopt a sub-particle multiphase ISM model and form star particles out of a cold phase gas that co-exist with a hot phase gas with high effective temperature. Such sub-particle SF models are designed to reproduce the observed Kennicutt (1998) law in a statistical sense (e.g., Springel & Hernquist, 2003; Schaye & Dalla Vecchia, 2008; Choi & Nagamine, 2010), even though the simulations cannot resolve the collapse of the cold gas clouds.
Despite of all the discussion above, one would have to perform a convergence study to be really confident that the results are unaffected by the resolution limit. To show that our results on the faint-end are not resolution dependent, we compare the LFs created from runs with seed particle counts of (long dash black), (dash-dot red) and (short dashed blue) all in a box size in Figure 14. Here we see that the faint-end slope from to remains largely unchanged in the N400L10 run.
8.2 UVB Radiation and H Cooling
Among the results presented in this paper, the most striking one is the very steep faint-end slopes in both galaxy mass and luminosity functions at , and we find a greater number of low-mass galaxies than what the current observations suggest. While it is possible that the current observations are missing these faint galaxies at high-redshift and that the upcoming JWST will detect these faint sources, we should also consider the uncertainties and limitations in the treatment of star formation in our numerical simulations.
In our simulations, stars are allowed to form once the gas density exceeds the threshold density of cm (see Nagamine et al., 2010, for justification of this value), and star particles are generated according to the SF law matched to the local Kennicutt (1998) law based on the cooling functions with metal line contributions (Choi & Nagamine, 2010). If the star formation efficiency decreases with increasing redshift in the real Universe, it is possible that our simulations are overproducing stars at high redshift. This could happen, for example, if star formation is controlled primarily by the amount of molecular gas, which is dependent on the gas metallicity. In the early Universe, ISM and IGM may not be enriched with metals and dust sufficiently yet, which leads to a lower molecular mass densities and lower SFRD (Gnedin & Kravtsov, 2010; Kuhlen et al., 2011).
To investigate the possible impact of H cooling on our result, we experiment with a simple model varying star formation threshold (VTH), in which the value of is increased by a factor of with increasing redshift. The result of this run is shown by the dash-dot blue line in Figure 15. As the SF threshold density is increased and the star formation is limited to regions with higher and higher density, the number of low mass galaxies is decreased significantly at , showing a much flatter faint-end slope than our fiducial runs.
In Figure 15 we also compare our fiducial and VTH model to the KMT09 simulation by Kuhlen et al. (2011), which utilizes an SF model based on H mass (magenta triangles). Although covering only small part of the total LF, we can clearly see that without dust extinction correction, this model overpredicts the observed rest-frame UV LF at , similarly to our simulations. When we apply a dust extinction effect with approximately to the Kuhlen et al.’s data (a shift of 1 mag similarly to our simulations), their result (filled magenta circles) agree with observations very well just like our result. This suggests that the simulations utilizing a SF model based on H mass produce a very similar LF to that presented here, except the sharp turnover at seen for the Kuhlen et al.’s data. We are currently in the process of implementing a similar SF model based on H mass in our simulations, and we will report the results in a subsequent paper (Thompson, et al., in preparation).
8.3 Origin of Steep Faint End Slope
In order to compare our result with some of the semi-analytic models based on the halo occupation model, in Figure 16 we examine the relationship between the rest-frame, monochromatic flux at 1500 Å and the halo mass for all the halos in our simulations at . The dashed blue line is a least-square fit to all the data points shown by the contour, the solid red line is the result of semi-analytic modeling by Trenti et al. (2010), and the yellow shade is another similar model result by Lee et al. (2009). Trenti et al. calibrated their model to the observed LF at , and Lee et al. to the observed LF and correlation functions at , with both assuming single galaxy occupation per halo. Lee et al.’s model implements a mass and luminosity threshold for galaxies, thus the yellow shade does not cover the entire range of our simulation.
We see good agreement between these model results at . Below our relationship begins to deviate from Trenti’s model, being higher by dex. This result implies that, in our simulations the low mass halos with are more efficient at producing stars than in the model of Trenti et al., making our faint-end slope steeper than theirs and Bouwens et al.’s.
To further validate our simulation results, in Figure 17 we compare the dark matter halo mass function (HMF) in our simulations (data points) to the analytic model of Sheth & Tormen (2002, dashed line). The HMF was assembled from our three runs in the same manner as the LF and GSMF. Here we find very good agreement with the model, indicating that we are producing proper number of halos in the mass range of .
Since we are producing correct number of halos, the above results suggest that our steeper than observed faint-end slope can be attributed to the fact that the low mass halos are producing stars more efficiently, rather than too many galaxies in a particular bin. This would have the effect of steepening the faint-end slope, since the deviation occurs primarily in the N400L10 run which controls the faint-end of the LF. Further work must be done to explain the exact baryonic physics of this effect.
Using cosmological SPH simulations, we examined the colors of high-redshift galaxies, quantified the evolution of luminosity and mass functions at , compared our results to the latest WFC3 observational results, and examined the implications for the reionisation of the Universe. Following are the main conclusions of the present work:
We find good agreement with observations in color-color space at , and 8. Most of the simulated galaxies with are consistent with the color selection criteria with a relatively tight distribution. This suggests that the galaxies selected by the current color selection criteria could have a wide range of extinction values with . We also find that the scatter in the observed data on the color-color plane is more dominated by the redshift scatter and photometric errors, rather than the variance in extinction.
The rest-frame UV LF of simulated galaxies agree well with the current WFC3 observations when we assume uniform dust extinctions of for and for . These extinction values are consistent with those obtained by Schaerer & de Barros (2010), who performed SED fitting including nebular emission lines. See Section 5 for more details on the extinction treatment.
We performed least fits to simulation LFs using the three parameter Schechter luminosity and mass functions. The best-fit Schechter parameters are summarised in Tables 2 and 3. Each of the three Schechter parameters, (, , and ) or (, , and ), are evolving with redshift. The faint-end slope is steeper than the current observational estimates with , and becomes steeper with and . The characteristic mass decreases towards higher redshift with . The characteristic magnitude does not evolve very much, but it does get dimmer from to by about 0.3 magnitude. The normalization decreases towards higher redshift as expected in the hierarchical structure formation model, with and .
We decomposed the total SFRD into three separate contributions from different galaxy mass ranges of and , in a similar fashion to the work by Choi & Nagamine (2011). We find that galaxies with are the primary contributor to the total SFRD at as shown in Figs. 11, 12, and Table 4, and therefore to the ionizing photon budget as well. Our simulation suggests that these faint galaxies can reionise the Universe by as long as the clumping factor is . This conclusion is consistent with that of our earlier direct radiative transfer calculations by Yajima et al. (2010).
We find that the contribution from low-mass galaxies with to the total SFRD (and hence the ionizing photon budget) decreases with decreasing redshift, from at to at , which is still dominant. This decreasing trend can also be seen in the contribution from faint-end galaxies () to the total UV luminosity density and stellar mass density (see Table 5 and Figure 13). Whereas the fractional contribution from galaxies with and are increasing gradually during the same epoch with much lower fraction (see Table 4 and Figure 12). The fractional contribution from galaxies with is particularly small with less than 10% at all redshifts of .
In summary, the steep faint-end slopes of high-redshift galaxies at seems to be a generic prediction of CDM model, if one calibrates the SF model to match the locally observed Kennicutt (1998) law or to the lower redshift LF results. In addition to the results of our simulations, the two simple semianalytic models obtain similarly steep faint-end slopes (Muñoz & Loeb, 2010; Trenti et al., 2010), as well as the more detailed semianalytic model of galaxy formation by Lo Faro et al. (2009). The existence of these low-mass galaxies below the current detection limit can be tested soon by the JWST, and we will continue to investigate the possibility of decreasing SF efficiency towards high-redshift by implementing more sophisticated physical processes (e.g., metallicity dependent SF model based on H masses) in our future simulations. Given the significant contribution by the faint-end galaxies to the overall budget of SFRD, ionizing photons, luminosity and stellar mass densities, it is critical to determine the number density and the nature of this population more accurately.
We are grateful to V. Springel for allowing us to use the original version of GADGET-3 code, on which the Choi & Nagamine (2010, 2011) simulations are based. We would also like to thank M. Trenti for communications regarding comparisons made in Figure 16. This work is supported in part by the NSF grant AST-0807491, NASA grant HST-AR-12143-01-A, National Aeronautics and Space Administration under Grant/Cooperative Agreement No. NNX08AE57A issued by the Nevada NASA EPSCoR program, and the President’s Infrastructure Award from UNLV. This research is also supported by the NSF through the TeraGrid resources provided by the Texas Advanced Computing Center. Some numerical simulations and analyses have also been performed on the UNLV Cosmology Cluster.
- Andrae (2010) Andrae R., 2010, ArXiv e-prints
- Bate & Burkert (1997) Bate M. R., Burkert A., 1997, MNRAS, 288, 1060
- Beckwith et al. (2006) Beckwith S. V. W., Stiavelli M., Koekemoer A. M., Caldwell J. A. R., Ferguson H. C., Hook R., Lucas R. A., Bergeron L. E., et al., 2006, AJ, 132, 1729
- Bouwens et al. (2006) Bouwens R. J., Illingworth G. D., Blakeslee J. P., Franx M., 2006, ApJ, 653, 53
- Bouwens et al. (2009) Bouwens R. J., Illingworth G. D., Franx M., Chary R., Meurer G. R., Conselice C. J., Ford H., Giavalisco M., et al., 2009, ApJ, 705, 936
- Bouwens et al. (2008) Bouwens R. J., Illingworth G. D., Franx M., Ford H., 2008, ApJ, 686, 230
- Bouwens et al. (2010) Bouwens R. J., Illingworth G. D., Oesch P. A., Labbe I., Trenti M., van Dokkum P., Franx M., Stiavelli M., et al., 2010, ArXiv e-prints
- Bouwens et al. (2010) Bouwens R. J., Illingworth G. D., Oesch P. A., Stiavelli M., van Dokkum P., Trenti M., Magee D., Labbé I., et al., 2010, ApJL, 709, L133
- Bouwens et al. (2011) Bouwens R. J., Illingworth G. D., Oesch P. A., Trenti M., Labbe I., Franx M., Stiavelli M., Carollo C. M., van Dokkum P., Magee D., 2011, ArXiv e-prints
- Bouwens et al. (2004) Bouwens R. J., Illingworth G. D., Thompson R. I., Blakeslee J. P., Dickinson M. E., Broadhurst T. J., Eisenstein D. J., Fan X., et al., 2004, ApJL, 606, L25
- Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
- Bunker et al. (2004) Bunker A. J., Stanway E. R., Ellis R. S., McMahon R. G., 2004, MNRAS, 355, 374
- Calzetti (1997) Calzetti D., 1997, AJ, 113, 162
- Castellano et al. (2010) Castellano M., Fontana A., Boutsia K., Grazian A., Pentericci L., Bouwens R., Dickinson M., Giavalisco M., et al., 2010, A&A, 511, A20
- Choi & Nagamine (2009) Choi J., Nagamine K., 2009, MNRAS, 395, 1776
- Choi & Nagamine (2010) Choi J.-H., Nagamine K., 2010, MNRAS, 407, 1464
- Choi & Nagamine (2011) Choi J.-H., Nagamine K., 2011, ArXiv e-prints
- Davé et al. (1999) Davé R., Hernquist L., Katz N., Weinberg D. H., 1999, ApJ, 511, 521
- Dickinson et al. (2004) Dickinson M., Stern D., Giavalisco M., Ferguson H. C., Tsvetanov Z., Chornock R., Cristiani S., Dawson S., et al., 2004, ApJL, 600, L99
- Drory et al. (2009) Drory N., Bundy K., Leauthaud A., Scoville N., Capak P., Ilbert O., Kartaltepe J. S., Kneib J. P. o., 2009, ApJ, 707, 1595
- Efstathiou (1992) Efstathiou G., 1992, MNRAS, 256, 43P
- Eisenstein & Hu (1999) Eisenstein D. J., Hu W., 1999, ApJ, 511, 5
- Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
- Finlator et al. (2006) Finlator K., Davé R., Papovich C., Hernquist L., 2006, ApJ, 639, 672
- Gnedin & Kravtsov (2010) Gnedin N. Y., Kravtsov A. V., 2010, ApJ, 714, 287
- Gnedin & Ostriker (1997) Gnedin N. Y., Ostriker J. P., 1997, ApJ, 486, 581
- Gonzalez et al. (2010) Gonzalez V., Labbe I., Bouwens R., Illingworth G., Franx M., Kriek M., 2010, ArXiv e-prints
- Haardt & Madau (1996) Haardt F., Madau P., 1996, ApJ, 461, 20
- Iliev et al. (2006) Iliev I. T., Ciardi B., Alvarez M. A., Maselli A., Ferrara A., Gnedin N. Y., Mellema G., Nakamoto T., et al., 2006, MNRAS, 371, 1057
- Katz et al. (1996) Katz N., Weinberg D. H., Hernquist L., 1996, ApJS, 105, 19
- Kennicutt (1998) Kennicutt Jr. R. C., 1998, ARA&A, 36, 189
- Komatsu et al. (2010) Komatsu E., Smith K. M., Dunkley J., Bennett C. L., Gold B., Hinshaw G., Jarosik N., Larson D., et al., 2010, ArXiv e-prints
- Kuhlen et al. (2011) Kuhlen M., Krumholz M., Madau P., Smith B., Wise J., 2011, ArXiv e-prints
- Lee et al. (2009) Lee K.-S., Giavalisco M., Conroy C., Wechsler R. H., Ferguson H. C., Somerville R. S., Dickinson M. E., Urry C. M., 2009, ApJ, 695, 368
- Lo Faro et al. (2009) Lo Faro B., Monaco P., Vanzella E., Fontanot F., Silva L., Cristiani S., 2009, MNRAS, 399, 827
- Madau (1995) Madau P., 1995, ApJ, 441, 18
- Madau et al. (1999) Madau P., Haardt F., Rees M. J., 1999, ApJ, 514, 648
- Malhotra et al. (2005) Malhotra S., Rhoads J. E., Pirzkal N., Haiman Z., Xu C., Daddi E., Yan H., Bergeron L. E., et al., 2005, ApJ, 626, 666
- Mannucci et al. (2007) Mannucci F., Buttery H., Maiolino R., Marconi A., Pozzetti L., 2007, A&A, 461, 423
- McLure et al. (2010) McLure R. J., Dunlop J. S., Cirasuolo M., Koekemoer A. M., Sabbi E., Stark D. P., Targett T. A., Ellis R. S., 2010, MNRAS, 403, 960
- Muñoz & Loeb (2010) Muñoz J. A., Loeb A., 2010, ArXiv e-prints
- Nagamine et al. (2005a) Nagamine K., Cen R., Hernquist L., Ostriker J. P., Springel V., 2005a, ApJ, 627, 608
- Nagamine et al. (2005b) Nagamine K., Cen R., Hernquist L., Ostriker J. P., Springel V., 2005b, ApJ, 618, 23
- Nagamine et al. (2010) Nagamine K., Choi J.-H., Yajima H., 2010, ApJL, 725, L219
- Nagamine et al. (2004) Nagamine K., Springel V., Hernquist L., Machacek M., 2004, MNRAS, 350, 385
- Night et al. (2006) Night C., Nagamine K., Springel V., Hernquist L., 2006, MNRAS, 366, 705
- Oesch et al. (2010) Oesch P. A., Bouwens R. J., Illingworth G. D., Carollo C. M., Franx M., Labbé I., Magee D., Stiavelli M., et al., 2010, ApJL, 709, L16
- Oesch et al. (2009) Oesch P. A., Carollo C. M., Stiavelli M., Trenti M., Bergeron L. E., Koekemoer A. M., Lucas R. A., Pavlovsky C. M., et al., 2009, ApJ, 690, 1350
- Ouchi et al. (2009) Ouchi M., Mobasher B., Shimasaku K., Ferguson H. C., Fall S. M., Ono Y., Kashikawa N., Morokuma T., Nakajima K., Okamura S., Dickinson M., Giavalisco M., Ohta K., 2009, ApJ, 706, 1136
- Pawlik et al. (2009) Pawlik A. H., Schaye J., van Scherpenzeel E., 2009, MNRAS, 394, 1812
- Reddy & Steidel (2009) Reddy N. A., Steidel C. C., 2009, ApJ, 692, 778
- Richard et al. (2006) Richard J., Pelló R., Schaerer D., Le Borgne J., Kneib J., 2006, A&A, 456, 861
- Robertson et al. (2007) Robertson B., Li Y., Cox T. J., Hernquist L., Hopkins P. F., 2007, ApJ, 667, 60
- Schaerer & de Barros (2010) Schaerer D., de Barros S., 2010, A&A, 515, A73+
- Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
- Schechter (1976) Schechter P., 1976, ApJ, 203, 297
- Sheth & Tormen (2002) Sheth R. K., Tormen G., 2002, MNRAS, 329, 61
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
- Trenti et al. (2010) Trenti M., Stiavelli M., Bouwens R. J., Oesch P., Shull J. M., Illingworth G. D., Bradley L. D., Carollo C. M., 2010, ApJL, 714, L202
- Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
- Wilkins et al. (2011) Wilkins S. M., Bunker A. J., Lorenzoni S., Caruana J., 2011, MNRAS, 411, 23
- Yajima et al. (2010) Yajima H., Choi J., Nagamine K., 2010, ArXiv e-prints
- Yan et al. (2009) Yan H., Windhorst R., Hathi N., Cohen S., Ryan R., O’Connell R., McCarthy P., 2009, ArXiv e-prints
- Yan & Windhorst (2004) Yan H., Windhorst R. A., 2004, ApJL, 612, L93