MOIRCS DEEP SURVEY. IV. EVOLUTION OF GALAXY STELLAR MASS FUNCTION BACK TO z ∼ 3

In order to estimate the photometric redshift and stellar mass of the sample galaxies, we performed SED fitting of the multiband photometry described above (UBVizJHK, 3.6 μm, 4.5 μm, and 5.8 μm) with population synthesis models. We adopted the standard minimum χ² method for the fitting procedure. The resulting best-fit redshift (i.e., photometric redshift) is used for objects without spectroscopic identifications and the best-fit stellar mass-to-luminosity (M/L) ratio is used to calculate the stellar mass. We adopted spectroscopic redshifts (if available) from the literature (Cohen et al. 2000; Cohen 2001; Dawson et al. 2001; Wirth et al. 2004; Cowie et al. 2004; Treu et al. 2005; Chapman et al. 2005; Reddy et al. 2006a; Barger et al. 2008), including those from our NIR spectroscopic observation of ∼ 20 star-forming BzK galaxies (Daddi et al. 2004) at z ∼ 2 with Subaru/MOIRCS (T. Yoshikawa et al. 2009, in preparation). For these objects, the SED fitting was performed fixing the redshift to each spectroscopic value.

In this study, three models, i.e., GALAXEV (Bruzual & Charlot 2003), PEGASE version 2 (Fioc & Rocca-Volmerange 1997), and Maraston (2005) model, were used as SED templates. Comparisons among the results with different SED template sets allow us to check the systematic effects of these models on the estimate of redshift and stellar mass. We also used a public photometric redshift code, EAZY (Brammer et al. 2008) for an independent check of the photometric redshift. In all SED models, Salpeter IMF (Salpeter 1955) with lower and upper mass limits of 0.1 and 100 M_☉ is adopted for easy comparison among the results with the models and those in other studies. The details of the model templates from each population synthesis model are given in the following.

• 1.  
  GALAXEV. We used the similar template sets described in Kajisawa & Yamada (2005). Exponentially decaying SFHs with the decaying timescale τ ranging between 0.1 and 20 Gyr are assumed. Calzetti extinction law (Calzetti et al. 2000) in the range of E(B − V) = 0.0–1.0 is adopted. Metallicity is changed from 1/50 to 1 solar metallicity.
• 2.  
  PEGASE2. We used the same template sets as those described in Grazian et al. (2006). Eight exponentially decaying SFHs with various gas infalling and SFR timescales are assumed (τ_gas = 0.1–5 Gyr and τ_SFR = 0.1–20 Gyr; see Table 4 in Grazian et al. 2006 for details). The model includes the metallicity evolution and dust extinction in a self-consistent way. Following Grazian et al. (2006), we also added passively evolving templates with a constant SFR and several truncation ages, and dusty star-forming templates with a constant SFR and the Calzetti extinction with the range of E(B − V) = 0.5–1.1.
• 3.  
  Maraston (2005). The same SFHs and the same range of the Calzetti extinction as in the case of GALAXEV are assumed (τ = 0.1–20 Gyr and E(B − V) = 0.0–1.0). Metallicity is changed from 0.5 to 2.0 solar metallicity. This model increases emphasis on contributions of thermally pulsating AGB stars to the SEDs. In all these 1–3 models, the model age is changed from 50 Myr to the age of the universe at the observed redshifts. We also added to the model SEDs the Lyman series absorption produced by the intergalactic medium following Madau (1995).
• 4.  
  EAZY. We used the default parameter setting and the default template set of eazy_v1.0. The default template set consists of five principal component templates constructed from a large number of PEGASE models that cover a distribution of the SFH of galaxies in the semi-analytic model by De Lucia & Blaizot (2007). In addition to these principal components, a dusty star-forming SED is added as a complementary template. EAZY fits the multiband photometry of the observed galaxies with a linear combination of these templates (see Brammer et al. 2008 for details). Since EAZY fits with the templates based on the principal components, we cannot derive the precise stellar M/L ratio from the fitting procedure. In order to estimate the stellar mass, therefore, we performed SED fitting with the 1–3 models, fixing the redshift to the output from EAZY. The results with the different population synthesis models at the same redshifts also enable us to check the systematic effect of these models on the estimate of stellar M/L ratio independently.

In Figure 1, we compare the photometric redshifts with a sample of 2102 spectroscopic redshifts from the literature described above. The photometric redshift accuracy and the fraction of the catastrophic failure with δz/(1+z_spec)>0.5 for each SED model are (δz/(1+z_spec) = −0.009 ± 0.077, 4.0% outliers) for GALAXEV, (−0.002 ± 0.103, 4.0%) for PEGASE2, (−0.007 ± 0.082, 3.9%) for the Maraston model, and (+0.003 ± 0.098, 4.7%) for EAZY, respectively. Although the photometric redshift accuracy is relatively good in all cases, there are some systematic differences among the models especially at 1.5 < z < 3.0 in Figure 1. We check the effects of these differences of the photometric redshifts on the SMF in Section 3.4.

