The Rest-frame Optical (900 nm) Galaxy Luminosity Function at z ∼ 4–7: Abundance Matching Points to Limited Evolution in the M_STAR/M_HALO Ratio at z ≥ 4

Our LF measurements are based on a composite sample of galaxies at $4\lesssim z\lesssim 7$, selected in the rest-frame optical (z' band, ${\lambda }_{\mathrm{eff}}\sim 0.9\,\mu {\rm{m}}$; see Figure 1 for its transmission efficiency), with the bulk of our sample formed by Lyman break galaxies (LBGs) from the CANDELS/GOODS-N, CANDELS/GOODS-S—ERS fields. The z ∼ 4 bin is complemented by a sample of galaxies from a catalog based on the UltraVISTA DR2 mosaics.

Zoom In Zoom Out Reset image size

The LBG samples in the CANDELS GOODS-N/S and ERS fields rely on the multi-wavelength photometric catalogs of Bouwens et al. (2015). These are based on the re-reduction of public HST imaging, and are enhanced by proprietary full-depth Spitzer/IRAC mosaics. Specifically, they benefit from novel full depth IRAC 5.8 μm and 8.0 μm mosaics, not available in the original catalog of Bouwens et al. (2015). The UltraVISTA DR2 catalog is based on the most recent publicly available mosaics at UV-to-NIR wavelengths (including UltraVISTA DR2 data sets), and complemented by an internal release of full-depth Spitzer/IRAC mosaics.

In the following sections, we briefly describe these two parent catalogs and detail the criteria we adopted to assemble our final sample of galaxies.

2.1.1. GOODS-N/S and ERS

2.1.2. UltraVISTA DR2

The first step consists of applying a cut in the K_s flux of the galaxies from the UltraVISTA sample, in order to control the detection completeness. The DR2 data is ∼1 mag deeper than DR1. Therefore, we set the threshold to ${K}_{{\rm{s}}}=24.4$ mag, corresponding to the 90% completeness in detection of point-sources (Muzzin et al. 2013a). Instead, we did not apply any cut in the detection band of the GOODS-N/S sample, as the method adopted to estimate the co-moving volumes already takes into account the incompleteness from the detection stage.

Successively, we excluded from our sample those galaxies with poor flux measurements in those IRAC bands used to compute the rest-frame z' luminosity. The variation in depth across each IRAC mosaic prevented us from applying a single value of flux threshold at this stage. Instead, we applied a cut in S/N to the flux in the IRAC band closest to the rest-frame z' band (i.e., IRAC 4.5 μm at z ∼ 4 and IRAC 5.8 μm and 8.0 μm at $z\gtrsim 5$).

Considering the gap in the photometric depth probed by the 4.5 μm mosaics compared to that reached by the 5.8–8.0 μm data, we opted for applying a distinct S/N cut depending on the considered redshift bin. The sample at z ∼ 4 was selected by applying the cut of ${\rm{S}}/{\rm{N}}\gt 5$ to the 4.5 μm flux; the samples at z ∼ 5, 6, 7 were assembled by considering the cumulative flux in the 5.8 and 8.0 μm bands as the inverse-variance weighted sum of the flux in these two bands. We then applied a cut to the corresponding S/N such that:

$$\begin{eqnarray}&&{\rm{S}}\,{/{\rm{N}}}_{5.8+8.0}\equiv \displaystyle \frac{{{\rm{S}}}_{5.8}{w}_{5.8}+{{\rm{S}}}_{8.0}{w}_{8.0}}{\sqrt{{w}_{5.8}+{w}_{8.0}}}\gt 4,\end{eqnarray} \tag{ 1 }$$

where ${S}_{\kappa }$ is the flux measurement in band κ and w_κ is the weight defined as $1/{\sigma }_{\kappa }^{2}$, with ${\sigma }_{\kappa }$ the corresponding flux uncertainty. The application of the S/N cut reduced the number of galaxies to 2644 (2040/604 for GOODS/UltraVISTA, respectively), 96, 17 and 4 at z ∼ 4, z ∼ 5, z ∼ 6, and z ∼ 7, respectively.

We further cleaned our sample, excluding those objects satisfying any of the following conditions: (1) the contribution to the 5.8 and 8.0 μm flux from neighboring objects is excessively high; (2) the source morphology is very uncertain or confused making IRAC photometry undetermined; (3) the source is detected at X-rays wavelengths, suggesting it is a lower redshift AGN; (4) the source is at higher redshift, but its SED is dominated by AGN light; (5) LBGs with a likely $z\lt 3.5$ solution from photometric redshift analysis. In Appendix A, we detail our application of these additional criteria in cleaning our sample.

Our final sample consists of 2098 galaxies at z ∼ 4 (1680 from the LBG sample and 418 from the UltraVISTA sample), 72 at z ∼ 5, 10 at z ∼ 6, and 2 objects at z ∼ 7. The distribution of the absolute magnitudes in the z' band for the sample is presented in Figure 2, for the four different redshift bins. It is noteworthy how the GOODS-N/S and UltraVISTA samples complement each other at z ∼ 4, allowing to fully exploit these data with little redundancy.

Zoom In Zoom Out Reset image size

The samples adopted in this work rely on LBG selection criteria, complemented at z ∼ 4 by a photometric redshift selected sample based on the UltraVISTA DR2 catalog.

The criteria adopted for the assembly of our samples introduce two potential biases to our estimates of LF, ${M}_{* }/L$ ratio and SMF (Fontana et al. 2006; Grazian et al. 2015). The Lyman Break criteria select, by construction, blue star-forming galaxies, and may thus exclude a greater fraction of red objects compared to photometric-redshift selections. Furthermore, even samples based on photometric redshifts can suffer incompleteness from very red sources, too faint to appear in the detection bands (usually H or K_s or a combination of NIR bands), but that emerge at redder wavelengths (e.g., IRAC). In the following we attempt to evaluate the impact of these biases on our sample.

From the stellar mass catalog of CANDELS/GOODS-S (Santini et al. 2015), we extracted those galaxies with photometric redshift $3.5\lt {z}_{\mathrm{phot}}\lt 4.5$. Successively, we applied the LBG criteria to their flux measurements from the multi-wavelength catalog of Guo et al. (2013). In Figure 3 we present, as a function of stellar mass, the ratio between the number of galaxies recovered through the LBG criteria and the number of galaxies in the photo-z sample. The plot shows that the LBG selection is able to recover $\gtrsim 75 \%$ of galaxies with stellar mass $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10$ (corresponding to $\sim -0.12$ dex offset in number density measurements); at $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10.5$ the galaxies missed by the LBG criteria amount to about 35% (∼0.2 dex). At higher stellar masses, the fraction of galaxies not entering the Lyman Break selection increases to $\gtrsim 60 \% \mbox{--}70 \%$ ($\sim 0.5\mbox{--}0.6$ dex) consistent with Grazian et al. (2015). Duncan et al. (2014) showed that photometric uncertainties scatter a large fraction of the measurements outside the LBG selection box; specifically, the LBG criteria recover only $\sim 1/4$ (equivalent to $\sim -0.6$ dex offset) of the galaxies recovered through photo-z (see also Dahlen et al. 2010). However, once selection criteria on the redshift probability distribution are introduced, excluding from the sample poorly constrained photometric redshifts, the resulting photo-z sample largely overlaps with the LBG one, as demonstrated by the fact that the resulting photometric redshift UV LFs agree well, usually within $1-\sigma$, with the LBG UV LF (Duncan et al. 2014; Finkelstein et al. 2015a).

Zoom In Zoom Out Reset image size

The photometric depth of the UltraVISTA DR2 catalog, ${K}_{{\rm{s}}}=24.4$ mag, corresponds to a stellar mass completeness limit for a passively evolving simple stellar population of $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10.6$ at z ∼ 4 and 11.2 at z ∼ 5. The depth in the GOODS-Deep fields correspond to limits in stellar mass of $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10.3$ at z ∼ 4 and 10.6 at z ∼ 5, respectively. We would like to remark that our analysis for the z ∼ 4 sample at stellar masses $\mathrm{log}({M}_{* }/{M}_{\odot })\gtrsim 10\mbox{--}10.5$ is dominated by the photometric redshift sample from UltraVISTA, covering a larger volume for bright sources than the GOODS fields. Therefore, our composite z ∼ 4 sample is only marginally affected by the LBG selection bias.

A number of works have studied the so called extremely red objects, characterized by very red ($\gtrsim 2\mbox{--}3$ mag) rest-UV/optical colors, making them more elusive in high-z samples (e.g., Yan et al. 2004; Huang et al. 2011; Caputi et al. 2012, 2015; Stefanon et al. 2015; Wang et al. 2016). Samples detected in IRAC bands suggest that many of these objects could be consistent with being $z\gtrsim 3$ massive galaxies.

Recently, Wang et al. (2016) analyzed the properties of $H-[4.5]\gt 2.25$ mag over the CANDELS/GOODS-N and GOODS-S fields. Interestingly they identified 18 sources not present in the ${H}_{160}$-band catalog, but included in the IRAC catalog of Ashby et al. (2013). Of these, 5 sources have an estimated photometric redshift $3.5\lt {z}_{\mathrm{phot}}\lt 4.5$ and have a stellar mass $10.5\lesssim \mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 11$. Since their analysis refers to the same fields we consider in our work (although likely the configurations at the detection stage are different), we can use their result to obtain a rough estimate of the fraction of objects missed by our selection. Assuming a ${M}_{* }/{L}_{z^{\prime} }\sim 0.2{M}_{\odot }/{L}_{\odot }$, quite typical for these masses and redshifts (as we show in Section 3.2), the stellar mass range of these galaxies would correspond to luminosities $-24.7\lesssim {M}_{z^{\prime} }\lesssim -23.4$ AB. This sample would then constitute $\sim 65 \%$ of the objects in our z ∼ 4 LBG sample with similar luminosities. This fraction drops to $\sim 8 \%$ when comparing the 5 sources to the ∼60 galaxies with the same photometric redshift and stellar mass over the CANDELS/GOODS fields.

Caputi et al. (2015) presented SMF measurements at $z\sim 3\mbox{--}5$ obtained complementing the SMF from a photometric-redshift, K_s-selected sample based on UltraVISTA DR2 data to SMF measurements from photometric redshift samples of K_s-dropouts detected in IRAC bands. The main result is that K_s dropouts can account for as high as ∼0.5 dex in number densities. However, Stefanon et al. (2015) showed that samples similar to those of Caputi et al. (2015) likely suffer from degeneracies in the measurement of photometric redshifts (and consequently stellar masses), and therefore the above estimate is still uncertain.

The depth of current IRAC data, however, is not sufficient to systematically inspect passively evolving stellar population with stellar masses below $\sim {10}^{10-10.5}{M}_{\odot }$ at $z\gtrsim 4$. We therefore caution the reader that any sample currently available dealing with stellar mass below the $\sim {10}^{10}{M}_{\odot }$ limit may still be biased against dusty and/or old galaxies. Forthcoming projects, like Spitzer/GREATS (I. Labbé et al. 2017, in preparation) and the JWST will allow us to obtain a more complete picture by probing the lower mass regime.

We implemented a Monte Carlo simulation based on real data to estimate the effects that our selection criteria in S/N and contamination polishing have on the sample of galaxies used in this work. For this simulation, we did not consider the effects of selection in the detection band, because the UltraVISTA sample is 90% or more complete in K_s by construction, whereas the effects of detection completeness in the GOODS-N/S sample have been taken into account when estimating the co-moving volumes adopted for the LF measurements.

We first defined a grid in apparent magnitude of width 0.20 mag. Given the small sizes of the galaxies compared to the IRAC PSF, for each magnitude value in the grid, we injected 100 point sources randomly distributed across a region of uniform depth in the 4.5, 5.8, and 8.0 μm mosaics of the GOODS-N field. We chose the GOODS-N because this field is characterized by the deeper Spitzer/IRAC 4.5, 5.8, and 8.0 μm band data among the fields considered for this work. We successively replicated the flux measurement, using the same procedure adopted for the actual photometry.

The completeness fraction in each magnitude bin was computed by comparing the number of objects satisfying our selection criteria (Section 2.2) to the number of objects initially injected into the simulation. For the completeness of the z ∼ 4 sample, the above process was applied to the 4.5 μm mosaic only. For the completeness of the samples at $z\geqslant 5$, the point sources were added at matching positions in the 5.8 and 8.0 μm mosaics. The selection of the S/N and contamination was then recovered by applying the corresponding criteria and assuming the SED to be flat in f_ν in the observed 5–9 μm region. This is a reasonable approximation because, as we show in Section 3.2, the median SEDs do not substantially deviate from a flat f_ν SED in the wavelength range covered by IRAC observations. The whole process was repeated $10\times$ in each band, to increase their statistical significance. The global completeness (i.e., the cumulative effects of S/N and contamination selection) at the different depths of the IRAC mosaics was obtained by rescaling the completeness in S/N selection to match the depth of the relevant region.

The results from our completeness simulation for the 4.5 μm and the 5.8 μm + 8.0 μm samples are presented in Figure 4. We first discuss the recovery of the contamination fraction; we then consider these results in the budget of our global completeness estimates.

Zoom In Zoom Out Reset image size

The fraction of objects in the 4.5 μm band contaminated by neighbors⁶ is negligible for objects brighter than ∼24 AB, and increases exponentially up to ∼27 AB, where it starts to flatten out. A similar behavior is observed for the 5.8 μm + 8.0 μm simulation, although shifted at brighter magnitudes due to the shallower depth of the 5.8 and 8.0 μm compared to the 4.5 μm. The flattening at the faint end is caused by a strong incompleteness in the data at such faint magnitudes, and likely does not reflect the true behavior. In what follows and our analysis, we do not consider the completeness for magnitudes fainter than those corresponding to the onset of the flattening, i.e., ∼27 AB and ∼25 AB for the 4.5 μm and 5.8 μm + 8.0 μm data, respectively.

As it could intuitively be expected, the bright end of the global completeness curve is dominated by the (small) fraction of purged objects. This effect becomes less and less pronounced at fainter magnitudes, corresponding to lower S/N, where the effective selection is driven by the S/N itself.

The random noise from the background can positively combine at the location of a given source, introducing an increase in the measured flux (flux boosting, Eddington 1913). The amount of this boost is inversely correlated to the S/N of the source. The flux for sources with very high S/N will mostly be the result of the photons emitted by the source itself, with reduced contribution from the background; on the other hand, for sources with low S/N, the background level can be typically just few factors smaller than the intrinsic signal from the source, making it sensitive to (positive) fluctuations of the background. Furthermore, sources do not uniformly distribute with flux, but rather follow an approximate power law, with fainter sources more numerous than brighter ones. Therefore, it is intrinsically more probable that fainter sources scatter to brighter fluxes than the reverse, giving origin to a net flux bias.

A second potential source of flux boosting comes from confusion noise: faint sources at apparent positions close to a brighter one are more likely to be blended into the brighter source, increasing the flux and decreasing the number of fainter objects. This effect is larger for flux measurement in those bands with wide PSF, like Spitzer/IRAC. However, in our case, the photometry in the IRAC bands was performed by adopting a higher-resolution morphological prior from HST mosaics (see Section 2.1). Furthermore, we applied a selection in flux contamination (see Section 2.2). Because these two factors drastically limit the potential contribution of confusion noise to the IRAC fluxes in our sample, we do not further consider its effects to the flux-boosting budget.

For each source, we estimated the flux bias as the ratio between the expected intrinsic flux f_i and the measured flux f_o. Because no direct measurement is possible, the intrinsic flux was recovered as the average flux obtained from an estimate of its probability distribution. This was constructed considering two distinct contributions: (1) the probability ${ \mathcal P }({f}_{i},{f}_{o})$ that the observed flux f_o is drawn from the distribution of intrinsic flux f_i, given the noise ${\sigma }_{o}$, and (2) the frequency ${ \mathcal P }({f}_{i})$ of occurrence of the intrinsic flux f_i. Assuming each probability is normalized to 1, the final probability distribution would then be ${{ \mathcal P }}_{\mathrm{tot}}({f}_{i})\,={ \mathcal P }({f}_{i},{f}_{o}){ \mathcal P }({f}_{i})$.⁷

Assuming a Gaussian noise, ${ \mathcal P }({f}_{i},{f}_{o})$ can be written as:

$$\begin{eqnarray}&&{ \mathcal P }({f}_{i},{f}_{o})=\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{o}^{2}}}\exp [-({f}_{i}-{f}_{o})/2{\sigma }_{o}^{2}],\end{eqnarray} \tag{ 2 }$$

