# Spitzer bright, ultravista faint sources in cosmos: the contribution to the overall population of massive galaxies at z=3-7

K. I. Caputi, O. Ilbert, C. Laigle, H. J. Mccracken, O. Le Fevre, J. Fynbo, B.

Milvang-Jensen, P. Capak, M. Salvato, Y. Taniguchi

► To cite this version:

K. I. Caputi, O. Ilbert, C. Laigle, H. J. Mccracken, O. Le Fevre, et al. Spitzer bright, ultravista faint sources in cosmos: the contribution to the overall population of massive galaxies at z=3-7. The Astrophysical Journal, 2015, 810 (1), pp.73. 10.1088/0004-637X/810/1/73. hal-01434385

HAL Id: hal-01434385

https://hal.science/hal-01434385

Submitted on 15 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SPITZER BRIGHT, ULTRAVISTA FAINT SOURCES IN COSMOS: THE CONTRIBUTION TO
THE OVERALL POPULATION OF MASSIVE GALAXIES AT z = 3–7

K. I. CAPUTI¹, O. ILBERT², C. LAIGLE³, H. J. MCCRACKEN³, O. LE FÈVRE², J. FYNBO⁴, B. MILVANG-JENSEN⁴, P. CAPAK⁵,
M. SALVATO⁶,⁷, AND Y. TANIGUCHI⁸

¹ Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700 AV Groningen, The Netherlands; karina@astro.rug.nl

² Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille), UMR 7326, F-13388, Marseille, France

³ Sorbonne Universités, UPMC Université Paris 6 et CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France

⁴ Dark Cosmology Center, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark

⁵ Spitzer Science Center, California Institute of Technology, Pasadena, CA 91125, USA

⁶ Max-Planck-Institut für Extraterrestrische Physik, Postfach 1312, D-85741 Garching bei München, Germany

⁷ Excellence Cluster, Boltzmann Strasse 2, D-85748 Garching, Germany

⁸ Research Center for Space and Cosmic Evolution, Ehime University, 2-5 Bunkyo-cho, Matsuyama 790-8577, Japan

ABSTRACT

We have analyzed a sample of 574 Spitzer 4.5 µm selected galaxies with [4.5] < 23 and $K_s^{auto} > 24$ (AB) over the UltraVISTA ultradeep COSMOS field. Our aim is to investigate whether these mid-infrared (mid-IR) bright, near-infrared (near-IR) faint sources contribute significantly to the overall population of massive galaxies at redshifts $z \ge 3$. By performing a spectral energy distribution (SED) analysis using up to 30 photometric bands, we have determined that the redshift distribution of our sample peaks at redshifts $z \approx 2.5-3.0$, and ~32% of the galaxies lie at $z \ge 3$. We have studied the contribution of these sources to the galaxy stellar mass function (GSMF) at high redshifts. We found that the [4.5] < 23, $K_s^{auto} > 24$ galaxies produce a negligible change to the GSMF previously determined for $K_s^{auto} < 24$ sources at $3 \le z < 4$, but their contribution is more important at $4 \le z < 5$, accounting for >50% of the galaxies with stellar masses $M_{st} \ge 6 \times 10^{10} M_\odot$. We also constrained the GSMF at the highest-mass end ($M_{st} \ge 2 \times 10^{11} M_\odot$) at $z \ge 5$. From their presence at $5 \le z < 6$ and virtual absence at higher redshifts, we can pinpoint quite precisely the moment of appearance of the first most massive galaxies as taking place in the ~0.2 Gyr of elapsed time between $z \sim 6$ and $z \sim 5$. Alternatively, if very massive galaxies existed earlier in cosmic time, they should have been significantly dust-obscured to lie beyond the detection limits of current, large-area, deep near-IR surveys.

Key words: galaxies: evolution – galaxies: high-redshift – infrared: galaxies

1. INTRODUCTION

Constraining the number density of massive ($M_{st} \gtrsim 5 \times 10^{10} M_\odot$) galaxies at different redshifts is very important to understand when galaxy buildup proceeded most efficiently in cosmic time. Over the past decade, multiple studies have shown that a significant fraction of the massive galaxies that we know today were already in place and massive at $z \sim 2$ (e.g., Saracco et al. 2004, 2005; Caputi et al. 2005, 2006; Labbé et al. 2005; Pozzetti et al. 2007). The search for massive galaxies has also been extended to higher redshifts $z \sim 3-4$, usually through the analysis of sources with red near-infrared (near-IR) colors that were mostly undetected at optical wavelengths (e.g., Franx et al. 2003; Kodama et al. 2007; Rodighiero et al. 2007). Until recently, however, the typical depths of deep near-IR surveys ($K_s \sim 23.5-24.0$ AB mag) and small areas covered have prevented a systematic search for the rare massive galaxies present at $z \gtrsim 4$.

The first epoch of appearance of massive galaxies in the early universe constitutes an important constraint for galaxy formation models. These models predict that massive galaxies are formed in the high-density fluctuation peaks of the matter density field (Cole & Kaiser 1989; Mo & White 1996). Reproducing the number density of massive galaxies since redshift $z \sim 7$ is critical to explain how massive galaxy formation proceeded after the epoch of reionization and the subsequent galaxy buildup until today.

The galaxy stellar mass function (GSMF) is a very important tool to investigate the number density evolution of galaxies of different stellar masses through cosmic times. Optical galaxy surveys have enabled the study of the GSMF up to redshift $z \sim 1$ (e.g., Pozzetti et al. 2010; Baldry et al. 2012; Davidzon et al. 2013; Moustakas et al. 2013), while near-IR surveys have extended these studies up to $z \sim 3-4$ (Fontana et al. 2006; Kajisawa et al. 2009; Marchesini et al. 2009; Bielby et al. 2012; Ilbert et al. 2013; Muzzin et al. 2013). Recent ultradeep near-IR surveys over small areas of the sky have allowed for GSMF studies at higher redshifts, providing constraints at the intermediate stellar mass regime, i.e., $M_{st} \sim 10^9-10^{10} M_\odot$ (e.g., González et al. 2011; Santini et al. 2012; Duncan et al. 2014; Grazian et al. 2015). However, the small surveyed areas make it difficult to properly sample the GSMF high-mass end, as massive galaxies are rare at high $z$. Alternatively, galaxy selections based on mid-IR images taken with the *Spitzer Space Telescope* Infrared Array Camera (IRAC; Fazio et al. 2004), some of which were conducted over larger areas of the sky, have offered the possibility of exploring the GSMF high-mass end up to $z \sim 5$ (Pérez-González et al. 2008; Mancini et al. 2009; Ilbert et al. 2010; Caputi et al. 2011; Domínguez-Sánchez et al. 2011; Stefanon et al. 2015).

The depth of the near-IR maps typically available in these larger fields, namely, $K_s \sim 23.5-24.0$, was sufficient to identify the vast majority (~95%-98%) of the IRAC sources to [4.5] = 23 AB mag. The small percentage of sources that
remained unidentified were usually neglected, as it was impossible to derive any of their properties—including their redshifts—without any detection beyond the IRAC bands. Indeed, if these unidentified sources had a similar redshift distribution to those that are identified, then taking them into account would not make any significant difference to the already-derived results, in particular those on the number density of massive galaxies at high z. However, if these missing galaxies had a *biased redshift distribution* toward high redshifts, their contribution in the early universe could be more important than previously assumed.

Hints on the importance of mid-IR bright, near-IR faint sources toward the study of massive galaxies at high redshifts have been provided by the analysis of IRAC extremely red sources (Wiklind et al. 2008; Huang et al. 2011; Caputi et al. 2012), which showed that most of them are massive galaxies at $3 \le z \le 5$. However, these examples are extreme cases of the mid-IR bright, near-IR faint galaxies that are present in ultradeep near-IR maps. These sources are not fully representative of the entire population of bright IRAC sources that remain beyond the typical identification limits of wide-area, near-IR surveys (i.e., $K_s \sim 24.0$).

In this paper we investigate a more representative sample of the IRAC bright sources that have been unidentified so far in large-area, deep near-IR surveys. Our aim is understanding their importance within the overall population of massive galaxies at high redshifts ($z \ge 3$). This study is possible thanks to the ongoing UltraVISTA near-IR survey (McCracken et al. 2012), which at the current stage already has a unique combination of area and photometric depth. The layout of this paper is as follows. In Section 2 we describe the photometric data sets used. In Section 3, we explain our sample selection and photometric redshift determinations through a spectral energy distribution (SED) analysis. In Sections 4 and 5, we present our results for galaxies at $3 \le z < 5$ and $5 \le z < 7$, respectively, including their contributions to the GSMF at these redshifts. We conduct an updated analysis of the cosmic stellar mass density evolution in Section 6. Finally, in Section 7, we summarize our findings and give some concluding remarks. We adopt throughout a cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1} \Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. All magnitudes and fluxes are total and refer to the AB system (Oke & Gunn 1983), unless otherwise stated. Stellar masses correspond to a Salpeter initial mass function (IMF) over 0.1–100 $M_\odot$.

## 2. DATA SETS

The Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) comprises a wealth of multi-wavelength imaging and spectroscopic data covering $1.4 \times 1.4 \text{ deg}^2$ of the sky, on a field centered at R.A. = 10:00:28.6 and decl. = +02:12:21.0 (J2000). The field has been defined by the original coverage in the optical $i_{814w}$ band with the Advanced Camera for Surveys on the Hubble Space Telescope. More extensive optical imaging for COSMOS has been taken with broadband, narrowband, and intermediate-band filters for SuprimeCam on the Subaru Telescope (Taniguchi et al. 2007). In the $u^*$ band, COSMOS has been observed with Megacam on the Canada-France-Hawaii Telescope (CFHT).

In addition, COSMOS has been targeted at mid- and far-IR wavelengths by the Spitzer Space Telescope (Werner et al. 2004) during the cryogenic mission, as part of the Spitzer Legacy Program S-COSMOS (Sanders et al. 2007). The four-

band observations with the IRAC (Fazio et al. 2004) resulted in source catalogs that are ~80% and ~70% complete at mag = 23 for 3.6 µm and 4.5 µm, respectively. The observations with the Multiband Photometer for Spitzer (MIPS; Rieke et al. 2004) achieved a depth of ~70 µJy at 24 µm (~80% catalog completeness).

The ongoing UltraVISTA survey (McCracken et al. 2012) is providing ultradeep near-IR images of the COSMOS field in four broad bands ($Y, J, H$, and $K_s$) and a narrow band NB118 (Milvang-Jensen et al. 2013). The survey strategy is such that it produces alternate deep and ultradeep stripes oriented N–S, covering a total of ~1.5 × 1.23 deg$^2$. The data used here correspond to the data release version 2 (DR2), which achieves depths of $K_s \approx 24.8$, $H \approx 24.7$, $J \approx 25.1$, and $Y \approx 25.4$ (5$\sigma$; 2 arcsec diameter apertures) on the ultradeep stripes (C. Laigle et al. 2015, in preparation). These values are 0.7–1.1 mag fainter than the characteristic depths of the UltraVISTA DR1 release discussed by McCracken et al. (2012).