Zoom In Zoom Out Reset image size

Table 1 lists the number of objects in each redshift bin for different SED models. These redshift bins, of which width is sufficiently larger than the typical photometric redshift errors, were defined so as to include a reasonable number of galaxies in the bin for calculating the SMF. The comoving volumes are 8.5 × 10⁴ Mpc³ (2.5 × 10⁴ Mpc³ for the deep field) at 0.5 < z < 1.0, 1.4 × 10⁵ Mpc³ (3.7 × 10⁴ Mpc³) at 1.0 < z < 1.5, 3.4 × 10⁵ Mpc³ (9.2 × 10⁴ Mpc³) at 1.5 < z < 2.5, and 3.4 × 10⁵ Mpc³ (9.3 × 10⁴ Mpc³) at 2.5 < z < 3.5, respectively. The parenthesis in Table 1 represents the number of objects with spectroscopic redshift. The fraction of spectroscopically identified sources is relatively high, thanks to the extensive spectroscopic surveys in this field.

As discussed in previous studies (e.g., Papovich et al. 2001; Kajisawa & Yamada 2005; Shapley et al. 2005), the uncertainty of the stellar mass is smaller than other parameters such as stellar age, star formation timescale, metallicity, and dust extinction. While the broadband SEDs are degenerated with respect to these parameters of the stellar populations (stellar age, star formation timescale, metallicity, and dust extinction), the stellar M/L ratio (and therefore stellar mass) is much less affected by the degeneracy because the effects of these parameters on the stellar M/L ratio tend to be canceled out with each other. Figure 2 shows the uncertainty of the stellar mass estimated from a χ² map in SED fitting with the GALAXEV model for each object. For sources without spectroscopic identifications, we also varied redshift as a free parameter in the calculation of the χ² map to take into account the photometric redshift error. The uncertainty increases with decreasing stellar mass and increasing redshift. The stellar mass errors become ∼ 0.3–0.4 dex at the limiting stellar mass (the vertical dashed lines in the figure) described in the following subsection.

Zoom In Zoom Out Reset image size

The K-band magnitude-limited sample does not have a sharp limit in stellar mass even at a fixed redshift, because the stellar M/L ratio at the observed K band varies with different stellar populations. We used the rest-frame U − V color distribution as a function of stellar mass in each redshift bin to estimate the limiting stellar mass above which most of the galaxies are expected to be brighter than the magnitude limits and detected in the K-band image. Since rest-frame U − V is correlated with stellar M/L ratio (e.g., Rudnick et al. 2003; Marchesini et al. 2007), one can predict the mass dependence of stellar M/L and then estimate the effect of the magnitude limit on the stellar mass distribution.

Figure 3 shows the rest-frame U −V distribution in each redshift bin for the wide sample (top panels) and the deep sample (middle panels). A dashed line in each panel represents the K-band magnitude limit (K = 23 for the wide sample and K = 24 for the deep one). All objects with stellar mass larger than this line (on the right side of the line in the figure) at each U −V value are brighter than the magnitude limit. In order to calculate the line, we used the stellar M/L ratio and the rest-frame color of the GALAXEV model with various SFHs, dust extinction, and metallicity. The maximum mass was selected from the possible range of the stellar mass for the models with K = 23 (or 24 for the deep sample) and each U −V color in order to depict the dashed line in the figure. Dashed-dot lines in Figure 3 represent the 90 percentile of U −V color at each stellar mass. We adopted the point where the lines of the magnitude limit and 90 percentile of U −V color cross each other as the limiting stellar mass. Above this limiting mass, more than 90% of objects are expected to be brighter than the K-band magnitude limit.

Zoom In Zoom Out Reset image size

The bottom panels in Figure 3 show the comparison of the 90 percentiles of U −V color for the wide and deep samples. The 90 percentiles for both the samples agree well with each other even near the limiting mass of the wide sample. This suggests that the incompleteness near the magnitude limit of the wide sample does not strongly affect the U −V color distribution, although some red galaxies might be missed on the left-hand side of the long-dashed line.

Figure 4 shows the calculated limiting stellar mass as a function of redshift for the wide and deep samples. We also estimated these mass limits for the cases with PEGASE2 and Maraston models as well as the case with GALAXEV. As seen in Figure 3 and other previous studies (e.g., Kajisawa & Yamada 2005; Kajisawa & Yamada 2006; Labbé et al. 2005; Taylor et al. 2009), less massive galaxies tend to have bluer rest-frame color even at high redshift. Such a mass-dependent color distribution can be seen well above the K-band magnitude limit up to at least z ∼ 2.5 in Figure 3. Since the bluer color of low-mass galaxies indicates a lower M/L ratio, we can detect galaxies down to the relatively lower mass limit with high completeness compared with, for example, the mass limit based on the M/L ratio of the passively evolving models, which is used in other previous studies (e.g., Dickinson et al. 2003; Fontana et al. 2004). On the other hand, our "wide" field data are relatively shallow for galaxies at z > 3 (top right panel in Figure 3) and the completeness is relatively low even at high mass, where galaxies tend to have red rest-frame colors (high M/L ratios). We use objects with the stellar mass larger than these mass limits at each redshift to estimate and discuss the SMF in the following.