normalized to a total probability of one.

The frequency associated with the intrinsic flux can be recovered from the (intrinsic) differential number count of sources, ${dN}({f}_{i})/{{df}}_{i}$. This can usually be described by a power-law form with negative index, which is thus divergent for ${f}_{i}\to 0$, preventing it from being normalized⁸ (see, e.g., Hogg & Turner 1998, who also discuss possible reasons why the divergence at ${f}_{i}\to 0$ is likely non-physical).

We therefore followed the formalism of Crawford et al. (2010), who introduced, as further constraint, the Poissonian probability that no sources brighter than f_i exist at the same location of the observed object. The expression for ${ \mathcal P }({f}_{i})$ then becomes:

$$\begin{eqnarray}&&{ \mathcal P }({f}_{i})=\displaystyle \frac{{dN}({f}_{i})}{{{df}}_{i}}\times \exp \left(-{\rm{\Delta }}{{\rm{\Omega }}}_{o}{\int }_{{f}_{i}}^{+\infty }{dN}/{df}\right),\end{eqnarray} \tag{ 3 }$$

where ${\rm{\Delta }}{{\rm{\Omega }}}_{o}$ corresponds to the area occupied by the source. In the left panel of Figure 5, we show examples of reconstructed ${({ \mathcal P })}_{\mathrm{tot}}$, for the cases of S/N = 6 and S/N = 3.8, where it is evident that the contribution of the faint source population to the expected intrinsic flux increases as the S/N decreases. The right panel of Figure 5 shows the expected flux boost as a function of S/N. For ${\rm{S}}/{\rm{N}}\gtrsim 4.5$, the flux boost is roughly the same amount as the flux uncertainty. However, for lower S/N, the estimated flux boost increases abruptly. For S/N ≲ 2, the expected flux boost is $\gtrsim 1.5$ mag, meaning that the recovery of the intrinsic flux for such low S/N data is highly uncertain.

Zoom In Zoom Out Reset image size

The S/N in the HST bands for the galaxies in our sample is $\gt 10$. At z ∼ 4, the S/N in the 4.5 μm band, adopted for the selection of the z ∼ 4 sample, is $\geqslant 5$ by construction; the S/N in the 3.6 μm band is $\geqslant 5$ as well, consistent with the nearly flat SEDs in that wavelength range. At $z\geqslant 5$, the selection was performed in a combination of 5.8 μm and 8.0 μm fluxes, adopting a ${\rm{S}}/{\rm{N}}\gt 4$ threshold. Figure 5 shows that the expected flux boost for ${\rm{S}}/{\rm{N}}\gt 5$ is $\lesssim 0.1$ mag. However, for lower S/N, typical of the selection of samples at $z\geqslant 5$, the correction can be as high as ∼0.8–1.0 mag.

We therefore applied the above correction to the fluxes in the 5.8 and 8.0 μm of the $z\geqslant 5$ samples. The average flux boost was ∼0.19 mag and ∼0.25 mag in the IRAC 5.8 μm and 8.0 μm bands, respectively.

In Appendix B, we present the SEDs of the 12 most luminous galaxies in the z ∼ 5 sample, along with the SEDs of the z ∼ 6 and z ∼ 7 samples, before and after applying the flux boost correction.

In the last few years, a number of works have studied the relation between the rest-frame UV luminosity and the stellar mass of high-redshift galaxies (e.g., Stark et al. 2009; González et al. 2011; Lee et al. 2011; McLure et al. 2011; Duncan et al. 2014; Spitler et al. 2014; Grazian et al. 2015, V. González et al. 2017, in preparation; Song et al. 2016). Indeed, a relation between the stellar mass and the UV luminosity is to be expected when considering a continuous star formation. Deviations from such a relation would then provide information on the age and metallicity of the stellar population and dust content of the considered galaxies.

The emerging picture is that, at z ∼ 4 and for stellar masses $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10$, the stellar mass increases monotonically with increasing UV luminosity; however, at stellar masses higher than $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10$, the trend becomes more uncertain: Spitler et al. (2014), using a sample of K_s-based photometric redshift selected galaxies, found indication of a turnover of the UV luminosity, with the more massive galaxies ($10.5\lesssim \mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 11$) spanning a wide range in UV luminosities (see also Oesch et al. 2013). Lee et al. (2011), on the other hand, using an LBG sample, found a linear relation between UV luminosities and stellar masses up to $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 11$. Considering the different criteria adopted by the two teams for the assembly of their samples, selection effects might be the main reason for the observed tension.

This observed discrepancy could, however, just be the tip of the proverbial iceberg. Indeed, current high-z surveys might still be missing lower-mass, intrinsically red galaxies (dusty and/or old), which could populate the ${M}_{* }\mbox{--}{M}_{\mathrm{UV}}$ plane outside the main sequence (Grazian et al. 2015 and our discussion in Section 2.3). The depth of the current NIR surveys does not allow us to further inspect this, which will likely remain an open issue until JWST.

A monotonic relation between the UV luminosity and the stellar mass has also been found at z ∼ 5 and z ∼ 6 for $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10$ (e.g., Stark et al. 2009; González et al. 2011; Duncan et al. 2014; Salmon et al. 2015; Song et al. 2016), with approximately the same slope and dispersion, but with an evolving normalization factor (but see, e.g., Figure 5 of Song et al. 2016 for further hints as to the existence of massive galaxies with faint UV luminosities).

Figure 6 presents the absolute magnitude in the rest-frame z' band (${M}_{z^{\prime} }$) as a function of the absolute magnitude in the UV (${M}_{{\mathrm{UV}}_{1600}}$), for our sample in the four redshift bins (z ∼ 4, 5, 6 and 7). In Figure 7, we present the binned median in the ${M}_{\mathrm{UV}}\mbox{--}{M}_{z^{\prime} }$ plane for the z ∼ 4 and z ∼ 5 samples.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The z ∼ 4 sample shows a clear correlation between the luminosities in the rest-frame UV and z' bands for ${M}_{z^{\prime} }\gtrsim -22$ mag, which can be described by the following best-fitting linear relation:

$$\begin{eqnarray}&&{M}_{\mathrm{UV}}=(-3.58\pm 1.49)+(0.79\pm 0.07)\times {M}_{z^{\prime} }.\end{eqnarray} \tag{ 4 }$$

The above best fit is marked by the magenta line in Figure 7, where we also indicate the $3\sigma$ limits corresponding to our 5.8 μm + 8.0 μm selection. At z ∼ 4 and z ∼ 5, the depth of the IRAC data allows us to not only probe the bright end, where the relation between M_UV and ${M}_{z^{\prime} }$ breaks, but to also explore the regime of the linear correlation expressed by Equation (4). Slopes of 0.4–0.5 in the $\mathrm{log}({M}_{* })\mbox{--}{M}_{\mathrm{UV}}$ plane (with nominal $1\sigma$ uncertainties of $\sim 0.05\mbox{--}0.1$) have been reported by, e.g., Duncan et al. (2014) and Song et al. (2016). Assuming a constant ${M}_{* }/{L}_{z^{\prime} }$ ratio (see Section 3.2), our measurements correspond to a slope of 0.5 in the $\mathrm{log}({M}_{* })\mbox{--}{M}_{\mathrm{UV}}$ plane, consistent with previous measurements.

Assuming that the SFR mostly comes from the UV light, and that ${M}_{z^{\prime} }$ is a good proxy for stellar mass measurements, we can also compare the slope we derived for the ${M}_{\mathrm{UV}}\mbox{--}{M}_{z^{\prime} }$ relation to that of the log(SFR)–$\mathrm{log}({M}_{* })$ from the literature. Indeed, the observed UV slopes of ${M}_{z}^{\prime} \gtrsim -22$ galaxies in our sample are $\beta \sim -2$ (see also Bouwens et al. 2010), consistent with star-forming galaxies and little-to-no dust extinction. Our measurement is perfectly consistent with the log(SFR)–$\mathrm{log}({M}_{* })$ slope of $\sim 0.8\pm 0.1$ recently measured by, e.g., Whitaker et al. (2015) for $z\lesssim 2.5$ star-forming galaxies with $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10.5$, and it has been shown to evolve little over the redshift range $0.5\lt z\lt 2.5$.

For absolute magnitudes brighter than ${M}_{z^{\prime} }\sim -22$ mag, the linear relation expressed by Equation (4) breaks as we observe the beginning of a turnover in the absolute UV–z' magnitude relation. Remarkably, this behavior is visible also for the z ∼ 5 sample, which is entirely based on LBG selection. This fact has important consequences for, e.g., SMF measurements; samples of galaxies selected at fixed rest-frame UV luminosity are potentially characterized by a wide range of stellar mass.

The absolute UV magnitude of galaxies with ${M}_{z^{\prime} }\sim -23.5\,\pm 0.8$ mag spans the full range of values observed for ${M}_{z^{\prime} }\lesssim -22.5$ mag. However, the bulk of values aggregates around the $[{M}_{\mathrm{UV}},{M}_{z^{\prime} }]\sim [-21.4,-23.5]$ mag region, and it is characterized by a large dispersion in M_UV ($\gtrsim 3$ mag). This result is qualitatively consistent with the findings of Spitler et al. (2014), assuming a correlation between the absolute magnitude ${M}_{z^{\prime} }$ and the stellar mass. Most of the galaxies with ${M}_{z^{\prime} }\sim -23.5$ come from the photometric redshift sample selected from the UltraVISTA catalog. As we will present in Section 3.2, our measurements of the mass-to-light ratios from stacking analysis show that galaxies with ${M}_{z^{\prime} }\lesssim -23.5$ statistically have stellar masses $\mathrm{log}({M}_{* }/{M}_{\odot })\gtrsim 10.6$. The above result then underlines the bias that LBG selections may introduce against massive systems.

At z ∼ 5, our data allow us to inspect the relation only for ${M}_{z^{\prime} }$ fainter than −23 and ${M}_{\mathrm{UV},1600}$ fainter than $\sim -22$. In this range of luminosities, our z ∼ 5 measurements are roughly consistent with the z ∼ 4 measurements in the same range of luminosities. The measurements for the z ∼ 6 and z ∼ 7 samples are still consistent with the trends observed at z ∼ 4, although the low number of objects does not allow us to derive any statistically significant conclusion.

Thus far, determinations of the stellar mass-to-light ratios (${M}_{* }/L$) for galaxies at $z\gt 4$ have involved the ${M}_{* }/{L}_{\mathrm{UV}}$ ratio. This quantity is fundamental to our understanding of galaxy formation and evolution, as it combines information on the recent (through the UV luminosity) and integrated (through the stellar mass) SFH (e.g., Stark et al. 2009). Nonetheless, the ${M}_{* }/{L}_{\mathrm{UV}}$ has been used to recover the stellar mass and SMF of high redshift galaxies with alternating success (see, e.g., González et al. 2011; Song et al. 2016). In this section, instead, we explore for the first time the ${M}_{* }/{L}_{z^{\prime} }$ properties of galaxies at $z\gtrsim 4$. The rest-frame z' luminosity is more sensitive to the stellar mass, compared to the UV luminosity, for two reasons. Although the UV light is emitted by massive, short-lived stars, and thus traces the SFH in the past few hundred Myr, the luminosity in the rest-frame optical region mostly originates from lower-mass, longer-lived stars, which constitute most of the stellar mass of galaxies. Furthermore, it is less sensitive to the dust extinction, and hence to the uncertainties in its determination, compared to the UV; for a Calzetti et al. (2000) extinction curve, an A_V = 1 mag gives ${A}_{\lambda 1600}\sim 2.5$ mag, compared to ${A}_{\lambda 9000}\sim 0.5$ mag.

Because we are interested more on average trends in the ${M}_{* }/{L}_{z^{\prime} }$ ratios here, rather than studying it for specific galaxies, we performed our analysis using the median stacked SEDs constructed from our composite sample. Due to the different photometric bands in the catalogs, we performed the stacking of sources separately for sources in the GOODS-N/S and UltraVISTA samples.

The stacked SEDs were constructed as follows (see also González et al. 2011). At each redshift interval, we divided the galaxies into sub-samples according to their ${M}_{z^{\prime} }$. The different depths reached by the 4.5 μm and 5.8 μm + 8.0 μm samples resulted in different numbers of subsamples across the redshift bins. Under the working assumption of limited variation in both redshift and SED shape in each bin of ${M}_{z^{\prime} }$, as well as for each HST band, we took the median of the individual flux measurements. Our assumption is also supported by the fact that the SEDs from stacking are generally characterized by a flat f_ν continuum at both rest-frame UV and optical regimes. Uncertainties on the median were computed from bootstrap techniques, drawing with replacement the same number of flux measurements as the number of galaxies in the considered absolute magnitude bin. Before median-combining, the fluxes were perturbed according to their associated uncertainty. The process was repeated 1000 times, and the standard deviation of the median values was taken as the final uncertainty. For the IRAC bands, median stacking was performed on the mosaic cutouts centered at the position of each source, previously cleaned from neighbors. Photometry was performed on the median of the images in apertures of $2\buildrel{\prime\prime}\over{.} 5$ diameter. Total fluxes were then recovered through the PSF growth curve. Uncertainties were computed by applying to the image cutouts the same bootstrap technique adopted for the median stacking of the fluxes, as described above. In randomly drawing the image cutouts, we preserved the total exposure time.

Photometric redshifts and z'-band luminosities were obtained from EAzY (Brammer et al. 2008) on the stacked SEDs. Briefly, EAzY initially selects the two SED templates that provide the closest match to the observed color in the two filters bracketing the rest-frame band of interest. The luminosity is then computed from the interpolation of the two colors, relative to the rest frame band of interest, obtained from the two selected SEDs (for full details, see Appendix C of Rudnick et al. 2003). Stellar masses were computed by running FAST (Kriek et al. 2009), adopting the Bruzual & Charlot (2003) template SEDs, a Chabrier (2003) IMF, solar metallicity, and a delayed-exponential SFH. The bands potentially contaminated by nebular emission were excluded from the fit. Because we performed the stacking in each band individually, assuming the same redshift for all sources, the flux in those bands close to the Lyman and the Balmer breaks potentially suffers from high scatter introduced by the range of redshifts of the galaxies in each sub-sample, depending on whether the break enters the band. Fluxes in these bands were therefore excluded from the fit with FAST. Specifically, for the z ∼ 4 LBG stacks, we excluded the ${B}_{435},{V}_{606}$ and H₁₆₀ bands; for the z ∼ 4 UltraVISTA stacks, we excluded the B, IA427, IA464, IA484, IA505, IA527, IA574, IA624, IA679, $g^{\prime} ,{g}^{+},V,H$ bands; for the z ∼ 5 stacks, we excluded the V₆₀₆ and the i₇₇₅ bands; for the z ∼ 6 stack, we excluded the i₇₇₅ band; for the z ∼ 7 stack, we excluded the I₈₁₄ band. The stacked SEDs, along with the best-fit templates from FAST, are presented in Figure 8.

Zoom In Zoom Out Reset image size

The photometric redshifts measured from the stacked SED are all consistent with the values of the corresponding redshift bin; the difference of the photometric redshifts of the median stacked SEDs and the median of the photometric redshifts of the individual sources in each subsample is ${\rm{\Delta }}z/(1+z)\lesssim 0.05$, i.e., within the uncertainties expected for photometric redshifts.

The inset in the left panel of Figure 8 presents the slope of the UV continuum (β), measured on the stacked photometry, as a function of absolute magnitude ${M}_{z^{\prime} }$. The stacked SEDs at z ∼ 4 and z ∼ 5 are characterized by a trend in the UV slope, with bluer slopes for low-luminosity galaxies, particularly evident for the z ∼ 4 stacks, and qualitatively consistent with the results of, e.g., González et al. (2011) and Oesch et al. (2013). Our measurements do not present evidence for evolution of β with redshift at fixed luminosity ($\approx {M}_{* }$), although the large uncertainties in β (especially for the higher redshift bins) may be blurring trends. Furthermore, the SEDs in the rest-frame UV wavelength region of the four or five brightest z ∼ 4 stacks do not differ too much one from each other, whereas they differ substantially at wavelengths redder than the Balmer/4000 Å break.

Recently, Oesch et al. (2013) presented stacked SEDs at z ∼ 4 in bins of z-band absolute magnitudes from a sample of galaxies based on CANDELS/GOODS-S, HUDF, and HUDF09-2. The sample benefits from deep IRAC 3.6 μm and 4.5 μm imaging from the IUDF program. The stacked SEDs (see, e.g., their Figure 2) show a clear trend of redder colors with increasing rest-frame z-band luminosity, in particular for ${M}_{z}\lt -21.5$. Our stacked SEDs are in qualitative agreement with those of Oesch et al. (2013), confirming the observed trend with luminosity. Furthermore, thanks to the wide area offered by UltraVISTA, which provides coverage for even brighter sources, we are able to extend the trend to even more luminous galaxies.

