UltraVISTA SMF at 4<z<7

Stellar mass functions of galaxies at from an IRAC-selected Sample in COSMOS/UltraVISTA: Limits on the abundance of very massive galaxies

Mauro Stefanon11affiliation: Physics and Astronomy Department, University of Missouri, Columbia, MO 65211, USA , Danilo Marchesini22affiliation: Physics and Astronomy Department, Tufts University, Robinson Hall, Room 257, Medford, MA, 02155, USA , Adam Muzzin33affiliation: Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands , Gabriel Brammer44affiliation: Space Telescope Science Institute, Baltimore, MD 21218, USA , James S. Dunlop55affiliation: SUPA - Scottish Universities Physics Alliance - Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK , Marijin Franx33affiliation: Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands , Johan P. U. Fynbo66affiliation: Dark Cosmology Centre, Niels Bohr Institute, Copenhagen University, Juliane Maries Vej 30, DK-2100 Copenhagen O, Denmark , Ivo Labbé33affiliation: Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands , Bo Milvang-Jensen66affiliation: Dark Cosmology Centre, Niels Bohr Institute, Copenhagen University, Juliane Maries Vej 30, DK-2100 Copenhagen O, Denmark , Pieter G. van Dokkum77affiliation: Department of Astronomy, Yale University, New Haven, CT 06511, USA Email: stefanonm@missouri.edu
Abstract

We build a Spitzer IRAC complete catalog of objects, obtained by complementing the -band selected UltraVISTA catalog with objects detected in IRAC only. With the aim of identifying massive (i.e., ) galaxies at , we consider the systematic effects on the measured photometric redshifts from the introduction of an old and dusty SED template and from the introduction of a bayesian prior taking into account the brightness of the objects, as well as the systematic effects from different star formation histories (SFHs) and from nebular emission lines in the recovery of stellar population parameters. We show that our results are most affected by the bayesian luminosity prior, while nebular emission lines and SFHs only introduce a small dispersion in the measurements. Specifically, the number of galaxies ranges from 52 to 382 depending on the adopted configuration. Using these results we investigate, for the first time, the evolution of the massive end of the stellar mass functions (SMFs) at . Given the rarity of very massive galaxies in the early universe, major contributions to the total error budget come from cosmic variance and poisson noise. The SMF obtained without the introduction of the bayesian luminosity prior does not show any evolution from to , implying that massive galaxies could already be present when the Universe was  Gyr old. However, the introduction of the bayesian luminosity prior reduces the number of galaxies with best fit masses by 83%, implying a rapid growth of very massive galaxies in the first 1.5 Gyr of cosmic history. From the stellar-mass complete sample, we identify one candidate of a very massive (), quiescent galaxy at , with MIPS m detection suggesting the presence of a powerful obscured AGN. Finally, we show that the number of massive galaxies at measured in this work matches the number of massive galaxies at predicted by current models of galaxy formation.

Subject headings:
galaxies: high-redshift, galaxies: evolution, galaxies: fundamental parameters, galaxies: luminosity function, mass function

1. Introduction

Understanding the formation and evolution of galaxies has remained a central topic in astrophysics for many decades. Since the pioneering works (e.g., Hoyle 1953), this field has seen a substantial step forward, especially in the last decade, thanks to the developments in data acquisition and analysis techniques. Specifically, this has allowed the community to study a large number of galaxies up to , with the frontier gradually shifting to (see e.g., Fontana et al. 2010; Pentericci et al. 2011; McLure et al. 2011; Yan et al. 2012; Caruana et al. 2012; Ono et al. 2012; Bouwens et al. 2013; Curtis-Lake et al. 2013; Labbé et al. 2013a; Finkelstein et al. 2013; Coe et al. 2013; Ellis et al. 2013; Bradley et al. 2014; Zitrin et al. 2014; Zheng et al. 2014).

Different measurements of galaxy masses provide information on different physical aspects. Dynamical masses measure the associated gravitational potential allowing for a closer insight into the dark matter halo properties and their evolution across cosmic time (see e.g., Gerhard et al. 2001; Thomas et al. 2011; Cappellari et al. 2013; Beifiori et al. 2014). However, recovering dynamical masses requires spectroscopy, more difficult to obtain for galaxies at high redshift, given their lower brightness. Stellar masses, on the other hand, reflect the information on the luminous matter content. They do not require spectroscopy, as they can also be recovered from modelling of the multi-wavelength photometry, and are thus available for a much larger sample of high redshift galaxies. The observed profile of the stellar mass function (SMF) across cosmic time is the result not only of the effects of the gravitational potential governing the assembly of the dark matter halos, but also of the physical processes which govern the baryonic matter, such as the formation of new stars or the quenching of star formation.

The SMF is then one of the statistical tools that are most commonly used to trace the evolution of the galaxy populations across cosmic time and is one of the main observables whose reproduction is a necessary step for the validation of galaxy formation models. Several measurements of the SMF exist up to (see e.g., Muzzin et al. 2013b and Ilbert et al. 2013 and references therein). The emerging picture is that some of the very massive galaxies were already in place at and their associated number density has rapidly evolved between and . From a complementary perspective, so-called archeological studies have shown that local most massive galaxies formed most of their stars during a short burst at (e.g., Thomas et al. 2010). This picture is however complicated by the fact that the evolution of dark matter halos follows a rigid hierarchical structure, with more massive haloes forming at later times only (see e.g., Springel et al. 2005; Baugh 2006 and references therein). It is therefore important to be able to track the formation and evolution of the most massive galaxies at earlier cosmic times.

The SMF measurements at are typically constructed from samples selected in the near-infrared (NIR) band. The advantage of the single-band photometric selection is that it allows for the assembly of stellar mass complete samples of galaxies. However, the rest-frame optical Balmer/4000Å break begins to enter the K band at and it is completely included by . Because of this effect, galaxies at with a given stellar mass will then be characterized by even fainter K-band fluxes than galaxies at with the same stellar mass and mass-to-light ratio (M/L). The dimmed K-band fluxes drop below the detection threshold, preventing these objects from being detected. Conversely, detected objects will be characterized by K-band fluxes equal or higher than the detection threshold, and thus correspond to higher values of stellar mass (for the same M/L value) required for the selection of stellar-mass-complete samples of galaxies.

Furthermore, the K-band selection becomes more and more sensitive to the extinction by dust with increasing redshift, as the rest-frame optical and ultra-violet (UV) enter the K-band. Given the current photometric depth of K-band selected samples, these factors thus limit our ability to perform statistical studies at higher redshifts. For this reason, SMF at are generally computed using samples of galaxies selected in the rest-frame UV via color selection criteria, through the so-called dropout technique (exceptions exist, such as Caputi et al. 2011 who measure the SMF up to using photometric redshifts from an IRAC m-selected sample in the UDS field). However, by construction, dropout selections are biased against evolved (i.e., quiescent) galaxies and/or dust-extincted systems as they preferentially pick those galaxies with brighter rest-frame UV luminosities, indicative of those systems with a recent burst or still ongoing star-formation and low dust content. This is even more important for massive galaxies, as recent studies have shown an increase of dust extinction with redshift (see e.g., Whitaker et al. 2010; Marchesini et al. 2010, 2014). Specifically, this prevents from searching for massive and massive and quiescent galaxies at (see e.g., Mobasher et al. 2005; McLure et al. 2006; Dunlop et al. 2007; Lundgren et al. 2014), potentially introducing a bias in our knowledge of the high-mass end population of galaxies and, in general, in the first 2 Gyr of cosmic history.

On the other side, as we will show in Sect. 3.8, the stellar mass completeness limit in bands at wavelength m (e.g., IRAC m and m bands) is roughly constant for , and, with current imaging depth, corresponds to (see also e.g., Fontana et al. 2006; Ilbert et al. 2010; Caputi et al. 2011). However, the larger point-spread functions (PSF) of IRAC bands, compared to those in the optical/NIR imaging, enhance the problem of source blending in the measurements of the fluxes; these effects become even more important when IRAC bands are considered for the detection of sources. The most commonly adopted solution consists in performing the photometry using positional and morphological information from higher resolution imaging, usually in bands different from that of interest, and under the assumption that the morphological properties of the objects in the lower resolution band do not significantly differ from those in the higher resolution one (see e.g., Fernández-Soto et al. 1999; Labbé et al. 2005; Laidler et al. 2007). If this approach solves the problem of contamination in the flux measurement, it can not be directly adopted for the detection of sources in IRAC bands, as it relies on an already existing catalog.