## 3. A SAMPLE OF SPITZER BRIGHT ([4.5] < 23) SOURCES WITH ULTRAVISTA FAINT ($K_s^{\text{auto}} > 24$) COUNTERPARTS IN THE COSMOS FIELD

### 3.1. Sample Selection and Multiwavelength Photometry

We used the publicly available S-COSMOS IRAC catalog<fnref n="9" /> (Ilbert et al. 2010) to select a sample of bright 4.5 µm sources with $[4.5] < 23$. To ensure the most reliable source detection and photometry, we only considered sources detected also at 3.6 µm and with an extraction flag equal to 0 in both bands (indicating that the sources are in unmasked areas).

Independently, we extracted a source catalog from the UltraVISTA DR2 $K_s$-band mosaic using the software SExtractor (Bertin & Arnouts 1996), imposing a detection threshold of $2.5\sigma$ over 5 contiguous pixels. Such a low threshold is appropriate in this case, as we are only aiming to find counterparts for the robustly detected sources in the IRAC bands. We measured aperture photometry in 3-arcsec diameter circular apertures for all the extracted sources in the $K_s$ band and derived aperture corrections from the curves of flux growth of isolated stars in the field. To perform photometric measurements on all the other UltraVISTA broad bands, we used SExtractor in dual-image mode, using the $K_s$-band mosaic as the detection image. The aperture corrections vary between −0.11 and −0.19 mag, depending on the band.

We preferred using 3-arcsec diameter apertures rather than smaller apertures because a priori we do not know the redshifts of our sources, and too small apertures can be inadequate to recover the photometry of low-$z$ galaxies. Nevertheless, we tested the use of 2-arcsec diameter apertures and found no significant difference in our results for high-$z$ galaxies, which we describe below.

The next step was cross-correlating the Spitzer IRAC $[4.5] < 23$ sources with the UltraVISTA $K_s$-band sources. We limited our study to the UltraVISTA ultradeep stripes (~0.8 deg$^2$; Figure 1) to ensure that we worked with near-IR data of homogeneous depth. The overall identification completeness of our IRAC $[4.5] < 23$ sources over these ultradeep stripes is ~99%.

As the aim here is to study the bright IRAC sources that could not be identified before with the UltraVISTA DR1 data

<footnote n="9">⁹ http://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-scan?projshort=COSMOS</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of our [4.5] < 23 sample with $K_s^{\text{auto}} > 24$ counterparts in the UltraVISTA/COSMOS field, where the red rectangles indicate the UltraVISTA ultradeep stripes. The dashed area at the bottom of the left-hand side stripe indicates a region with no data, owing to a faulty detector in the VISTA telescope wide-field camera VIRCAM (Sutherland et al. 2015).</figcaption></figure>

set, we looked for [4.5] < 23 sources with a $K_s^{\text{auto}} > 24$ counterpart. Exclusively for the purpose of identifying sources with faint $K_s$ counterparts, we used a $K_s$ magnitude cut based on the SExtractor “mag_auto,” as this guarantees that we are dealing with a sample complementary to those obtained from the UltraVISTA DR1 images (Ilbert et al. 2013). This selection makes this new sample also complementary to the IRAC galaxy sample studied by Caputi et al. (2011), in which the identification completeness of the [4.5] < 23 galaxies was ~96%, and less than 1% of the identified [4.5] < 23 sources had $K_s^{\text{auto}} > 24$ at all redshifts, and <4% at $z \ge 3$. This complementarity is important to assess the corrections that the newly identified bright IRAC galaxies introduce to the GSMF at $z \ge 3$ (see Sections 4.2 and 5.2). For all other purposes, we used total magnitudes obtained from aperture magnitudes and corresponding aperture corrections.

We identified 604 IRAC [4.5] < 23 sources with $K_s^{\text{auto}} > 24$ counterparts within a 1-arcsec matching radius. To avoid dealing with cases with severely contaminated photometry, we excluded from our sample 30 sources that also have a brighter ($K_s^{\text{auto}} < 24$) neighbor within 2 arcsec radius. Our final, clean, [4.5] < 23 and $K_s^{\text{auto}} > 24$ sample contained 574 sources. Note that, in spite of imposing a matching radius of 1 arcsec, the median separation between the IRAC and $K_s$-band centroids in our sample is only 0.26 arcsec.

We also compiled the CFHT *u*-band and SUBARU optical ancillary data available for our sources, including the COSMOS broadband, narrowband, and intermediate-band photometry. To obtain this multiwavelength photometry, we ran SExtractor once again in dual-image mode, using the UltraVISTA $K_s$-band mosaic as the detection image, with the same source extraction parameters as before. As for the UltraVISTA data, we measured aperture photometry in 3-arcsec diameter apertures and derived aperture corrections from the curves of flux growth of isolated stars in the field.

Our finally compiled catalog contains photometric information for our 574 sources in 30 photometric bands. In total, more

than 75% of our 574 IRAC [4.5] < 23 sources with $K_s^{\text{auto}} > 24$ counterparts are detected in the z, $z^{++}$, and different optical bands (in addition to being detected in the near- and mid-IR). All the remaining sources are detected in at least six IRAC/ UltraVISTA bands. All the magnitudes in our final catalog are total (obtained from corrected aperture magnitudes) and have been corrected for galactic extinction. As in Ilbert et al. (2013), we have multiplied all error bars by a factor of 1.5 to account for underestimated photometric errors in SExtractor’s output, mainly owing to correlations between pixels produced by image resampling.

To identify galactic stars, we used a ($B-J$) versus ($J-[3.6]$) color–color diagram, similarly to Caputi et al. (2011). We found that none of our sources display the characteristic blue ($J-[3.6]$) colors of galactic stars. As we explain in Section 3.2, we additionally checked for the presence of red dwarf contaminants but found none either within our sample.

## 3.2. SED Modeling and Photometric Redshifts

We performed the SED modeling of our sources with a customized $\chi^2$-minimization fitting code, using Bruzual & Charlot (2003) synthetic stellar templates for different star formation histories: a single stellar population, and exponentially declining models with characteristic times ranging between $\tau = 0.1$ and 5.0 Gyr, all with solar metallicity. For each galaxy, we tested all redshifts between $z = 0$ and 7, with a step $dz = 0.02$. To account for the galaxy internal extinction, we convolved the stellar templates with the Calzetti et al. (2000) reddening law, allowing for V-band extinction values $0 \le A_V \le 6$. We used the Madau (1995) prescription to include the effects of the intergalactic medium attenuation at $\lambda_{\text{rest}} < 1216 \AA$. As an output of our code, we obtained each galaxy best-fit photometric redshift ($z_{\text{phot}}$) and derived parameters, including the stellar mass ($M_{\text{st}}$).

In a first step, we performed the SED modeling of our sources over 28 photometric bands, from $u^*$ through the IRAC 4.5 µm band. The reason to exclude the longest-wavelength IRAC bands in this step is that, for low-$z$ sources, dust emission (including that from polycyclic aromatic hydrocarbons) can significantly contribute to the photometry beyond ~5 µm, so the SED fitting with pure stellar templates may not be adequate. As a priori we did not know the redshifts of our sources, we excluded the longest-wavelength data in our first photometric redshift run. For the SED modeling, when a source was nondetected in a given broad band, we rejected any template that produced a flux above the $3\sigma$ detection limit in that band. Narrow and intermediate bands were ignored in the cases of nondetections.

We imposed a maximum redshift for each source, according to its detections at short wavelengths, following the criteria explained in Caputi et al. (2011): for sources detected in the $u^*$ band (with $>2\sigma$ significance), the maximum accepted redshift was $z_{\text{max}} = 3.6$. For sources significantly detected in *Bj*, *Vj*, and *r*⁺, the maximum imposed redshifts were 4.6, 5.6, and 6.4, respectively. We note that, in practice, these maximum redshift constraints had to be applied for a small amount of sources (~2% of the sample), as our code directly produced best-fit redshifts consistent with the short-wavelength photometric detections in the vast majority of cases.

For all the $z_{\text{phot}} \ge 3$ galaxy candidates obtained in the SED modeling run (a total of 252 sources), we performed a second
SED modeling considering the 30 photometric bands (i.e., with the full IRAC photometry). This run with full IRAC photometry was only performed in the cases of reliable IRAC 5.8 and 8.0 µm detections, i.e., only for magnitudes <22.5 and <22.0 AB mag, at 5.8 and 8.0 µm, respectively (which corresponds to >3σ detections). We analyzed the SED fitting results for these sources considering the results of the 30-band SED modeling, and also the full χ² map in all cases, in order to investigate degeneracies in $A_V-z$ space.

To accept a $z_{phot} \ge 3$ galaxy candidate, we imposed the criterion that the best-fit solution should be $z_{phot} \ge 3$ both in the 28- and the 30-band SED fitting, when both were available. For the $z_{phot} \ge 5$ candidates, we also requested two additional conditions, namely, (i) no secondary χ² local minimum should exist at $z_{phot} < 5$ within $1\sigma$ of the best-fit solution at $z \ge 5$ (i.e., reduced χ² (secondary $z_{phot}$)-χ² (primary $z_{phot}$) > 1); (ii) the median of the marginalized $P(z)$ versus $z$ distribution should also be at $z_{phot} \ge 5$. The latter is especially helpful in cases of broad $P(z)$ distributions, where no secondary χ² local minimum is identified within $1\sigma$, but $P(z)$ significantly extends to lower $z$. These are conservative criteria that ensure that we keep only reliable high-$z$ sources in our sample.

Overall, we found that 188 out of 252 $z_{phot} \ge 3$ candidates satisfy all these criteria. For the remaining candidates, we replaced the best-fitting solutions by the 30-band best solutions, when available, the secondary lower $z_{phot}$ within $1\sigma$ confidence, or the lower median $z_{phot}$ values, with corresponding best-fit parameters.

In the analysis of the $z_{phot} \ge 3$ galaxy candidates with full IRAC photometry, we searched for galaxies with a plausible IR power-law SED component, following the power-law subtraction methodology explained in Caputi (2013). Such a power-law SED component is a clear indication of active galactic nuclei (AGNs) at such high $z$, but only the warmest sources at $z \ge 3$ can be identified using data up to 8 µm (the maximum contribution of the AGN-associated IR power law occurs at rest wavelengths 1–2 µm, and thus this maximum is shifted beyond observed ~8 µm at $z_{phot} \ge 3$; see discussion in Caputi 2013, 2014). Nevertheless, our aim here is not to make a complete AGN census among our galaxies, but rather to identify cases where the $z_{phot}$ and stellar masses could be affected by a nonstellar component in the SED. So, it is interesting to identify which of our high-$z$ galaxies do manifest an IR power-law signature. Among our $z_{phot} \ge 3$ candidate galaxies, only three display a mid-IR excess that suggests an AGN presence (and one of these has a revised redshift $z_{phot} < 3$).

