© ESO 2019

1. Introduction

The formation of galaxies and their properties are driven by the evolution of both dark matter and baryons. While dark matter halos assembled most of their mass through sequential merging, some fundamental questions remain still open (Naab & Ostriker 2017), particularly those regarding the formation of massive galaxies and their transformation into quiescent systems.

Recent developments support a scenario where galaxies assembled most of their mass not only through sequential merging, but also through smooth accretion by cold gas streams that penetrate the shock-heated media of massive dark matter halos and grow dense, unstable, turbulent discs with bulges, and trigger rapid star formation (e.g., Dekel et al. 2009).

Different models of galaxy formation and evolution have been developed in order to explain the emergence of massive objects (e.g., Sheth et al. 2001; Henriques et al. 2015) and their transformation into quiescent systems (e.g., Guo et al. 2011, 2013; Naab et al. 2014; Henriques et al. 2017). However, when these models are applied to simulations (e.g., the Munich Simulation: Kauffmann et al. 1999; Springel et al. 2001; Croton 2006; De Lucia & Blaizot 2007; Henriques et al. 2015; Croton et al. 2016), they predict results which are not fully consistent with the observed evolution of galaxy properties. In particular, the presence of a population of high-redshift massive (and sometimes quiescent) galaxies can create tension between models and observations (e.g., the impossibly early galaxy problem, as in Steinhardt et al. 2016, referring to the problem of such rapid assembly of massive systems). Therefore, a robust determination of the abundance of massive quiescent galaxies at high redshift is a powerful test bench for galaxy formation models that also have to reproduce the mechanisms for quenching star formation in order to produce quiescent galaxies.

Many recent studies (e.g., Daddi et al. 2004; Ilbert et al. 2006; Fontana et al. 2006; Wiklind et al. 2008; Caputi et al. 2012; Castro-Rodríguez & López-Corredoira 2012; Stefanon et al. 2013; Muzzin et al. 2013; Nayyeri et al. 2014; Straatman et al. 2014; Mawatari et al. 2016; Wang et al. 2016; Merlin et al. 2018, 2019; Deshmukh et al. 2018) have identified a population of massive quiescent galaxies at high redshift (z  ≳  2.5). Due to the current observational limits, these objects are just photometric candidates based on spectral energy distribution (SED) fitting analysis of broad-band photometry with population synthesis models. There are few notable exceptions for which deep near-infrared (NIR) spectroscopy confirms the passive nature of the candidates and determines their redshifts (at z  =  3 in Gobat et al. 2012, and z  =  3.7 in Glazebrook et al. 2017). Moreover, in a recent work Schreiber et al. (2018) obtain Keck–MOSFIRE H and K-band spectra for 24 candidate quiescent galaxies at 3 <  z <  4 and confirm the redshift and passive nature for eight of them, adding further evidence to the existence of a significant population of quiescent galaxies at high-redshift. This is a strong evidence that star formation in these galaxies occurred very fast within the first billion years of the universe, followed by a passive evolution (Mancini et al. 2009; Stefanon et al. 2013; Straatman et al. 2014; Oesch et al. 2013). The common technique adopted to select and identify quiescent galaxies at z >  3 consists of two steps: the selection of likely candidates based on colors (observed or rest-frame) and the subsequent identification of the most likely old and high redshift galaxies among these candidates by means of a SED fitting analysis with population synthesis models (e.g., Wiklind et al. 2008; Nayyeri et al. 2014; Wang et al. 2016). This two-step process is required because, as it is discussed further in this paper, the colors of these high redshift quiescent galaxies are to some extent degenerate with regard to those of dusty star-forming galaxies at the same or lower redshift (as is also known for quiescent galaxies at lower redshift, e.g., Pozzetti & Mannucci 2000). A recent example is the work by Merlin et al. (2018) in which the authors search for passive galaxies at z >  3 in the GOODS-South field using photometric data and then select galaxies with star formation rate SFR = 0 with the aid of a SED fitting analysis adopting top-hat star formation histories. All these results demonstrate the increasing importance of the quest for identifying the population of passive and quiescent galaxies at high redshift.

The aim of this paper is to devise observed color-selection diagrams that can efficiently identify the bulk of the population of z ≥ 2.5 massive quiescent galaxies, including those that could be missed by other selection techniques. In Sect. 2, we present new observed color selections to identify quiescent galaxies at z ≳ 2.5. We carefully studied the effects of all the parameters characterizing evolutionary tracks, which is presented in Appendix A. The data used to test the new criteria are given in Sect. 3. In Sect. 3.3, we present the selection of the quiescent candidates through color selections and in Sect. 4, the SED fitting analysis needed to ascertain their passive nature. In Sect. 5, we revise the number density of quiescent objects at high redshift and compare our results to recent galaxy formation models. Moreover, in Sect. 6, we present our estimate of the mass functions for the quiescent population. A discussion on our results and a comparison with previous observations and current semi-analytic models is presented in Sect. 7. A summary of the methods adopted and the results is presented in Sect. 8. Throughout this paper, a standard cold dark matter (ΛCDM) cosmology is adopted with Hubble constant H₀ = 70 km s⁻¹ Mpc⁻¹, total matter density Ω_M = 0.3 and dark energy density Ω_Λ = 0.7. All magnitudes are expressed in the AB system and log is base 10 logarithm, if not otherwise specified.

2. Color selections at z > 2.5

One of the most efficient ways of selecting quiescent galaxies is based on colors, a method which relies on the identification of spectral features via photometry.

Color–color selection can be performed using rest-frame or observed quantities. The approach with rest-frame colors, such as UVJ (Williams et al. 2009) and NUVRJ (see Laigle et al. 2016 or Ilbert et al. 2013 for details), requires an estimate of the redshift, which for large and deep surveys is usually derived from photometry, with rest-frame colors evaluated using the observed SED. High-redshift quiescent galaxies are characterized by SEDs which make them often too faint to be detected at optical wavelengths, and therefore rest-frame colors are usually extrapolated (Ilbert et al. 2013). Instead, apparent colors (e.g., Nayyeri et al. 2014) have the advantage of relying only on observations, although the lack of redshift measurement can add a further level of degeneracy. In this paper we will follow this second approach, in order to pre-select quiescent candidates based only upon observed data, and then an ad-hoc SED fitting analysis is performed in order to confirm their nature.

To identify quiescent galaxies by means of observed colors we use the strength of the Balmer and D4000 breaks (e.g., Bica et al. 1994), redshifted to NIR wavelengths for 2.5 ≲ z ≲ 4.5. The Balmer break at λ_rest = 3646 Å is due to the absorption by Balmer series down to the Balmer limit (strongest in A-type stars). The so-called D4000 break at 4000 Å is mostly produced by H and K absorption lines of calcium at λ_rest = 3968 Å and λ_rest = 3933 Å respectively, characteristic of solar-type stars, along with several metallic absorption lines. Both discontinuities can be used to select galaxies which ceased star formation (post-starburst, quiescent passive) since their breaks are much more prominent than in star forming galaxies.

2.1. Color predictions with evolutionary tracks

To identify the locus populated by high-z passive galaxies on a color–color diagram we used evolutionary tracks derived with stellar population synthesis (SPS) models. Some previous work on color selection determination (e.g., Nayyeri et al. 2014; Wiklind et al. 2008; Wang et al. 2016; Mawatari et al. 2016) adopts a fixed age for the stellar populations and follows their color evolution with redshift. In our work, we decided to follow a more physically-motivated approach, letting colors evolve also with age after having fixed a formation redshift in order to avoid ages larger than the age of the universe at the considered redshift. Therefore, in our approach we have to consider several values of formation redshift to explore the whole parameter space of the evolutionary tracks, as done in Appendix A.7.

To parametrize quiescent objects, we computed color tracks from BC03 (Bruzual & Charlot 2003) models, adopting a Chabrier IMF (Chabrier 2003), solar metallicity (Z  =  Z_⊙  =  0.02), and exponentially declining star formation histories with different e-folding times τ. Star-forming galaxies are generally assumed to be reproduced by a constant star formation template: we built evolutionary tracks of this galaxy template using the code fsps (Conroy et al. 2009, 2010), including the contribution of nebular emission lines that are evaluated with the results of Byler et al. (2017) based on CLOUDY (Ferland et al. 2013). Among the free parameters available in the model, we chose to adopt the Chabrier IMF, Padova’s stellar isochrones, and BaSeL libraries for better consistency with BC03 models. We considered the aging of both the star-forming and the quiescent templates. Therefore, a redshift of formation must be fixed when computing the color tracks as a function of redshift, and we chose as our reference z_form  =  6. We also allowed for different values of dust extinction A_V applied through the Calzetti’s attenuation law (Calzetti et al. 2000).

We tested the robustness of the selection criteria against different choices of parameters for extinction laws, initial mass functions (IMFs), star formation histories (SFHs), metallicities, stellar population synthesis (SPS) models (e.g., Maraston 2005; hereafter M05), redshift of formation, and also considering whether the inclusion (or not) of emission lines changes the selection criteria in Appendix A.

We systematically looked for optimal selections using near- and mid-infrared observed colors. It is indeed necessary to use red bands and/or near/mid-infrared bands for three main reasons: first, because quiescent galaxies are characterized by intrinsically red SEDs; second, because the effect of redshift and k-correction makes these objects almost undetected in optical blue bands due to their faint UV emission; and finally, because at z >  3, the radiation blue-ward the Lyman-limit is largely suppressed by the inter-galactic medium (IGM) absorption and, therefore, all the galaxies start to become suppressed in the u filter (Guhathakurta et al. 1990; Steidel et al. 1996, 1999). In particular, we used the photometry in the bands i, z, Y, J, H, K_s, IRAC[3.6] and IRAC[4.5], with transmission curves of the instruments which have observed the COSMOS field¹, and whose characteristics are listed in Table 1.

Based on the vast exploration of possible colors, we selected two criteria that perform well in selecting quiescent objects at 2.5 ≲ z ≲ 4.5. We stress that the purpose of this work is to be as inclusive as possible: we decided to maximize the completeness of the quiescent galaxies selection by including objects that can be missed by other criteria. The consequence of this choice will be the lower purity of the sample; contamination from other types of objects will have to be removed with the aid of an SED fitting.

2.2. JK_s[3.6][4.5] selection