Zoom In Zoom Out Reset image size

The SMF of galaxies was derived with the non-parametric 1/V_max formalism and the parametric STY method (Sandage et al. 1979). Both methods are commonly used to estimate the luminosity function and SMF of galaxies.

In the V_max method, V_max was calculated with the best-fit model SED template for each galaxy. For each best-fit SED, we estimated the K-band apparent magnitude as a function of redshift, taking into account both the luminosity distance and K correction. Then we determined the maximum redshift, z_max above which the object becomes fainter than the K-band magnitude limit (K = 23 for the wide sample or K = 24 for the deep sample). V_max is a comoving volume integrated from the lower limit of each redshift bin to z_max or the upper limit of the bin (the smaller of these two). Then 1/V_max estimates were used to calculate the number density of galaxies in each mass bin.

In the STY method, assuming the Schechter function form (Schechter 1976) for the SMF, we estimated best-fit values of the Schechter parameters (α, M*, ϕ*). The limiting stellar mass M_lim(z) (described in the previous subsection) for the redshift of each object was used to calculate the probability that the object has the observed stellar mass M as

$$\begin{equation} p=\frac{\phi (M)}{\int ^{\infty }_{M_{\rm lim}(z)}\phi (M) dM}. \end{equation} \tag{ 1 }$$

Here ϕ(M) is the SMF represented by the Schechter function. We searched for the values of the Schechter parameters (α, M*, ϕ*) which maximize the likelihood L = ∏p, the products of the probability densities for the objects with the stellar mass larger than M_lim(z) in each redshift bin. Both the wide and deep samples were used simultaneously in the maximum likelihood technique.

Figure 5 shows the SMF of galaxies in the different redshift bins for the different SED models. The results of the V_max method and the best-fit Schechter function estimated with the STY method are plotted in each panel. Error bars are based on the Poisson statistics. Dashed lines show the local SMF derived from the 2dF and Two Micron All Sky Surveys (Cole et al. 2001) with the small correction for the "maximum age" method as described in Fontana et al. (2004). In Figure 6, we plot the combined wide and deep complete data (same as the solid symbols in Figure 5) for the different SED models with different symbols in the same panel.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 5 and 6 show that the SMFs obtained from different samples and SED models are in good agreement, although there are some systematic differences among the SMFs. The number densities for the deep sample are systematically larger (by ∼ 0.1 dex) than those for the wide sample at 0.5 < z < 1.0 for all SED models. Our deep field is centered at the HDF-N field where the extensive spectroscopic surveys revealed the large-scale structures at z = 0.85 and z = 1.02 (Cohen et al. 2000; Wirth et al. 2004). These filaments or clumps around the HDF-N could cause the slightly larger number density in our deep field, and such differences can be considered as the possible field-to-field variance. Slighter excess in the deep field can also be seen at 2.5 < z < 3.5 in most SED-model cases, and this may also be due to large-scale structures. On the other hand, the systematic differences among the different SED models in Figure 6 seem to be larger. The number density of galaxies for the Maraston model is systematically smaller by ∼ 0.15–0.2 dex (up to ∼ 0.5 dex at > 10¹¹ M_☉) than those for the GALAXEV and PEGASE2 models especially at z > 1.5. Since the differences are also seen between the EAZY+Maraston and EAZY+GALAXEV/PEGASE2, where the same redshifts of EAZY are used, the differences of the estimated stellar M/L ratio among the SED models cause the systematic differences of the number density. For example, Maraston et al. (2006) performed the broadband SED fitting of relatively young (∼ 0.2–2 Gyr stellar age) galaxies at 1.4 < z < 2.7 with the Maraston and GALAXEV models and reported that the Maraston model gives systematically younger age and lower stellar mass (∼ 60%) than the GALAXEV model. Such difference of the estimated stellar mass seems to explain the differences seen in Figure 6. If the fraction of young galaxies becomes larger at high redshift, the larger differences in the SMFs at z > 1.5 could also be explained because the contribution of TP-AGB stars is expected to be significant in the relatively young ages (∼ 0.2–2 Gyr, Maraston 2005; Maraston et al. 2006).