In the top panel of Figure 9, we present the ${M}_{* }/{L}_{z^{\prime} }$ measured⁹ from the stacked SEDs as a function of z' absolute magnitude. Total uncertainties were obtained by propagating the 68% confidence intervals in stellar mass generated by FAST and the uncertainties in luminosity, taken as the flux uncertainties from stacking. At z ∼ 4, the ${M}_{* }/{L}_{z^{\prime} }$ ratio is consistent with being constant for ${M}_{z^{\prime} }$ fainter than $\sim -22.5$ mag. We find:

$$\begin{eqnarray}&&\mathrm{log}({M}_{* }/{L}_{z^{\prime} })=-0.87\pm 0.09,\mathrm{for}\ \ {M}_{z^{\prime} }\gtrsim -22.5.\end{eqnarray} \tag{ 5 }$$

For ${M}_{z^{\prime} }\lesssim -22.5$ mag, there is indication of ${M}_{* }/{L}_{z^{\prime} }$ increasing with the luminosity, although the error bars are large. The best-fit SEDs of the most luminous stacks have a nearly constant age ($\sim {10}^{8.8}$ yr), and show an A_V slightly increasing with stellar mass (from 1.0 to 1.2 mag). A linear fit of the $\mathrm{log}({M}_{* }/{L}_{z^{\prime} })$ values for ${M}_{z^{\prime} }\lesssim -22.5$ mag resulted in the following relation:

$$\begin{eqnarray}&&\mathrm{log}({M}_{* }/{L}_{z^{\prime} })=(-5.1\pm 4.7)-(0.19\pm 0.18)\times {M}_{z^{\prime} },\\ &&\quad \mathrm{for}\ \ {M}_{z^{\prime} }\lesssim -22.5.\end{eqnarray} \tag{ 6 }$$

Our linear relation recovers a constant value of $\mathrm{log}({M}_{* }/{L}_{z^{\prime} })\,=-0.87$ at ${M}_{z^{\prime} }\sim -22.5$. However, the uncertainties on the fit parameters make the above relation also consistent with a constant value.

Zoom In Zoom Out Reset image size

The ${M}_{* }/{L}_{z^{\prime} }$ ratios for the z ∼ 5 and z ∼ 6 are consistent with ${M}_{* }/{L}_{z^{\prime} }$ measurements of the z ∼ 4, ${M}_{z}\gt -22.5$ AB stacks. The ${M}_{* }/{L}_{z^{\prime} }$ for the z ∼ 7 bin is consistent with the average ${M}_{* }/{L}_{z^{\prime} }$ only at the $\sim 3\sigma$ level. We note, however, that the z ∼ 7 sample only includes two sources, thus reducing the statistical significance of the observed disagreement.

The above results are consistent with what is observed in Figure 6. In Section 3.1, we showed that galaxies more luminous than ${M}_{z^{\prime} }\sim -22.5$ mag form a cloud in the rest-frame UV–z' plane around ${M}_{\mathrm{UV}}\sim -21.4$ mag. From our stacking analysis, the average apparent magnitude at ${\lambda }_{\mathrm{obs}}\sim 8000$ Å (i.e., the rest-frame UV₁₆₀₀) for the stacked SEDs with $[4.5\,\mu {\rm{m}}]\lt 23.5$ AB is ∼24.7 AB, which corresponds at z ∼ 4 to an absolute magnitude ${M}_{\mathrm{UV},1600}\sim -21.4$. According to the above relation, the stellar mass corresponding to ${M}_{z^{\prime} }\sim -23$ mag is $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10.3$. This behavior raises concern about potential biases that can occur when adopting the UV luminosity and ${M}_{* }/{L}_{\mathrm{UV}}$ in the measurement of stellar masses, particularly for massive galaxies.

A constant ${M}_{* }/L$ is equivalent to a slope of −0.4 in the log(stellar mass)—absolute magnitude plane. Our result at z ∼ 4, obtained for galaxies with ${M}_{z^{\prime} }\lt -23$ mag, is consistent with the $\sim -0.4$ slopes found in the stellar mass—M_UV plane (see, e.g., Duncan et al. 2014; Grazian et al. 2015). Steeper slopes, such as those found by Stark et al. (2009), González et al. (2011), Lee et al. (2011), McLure et al. (2011), and Song et al. (2016), require the $\mathrm{log}({M}_{* }/L)$ to decrease for fainter galaxies, increase for brighter galaxies, or a combination of both effects. The origin for this is still unclear, as it could be a mix between selection effects (see, e.g., Grazian et al. 2015) and nebular emission contamination, which could boost the stellar masses of the more luminous galaxies (e.g., Song et al. 2016).

In the bottom panel of Figure 9, we compare the measurements of stellar mass and luminosity recovered from our stacking analysis to a control sample. This sample is composed of individual measurements selected from the COSMOS/UltraVISTA and CANDELS/GOODS-S catalogs to have photometric redshifts $3.5\lt {z}_{\mathrm{phot}}\lt 4.5$.

The individual measurements of the control sample distribute according to a monotonic relation defining a main sequence, with a scatter of about ∼0.7 dex. This correlation holds over a wide range of values, both in stellar mass and luminosity. In particular, we notice the absence of any turnover, as instead is observed when considering ${M}_{\mathrm{UV}}\mbox{--}{M}_{* }$, e.g., Spitler et al. (2014) or equivalently, our Figure 6. Furthermore, there are virtually no measurements outside the main sequence.

The relation between the stellar mass and luminosity recovered from the stacking analysis is in excellent agreement with the values of the control sample. We remark here that the control sample was selected based exclusively on photometric redshifts criteria. The agreement between our stacking results in the GOODS field, then, indicates that the stacking does not suffer any major bias from the LBG selection. This result is not unexpected, however, given the low fraction of objects missed by the LBG selection for stellar masses $\mathrm{log}({M}_{* }/{M}_{\odot })\lesssim 10.2$, as we showed in Section 2.3.

Together, these two results increase our confidence in the reliability of the z' band as a proxy for the stellar mass, s well as the robustness of our stacking analysis.

The LFs were measured adopting the $1/{V}_{\max }$ estimator (Schmidt 1968). Although this method is intrinsically sensitive to local overdensities of galaxies, the clustering is expected to be negligible at $z\gt 4$. On the other hand, the $1/{V}_{\max }$ method directly provides the normalization of the LF. Furthermore, and most importantly, the coherent analysis extension developed by Avni & Bahcall (1980) is key to this work.

As we showed in Section 2.1, our composite sample is based on a dual-band flux selection, corresponding to a double flux threshold. The detection process introduces the first flux cut in the corresponding band (K_s or ${\chi }^{2}$ image built from the HST NIR bands, for the UltraVISTA and GOODS-N/S sample, respectively). The S/N cut on the flux in the IRAC band closest to the rest-frame z' is responsible for the second flux threshold in the relevant IRAC band.

For each galaxy in the sample, each flux threshold generates an upper limit to the redshift the specific galaxy can have and still be included in the sample. These different upper limits in redshift correspond to different comoving volumes for each object that could potentially enter the V_max computation. The coherent approach allowed us to take this double selection into consideration in a consistent way: the upper limit in redshift, used to compute the comoving volume, was taken to be the smaller of the two redshift upper limits, computed based on the threshold in the corresponding selection band. Furthermore, as we showed in Figure 1, the depth of the IRAC mosaics is highly inhomogeneous. Therefore, for the computation of the comoving volumes in each field, we divided the IRAC footprint into a number of sub-fields, such that each sub-field was characterized by nearly homogeneous depths in both the detection¹⁰ and relevant IRAC band. Again, the Avni & Bahcall (1980) prescription allowed us to coherently analyze the different sub-samples.

Comoving volumes were computed differently, depending on the field and on the band driving the selection. For the galaxies in the GOODS-N/S fields, we used the comoving volumes computed by Bouwens et al. (2015). These volumes were estimated using an extensive Monte Carlo simulation based on real data. Sources were added to the different mosaics and recovered following the same procedure applied for the assembly of the LBG sample. Such volume estimates natively take into account the selection effects at the detection stage, and correct for flux-boosting effects and contamination by lower redshift interlopers and brown dwarfs. The volumes V_i for those objects i in the GOODS-N/S fields whose redshift upper limit z_up was driven by the IRAC S/N threshold (${z}_{\mathrm{up}}\equiv {z}_{\mathrm{up},\mathrm{IRAC}}$) were rescaled by the ratio between the volume associated with the redshift upper limit of the IRAC band ${V}_{i}({z}_{\mathrm{up},\mathrm{IRAC}})$ and that of the ${\chi }^{2}$ image ${V}_{i}({z}_{\mathrm{up},{\chi }^{2}})$:

$$\begin{eqnarray}&&{V}_{i,\mathrm{IRAC}}={V}_{i,\mathrm{GOODS}}\times \displaystyle \frac{{V}_{i}({z}_{\mathrm{up},\mathrm{IRAC}})}{{V}_{i}({z}_{\mathrm{up},{\chi }^{2}})}.\end{eqnarray} \tag{ 7 }$$

For the UltraVISTA sample, the volumes were computed directly from the limits in redshift that correspond to the flux limits in the K_s and 4.5 μm bands.

We computed the LF in four redshift bins centered at z ∼ 4, z ∼ 5, z ∼ 6, and z ∼ 7. Although the IRAC data potentially allowed us to consider galaxies at z ∼ 8, we did not find any candidate with reliable flux measurement in the IRAC 5.8 and 8.0 μm. Uncertainties on the LF measurements were derived by combining in quadrature the Poisson noise in the approximation of Gehrels (1986) to an estimate of cosmic variance from the recipe of Moster et al. (2011). The average cosmic variance value obtained for the z ∼ 4 UltraVISTA sample was ∼0.43; the average cosmic variance estimates for the GOODS-N/S sample were $\sim 0.27,\sim 0.41,\sim 0.58$, and ∼0.80, respectively, for the z ∼ 4, 5, 6, and z ∼ 7 redshift bins. The high values of the cosmic variance registered for all redshifts and luminosities are the dominant source of stochastic uncertainties in our LF measurements.

Our LF measurements are presented in Figure 10 and Table 1. The LF at z ∼ 7 consists of a single-bin measurement, and is characterized by large uncertainties that do not allow us to properly constrain its shape. The absolute magnitude range of the z ∼ 4 LF spans ∼5 magnitudes, $\gtrsim 3\times$ more than the magnitude range of the $z\geqslant 5$ LFs. The larger absolute magnitude range available at z ∼ 4 is the result of a number of distinct factors. First, the increased depth in the 4.5 μm band from the combination of Spitzer/IRAC cryogenic and post-cryogenic epochs enables us to reach fainter absolute magnitudes than the cryogenic $5.8+8.0\,\mu {\rm{m}}$ data alone at $z\geqslant 5$. Second, the smaller PSF size of the 4.5 μm data, compared to the 5.8 μm and 8.0 μm bands, allows to reach fainter fluxes for the same exposure time and detector efficiency. Third, the availability of COSMOS/UltraVISTA data over an area $\sim 4\times$ larger than the GOODS-N/S footprint allowed us to recover the exponential decline of the bright end of the z ∼ 4 LF, which is otherwise inaccessible by the small footprint of the GOODS-N/S mosaics.