Recently, it has been pointed out that the presence of emission lines can produce contamination in the selection of high-$z$ galaxies (Zackrisson et al. 2008; Schaerer & de Barros 2009). However, this problem typically arises for galaxies with blue, rather than red, near-/mid-IR colors (de Barros et al. 2014). We investigated the possible presence of low-redshift interlopers in our $z_{phot} \ge 3$ galaxy sample produced by line emitters by doing an independent photometric redshift run using the public code LePhare (Arnouts et al. 1999; Ilbert et al. 2006), including emission lines. We found that only two of our $z_{phot} \ge 3$ galaxies have a significantly lower best-fit redshift (1.25 < $z$ < 2.6), while all the rest are confirmed to be at $z_{phot} \ge 3$. We adopted the lower redshifts for the two interlopers found with SED fitting with emission lines.

As a summary, we have 185 galaxies at $z_{phot} \ge 3$, which make ~32% of our total [4.5] < 23, $K_s^{auto} > 24$ sample.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Redshift distribution of our [4.5] < 23, $K_s^{auto} > 24$ sample (thick solid histogram). For a comparison, we also show the renormalized redshift distributions of the [4.5] < 23, $K_s^{auto} < 24$ galaxies (thin solid histogram), and all $K_s^{auto} < 24$ galaxies (thin dashed histogram).</figcaption></figure>

Among these high-$z$ sources, we identified nine reliable galaxy candidates at $z_{phot} \ge 5$.

None of our $z_{phot} \ge 3$ sources appear to be a red dwarf contaminant. We modeled the SEDs of these galaxies using the characteristic stellar templates of M and L stars (Rayner et al. 2009), but in no case do these templates produce a good fitting of the observed light. Our $z_{phot} \ge 3$ sample does not contain later-type dwarf contaminants either. Most T dwarfs and the Y0 dwarfs recently discovered by the *Wide-field Infrared Survey Explorer* are bright in the mid-IR, but bright in the near-IR as well, and typically have [3.4] – [4.6] $\gtrsim$ 1.5 (Eisenhardt et al. 2010; Kirkpatrick et al. 2011). Our sources are near-IR faint and all have [3.6] – [4.5] $\lesssim 1$.

### 3.3. Redshift Distribution and IR Colors

Figure 2 shows the resulting redshift distribution of our [4.5] < 23, $K_s^{auto} > 24$ sample. Virtually all our galaxies lie at redshifts $z_{phot} > 1$, with a clear peak at redshifts 2.5 < $z_{phot} < 3.0$. This bias toward high redshifts is a consequence of the red colors imposed by our double magnitude-limited selection. This effect is similar to that observed in classical “extremely red galaxies,” whose redshifts typically are $z \gtrsim 1.0$ (Caputi et al. 2004; Georgakakis et al. 2006), and dust-obscured galaxies, which mainly lie at $z \gtrsim 2$ (Dey et al. 2008; Pope et al. 2008).

In addition, the redshift distribution has a significant tail at $z_{phot} \ge 3$, containing ~32% of the galaxies in our sample. In this high-$z$ tail, there are nine galaxies at redshifts $z_{phot} \ge 5.0$.

About 10% of our galaxies have no $z_{phot}$ determination. There are two reasons for this: for 2% of our sample, the minimum χ² value is too high to trust the resulting $z_{phot}$. For the other 8% of sources, the resulting probability distribution in redshift space $P(z)$ is basically flat at $z \gtrsim 2$, so it is not possible to decide on any redshift.

As our galaxies are faint at optical and near-IR wavelengths, it is currently very difficult to obtain spectroscopic confirmations for our redshifts. We have cross-correlated our catalog with a compilation of all available sources with spectroscopic redshifts in COSMOS, including those in the zCOSMOS (Lilly
<figure><img src="image_4.png" /><figcaption>**Figure 3.** $K_s$ – [4.5] vs. $H - K_s$ color–color diagram for our [4.5] < 23, $K_s^{\text{auto}}$ > 24 sources at different redshifts. The dashed line delimits the region corresponding to $H - [4.5] > 4$ colors, which is the color cut imposed for the sample selection in Caputi et al. (2012). Median error bars on the colors are shown in the bottom left corner of the plot.</figcaption></figure>

et al. 2007) and VUDS surveys (Le Fèvre et al. 2015), and data taken with the DEIMOS and FMOS spectrographs (P. Capak et al. 2015, in preparation; J. Kartaltepe et al. 2015, in preparation), among others. We found a total of six matches, and in four of these cases our photometric redshifts are in very good agreement with the spectroscopic redshifts within the error bars (these redshifts range between $z = 2.10$ and $z = 3.88$). The other two sources have low signal-to-noise spectra, with (only tentative) redshift measurements that do not agree with our photometric determinations (which are $z_{phot} = 1.68$ and 4.08). As the spectral quality of these two sources is not really useful for a diagnostic, we considered that our photometric redshift estimates are likely correct.

In Figure 3, we show a $K_s$ – [4.5] versus $H - K_s$ color–color plot, where we have segregated our galaxies by redshift. We see that our galaxies span a wide range of $K_s$ – [4.5] colors, including a minority with $K_s$ – [4.5] $\le 1$. This is because our colors are based on total magnitudes (obtained from aperture magnitudes, plus corrections), while the $K_s^{\text{auto}}$ > 24 cut applied in our sample selection refers to the SExtractor $K_s$ “mag_auto.” (The median difference between the $K_s$ total and the $K_s^{\text{auto}}$ magnitudes in our sample is about –0.07 mag.)

The $K_s$ – [4.5] color, in particular, is only mildly related to the redshift for the [4.5] < 23 galaxies. About two-thirds of the sources with $K_s$ – [4.5] > 2 lie at $z_{phot} \ge 3$, and the other third are mostly $2 < z_{phot} < 3$ galaxies with extinction values $A_V \approx 3.0$. The $H - K_s$ color is potentially a better redshift discriminator: at least two-thirds of the $H - K_s \gtrsim 1$ sources lie at $3 \le z_{phot} < 5$, and only 13% are confirmed at lower redshifts (there are about 20% of sources with these colors that have no redshift estimates in our sample). The $z_{phot} \ge 5$ galaxy candidates, instead, are characterized by flat near-/mid-IR

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Estimated stellar masses vs. photometric redshifts for our [4.5] < 23, $K_s^{\text{auto}}$ > 24 galaxies at $3 \le z < 5$ (red filled circles). The data points in the background indicate [4.5] < 23, $K_s^{\text{auto}}$ < 24 galaxies from two different samples: the UDS sample analyzed by Caputi et al. (2011) (small black dots) and the UltraVISTA DR1 sample analyzed by Ilbert et al. (2013) (crosses). The stellar masses derived by Caputi et al. (2011) have been multiplied by a factor of 1.24 to obtain a crude conversion from BC07 to BC03 templates, while the Ilbert et al. (2013) stellar masses have been multiplied by a factor of 1.7 to convert from a Chabrier to a Salpeter IMF over (0.1–100) $M_{\odot}$. The horizontal lines indicate the stellar mass completeness limits imposed by the [4.5] < 23 mag cut at $z \sim 3.5$ and 4.5.</figcaption></figure>

spectra, with $H - K_s$ colors around 0, similarly to most $z_{phot} < 3$ galaxies.

Only a small fraction of our sources display the extremely red $H - [4.5] > 4$ colors characterizing the galaxies studied by Caputi et al. (2012). This is mainly a consequence of the insufficient depth of the UltraVISTA DR2 images to probe such red sources. The bulk of these extremely red galaxies are expected to be among the ~1% of [4.5] < 23 sources that remain unidentified in the UltraVISTA DR2. The sources that do satisfy $H - [4.5] > 4$ in our current sample appear to be at $3 \le z < 5$, consistently with the results of Caputi et al. (2012).

The majority of the ~10% of sources with no redshift estimate in our sample display similarly red $H - K_s$ colors to the $3 \le z < 5$ sources, so they also likely lie at these redshifts. Only a few sources with no redshift estimates display similar colors to the $z = 5-6$ candidates, but in this case we cannot conclude on a plausible high $z$, as many sources at $z < 3$, and a minority at $3 \le z < 5$, also have similar colors.

# 4. MASSIVE GALAXIES AT $3 \le z < 5$

## 4.1. Properties and Comparison with Near-IR Brighter, Massive Galaxies at Similar Redshifts

A total of 176 galaxies in our [4.5] < 23, $K_s^{\text{auto}}$ > 24 sample lie at redshifts $3 \le z < 5$. Figure 4 shows the derived stellar masses versus photometric redshifts for all these galaxies. We see that the stellar masses range from ~10⁹ to ~6 × 10¹¹ $M_{\odot}$, and about 66% of our galaxies are quite massive, with $M_{\text{st}} > 5 × 10^{10} M_{\odot}$. At redshifts $z > 4$, almost 90% of our galaxies display such large stellar masses. The horizontal lines centered at $z = 3.5$ and 4.5 in Figure 4 indicate the stellar mass completeness limits imposed by the [4.5] < 23 mag cut at these redshifts. These limiting stellar masses correspond to galaxies described by a single-stellar-population template, with no dust,
and with an age equal to the age of the universe at those redshifts.

For a comparison, in the same plot we show the stellar masses obtained for [4.5] < 23 galaxies with brighter $K_s$ counterparts in COSMOS and the Ultra Deep Survey (UDS) field. By simple inspection, it becomes clear that the new galaxies analyzed here contribute significantly to the popula- tion of massive galaxies at $z \ge 3$. This is especially true at $4 \le z < 5$, where our newly analyzed [4.5] < 23, $K_s^{\text{auto}} > 24$ galaxies increase by a factor $\gtrsim 2$ the number density of known massive galaxies. In the next section, we investigate the effects that this substantial increment has on the GSMF at these redshifts.

The maximum extinction value associated with the best-fit solutions of our $3 \le z < 5$ sources is $A_V = 2.7$, and the median value is of only $A_V = 1.2$. Thus, we conclude that the extinction values keep moderate even for the reddest sources in our sample (with respect to the maximum $A_V$ value allowed in our SED fitting grid, i.e., $A_V = 6$).

It is interesting to compare these extinctions and other average SED properties with those obtained previously for massive galaxies in other sample selections at $3 \le z < 5$. For the extremely red sources with $3 \le z < 5$ in the Caputi et al. (2012) sample, the median derived extinction was $A_V = 2.2$, which is significantly higher than the median value found for our sample here, but still relatively moderate. In addition, the $3 \le z < 5$ galaxies in Caputi et al. (2012) had associated best-fit ages $\gtrsim 0.5$ Gyr in all cases, while around one-quarter of the galaxies in our current $3 \le z < 5$ sample have best-fit younger ages. On the other hand, in the case of the [4.5] < 23 galaxies at $3 \le z < 5$ included in the Caputi et al. (2011) sample, which had no color restriction, the median extinction was $A_V = 0.8$, and 60% of the sample had best-fit ages younger than $\sim 0.5$ Gyr.

Thus, we conclude that the colors of massive galaxies at $3 \le z < 5$ are the result of a combination of age and dust extinction: redder galaxies are typically both older and dustier than bluer ones.

## 4.2. The GSMF at $3 \le z < 5$