We designed the color selection using the observed colors J − K_s and IRAC[3.6]−IRAC[4.5] (hereafter [3.6]−[4.5]) to select quiescent galaxies in the redshift range z ≈ 2 − 4. At z ≈ 2 the Balmer break is redshifted at λ ≈ 1.09 μm, while for z ≈ 4 to λ ≈ 1.82 μm; the D4000 break is redshifted at λ ≈ 1.2 μm at z ≈ 2 and at λ ≈ 2.0 μm at z ≈ 4. The color J − K_s can therefore identify the breaks, while the color [3.6]−[4.5] has been chosen since it best complements the other color to avoid overlaps between different tracks.

In the top panels of Fig. 1, two SEDs of quiescent galaxies at fixed ages and redshifted to different z with no dust extinction are shown to illustrate the location of the D4000 and Balmer breaks for the objects of interest with respect to the used filters. The cyan line shows a SED of a quiescent galaxy of 0.5 Gyr redshifted to z = 4 while the magenta line shows a SED of a quiescent galaxy of 1 Gyr at z = 3. Filter transmission curves are also shown in order to indicate the location of the Balmer and D4000 breaks with respect to the chosen colors.

Fig. 1. Top panels: two SEDs of quiescent galaxies. SED in magenta represents a population of 1 Gyr redshifted to z = 3 (i.e., zform ≈ 6), while the cyan line shows the SED of a population of 0.5 Gyr redshifted to z = 4 (i.e., again zform ≈ 6). Filters transmission curves used in COSMOS field and in the tracks evaluations are also shown. Bottom panels: JKs[3.6][4.5] and HKs[3.6] color–color diagrams with different evolutionary tracks both for star-forming galaxies and passive galaxies. Tracks representing star-forming galaxies are shown in cyan, blue and black with AV = 0, AV = 3, and AV = 5 respectively. Three tracks for quiescent galaxies (with exponentially declining star formation history τ = 0.1 Gyr and solar metallicity) are shown in magenta, orange and red with different extinction values (AV = 0, AV = 1, and AV = 2 respectively). Gray shaded areas represent the selection region for quiescent galaxies at 2.5 ≲ z ≲ 4 in the left panel and at 3 ≲ z ≲ 4.5 in the right panel. The numbers near the tracks of the same color represent the redshift. A vector corresponding to a magnitude extinction of AV = 1 using Calzetti’s law is also shown.

Fig. 2. Stellar mass as a function of redshift for different theoretical parametrizations of quiescent galaxies: we adopted an exponentially declining star formation history with characteristic time τ = 0.1 Gyr, Z = Z⊙ = 0.02, AV = 0, zform = 6. Red and black lines correspond to magnitudes m[4.5] = 22.27 and 24, respectively. Orange, green, and magenta lines represent different models of a quiescent galaxy differing in redshift of formation, metallicity, and e-folding time of the SFH, respectively.

We identified the following criteria to select quiescent galaxies in the redshift range z = 2 − 4 and limit the contamination from star-forming galaxies:

$$\begin{array}{c}\hfill \{\begin{array}{c}([3.6]-[4.5])\ge 0\hfill \\ (J-K_s)\ge 0.975([3.6]-[4.5])+1.92\text{for}[3.6]-[4.5]<0.4\hfill \\ (J-K_s)\ge 9.0([3.6]-[4.5])-1.33\text{for}[3.6]-[4.5]\ge 0.4\hfill \end{array}\end{array}$$(1)

which correspond to the gray shaded area in the left panel of Fig. 1. In the same panel, we show the color–color diagram and different evolutionary tracks built as described in the previous section. More precisely, we show three different tracks parametrizing quiescent galaxies with different extinction values. In particular, magenta tracks show a track with A_V = 0, orange a track with A_V = 1 and the red a track with A_V = 2. The choice of A_V = 2 is an extreme case: in fact, we do not expect to find any quiescent candidate with such high extinction value. However, in order to be as inclusive as possible, we do not exclude the region covered by this track. High-extinction star-forming galaxies can, in principle, contaminate the sample, but even though the black track in Fig. 1, which reproduces the color evolution of a galaxy with constant star formation history and A_V = 5, falls within the selection region, these interlopers are at lower redshift (z ≲ 2) and can be removed from the sample once SED fitting is performed.

2.3. HK_s[3.6] selection

This color selection is designed to identify quiescent galaxies at slightly higher redshift with respect to the previous color selection, that is, in the range z = 3 − 4.5, using filters H, K_s and IRAC 3.6 μm. From z ≈ 3 to z ≈ 4.5, both the Balmer and the D4000 break are located in the wavelength range bracketed by H and K_s filters (from λ = 1.64 μm to λ = 2.15 μm). Since the D4000 and Balmer breaks develop on different timescales, the color H − K_s can be used to select quiescent galaxies for broad ranges of age. In the top-right panel of Fig. 1, two SEDs of quiescent galaxies at fixed ages and redshifted to different z are shown to illustrate the location of the breaks with respect to the filters used.

In the bottom-right panel of Fig. 1, we show the color–color diagram and different evolutionary tracks built as described in the previous section. As for the previous color selection, the color K_s − [3.6] has been chosen in order to best complement the information of the other color.

To select quiescent galaxies in the redshift range z = 3 − 4.5, and to exclude most star-forming galaxies, we designed the following selection criteria:

$$\begin{array}{c}\hfill \{\begin{array}{c}{K}_s-[3.6]\ge 0.2\hfill \\ H-K_s\ge 0.617(K_s-[3.6])+0.58\hfill \end{array}\end{array}$$(2)

which correspond to the gray shaded area in the right panel of Fig. 1. Star-forming galaxies may enter the color selection region for extreme values of extinction (A_V ≳ 5) at lower redshift: the contamination by dusty and star-forming galaxies cannot be completely excluded, but they can be removed with the SED fitting analysis.

3. Application to real data

We applied the proposed color–color diagrams to real data. To this aim, we used the multi-waveband deep photometric observations available in the COSMOS field (Scoville et al. 2007).

3.1. The COSMOS2015 catalog and colors

We used the COSMOS2015 catalog (Laigle et al. 2016), which is the latest public release of data in the COSMOS field. This catalog provides the deepest optical and infrared observations of the field and the photometry has been obtained using a detection image evaluated by combining NIR images of UltraVISTA (YJHK_s) with the optical z⁺⁺-band data from Subaru. Since deep optical and IR data are especially important for the purpose of this paper, the characteristics of the catalog are particularly suitable for our work. In COSMOS2015, optical and NIR photometric data are provided in different apertures: 2″, 3″ (denoted APER2 and APER3 respectively), ISO (isophotal), and AUTO apertures. AUTO apertures are object to object variable apertures, automatically evaluated by SExtractor (Bertin & Arnouts 1996) in order to approximately measure the total flux of the object. IRAC data are instead provided only as total magnitudes and fluxes.

In the color–color diagrams presented in Sects. 2.3 and 2.2, we combine NIR and IRAC data and, therefore, total magnitudes are needed for all the filters. Following Laigle et al. 2016, instead of using AUTO magnitudes in NIR filters, we rescaled APER3 magnitudes to reproduce the total flux normalization for each object since these aperture magnitudes provide better results in SED fitting for photometric redshifts. We computed total magnitudes in optical and NIR bands following the prescriptions in Eqs. (9) and (10) of Laigle et al. 2016. In particular, we applied to each object a photometric offset that normalizes the aperture magnitudes to the total ones, preserving the colors that have been better determined in aperture photometry. We also corrected magnitudes for foreground galactic extinction using the values given in the catalog, which are computed at each object position using the Schlegel et al. (1998) dust maps. Finally, we applied corrections to both NIR and IRAC magnitudes using the small systematic offsets given in the catalog, estimated by matching the predicted magnitudes and the observed ones using the spectroscopic sample.

A census of the photometric bands adopted in this paper and the depth of observations at 2σ and 3σ are presented in Table 1. The 3σ values are provided in Laigle et al. (2016), while 2σ, that will be used in the following, have been evaluated as the magnitude of objects with photometric error Δm ∼ 0.54 (corresponding to S/N = 2).

3.2. Parent sample

In this work we focus on massive galaxies and in order to build a mass-selected sample, we used as a proxy NIR apparent magnitudes in filter IRAC[4.5]. In fact, since stars with small masses (M ≲ 1 M_⊙), which make up the bulk of the mass of a galaxy, mostly emit in the redder part of the spectrum, a selection in NIR magnitudes roughly corresponds to a selection in stellar mass. Moreover, as shown e.g., in Bell & de Jong (2001), the stellar mass-to-light ratios in the rest-frame NIR bands vary only by a factor of 2 or less over a wide range of star-forming histories, in contrast with a factor of 10 in blue bands. This means that the luminosity of a galaxy in the NIR is a good tracer of its stellar mass.

Using population synthesis models, we derived the relations between stellar mass and magnitudes in filter IRAC[4.5] as a function of redshift: at z = 3 the expected magnitude of a passive galaxy with mass M = 10^10.6 M_⊙ is m_[4.5] ≃ 22.3, that is, passive galaxies brighter than this limit should be characterized by stellar masses larger than 10^10.6 M_⊙. The same type of object with the same stellar mass can be observed with a magnitude of m_[4.5] ≃ 22.6 at z = 4 and m_[4.5] ≃ 22.3 at z = 5. Assuming a different formation redshift z_form = 8, the expected magnitude of a passive galaxy with M = 10^10.6 M_⊙ at z = 3 is m_[4.5] ≃ 22.4, while considering models with low metallicity (Z = 0.008) or a longer timescale of star formation (τ = 0.3 Gyr), the corresponding magnitude is m_[4.5] ≃ 22.2.

Based on previous considerations, we first defined and selected a parent sample with a cut at m_[4.5] ≤ 24 to guarantee both the sensitivity to stellar masses log(M/M_⊙)≳10 for quiescent galaxies and reliable photometry. In terms of the signal-to-noise ratio, this magnitude limit ensures that more than 99% of the considered data have an S/N >  3 (also see Table 1). Moreover, we considered only sources with S/N ≥ 3 (i.e., Δm ≤ 0.36) at 4.5 μm. Spectroscopically confirmed stars and active galactic nuclei (AGN), for a total of 743 objects, have also been removed.

All the above selections were performed inside the region with the best quality photometry: FLAG_PETER, which defines good optical areas, that is, masking regions where bright stars may contaminate the photometry of nearby objects, while FLAG_COSMOS and FLAG_HJMCC define the area in COSMOS field covered by UltraVISTA. UltraVISTA data are primarily important in this work as they are the basis of our selections. The final sky area considered after applying the cited flags is 1.38 deg².

The parent sample is then composed of 212 897 objects (out of the 1 182 108 entries in the original COSMOS2015 catalog). The whole parent sample is shown in gray in Fig. 3 in our color–color diagrams.

Fig. 3. Observed color–color diagrams. Left: JKs[3.6][4.5]. Right: HKs[3.6]. Gray shaded areas represent the color selection regions defined in Sect. 2 and the evolution of the colors of a quiescent and a star-forming galaxies are plotted as a reference. Gray points represent the whole parent sample selected as in Sect. 3.2, black points show the color-selected-only quiescent candidate. Red points are “bona fide” objects, while orange points belong to the “extended” sample, both defined in Sect. 4. The mean photometric errors of the color-selected candidates (black points) are also shown. Objects that are non-detected in one or more filters in the color-selected sample are not shown in order to improve clarity given their large number. Only undetected objects belonging to extended or bona fide samples are shown with arrows indicating the lower limits of the involved color. Non-detections correspond to < 2σ.

3.3. Photometric requirements

To select quiescent candidates through color–color diagrams, we impose the following photometric requirements:

1. We defined as the detection limit in each filter of interest a conservative value of S/N = 2, corresponding to a maximum photometric error of Δm = 0.54 mag, in order to include in the color selection any possible quiescent candidate. In addition, we also considered galaxies non-detected (i.e., with S/N <  2) in some of the filters used in the diagrams; when a non-detection occured, the magnitude limit in that filter was considered from Table 1 to derive a lower or upper limit of the color.

2. To minimize the fraction of lower redshift and star-forming objects, we removed from our sample those objects detected in the u band with S/N >  2; given our parent sample cut at m_[4.5] ≤ 24, this choice corresponds to removing from the sample those objects with blue u − m_[4.5] colors, which are characteristic only of star-forming galaxies or very-low-redshift quiescent galaxies.

3. We imposed a non-detection condition at 24 μm measured by MIPS (Multiband Imaging Photometer on Spitzer) to minimize the fraction of dusty star-forming galaxies that would significantly emit at 24 μm because of the dust re-emission; a detection limit of S/N = 2 was chosen.

3.4. Candidates identification through color selections

Candidates are identified from the parent sample using color–color diagrams defined in Sect. 2. As mentioned above, we also include in the selection extremely faint quiescent objects that may not be detected in some of the filters involved in the diagram. Moreover, sources which satisfy the selection criteria but do not meet the u and/or the 24 μm non-detection conditions were removed from the sample (a total of 182 743 objects in the parent sample were rejected due to these cuts).

A total of 1047 objects (0.49% of parent sample) were selected using at least one of the two criteria (183 objects were selected by both color selections simultaneously. See Table 2 for details). In Fig. 3, we show these candidates as black dots, together with two evolutionary tracks for comparison with previous plots. Objects that are non-detected in one or more filters in the color-selected sample are not shown for the purposes of clarity given their large number. We notice that a large number of objects is located just below the selection boxes: the color errors will therefore result in a larger number of objects scattered upwards into the selection box than downwards, producing an overestimate of the number of selected objects. This bias and the other possible sources of contamination are the object of the SED fitting refinement illustrated in Sect. 4.

Figure 4 illustrates the properties of the selected objects in gray: the redshift distribution of the selected sample peaks at z_phot ≈ 3.5 and appears to be bimodal, as a result of the combination of the two adopted color diagrams. A non-negligible fraction of objects exhibits photometric redshifts in the range 2 <  z_phot <  2.5: a contamination by non-quiescent and/or z <  2.5 objects was expected. In fact, the chosen color selections (with generous selection boxes), the parent sample (with a deep cut at 24 magnitudes at 4.5 μm), and the photometric requirements (with an undetection requirement at 2σ in u and 24 μm bands) have been designed in order to include any possible quiescent candidate, given their rareness at z ≥ 2.5. We also expect a certain degree of contamination by star-forming galaxies with prominent emission lines. With the aid of the full photometric data, the SED fitting procedure is expected to be able to identify the truly quiescent objects at the highest redshifts and remove the interlopers belonging to different redshift ranges or galaxy types.

Fig. 4. Distribution of physical properties of the color selected sample (gray) of the extended (orange) and bona fide (red) sub-samples.

4. SED fitting and physical properties

The SED fitting analysis of the selected objects allows us to select the most massive (log(M_*/M_⊙) > 10.6) high redshift (z ≳ 2.5 − 3) quiescent galaxies, hereafter defined as galaxies characterized by log(sSFR [yr⁻¹]) <  − 10.5, where sSFR is the specific star formation rate (sSFR = SFR/M_*). Since there are no spectroscopic redshifts measured for the candidates selected in the previous section, to evaluate the physical properties of the candidates we assumed the photometric redshifts that have been carefully optimized by Laigle et al. (2016) in the COSMOS2015 catalog. We then carried out the SED fitting to estimate physical properties, such as the stellar mass, extinction, age, and specific star formation rate (sSFR). The exploration of the parameter space has not been designed to consider all the possible SEDs, but only to break the degeneracy between the two specific galaxy populations we expect from the color–color selection, that is, high-redshift quiescent objects and dusty star-forming contaminants.

We used the hyperzmass code (Bolzonella et al. 2000, 2010) to evaluate the best fit SED, corresponding to the minimum χ² derived from the comparison between observed and model photometry at fixed z = z_phot from COSMOS2015. From the best fit SED we derive the stellar mass, star formation rate, age, and extinction. In the following we describe the relevant parameters of the fitting procedure, which makes use of the input photometric catalog including the information of photometric redshift.

To avoid degeneracies between similarly probable solutions, we tailored the set of templates considering only the two populations we expected to be included in our selection. We used two evolving templates with solar metallicity and Chabrier IMF: an exponentially declining star formation history built with BC03 SPS, with timescale of 0.1 Gyr, representing a quiescent galaxy, and a template with constant star formation including emission lines, built with fsps (Conroy et al. 2009, 2010), appropriate for star-forming galaxies. Each template contains 221 spectra for evolving ages from t = 0 to t = 20 Gyr. No formation redshift was imposed but only ages smaller than the age of the universe at the photometric redshift were considered. Other population synthesis models are available from the literature (e.g., M05); different choices could lead to variation in the stellar masses of 0.10 − 0.15 dex depending on the considered SPS model (e.g., Walcher et al. 2011; Ilbert et al. 2013). For consistency and continuity with previous works (e.g., Laigle et al. 2016; Ilbert et al. 2013; Davidzon et al. 2017; Muzzin et al. 2013), BC03+fsps models were adopted as reference templates here.

Moreover, we performed several SED-fitting runs using different extinction laws: in particular, those of Calzetti et al. (2000), which are characteristic of starburst galaxies, Fitzpatrick (1986) for the Large Magellanic Cloud, and Seaton (1979) for the Milky Way (see Appendix A.2 for the expected differences in the color selection). We chose as our reference the results obtained using the law of Calzetti et al. (2000). We set the range of extinction between A_V  =  0.0 and A_V  =  5.0 both for the τ model and for star forming galaxies. Although high values of extinction are not expected in quiescent galaxies, we let the parameter assume all the values in the range to avoid any bias on the age, and consequently, on the stellar mass because of the well-known degeneracy among these parameters (e.g., Papovich et al. 2001). As a result of making this choice, we expected to select galaxies that would be more consistent with being quiescent than dusty star-forming.

The filters used for the fit are u, B, V, r, i⁺, z⁺⁺,Y, J, H, K_s, IRAC[3.6] and IRAC[4.5]. The fluxes in the IRAC[5.8] and IRAC[8.0] bands have not been included in the final fit because of the shallower depth and worse PSF (point spread function) of these bands, which makes them much less suitable for constraining the SED shape. The filters used in the fit sample the optical/NIR wavelengths well enough even at the high redshifts considered in this work, so that the stellar mass can still be reliably determined (see e.g., Taylor et al. 2011; Pacifici et al. 2012; Mitchell et al. 2013; Conroy 2013). To avoid the problem of dealing with non-detections that should require a statistical treatment as in Sawicki (2012), we directly used fluxes (corrected to total) with their face values and errors. In this way, the χ² can be estimated also for very small or negative values of the measured fluxes, as explained in Laigle et al. (2016), and the number of data points used to constrain the fit is always ≥9.

Moreover, the intergalactic medium average opacity of intervening systems along the line of sight of high redshift objects has been taken into account following Madau (1995).

4.1. Results of SED fitting procedure

We divided our color-selected sample in two classes (the extended and the bona fide samples), following the results of the SED fitting.

The extended sample is defined to have z_phot  ≥  2.5, log(M_*/M_⊙)  ≥  10.6, log(sSFR [yr⁻¹])≤ − 10.5 and the reduced χ² of the best fit with the model of quiescent galaxies smaller than the one using the model with constant star formation (${\chi }_{\mathrm{q}}^2<{\chi }_{\mathrm{sf}}^2$).

The bona fide sample, characterized by more restrictive criteria, i.e., z_phot ≥ 3.0, log(M_*/M_⊙)≥10.6, log(sSFR [yr⁻¹]) ≤ −11, ${\chi }_{\mathrm{sf}}^2-{\chi }_{\mathrm{q}}^2>2$, and ${\chi }_{\mathrm{q}}^2<2$.

In particular, we choose the stellar mass cut at log(M_*/M_⊙)≥10.6 in order to guarantee the completeness of our samples and facilitate the comparison with results and mass functions that have previously been published in the literature. We present the results of all the selections in Table 2. The total number of objects in the bona fide sample is nine. The extended sample is made up of 128 objects. By using the different extinction laws mentioned in the previous section (i.e., Fitzpatrick 1986; Seaton 1979), we obtain almost identical results. In further detail, we find that eight objects of the bona fide sample were selected independently from the extinction law adopted in the fit, while the object left out from our reference sample is the one with the largest amount of dust (A_V = 0.9) in Fig. 4. This object is fitted with a value of A_V = 1.0 by Fitzpatrick and Seaton extinction laws, but is excluded from the bona fide sample due to the higher sSFR (log(sSFR [yr⁻¹]) ≈ −10.9 both for Fitzpatrck and Seaton laws). Also, the extended sample shows little variation compared to our reference results derived using Calzetti et al. (2000) extinction law: we select 131 objects assuming Fitzpatrick’s law, and 137 with Seaton’s law. In the color–color diagrams illustrated in Fig. 3, together with the parent sample in gray and candidates selected in Sect. 3.3 in black, we show the bona fide and extended samples with red and orange points, respectively. Red and orange points below the selection regions correspond to objects that are non-detected in one filter (H or J): they have been fully considered in the samples as their color is consistent with the selection criteria when considering it as a lower limit.

We present in Fig. 5, the SED fitting results for all the galaxies in the bona fide sample, showing both the best-fit quiescent template, and the best-fit using only the star-forming template. The star-forming model at the same redshift is not able to reproduce the break strength located between the bands H and K_s, nor the faint fluxes at blue/near UV rest frame wavelengths. Similar results for the fit with the two classes of SEDs are valid for all the other objects: on one side, the inclusion of emission lines is not sufficient to mimic the D4000 break, while on the other, the fluxes blue-ward of the break are too high despite the large values allowed for dust extinction.

Fig. 5. Best-fit models for all bona fide objects: colored points represent photometric points with error-bars, red lines are best-fit SEDs. The best fit using the star-forming template is shown with a blue line. Error bars in the wavelength scale refer to the width of the considered filter. Physical parameters evaluated through the SED fitting are shown in inserted labels. Black tilted crosses and stars show the synthetic photometry, evaluated by integrating the best-fit templates in the several bands, for the best-fit quiescent and star-forming templates respectively. All filters used in the fit are listed in inserted legends for each plot.

Figure 4 summarizes the main physical properties of the extended and bona fide samples derived through the SED fitting analysis, taking as a reference fit the one using Calzetti’s extinction law. In particular, the bona fide sample is shown in red while the extended sample is in orange. The properties shown are the photometric redshift and the observed magnitude at 4.5 μm from the COSMOS2015 catalog, the stellar mass, the sSFR, the age, and the dust extinction from the SED fitting analysis. Many of the extended sample objects are located at 2.5 ≲ z ≲ 3.0, whereas for the bona fide sample, the median redshift is z ≈ 3.38. Compared to the total sample selected in Sect. 3.3 (shown in gray), the distributions of m_[4.5] peak at brighter magnitudes. The other panels of the figure show how the bona fide and extended samples represent the most massive and quiescent galaxies chosen through color selections. We notice that three galaxies of the bona fide sample have best fit with 0.5 ≤ A_V ≤ 0.9, which are extinction values approaching those of Lyman-break galaxies at the same redshifts (Shapley et al. 2005). However, recent works (Gobat et al. 2018; Martis et al. 2019) claim that quiescent galaxies at high redshifts (1 <  z <  4) can contain at least two orders of magnitude more dust at a fixed stellar mass compared with local early-type galaxies (ETGs). They found that these dusty (A_V ≥ 1.0) high-z quiescent galaxies can comprise up to ∼20 − 25% of the population of quiescent galaxies, which is consistent with our findings. Moreover, we stress that the results we find are dependent on the population synthesis model we adopt. In fact, if we adopt a different SPS model (e.g., Maraston 2005), the best-fit A_V value may change, along with other physical parameters degenerate with dust extinction, such as the age of the galaxy. We explore this possibility, using the SPS model of Maraston (2005) in the following section.

Out of our 128 quiescent candidates (extended sample) we find that 98 of them are also classified as quiescent in the COSMOS2015 catalog, where quiescent galaxies are identified using the locations of galaxies in the color–color plane NUV − r/r − J (Williams et al. 2009).

4.2. M05 vs. BC03 models

We derive the best-fit physical parameters also using M05 models (Maraston 2005) to evaluate the robustness of our selections. The two models (i.e., BC03 and M05) differ in the treatment of the thermally-pulsing asymptotic giant branch (TP-AGB) phase. The contribution of the TP-AGB stars to the integrated light of a synthetic stellar population critically depends on what is adopted for the stellar mass loss during this phase. The higher the mass loss, the sooner the star loses its envelope and the sooner the TP-AGB phase is terminated. In M05 the TP-AGB phase contribution is much higher than in BC03 models. The result is that the M05 models are brighter and redder than the BC03 models for ages between ∼0.2 and ∼2 Gyr at λ >  2 − 2.5 μm (Maraston et al. 2006). This implies that the M05 models give, in general, younger ages and lower stellar masses compared to BC03 models. At older ages, this tendency is reversed.

In order to perform a fair comparison we apply the same criteria used to define the bona fide and extended samples described at the beginning of the section. While we adopt for BC03 models a Chabrier IMF for consistency and continuity with other works in the literature (e.g., Ilbert et al. 2013; Laigle et al. 2016; Muzzin et al. 2013; Davidzon et al. 2017), for M05 models, a Chabrier IMF is not available, so we adopt a Kroupa IMF (Kroupa 2001) and apply a statistical offset to stellar masses, such as log M_{*, Chabrier} = log M_*, Kroupa − 0.04, to take into account the different IMF. All other parameters are consistent between the two models. Moreover, for the M05 SED fitting run we also adopted the star-forming template that includes emission lines built with fsps.

Using M05 models to fit the data, 20% fewer objects are included in the bona fide sample with respect to BC03 results (i.e., eight objects are selected in the M05 bona fide sample), while ≈33% fewer objects are included in the extended sample. The case of having different numbers of objects is expected since M05 models predict for ages between ∼0.2 and ∼2 Gyr (i.e., the typical ages for galaxies at z  ≳  3) a higher flux per unit mass with respect to BC03 models (Maraston et al. 2006) at λ ≳ 2 μm, due to a higher contribution from stars in the TP-AGB phase. This reflects in a lower normalization needed for the theoretical SED to fit the observed data, and ultimately to a lower mass estimate. As a consequence, fewer objects satisfy the log(M_*/M _⊙ )≥10.6 condition.

To have a fair comparison between the number of objects selected with BC03 and M05 models, we compute the median offset of the mass estimation between BC03 and M05 models of the bona fide and extended samples. We find that the median value of log(M_*, BC03/M_*, M05) is ≈0.24 when considering the extended sample, which is larger than the value of 0.14 used by Henriques et al. (2015). However, the value adopted by Henriques et al. (2015) was derived by Domínguez Sánchez et al. (2011) and Pozzetti et al. (2010) for a mix of different galaxy populations, while the difference is expected to be larger for objects that quenched their star formation. Assuming a lower mass cut of log(M_*/M_⊙)≥10.6 − 0.24 = 10.36, the total number of objects in bona fide and extended samples is more in agreement with the results obtained using BC03 models, which are, respectively, 8% smaller and 33% bigger than the BC03 results.

We also investigated the dependence of the derived age on the adopted model. If we leave the dust extinction free to change in the process of SED-fitting, then similar ages are obtained either with M05 or BC03 models (both for the extended and bona fide samples). However, the BC03 model predicts a mean extinction value of 0.81 mag for the extended sample and ⟨A_V⟩ = 0.31 for the bona fide, while M05 models predict a mean extinction of 0.17 and 0.09 mag for the extended and bona fide samples, respectively. If, instead, we do not allow for dust extinction (i.e., we fix A_V = 0), the mean values of the ages predicted by the two models differ, with M05 models predicting ≈0.4 Gyr younger ages with respect to BC03 models. These results are in agreement with the work presented in Maraston et al. (2006), where similar conclusions are found for a sample of seven passive galaxies at z >  2 when dust extinction is forced to assume only very small values.

There is no consensus, yet, in the literature about which of the two models better represents high redshift galaxies, in other words, about the relative importance of the TP-AGB phase in high redshift objects. In fact, while, for example, Maraston et al. (2006) finds a better fit of high-z galaxies using M05 models, other studies, such as that of Kriek et al. (2010), find better fits using BC03 models. In this paper, we do not find any strong evidence in favor or against any of the two considered SPS models. In the following of this paper, we will use BC03 results to ease the comparison with previous work (e.g., Davidzon et al. 2017; Ilbert et al. 2013; Muzzin et al. 2013).

5. Number and mass densities

We use the results obtained in the previous section to estimate the number and stellar mass densities of quiescent galaxies at z >  2.5 and compare them to observations and state-of-the-art semi-analytical models. In Tables 3 and 4, the number and mass densities of the bona fide and extended samples are presented in redshift bins and shown in Fig. 6. To evaluate the number and mass densities, we consider only massive galaxies with M ≥ 10^10.6 M_⊙, both for the extended sample (red diamonds) and for the bona fide sample (red squares).

Fig. 6. Number and stellar mass densities of galaxies with M ≥ 1010.6 M⊙ as function of redshift for total population of galaxies (in gray) or quiescent galaxies (colored points) in literature, compared to those obtained in this work. Our results are shown as red diamonds and squares for the BC03 extended and bona fide samples, respectively, with error bars representing the Poissonian errors, cosmic variance, and photometric redshift errors added in quadrature. The forecasts of two semi-analytic models are also shown: the black dashed and continuous lines represent the Guo et al. (2013) and the H15 models, respectively, estimated for the entire population of galaxies and convolved for the Eddington bias. The black dotted line represents the H15 model for the entire population not convolved for the Eddington bias. Blue and cyan shaded regions represent the H15 model with Poissonian error for the selected quiescent populations characterized by log(sSFR [yr−1])≤ − 10.5 and ≤ − 11, respectively.

The error bars we show in Fig. 6 have been evaluated by adding in quadrature Gaussian errors for the object counting (we approximated Poissonian statistic with the Gaussian one given that for large counts, the two statistics match), cosmic variance, and the scatter between different redshift bins due to photometric redshift errors. To take into account cosmic variance, we follow the prescriptions of Moster et al. (2011). In particular, we choose the field size to be consistent with the COSMOS field we are considering, we impose the same mass cut we adopt (i.e., log(M_*/M_⊙)≥10.6) and for each redshift bin we choose the median redshift of the objects grouped in that bin. We find that, for our mass range, the cosmic variance ranges from a value of σ_CV = 0.15 at 2.5 ≤ z <  3.0 to σ_CV = 0.25 from at z ≥ 4.0. To take into account the scatter of objects between different redshift bins due to photometric errors, we checked how much the number of objects in the adopted bins can change by applying a noise of σ_z ∼ 0.03(1 + z) to the redshift distribution of a sample selected with m_4.5 <  24. We used the COSMOS catalog itself and a lightcone from H15 model to estimate the variation in the number of objects in high redshift bins, and we find that it is always less than ∼7%. We adopted an additional error of 10% on our measurements of number and mass densities. We notice that the assumption on the scatter can be optimistic given the lack of spectroscopic control samples for this class of objects, but even assuming a scatter as large as σ_z ∼ 0.10(1 + z), we expect a contamination of 11% in the wide bin 3.0 ≤ z <  4.5 adopted in Sect. 6, and a maximum contamination of ∼20% at 3.5 ≤ z <  4.0 and ≲50% at 4.0 ≤ z <  4.5.

We find the following results. For the extended sample: we find a decrease by a factor of $\sim 14{}_{-7}^{+3}$ in the number of quiescent object from the lowest redshift bin (i.e., 2.5  ≤  z <  3.0) to the highest one (i.e., z  ≥  4.0), and the same decrease is found by considering the mass densities between the same redshift intervals. In addition, we find a decrease by a factor $\sim 66{}_{-27}^{+28}$ and $\sim 100{\pm }_{-35}^{+29}$ in the number density and mass density, respectively, between Davidzon et al. (2017) data for the quiescent population at 0.2 <  z <  0.5 and our estimate in the redshift bin 2.5 ≤ z <  3.0. Between the same redshift bins, considering the results of Moutard et al. (2016), which are evaluated on the VIPERS multi-lambda catalog (on 22 deg²) and are, therefore, more statistically significant and less subject to cosmic variance with respect to results on COSMOS, we find a decrease by a factor $\sim 57{}_{-19}^{+14}$ and $\sim 78{\pm }_{-19}^{+15}$ in the number density and mass density, respectively. Concerning the bona fide sample: between the lowest redshift bin (i.e., 3.0 ≤ z <  3.5) to the highest (i.e., 3.5 ≤ z <  4.0) we find an evolution by a factor of ∼1.7 in the number density and almost no evolution in the mass densities.

The errors on the estimated factors have been evaluated by propagating the errors of our measurements. When considering the Davidzon et al. (2017) data, the errors on their measurements have also been considered and propagated. Interestingly, the 10% uncertainty we add to our error budget in order to take photometric redshift errors into account, is similar to what Ilbert et al. (2013) quoted as error relative to the template fitting procedure, including photometric redshift error and stellar mass estimate uncertainties.

5.1. Comparison with previous results

In Fig. 6, we compare the evolution of number and stellar mass densities with results from the literature, specifically, Pozzetti et al. (2010), Davidzon et al. (2013, 2017), Santini et al. (2012), Ilbert et al. (2013), Muzzin et al. (2013), Mortlock et al. (2011), Fontana et al. (2006), Caputi et al. (2015), Stark et al. (2009), Song et al. (2016), and Grazian et al. (2015). The datapoints have been rescaled where necessary to the cosmology we adopted (see Sect. 1) and to a Chabrier IMF. We consider only massive galaxies with M ≥ 10^10.6 M_⊙ both for the full sample of galaxies (gray points) and for the sub-sample of quiescent galaxies (colored points). The literature we used featured data which focused on quiescent galaxies (Davidzon et al. 2017; Ilbert et al. 2013; Muzzin et al. 2013; Moutard et al. 2016). All of them, like the present work, collect a galaxy sample where photometric redshifts and stellar masses are derived via SED fitting. The classification of quiescent galaxies is based on the sSFR in the present work, whereas it is based on the NUVrJ diagram for Davidzon et al. (2017), Moutard et al. (2016), and Ilbert et al. (2013), and on the UVJ diagram for Muzzin et al. (2013). We have been able to evaluate the number and mass densities of quiescent galaxies at the highest redshifts from the stellar mass functions that have been published. Data from Davidzon et al. (2017) are estimated by means of the integral of their stellar mass function for the quiescent population, obtained from COSMOS2015 catalog. The same recipe has been applied to the stellar mass function of quiescent galaxies of Muzzin et al. (2013) selected from a NIR-selected sample of galaxies out to z ∼ 4, and to the mass functions of Moutard et al. (2016) at z ≤ 1.5 derived in a 22 deg² field. The results of the present work are in broad agreement with the results obtained by Ilbert et al. (2013), Muzzin et al. (2013), and Davidzon et al. (2017).

5.2. Comparison with semi-analytic models

We compared our results with the number and stellar mass densities obtained using semi-analytic models by Guo et al. (2013; hereafter G13) and Henriques et al. (2015; hereafter H15) for the total and the quiescent populations. Both models are based on the Munich galaxy formation model (Kauffmann et al. 1999; Springel et al. 2001; Croton 2006; De Lucia & Blaizot 2007; Guo et al. 2011, 2013; Henriques et al. 2015) which has been implemented in the Millennium (Springel et al. 2005) simulation of dark matter in a box with comoving side of 500 Mpc h⁻¹, and with cosmological parameters adopted from the Wilkinson Microwave Anisotropy Probe (Komatsu 2006) and from the Planck Collaboration I (2014) for G13 and H15, respectively. We used the data from the Millennium database², selecting galaxies directly from the snapshots of the simulation. Before applying the same selection criteria adopted for observed galaxies, we convolved stellar masses with a Gaussian in log(M_*/M_⊙), with width increasing with redshift σ_{log M*} = 0.08(1 + z), in order to account for the Eddington bias, as done in Henriques et al. (2015). We set the Hubble parameter as in our reference ΛCDM cosmology, and we selected galaxies with log(M_*/M_⊙)≥10.6 at all sSFR for the total population shown in Fig. 6. For quiescent galaxies, we selected massive galaxies (i.e., log(M_*/M_⊙)≥10.6) with two different cuts in log(sSFR [yr⁻¹]): one at −10.5 to match our extended sample definition and one at −11 to match the bona fide selection (shown in Fig. 6 as blue and cyan lines respectively). It can be argued that the different star formation histories assumed in our SED fit and in the models may bias the comparison between the values. However, it has been shown by Laigle et al. (2019) that the instantaneous or short time-scale sSFR derived for a COSMOS-like photometric sample using exponentially declining SFHs reproduces the intrinsic values of a simulation with stochastic SFHs reasonably well.

In Fig. 6, we also show the effect of the correction for the Eddington bias for the total population in the H15 model: the convolution with the Eddington bias tends to increase the number counts of galaxies especially at high redshift, since more galaxies with small masses (which are more numerous than galaxies at high masses) are scattered upwards of the mass threshold than high mass galaxies scattered downwards of the same limit. Since the observed values are naturally the result of the convolution of intrinsic properties and observational errors, it is fundamental to take this bias into account in order to make a fair comparison between models and observations. Figure 6 shows that the amount of the convolved error is what is needed to match the densities of massive observed galaxies at high redshift.

While the total number and mass densities by H15 and G13 are in good agreement with the literature data, Fig. 6 shows that the same quantities for passive galaxies derived by H15 are underestimated when compared with results from the present work and from the literature at high redshift: even considering a large scattering in photometric redshifts as mentioned above, the decrease in number and mass densities is not sufficient to fully reconcile the results with the models, at least for the extended sample, as we discuss in Sect. 7.

6. Stellar mass functions

The stellar mass function (SMF) of galaxies has been studied extensively over the past years out to z ∼ 4 − 5, for both star-forming and quiescent galaxies (e.g., Fontana et al. 2006; Ilbert et al. 2013; Davidzon et al. 2017; Muzzin et al. 2013; Grazian et al. 2015). In this section, we estimate the stellar mass function for quiescent galaxies at z >  2.5 and compare it to other observations and models.

We estimate the stellar mass function of the extended sample into two redshift bins: 2.5 <  z <  3 and 3.0 <  z <  4.5 with the 1/V_max method (considering, for all the galaxies, the V_max equal to the volume defined by the redshift bin since we assume it to be complete for masses above ∼10¹⁰ M_⊙, as shown in Fig. 2). We fit our data (3 mass bins for each redshift interval) with a Schechter function (Schechter 1976):

$$\begin{array}{c}\hfill \mathrm{\Phi }(M)\mathrm{d}M={\mathrm{\Phi }}^{\ast }{\left(\frac{M}{M^{\ast }}\right)}^{\alpha }exp(-\frac{M}{M^{\ast }})\frac{\mathrm{d}M}{M^{\ast }}\end{array}$$(3)

using the software package emcee (Foreman-Mackey et al. 2013), an MIT-licensed pure-python implementation of affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare 2010). Since we are probing only the massive end of the function, the slope of the power law part of the function is not constrained. We thus fix the parameter α = 1.15, following Davidzon et al. (2017) and assuming the same slope they find at lower redshift. We report in Table 5 the Schechter parameters fitting the data points in the two redshift bins along with the 1σ errors.

In Fig. 7, the data points along with the best-fit Schechter functions are shown. Error bars include the contribution of Poissonian errors, cosmic variance, and photometric redshift errors added in quadrature, while in the abscissa, they represent the mass bin Δlog M/M_⊙ = 0.2. For the cosmic variance, we once again used the prescriptions of Moster et al. (2011), deriving for each bin the cosmic variance on the COSMOS field, at the median redshift of the objects grouped in each bin and in the mass range covered by the bin. To account for photometric redshift errors, we considered a 10% error, as explained in Sect. 5.

Fig. 7. Stellar mass function of extended galaxies, in two redshift bins between z = 2.5 and 4.5. In each panel, the data points are shown as red circles in bins of Δlog M/M⊙ = 0.1. Error bars include Poisson noise, cosmic variance, and photometric redshift errors. The data points are fitted by a single Schechter function, shown by a red solid line, while the red shaded area is its 1σ uncertainty. Data points (with their error bars) from other works are also shown: violet squares are taken from Davidzon et al. (2017), in magenta triangles we plot the mass functions by Ilbert et al. (2013), and in orange circles the one by Muzzin et al. (2013). The blue lines shows the H15 model for quiescent galaxies at z = 2.68 (left) and z = 3.44 (right).

We compared our SMF with the literature data in the same redshift bins (rescaling masses to Chabrier IMF and ΛCDM cosmology when required). We selected previous works that focused on quiescent galaxies (Davidzon et al. 2017; Ilbert et al. 2013; Muzzin et al. 2013) and whose characteristics are explained in Sect. 5.1. In Fig. 7, we also show in blue the SMF for the H15 model (described in Sect. 5.2) that is derived in the snapshot located at z = 2.68 for the redshift bin 2.5 ≤ z <  3.0 and at z = 3.44 for the redshift bin 3.0 ≤ z <  4.5. We chose the snapshots located at redshifts that were the nearest to the median redshift of our observed galaxies in the two bins, which are z = 2.60 and z = 3.35. We selected quiescent galaxies in the snapshots in the same way we selected the extended sample, that is, log(M_*/M_⊙)≥10.6 and log(sSFR [yr⁻¹])≤ − 10.5, after convolving the masses with a Gaussian of width 0.08(1 + z) in log(M_*). In Sect. 7, we present an extensive discussion of our results and a comparison with results in the literature and semi-analytic models.

7. Discussion

7.1. Comparison with previous observations

In Fig. 6, we compare our estimates of the number and mass densities with several results presented in the literature. The observed densities of the quiescent population are characterized by a rather small scatter: our results are in good agreement with Davidzon et al. (2017) at all redshifts, while larger differences are visible with Muzzin et al. (2013) and Ilbert et al. (2013), especially at z <  3. This can be due to different input photometric datasets, SED fitting analysis performed on the data, and the different criteria adopted to define a quiescent galaxy (see Sect. 5.1).

From the comparison of the observed SMFs in Fig. 7, it is evident that at z >  3.0 Muzzin et al. (2013) estimate a greater number of very massive (> 10¹¹ M_⊙) quiescent galaxies than all the other works. The estimate by Ilbert et al. (2013) is, instead, only slightly higher than ours at 2.5 <  z <  3, while we obtain the best agreement with Davidzon et al. (2017) in both redshift bins. Since all the SMFs are derived in the same field, differences cannot be simply ascribed to observations on different fields.

Based on Fig. 7, there is a noticeable evolution between the two redshift bins. This effect can be due to the fast evolution of colors. We note once again that our classification is based on apparent color selections designed to identify the Balmer and D4000 breaks. A passively evolving stellar population develops pronounced Balmer and D4000 breaks in t ≈ 0.3 Gyr (Bica et al. 1994), while the redshift interval from z = 4.5 to z = 2.5 corresponds to 1.2 Gyr. This means that a galaxy that quenched its star formation at z ≈ 4 − 4.5, and since then has been passively evolving, had enough time to develop the pronounced Balmer and D4000 breaks needed to be identified as quiescent at z ≈ 2.5 − 3 through our color selections. Therefore, if the quenching of the star formation occurs at around z ≈ 4 − 4.5 for a large number of objects, it is reasonable to expect such a fast evolution in the number of passive galaxies. Another possible explanation for the evolution between the two redshift bins is the evolution of the stellar mass: the rapid increase of objects with stellar masses M_* >  10^10.6 M_⊙ may reflect the emergence of a high merging rate, assembling more and more massive quiescent galaxies as they move down with redshift, although the lack of evolution in M^* of the Schechter function can challenge this interpretation.

7.2. Comparison with models

According to our results and other results in the literature, the H15 model tends to underpredict quiescent galaxies at z ≳ 2.5, both in number and in the content of their stellar mass, as visible in Fig. 6. In particular, when considering the bona fide sample, we find that the H15 model underpredicts the number of quiescent objects in the redshift bin 3.0 ≤ z <  3.5 by a factor of $\sim 4.9{}_{-2.4}^{+2.3}$ up to a value of $\sim 6.0{}_{-3.5}^{+3.1}$ at z ∼ 4. By checking the differences in the mass densities, we find that the model underpredicts observed values at 3.0 ≤ z <  3.5 by a factor of $\sim 4.2{}_{-2.4}^{+2.0}$ and it shows a difference by a factor of $\sim 6.2{}_{-2.7}^{+2.5}$ at z ∼ 4. The differences between H15 model and our results are all the more evident when considering the extended sample. In this case, in the redshift bin 2.5 ≤ z <  3.0 we find a difference by a factor of $\sim 11.9{}_{-2.9}^{+2.7}$ in the number densities, while the mass densities differ by a factor of $\sim 11.1{}_{-1.2}^{+1.0}$. Moreover, at z ∼ 4 we find a difference by a factor of $\sim 10.1{}_{-6.1}^{+5.8}$ in the number densities and by a factor of $\sim 19.2{}_{-4.1}^{+3.9}$ when considering the mass densities (also in this case, we simply propagated the errors in our measurements).

Considering the SMFs in the left panel of Fig. 7 (i.e., at 2.5 ≤ z <  3.0), the H15 model appears not to reproduce the shape of the observed SMFs. In particular, it emerges that the model overpredicts low-mass (log(M_*/M_⊙)≲10.4) and high-mass (log(M_*/M_⊙)≳11.4) systems while objects in the intermediate mass range appear to be underpredicted (by a factor of ∼8.7 at log(M_*/M_⊙) = 10.9). Also, in the right panel of Fig. 6 (i.e., at 3.0 ≤ z <  4.5), it is evident that the model seems to under-predict objects in the whole mass range (by a factor of ∼4.4 at log(M_*/M_⊙) = 10.9) with the exception of low-mass systems, i.e., log(M_*/M_⊙)≲10.2. These discrepancies reflect the intrinsic difficulty in treating the processes involved in galaxy formation, and in particular, the processes related to the transformation of star-forming galaxies into quiescent objects and their mass assembly at z ≳ 2.5. In Cecchi et al. (2019) we compare high-redshift observed quiescent galaxies to different semi-analytic models (SAMs) and discuss the quenching mechanisms that led to their formation.

In the current scenario, such massive quiescent galaxies are the result of a strong active galactic nuclei (AGN) feedback from the central supermassive black holes (BHs) which significantly affects galaxy formation processes. The AGN feedback can take place in two main ways: the jet (or radio) mode, and the radiative mode (sometimes called quasar mode or bright mode). The mechanism that is primarily, or entirely, responsible for quenching in the jet mode is connected with highly collimated jets of relativistic particles, in which star formation dies out because the hot gas halo is continually heated, and the supply of new cold gas is cut off (Bower et al. 2006; Croton 2006; Somerville et al. 2008; Kimm et al. 2009). In the radiative mode, quenching is associated with mergers and rapid BH growth, followed by a quasar wind which expels gas from galaxy’s center, explaining the growth of the quiescent population (Hopkins et al. 2008a,b). Cosmological zoom-in simulations, which include fast momentum-driven AGN winds, also appear to be able to quench star formation. Moreover, it appears that it can also maintain quiescence over long timescales without any explicit jet mode type feedback (Choi et al. 2014). Conversely, Gabor & Davé (2015) suggested that the presence of a gaseous halo kept hot by AGN feedback is sufficient to quench a galaxy without the need for additional radiative mode feedback. What the relative importance is of these AGN feedback mechanisms that are responsible for the early appearance of the population of quiescent galaxies at high redshift is still a matter of debate. However, a feedback mechanism is fundamental in order to quench star formation and form massive quiescent galaxies.

In the H15 model, massive galaxies can be quenched by AGN feedback depending on black-hole and hot-gas mass and, therefore, indirectly on stellar mass. In addition, galaxies of any mass can be quenched by ram-pressure or tidal stripping of gas and through the suppression of gaseous infall. It is argued by Henriques et al. (2017) that this combination of processes produces quenching efficiencies which depend on stellar mass, host halo mass, environment density, distance to group centre, and group central galaxy properties. In the case of massive galaxies, the quenching is likely due to AGN feedback. As explained in Henriques et al. (2017), both the quasar and the radio accretion modes on the black hole are considered in the model; while the quasar mode produces no feedback on the galaxy, the radio mode produces a strong feedback by avoiding further hot gas condensation and, therefore, star formation is suppressed once the cold gas is exhausted.

If our result (that is, the under-prediction of quiescent objects at 2.5 ≲ z ≲ 3 in H15) is confirmed, it would imply that the considered model is not efficient enough in producing massive quiescent galaxies at high redshifts. Considering the ingredients of the model, this may be due to the parametrization of the radio mode AGN feedback on the galaxy, or the timescale of this feedback process is too long to produce massive quiescent objects at such redshifts or, in addition, a different accretion mechanism on the black hole should be considered.

On the observations side, the scatter between different observed samples of quiescent objects makes it evident that there is a need for deeper observations on different fields (all the data reported in the plot for quiescent galaxies belong to the COSMOS field) and spectroscopic confirmation to better constrain the population of massive quiescent galaxies at z ≳ 2.5 − 3.0, whose presence should be explained by models of galaxy formation but still aren’t (see Somerville & Davé 2015; Naab & Ostriker 2017 for two extensive reviews on theoretical state-of-the-art models of galaxy formation).

8. Summary and conclusions

The results of the present study can be summarized as follows:

1. We identified two new color selections using NIR bands to select quiescent galaxies at 2.5 ≲ z ≲ 4.5. The color selections are based on the identification of strong spectral features characterizing these evolved objects, that is, the D4000 and Balmer breaks.

2. We studied the effectiveness of our color selections by exploring all the parameters characterizing evolutionary tracks, as well as the effects of emission lines.

3. We applied the color selection to the COSMOS2015 catalog, selecting a parent sample with a cut at m_[4.5] ≤ 24. Through a SED fitting analysis we tightened the selection considering only the most massive (log(M_*/M_⊙) > 10.6) high redshift (z_phot ≳ 2.5) quiescent (log(sSFR [yr⁻¹]) <  − 10.5) galaxies. The proposed color selections coupled with SED fitting analysis allowed us to build a reliable sample of quiescent candidates, aimed at maximizing its completeness.

4. The objects consistent with being the most quiescent massive galaxies (i.e., the bona fide sample) have observed number densities that decrease by a factor of ∼1.7 from the redshift bin 3.0 ≤ z <  3.5 to the bin z ≥ 4 while mass densities show almost no evolution with redshift. Considering the extended sample, we found a decrease by a factor of $\sim 12.4_{-7}^{+3}$ in the number of quiescent object and by a factor of $\sim 10.9_{-6}^{+3}$ in their mass densities from the lowest redshift bin (i.e., 2.5 ≤ z <  3.0) to the highest one (i.e., z ≥ 4.0). In addition, we find a difference in the number density (mass density) by a factor of $\sim 66{}_{-27}^{+28}$ ($\sim 100{}_{-35}^{+29}$) between our lowest redshift bin (i.e., 2.5 ≤ z <  3.0) and Davidzon et al. (2017) data for the quiescent population in the redshift bin 0.2 <  z <  0.5.

5. We estimated the stellar mass functions for our sample of quiescent galaxies and their fit with a Schechter function.

6. According to our results, the semi-analytical model by Henriques et al. (2015) is not able to fully account for the number and mass densities and the stellar mass functions of candidate quiescent objects unless a severe effect of contamination is what is affecting our bona fide sample, in which case the significance of the disagreement with observed data can decrease. In particular, considering the bona fide sample, the model seems to underpredict the number of quiescent objects in the redshift bin 3.0 ≤ z <  3.5 by a factor of $\sim 4.9{}_{-2.3}^{+2.4}$ that grows to a value of $\sim 6.0{}_{-3.5}^{+3.1}$ at z ≥ 4.0. In considering the mass densities, we find that the model underpredicts our estimates in the redshift bin 3.0 ≤ z <  3.5 by a factor of $\sim 4.2{}_{-2.0}^{+2.4}$ and it shows a difference by a factor of $\sim 6.2{}_{-2.5}^{+2.7}$ at z ≥ 4.0. The differences between H15 model and our results is notably evident when considering the extended sample. In this case, the redshift bin 2.5 ≤ z <  3.0 demonstrates a difference by a factor of $\sim 11.9{}_{-2.9}^{+2.7}$ in the number densities, while the mass densities differ by a factor of $\sim 11.1{}_{-1.2}^{+1.0}$. Moreover, at z ≥ 4.0 we find a difference by a factor of $\sim 10.1{}_{-6.1}^{+5.8}$ in the number densities and by a factor of $\sim 19.2{}_{-4.1}^{+3.9}$ considering the mass densities. Moreover, the shape of the SMF is not fully reproduced by the model.

7. This method can be used to select targets for spectroscopic follow-up, especially with future facilities such as the James Webb Space Telescope. If future spectroscopic observations confirm the presence of a population of quiescent galaxies at high redshift, as proposed in the present work, some further internal mechanisms of quenching will be needed to explain the presence of such galaxies at high redshift that are becoming more and more numerous, given the many studies devoted to this topic. In semi-analytic models, the AGN-feedback needed to quench star-formation at high stellar masses could be an ingredient that still needs to be fully understood and parametrized differently.

Acknowledgments

We thank Bruno Henriques, Gianni Zamorani and Lucia Pozzetti for fruitful discussions on the topic. Data compilations the studies used in this paper were made much more accurate and efficient by the online WEBPLOTDIGITALIZER code. We acknowledge the support from the grants PRIN-MIUR 2015 and ASI n.2018-23-HH.0.

References

Appendix A: Exploring evolutionary tracks’s parameters space

We explored the effects of different evolutionary tracks on the choice of color criteria. We analyzed, in particular, the impact of considering star forming galaxies with emission lines, different extinction laws, IMFs, stellar population synthesis (SPS) models, metallicities, SFHs, and redshifts of formation. In all following plots: gray shaded areas represent the selection region for quiescent galaxies, colored numbers represent the redshift of the nearby point for the evolutionary track with the same color and a vector corresponding to a magnitude extinction of A_V = 1 using Calzetti’s law is shown.

A.1. Emission lines

We built the evolution of a galaxy template of constant star formation with the contribution of nebular emission lines using the code fsps (Conroy et al. 2009, 2010). As shown in Fig. A.1, considering a template with emission lines in the HK_s[3.6] diagram does not affect the selection criteria. Instead, considering the JK_s[3.6][4.5] selection, some contamination by star forming galaxies with emission lines is expected at 2.2 ≲ z ≲ 3.0 for extinction values around A_V ∼ 3. In Fig. A.1, it can be seen that higher (or lower) values of A_V would not affect the selection box contaminating the candidates. Considering the SED fitting procedure for a star-forming template with emission lines, this contamination may be kept under control with the aid of optical bands where the emission-line galaxies should be characterized by a larger flux with respect to a quiescent object (see also the top panels of Fig. A.1).

Fig. A.1. Effect of emission lines on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies (magenta and cyan lines), representing a population 1 Gyr old redshifted to z = 3 and a population 0.5 Gyr old redshifted to z = 4 (i.e., zform ≈ 6) built with BC03 models. A star forming SED with emission lines built from fsps, 1 Gyr old and redshifted to z = 2.5 is shown in dark-cyan. Also the filters transmission curves, used in COSMOS field and in the tracks computation, are shown. Bottom: evolution in color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue, black dark-cyan, and cyan, and their characteristics are shown in the legend inside the plot. A track for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) is shown in magenta.