Zoom In Zoom Out Reset image size

Table 1.  V_max Measurements of the Luminosity Functions

z	${M}_{z^{\prime} }$	${\rm{\Delta }}{M}_{z^{\prime} }$	Φ	#
bin	(mag)	(mag)	(10−5 Mpc−3 mag−1)	gal.
z ∼ 4	−25.25	0.25	${0.044}_{-0.041}^{+0.102}$	1
	−25.00	0.25	${0.088}_{-0.068}^{+0.121}$	2
	−24.75	0.25	${0.22}_{-0.13}^{+0.18}$	5
	−24.50	0.25	${0.70}_{-0.35}^{+0.37}$	16
	−24.25	0.25	${1.27}_{-0.59}^{+0.62}$	29
	−24.00	0.25	${1.71}_{-0.78}^{+0.80}$	39
	−23.75	0.25	${2.4}_{-1.1}^{+1.1}$	55
	−23.50	0.25	${3.9}_{-1.7}^{+1.8}$	78
	−23.25	0.25	${6.6}_{-2.9}^{+3.0}$	100
	−23.00	0.25	${10.8}_{-4.8}^{+4.8}$	86
	−22.75	0.25	${16.8}_{-7.5}^{+7.6}$	64
	−22.50	0.25	$25{.}_{-11.}^{+11.}$	72
	−22.25	0.25	${28.6}_{-7.2}^{+7.5}$	64
	−22.00	0.25	$44{.}_{-11.}^{+11.}$	87
	−21.75	0.25	$62{.}_{-15.}^{+15.}$	123
	−21.50	0.25	$65{.}_{-15.}^{+16.}$	128
	−21.25	0.25	$95{.}_{-22.}^{+22.}$	186
	−21.00	0.25	$101{.}_{-23.}^{+24.}$	194
	−20.75	0.25	$103{.}_{-24.}^{+24.}$	196
	−20.50	0.25	$126{.}_{-29.}^{+29.}$	207
	−20.25	0.25	$162{.}_{-37.}^{+38.}$	193
	−20.00	0.25	$191{.}_{-45.}^{+46.}$	119
	−19.75	0.25	$215{.}_{-56.}^{+59.}$	49
z ∼ 5	−23.75	0.50	${0.56}_{-0.41}^{+0.76}$	2
	−23.25	0.50	${2.6}_{-1.2}^{+1.5}$	9
	−22.75	0.50	${4.6}_{-2.0}^{+2.2}$	16
	−22.25	0.50	${8.2}_{-3.4}^{+3.7}$	20
	−21.75	0.50	$36{.}_{-15.}^{+17.}$	17
	−21.25	0.50	$83{.}_{-41.}^{+50.}$	8
z ∼ 6	−23.40	0.80	${0.20}_{-0.19}^{+0.46}$	1
	−22.50	1.00	${3.2}_{-2.1}^{+2.9}$	4
	−21.65	0.70	$26{.}_{-16.}^{+21.}$	5
z ∼ 7	−22.25	0.50	${2.9}_{-2.4}^{+4.1}$	2

Download table as:  ASCIITypeset image

In order to verify the consistency of the z ∼ 4 LF with respect to the GOODS-N/S and UltraVISTA data, we also computed the LF separately on each one of these two data sets. The resulting LFs are marked in Figure 10 by the open squares, and show a good agreement with the LF from the composite sample.

The large uncertainties associated to the number density measurements at z ∼ 5, 6, and 7 do not allow us to disentangle whether the evolution is in luminosity, number density, or both. In Section 3.5, we attempt to analyze this in a more quantitative way.

The dashed gray curve in Figure 10 marks the $z\sim 0$ LF of Kelvin et al. (2014), measured with data from the Galaxy and Mass Assembly (GAMA) survey (Driver et al. 2009). Compared to our lowest redshift LF measurement, the $z\sim 0$ LF is characterized by a steeper decay for ${M}_{z^{\prime} }\lesssim -24.5$. Fully understanding the evolution of the LF from z ∼ 4 to $z\sim 0$ goes beyond the scope of the present work. However, qualitatively, a decrease in luminosity at $z\sim 0$ compared to z ∼ 4 is expected, considering the lower values of the star formation rate density at $z\sim 0$ than at z ∼ 4 (e.g., Madau & Dickinson 2014), and that the z' band may retain the effects of the recent SFH. This is particularly true for the bright-end of the LF: indeed, the brightest (i.e., most massive) galaxies that are still forming stars at z ∼ 4 are likely to become quenched by $z\sim 0$ (e.g., Muzzin et al. 2013b).

We generated SMF measurements by taking advantage of the mass-to-light ratios we measured in Section 3.2. The fact that ${M}_{* }/{L}_{z^{\prime} }$ does not decrease with luminosity makes the shape of the LF resemble that of the SMF, allowing us to attempt a simple and straightforward conversion of the LF into the SMF. Other M/L relations allow for the recovery of the SMF from the LF, although in a less straightforward way. We further discuss this in Section 4.1.

We adopt the following very simple procedure. We assume that the constant $\mathrm{log}({M}_{* }/{L}_{z^{\prime} })$ and the linear relation observed at z ∼ 4 (Equations (5) and (6)) are valid at all redshifts. The absolute magnitudes corresponding to the bin centers of the LFs are converted into stellar mass by applying the relevant $\mathrm{log}({M}_{* }/{L}_{z^{\prime} })$ relation, depending on the ${M}_{z^{\prime} }$ value (see Equations (5) and (6)). We then differentiate the two relations, solving for dM_*. The obtained values, specific for each M_* bin are used to rescale the LF normalization, to take into account the change in units from mag⁻¹ to dex⁻¹.

Our SMF measurements are presented in Figure 11 and Table 2. Unsurprisingly, the z ∼ 4 SMF covers a range in stellar mass wider than the z ∼ 5, 6, and 7 SMFs, for the same reasons we described for the LF.

Zoom In Zoom Out Reset image size

Table 2.  Stellar Mass Function Measurements

z	Stellar mass	Φ
bin	$\mathrm{log}({M}_{* }/{M}_{\odot })$	(10−5 Mpc−3 dex−1)
z ∼ 4	11.55	${0.075}_{-0.070}^{+0.174}$
	11.41	${0.15}_{-0.12}^{+0.21}$
	11.26	${0.37}_{-0.23}^{+0.30}$
	11.12	${1.20}_{-0.59}^{+0.64}$
	10.97	${2.2}_{-1.0}^{+1.0}$
	10.82	${2.9}_{-1.3}^{+1.4}$
	10.68	${4.1}_{-1.9}^{+1.9}$
	10.53	${6.7}_{-3.0}^{+3.0}$
	10.38	${11.3}_{-5.0}^{+5.0}$
	10.24	${18.5}_{-8.2}^{+8.3}$
	10.09	$29{.}_{-13.}^{+13.}$
	9.94	$61{.}_{-27.}^{+28.}$
	9.84	$71{.}_{-18.}^{+19.}$
	9.74	$109{.}_{-27.}^{+27.}$
	9.64	$155{.}_{-37.}^{+37.}$
	9.54	$162{.}_{-38.}^{+39.}$
	9.44	$238{.}_{-55.}^{+56.}$
	9.34	$252{.}_{-58.}^{+59.}$
	9.24	$258{.}_{-60.}^{+60.}$
	9.14	$316{.}_{-73.}^{+73.}$
	9.04	$405{.}_{-94.}^{+94.}$
	8.94	$477{.}_{-114.}^{+115.}$
	8.84	$537{.}_{-141.}^{+147.}$
z ∼ 5	10.68	${0.95}_{-0.70}^{+1.29}$
	10.38	${4.4}_{-2.1}^{+2.5}$
	10.09	${7.8}_{-3.4}^{+3.7}$
	9.84	${20.5}_{-8.5}^{+9.2}$
	9.64	$91{.}_{-39.}^{+42.}$
	9.44	$208{.}_{-102.}^{+125.}$
z ∼ 6	10.47	${0.34}_{-0.32}^{+0.79}$
	9.94	${8.1}_{-5.3}^{+7.4}$
	9.60	$65{.}_{-40.}^{+53.}$
z ∼ 7	9.84	${7.4}_{-6.0}^{+10.4}$

Download table as:  ASCIITypeset image

In Figure 11, we also plot a compilation of SMF measurements from the literature (Pérez-González et al. 2008; Marchesini et al. 2009, 2010; Stark et al. 2009; González et al. 2011; Lee et al. 2012; Santini et al. 2012; Ilbert et al. 2013; Muzzin et al. 2013b; Duncan et al. 2014; Stefanon et al. 2015; Caputi et al. 2015; Grazian et al. 2015; Song et al. 2016; Davidzon et al. 2017). At z ∼ 4, starting from the low-mass end where the measurements are generally quite consistent with each other, the discrepancies increase with increasing stellar mass. One possible reason for the increased dispersion at higher masses is that galaxies constituting the low-mass end are mostly star-forming. Their redshift can then be assessed through the location of the observed Lyman break (either from dropouts or photometric redshift selections). The massive end, instead, possibly also includes more evolved and/or dusty systems, and it is therefore more sensitive to the degeneracy in identifying the observed break as either the Balmer/4000 Å break or the Lyman break.

Our z ∼ 4 SMF determination is in good agreement with the SMFs of Stark et al. (2009), Lee et al. (2011), Stefanon et al. (2015), Caputi et al. (2015), Grazian et al. (2015), Song et al. (2016), and Davidzon et al. (2017). This is quite remarkable, because these SMFs have been recovered from different selection techniques. Specifically, Stark et al. (2009) and Lee et al. (2011) have measurements based on dropouts samples from HST/WFC3 data; the SMF of Grazian et al. (2015) was built from a H₁₆₀-detected photometric redshift sample over the CANDELS/GOODS-S and CANDELS/UDS fields; Stefanon et al. (2015) assembled a composite sample that complements a K_s-detected catalog from UltraVISTA data with detections in IRAC 3.6 and 4.5 μm bands; Caputi et al. (2015) use measurements based on a K_s-detected SMF complemented by SMF measurements from detections in IRAC 4.5 μm. The sample selection of both Stefanon et al. (2015) and Caputi et al. (2015) relies on photometric-redshift measurements. Song et al. (2016) recovered the SMF measurements by converting the UV LF into SMF through a linear ${M}_{* }\mbox{--}{M}_{\mathrm{UV}}$ relation for ${M}_{* }\lt {10}^{10}{M}_{\odot }$, complemented by bootstrapped estimates at ${M}_{* }\gt {10}^{10}{M}_{\odot }$. Finally, the SMF of Davidzon et al. (2017) was based on a photometric-redshift sample from K_s-detection in COSMOS/UltraVISTA DR2 mosaics. On the other side, the normalization of our SMF is higher than González et al. (2011) SMFs; these measurements were obtained by converting the observed UV LF into SMF through ${M}_{* }/{L}_{\mathrm{UV}}$ measurements. The discrepancy between our SMF (and the bulk of the other SMF determinations) could be due to a steeper ${M}_{* }/{L}_{\mathrm{UV}}$ relation found by González et al. (2011), and consequent lower normalization term.