We can see general evolutionary features in Figures 5 and 6 in spite of the field variance and the systematic differences among the SED models mentioned above. We note, first, that the overall number density decreases gradually with redshift in all cases. While the number density of galaxies at 0.5 < z < 1.0 is similar with that in the local universe, the number density at 2.5 < z < 3.5 is about an order of magnitude smaller than the local value. Figure 7 and Table 2 show the best-fit Schechter parameters (α, M*, ϕ*) estimated with the STY method. The redshift evolution of the overall number density can be seen as the decrease of the normalization of the SMF ϕ*. ϕ* decreases down to ∼ 50% of that in the local universe at z ∼ 0.75, ∼ 16% at z ∼ 1.25, ∼ 9% at z ∼ 2, and ∼ 3% at z ∼ 3. Similar evolution of the SMF is also seen in previous studies of general fields (e.g., Fontana et al. 2006; Pérez-González et al. 2008; Marchesini et al. 2009).

Zoom In Zoom Out Reset image size

Second, we found the mass-dependent evolution of the SMF. The evolution of the number density of low-mass galaxies with M_star ∼ 10⁹–10¹⁰ M_☉ is smaller than that of the massive galaxies with stellar mass ∼10¹¹ M_☉. While the number density of galaxies with ∼10¹¹ M_☉ at 2.5 < z < 3.5 is smaller by a factor of ∼ 20 than the local value, that of the galaxies with ∼5 × 10⁹ M_☉ is smaller by only a factor of ∼ 6 at the same redshift. The shape of the SMF at 0.5 < z < 1.0 is similar with that in the local universe, and it becomes steeper with redshift at z > 1. This can be seen in Figure 7 and Table 2 as a steepening of the low-mass slope α with redshift. α decreases with redshift gradually from α = −1.29 ± 0.03(±0.04) at 0.5 < z < 1.0 to α = −1.48 ± 0.06(±0.07) at 1.5 < z < 2.5 and α = −1.62 ± 0.14(±0.06) at 2.5 < z < 3.5. The quoted errors are statistical errors estimated from the maximum likelihood method. The values in parenthesis show the uncertainty due to the different SED models.

The uncertainty of the Schechter parameters becomes larger with redshift, especially at 2.5 < z < 3.5, because of the larger limiting mass at higher redshift and of the small number of galaxies even at relatively high mass due to the evolution of the overall number density mentioned above. Nonetheless, since the uncertainty at z < 2.5 is rather small by virtue of our deep and wide NIR data, the evolution of the low-mass slope of the SMF between 0 ≲ z ≲ 3 is found to be significant. Figure 8 shows the best-fit Schechter parameters and their uncertainty for the different SED models in the M*–α plane, which represents the evolution of the shape of the SMF. In all cases, the evolution of α is significant, although there is degeneracy between M* and α especially at 2.5 < z < 3.5, where we can reach only to the relatively high stellar mass.

Zoom In Zoom Out Reset image size

On the other hand, the characteristic mass M* shows no significant evolution except for the results with the Maraston model. The M* values at 0.5 < z < 3.5 are similar with or slightly larger than those in the local universe. No significant evolution for M* is also seen in previous studies (Fontana et al. 2006; Pozzetti et al. 2007; Marchesini et al. 2009). For the Maraston model, M* becomes smaller by a factor of ∼ 2–2.5 at z > 1.5. As in the above discussion of the overall number density, this can be explained by the systematically lower stellar M/L ratio of the Maraston model because the same result is also seen when the same photometric redshifts (EAZY) are used. While passively evolving galaxies dominate the massive end (≳10¹¹ M_☉) of the SMF at z ≲ 1 (e.g., Juneau et al. 2005; Borch et al. 2006; Vergani et al. 2008; Ilbert et al. 2009), many massive star-forming (i.e., relatively young) galaxies have been found at z ∼ 2 (e.g., Daddi et al. 2007; Papovich et al. 2006; Borys et al. 2005; Shapley et al. 2004). It is possible that TP-AGB stars contribute to the SED of massive galaxies significantly only at z > 1.5. Therefore, the Maraston model would give systematically lower stellar mass for these massive galaxies.

Here we investigate possible biases for the results in the previous subsection, in particular, systematic effects which could cause the steepening of the low-mass slope at high redshift.