### 4.2.1. An Updated Calculation of the GSMF
High-mass End at $3 \le z < 5$

A main goal of this work is to assess the importance of the newly identified $3 \le z < 5$ galaxies with respect to the population of massive galaxies known at these redshifts. For this, we analyze their contribution to the high-mass end of the GSMF. We note that for this analysis we did not exclude any source as a result of the presence of an IR AGN power-law component. We rather considered corrected stellar masses, computed after power-law subtraction from the SED (see Caputi 2013), in the few necessary cases discussed before.

Figure 5 shows the GSMF at $3 \le z < 4$ and $4 \le z < 5$. The red circles in Figure 5 and values in Table 1 indicate the updated $V/V_{\text{max}}$ GSMF values, obtained by adding the contribution of our [4.5] < 23, $K_s^{\text{auto}} > 24$ galaxies to the average GSMF values obtained by other recent analyses of large-area, near-/mid-IR surveys. These previous GSMF determinations are based on the IRAC survey of the UDS field (Caputi et al. 2011) and UltraVISTA DR1 in COSMOS (Ilbert et al. 2013). From all these studies, the analyzed sources are complementary to our current sample (i.e., they have

$K_s^{\text{auto}} < 24$), and the differently obtained GSMFs are in very good agreement among themselves. In the case of the UDS data, we have recomputed the GSMF values in the redshift bins considered here after excluding a minor fraction of sources with $K_s^{\text{auto}} > 24$ in the Caputi et al. (2011) sample.

For a correct joint analysis and comparison, we have converted the UltraVISTA previous determinations to the same IMF used here, i.e., Salpeter over stellar masses (0.1–100) $M_{\odot}$, and the stellar mass values from Caputi et al. (2011) have been multiplied by a factor of 1.24 to provide a crude correction from the 2007 to the 2003 version of the Bruzual & Charlot templates. To obtain the average GSMF of these previous surveys, we have linearly interpolated the UltraVISTA DR1 values at the stellar mass bin centers considered here.

To consider the $V/V_{\text{max}}$ contribution of each of our new galaxies, we took into account the fact that they are selected with a magnitude limit [4.5] < 23. This criterion imposes a maximum redshift at which the galaxy would be in the sample. On the other hand, the $K_s^{\text{auto}} > 24$ criterion imposes a minimum redshift at which the galaxy would be included. However, since here we are analyzing jointly previously selected $K_s^{\text{auto}} < 24$ galaxy samples and our new sample with $K_s^{\text{auto}} > 24$, this extra correction is not necessary.

The red solid lines in Figure 5 and parameter values in Table 2 indicate the result of an independent computation of the GMSF, obtained by applying the STY (Sandage et al. 1979) maximum likelihood analysis, and assuming that the GSMF has the shape of a single Schechter function (Schechter 1976). For this GSMF computation, we considered the complementary IRAC galaxy sample in the UDS from Caputi et al. (2011), along with the new UltraVISTA COSMOS sample analyzed here (each with their corresponding magnitude limits and surveyed areas). The STY method is a parametric technique that, in contrast to the $V/V_{\text{max}}$ method, involves no data binning and has no implicit assumption on a uniform galaxy spatial distribution. However, it has the disadvantage that it does not automatically provide the normalization of the GSMF, which rather has to be provided as an input parameter. Here we considered the number density of galaxies with $M_{\text{st}} > 10^{11} M_{\odot}$ in our sample at different redshifts to compute this normalization. It is worth reminding the reader that the curves resulting from the STY method do not constitute a fitting to the $V_{\text{max}}$ points. We refer the reader to Caputi et al. (2011) for more details about the GSMF computation using the STY method.

Table 2 contains the Schechter function parameter values obtained with the maximum likelihood STY analysis. In this analysis, we have left both the $\alpha$ and $M^*$ parameters free. The values of $\Phi_M^*$ and $\rho_M$ carry error bars that include the largest errors among the $1\sigma$ errors of the maximum likelihood analysis and the mock realizations described in Section 4.2.3, and a fixed 20% fiducial error to account for cosmic variance (following the determinations of Ilbert et al. 2013 for massive galaxies).

We show the $1\sigma$, $2\sigma$, and $3\sigma$ confidence levels on the $\alpha-M^*$ plane in Figure 6. There is the misconception that $\alpha$, commonly referred to as the faint-end slope, cannot be left as a free parameter unless the low-mass end of the GSMF can be well constrained. However, this is not the case, for two reasons: first, the value of $\alpha$ affects exclusively the GSMF low-mass end only when it is close to 1, as the Schechter function is proportional to $(M/M^*)^{(1-\alpha)} \times \exp(-M/M^*)$. In this case, the effect on the high-mass end is basically negligible. Instead, when $\alpha$ is
<figure><img src="image_6.png" /><figcaption>**Figure 5.** GSMF at 3.0 ≤ z < 4.0 (left) and 4.0 ≤ z < 5.0 (right), obtained by considering the results of previous $K_s^{auto} < 24$ surveys and the new [4.5] < 23, $K_s^{auto} > 24$ galaxies analyzed here. The large filled circles correspond to the total GSMF, while the small open circles show the contribution of only our [4.5] < 23, $K_s^{auto} > 24$ galaxies, as computed with the $V/V_{max}$ method. Other symbols correspond to previous $K_s^{auto} < 24$ surveys: diamonds (Caputi et al. 2011) and squares (Ilbert et al. 2013), with the dot-dashed line on the left-hand side panel indicating the average between the two determinations at 3.0 ≤ z < 4.0. The GSMF determination by Muzzin et al. (2013; asterisks) is also shown for a comparison. All $V/V_{max}$ data points correspond to a Salpeter IMF over (0.1–100) $M_{\odot}$, and the stellar masses from Caputi et al. (2011) have been multiplied by a factor of 1.24 to obtain a crude conversion from BC07 to BC03 templates. The solid red line in each panel indicates the resulting GSMF obtained with the STY maximum likelihood analysis, assuming a single Schechter function. This maximum likelihood analysis has been performed on a combination of the Caputi et al. (2011) sample and the new galaxies analyzed in this paper. The vertical dotted line in each panel indicates the stellar mass completeness limit imposed by the IRAC [4.5] = 23 mag cut. The dashed line indicates the local GSMF (Cole et al. 2001).</figcaption></figure>

<table><caption>**Table 1**
GSMF Values Obtained with the $V/V_{max}$ Method at $3 \le z < 5$</caption><thead><tr><th rowspan="2">log<sub>10</sub>(M/M<sub>☉</sub>)</th><th colspan="2">log<sub>10</sub>(Φ<sub>M</sub>/Mpc<sup>-3</sup> dex<sup>-1</sup>)</th></tr><tr><th>3.0 ≤ z < 4.0</th><th>4.0 ≤ z < 5.0</th></tr></thead><tbody><tr><td>10.40</td><td>&gt; -3.56</td><td>&gt; -4.44</td></tr><tr><td>10.60</td><td>&gt; -3.65</td><td>&gt; -4.46</td></tr><tr><td>10.80</td><td>-3.82<sup>+0.17</sup><sub>-0.16</sub></td><td>&gt; -4.52</td></tr><tr><td>11.00</td><td>-3.97<sup>+0.19</sup><sub>-0.21</sub></td><td>-4.51<sup>+0.17</sup><sub>-0.21</sub></td></tr><tr><td>11.20</td><td>-4.25<sup>+0.20</sup><sub>-0.25</sub></td><td>-4.65<sup>+0.19</sup><sub>-0.25</sub></td></tr><tr><td>11.40</td><td>-4.51<sup>+0.25</sup><sub>-0.31</sub></td><td>-4.91<sup>+0.20</sup><sub>-0.27</sub></td></tr><tr><td>11.60</td><td>-4.94<sup>+0.30</sup><sub>-0.53</sub></td><td>-5.20<sup>+0.31</sup><sub>-0.47</sub></td></tr><tr><td>11.80</td><td>-5.34<sup>+0.30</sup><sub>-0.61</sub></td><td>-5.88<sup>+0.35</sup><sub>-0.70</sub></td></tr></tbody></table>

sufficiently larger than 1, it has two effects: it governs the low-mass end slope and also modulates the exponential decline in the high-mass end.

Second, the STY method takes into account the survey limiting magnitude, providing a suitable extrapolation to stellar masses below completeness when the regime around and below $M^*$ is reasonably well sampled (see Figure 5). This does not happen when a Schechter function is used to do a simple fitting of the $V/V_{max}$ data points, which is the most typical case in the literature. Of course, having a galaxy sample that represents significantly both the low- and high-mass end would be the ideal situation for computing the Schechter function parameter values. Sampling a wide range of stellar masses is also necessary to probe, for instance, whether the GSMF follows a single or double Schechter function, but this is beyond the scope of this paper. In spite of not reaching the low stellar mass

end, our sample stellar mass completeness limits are below $M^*$ at $3 \le z < 5$, so the computation of the Schechter $\alpha$ parameter is meaningful. The agreement of our maximum likelihood analysis results with previous literature values suggests that this is indeed the case.

In Figure 7, we show the relative contribution of our [4.5] < 23, $K_s^{auto} > 24$ galaxies to the total GSMF. At $3 \le z < 4$, we see that our new galaxies produce a basically negligible correction to the high-mass end of the GSMF, indicating that a $K_s^{auto} < 24$ survey contains the vast majority of massive galaxies at these redshifts. In fact, as we discuss below, deeper near-IR surveys appear to have a relatively small additional contribution even at intermediate stellar masses ($10^{10} \lesssim M_{\text{st}} \lesssim 10^{11} M_{\odot}$) at these redshifts.

Instead, the contribution of the [4.5] < 23, $K_s^{auto} > 24$ galaxies to the GSMF high-mass end at $4 \le z < 5$ appears to be more significant. Notably, it allows us to constrain the GSMF at very high stellar masses ($M_{\text{st}} \gtrsim 3 \times 10^{11} M_{\odot}$), which was not possible with the sample of IRAC galaxies with brighter $K_s$ counterparts in Caputi et al. (2011). We have only 11 galaxies with such high stellar masses within our sample at $4 \le z < 5$ in the UltraVISTA ultradeep stripes, which indicates that these galaxies are indeed rare. Overall, our new galaxies account for $\gtrsim 50\%$ of the galaxies with stellar masses $M_{\text{st}} \gtrsim 6 \times 10^{10} M_{\odot}$ at $4 \le z < 5$. Therefore, we conclude that finding the bulk of massive galaxies at $z \gtrsim 4$ requires ultradeep near-IR surveys covering large areas of the sky.

Our present STY analysis based on the combination of the UDS IRAC galaxy sample and our current UltraVISTA [4.5] < 23, $K_s^{auto} > 24$ galaxies yields $\alpha = 1.72 \pm 0.06$ at $3 \le z < 4$ and $\alpha = 1.88_{-0.18}^{+0.12}$ at $4 \le z < 5$, confirming that
<table><caption>**Table 2**