At the massive end ($\mathrm{log}({M}_{* }/{M}_{\odot })\gtrsim 11.2$), we observe a discrepancy between our z ∼ 4 SMF and some of the corresponding measurements from the literature (e.g., Ilbert et al. 2013; Muzzin et al. 2013b). This discrepancy could, at least in part, be explained by our SMF lacking any scatter in ${M}_{* }/{L}_{z^{\prime} }$ for a given ${L}_{z^{\prime} }$. Our stacking analysis, by construction, recovers median ${M}_{* }/{L}_{z^{\prime} }$ ratios, potentially excluding such extreme cases as very dusty/old systems with very high stellar masses. However, the bottom panel of Figure 6, and our discussion in Section 2.3, show that our sample is not strongly biased against this class of objects in the limits of current data. Nonetheless, they do not allow us to properly ascertain their existence at lower stellar masses. Aside from that, because of the Eddington bias, a distribution in the observed ${M}_{* }/{L}_{z^{\prime} }$ values for a specific luminosity would introduce a higher fraction of lower stellar mass objects that are scattered to higher stellar masses than the opposite, increasing the number density of the massive objects. One additional potential reason for this discrepancy is that the current and previous estimates of photometric redshifts do not agree. Compared to Muzzin et al. (2013b) or Ilbert et al. (2013), the DR2 version of the UltraVISTA catalog benefits from deeper NIR and IRAC data, providing improved photometric redshift constraints. The agreement between our estimates and those of Grazian et al. (2015) and Davidzon et al. (2017), both obtained from photometric redshift samples, is noteworthy—suggesting that our composite sample suffers from small selection bias, even at higher masses.

At z ∼ 5, the measurements of the low-mass end ($\mathrm{log}({M}_{* }/{M}_{\odot })\sim 9\mbox{--}9.5$) of the SMF are characterized by a larger scatter than for $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10\mbox{--}10.5$. At higher stellar masses, measurements are broadly consistent with each other—mainly because of the large uncertainties on the number densities. Our SMF determination overlaps with the measurements of Stark et al. (2009), González et al. (2011), Lee et al. (2011), Grazian et al. (2015), Song et al. (2016), and Davidzon et al. (2017). However, it lies below the SMF measurements of Duncan et al. (2014).

At even higher redshift, the Poisson uncertainties on the number densities start to be of the same order as those from cosmic variance. The small number of galaxies at z ∼ 6 and z ∼ 7 (10 and 2 galaxies, respectively) generates large Poisson uncertainties, resulting in a broad agreement among the different SMF determinations. Our SMF overlaps with the measurements of Stark et al. (2009), González et al. (2011), Grazian et al. (2015), and Song et al. (2016). The large uncertainties also make our SMF measurements roughly consistent with those of Duncan et al. (2014).

The shape of the LF and SMF of galaxies at high redshift is usually well-described by a Schechter (1976) function over a wide range of redshifts ($2\lesssim z\lesssim 10$, but see also, e.g., Bowler et al. 2014). We therefore fitted our LF measurements with a Schechter function:

$$\begin{eqnarray}&&{\rm{\Phi }}(M)={{\rm{\Phi }}}^{* }\displaystyle \frac{\mathrm{ln}(10)}{2.5}{10}^{-0.4(M-{{ \mathcal M }}^{* })(\alpha +1)}\exp (-{10}^{-0.4(M-{{ \mathcal M }}^{* })}),\end{eqnarray} \tag{ 8 }$$

where ${{ \mathcal M }}^{* }$ is the characteristic magnitude, corresponding to the knee of the density distribution. Here, α is the faint end slope and ${{\rm{\Phi }}}^{* }$ is a global normalization factor. The single-bin measurement at z ∼ 7 does not allow us to properly inspect the shape of the LF. However, LF estimates at similar redshift, although in different bands, show that it is reasonably well-determined and consistent with a single Schechter form (Bouwens et al. 2015; Finkelstein et al. 2015a), or a double power law (Bowler et al. 2014, 2017).

We performed a best  fit to our LFs' measurements, using the Levenberg-Marquardt method. For the z ∼ 4 LF, we left the three Schechter parameters free to vary. Because the faint-end slope of the $z\geqslant 5$ LFs is poorly sampled, the best fits to the $z\geqslant 5$ LFs were done by fixing α to the value obtained for the z ∼ 4 LF. The fit to the z ∼ 5 and z ∼ 6 LFs were done in three different configurations: both ${{ \mathcal M }}^{* }$ or ${{\rm{\Phi }}}^{* }$ as free parameters, with ${{ \mathcal M }}^{* }$ as unique free parameter, or with ${{\rm{\Phi }}}^{* }$ as the only free parameter. We did not fit the z ∼ 7 LF, leaving both ${{ \mathcal M }}^{* }$ or ${{\rm{\Phi }}}^{* }$ as free parameters. The resulting best-fit Schechter functions for the case of pure luminosity evolution are marked in Figure 10 with solid curves. Table 3 lists the recovered parameters and their uncertainties.

Table 3.  Luminosity Function Best-fit Schechter Parameters

Redshift	${{\rm{\Phi }}}^{* }$	${{ \mathcal M }}^{* }$	α
bin	(10−5 Mpc−3 mag−1)	(mag)
4	19.4 ± 7.1	−23.38 ± 0.19	−1.79 ± 0.09
5	9.1 ± 7.7	−22.99 ± 0.51	−1.79
	19.4	−22.62 ± 0.12	−1.79
	6.85 ± 1.36	−23.38	−1.79
6	51.1 ± 36.0	−21.77 ± 0.42	−1.79
	19.4	−22.09 ± 0.17	−1.79
	8.59 ± 5.40	−23.38	−1.79
7	19.4	−21.82 ± 0.39	−1.79
	1.99 ± 1.62	−23.38	−1.79

Download table as:  ASCIITypeset image

Visual inspection of the best-fitting Schechter function shows that, overall, there is a preference for a pure luminosity evolution against a pure density evolution. Moreover, when both ${{ \mathcal M }}^{* }$ and ${{\rm{\Phi }}}^{* }$ were left free to vary, the values of ${{\rm{\Phi }}}^{* }$ were characterized by large uncertainties, making them consistent with the value of ${{\rm{\Phi }}}^{* }$ at z ∼ 4, i.e., corresponding to no evolution with redshift. On the contrary, the values of ${{ \mathcal M }}^{* }$ showed a more clear trend with redshift, increasing our confidence in the luminosity evolution. We note, however, that these results potentially suffer from the limited coverage of ${{ \mathcal M }}^{* }$.

The luminosity evolution registered here is in contrast with the most recent measurements of the evolution of the UV LF at $z\gt 4$, where indication is found that the characteristic magnitude evolves very little with redshift (likely constrained by the impact of dust extinction at high masses—Bouwens et al. 2009; Reddy et al. 2010), as most of the evolution seems to be driven by variation in the overall density (see, e.g., van der Burg et al. 2010; Bouwens et al. 2015; Finkelstein et al. 2015a).

The LFs presented in our work are the first determination of the rest-frame optical LF at $z\gtrsim 4$, and therefore direct comparisons to previous estimates are not possible. However, recent works have recovered the measurement of the characteristic magnitude in the rest-frame z'-band at z ∼ 4 (Oesch et al. 2013), or have studied the evolution of the LF up to $z\lesssim 4$ in rest-frame optical bands close to z' (Marchesini et al. 2012; Stefanon & Marchesini 2013). In the following paragraphs, we will compare our determination of the z ∼ 4 characteristic magnitude to the estimates from the above three works.

Oesch et al. (2013) estimated the characteristic magnitude of the z ∼ 4 LF by applying a correction based on the average ${i}_{775}-[4.5]$ color to the characteristic magnitude of the z ∼ 4 UV LF. The obtained value, ${{ \mathcal M }}_{z}^{* }=-21.7$ AB, is ∼1.7 mag fainter than what we found in this work. A possible reason for this large discrepancy is the lack of galaxies brighter than ${M}_{z}\lt -23$ from the sample of Oesch et al. (2013), which instead have become accessible through the deep and wide area UltraVISTA data.

Marchesini et al. (2012) presented estimates of the rest-frame V-band LF at $z\lesssim 4$, obtained from a composite sample that includes wide-area data from the NMBS, FIRES, FIREWORKS, HDFN, HUDF, and GOODS-S programs. The characteristic magnitude recovered from a maximum likelihood analysis is ${{ \mathcal M }}_{V}^{* }=-{22.76}_{-0.63}^{+0.40}$. This value is brighter than that found by Oesch et al. (2013), although it is still ∼0.6 mag fainter than the estimate in our work. Nonetheless, it is consistent at $1\sigma$ with our estimate, considering the associated uncertainties. Our stacking analysis (see Figure 8) shows that the brighter galaxies have redder $(V-z^{\prime} )$ colors, reaching $(V-z^{\prime} )\sim 1$ mag for the brightest stacks. Using our stacked SED at ${M}_{z^{\prime} }\sim -23.5$, close to the value of the characteristic magnitude of our z ∼ 4 LF, we find a rest-frame color $(V-z^{\prime} )=+0.53$ mag; applied to the ${{ \mathcal M }}_{V}^{* }$, it gives a ${{ \mathcal M }}_{z^{\prime} }^{* }=-23.29$ mag, very close to our best-fit ${{ \mathcal M }}_{z^{\prime} }^{* }\,=-23.38$ mag.

Stefanon & Marchesini (2013) measured the evolution of the rest-frame J-and H-band LF up to $z\sim 3.5$ by using a composite sample of galaxies from the MUSYC, FIRES, and GOODS-CDFS programs. From a maximum likelihood analysis, the characteristic magnitude of the $z\sim 3.25$ rest-frame J-band LF was estimated to be ${{ \mathcal M }}_{J}^{* }=-{23.28}_{-0.29}^{+0.33}$ mag. Applying the same analysis adopted for the comparison to Marchesini et al. (2012), we find a rest-frame color $(z^{\prime} -J)=+0.11$ mag, and a corresponding ${{ \mathcal M }}_{z^{\prime} }^{* }=-23.17$ mag, still consistent with our estimate at $1\sigma$ level.

We also performed a Schechter fit to the SMFs. However, the differential ${M}_{* }/{L}_{z^{\prime} }$ we measured at z ∼ 4 and applied to the $z\geqslant 5$ LFs has the effect of stretching the original LFs. The consequence of this stretching, along with the limited number of measurements available at each redshift, is that the $z\geqslant 5$ SMFs do not show any robust evidence for the exponential cut at the massive end, preventing any reliable estimate of the characteristic stellar mass. A similar result was found by Grazian et al. (2015), even though the issue is not yet settled (see, e.g., Caputi et al. 2015). Therefore, we performed the Schechter fit only on the z ∼ 4 SMF, and obtained the following results: $\alpha \,=-1.93\pm 0.24$, ${{\rm{\Phi }}}_{* }^{* }=(2.72\pm 1.81)\times {10}^{-5}$ Mpc⁻³ dex⁻¹ and $\mathrm{log}({M}_{* }^{* }/{M}_{\odot })=10.96\pm 0.33$. The corresponding best-fit Schechter function is presented in Figure 11.