A.2. Extinction laws

We also tested our color–color diagrams using three different extinction laws. In particular, as shown in Fig. A.2, we computed the evolutionary tracks using the Calzetti law (Calzetti et al. 2000), characteristic of starburst (SB) galaxies, the Fitzpatrick law (Fitzpatrick 1986) derived for the Large Magellanic Cloud (LMC), and the Seaton law (Seaton 1979) obtained for the Milky Way (MW). The main differences among them are the bump at 2200 Å, which is absent in starburst galaxies, and the slope in the UV, which is steeper for the LMC extinction law.

Fig. A.2. Effect of different extinction laws on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies (magenta and cyan lines), that is, a population 1 Gyr old redshifted to z = 3 and a population 0.5 Gyr old redshifted to z = 4 (i.e., zform ≈ 6) built with BC03 models. Bottom: evolution in the color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in solid blue, solid black and dashed cyan lines, for the Calzetti, Seaton and Fitzpatrick laws, respectively, and their characteristics are shown in the inserted legends. Tracks for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) are shown in solid magenta, solid orange, dashed red and brown for the Calzetti (AV = 0 and AV = 1), Seaton and Fitzpatrick laws, respectively.

As shown in Fig. A.2, tracks built with MW and LMC laws are, in general, very similar and differ slightly from the SB law we used as a reference in the present paper. Adopting different extinction laws does not affect the selection criteria of quiescent galaxies. The tracks of color evolution for star-forming galaxies can cross our selection boxes in the JK_s[3.6][4.5] diagram when adopting LMC and MW extinction laws, but only at z ≲ 3.