The aim of this work is to search for a population of massive () galaxies at and to study their evolution through the analysis of the SMF in three redshift bins: , and . To this aim, we used an IRAC m-complete sample of galaxies assembled complementing the UltraVISTA DR1 -band selected catalog by Muzzin et al. (2013a) with detections on the residual maps of IRAC m and m bands. Indeed, the IRAC flux measurements included in the UltraVISTA -band selected catalog were obtained using the technique described above, which relies on position and morphological information from the K-band map. These measurements are less sensitive to the contamination from neighbours than those performed directly on the IRAC maps. Performing the detections on the IRAC residual maps allows for the inclusion in the catalog of all those sources with K-band fluxes below the K-band detection threshold. The advantage of an IRAC-based sample is that it is possible to detect galaxies at with stellar masses lower than with the current UltraVISTA -band based sample, for the same M/L value, since galaxies in the IRAC m and m do not suffer from the dimming due to the Balmer/4000Å break, which instead affects the K-band data. This approach then allowed us to exploit the deeper completeness in stellar mass associated to the IRAC bands, while solving the issues in the detection and flux measurements from source blending in purely IRAC-detected catalogs.

The still large uncertainties in the knowledge of the SEDs of galaxies in this range of redshift can introduce systematic effects which potentially bias the measurements of photometric redshifts and stellar population parameters. We approach this problem by exploring different configurations for the measurements of photometric redshifts and stellar population parameters and treat the results as systematic effects. In order to make the presentation more organized, we assume the sample obtained from the most relaxed configuration to be the default sample and consider the statistical differences arising from the other configurations as systematic effects. One of the main results of this approach, presented in Section 4.1, is a sample of galaxies at and with irrespective of the configuration adopted for the measurement of photometric redshifts and stellar population parameters.

This paper is organized as follows. In Section 2 we describe how we build the IRAC-complete sample, complementing the UltraVISTA DR1 -band selected catalog by Muzzin et al. (2013a) with detections on the residual maps of IRAC m and m bands. In Section 3 we present the step-by-step process we adopted to obtain a clean sample of massive () galaxies at , removing galaxies potentially contaminated by AGN and potential lower redshift interlopers, and addressing different factors which can affect the measurements of photometric redshifts and stellar population parameters. In Section 4 we present the robust sample of galaxies with at irrespectively of the systematic effects discussed in Section 3. Specifically, in Section 4.3 we present the SMF measurements in the redshift range and show the systematic effects arising from the factors presented in Section 3 on the SMF measurements. Our results are discussed in Section 5, while we summarize and conclude in Section 6.

Throughout this work we adopt a concordance cosmology, with  km/s/Mpc, and . Unless otherwise specified, all magnitudes are referred to the AB system, while stellar masses were computed using the Chabrier (2003) initial mass function (IMF).

2. The photometric sample

For this work, we complemented the public catalog from Muzzin et al. (2013a) with detections on the IRAC residual maps from S-COSMOS (Sanders et al., 2007; Ilbert et al., 2010). In the following paragraphs we present the details of the procedure we followed to construct the IRAC-complete sample from the -band selected catalog of Muzzin et al. (2013a). The -band selected catalog of Muzzin et al. (2013a) is based on the first data release (DR1) of the UltraVISTA survey (McCracken et al., 2012) and delivers 30-bands flux information for 262615 objects. The DR1 of the UltraVISTA survey is characterized by deep (90% point source completeness AB) imaging in four broad-band NIR filters () as well as one narrow-band filter centered on H at (NB118). Images in each band cover an area of degrees centered on the COSMOS field (Scoville et al. 2007). This field is characterized by extensive deep multi-wavelength coverage, ranging from the X-rays (Hasinger et al. 2007; Elvis et al. 2009) to the radio (Schinnerer et al. 2007). Specifically, Spitzer IRAC imaging data were collected in the framework of the S-COSMOS project (Sanders et al., 2007; Ilbert et al., 2010).

The flux measurements in the S-COSMOS IRAC for the UltraVISTA -band selected catalog of Muzzin et al. (2013a) were obtained using a source-fitting code specifically developed to recover fluxes in heavily confused images (Labbé et al., 2005). The code is based on the assumption that the morphology of each object in the higher-resolution band ( for our case) does not differ too much from the intrinsic morphology in the lower resolution image (i.e., IRAC). It is then possible to use the position and brightness profile of the source on the high resolution image as a prior for the same object in the lower resolution image. The brightness profile of the source in the high-resolution image is convolved with the kernel required to match the low-resolution image PSF. The result is then a template of the object in the low-resolution image modulo its total flux which is obtained via best fit. The actual flux measurements for each object is however performed in apertures on a per-object basis, after all the neighbouring objects have been removed using the information from the fitting procedure. One of the diagnostic outputs from this procedure is the image resulting from subtracting all the fitted sources from the input science frame, i.e., a residual image. Since no new detection is performed during the source-fitting measurements, the residual image will contain those sources not detected in the high-resolution image. An example of this process is shown in Figure 1. We complemented the UltraVISTA DR1 -band selected photometric catalog of Muzzin et al. (2013a) with objects detected on the IRAC-band residual images. These objects are too faint to be robustly detected in the DR1 K-band image and therefore were not included in the K-selected UltraVISTA DR1 catalogIRAC bands, whose -band counterparts are too faint to be included in the original -band selected UltraVISTA catalog. This was achieved by performing an additional detection on the residual images resulting from the template fitting procedure as described below.

Figure 1.— Example of the steps involved in the detection of sources in the IRAC bands. The left and right columns refer to the same procedure in the IRAC m and m, respectively. In each column, the top panel shows the cutout of the UltraVISTA -band frame centered on the position of an IRAC source randomly picked among those originally not included in the -band catalog of Muzzin et al. (2013a). Its position is also marked across all panels by the green circle to facilitate its identification. The second panel reproduces the brightness profiles of all the objects, for the same region of sky, from the -band catalog, after convolving them with the kernel required to match the PSF of the corresponding IRAC channel and after applying to each source the flux scale factor from the best fit process. This image, then, constitutes the best-fit model of the science frame based on the information available from the band. Sources not included in the -based catalog will not be reproduced in the model image. The third panel panel shows the original IRAC science frame. The panel at the bottom presents the residual image obtained subtracting the model from the science frame. The additional IRAC source is still present, cleaned from neighbors. We adopted the residual images as input for the detection with SExtractor.

At first, we run SExtractor (Bertin & Arnouts 1996) on S-COSMOS IRAC m and m bands residual images independently. Only those objects with a matching position within a radius of 5 pixels (equivalent to ), and with a SNR in the two IRAC bands were kept. This approach limits the number of spurious sources, although it introduces a completeness in stellar mass shifted towards higher stellar masses, given the shallower depth of the IRAC m map compared to the IRAC m one. Objects within from the border of the UltraVISTA -band map were purged to remove detection in those regions with lower signal-to-noise ratio. In total, 408 new sources over an area of square degree were added to the existing UltraVISTA catalog, representing an increase up to 39% (depending on the configuration adopted for the recovery of photometric redshifts - see Section 3) to the number of galaxies at in the sample.

Spectral energy distributions (SEDs) were then built using matched aperture photometry from CFHTLS (Cuillandre et al., 2012), Subaru (Taniguchi et al., 2007), UltraVISTA (McCracken et al., 2012) and S-COSMOS (Sanders et al., 2007). Photometry was performed with SExtractor in dual mode, with IRAC m residual image as the detection image, to take advantage of its slightly better image quality compared to IRAC m map. We adopted apertures of diameter for CFHTLS, Subaru and UltraVISTA frames, for IRAC bands and for MIPS 24m as a compromise between the highest signal-to-noise ratio measurements and the possible loss of flux due to the uncertainty in the position of the source on the IRAC frame. The sizes of the aperture adopted for the photometry were chosen taking into consideration potential contamination from neighboring objects. The FWHM of the broadest PSF in ground-based data is about , so that the chosen aperture corresponds to FWHM. The fraction of pairs closer than the adopted aperture radius is about , a value which still ensures a limited contamination fraction from neighbors. Finally, objects close to bright stars and those with contamination by bright saturated nearby objects in more than three filters were excluded from the selection (flag use=1 in the UltraVISTA catalog). For IRAC and MIPS data, instead, the flux measurement is done for each object after cleaning from neighbors, minimizing potential contamination effects before the actual measurement is done.

Total fluxes were finally computed applying an aperture correction obtained from the curve of growth of a sample of bright and isolated point sources in each band.