The larger limiting mass at higher redshift causes the degeneracy between M* and α as mentioned in the previous section. Since low-mass galaxies near the limiting mass tend to be faint in each band, the photometric errors are larger, which results in the large uncertainty in their photometric redshift. The large errors of the photometric redshifts of faint objects might lead to the systematic increase of the low-mass galaxies at high redshift because the redshift distribution of the K-selected sample has a peak around z ∼ 1 and a tail to higher redshift. In order to evaluate the effect on the low-mass slope, we performed a Monte Carlo simulation, assuming no evolution of the shape of the SMF (i.e., constant M* and α). At first, we constructed mock catalogs with the observed M* and α at 0.5 < z < 1.0. For the normalization of the SMF, we assumed ϕ* evolving as ϕ*(z) ∝ (1 + z)⁻² so that the redshift distribution of the mock sample is consistent with our observation. Even if we assume ϕ*(z) ∝ (1 + z)⁻¹ or ϕ*(z) ∝ (1 + z)⁻³, the results shown in the following do not change significantly. Stellar mass and redshift of mock objects were randomly selected from the ranges of 10⁸ M_☉< M_star < 10¹² M_☉ and 0 < z < 6, using the probability distribution estimated from the assumed SMF and the corresponding comoving volume of our survey at each redshift. For each mock object, we randomly extract an object from the observed sample with similar mass and redshift (allowing duplicate) and adopted its observed multiband photometry. The observed multiband photometry was extracted from the catalog that contains all sources detected on the K-band image (including K > 23 or K > 24 objects) in order to take into account the scattering of objects fainter than the magnitude limits. Then we added random offsets to the multiband photometry according to the measured photometric errors, and adopted the mock object if the resulting K-band magnitude of the object was brighter than the magnitude limits (K < 23 for the wide sample and K < 24 for the deep sample). We repeated this procedure and made the mock catalogs with the same sample sizes as the observed wide and deep samples. The same SED-fitting procedure as for the observed one was performed in order to estimate the photometric redshift and stellar mass of the mock objects. For objects with spectroscopic identification, redshifts are fixed to the spectroscopic values. We performed 200 simulations and calculated the best-fit Schechter parameters in each simulation with the STY method.

Figure 9 shows the results of the simulation in the case with GALAXEV (the results for the other models are similar). A relatively large scatter of α is seen in the highest redshift bin, which probably reflects the large degeneracy between M* and α due to the large limiting mass. Furthermore, the simulated α distributes around a systematically steeper (by ∼ 0.1–0.15) value than the assumed one at 2.5 < z < 3.5, while the simulated values tend to be slightly flatter (by ∼ 0.1) in lower redshift bins. However, since the observed evolution of α is much stronger than the systematic effects in Figure 9, the observed steepening of the low-mass slope at high redshift is significant, especially at 1 < z < 2.5.

Zoom In Zoom Out Reset image size

In the simulation, the random offsets added to the multiband photometry and the recalculation of the photometric redshift include the effect of the catastrophic failure. We extracted the mock objects whose photometric redshift was changed catastrophically (δz/(1 + z)>0.5) by the random offsets and checked the effect of these objects on the resulting SMF. Figure 10 shows the fractional increase and decrease of the number of galaxies due to the catastrophic failure of the photometric redshift as a function of stellar mass in each redshift bin. At 0.5 < z < 1.5, there is ∼10%–20% decrease near the limiting stellar mass, while there is only negligible effect of the catastrophic failure at M_star ≳ 10¹⁰ M_☉. Most contamination from z < 0.5 or z > 1.5 occurs only at the stellar mass lower than limiting mass and it is not plotted in the figure. About 10%–20% decrease near the limiting stellar mass might cause a slightly flatter low-mass slope seen in Figure 9, but the effect is relatively small (∼0.1 dex decrease in the number density). At z > 1.5, the effect of the catastrophic failure is similar or even smaller than that at low redshift. Furthermore, the contamination from lower redshift and the dropout from the redshift bin tend to be canceled out with each other, which results in the negligible effect on the SMF.

Zoom In Zoom Out Reset image size

In Figure 11, we also show the results of the same simulation on the α–ϕ* plane to investigate the effect on the normalization. The systematic effect on the normalization of the SMF due to the photometric redshift errors seems to be relatively small, although the degeneracy between the parameters makes direct comparison difficult.

Zoom In Zoom Out Reset image size

Next, we discuss a possibility of the over-deblending of faint objects with relatively low S/N ratio near the detection limit. Although the K band, where we performed the source detection, corresponds to the rest-frame B band even at z = 3.5, the morphological K correction could enhance the over-deblending at high redshift because galaxies tend to show the patchy appearance in shorter wavelengths due to the dominance of young stars and the dust extinction (e.g., Kuchinski et al. 2001; Rawat et al. 2009). Figure 12 shows the fraction of the objects with the deblending flag by SExtractor as a function of stellar mass in each redshift bin for the wide and deep samples. We cannot see the mass nor redshift dependence of the fraction of the deblended objects. We conclude that the effect of deblending does not cause the evolution of the low-mass slope.

Zoom In Zoom Out Reset image size

Finally, we estimated the stellar mass of galaxies in the highest redshift bin with the GALAXEV templates with two-component (old and young) SFHs. If the old stellar population is hidden by recent star formation, the stellar mass could be underestimated especially for relatively blue low-mass galaxies at high redshift, where the S/N ratio tends to be low (e.g., Papovich et al. 2001; Drory et al. 2005). We used the exponentially decaying star formation models (young component) with an old population component. For the young population, free parameters are stellar age, star formation timescale τ, color excess, metallicity (same as for one-component SFH). For the old component, we limited the star formation timescale and stellar age to shorter and older values than the young population, respectively, and assumed no dust extinction. Figure 13 shows the comparison of the stellar masses estimated with one- and two-component models for galaxies at 2.5 < z < 3.5. No significant systematic difference of the stellar mass can be seen. The scatter is consistent with the uncertainty of the stellar mass of these galaxies shown in Figure 2, although the stellar mass with the two-component model is slightly larger for a small fraction of galaxies at low-mass end. The long wavelength data with Spitzer/IRAC, which sample the rest-frame NIR region even at high redshift, could make the systematic uncertainty relatively small (Fontana et al. 2006; Elsner et al. 2008).