A.3. IMF

We tested our selections with the IMFs of Salpeter (1955), Chabrier (2003), and Kroupa (2001) and we found no appreciable differences in evolutionary tracks. This was expected because the SED of a passively evolving galaxy is mainly determined by the stars at the turnoff and considering the possible ages of a galaxy at z ≳ 3, the stars dominating the emission, and, therefore, the colors, are in the portion of the IMF that is similar for all the mentioned IMFs.

A.4. SPS models

The difference between BC03 and M05 models (Bruzual & Charlot 2003; Maraston 2005) for the building of evolutionary tracks is explored in Fig. A.3. As explained in Sect. 4.2, for λ >  2 − 2.5 μm and ages between ∼0.2 and ∼2 Gyr, a notable difference in flux between the two models is expected (Maraston et al. 2006) which means that a difference in K_s − [3.6] color is also expected. However, both models generate evolutionary tracks for quiescent galaxies which fall in the selection region for approximately the same redshift range, and leave tracks for star-forming galaxies outside of it. Therefore, the choice of the model adopted does not affect the color selection criteria. The difference between BC03 and M05 results in SED fitting for observed galaxies is also explored in Sect. 4.2, where no appreciable difference has been found.

Fig. A.3. Effect of the choice of SPS models on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies, that is, a population 1 Gyr old redshifted to z = 3 (i.e. zform ≈ 6) built with BC03 models (magenta line) and with M05 models (orange line). Bottom: evolution in color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan, while their characteristics are shown in the inserted legends. Tracks for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) are shown in magenta and orange for BC03 and M05, respectively.