Schechter Function Free Parameter Values Obtained for the GSMF Computed with the Maximum Likelihood STY Analysis</caption><thead><tr><th>Redshift</th><th>α</th><th>M*(M<sub>☉</sub>)</th><th>Φ<sup>*</sup><sub>M</sub>(Mpc<sup>-3</sup> dex<sup>-1</sup>)</th><th>ρ<sup>compl.</sup><sub>M</sub>(M<sub>☉</sub> Mpc<sup>-3</sup>)<sup>a</sup></th><th>ρ<sup>TOTAL</sup><sub>M</sub>(M<sub>☉</sub> Mpc<sup>-3</sup>)<sup>b</sup></th></tr></thead><tbody><tr><td>3.0 &le; z &lt; 4.0</td><td>1.72 &plusmn; 0.06</td><td>(1.82<sup>+0.18</sup><sub>-0.16</sub>) &times; 10<sup>11</sup></td><td>(1.62<sup>+0.85</sup><sub>-0.48</sub>) &times; 10<sup>-4</sup></td><td>(1.11<sup>+0.58</sup><sub>-0.33</sub>) &times; 10<sup>7</sup></td><td>(3.56<sup>+1.87</sup><sub>-1.06</sub>) &times; 10<sup>7</sup></td></tr><tr><td>c</td><td>1.58 &plusmn; 0.06</td><td>(1.79<sup>+0.16</sup><sub>-0.21</sub>) &times; 10<sup>11</sup></td><td>(1.23<sup>+0.95</sup><sub>-0.26</sub>) &times; 10<sup>-4</sup></td><td>(7.86<sup>+6.03</sup><sub>-1.65</sub>) &times; 10<sup>6</sup></td><td>(1.92<sup>+1.47</sup><sub>-0.40</sub>) &times; 10<sup>7</sup></td></tr><tr><td>4.0 &le; z &lt; 5.0</td><td>1.88<sup>+0.12</sup><sub>-0.18</sub></td><td>(2.40<sup>+0.48</sup><sub>-0.74</sub>) &times; 10<sup>11</sup></td><td>(4.39<sup>+2.49</sup><sub>-1.09</sub>) &times; 10<sup>-5</sup></td><td>(3.41<sup>+1.93</sup><sub>-0.84</sub>) &times; 10<sup>6</sup></td><td>(2.10<sup>+1.19</sup><sub>-0.52</sub>) &times; 10<sup>7</sup></td></tr><tr><td>c</td><td>1.58<sup>+0.12</sup><sub>-0.16</sub></td><td>(2.00<sup>+1.16</sup><sub>-0.59</sub>) &times; 10<sup>11</sup></td><td>(6.14<sup>+3.55</sup><sub>-1.52</sub>) &times; 10<sup>-5</sup></td><td>(3.32<sup>+1.92</sup><sub>-0.81</sub>) &times; 10<sup>6</sup></td><td>(1.07<sup>+0.62</sup><sub>-0.26</sub>) &times; 10<sup>7</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Stellar mass density values obtained by integrating the resulting Schechter functions above stellar mass completeness.

<sup>b</sup> Stellar mass density values obtained by integrating the resulting Schechter functions for stellar masses $M_{st} > 10^8 M_{☉}$.

<sup>c</sup> The second row for each redshift bin provides the values obtained after correction for Eddington bias.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** 1σ, 2σ, and 3σ confidence levels for the maximum likelihood free parameters at 3.0 ≤ z < 4.0 (top) and 4.0 ≤ z < 5.0 (bottom). The dashed curves correspond to the original maximum likelihood analysis performed considering a Schechter function, i.e., with no correction for Eddington bias. The solid curves correspond to the Eddington-bias-corrected maximum likelihood analysis, i.e., taking into account the convolution of the Schechter function with Gaussian kernels.</figcaption></figure>

α appears to be significantly higher at high z than at z = 0, as was found in several previous studies (e.g., Kajisawa et al. 2009; Caputi et al. 2011; Santini et al. 2012). The α value that we obtained at 3 ≤ z < 4 is slightly lower (in absolute value) than that obtained in Caputi et al. (2011). This is mainly because we have improved the statistics here

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Relative contribution of our new [4.5] < 23, $K_s^{auto} > 24$ galaxies to the total GSMF at different stellar masses, at both 3.0 ≤ z < 4.0 and 4.0 ≤ z < 5.0. At 3.0 ≤ z < 4.0, the contribution is small at all stellar masses, while at 4.0 ≤ z < 5.0 it becomes important at stellar masses above ~6 × 10<sup>10</sup> M<sub>☉</sub>.</figcaption></figure>

by considering broader redshift bins (in Caputi et al. 2011 we used three different redshift bins, i.e., 3.0 ≤ z < 3.5, 3.5 ≤ z < 4.25, and 4.25 ≤ z < 5.0, for the GSMF computation). At 4 ≤ z < 5, our obtained α value (α = 1.88<sup>+0.12</sup><sub>-0.18</sub>) is very similar to that obtained by Caputi et al. (2011) at similar redshifts, even when the new galaxies make a significant contribution to the GSMF at these redshifts. We do obtain a higher M* value than Caputi et al. (2011), although still consistent within the error bars.

Caputi et al. (2011) found that a single power-law shape and a Schechter function could not be differentiated as a functional form for the GSMF at 4.25 ≤ z < 5.0 in the maximum likelihood STY analysis. And they argued that this could be due to the insufficient sampling of the GSMF at these redshifts. In our present analysis of the GSMF at 4.0 ≤ z < 5.0, a single power-law shape is discarded with >5σ confidence with respect to the Schechter function, indicating that the latter functional form produces a better representation of the GSMF shape up to at least z = 5.

### 4.2.2. Comparison with Other Recent Studies

It is interesting to compare our new results for the GSMF at 3.0 ≤ z < 5.0 with the latest determinations found in the literature. Figure 8 shows our obtained GSMF along with that
<figure><img src="image_9.png" /><figcaption>**Figure 8.** GSMF obtained in this paper at 3.0 ≤ z < 4.0 and 4.0 ≤ z < 5.0, compared to those obtained with data sets from the CANDELS survey (Duncan et al. 2014; Grazian et al. 2015) and other recent works (Santini et al. 2012; Stefanon et al. 2015). No GSMF determination is available for CANDELS and Stefanon et al. (2015) at z < 4, so we only show their results at z ≥ 4. As in Figure 5, the vertical dotted line in each panel indicates the stellar mass completeness limit imposed by our IRAC magnitude cut. Stellar masses correspond to a Salpeter IMF over (0.1–100) M⊙ in all cases.</figcaption></figure>

obtained using ultradeep near-IR data sets from the CANDELS survey (Duncan et al. 2014; Grazian et al. 2015) and other recent works (Santini et al. 2012; Stefanon et al. 2015). We show the computed V/Vmax points in the different cases, as this is the clearest and most direct way of comparing two GSMF determinations.

We see that our GSMF V/Vmax values above our stellar mass completeness thresholds are in good agreement with the values obtained with CANDELS at similar redshifts, and also the Stefanon et al. (2015) values, within the error bars. Note that, in the upper panel of Figure 8, our V/Vmax points are systematically higher than the CANDELS values above our

stellar mass completeness limit because the mean redshifts of the samples are slightly different (z ~ 3.5 and 4, respectively).

In contrast, at high stellar masses, the V/Vmax points of Santini et al. (2012) are systematically higher than all other determinations (although still consistent within the error bars at 3 ≤ z < 4). Their GSMF has been computed analyzing a pencil-beam survey on the Chandra Deep Field South Early Science Release area, which is too small to study the GSMF high-mass end. At Mst < 10^{11} M⊙, Santini et al. V/Vmax points appear to be in better agreement with the other results, including our own.

The comparison of the individual Schechter function free parameters is also a common practice in the literature, but one should always keep in mind that α and M* are coupled, so the discussion of each of them separately has to be taken with care. Often, a comparison of these parameters separately is misleading, resulting in wrong conclusions about apparent discrepancies between different studies. With this caveat in mind, we compare the α values derived here with those obtained in CANDELS.

We obtained α = 1.72 ± 0.06 at z ~ 3.5 and α = 1.88^{+0.12}_{-0.18} at z ~ 4.5. These values are consistent with those derived by Duncan et al. (2014) and Grazian et al. (2015), within the error bars. At z ~ 4, Duncan et al. obtained α = 1.89^{+0.15}_{-0.13}, while Grazian et al. derived α = 1.63 ± 0.05. At z ~ 5, the values derived by these authors are 1.74^{+0.41}_{-0.29} and 1.63 ± 0.05, respectively. Note that Grazian et al. quoted values are those obtained after correction for Eddington bias, which explains why they are lower than the others. As we explain in Section 4.2.3, our own α values also become similarly lower after considering this correction.

Therefore, all these results taken together confirm the now well-known fact that the derived α (absolute) values are significantly larger at high redshifts than in the local universe. The dispersion between the mean α values derived from different data sets is still at the ~10%–20% level. This is due to the selection effects intrinsic to the different analyzed samples and also the different ways of deriving the Schechter function characteristic parameter values. Here we derive these values through an STY maximum likelihood analysis, while Duncan et al. (2014) and Grazian et al. (2015) simply fitted the V/Vmax points. In the latter, the highest stellar mass bins, with large error bars, are basically ignored, while the STY method provides a fairer weight for galaxies of different stellar masses.

### 4.2.3. Error Analysis and the Eddington Bias in the GSMF Determination

The errors in the photometric redshifts and stellar masses of our galaxies propagate into uncertainties in the GSMF determination. To analyze this effect, we repeated the STY maximum likelihood 100 times on different mock galaxy catalogs. Each of these mock catalogs contains the “same” galaxies as the original catalog, but with randomly generated redshifts and stellar masses within the error bars. To assign a photometric redshift for each galaxy in the mock catalog, we computed a random value following a Gaussian distribution around the real, best-fit photometric redshift, with an rms given by the 1σ confidence interval of the P(z) distribution resulting from each galaxy SED fitting. This treatment is adequate for red sources, as the typical P(z) distribution is broad around a single best-redshift peak, rather than having multiple peaks with similar probability.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Total stellar mass errors vs. stellar mass for the $3.0 \le z < 5.0$ galaxies in our sample. These errors are the sum in quadrature of the stellar mass errors produced by the $z_{phot}$ uncertainties and the errors at fixed $z_{phot}$. In most cases, the total errors are dominated by the component due to $z_{phot}$ uncertainties. The stellar mass errors at fixed redshift only have a contribution of ~0.06 dex, while most total values are larger than this number.</figcaption></figure>

The $z_{phot}$ errors automatically produce variations on the derived stellar masses, which we recomputed consistently at each new mock redshift. Hence, we considered that the $\pm 1\sigma$ error of each stellar mass, due to the redshift variations, was given by the values encompassing 68% of the resulting stellar masses for the mock $z_{phot}$ of each galaxy, around the “real” stellar mass value.