This sample was finally cleaned from potential brown dwarfs (BD) and variable objects. Current observations of BD SEDs do not yet cover the wavelength region redder than m. The molecular bands in the photosphere of the cool stars result in peaks and plateaux which largely overlap with the NIR broad bands (see e.g., Figure 5 in Bowler et al. 2012). These features do not allow for a reliable removal of the degeneracy between BD and red/high-z galaxies. The identification of candidate brown dwarfs was therefore carried out in two steps: at first we pre-selected the sample of candidate BDs from color-color plots; the PSF of each candidate was then measured on the public HST/ACS frames (Koekemoer et al. 2007 - where ACS was not available, we considered those bands with the best compromise between resolution and S/N). We thus excluded from our sample those pre-selected objects with a full-width-half-maximum (FWHM) smaller than for HST ( for ground-based), being consistent with point sources. We further removed sources which showed clear signs of temporal variability. Through this selection we removed 54 objects from the original sample, leading to a sample of 502 galaxies at .

Figure 2.— The five panels track the effect of the selections applied to the sample and of the different configurations adopted for the measurements of photometric redshifts. From top to bottom, the panels refer to (a) the initial sample, obtained after cleaning from brown dwarfs, (b) the initial sample after eliminating AGN candidates; (c) the AGN-cleaned sample after removing those objects with an emission blueward of the Lyman limit; (d) the sample obtained after introducing an old and dusty template for the measurements of photometric redshift and (e) the sample obtained after applying the bayesian luminosity prior in the measurements of the photometric redshifts. In each panel, the filled blue circles correspond to the object properties after applying the corresponding selection, while filled grey circles represent the objects from the progenitor sample, i.e., the sample identified by filled blue circles in the previous panel. The orange curve delimits the stellar-mass complete sample (see Sect. 3.8). The percentages of objects excluded from the sample because potential AGN (18%) or because showing a detection blue-ward of the Lyman limit (8%) or after introducing the dusty template (19%) are relatively small, the percentage of objects excluded from the sample when the luminosity prior is introduced in the measurements of photometric redshifts is very high, reaching 83%.

3. Photometric redshifts and stellar population parameters

The selection of a sample of reliable (massive) galaxies from a photometric catalog requires taking into account many different aspects that can possibly taint the sample and/or introduce systematic effects in the measurements of their physical parameters (see e.g., Dahlen et al. 2010; Marchesini et al. 2010; Lee et al. 2012). Specifically, when the measurements of photometric redshifts are involved, the lack of extensive spectroscopic data on galaxies limits our knowledge on the SEDs of such objects, making the choice of suitable SED templates less straightforward, thus affecting the reliability of the measurements of photometric redshifts. A second source of uncertainty in the estimation of photometric redshifts comes from the adoption of a flux- and redshift-based bayesian prior. As originally shown in Benítez (2000), the adoption of prior information through the Bayesian formalism can drastically limit the fraction of outliers and reduce systematic biases in the measurement of photometric redshifts. Specifically, the prior adopted in this work consists in the distribution of galaxies as a function of redshift and of flux density as it was recovered from semi-analytic models (see Section 3.5). Furthermore, the presence, in a galaxy core, of an Active Galactic Nucleus (AGN) can bias the measurements of the stellar population parameters; finally, as it has recently been shown, biases in the estimation of the stellar masses can also be introduced by nebular emission lines. In the following sections we will systematically address these uncertainties. The involved steps are briefly summarized hereafter. Figures 2 through 10 track the discussed systematic effects as this sample is polished. We start from the full sample of galaxies obtained after cleaning the composite sample from point-source objects and potential brown dwarfs (the preparatory sample - Section 3.1). We then identify and remove from this sample galaxies potentially contaminated by AGN (Section 3.2). Successively, we purge those objects with a non-zero flux in those bands bluer than the Lyman limit at the redshift measured for each galaxy, as such objects are inconsistent with galaxies (Section 3.3). The sample obtained after this multi-step polishing process constitutes our default sample of galaxies. Panels a) through c) in Figure 2 give a graphical representation of these steps. Successively, we study the systematic effects on the redshift measurement of the inclusion in the set of SED templates of a maximally red SED template (Section 3.4) and of the bayesian luminosity prior (Section 3.5) - see also panels d) and e) in Figure 2, respectively. Finally we analyze the systematic effects of contamination by nebular emission lines and of different star-formation histories (SFHs) in the recovery of the stellar population parameters (Sections 3.6 and 3.7).

In total, from the combination of the above systematic effect analysis, this approach generated 32 samples of massive galaxies, each one corresponding to a different configuration: two configurations are the result of the introduction/exclusion of the old and dusty template; two configurations result from the activation or not of the bayesian prior on the flux, four configurations proceed from the different treatment of the contamination by nebular lines and two more configurations derive from the considered SFHs. Given the current uncertainties in the measurement of the physical parameters of galaxies from broad-band photometry, we would like to stress here that each one of these samples is potentially a consistent measurement of the properties of galaxies. In order to make the presentation of the results clearer, in this work we assume as a reference the sample obtained from the configuration with the smaller set of assumptions, i.e., photometric redshifts computed without the inclusion of the old and dusty template and without activating the bayesian prior on the flux, and stellar population parameters measured with the delayed exponential SFH and without applying any correction for the contamination by nebular lines. We therefore consider all the other configurations as systematic effects.

3.1. The preparatory sample

Figure 3.— Examples of two AGN candidates. The measured photometry in observer frame is represented by the filled colored squares with error bars (yellow for optical bands, orange for UltraVISTA and bands, red for IRAC to m). The best fitting SED from EAzY and from FAST are shown as the blue and pink curves respectively. The main physical parameters are also listed at the top-left corner. Top to bottom they are: the photometric redshift (), the , the and the extinction expressed in magnitudes. The inset shows a cutout of the object in the HST ACS F814W band and UltraVISTA for the object on the left and right respectively. The HST ACS F814W cutout shows a point-source morphology. Objects like the two presented here were excluded from the sample.

Photometric redshifts were initially computed for all the 502 objects in the sample presented on Sect. 2, using EAzY (Brammer et al., 2008). EAzY performs the redshift measurements by fitting the observed spectral energy distributions to linear combinations of a number of templates. The template set adopted for the preparatory sample is constituted by six templates from the PEGASE models (Fioc & Rocca-Volmerange, 1999) which also include emission lines and was presented in Muzzin et al. (2013a). Although EAzY allows for the use of a bayesian prior, we did not activate this option at this stage. The systematic effects of the introduction of the bayesian prior will be discussed in Sect. 3.5

Since EAzY does not deliver information on the stellar population parameters, the photometric redshifts were then used as an input to FAST (Kriek et al., 2009b) for the measurements of stellar masses, star-formation rates and ages. Indeed FAST can measure photometric redshifts using the minimization procedure, although it does not provide the possibility of introducing any bayesian luminosity prior, as instead is the case for EAzY. The uncertainties associated to the stellar population parameters are natively computed by FAST through Monte Carlo simulations. The errors on the stellar population parameters quoted in the figures of this paper refer to the upper and lower 68% confidence intervals produced by FAST, unless otherwise specified. The synergy between EAzY and FAST allows then for an optimal measurement of photometric redshifts and stellar population parameters.

We adopted the Bruzual & Charlot (2003) models, a Chabrier (2003) initial mass function, solar metallicity, and a delayed-exponential (-model) SFH. We also note that the set of templates adopted for the measurements of the stellar population parameters does not include any SED template with AGN emission. These templates were not included since a robust fit of such models to the observed data for galaxies would also require coverage in bands red-ward of IRAC/MIPS, unavailable for this work. The distribution of the stellar mass with redshift for this sample is presented in panel a) of Figure 2.

3.2. AGN contamination

The presence, in a galaxy core, of an active galactic nucleus (AGN) can bias the measurements of the stellar population parameters obtained thorough SED modelling when an AGN component is not properly included in the set of SED templates. However, a reliable estimation of the contribution from the AGN requires the modelling of the rest-frame infra-red region of the SED, which would provide information on the dusty torus around the central black hole. Spectroscopy would certainly be the preferred tool, as it allows to recover information on the main emission lines characteristic of AGNs, such as Ly, Ly, H, H, NV Å, SiIVÅ, CIII] Å, CIV Å, NIII] Å, OIII] Å, MgII Å and CaII Å. However, this kind of information is extremely difficult to obtain for high-z galaxies, as the expected flux is very low.

Given the above uncertainties, we opted for removing from the sample all potential Type-1 AGNs. When viewed face-on, AGNs are characterized by extremely compact or point-source morphologies, corresponding to the compact region around the central black hole. On the other side, emission from the AGN can strongly contaminate the rest-frame optical region of the SED resulting in a characteristic excess compared to the SEDs of galaxies not hosting an AGN.