Zoom In Zoom Out Reset image size

In this subsection, we compare our results in the MODS field with those in other general fields.

Figure 14 shows a comparison of the integrated stellar mass density as a function of redshift. We calculated the stellar mass density by integrating the best-fit Schechter function over 10⁸–10¹³M_☉ in each redshift bin. Different symbols represent the results with the different SED models. Compilations of Wilkins et al. (2008) and Marchesini et al. (2009) are also plotted. All data points are the results with the Salpeter IMF. The stellar mass density in the MODS field decreases with redshift as seen in other previous studies. When compared with those in the local universe, the stellar mass density is ∼53%–72% of the local value at 0.5 < z < 1.0, ∼22%–34% at 1.0 < z < 1.5, ∼8%–18% at 1.5 < z < 2.5, and ∼4%–9% at 2.5 < z < 3.5. For the GALAXEV model, which is widely used in other previous studies (e.g., Drory et al. 2005; Fontana et al. 2006; Pozzetti et al. 2007; Elsner et al. 2008; Marchesini et al. 2009), our results are consistent with those of previous studies at each redshift (Figure 14). A slightly large value at 0.5 < z < 1.0 is probably due to the large-scale structures around the HDF-N mentioned above. As discussed in Section 3.4, the Maraston model gives systematically lower stellar mass density, which is about three-quarters of those with the other two models (GALAXEV and PEGASE2) at 0.5 < z < 1.0, about two-thirds at 1.0 < z < 1.5 and about a half at z > 1.5.

Zoom In Zoom Out Reset image size

In Figure 15, we compare the SMF of galaxies in the MODS field with those in other studies. For simplicity, we plotted the result with the GALAXEV model. At the massive end (≳10¹¹ M_☉), our result is consistent with those in other general fields, although the uncertainty of our data and the variance among the different surveys are relatively large because of the small number of observed massive galaxies and probably their strong clustering (e.g., van Dokkum et al. 2006; Conselice et al. 2007; Stringer et al. 2009). The agreement of the SMF at the massive end is consistent with the comparison of the integrated stellar mass density because galaxies around M* (∼10¹¹ M_☉) dominate the stellar mass density (e.g., Gwyn & Hartwick 2005, see also next subsection). No significant evolution of M* seen in the MODS field at 0.5 < z < 3.5 is also consistent with results of previous studies (Fontana et al. 2006; Bertone et al. 2007; Pozzetti et al. 2007; Elsner et al. 2008; Pérez-González et al. 2008; Marchesini et al. 2009).

Zoom In Zoom Out Reset image size

On the other hand, at the lower mass of the SMF, there is some variance in different fields including the MODS field in spite of the relatively small statistical errors. In addition to the differences in the normalization, which are similar offsets with those at the massive end, there is differences in the low-mass slope especially at high redshift. One possible reason for this is the different depths of the surveys and therefore the different limiting stellar masses. The SMF in the MODS field shows a slight upturn around ∼10¹⁰ M_☉ and a steeper slope at <10¹⁰ M_☉, although this is marginal in the highest redshift bin. This indication can be seen in the results with all SED models (Figure 5) and does not depend on the systematic error of the photometric redshift and stellar M/L ratio. Similar upturns in the SMF are also seen in field galaxies in the local universe (Baldry et al. 2008) and at intermediate and high redshift (Fontana et al. 2006; Pozzetti et al. 2007; Elsner et al. 2008). Therefore, the relatively shallow survey with the limiting stellar mass larger than ∼10¹⁰ M_☉ could miss the upturn and result in the systematically flatter low-mass slope even if the completeness to the limit is sufficiently high or the incompleteness correction is properly done. For example, the SMFs from the different surveys show consistent low-mass slopes in the 0.5 < z < 1.0 bin in Figure 15, where all the surveys reach below 10¹⁰ M_☉. At 1.0 < z < 1.5, several surveys reach only to ∼10¹⁰ M_☉. While the slope at 10¹⁰–10¹¹ M_☉ is similar with those in other surveys, the SMF in the MODS field shows the steeper slope at 10⁹–10¹⁰ M_☉ than that at larger stellar mass. This could be the case in the higher redshift bins.