In addition to this, we considered that the stellar mass of each galaxy was prone to an additional error, at fixed redshift, given by a Gaussian distribution with a fixed width of $0.30 \times M_{st}$, i.e., $\sigma = 0.15 \times M_{st}$. This is the maximum error on $M_{st}$ observed for $\ge 68\%$ of our galaxies at fixed redshift, similarly to the findings in the Caputi et al. (2011) sample. In fact, this additional error constitutes a relatively minor contribution to the total stellar mass error budget, which is dominated by the errors produced by the photometric redshift uncertainties in most cases, as can be seen in Figure 9.

At $3 \le z < 4$, we found that the resulting scatter in the GSMF determination is very small. Instead, at $4 \le z < 5$, the scatter observed in the GSMF is more important (Figure 10). The $\alpha$ values obtained in the different realizations of the GSMF vary between $\alpha = 1.64$ and $2.02$. We adjusted the error bars of our $\alpha$ value quoted in Table 2 to account for this possible range (which is larger than the range given by the $1\sigma$ errors of the maximum likelihood fitting on the real catalogs).

On the other hand, the GSMF determination can potentially be affected by the so-called Eddington bias (Eddington 1913). This is a consequence of the errors in the photometric redshifts and stellar mass estimates, which introduce non-negligible scatter in the GSMF. This effect becomes significant when the errors are very large and/or the number of galaxies considered to compute the GSMF is low. Because of the latter, the Eddington bias affects mostly the GSMF high- and low-mass

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Full range of GSMF at $4.0 \le z < 5.0$ obtained from 100 mock galaxy catalogs (filled orange dots). The red solid line indicates the GSMF obtained applying the STY maximum likelihood analysis on our real sample. This plot reflects the uncertainties on the GSMF produced by the $z_{phot}$ and stellar mass errors. The long-dashed line corresponds to the Eddington-bias corrected GSMF.</figcaption></figure>

ends, typically producing a flattening, which may result in an artificially high (absolute) $\alpha$ value.

To investigate the effect of the Eddington bias, we considered that the observed function describing the GSMF is in fact the convolution of the “real” GSMF Schechter function with a Gaussian kernel, which is characterized by an rms given by the stellar mass error (see, e.g., Zucca et al. 1997; Teerikorpi 2004; Ilbert et al. 2013 for more details). Here we actually considered two Gaussian kernels, whose characteristic rms values correspond to the two stellar mass errors described above (one produced by the $z_{phot}$ and another at fixed redshift), taking into account these errors on each galaxy on an individual basis. We then repeated the STY maximum likelihood analysis on our real galaxy sample, adopting this convolution as the STY functional form.

We found that the best-fit parameters of the “real” Schechter function are somewhat different from those obtained for the STY analysis with a plain Schechter function without error convolution (see Table 2). This is especially the case for the resulting $\alpha$ value: for the “real” Schechter functions, we get $\alpha = 1.58 \pm 0.06$ and $\alpha = 1.58^{+0.12}_{-0.16}$, at $3.0 \le z < 4.0$ and $4.0 \le z < 5.0$, respectively. These values are consistent with those derived by Grazian et al. (2015), after a similar Eddington bias analysis. At the same time, the change we found in the derived value of $M^*$ is small, which implies that the resulting best Schechter curve after error deconvolution does not differ much with respect to the original Schechter curve at the highest-mass end (Figure 10).