A.5. Metallicity

The effect of different metallicities is explored in Fig. A.4 (for Z = 0.02,  0.008,  0.0004): we found that all the derived tracks enter the color selection for z ≈ 3 − 3.5 and exit for z ≈ 4 − 4.5, suggesting that the color selection is valid for a wide range of metallicities. Assuming a metallicity in the range of Z = Z_⊙ to 0.2 Z_⊙ is reasonable even for high redshift galaxies (e.g., Sawicki & Yee 1998; Papovich et al. 2001; Shapley et al. 2001; Mannucci et al. 2010; Finkelstein 2016). For instance, according to Glazebrook et al. (2017), a quiescent galaxy at z = 3.7 is consistent with having a solar metallicity.

Fig. A.4. Effect of metallicity on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown. In particular, the magenta SED represents a population 1 Gyr old redshifted to z = 3 (i.e. zform ≈ 6) with Z = Z⊙, while the green line shows the SED of a population of a 1 Gyr redshifted to z = 3 (i.e. again zform ≈ 6) with Z = 0.2 Z⊙. Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan, while their characteristics are shown in the inserted legends. Three tracks for quiescent galaxies are shown in magenta, dark-green and light-green for different metallicities (Z = Z⊙, Z = 0.4 Z⊙, Z = 0.2 Z⊙ respectively).