The large uncertainties associated with the Schechter parameters make them consistent with most of the measurements from the literature. The low number of massive galaxies from the exponential part of the Schechter function, at the massive end, suffer from high uncertainties from both (relative) Poisson noise and cosmic variance. The massive end of our SMF is based on the UltraVISTA sample, suggesting a more adequate coverage of massive galaxies.

One of the main aims of this paper is to test the LF in rest-frame optical bands, as a proxy for SMF measurements at high redshifts ($z\gtrsim 4$). The compilation of SMF measurements presented in Figure 11 shows that, for $z\geqslant 4$, they are characterized by a scatter that can be as large as $\gtrsim 1$ dex. This has immediate consequences for our understanding of more fundamental and global properties, like the evolution with cosmic time of the stellar mass density or the stellar-to-halo mass relation. The likely main reason for this large scatter is sample selection. In addition to this, further systematics may arise from our limited knowledge of some of the specific aspects characterizing the stellar population of each galaxy, including contamination from emission lines (see, e.g., Conroy & Wechsler 2009; Stark et al. 2009; Behroozi et al. 2010).

Reliable stellar mass measurements require coverage of rest-frame optical/NIR wavelengths, probed by Spitzer/IRAC at $z\gtrsim 4$. In particular, current measurements of the low-mass end of the $z\gt 4$ SMFs rely on IRAC data over the GOODS-N/S fields, as these are the only fields with sufficient depth that provide photometric coverage over a sufficient area. The 3.6 and 4.5 μm bands, although characterized by deeper data than the 5.8 and 8.0 μm bands at $z\gtrsim 4$, are potentially contaminated by nebular lines. Previous work has attempted to estimate the impact that nebular emission can have on the recovery of stellar mass, yet no concordance has been found so far. Specifically, estimates based on relations between the EW of Hβ and the ionizing properties of the best-fit SEDs generally predict systematics $\lesssim 0.2$ dex, and are nearly independent of redshift (e.g., Duncan et al. 2014; Grazian et al. 2015; Salmon et al. 2015). On the other side, estimates based on the measurements of the evolution of the Hα EW with redshift indicate an increasing contribution of nebular lines to the stellar mass measurements, ranging from ∼0 dex at z ∼ 4, to z ∼ 0.6 dex at z ∼ 7, as seen in Stark et al. 2013. Those authors presented, among other estimates, an extrapolation to $z\sim 7$, later supported by other observations (e.g., Smit et al. 2014; Faisst et al. 2016; Rasappu et al. 2016). Such systematics in the measurement of stellar mass introduce up to ∼0.5 dex offset in the number densities for the higher redshift bins. To further complicate the picture, recent works have shown that, when the correction for nebular line contamination is applied on a statistical basis, irrespective of the specific SEDs, it can even boost the stellar mass measurements by up to ∼0.2–0.3 dex (Stefanon et al. 2015; Nayyeri et al. 2017; Stefanon et al. 2017).

One possible way to circumvent this problem is estimating the stellar mass exclusively from the 5.8 and 8.0 μm bands. These bands cover a region of the SEDs of $z\gtrsim 4$ galaxies free from contamination by strong nebular lines. However, the larger PSF FWHM of Spitzer/IRAC 5.8 and 8.0 μm bands, compared to the 3.6 and 4.5 μm bands, may introduce blending effects in the flux measurements. Labbé et al. (2015) have shown that mophongo, the software we adopted to measure the fluxes in the IRAC bands, does not introduce any substantial systematics in the flux measurement, even in very crowded regions. Nonetheless, given the larger FWHM and the lower S/N characterizing the 5.8 μm and 8.0 μm data, one could expect an increased scatter in the flux measurements, compared to the bluer IRAC bands.

Through the ${M}_{* }/{L}_{z^{\prime} }$ we derived in the present work, the S/N cuts we applied to the flux in the IRAC 5.8 and 8.0 μm bands identify a range in stellar mass where their measurements can be considered reliable. The LF and SMF from this work can thus be regarded as a indicative of most of the current SMF measurements at $z\geqslant 4$, as they are based on the subsample of objects with the highest S/N measurements in those bands more sensitive to the stellar mass and with reduced contamination from nebular emission. This is visible in Figure 11; at z ∼ 5 and above the lowest stellar mass over which our SMFs are defined, we see S/N $\approx \,5\mbox{--}10\times$ higher than most current SMF determinations. Specifically, this also means that stellar masses below our low-mass limits are necessarily based on very low S/N measurements in the IRAC 5.8 and 8.0 μm bands, (still uncertain) correction for nebular emission contamination, or a combination of the two.

Our analysis showed that SMFs consistent with the average determinations from the literature could be recovered by applying a simple ${M}_{* }/{L}_{z^{\prime} }$ relation to the observed z'-band LF, with stellar masses measured from common stellar population parameters (e.g., delayed exponential SFH, solar metallicity, Chabrier IMF). Our simple transformation of the LF into SMF was supported by the non-decreasing ${M}_{* }/{L}_{z^{\prime} }$ ratio for increasing luminosities observed at z ∼ 4. Relations between the ${M}_{* }/L$ and L other than that (e.g., if the ${M}_{* }/L$ presented a minimum for some value of L) would still allow the conversion, but would require us to consider, for specific bins of M_*, the contributions to the density originating from different bins of L. A non-decreasing ${M}_{* }/{L}_{z^{\prime} }$, instead, constitutes an injective mapping between luminosity and stellar mass; galaxies with higher luminosity will always have higher stellar mass, and galaxies with lower luminosity will always have lower stellar mass.

A different scenario arises when multiple ${M}_{* }/L$ exist in correspondence with a single value of L, as can be the case, for instance, for UV luminosity versus stellar mass. Figure 6 (along with the bottom panel of Figure 9) shows that, for ${M}_{\mathrm{UV}}\gtrsim -20$ mag, there are broadly two very different ${M}_{* }/{L}_{\mathrm{UV}}$ values. More specifically, this means that, if galaxies with a higher ${M}_{* }/{L}_{\mathrm{UV}}$ ratio are included in the sample, those bins of luminosity will include (and mix) the contribution from both high- and low-mass galaxies. If this effect is not properly taken into account, it introduces an overestimate of the low-mass-end slope and an underestimate of the massive end of the SMF. The only way to deal with this problem is to directly count the number of galaxies in each of the two ${M}_{* }/L$ bins.

Finally, the ideal rest-frame band for this kind of study is probably one for which the ${M}_{* }/L$ ratio would not depend on the luminosity, as any luminosity dependence could potentially hide effects from, e.g., SFH. Our measurements of the ${M}_{* }/{L}_{z^{\prime} }$ relation suggest that, for a wide range in luminosity, they are consistent with a constant value. A log-linear relation arises at the bright (massive) end of the LF (SMF), suggesting that the z'-band at these luminosities holds the signature of the stellar population age and/or the dust content, at some level.

Given the potential systematics on the SMF measurements discussed above, LF estimates can provide a valid alternative for recovering the halo masses (M_h) for high-redshift galaxies through abundance matching techniques (e.g., Behroozi et al. 2013b; Finkelstein et al. 2015b; Steinhardt et al. 2016). To date, measurements of $z\geqslant 4$ LFs are mostly available in the rest-frame UV. The adoption of UV LFs in the M_h/L estimates provides information on the relative importance of star formation processes (e.g., gas cooling, stellar ejections, SFR timescales) versus the hierarchical growth of the dark matter halos (Bouwens et al. 2015). However, stellar masses are likely to be more strongly correlated with the halo masses than UV luminosities.

The Schechter fits performed in Section 3.5 suggest that the evolution of the z'-band LFs can be accounted for by an increase of luminosity with cosmic time uniformly across luminosities at a given redshift. Because we have shown that the z'-band LF is a reasonable proxy for the SMF, it is tempting to analyze the evolution of the rest-frame optical LF obtained in the present work in terms of the evolution of the dark-matter HMF.

We perform a first analysis as follows. We apply a simple abundance-matching technique (Vale & Ostriker 2004) consisting of matching the cumulative number density of the LF to that of Behroozi et al. (2013b) HMF obtained from HMFcalc ¹¹ (Murray et al. 2013) and recover the evolution in mass of the HMF at $z\geqslant 5$ relative to z ∼ 4. Given the rapid evolution of the HMF in this range of redshift, we adopt the HMFs at $z=3.78,4.95,5.76$, and 6.87, corresponding to the (median) photometric redshifts of the stacked SEDs. The HMFs assume ${\sigma }_{8}=0.81$. We find a relative displacement in halo mass of 0.76, 0.52, and 0.35 dex, corresponding to ∼1.9, ∼1.3, and ∼0.9 mag from $z\sim 7,6$ and 5 to z ∼ 4, respectively. Figure 12 shows the result of applying the above offsets to the ${{ \mathcal M }}_{z^{\prime} }^{* }$ of the z ∼ 4 Schechter parameterization. The solid blue curve marks the best-fit Schechter function at z ∼ 4, whereas the dashed curves represent the z ∼ 4 LF rigidly shifted by the corresponding amount at $z\geqslant 5$. The agreement between the predicted and observed LF is very good at all redshifts, suggesting that the z'-band LF could trace the evolution of the HMF.

Zoom In Zoom Out Reset image size

In this section, we discuss the evolution of those halos associated with a constant cumulative number density of $3.1\times {10}^{-5}$ Mpc⁻³ over $4\lesssim z\lesssim 7$. This value corresponds to the cumulative number density of ${{ \mathcal M }}_{z^{\prime} }^{* }$ galaxies at z ∼ 4 (${M}_{z^{\prime} }=-23.38$ mag), and allows us to recover the corresponding absolute magnitude up to z ∼ 7 with strongly reduced dependence (i.e., $\lesssim 0.5$ mag) on the extrapolation of the LFs to magnitudes fainter than actually observed. Table 4 lists a compilation of the values of the main parameters recovered with our analysis.

Table 4.  Values of the Main Observables from our Abundance Matching Analysis

	z	Cum. Den.a	${M}_{z^{\prime} }$	Mh	${L}_{z^{\prime} }/{M}_{h}$	M*	${M}_{* }/{M}_{h}$	ISFEb
	bin	$[\mathrm{log}({\mathrm{Mpc}}^{-3})]$	[mag]	$[\mathrm{log}({M}_{h}/{M}_{\odot })]$	[L⊙/M⊙]	$[\mathrm{log}({M}_{* }/{M}_{\odot })]$
Fixed den.	3.78	−4.51	−23.39 ± 0.27	12.32	${0.069}_{-0.015}^{+0.019}$	10.46 ± 0.16	${0.014}_{-0.004}^{+0.006}$	${0.087}_{-0.027}^{+0.038}$
	4.95	−4.51	−22.63 ± 0.23	11.91	${0.089}_{-0.017}^{+0.021}$	10.02 ± 0.13	${0.013}_{-0.003}^{+0.005}$	${0.081}_{-0.021}^{+0.029}$
	5.76	−4.51	−22.09 ± 0.26	11.67	${0.095}_{-0.020}^{+0.025}$	9.77 ± 0.10	${0.013}_{-0.003}^{+0.003}$	${0.082}_{-0.017}^{+0.022}$
	6.87	−4.51	−21.82 ± 0.43	11.40	${0.138}_{-0.045}^{+0.067}$	9.67 ± 0.17	${0.019}_{-0.006}^{+0.009}$	${0.118}_{-0.039}^{+0.058}$