In order to select all potential AGNs, we visually inspected both the image stamps and the observed SEDs of all the objects in the preparatory sample and removed those showing a point-source morphology and/or whose SED presented the characteristic excess in the observed SED (occurring in the IRAC bands, given the redshift range considered here), signature of potential contamination by emission from the AGN. The right panel of Figure 3 shows the SED and image stamp of an object randomly chosen among those whose SED presents an excess in the rest-frame optical (corresponding to the IRAC bands in the observer frame). The visual inspection on morphology was performed on the ACS F814W band, and adopting optical/NIR bands where no ACS coverage was available. In the left panel of Figure 3 we show an example of such objects, with an extremely compact morphology and a strong and broad emission possibly corresponding to the rest-frame Lyman . We further improved this selection by matching our sample to public XMM-Newton (Cappelluti et al., 2009), Chandra/ACIS catalogs (Elvis et al., 2009) and to the SDSS Quasar Catalog (Pâris et al., 2012).

With this procedure, we flagged as potential AGNs 91 sources, corresponding to 18% of the preparatory sample; of these, 12 sources were already detected in X-ray data, while 79 are candidate new AGN. The objects flagged as AGN are visible as grey filled circles in panel b) of Figure 2. As it can be seen from panel b), several objects with high stellar mass turn out to actually be potential AGNs. The high stellar mass recovered from the SED fitting for theses objects is likely the result of the fact that when the stellar population parameters are measured adopting an SED template set which does not include any AGN-specific SED, the excess in flux introduced by the existing AGN component biases towards higher values the stellar mass measurements. This, then, emphasizes the importance of this kind of selection in such works, as it can significantly bias the statistical analysis. All the 91 sources flagged as potential AGN were removed from the preparatory sample.

3.3. Galaxies at : the default sample

The neutral hydrogen clouds which constitute the inter-galactic medium absorb the light emitted by distant sources when the wavelength of the photon is shorter than that corresponding to the 912Å Lyman limit at the redshift of the emitting source measured from the rest-frame of the HI cloud. Although the amount of absorption depends on the redshift of the emitting object, for objects at , the absorption can be considered total (see e.g., Moller & Jakobsen 1990, Madau 1995). From an observational point of view, this physical effect manifests itself as an absence of flux in bands covering the wavelength region blueward of the 912Å Lyman limit at the observer-frame. Specifically, this means that any object with a non-zero flux measurements in bands covering the wavelength region Å (with the presumed redshift of the source) is instead likely to be at lower redshift. We note here that the Lyman forest for sources at still transmits of the photons (see e.g., Madau 1995), meaning that we can only consider as effective limit the Lyman 912Å limit, but not the 1216Å Lyman-. Although by the transmission has dropped to 2%, for consistency we apply the same selection criteria over the full redshift range.

Figure 4.— Examples of objects excluded from the sample of galaxies because presenting an excess in the flux measured in bands blue-ward of the observer-frame Lyman limit. The colored squares mark the flux measurements in observer frame with the associated error bars; the blue and pink curves represent the best-fit SED from EAzY and FAST respectively. The main physical properties are listed at the top-left corner (see Figure 3 for further details). The inset reproduces the B-band cutout centered on the object position. The photometry in some of the bands bluer than the Lyman limit present an excess of flux. The presence of emission associated to the object is confirmed by the cutout. Objects like the two shown here were excluded from the sample of galaxies.

As a consistency check, for each object we stacked the cutouts in those bands bluer than the observer frame Lyman limit and excluded from the sample those objects whose stacked image showed a clear excess at visual inspection. We also visually inspected the SEDs of the AGN-purged sample to identify those objects with a non-zero flux (at 1) in those bands bluer than the Lyman limit. The image stamps of the selected objects in the bands blueward of the observer-frame Lyman limit were successively visually inspected in order to disentangle whether the observed flux excess was the result of contamination from bright and/or nearby objects or it was a genuine emission from the object. In this latter case, the object was removed from the sample. Following this procedure, we removed from the AGN-cleaned sample 33 sources with non-zero flux measurements blueward of the Lyman limit; these sources are marked by grey filled circle in panel c) of Figure 2, while in Figure 4 we show the SEDs and cutouts for two objects randomly chosen among those removed from the sample. The low number of objects excluded from the sample in this step supports the robustness of our photometric redshift measurements.

In our analysis, we considered as default sample the preparatory sample cleaned from AGN and from sources with clear detection in bands bluer than the Lyman limit at the redshift of each source. This sample includes 382 galaxies, which correspond to a fraction of 75% with respect to the preparatory sample, or 92% of the AGN-cleaned sample.

3.4. SED templates: the maximally red SED

We investigated the systematic effects on the measurement of photometric redshifts of the inclusion of a maximally red template, i. e. a passively evolving 1.5 Gyr old and dusty ( mag) galaxy from Muzzin et al. (2013a). A non-negligible fraction of objects whose redshift was measured to be at with the default template set, have when the maximally red template is introduced. This is graphically shown in panel d) of Figure 2. The grey filled circles mark the sample of pure galaxies, from the previous section, while the blue circles identify the objects whose redshifts have been computed with the template set containing the maximally red template and whose stellar masses have been recomputed according to the new photometric redshifts. The fraction of objects which have after the introduction of the old and dusty template with respect to the default population is 81% (61% when considering the initial sample). The SED and of two objects randomly chosen from the objects with without the dusty template but with when the old and dusty template is introduced, are shown in Figure 5. The introduction of the dusty template generally provides better fits to the data, in agreement with what found in e.g., Muzzin et al. (2013a).

Figure 5.— Examples of photometric redshift measurements with and without including the old and dusty template into the template set. Each one of the two rows refers to a distinct object, randomly chosen among those whose photometric redshift was without the inclusion of the old and dusty template, but became after including it in the template set. For each object, the panel on the left shows the measured photometry in observer frame (colored points) together with the high-redshift and low-redshift best-fit templates (magenta and blue curves respectively). The panel on the right shown the probability distribution of the redshift () for the two cases of no old and dusty template (magenta region) and with the old and dusty template (blue region). The introduction of the dusty template favors the lower redshift solutions.

3.5. Effect of bayesian luminosity prior

The core of most photometric redshifts codes is a minimization on the observed photometry of fluxes from a set of templates. However there are degenerate cases where the minimisation process alone can lead to an incorrect value. The most typical case is the inability, usually for faint objects, to distinguish if a break observed in the photometry is the result of the Balmer/4000Å break or of the Lyman break at 912Å. In order to provide a way to remove such degeneracies, Benítez (2000) introduced the bayesian analysis to the measurements of photometric redshifts. The idea is to complement the flux measurements with information such as the distribution of apparent magnitude of objects with redshift. This technique has proven to be reliable at . However, the prior is generally built on the basis of either properties of galaxies at lower redshift or from semi-analytic models (as is the case for EAzY). This makes its robustness questionable when it comes to measure photometric redshifts at higher z. For this reason, we considered the introduction of the bayesian prior in the measurements of the photometric redshift as a systematic effect. For the subsample selected from the UltraVISTA catalog, we needed to extend up to the original EAzY K-band prior, which reaches . This was done by fitting the standard EAzY prior with a functional form , with and free parameters (Benítez, 2000; Brammer et al., 2008) and extrapolating its values to .

By construction the sample recovered from the IRAC residual images is characterized by very faint (when not absent) fluxes in the K-band, making unsuitable the adoption of the prior in the K-band. We therefore built a prior for the IRAC m band. This was achieved by fitting the same functional form as before to the distribution of redshifts in bins of apparent magnitudes for objects from the simulation data by Henriques et al. (2012), based on the semi-analytic model presented in Guo et al. (2011), and accessed through the Virgo - Millennium database (Lemson & Virgo Consortium 2006, Springel et al. 2005).

The activation of the prior results in 52 residual objects at redshift , i.e., 83% of the objects from the sample at after the introduction of the old and dusty template have a redshift when the bayesian prior is also activated (the percentage becomes 90% when comparing the 52 objects to the 502 objects in the initial sample). The introduction of the prior has then the largest systematic effect, among those considered in this work, on the sample selection. A graphical representation of this substantial selection effect is presented in panel e) of Figure 2, while the SED and of two objects randomly chosen from the objects with without applying the bayesian luminosity prior but with when the bayesian luminosity prior is introduced, are shown in Figure 6.

We would like to note here that, since the prior is based on semi-analytic models, which are still very uncertain in the redshift range considered in this work, the sample of galaxies selected using photometric redshifts obtained with the adoption of the prior should be considered as one of the sources of systematic effect.