Fig. A.5. Effect of SFHs on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown in magenta SED and cyan, representing a population of 1 Gyr redshifted to z = 3 (i.e. zform ≈ 6), and a population of 0.5 Gyr redshifted to z = 4 (i.e., again zform ≈ 6), respectively. Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan and their characteristics are shown in the inserted legends. Three tracks for quiescent galaxies are shown in magenta, orange and red: the magenta track is parametrized with an exponentially declining star formation with τ = 0.1 Gyr, orange track has τ = 0.3 Gyr, and red track has been parametrized with a delayed star formation history where τ = 0.1 Gyr.

Fig. A.6. Effect of formation redshift on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown. In particular, magenta SED represents a population of 1 Gyr redshifted to z = 3 (i.e. zform ≈ 6), while cyan line shows SED for a population of 1.5 Gyr redshifted to z = 3 (i.e. zform ≈ 8). Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan and their characteristics are shown in the inserted legends. Two tracks for quiescent galaxies with redshift of formation zform = 6 and zform = 8 are shown in magenta and orange respectively. Their e-folding time is τ = 0.1 Gyr and they have solar metallicity.

A.6. Star formation histories

As already mentioned, star formation histories for quiescent galaxies have been parametrized with an exponentially declining star formation with short e-folding times in order to guarantee a negligible SFR after a few hundreds Myr. Star formation histories follow the formula SFR(t)  ∝  τ⁻¹e^−t/τ where τ is the e-folding time. High-redshift quiescent galaxies need to have formed their stellar mass quickly since they are observed at z ≥ 3 when the universe was only 2 Gyr old; therefore, a timescale of τ  =  0.1 Gyr is often assumed. By considering different e-folding times and assuming the same redshift of formation z_form  =  6, models with a longer e-folding time will be characterised by bluer colors with respect to the model, with a shorter τ at the same redshift. The color selection is still valid assuming an e-folding time of 0.3 Gyr, even though at fixed ages (or redshift) tracks with longer e-folding times are bluer. Assuming even larger e-folding times will not allow for the selection of passive galaxies at high redshift since the stellar population do not have enough time to evolve and develop the D4000 and Balmer breaks.