This result is indicating that our originally derived GSMF is affected by the Eddington bias. After correction, it still holds that the real Schechter function $\alpha$ value is higher at $z \ge 3$ than in the local universe, but the differences are less dramatic than those resulting from the GSMF analysis taken at face value.
<table><caption>**Table 3**
IRAC Coordinates and Properties of the $z \ge 5$ Galaxy Candidates</caption><thead><tr><th>ID</th><th>R.A.(J2000)</th><th>Decl.(J2000)</th><th>z<sub>phot</sub></th><th>A<sub>V</sub> (mag)</th><th>M<sub>st</sub> (×10<sup>11</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>66110</td><td>09:57:38.80</td><td>+01:44:37.9</td><td>5.48<sup>+1.04</sup><sub>-0.92</sub></td><td>1.80 ± 0.45</td><td>2.95 ± 0.95</td></tr><tr><td>120769</td><td>09:57:41.32</td><td>+01:59:41.5</td><td>5.56<sup>+1.40</sup><sub>-0.68</sub></td><td>1.20 ± 0.30</td><td>1.23 ± 0.39</td></tr><tr><td>144397</td><td>10:00:47.88</td><td>+02:06:09.2</td><td>5.04<sup>+1.96</sup><sub>-2.28</sub></td><td>2.40 ± 0.60</td><td>1.85 ± 0.59</td></tr><tr><td>154286</td><td>10:00:56.69</td><td>+02:08:49.4</td><td>5.00<sup>+1.24</sup><sub>-1.98</sub></td><td>0.30 ± 0.30</td><td>3.35 ± 1.10</td></tr><tr><td>177325</td><td>09:57:28.92</td><td>+02:14:46.6</td><td>5.04<sup>+1.96</sup><sub>-1.80</sub></td><td>2.10 ± 0.60</td><td>1.37 ± 0.44</td></tr><tr><td>185047</td><td>10:01:56.87</td><td>+02:16:51.2</td><td>5.40<sup>+0.80</sup><sub>-0.40</sub></td><td>0.30 ± 0.30</td><td>0.63 ± 0.20</td></tr><tr><td>231699</td><td>09:57:27.65</td><td>+02:29:06.6</td><td>6.04<sup>+0.96</sup><sub>-0.24</sub></td><td>0.50 ± 0.30</td><td>1.77 ± 0.57</td></tr><tr><td>234978</td><td>10:02:12.64</td><td>+02:30:01.0</td><td>5.00<sup>+1.80</sup><sub>-0.96</sub></td><td>1.20 ± 0.45</td><td>2.76 ± 0.88</td></tr><tr><td>260169</td><td>10:01:57.76</td><td>+02:36:48.2</td><td>5.68<sup>+0.16</sup><sub>-0.72</sub></td><td>0.70 ± 0.30</td><td>1.00 ± 0.32</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** $K_s$ and 4.5 µm postage stamps of our $z \ge 5$ galaxy candidates. The field shown in each stamp is of ~12 × 12 arcsec².</figcaption></figure>

## 5. MASSIVE GALAXY CANDIDATES AT $z \ge 5$

### 5.1. Sample Properties

Our [4.5] < 23, $K_s^{auto} > 24$ sample contains a total of nine $z \ge 5$ galaxy candidates, out of which only one is at $z > 6$ (with best $z_{phot} = 6.04$). As explained in Section 3.2, this is a conservative list of objects, for which both the $\chi^2_{min}$ photometric redshifts and median of the $P(z)$ distribution indicate that these are very likely $z \ge 5$ galaxies. We list the coordinates and properties of these objects in Table 3 and show their near- and mid-IR images in Figure 11.

In addition, Figure 12 shows the best-fitting SEDs and marginalized $P(z)$ versus $z$ distributions of these nine galaxies. From these plots, we see that the presence of the 4000 Å break shifted beyond the $K_s$ band, in combination with dust extinction, produces the near- to mid-IR red colors of these sources. In one case, instead, the red color is mainly produced by a 4.5 µm flux excess, which is likely due to the presence of an emission line. For this source, our best $z_{phot}$ value (obtained only with continuum SED fitting) is consistent with the $H\alpha$

emission line ($\lambda_{rest} = 6563$ Å) being shifted into the 4.5 µm band. Note that, as expected, the $P(z)$ distributions have well-localized peaks when the sources are detected in many bands, but become broader when there are multiple nondetections.

The best-fitting extinction values of our $z_{phot} \ge 5$ candidates span a range between $A_V = 0.30$ and 2.40 mag. We have identified one source with an IR power-law excess (id #144397), which, at these redshifts, indicates the presence of a very warm AGN. As before, this has been determined through the SED analysis extended to all IRAC bands, performed following the prescription of Caputi (2013). This suggests that dust-obscured AGNs have been present in the universe since the first few billion years of cosmic time. Indeed, another $z \ge 5$ candidate, source #260169, is detected at 24 µm with flux density $S_ν(24~µm) = (133 \pm 19)~µJy$. The redshift estimate of this source is one of the most secure among our $z \ge 5$ candidates—note that it is detected in most of the near-IR and optical wavebands, except at the shortest wavelengths (Figure 12). Hence, this 24 µm emission is very likely due to an AGN.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Best SED fittings and marginalized probability density distributions $P(z)$ vs. $z$ for our $z \ge 5$ galaxy candidates. For clarity, only the broadband photometric data points are shown in the SEDs.</figcaption></figure>

On the other hand, source #154286 may also have a (marginal) 24 µm detection, but the exact flux density is difficult to determine precisely as it is blended with a bright neighbor at 24 µm (which can be seen at the bottom of the IRAC stamps shown in Figure 11). Hence, for this source, there could be two possibilities: either it has an AGN mid-IR excess, as source #260169, or the real redshift is lower than our best

value quoted here ($z_{phot} = 5.00$). Given the current information, we cannot decide among these two possibilities, but we decided to keep this source in our $z \ge 5$ candidate list, as it complies with all our other selection criteria. Neither this source nor any other of our $z \ge 5$ candidates are detected in the latest SCUBA2 maps of the COSMOS field (C. Chen et al. 2015, in preparation).
<figure><img src="image_14.png" /><figcaption>**Figure 13.** GSMF high-mass end at $5.0 \le z < 6.0$ (left) and $6.0 \le z < 7.0$ (right). For the GSMF at $5.0 \le z < 6.0$, we considered the combination of all UltraVISTA sources in the DR1 release and our new sources with $[4.5] < z < 23$, $K_s^{auto} > 24$ over the UltraVISTA ultradeep stripes. Symbols are the same as in Figure 5. The GSMF constraints at $6.0 \le z < 7.0$ are based on our single $z > 6$ candidate (as no $z \ge 6$ galaxies were identified in the UltraVISTA DR1) and upper limits obtained from the absence of very massive galaxies in the UltraVISTA ultradeep area. Note that with our sample selection cut at $[4.5] < 23$, we only aim at constraining the highest-mass end of the GSMF at these high redshifts.</figcaption></figure>

## 5.2. The GSMF at $z \ge 5$

### 5.2.1. Constraints to the GSMF High-mass End at $5 \le z \le 7$

Although here we only find a small number of galaxies at $5 \le z < 7$, it is still important to understand their contribution to the GSMF. Even the fact that we find only one galaxy at $z > 6$ provides important constraints at high $z$, given our large surveyed area. Note that we only aim at constraining the highest-mass end of the GSMF at these high redshifts with our sample selection cut at $[4.5] < 23$. In this sense, our study is complementary to those selecting, e.g., Lyman break galaxies at these redshifts, which sample the intermediate-mass regime. Our possibility of investigating the GSMF highest-mass end at $z_{phot} \ge 5$ is unique, thanks to the combination of the COSMOS large area and depth of the UltraVISTA survey.

Figure 13 shows our GSMF computed with the $V/V_{max}$ method at redshifts $5 \le z_{phot} < 6$ and $6 \le z_{phot} < 7$ (see also Table 4). For the GSMF at $5 \le z < 6$, we considered a combination of the UltraVISTA galaxies with $K_s^{auto} < 24$ (comprising 70 galaxies from the DR1 over $\sim 1.5$ deg$^2$) and the new galaxies with $[4.5] < z < 23$, $K_s^{auto} > 24$ studied here. At $6 \le z_{phot} < 7$, instead, we only considered our single $z > 6$ galaxy in our new sample, as no robust galaxy candidate with $K_s^{auto} < 24$ has been identified in the UltraVISTA DR1 at these high redshifts. The upper limits have been computed considering the fact that we have no galaxies with $M > 2 \times 10^{11} M_\odot$ at $z \ge 6$ in the UltraVISTA ultradeep area.

From our GSMF determination at $z_{phot} \ge 5$, it may seem apparent that the highest-mass end at $M_{st} > 1.6-2.0 \times 10^{11} M_\odot$ had only a modest evolution from $z \sim 5.5$ to $z \sim 3.5$. However, note that the number density of such massive galaxies still rose by a factor of $\sim 4$ in a time period of only $\sim 0.8$ Gyr.

Instead, the number density of such very massive galaxies drops sharply at $z > 6$. We find almost 1 dex difference in the

<table><caption>**Table 4**
GSMF Constraints at $5 \le z \le 7$</caption><thead><tr><th rowspan="2">log<sub>10</sub>(M/M<sub>☉</sub>)</th><th colspan="2">log<sub>10</sub>(Φ<sub>M</sub>/Mpc<sup>-3</sup> dex<sup>-1</sup>)</th></tr><tr><th>5.0 &le; z &le; 6.0</th><th>6.0 &le; z &le; 7.0</th></tr></thead><tbody><tr><td>10.40</td><td>&gt; -5.19</td><td>...</td></tr><tr><td>10.60</td><td>&gt; -5.44</td><td>...</td></tr><tr><td>10.80</td><td>&gt; -5.33</td><td>...</td></tr><tr><td>11.00</td><td>&gt; -5.17</td><td>...</td></tr><tr><td>11.20</td><td>&gt; -5.41</td><td>&gt; -6.01</td></tr><tr><td>11.40</td><td>-5.10<sup>+0.25</sup><sub>-0.39</sub></td><td>&lt; -6.14</td></tr><tr><td>11.60</td><td>-5.36<sup>+0.27</sup><sub>-0.42</sub></td><td>&lt; -6.14</td></tr></tbody></table>

number density of very massive galaxies between $z \sim 6.5$ and $z \sim 5.5$. Taking into account that the cosmic time elapsed between these two redshifts is of only 0.2 Gyr, this result is striking, as we could be pinpointing quite precisely in cosmic time the moment in which the first significant population of very massive galaxies appear. We discuss this further in Section 5.3.

### 5.2.2. Comparison with Other Recent Studies

In Figure 14, we show our GSMF constraints in the context of other recent results based on the study of near-/mid-IR-selected galaxies. Some of these GSMF studies are based on data sets from the CANDELS survey (Duncan et al. 2014; Grazian et al. 2015). These authors used the ultradeep ($H \sim 27-28$ mag) CANDELS data to select galaxies over an area that is at least nine times smaller than that covered by the UltraVISTA ultradeep survey. The other results considered for comparison here are those obtained by Stefanon et al. (2015), which are based on the IRAC S-COSMOS data and the UltraVISTA DR1 release.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** GSMF constraints provided by our sample at $5.0 \le z < 7.0$, in the context of the results obtained in the CANDELS survey (Duncan et al. 2014; Grazian et al. 2015), and the upper limits obtained by Stefanon et al. (2015), at similar redshifts. As in Figure 5, the vertical dotted line indicates the stellar mass completeness limit imposed by our IRAC magnitude cut. Stellar masses correspond to a Salpeter IMF over (0.1–100) $M_{\odot}$ in all cases.</figcaption></figure>

From the results in CANDELS, we see that the GSMF covers a stellar mass regime that is complementary to that studied here: CANDELS probes mainly the intermediate-mass regime ($M_{\text{st}} \sim 10^{10}-10^{11} M_{\odot}$) at $z > 4$, while here we constrain the GSMF highest-mass end. Indeed, note that the CANDELS sample virtually contains no galaxy with $M_{\text{st}} \gtrsim 10^{11} M_{\odot}$ at $z > 5.5$. This is very likely an effect of the small volume sampled by CANDELS. In the UltraVISTA ultradeep stripes, we have a total of 22 galaxies with $M_{\text{st}} \gtrsim 10^{11} M_{\odot}$ at $z > 5.5$ (including all DR1 sources with $K_s^{\text{auto}} < 24$ and our new $K_s^{\text{auto}} > 24$ sources). If the distribution of these galaxies were perfectly homogeneous (i.e., if there were no cosmic variance), then at most one or two galaxies would be expected in the

CANDELS field area analyzed by Duncan et al. (2014) and Grazian et al. (2015).

In Figure 14, we also show the GSMF values and upper limits derived by Stefanon et al. (2015). These authors have used the same S-COSMOS IRAC sample that we used here, but in combination with the UltraVISTA DR1 data release. Therefore, in this work we have additional information to that presented by Stefanon et al., which allows us to make a better estimate of the GSMF highest-mass end at $z \ge 5$. We see that our GSMF determinations are consistent with the values derived by Stefanon et al. at $5 \le z \le 6$ within the error bars. At $6 \le z \le 7$, we place more stringent upper limits to the number density of galaxies with $M_{\text{st}} \gtrsim 2 \times 10^{11} M_{\odot}$, which are about 0.3–0.5 dex lower than their previous determination (see discussion in Section 5.3).

### 5.3. Discussion: Galaxies with $M_{\text{st}} > 10^{11} M_{\odot}$ at $z \sim 6$?

In this paper we find only one massive galaxy candidate at $z \ge 6$ over the UltraVISTA ultradeep stripe area (~0.8 deg²). The best-fit SED indicates that this is a 0.1 Gyr old galaxy at $z_{\text{phot}} = 6.04$ (the age of the universe is ~0.9 Gyr at $z=6$), with extinction $A_V = 0.30$. This source is not detected at 24 $\mu$m, as expected for such a distant source (unless it were an AGN). The derived stellar mass is $M_{\text{st}} \approx 1.8 \times 10^{11} M_{\odot}$. This single galaxy implies a minimum number density of ~1.3 × 10⁻⁷ Mpc⁻³ for $M_{\text{st}} > 10^{11} M_{\odot}$ at $z \sim 6$. Note that the stellar mass value is below the completeness limit imposed by our [4.5] < 23 cut, so we should consider this number density as a lower limit.

According to the Sheth & Tormen (1999) formalism, the expected number density of dark matter halos with $M > 10^{12.5} M_{\odot}$ is ~10⁻⁷ Mpc⁻³ at $z=6$, which is similar to the implied number density of $M_{\text{st}} > 10^{11} M_{\odot}$ galaxies that we obtain here. However, as our figure is a lower limit, we could plausibly expect some conflict between the observed number density of massive galaxies and the number density of dark matter halos that can host them. Of course, identifying galaxies with halos assumes that we know the baryon conversion fraction at high redshifts, which is in fact not known. And it is unclear whether massive halos at those redshifts would host massive galaxies or rather multiple lower-mass galaxies (both cases are possible depending on the model parameters).

In fact, current semianalytic models basically do not predict any galaxy with $M \gtrsim 5 \times 10^{10} M_{\odot}$ at $z=6$ (Y. Lu 2015, private communication). However, there is no reason to conclude that the different candidates found so far in CANDELS and in this work are not real and merely the consequence of errors in the photometric redshifts and stellar mass determinations.

Our results also set strong upper limits to the number density of galaxies with $M_{\text{st}} \gtrsim 2 \times 10^{11} M_{\odot}$ at $z > 6$, by finding no robust candidate with such high stellar masses at these redshifts. Remarkably, the number density is significantly higher at $5 \le z < 6$, which strongly suggests that the most massive galaxies in the universe only become that massive at $z < 6$.

As an alternative, one could wonder whether there could still be very massive galaxies with $M_{\text{st}} \gtrsim 2 \times 10^{11} M_{\odot}$ at $z \sim 6$, among the ~1% of [4.5] < 23 sources that remain unidentified in the UltraVISTA DR2 survey, and/or among the sources with unknown redshifts in our current sample. A few of the latter have similar near-/mid-IR colors to our $z \sim 6$ candidate
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Redshift evolution of the cosmic stellar mass density, given by our new GSMF determinations and a compilation of recent results from the literature. At $3 \le z < 5$, our stellar mass density values have been computed by integrating the best STY Schechter functions obtained here at stellar masses $log_{10}(M) > 10^8 M_{\odot}$ (after correction for Eddington bias). At $5 \le z < 6$, we only sum up the contributions of our galaxies in the stellar mass bins above the IRAC completeness limit (see Figure 13), so the stellar mass density constitutes a lower limit. The open circle indicates our best estimate of the total stellar mass density at $5 \le z < 6$, taking into account the results of the CANDELS survey and our own, in complementary stellar mass regimes. We include also different model predictions: solid and dotted lines correspond to the mean and 99% confidence interval in Lu et al. (2014); the dashed and dot-dashed lines correspond to mean values in Somerville et al. (2012) and Croton et al. (2006), respectively. All stellar mass density values in this plot refer to a Salpeter IMF over (0.1–100) $M_{\odot}$.</figcaption></figure>

(Figure 3), but unfortunately these colors are not conclusive for $z \sim 6$ objects, as some sources at $z < 3$ and a few at $3 \le z < 5$ display similar colors. Although in principle we cannot completely exclude the possibility that very massive galaxies exist at $z > 6$, we note that this could only happen if these sources had significant dust extinction. With no dust, the 4000 Å break of a maximal age galaxy at $z = 6$ would produce $K_s - [4.5] \lesssim 1$, so all the [4.5] < 23 with these characteristics should now be identified in the UltraVISTA survey. Hence, we conclude that, unless sources with significant amounts of dust extinction exist at $z > 6$, the appearance of the most massive galaxies only happens at $5 < z < 6$. Future studies with ALMA and the James Webb Space Telescope will allow us to confirm whether such very massive galaxies really exist at $z > 6$, but candidates will need to be found in advance with large-area ultradeep near-IR surveys.

## 6. THE EVOLUTION OF THE COSMIC STELLAR MASS DENSITY

We used our new GSMF determinations to set updated constraints on the cosmic stellar mass density up to $z = 6$. Figure 15 shows the evolution of the cosmic stellar mass density given by our new results and a compilation of recent results from the literature. Our values at $3 \le z < 5$ have been obtained by integrating the Schechter function after correction for Eddington bias, with the parameter values given in Table 2,

so they can be considered total values, and we also show the values obtained by integrating the GSMF only above the stellar mass completeness limits. At $5 \le z \le 6$, instead, we do not probe a wide enough stellar mass range to attempt a good Schechter function determination, so we only sum up the individual contributions of our galaxies (see Figure 13). Therefore, at these high redshifts, our sample provides only a lower limit of the cosmic stellar mass density.

As the stellar mass regime that we probe at $5 \le z \le 6$ is complementary to that probed by CANDELS, we can estimate the total stellar mass density estimates at these redshifts by summing up the contributions obtained here and the average of the CANDELS determinations (Duncan et al. 2014; Grazian et al. 2015). In this case, we only considered the galaxies in our sample above our stellar mass completeness limit (imposed by the IRAC completeness) and the CANDELS results at lower stellar masses. The resulting total stellar mass density estimate is indicated with a large open circle in Figure 15.

Our new stellar mass density determinations at $3 \le z < 5$ are in agreement with most recent results from the literature, within the error bars. However, at $z \sim 3-4$ we cannot reproduce the high stellar mass density values found by Mortlock et al. (2011) and Santini et al. (2012), even within the errors. These two studies have been based on very small areas of the sky, so cosmic variance is the most likely reason for this discrepancy. On the other hand, at $z \sim 4-5$, our new stellar mass density determination is significantly above the value obtained by Caputi et al. (2011). In this case, this is mainly due to the incompleteness in the near-IR identification of high-z IRAC sources (Caputi et al. 2011 used the UDS data release 5). This effect clearly illustrates the importance of the ultradeep near-IR data in identifying massive galaxies and determining the total cosmic stellar mass density, particularly at $4 \lesssim z \lesssim 5$.

Figure 13 also includes different model predictions of the cosmic stellar mass density evolution (Croton et al. 2006; Somerville et al. 2012; Lu et al. 2014). These models, and particularly the latest ones, reproduce reasonably well the stellar mass density values up to $z \sim 7$, as they are dominated by the contribution of intermediate and low stellar mass galaxies. However, they cannot reproduce the presence of massive galaxies with $M_{\text{st}} \gtrsim 5 \times 10^{10} M_{\odot}$ up to $z \sim 6$ (Y. Lu 2015, private communication).

## 7. SUMMARY AND CONCLUSIONS

In this paper we have studied a sample of 574 Spitzer IRAC bright ([4.5] < 23), near-IR faint ($K_s^{\text{auto}} > 24$) galaxies over ~0.8 deg$^2$ of the UltraVISTA ultradeep COSMOS field. This is the first time that the study of such galaxies can be conducted over a large area of the sky, and it is becoming possible thanks to the unique combination of area and depth that is being achieved by the UltraVISTA survey. Our galaxy sample constitutes a small (<1%) fraction of the overall [4.5] < 23 galaxy population. However, in this paper we have shown that these galaxies, previously unidentified in statistically large samples, provide a significant contribution to the total population of massive galaxies at high redshifts.

Indeed, from the SED analysis of our [4.5] < 23, $K_s^{\text{auto}} > 24$ galaxies, we have determined that their redshift distribution peaks at redshift $z \sim 2.5-3.0$, and ~32% of them lie at $3 \le z \le 6.04$. We found that colors $H - [4.5] > 4$ almost
exclusively select galaxies at $3 \le z < 5$, consistently with the result of theometry ~0.5 mag deeper than the current DR2 release, will Caputi et al. (2012).

We analyzed the contribution of our $z \ge 3$ galaxies to the GSMF, finding of massive galaxies at high $z$.

mass end at high redshifts. We found that our

galaxies make a very minor contribution to the GSMF based on data products from observations made with ESO 3 $\le z < 4$, previously determined with $K_s^{\text{auto}} < 24$ galaxies. Instead, telescopes at the La Silla Paranal Observatory under ESO have a more significant role within the $4 \le z < 5$ GSMF, accounting for program ID 179.A-2005 and on data products produced by 50% of the galaxies with stellar masses $> 10^{10} M_{\odot}$. We concluded that the Hubble and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. Based on observations considering

these IRAC bright, near-IR faint galaxies at $4 \le z < 5$ is excluded out with the *Spitzer Space Telescope*, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA; the NASA/ESA Hubble Space Telescope, obtained and archived at the Space Telescope Science Institute; and the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

previously found, after correcting for the effect of Eddington bias. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. O.I. acknowledges funding from the French Agence Nationale de la Recherche (ANR) for the project SAGACE. B.M.-J. and J.F. acknowledge support from the factor of about four between these redshifts. The number density at $z \sim 0.8$ is funded by the DNR. We thank Kenneth Duncan, Yu Lu, and increased by another factor of 10 later by $z \sim 1.5-2.0$, in the 2 billion years encompassing the star formation activity peak of the universe (Hopkins & Santini for providing us data in electronic format; Beacom 2006; Behroozi et al. 2013). Finally, that number density increased by another factor of about four until reaching the present value. Thus, in conclusion, the most massive galaxies were formed quite effectively after the

first billion years: almost as effectively as during the peak activity epoch, and much more effectively than over the past 10 billion years of cosmic time.

## REFERENCES

The presence of very massive ($M_{\odot} \gtrsim 10^{11} M_{\odot}$) galaxies

in our sample at $5 \le z < 6$ and the virtual absence of massive galaxies at $z \le 0.8$ by Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540

Balay, J. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621

Bell, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

end, which suggests that the appearance of such massive galaxies depends on, & Arnouts, S. 1996, A&AS, 117, 393

place in the few hundred million years of elapsed time between 1980 and 2000, by Bielby, R., Hudelot, P., McCracken, H. J., et al. 2012, A&A, 545, A23

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

surveys like CANDELS (Grogin et al. 2011; Koekemoer et al. 2011), by D. Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

galaxy surveys are necessary to sample this very massive galaxy at $z = 1.4$ by I. 2014, IJMPD, 23, 1430015

population at high $z$. The combination of area and depth of UltraVISTA by I. 2014, MNRAS, 413, 162

currently unique for this purpose. by Caputi, K. I., Dunlop, J. S., McLure, R. J., & Roche, N. D. 2004, MNRAS, 353, 30

The only alternative to this conclusion is that, among the [4.5] $z \le 3$ galaxies that remain unidentified and/or those that have no redshift, 607

determination in our current sample, there is a population of very massive galaxies at $z \le 0.8$ by Arnouts, S., Cristiani, S., Moscardini, L., et al. 2012, ApJL, 750, L20

galaxies at $z > 6$ that are significantly dust obscured. One possible source of MNRAS, 366, 609

candidate for such galaxies was discussed by Caputi et al. (2012) in one of & Kaiser, N. 1989, MNRAS, 237, 1127

the

CANDELS fields, but the level of dust obscuration ($A_V = 0.90$ mag) at $z \le 1.4$ by I., Bolzonella, M., Coupon, J., et al. 2013, A&A, 558, A23 de Barros, S., Schaerer, D., & Stark, D. P. 2014, A&A, 563, A81

see Oesch et al. 2015). Further studies in other fields, as well as by Soifer, B. T., Desai, V., et al. 2008, ApJ, 677, 943

follow-up with JWST and ALMA, are necessary to confirm whether the sources exist at $z > 6$. In any case, a substantial fraction of the sources are in the field of view of the HST.

unidentified IRAC [4.5] $< 23$ sources are more likely dust-obscured massive galaxies at $z \le 0.8$ by I., Bolzonella, M., Coupon, J., et al. 2013, A&A, 558, A23

galaxies at $z \le 0.8$ at $z \le 4$ by Bielby, R., Hudelot, P., McCracken, H. J., et al. 2012, ApJ, 770, 155

diagram shown in Figure 3 and our GSMF determinations. Further analysis by G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

of this problem after the UltraVISTA completion, which will achieve 600

Franx, M., Labbé, I., Rudnick, G., et al. 2003, ApJL, 587, L79

Georgakakis, A., Hopkins, A. M., Afonso, J., et al. 2006, MNRAS, 367, 331

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJL, 735, L34

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142

Houck, J. R., Soifer, B. T., Weedman, D., et al. 2005, ApJL, 622, L105

Huang, J. S., Zheng, X. Z., Rigopoulou, D., et al. 2011, ApJL, 742, L13

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, A&A, 457, 841
Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Ilbert, O., Salvato, M., Le Floc'h, E., et al. 2010, ApJ, 709, 644

Kajisawa, M., Ichikawa, T., Tanaka, I., et al. 2009, ApJ, 702, 1393

Kirkpatrick, J. D., Cushing, M. C., Gelino, C. R., et al. 2011, ApJS, 197, 19

Kodama, T., Tanaka, I., Kajisawa, M., et al. 2007, MNRAS, 377, 1717

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Labbé, I., Huang, J., Franx, M., et al. 2005, ApJL, 624, L81

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, A&A, 576, A79

Lilly, S. J., Le Fèvre, O., Renzini, A., et al. 2007, ApJS, 172, 70

Lu, Y., Wechsler, R. H., Somerville, R. S., et al. 2014, ApJ, 795, 123

Madau, P. 1995, ApJ, 441, 18

Mancini, C., Matute, I., Cimatti, A., et al. 2009, A&A, 500, 705

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

Milvang-Jensen, B., Freudling, W., Zabl, J., et al. 2013, A&A, 560, A94

Mo, H. J., & White, S. D. M. 1996, MNRAS, 282, 347

Mortlock, A., Conselice, C. J., Bluck, A. F. L., et al. 2011, MNRAS, 413, 2845

Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Oesch, P., van Dokkum, P., Illingworth, G., et al. 2015, ApJL, 804, L30

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Pope, A., Bussmann, R. S., & Dey, A. 2008, ApJ, 689, 127

Pozzetti, L., Bolzonella, M., Lamareille, F., et al. 2007, A&A, 474, 443

Pozzetti, L., Bolzonella, M., Zucca, E., et al. 2010, A&A, 523, A13

Rayner, J. T., Cushing, M. C., & Vacca, W. D. 2009, ApJS, 185, 289

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Rodighiero, G., Cimatti, A., Franceschini, A., et al. 2007, A&A, 470, 21

Salpeter, E. E. 1955, ApJ, 121, 161

Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Santini, P., Fontana, A., Grazian, A., et al. 2012, A&A, 538, 33

Saracco, P., Longhetti, M., Giallongo, E., et al. 2004, A&A, 420, 125

Saracco, P., Longhetti, M., Severgnini, P., et al. 2005, MNRAS, 357, L40

Schäerer, D., & de Barros, S. 2009, A&A, 502, 423

Schechter, P. 1976, ApJ, 203, 297

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Somerville, R. S., Gilmore, R. C., Primack, J. R., & Domínguez, A. 2012, MNRAS, 423, 1992

Stefanon, M., Marchesini, D., Muzzin, A., et al. 2015, ApJ, 803, 11

Sutherland, W., Emerson, J., Dalton, G., et al. 2015, A&A, 575, A25

Taniguchi, Y., Scoville, N., Murayama, T., et al. 2007, ApJS, 172, 9

Teerikorpi, P. 2004, A&A, 424, 73

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Wiklind, T., Dickinson, M., Ferguson, H. C., et al. 2008, ApJ, 676, 781

Zackrisson, E., Bergvall, N., & Leitet, E. 2008, ApJL, 676, L9

Zucca, E., Zamorani, G., Vettolani, G., et al. 1997, A&A, 326, 477