Evol. den.	3.78	−4.51	−23.39 ± 0.27	12.32	${0.069}_{-0.015}^{+0.019}$	10.46 ± 0.16	${0.014}_{-0.004}^{+0.006}$	${0.087}_{-0.027}^{+0.038}$
	4.95	−4.32	−22.41 ± 0.23	11.83	${0.088}_{-0.017}^{+0.020}$	9.90 ± 0.12	${0.012}_{-0.002}^{+0.004}$	${0.075}_{-0.015}^{+0.024}$
	5.76	−4.14	−21.65 ± 0.26	11.51	${0.090}_{-0.019}^{+0.024}$	9.60 ± 0.10	${0.012}_{-0.003}^{+0.003}$	${0.078}_{-0.016}^{+0.021}$
	6.87	−3.88	−21.02 ± 0.43	11.14	${0.119}_{-0.039}^{+0.058}$	9.34 ± 0.17	${0.016}_{-0.005}^{+0.008}$	${0.102}_{-0.034}^{+0.050}$

Notes.

^aCumulative number density adopted for the abundance matching. ^bIntegrated Star Formation Efficiency $\equiv \,({M}_{* }/{M}_{h})/({{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m})$, ${{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}\equiv 0.157$.

Download table as:  ASCIITypeset image

The abundance matching performed through cumulative number density implicitly assumes that each halo contains one and only one galaxy, and that halos of the same mass contain galaxies of the same z' luminosity (∼stellar mass). Indeed, recent measurements have shown that the scatter between halo mass-stellar mass relation is quite small, ∼0.15–0.20 dex (e.g., Conroy & Wechsler 2009; Moster et al. 2010; Tinker et al. 2017; but see also Gu et al. 2016, who found scatter of up to 0.32 dex).

The cumulative number densities of the LFs were computed adopting the Schechter parameterization presented in Section 3.5. For the LFs at $z\geqslant 5$, we adopted the best-fit Schechter functions obtained when the characteristic magnitude was assumed to be the only free parameter of the fit, coinciding with the case of pure luminosity evolution.

The values of the absolute magnitudes we obtain from our procedure are ${M}_{z^{\prime} }\sim -23.39,-22.63,-22.09,-21.82$ mag for the z ∼ 4, z ∼ 5, z ∼ 6, and z ∼ 7 cases, respectively. These values are consistent within $1\sigma$, with the characteristic magnitudes of our Schechter fit. This is not surprising, considering the pure luminosity evolution of the Schechter fits themselves.

The matches to the cumulative number density performed on the HMFs resulted in halo masses $\mathrm{log}({M}_{h}/{M}_{\odot })\,=12.32,11.91,11.67,11.40$ for the z ∼ 4, z ∼ 5, z ∼ 6, and z ∼ 7 cases, respectively. These displacements correspond to an evolution in the halo mass of ${{ \mathcal M }}_{z^{\prime} }^{* }$ galaxies of $\sim 0.9,\sim 0.6$, and ∼0.3 dex from $z\sim 7,6,5$ to z ∼ 4, respectively. We note, however, that a rigid displacement in mass of the HMF is sufficient to reproduce the HMF evolution at these redshifts only for halo masses $\mathrm{log}({M}_{h}/{M}_{\odot })\gtrsim 12;$ at lower halo masses, the displacement in mass must be coupled to a steepening with redshift of the low-mass end slope.

We can now use the above results on the evolution of the luminosity and of the halo mass to recover the evolution with redshift of the light-to-halo mass for galaxies at fixed cumulative number density. Combining the two, we obtain ${L}_{z^{\prime} }/{M}_{h}\sim 0.069,0.089,0.095,0.138$ in units of L_⊙/M_⊙. These values are also presented in the left panel of Figure 13 and suggest a mild increase with redshift (a factor $\lesssim 2\times$), although the large uncertainties make them consistent with a constant value across the 800 Myr of cosmic time, from z ∼ 4 to z ∼ 7.

Zoom In Zoom Out Reset image size

Using the results on the ${M}_{* }/{L}_{z^{\prime} }$ from Section 3.2, we can convert the ${L}_{z^{\prime} }/{M}_{h}$ into ${M}_{* }/{M}_{h}$. The result of this is shown in the right panel of Figure 13. The corresponding values are listed in Table 4. The stellar-to-halo mass ratio does not present any significant evolution with redshift. In the same panel, we convert the ${M}_{* }/{M}_{h}$ into the integrated star formation efficiency (ISFE), i.e., ${M}_{* }/{M}_{h}$ in units of ${{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}$ (Conroy & Wechsler 2009), which is equivalent to the stellar baryon fraction (Finkelstein et al. 2015b), using ${{\rm{\Omega }}}_{b}/{{\rm{\Omega }}}_{m}=0.157$ (Planck Collaboration et al. 2016). Our measurements are consistent with the ISFE being constant with redshift. We stress here that this result refers to ${M}_{h}\gtrsim {10}^{12}{M}_{\odot }$ and does not exclude the existence of evolution with redshift at lower halo masses. We defer a more complete analysis of the dependence of the ${M}_{* }/{M}_{h}$ ratio on halo mass to a future work. Furthermore, our samples at $z\sim 5\mbox{--}7$ are entirely based on LBG selection. If non-negligible numbers of redder (dustier/more evolved) galaxies exist at these epochs, they would affect SMF (and likely its massive end, e.g., Caputi et al. 2015; Stefanon et al. 2015) and, consequently, the recovered M_h.

Matching galaxies at a constant cumulative number density, however, does not consider the effect of major mergers in the galaxy ranking. We therefore repeated the same analysis using a cumulative number density that evolves with redshift, following the recipe of Behroozi et al. (2013a). The results are listed in Table 4, and plotted as red circles in Figure 13. No significant difference with the constant cumulative number density match is observed. We also note that the values for the z ∼ 7 bin rely on the extrapolation of the LF to luminosities below those currently probed by our sample. Those measurements should then be treated with caution.

In the right panel of Figure 13, we also plot recent estimates of the ${M}_{* }/{M}_{h}$ from the literature: Durkalec et al. (2015), Finkelstein et al. (2015a, converted to a Chabrier 2003 IMF by applying a factor 0.55) and Harikane et al. (2016). Durkalec et al. (2015) applied the measurements of the two-point correlation function to a halo occupation model to recover the halo mass of samples of galaxies at $z\sim 2\mbox{--}5$ with spectroscopic redshift from the VIMOS Ultra Deep Survey (Le Fèvre et al. 2015). Finkelstein et al. (2015b) measured the evolution of the ${M}_{* }/{M}_{h}$ by abundance matching the $z\sim 4\mbox{--}7$ UV LF. Harikane et al. (2016) recovered ${M}_{* }/{M}_{h}$ from the clustering of LBGs selected at $z\sim 4\mbox{--}7$ from a variety of programs, including CANDELS, the Hubble Frontier Fields (PI: J. Lotz), and Subaru Hyper-Suprime Cam Subaru Strategic Program (PI: S. Miyazaki).

At face value, our measurements are consistent with those of Durkalec et al. (2015) at z ∼ 4, and with those of Finkelstein et al. (2015b) at $z\sim 6\mbox{--}7;$ however, they are inconsistent with those of Harikane et al. (2016) over the full range of redshift, and with those of Finkelstein et al. (2015b) at $z\sim 4\mbox{--}5$. Recently, Mancuso et al. (2016), applying abundance matching to the evolution of the SFR function recovered from UV+far-IR data, found indication for a non-evolving ${M}_{* }/{M}_{h}$ ratio at $z\geqslant 4$.

The ${M}_{* }/{M}_{h}$ measurements presented in the right panel of Figure 13 were obtained from a variety of methods, and ultimately refer to different M_h estimates, making a straightforward comparison difficult to interpret. Specifically, our measurements of the halo mass are based on a constant cumulative number-density match with $\mathrm{log}({M}_{h}/{M}_{\odot })\sim 12.3,11.9,11.7,11.4$ at $z\sim 4,5,6,7$, respectively. The increase with redshift of the ${M}_{* }/{M}_{h}$ of Harikane et al. (2016) refers to a fixed $\mathrm{log}({M}_{h}/{M}_{\odot })=11$ across $z\sim 4\mbox{--}7$. Finkelstein et al. (2015b) report an increase with redshift of the ISFE $\propto \,(0.024\pm 0.07)\,\times z$. However, this trend is most likely driven by the point at z ∼ 4, and originates from limited evolution of the UV LF between z ∼ 4 and z ∼ 5 (${\rm{\Delta }}{{ \mathcal M }}_{\mathrm{UV}}^{* }\sim 0.08\pm 0.16$ mag) observed by Finkelstein et al. (2015a), in contrast to the large-luminosity evolution (${\rm{\Delta }}{{ \mathcal M }}_{z^{\prime} }^{* }\sim 0.76\pm 0.22$ mag) observed in our work over the same redshift interval. This, along with the constant characteristic magnitude of the UV LF (i.e., ≈constant SFR) assumed as a criterion for the abundance matching, generates a reference cumulative number density that decreases with redshift, and halo masses $\mathrm{log}({M}_{h}/{M}_{\odot })\sim 11.9,11.7,11.6,11.3$. We note here that our measurements are consistent with those of Finkelstein et al. (2015b) at $z\sim 6\mbox{--}7$, i.e., where the M_h recovered by the two teams are more similar, whereas they are inconsistent at $z\sim 4\mbox{--}5$, where the M_h differ. Finally Durkalec et al. (2015) estimates refer to a diversity of halo masses and redshift ranges: $\mathrm{log}({M}_{h}/{M}_{\odot })\sim 11.1,11.5,11.2$ for ${z}_{\mathrm{mean}}\sim 2.5,3.0,3.5$, respectively.

Our finding of ${M}_{* }/{M}_{h}$ independent of redshift is also qualitatively in agreement with recent estimates of the galaxy bias, observed to be nearly constant over $z\sim 4\mbox{--}6$ (Barone-Nugent et al. 2014), although a measurement at z ∼ 7 from the same work seems to suggest a potential change of the star formation efficiency at earlier epochs.

The picture is not settled, even from a theoretical perspective. Indeed, some of the models predict an increase with cosmic time of the ${M}_{* }/{M}_{h}$ ratio for a fixed M_h (e.g., Somerville et al. 2015). Other models find that ${M}_{* }/{M}_{h}$ decreases with cosmic time at fixed halo mass (Moster et al. 2013; Behroozi & Silk 2015). However, when considering the evolution of the same population of galaxies (through an evolving number density), Behroozi & Silk (2015) find that ${M}_{* }/{M}_{h}$ is nearly independent on redshift. Finally, other models, instead, find ${M}_{* }/{M}_{h}$ to be insensitive to redshift (e.g., O'Shea et al. 2015; Mutch et al. 2016). Some semi-empirical models have been able to reproduce the evolution of the UV LF from z ∼ 2 to z ∼ 10 under the assumption that, for star-forming galaxies, the ${M}_{* }/{M}_{h}$ depended on M_h but not on redshift (Trenti et al. 2010; Tacchella et al. 2013; Mason et al. 2015).