Figure 6.— Similar to Figure 5 but for the case of excluding (magenta curves and regions) or including (blue curves and regions) the bayesian luminosity prior. The bayesian luminosity prior favours the low-redshift solution.
Figure 7.— Excess in stellar mass measured after correcting the photometry for nebular emission lines. The top row refers to photometric redshifts computed without the old and dusty template and with no bayesian luminosity prior, the central row to the inclusion of the old and dusty template but not applying the luminosity prior while the bottom row refers to photometric redshifts obtained with the inclusion of the old and dusty template and with the adoption of the bayesian luminosity prior. Points are color coded according to their redshift, as specified by the legend in the top-left panel. These three configurations correspond to panels (c), (d) and (e) of Figure 2. For each row, panels from left to right correspond to one of three methods implemented to correct the photometry from contamination by emission lines. Specifically, the panels on the left refer to the EW computed from EAzY templates; the central panels to the procedure presented in Smit et al. (2014), while the panels on the right refer to the exclusion of the fluxes in those bands possibly contaminated by the nebular lines. The yellow circles and error bars mark the average values together with the associated standard deviation. The robust mean and associated standard deviation are quoted in parentheses in each panel. The largest excess correction is found with Smit et al. (2014) method, while the exclusion of the contaminated bands has the effect of introducing a scatter in the stellar mass measurement, although the average excess is nearly zero. The same behavior is observed for the case of including the old and dusty template and with the bayesian luminosity prior, although the lower number ob objects in the sample affects the statistics measurement.

3.6. Nebular line contamination

Recently, a number of works have investigated how the contamination to broad- and narrow-band photometry by nebular emission lines may introduce a systematic excess in the measurements of the stellar masses (and consequently reduce the sSFRs) of galaxies (see e.g., Atek et al. 2011; Shim et al. 2011; de Barros et al. 2014; Stark et al. 2013; Bowler et al. 2014; Straatman et al. 2014; Labbé et al. 2013b; Oesch et al. 2013; Ellis et al. 2013; Schenker et al. 2013; González et al. 2014), although a consensus on the effectiveness of this correction, in particular for galaxies at , has not yet been reached (see e.g., Labbé et al. 2013b; Bowler et al. 2014). Given the above uncertainties in the potential contamination from nebular lines at high redshift, we can not adopt a single and straightforward approach to correct for the nebular line contamination. Instead, the systematic effects that the contamination by nebular lines emission can have on the measurement of the stellar population parameters (in particular of the stellar mass) have been analyzed implementing three independent recipes to correct the observed photometry from potential contamination, with each recipe applied on a per-objects basis. Stellar masses and population parameters were then re-computed, adopting the same configuration used for the default sample. The result of this process are four different measurements of the stellar mass and population parameters associated to each sample of galaxies: one sample obtained with the original photometry and three new samples each one associated to one of the distinct methods of the nebular line contamination correction.

Figure 8.— Examples of stellar population parameters recovered from SED fitting on the photometric catalog after correction from nebular lines emission. The two rows refer each to a different object. In each panel, the photometry (in the observer frame) adopted for the fitting is represented by the colored squares; grey squares mark the photometry before applying the corrections. The pink curves represent the best-fit templates from FAST without applying any correction for emission line contamination; the violet curves in the remaining panels represent the best-fit template from FAST for each method implemented to correct for emission lines contamination. The best-fit template for the no-correction scenario is reported in each panel for comparison (pink curve). The blue curve in the leftmost panel marks the best-fit SED from EAzY. The main physical properties are listed at the top-left corner, together with the 68% confidence level uncertainties (see Figure 3 for further details). Left to right, the panels refer to the cases of original data, photometry corrected from the EW of lines in the best-fit EAzY template, photometry corrected following the procedure in Smit et al. (2014); and excluding from the photometry those bands being potentially contaminated by nebular emission.

In the first method, for each object we identified those bands which could be contaminated by the (redshifted) main nebular emission lines (Ly, H, H, [O II], [O III]). The potential contribution was computed from the line equivalent width (EW) recovered from the best-fit EAzY template, and the corresponding flux was then rescaled by the factor , where is the filter efficiency and is the redshifted wavelength of the emission line (see e.g., Eq. 1 in Smit et al. 2014). The emission lines in the EAzY templates are tuned for objects at redshifts , resulting in EW smaller than those observed at higher redshift. This method reflects then an optimistic scenario, as the contamination by nebular lines is likely larger than the values recovered through this configuration.

Figure 9.— Left panel: example of object for which the Smit et al. (2014) method for the correction of nebular emission contamination introduces an increase in stellar mass. The colored points correspond to the photometry (in the observer frame) after applying the correction of contamination. The original measurements are shown as filled grey squares. The best-fit FAST solution adopting the original photometry is shown by the pink curve, while the best-fit solution with the corrected photometry is represented by the violet curve. The main stellar population parameters are also reported, with the text colour matching the model they refer to. Specifically, top to bottom they are: the photometric redshift (), the , the SFR) in units of , the , the and the extinction expressed in magnitudes. Quoted errors refer to the 68% confidence intervals. Right panel: example of object for which the Smit et al. (2014) method for the correction of nebular emission contamination introduces an increase in stellar mass. Same plotting conventions as for the left panel. The increase in stellar mass is due to a significant increase in either the dust extinction or the age.