At z > 1.5, on the other hand, Elsner et al. (2008) and Marchesini et al. (2009) reach below 10¹⁰ M_☉ and show the flatter low-mass slopes than those in the present study. In Elsner et al. (2008), the low-mass slope is fixed to the weighted mean value α ∼ −1.36, which is constrained mainly by galaxies at relatively low redshift (z ≲ 1, their Figure 7). Marchesini et al. (2009) discussed the effects of various random and systematic uncertainties on the estimate of the Schechter parameters of the SMF and suggested that the uncertainty of the low-mass slope in their analysis is relatively large especially at high redshift due to the statistical error and the systematic errors of the different SED models. The difference in the low-mass slope between Marchesini et al. (2009) and the present study could be attributed to field-to-field variance. In the local universe, several studies suggest that the shape of the SMF depends on the environment of galaxies (Balogh et al. 2001; Baldry et al. 2006; Baldry et al. 2008), and this could be the case at higher redshift (e.g., Scodeggio et al. 2009). In Marchesini et al. (2009), the low-mass slope at high redshift is constrained mainly by the FIRES data, whose survey area is about a fifth of our study. If the shape of the SMF depends on the environment even at high redshift, the low-mass slope estimated from the small area surveys could be affected by the field-to-field variance. A larger area survey with a similar or fainter limiting magnitude in the NIR wavelength is needed to conclude whether there exists the field-to-field variance of the shape of the SMF at z∼ 2–3 or not.

The decrease of the overall number density and the steepening of the low-mass slope with redshift seen in the evolution of the SMF of galaxies in the MODS field suggest the possibility of the mass-dependent number density evolution. Figure 16 shows the evolution of the number density of galaxies in different stellar mass ranges for the different SED models. Relatively low-mass galaxies with the stellar mass of 10^9.5–10^10.5 M_☉ tend to show a weaker evolution of the number density especially at z ≳ 1.5 than M* (∼10¹¹ M_☉) galaxies, although the trend is relatively marginal in the results with the GALAXEV/PEGASE2 models. The relative number density of low-mass galaxies had been larger at high redshift than that in the local universe. In other words, the number density of galaxies around M* might have more rapidly increased at z ≳ 1.5 than that of low-mass galaxies. On the other hand, several studies based on wide area but shallow surveys suggest that for M_star > 10¹¹ M_☉, more massive (e.g., ≳ 10^11.5–10¹² M_☉) galaxies show weaker evolution in the number density at z ≳ 1.5 (Berta et al. 2007; Pérez-González et al. 2008). Although the survey volume of the present study is not large enough to constrain strongly the number density of very massive galaxies with 10^11.5–10¹² M_☉, galaxies around 10¹¹ M_☉ might have been formed effectively at 1 ≲ z ≲ 3 compared with higher and lower mass galaxies (e.g., Fontana et al. 2004; Franceschini et al. 2006).

Zoom In Zoom Out Reset image size

Figure 17 shows the contribution of galaxies in different stellar mass ranges to the cosmic stellar mass density. While the stellar mass density is dominated by massive galaxies with M_star > 10^10.5 M_☉ in the local universe, the contribution of lower mass galaxies increases with redshift and is significant at z ∼ 3. Reddy & Steidel (2009) reported a relatively steep faint-end slope (α ∼ −1.7) of the luminosity function for Lyman Break Galaxies (LBGs) at 1.9 < z < 3.4. They estimated the stellar mass of faint LBGs with the log-linear relation between SFR and stellar mass derived in Sawicki et al. (2007), and suggested that those faint LBGs, which typically have M_star ≲ 10¹⁰ M_☉ (Figure 12 in Reddy & Steidel 2009), show a significant contribution to the total stellar mass density. Since most low-mass galaxies at z∼ 2–3 have rather blue rest-frame colors (typically rest U − V ≲ 0, Figure 3; Kajisawa & Yamada 2005; Kajisawa & Yamada 2006) and satisfy the LBG criteria (Labbé et al. 2005), the steep low-mass slope of the SMF at high redshift in our study seems to be consistent with the result in Reddy & Steidel (2009). Note that the integrated stellar mass density in the MODS field is consistent with those of other previous studies as mentioned in the previous subsection, but about a half of the estimate in Reddy & Steidel (2009). Then, although the contribution of low-mass galaxies to the total stellar mass density in the MODS field is larger than those estimated in previous studies, the stellar mass density is still smaller than that expected from the evolution of the cosmic SFR density (e.g., Wilkins et al. 2008; Reddy & Steidel 2009).

Zoom In Zoom Out Reset image size