Delayed star formation histories, parametrized as ∝τ⁻²t e^−t/τ, have also been explored. In particular, choosing timescales of the SFH that are small compared to the time between z_form and z_obs (e.g., τ ≈ 0.1 Gyr), the color selections are still valid, although on a slightly narrower redshift interval. Longer e-folding times would instead represent star-forming galaxies in the redshift interval z ≈ 2.5 − 4.

A.7. Redshift of formation

Models with higher z_form are, in general, characterized by redder colors at the same redshift since the galaxy which forms at a higher z will exhibit prominent breaks at earlier epochs. This leads to a slightly wider redshift interval for the selection (for z_form = 8 we select quiescent galaxies in the range 3 ≲ z ≲ 4.8 for HK_s[3.6] and in the range 2.5 ≲ z ≲ 4.5 for JK_s[3.6][4.5]) since galaxies enter the color selections at higher redshifts than they did with our adopted reference value of z_form = 6.

All Tables

All Figures

Fig. 1. Top panels: two SEDs of quiescent galaxies. SED in magenta represents a population of 1 Gyr redshifted to z = 3 (i.e., zform ≈ 6), while the cyan line shows the SED of a population of 0.5 Gyr redshifted to z = 4 (i.e., again zform ≈ 6). Filters transmission curves used in COSMOS field and in the tracks evaluations are also shown. Bottom panels: JKs[3.6][4.5] and HKs[3.6] color–color diagrams with different evolutionary tracks both for star-forming galaxies and passive galaxies. Tracks representing star-forming galaxies are shown in cyan, blue and black with AV = 0, AV = 3, and AV = 5 respectively. Three tracks for quiescent galaxies (with exponentially declining star formation history τ = 0.1 Gyr and solar metallicity) are shown in magenta, orange and red with different extinction values (AV = 0, AV = 1, and AV = 2 respectively). Gray shaded areas represent the selection region for quiescent galaxies at 2.5 ≲ z ≲ 4 in the left panel and at 3 ≲ z ≲ 4.5 in the right panel. The numbers near the tracks of the same color represent the redshift. A vector corresponding to a magnitude extinction of AV = 1 using Calzetti’s law is also shown.

In the text

Fig. 2. Stellar mass as a function of redshift for different theoretical parametrizations of quiescent galaxies: we adopted an exponentially declining star formation history with characteristic time τ = 0.1 Gyr, Z = Z⊙ = 0.02, AV = 0, zform = 6. Red and black lines correspond to magnitudes m[4.5] = 22.27 and 24, respectively. Orange, green, and magenta lines represent different models of a quiescent galaxy differing in redshift of formation, metallicity, and e-folding time of the SFH, respectively.

In the text

Fig. 3. Observed color–color diagrams. Left: JKs[3.6][4.5]. Right: HKs[3.6]. Gray shaded areas represent the color selection regions defined in Sect. 2 and the evolution of the colors of a quiescent and a star-forming galaxies are plotted as a reference. Gray points represent the whole parent sample selected as in Sect. 3.2, black points show the color-selected-only quiescent candidate. Red points are “bona fide” objects, while orange points belong to the “extended” sample, both defined in Sect. 4. The mean photometric errors of the color-selected candidates (black points) are also shown. Objects that are non-detected in one or more filters in the color-selected sample are not shown in order to improve clarity given their large number. Only undetected objects belonging to extended or bona fide samples are shown with arrows indicating the lower limits of the involved color. Non-detections correspond to < 2σ.

In the text

	Fig. 4. Distribution of physical properties of the color selected sample (gray) of the extended (orange) and bona fide (red) sub-samples.
In the text

Fig. 5. Best-fit models for all bona fide objects: colored points represent photometric points with error-bars, red lines are best-fit SEDs. The best fit using the star-forming template is shown with a blue line. Error bars in the wavelength scale refer to the width of the considered filter. Physical parameters evaluated through the SED fitting are shown in inserted labels. Black tilted crosses and stars show the synthetic photometry, evaluated by integrating the best-fit templates in the several bands, for the best-fit quiescent and star-forming templates respectively. All filters used in the fit are listed in inserted legends for each plot.

In the text

Fig. 6. Number and stellar mass densities of galaxies with M ≥ 1010.6 M⊙ as function of redshift for total population of galaxies (in gray) or quiescent galaxies (colored points) in literature, compared to those obtained in this work. Our results are shown as red diamonds and squares for the BC03 extended and bona fide samples, respectively, with error bars representing the Poissonian errors, cosmic variance, and photometric redshift errors added in quadrature. The forecasts of two semi-analytic models are also shown: the black dashed and continuous lines represent the Guo et al. (2013) and the H15 models, respectively, estimated for the entire population of galaxies and convolved for the Eddington bias. The black dotted line represents the H15 model for the entire population not convolved for the Eddington bias. Blue and cyan shaded regions represent the H15 model with Poissonian error for the selected quiescent populations characterized by log(sSFR [yr−1])≤ − 10.5 and ≤ − 11, respectively.

In the text

Fig. 7. Stellar mass function of extended galaxies, in two redshift bins between z = 2.5 and 4.5. In each panel, the data points are shown as red circles in bins of Δlog M/M⊙ = 0.1. Error bars include Poisson noise, cosmic variance, and photometric redshift errors. The data points are fitted by a single Schechter function, shown by a red solid line, while the red shaded area is its 1σ uncertainty. Data points (with their error bars) from other works are also shown: violet squares are taken from Davidzon et al. (2017), in magenta triangles we plot the mass functions by Ilbert et al. (2013), and in orange circles the one by Muzzin et al. (2013). The blue lines shows the H15 model for quiescent galaxies at z = 2.68 (left) and z = 3.44 (right).

In the text

Fig. A.1. Effect of emission lines on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies (magenta and cyan lines), representing a population 1 Gyr old redshifted to z = 3 and a population 0.5 Gyr old redshifted to z = 4 (i.e., zform ≈ 6) built with BC03 models. A star forming SED with emission lines built from fsps, 1 Gyr old and redshifted to z = 2.5 is shown in dark-cyan. Also the filters transmission curves, used in COSMOS field and in the tracks computation, are shown. Bottom: evolution in color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue, black dark-cyan, and cyan, and their characteristics are shown in the legend inside the plot. A track for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) is shown in magenta.

In the text

Fig. A.2. Effect of different extinction laws on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies (magenta and cyan lines), that is, a population 1 Gyr old redshifted to z = 3 and a population 0.5 Gyr old redshifted to z = 4 (i.e., zform ≈ 6) built with BC03 models. Bottom: evolution in the color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in solid blue, solid black and dashed cyan lines, for the Calzetti, Seaton and Fitzpatrick laws, respectively, and their characteristics are shown in the inserted legends. Tracks for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) are shown in solid magenta, solid orange, dashed red and brown for the Calzetti (AV = 0 and AV = 1), Seaton and Fitzpatrick laws, respectively.

In the text

Fig. A.3. Effect of the choice of SPS models on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies, that is, a population 1 Gyr old redshifted to z = 3 (i.e. zform ≈ 6) built with BC03 models (magenta line) and with M05 models (orange line). Bottom: evolution in color–color plots for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan, while their characteristics are shown in the inserted legends. Tracks for quiescent galaxies (with e-folding time of τ = 0.1 Gyr and solar metallicity) are shown in magenta and orange for BC03 and M05, respectively.

In the text

Fig. A.4. Effect of metallicity on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown. In particular, the magenta SED represents a population 1 Gyr old redshifted to z = 3 (i.e. zform ≈ 6) with Z = Z⊙, while the green line shows the SED of a population of a 1 Gyr redshifted to z = 3 (i.e. again zform ≈ 6) with Z = 0.2 Z⊙. Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan, while their characteristics are shown in the inserted legends. Three tracks for quiescent galaxies are shown in magenta, dark-green and light-green for different metallicities (Z = Z⊙, Z = 0.4 Z⊙, Z = 0.2 Z⊙ respectively).

In the text

Fig. A.5. Effect of SFHs on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown in magenta SED and cyan, representing a population of 1 Gyr redshifted to z = 3 (i.e. zform ≈ 6), and a population of 0.5 Gyr redshifted to z = 4 (i.e., again zform ≈ 6), respectively. Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan and their characteristics are shown in the inserted legends. Three tracks for quiescent galaxies are shown in magenta, orange and red: the magenta track is parametrized with an exponentially declining star formation with τ = 0.1 Gyr, orange track has τ = 0.3 Gyr, and red track has been parametrized with a delayed star formation history where τ = 0.1 Gyr.

In the text

Fig. A.6. Effect of formation redshift on JKs[3.6][4.5] and HKs[3.6] color–color diagrams. Top panels: two SEDs of quiescent galaxies are shown. In particular, magenta SED represents a population of 1 Gyr redshifted to z = 3 (i.e. zform ≈ 6), while cyan line shows SED for a population of 1.5 Gyr redshifted to z = 3 (i.e. zform ≈ 8). Bottom: color–color plots with evolutionary tracks for star-forming and passive galaxies. Tracks representing star-forming galaxies are shown in blue and cyan and their characteristics are shown in the inserted legends. Two tracks for quiescent galaxies with redshift of formation zform = 6 and zform = 8 are shown in magenta and orange respectively. Their e-folding time is τ = 0.1 Gyr and they have solar metallicity.

In the text