Recent works (see e.g., Holden et al. 2014; Smit et al. 2014) have shown that star-forming galaxies are commonly characterized by high line ratios and large equivalent width which can even more bias the measurements of stellar masses. Specifically, Smit et al. (2014) measured the rest-frame EW([O III]+H for a sample of galaxies at , finding a lower limit value of 637Å for the full sample, and an EW=1582Å for the bluest sources. These values were shown to be consistent with the extrapolation to of the EW(H measured by Fumagalli et al. (2012), assuming an evolution of the EW as and converting the EW(H into EW([O III]+H using the line intensity ratios in Anders & Fritze-v. Alvensleben (2003) for . Following Smit et al. (2014), our second method was then implemented as follows. At first we computed the rest-frame EW([O III]+H at the redshift observed for each galaxy applying the evolution to the EW([O III]+H measured by Smit et al. (2014) at . We adopted for the rest-frame EW of [O III]+H at the value EW([O III]+H=1582Å, which coincides with the higher EW from Smit et al. (2014), as a way to consider the highest nebular emission contamination still consistent with observations. The EW of all the nebular emission lines in Table 1 of Anders & Fritze-v. Alvensleben (2003) where then computed adopting the line intensity ratios corresponding to . The new flux was finally computed from Eq. (1) in Smit et al. (2014).

Our third method consisted in removing from the photometric catalog the fluxes in those bands which could be contaminated by the same nebular emission lines considered for the first method. This method is driven by the idea of not imposing any constraint on the contribution of the nebular lines in each photometric band.

The effects of the above three recipes on the stellar mass measurements are summarized in Figure 7, while in Figure 8 we present examples of SED fitting to recover stellar population parameters before and after correcting the photometry for nebular emission contamination. For the default sample, with the first approach we find an excess in stellar mass with respect to the stellar masses obtained without applying any correction for nebular line contamination, with biweight mean value of dex for the full sample ( dex, dex and dex for the stellar masses in the , and bins respectively). These values are approximately one order of magnitude smaller than previous determinations. Specifically, Schenker et al. (2013) report an excess of 0.2 dex for galaxies at ; at González et al. (2014) find a marginal correction, while Stark et al. (2013) report an excess of 0.04 dex; at the excess by González et al. (2014) is marginal, while Stark et al. (2013) find 0.1 dex; at both Stark et al. (2013) and González et al. (2014) find an excess of about 0.26 dex; finally an excess of 0.48 dex is found in Labbé et al. (2013b) for galaxies. However, the small values found for the excesses in stellar mass are consistent with the working hypothesis that the EW of nebular emission lines in the EAzY templates are smaller than those observed for high-redshift galaxies.

The biweight mean excess in stellar mass from the Smit et al. (2014) method is  dex for sources at ( dex,  dex and  dex in the three redshift bin respectively) which appear to be more consistent with the literature, although they are also consistent with no excess. We note that there is a group of galaxies which experienced an increase in stellar mass, rather than a decrease. For stellar mass , these are almost entirely objects that were detected on IRAC residual images and with redshift . An example, randomly extracted from the sample, in shown in the left panel of Figure 9. The SEDs of these objects are characterized by absence of flux in the optical bands and, by construction, little to no flux also in the UltraVISTA NIR bands. At these redshifts, the H+N II and H+[O III] enter the IRAC m and m bands. Under our working hypotheses, the EW associated to these nebular lines are large, resulting in a large correction factors in these two bands. Since the redshift is not re-computed after the correction is applied, the stellar population parameters are best fitted by a template with higher dust extinction, which translates into higher stellar masses.

A different effect seems to be responsible most often for the larger stellar masses estimated when adopting the Smit et al. (2014) method to correct for emission line contamination in galaxies with . An example is shown in the right panel of Figure 9. In this case, the galaxies are at and are characterized by extreme SFR (SFR ) and an extinction of about 1 mag. Here, the strong contribution from the Lyman- line in our model assumptions translates into large EW for the intermediate optical bands, resulting in low flux (generally consistent with 0 erg s cm Å) in the bands potentially contaminated by Lyman- emission, after the correction is applied. The resulting best-fit template is then characterized by a strong Ly absorption line, which excludes the high-star-formation solution. The observed red color is then explained by the resulting older population of stars, rather than the effect of dust extinction, which translates into lower SFR and and higher age.

The major effect in the measurements of the stellar mass introduced by our third method is to increase the spread in stellar masses, likely the result of the lower number of bands available for the fit. Indeed the measured excesses are  dex for and  dex,  dex and  dex in the three redshift bins. The lack of flux information during the fit is particularly critical for objects at . In this range of redshift H and H enter IRAC m and m bands. The simultaneous exclusion of these two bands from the fitting process introduces a higher degree of freedom in the choice of the best-fitting template. For those sources with a K-band flux value larger (in ) than the flux value in IRAC m and a plateau in the observer-frame NIR wavelength shorter than m, the fitting code favors a solution without a strong Balmer break and with high values of dust extinction. The combination of these two effects generates best-fit templates characterized by a very unlikely high SFRs, with common values about fewyr (see Figure 10 for an example). If, instead, the observer-frame NIR region presents flux decreasing with wavelength, the best-fit template does not show such extreme values of SFR and (see e.g., the lower-right panel of Figure 8). The effectiveness of the results from this third prescription were then considered on a per-case basis.

Figure 10.— Example of the effect of measuring the stellar population parameters excluding the bands possibly contaminated by nebular emission. Same plotting conventions as for Figure 9. The exclusion of the fluxes in the IRAC m and m bands and the fact that the flux in K-band is higher than the flux in IRAC m disfavors solutions with pronounced Balmer break. This results in best-fit templates with extremely high SFR.

When the three methods are applied to the sample obtained from the introduction of the old and dusty template and to the sample with both the old and dusty template and applying the bayesian luminosity prior, we find average excesses of , and and , and respectively for the three methods, in qualitative agreement with the results found for the default sample. The method introducing the largest excess in stellar mass is that from Smit et al. (2014).

The contribution of the nebular emission lines to the stellar mass measurements derived with the above scenarios are to be intended in a statistical sense. A quantitative and accurate determination of the contamination by emission lines for galaxies at would require spectroscopic studies on large samples, which is currently missing.

3.7. Star formation histories

We finally analyzed the systematic effects of adopting different SFHs in the measurement of the stellar population parameters; specifically, in addition to the delayed-exponential SFH, which constitutes out default configuration,stellar population parameters were also computed adopting an exponential SFH. The comparison between the stellar masses from the two SFHs is shown in Figure 11. For most of the galaxies, we observe no systematic offset in the stellar mass measurements. We can identify, however, a small sample of galaxies characterized by a stellar mass larger by  dex when using the exponential SFH. The observed SED and best-fit templates for both SFHs for one object randomly chosen among those showing the increase in stellar mass when adopting the exponential SFH is shown in Figure 12. The observed SED of these objects are characterized by a plateau in the NIR bands shorter of m. During the fitting process, this plateau can be described by either a very young, highly star-forming galaxy with strong extinction by dust, or by an object with lower dust content and SFR, but older stellar population. The best fit SED template from the delayed-exponential SFH is characterized by a slowly increasing SFH (/yr) with young age ((age/yr)), high SFR (SFRfewyr) and high dust extinction ( mag). The exponential SFH, instead, provides a best fit SED with older age ((age/yr)), and, correspondingly, lower dust extinction ( mag), resulting in a larger stellar mass. Although the delayed exponential SFH can, in principle, mimic the same best fit SED of the exponential SFH, as it can be seen from Figure 12, the solution with large is finally preferred as it provides a slightly better fit to the data in the NIR wavelength region.

Figure 11.— Excess in stellar mass adopting the exponential SFH compared to the stellar mass from the delayed-exponential SFH. The top panel refers to photometric redshifts computed without the old and dusty template and with no bayesian luminosity prior, the central panel to the inclusion of the old and dusty template but not applying the prior while the bottom panel refers to photometric redshifts obtained with the inclusion of the old and dusty template and with the adoption of the bayesian luminosity prior. The yellow circles and error bars mark the average values together with the associated standard deviation. The blue points identify those galaxies with sSFRyr. The robust mean and associated standard deviation are quoted in parentheses in each panel. For most of the galaxies and configurations, the different SFHs do not alter significantly the measurement of the stellar mass. In the cases without prior, a few massive galaxies have stellar masses larger by dex than the corresponding from the delayed exponential case.
Figure 12.— Example of objects whose stellar mass obtained assuming an exponential SFH has an excess of  dex compared to the stellar mass from the delayed-exponential SFH. The colored points mark the observed photometry in the observer frame; the pink curve marks the best-fit template from the delayed-exponential SFH, while the violet curve represents the bet-fit template from the exponential SFH. Stellar population parameters are also reported in the labels, coded by the color of the corresponding best-fit template. Other plotting conventions as in Figure 9. The exponential SFH provides an SED with older age and hence a higher stellar mass.

3.8. Stellar mass completeness

Figure 13.— Stellar mass completeness as a function of redshift: the black points represent the galaxies from the sample obtained without bayesian luminosity prior. The red solid, dashed and dotted lines respectively mark the 90%, 50% and 5% completeness limits for a passively evolving galaxy with , obtained considering the corresponding detection level in IRAC m, as described in the text, while the solid yellow line marks the 90% completeness limit in IRAC m with an additional absorption of  mag, and corresponds to the limit in stellar mass we adopted in this work. The inset shows the detection completeness as a function of apparent magnitude for a point source in IRAC m (grey dotted curve) and IRAC m (red solid curve), with 90% confidence levels equal to 23.4 mag and 22.9 mag respectively. The depth of the IRAC maps allows us to consider only galaxies with .

In Figure 13 we mark with the solid red curve the completeness in stellar mass corresponding to the evolution from of a simple stellar population (SSP) model with stellar mass corresponding to that computed from the 90% completeness limit in IRAC m. However, as shown by the inset, the 90% completeness limit in IRAC m is fainter by 0.5 mag than in IRAC m (IRAC m and m 90% completeness limits are 23.4 mag and 22.9 mag respectively). Furthermore, recently there has been evidence of significant dust extinction at the high-mass end of galaxies (see e.g., Whitaker et al. 2010; Marchesini et al. 2010, 2014). Because of these two factors, we adopt for the completeness limit the orange solid curve in Figure 13, which includes the effect of 2 mag of extinction in the V-band to the completeness from IRAC m. As it is shown in Figure 13, the current depth of the IRAC coverage to the COSMOS field only allows for and galaxies at and respectively, corresponding to the high-mass end of the SMF.

Figure 14.— These plots show the SED of the stellar mass complete sample (one object per row) obtained adopting the most conservative configuration, i.e., with the inclusion of the old and dusty template in the set of templates used for the measurement of photometric redshifts, and with the application of the bayesian luminosity prior in the measurement of photometric redshifts. Each panel refers to different fluxes adopted for the measurement of the stellar population parameters, as explained in Figure 8.
Figure 14.— - Continued

3.9. Cosmic Variance

The total area covered by COSMOS/UltraVISTA is approximately square degrees in one single field, making the effects of cosmic variance not negligible for very massive galaxies. The contribution of cosmic variance to the total error budget was estimated through the recipe of Moster et al. (2011). The average uncertainties due to this effect are 21%, 28% and 37% for the , and redshift bins respectively for stellar masses . These values were added in quadrature to the poissonian error of the SMFs.

Figure 15.— The robust massive galaxy (corresponding to ID 43320 in Figure 14). The three panels on the left show the SED from the measured photometry in the observer frame (colored points with error bars), together with the best-fit template from EAzY (blue curve) and from FAST (pink curve), this latter obtained without applying any correction for nebular emission contamination to the photometry. The main physical properties are listed at the top-left corner (see Figure 3 for further details). Each one of the three panels refers to a different configuration adopted for the measurement of the photometric redshift. Left to right, these represent the cases of: excluding the bayesian luminosity prior and the old and dusty template; excluding the bayesian luminosity prior, but introducing the old and dusty template among the set of SED templates adopted for the photometric redshift measurements; activating the bayesian luminosity prior and including the old and dusty template. The panel on the right shows the redshift probability distribution for the three cases. The inset ( wide on each side), centred at the position of the object, shows the results from stacking the filters bluer than the Lyman limit. The best-fit templates well describe the photometric data. The are characterized by a narrow peak at , with a secondary peak at which appears when the luminosity prior is introduced, but whose probability is .

4. The population of galaxies

4.1. Robust massive galaxies at

In this work, we consider as massive a galaxy with , while quiescent a galaxy whose specific star formation rate (sSFR) is smaller than , with the Hubble time at redshift (Damen et al., 2009; Lundgren et al., 2014). Our stellar mass complete sample is then composed by massive galaxies only. However, the different configurations adopted for the computation of redshift and stellar masses produce samples which, in general, do not always contain the same galaxies. As we discussed in Section 3, systematic effects can arise in the measurements of photometric redshifts and/or stellar population parameters depending on the adopted set of SED templates, the inclusion of the bayesian luminosity prior, and the inclusion of the potential contamination from nebular lines.

Our analysis allows us to identify seven robust massive galaxies with redshift measured adopting the most restrictive configuration, i.e., with the adoption of both the bayesian luminosity prior and the maximally red SED template. Their SEDs are presented in Figure 14. For each object, the panel on the left presents the original photometric points together with the best-fit SED template from EAzY (solid blue curve) and the best-fit FAST template with the default assumptions (dash-dotted pink curve). The three panels on the right show the results of the FAST run on the photometric catalogs obtained after applying the correction from nebular emission lines contamination. We note that, out of the seven objects, only three of them (ID 43320, 196141 and 203033) would be included in the stellar-mass complete sample irrespectively of the assumptions adopted for the measurement of photometric redshifts and stellar population parameters.

The sample is characterized by stellar masses in the range , specific star-formation rates in the range , extinction in the range  mag and ages age/yr (obtained with the default assumptions: no luminosity prior, no old/dusty template and delayed exponential SFH), with median of the log-values of ,  yr and yr for the stellar mass, sSFR and age, respectively and consistent with what found by Bowler et al. (2014), while the median of the dust extinction is 0.8 mag. Specifically, the high value of the dust extinction is qualitatively in agreement with the trend of increasing with redshift found by e.g., Whitaker et al. (2010) and Marchesini et al. (2014). Most of the SEDs are characteristic of star-forming or post-starburst galaxies, and their redshift probability distributions show very pronounced peak, with a small (where non absent) secondary peak at . The redshift probability distribution are well constrained because most of these objects show both the Lyman break and the Balmer/4000Å break (although, given the young ages, this second break is mostly originated by the Balmer break).

4.2. Massive and quiescent galaxies at

The adoption of a set of different configurations for the measurements of photometric redshifts and stellar population parameters, as described in Sect. 3, in principle can produce different values of the physical properties from the same photometry. Specifically, the effect of both the old/dusty SED template and the introduction of the bayesian luminosity prior result in photometric redshifts generally lower than those obtained when excluding the luminosity prior and the old/dusty SED template.

Despite the above indetermination, we are able to identify one robust candidate for a massive galaxy at irrespectively of the configuration adopted for the measurements of photometric redshift and stellar population parameters, whose sSFR from SED fitting is consistent with being quiescent. Its SED is presented in the top row of Figure 14 labeled as ID 43320; in Figure 15 we show the best-fit SEDs from EAzY for the three cases of no luminosity prior and no old/dusty template, no luminosity prior and old/dusty template, and luminosity prior and old/dusty template (blue curves) and corresponding FAST best-fit templates (pink curves), together with the redshift probability distributions for the three above cases (filled regions in the panel on the right).

In the following sections we analyze in more detail its main physical parameters.

4.2.1 Observed SED

The photometric SED for this object is presented in the top row of Figure 14, labeled as ID 43320. The flux in the IA527 filter presents an excess compared to that in the adjacent filters. Although this would be consistent with a Lyman emission, a visual inspection of the science frame did not show any clear evidence for emission from an object. Instead, the same region is crossed by an horizontal band, possibly an instrumental defect, few pixels wide and characterized by an emission slightly above the local background, which is likely the origin of the measured flux excess.

The MIPS 24m map shows a clear detection at level, indicating either that a solution at could be more appropriate or that the rest-frame emission at m originates from the torus hot dust heated by a hidden AGN.

As we show in Section 4.2.4, the MIPS m emission measured for this source corresponds to high infra-red luminosity. If this excess originated from obscured star formation, the high IR luminosity would make the galaxy detected in Herschel maps. The HerMES SPIRE 250m map (Oliver et al. 2012; Viero et al. 2013) does not show any significant source centered at the position of this galaxy; according to the HerMES catalog (Smith et al. 2012; Roseboom et al. 2010; Wang et al. 2014), the closest source lays at an apparent distance of about . However, the large FWHM of SPIRE m (FWHM) does not allow us to completely rule out the absence of SPIRE 250m flux for this galaxy.

4.2.2 Redshift

In Figure 15 we show the best-fit SEDs from EAzY for the three cases of no luminosity prior and no old/dusty template, no luminosity prior and old/dusty template, and luminosity prior with old/dusty template (blue curves); the panel on the right shows the redshift probability distributions for the three above cases (filled regions). The galaxy has a photometric redshift consistent with , depending on the configuration adopted for the measurement. The best-fit SEDs well describe the photometric points, supporting the redshift measurement. In the right panel of Figure 15 we also show the stack of those bands bluer than the Lyman limit under the assumption that . The resulting image does not present any clear evidence of flux excess, increasing the confidence on the measured value for the redshift. According to the distribution of redshift resulting from the SED fitting, the probability of this galaxy to be at is smaller than 21% (this upper value corresponds to the introduction of the bayesian luminosity prior and of the dusty template). However, forcing the redshift to be produces a best-fit solution which is worse than the one.

4.2.3 Stellar mass

The recovered stellar population parameters vary both because of the different redshift measurement and because of the different recipes we adopted to take into account the potential contamination by nebular lines. Despite this, the stellar mass measurements for this object are all consistent with the value of . The FAST best-fit SEDs are marked by the pink curves in Figure 15. The stellar mass changes only marginally even when assuming that IRAC m and m are contaminated by strong emission lines, as from Smit et al. (2014) recipe for which the rest-frame EW([O III]+)1230Å at . We note that the measurements of the stellar population parameters obtained excluding from the fit those bands potentially contaminated by nebular emission when no luminosity prior is adopted for the measurement of photometric redshifts suffer from the issue presented in Section 3.6 (see also Figure 10): the exclusion of the IRAC m and m fluxes, together with the fact that the -band flux is higher than the IRAC m disfavors a solution with pronounced Balmer/4000Å break; the best-fit SED is instead characterized by extremely high SFR (SFRyr). The corresponding stellar population parameters were then excluded during the selection process. A definitive assessment of the intrinsic physical properties of this object will necessarily require spectroscopic observations.

4.2.4 sSFR

The values for the sSFR recovered from the SED analyses are all consistent with log(sSFR/yr)=, with and an extinction  mag. Similarly to the stellar mass measurement, also the sSFR change marginally under the assumption that IRAC m and m are contaminated by strong emission lines. At redshift the sSFR values satisfies the criterion of sSFR for the identification of quiescent galaxies. However, the value of the sSFR from the SED fitting analysis is potentially in contrast with the observed MIPS m flux.

Assuming that the observed MIPS 24m flux comes entirely from the dust-enshrouded star formation, at the measured redshift, it corresponds to a luminosity (adopting the recipe in Wuyts et al. (2008) with Dale & Helou (2002) template set; no significant discrepancy was obtained using the Chary & Elbaz (2001) recipe). Using this value of the infrared luminosity, we estimate the star-formation rate to be SFRyr for a Kroupa (2001) IMF (Kennicutt 1998; Bell et al. 2005; Muzzin et al. 2013a), an unlikely high physical value. If instead we assume for this galaxy a redshift , roughly corresponding to the secondary peak in the distribution, its total infra-red luminosity () and the SFR would be typical of hyper-luminous infra-red galaxies (HLIRGs). Indeed, using the same prescriptions used above, the luminosity recovered from the MIPS 24m flux would be , with a star-formation rate SFRyr for a Kroupa (2001) IMF.

If the MIPS emission of this galaxy originated from obscured star formation, the high infra-red luminosity at , and most likely even that at , would make the galaxy detected in Herschel data. However, as discussed in Section 4.2.1, current HerMES data do not show any clear evidence of flux which could be associated to this object. This fact, together with the unlikely high SFR recovered from the observed MIPS m when the galaxy is considered to be at redshift , suggest that the observed MIPS flux is likely not originated by an intense SFR, but, instead, by the torus of hot dust heated by a hidden AGN, supporting the SFR values recovered from the SED fitting analysis. However, a robust assessment of the presence of an AGN (and hence of the origin of the MIPS excess) requires spectroscopy, which is still lacking.

4.3. Stellar Mass Functions

We computed the SMFs in the three redshift bins , and using the 1/Vmax formalism (Schmidt, 1968; Avni & Bahcall, 1980), which provides both the shape and the normalization of the SMF at the same time. Upper and lower poissonian uncertainties were computed using the recipe of Gehrels (1986), valid for small samples. Uncertainties due to cosmic variance computed as described in Sect. 3.9 were added in quadrature. Estimating how the uncertainties in the stellar mass measurements propagate into the computation of the SMF is a non-trivial task since the stellar mass is not directly observed, but it is instead recovered through SED template fitting on multi-wavelength photometry. Similarly to what done in Muzzin et al. (2013b), we then estimated the effects that flux uncertainties have on the measurement of the SMF using Monte Carlo simulations. We implemented 100 Monte Carlo realizations of the multi-wavelength photometric catalog randomly perturbing the measured fluxes according to the flux uncertainties. Photometric redshifts and stellar population parameters were then re-computed using the same methods applied to the original unperturbed sample. For each one of the 100 catalogs, the SMFs were finally measured. The dispersion of the density values of the SMF in each stellar mass bin then gives an estimate of the uncertainties in the SMF from photometry errors propagating to photometric redshifts and stellar mass measurements. The uncertainties ranged from 8% to 40%. These uncertainties were added in quadrature to the uncertainties from poissonian noise and cosmic variance. Furthermore, given the high uncertainties in the dust content of the SED templates adopted for the recovery of the photometric redshift, we caution the reader that the lower stellar mass bin in each redshift range may still suffer from incompleteness.

Figure 16.— Stellar mass functions at (top panels), (middle panels), (bottom panels). The plots in the left column refer to the SMF of the robust sample of massive galaxies presented in Section 4.1, while in the column on the right we present the SMF measurements for different configurations adopted for the measurement of photometric redshifts and stellar population parameters. For the measurements in our work, the vertical error bars include the effect of cosmic variance and the systematic uncertainties arising from the adopted SFHs, the correction of the nebular lines in the photometry, the inclusion of the old and dusty template and the application of the bayesian luminosity prior for the measurements of the photometric redshifts. The horizontal error bars reflect the bin size adopted for the computation of the SMF. When the bayesian luminosity prior is introduced in the computation of the photometric redshifts, most of the resulting SMFs measurements turn into upper limits. The measurements from the luminosity prior case have been arbitrarily offset by -0.05 dex to help visualisation. The SMF at is consistent with previous measurements at level up to .
Figure 17.— Comparison of the SMF in the three redshift bins. The top panel refers to photometric redshifts obtained without the bayesian luminosity prior, while the lower panel refers to the case with the bayesian luminosity prior. In each panel, the points refer to the median value of the SMF measurements obtained with all the different configurations adopted to compute the stellar masses (i.e., using the measured flux as well as after applying the correction for nebular emission contamination, as described in Section 3.6. For comparison, the SMF measurements from Muzzin et al. (2013b) for are also plotted. Specifically, the blue shaded region corresponds to the SMF. No evidence for evolution is found in the SMFs from to .

In Table 1 and 2 we report the measurements with the total errors of the SMFs for all the combinations for the measurements of photometric redshifts and stellar population parameters adopted in this work for the delayed-exponential and exponential SFH, respectively. A graphical representation of the SMFs is shown in Figure 16. The column on the left shows the SMF of the robust sample of massive galaxies presented in Section 4.1, while in the right column we present the SMF measurements from different configurations adopted for the measurement of photometric redshifts and stellar population parameters. For each stellar mass bin, points correspond to the median of the measurements from the different configurations adopted for the measurement of photometric redshift and/or for the correction of nebular emission contamination; upper (lower) error bars correspond to the maximum (minimum) value. A direct comparison among the SMFs in the three redshift bins is presented in Figure 17. We however note that, given the large uncertainties in the SED fitting arising particularly at from the exclusion of the flux in those bands potentially contaminated by nebular emission (see discussion in Sect. 3.6), the values reported in Figures 16 and 17 do not include the measurements of the SMF obtained with this configuration, although, for completeness, these values are included in Table 1 and 2.

range Central Mpc dex
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Table 1SMF measurements in the three redshift bins, corresponding to the delayed-exponential SFH and for the four different configurations adopted for the measurements of the photometric redshifts, namely without luminosity prior and old/dusty template, without luminosity prior but including the old/dusty template, with luminosity prior but without old/dusty template, and with luminosity prior and the old/dusty template. For each photometric redshift measurement, the SMF corresponding to each of the four configurations adopted for the measurement of the stellar population parameters and the correction of nebular line contamination are reported. Uncertainties include poisson noise and cosmic variance.
range Central Mpc dex
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
No prior, No old/dusty No prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Prior, No old/dusty Prior+old/dusty
Default EAzY lines Smit+2014 Excl. bands Default EAzY lines Smit+2014 Excl. bands
Table 2SMF measurements in the three redshift bins, corresponding to the exponential SFH and for the four different configurations adopted for the measurements of the photometric redshifts, namely without luminosity prior and old/dusty template, without luminosity prior but including the old/dusty template, with prior but without old/dusty template, and with luminosity prior and the old/dusty template. For each photometric redshift measurement, the SMF corresponding to each of the four configurations adopted for the measurement of the stellar population parameters and the correction of nebular line contamination are reported. Uncertainties include poisson noise and cosmic variance.

The configurations we adopted for the measurements of redshift and stellar population parameters, described in Section 3, introduce a dispersion in our SMF measurements. The considered SFHs introduce the smaller dispersion in the SMF measurements; the dispersion of the stellar mass measurements from the inclusion of the nebular emission lines is smaller than the uncertainty from the combination of poissonian noise and cosmic variance. On the contrary, the inclusion or not of the bayesian luminosity prior in the measurements of the photometric redshifts influences the values of photometric redshifts themselves for a non-negligible fraction of objects, and, as a consequence, of the stellar masses. Specifically, the adoption of the bayesian luminosity prior transforms many measurements into upper limits.

In the following subsections we compare our SMF estimates with measurements from the literature.

4.3.1 The SMF at

In the redshift range , SMF measurements have already been published by Stark et al. (2009); Caputi et al. (2011); González et al. (2011); Santini et al. (2012); Lee et al. (2012). Of these, the SMFs of Stark et al. (2009), González et al. (2011) and of Lee et al. (2012) are based on dropouts selections, while those of Caputi et al. (2011) and Santini et al. (2012) are based on photometric redshift measurements. The measurements by Caputi et al. (2011) are based on observations over relatively wide fields (0.6 deg), while the other works over fields about one order of magnitude narrower than UltraVISTA. Since the total exposure time is roughly of the same order, this translates into detection of fainter sources, and allowing for a lower stellar mass completeness limit for the same stellar mass-to-light ratio. The smaller fields together with the selection effects introduced by the dropout technique, have reduced so far the chances of detecting the most massive objects. This is evident from the top-left plot of Figure 16. Points from the literature reach an upper limit in stellar mass which is at most , the only exception here being the measurements at by Santini et al. (2012), although the error bar is quite large. On the opposite side, González et al. (2011) measure the SMF up to , but their stellar mass range never overlaps with ours.

The stellar mass range over which our SMF measurements is defined overlaps with those by Stark et al. (2009), Caputi et al. (2011) and Lee et al. (2012) in the lowest stellar mass bin, while it overlaps with Santini et al. (2012) up to . The lowest stellar mass bin measurements with no luminosity prior is consistent with the measurements by Stark et al. (2009) and with Caputi et al. (2011) and at -level with the measurements by Lee et al. (2012), while it is consistent with Santini et al. (2012) at 3. We however note that Santini et al. (2012) do not include the effects of cosmic variance in the error budget, which for we estimate using the recipe of Moster et al. (2011) to be as high as 50% and 70% for and 11 respectively. At higher stellar mass, our SMF becomes consistent with both the Schechter parameterization and Vmax measurements by Santini et al. (2012) at 1 level. Under the caveat that extrapolations are characterized by a high degree of uncertainty and should be considered with care, we finally compare our measurements to the extrapolation of the Schecher fit of Lee et al. (2012). While the measurements obtained without the bayesian luminosity prior are