What is the origin of the evolution of the low-mass slope? In hierarchical structure formation scenarios, low-mass objects continue to grow by merging over cosmic time and result in more massive objects at later epochs. Therefore, one can expect the steeper low-mass slope at the higher redshift (e.g., Khochfar et al. 2007; Ryan et al. 2007). Figure 18 compares the observed SMF in the MODS field with the predictions of the publicly available semi-analytic models by Bower et al. (2006) and Bertone et al. (2007). These two models are based on the Millennium Simulation of the growth of dark matter structure in a Λcold dark matter cosmology (Springel et al. 2005) but are made by the different procedures in many aspects such as, e.g., the construction of halo merger trees, the implementation of the star formation, the feedback by supernova explosions and active galactic nuclei activities. We convolved these model predictions with the Gaussian representing the uncertainty of the stellar mass estimate shown in Figure 2. The median value in Figure 2 is chosen as the width of the Gaussian depending on galaxy's stellar mass and redshift. In the models, the low-mass slope similarly steepens with redshift in Figure 18. However, it should be noted that there are discrepancies between the observation and the models in the detailed shape of the SMF and that the evolution of the overall number density of model galaxies is significantly weaker (i.e., larger number density at high redshift) than that of the present observation especially at low mass as in previous studies (e.g., Fontana et al. 2006; Kitzbichler & White 2007; Bertone et al. 2007). On the other hand, several studies found that semi-analytic models tend to predict the smaller number density of galaxies at the massive end (M_star ≳ 10^11.5–10^11.75 M_☉) than observations (e.g., Conselice et al. 2007; Bertone et al. 2007). A similar indication is also seen in Figure 18, although our data cannot strongly constrain the number density of these very massive galaxies.

Zoom In Zoom Out Reset image size

Since the stellar mass of a galaxy is generally dominated by long-lived low-mass stars, it can be generally considered as the integral of the past SFR and expected to increase monotonically with time. Therefore, for example, only suppression of star formation in low-mass galaxies cannot cause the flattening of the low-mass slope of the SMF with time, while it could lead to the evolution of the faint-end slope of the rest-frame UV luminosity function, which directly reflects the SFR, as discussed in Reddy & Steidel (2009). On the other hand, the stellar mass growth by star formation activities in each galaxy could change the slope of the SMF. There are several observational studies that suggest that low-mass galaxies have higher SFR relative to their stellar mass (i.e., SFR/M_star, specific SFR) than massive galaxies even at z∼ 2–3 (e.g., Feulner et al. 2005; Papovich et al. 2006; Reddy et al. 2006b, see also Santini et al. 2009). These results for the mass dependence of SFR/M_star suggest that the expected growth rate of stellar mass by star formation activities in a galaxy is higher in lower mass galaxies. If we assume the observed trends of the low-mass slope and mass dependence of SFR/M_star continue down to lower mass, a net increase of the number of galaxies in a mass bin due to star formation activities (inflow from lower mass bins − outflow into higher mass bins) is expected to be larger at lower stellar mass, which leads to the steepening of the low-mass slope with time, i.e., opposite to that seen in Figures 5 and 6. For example, Drory & Alvarez (2008) calculated the expected change of the SMF due to star formation only from the average SFR at each stellar mass of galaxies in the FORS Deep Field (Feulner et al. 2005) and showed that the number density of lower mass galaxies is expected to increase more rapidly with time especially at z≲ 3. By subtracting the expected contribution due to star formation from the observed redshift evolution of the SMF, Drory & Alvarez (2008) evaluated the effect of the hierarchical merging on the evolution of the SMF (i.e., destroying small galaxies and building more massive ones). In their analysis, at M_star ≲ 10¹⁰ M_☉, lower mass galaxies tend to be consumed in mergers at higher rate to counterbalance the effect of the mass-dependent SFR. Since no evolution of the low-mass end slope is assumed in Drory & Alvarez (2008), the steepening with redshift indicates the effect of merging larger than that estimated in Drory & Alvarez (2008). Therefore, the evolution of the SMF of galaxies in the MODS field suggests that the hierarchical merging process was important for the stellar mass assembly of these low-mass galaxies at z∼ 1–3, when the cosmic stellar mass density had been increasing rapidly. The relatively high fraction of peculiar/irregular morphology in low-mass galaxies at z > 1 in the NIR wavelength (Kajisawa & Yamada 2005; Conselice et al. 2005) might also imply intense merging activity of these galaxies.

In this context, the fact that low-mass galaxies tend to have relatively younger stellar age and higher SFR/M_star than massive galaxies at low to intermediate redshift (i.e., "downsizing," e.g., Kauffmann et al. 2003; Heavens et al. 2004) might not necessarily mean relative scarcity of low-mass galaxies at high redshift. A significant fraction of old stars in relatively massive (e.g., ∼10¹¹ M_☉) galaxies at low redshift could be formed at high redshift in actively star-forming low-mass galaxies which had coalesced into more massive objects through the mergers. One of the means to address this is to investigate the distribution of the SFR and mean stellar age of galaxies as a function of stellar mass down to low mass at high redshift. The combination of the distribution of the SFR and mean stellar age as a function of stellar mass and the SMF shown in this study allows us to measure the contributions of galaxies in a certain mass range to the cosmic SFR density and stellar mass density and their age distribution as a function of redshift. We will be able to use these quantities in order to estimate what the mass range of high-z galaxies could have assembled into more massive galaxies at later epochs.