LARgE Survey. I. Dead Monsters: the Massive End of the Passive Galaxy Stellar Mass Function at Cosmic Noon

TL;DR

This study presents the largest survey of massive quiescent galaxies at redshift ~1.6, revealing the stellar mass function's shape and implications for galaxy evolution mechanisms during cosmic noon.

Contribution

It introduces a large, statistically robust catalog of z~1.6 quiescent galaxies and analyzes their stellar mass function, highlighting the early establishment of mass quenching mechanisms.

Findings

The stellar mass function follows a Schechter function shape.

Mass quenching mechanisms were already established by z~1.6.

Hints of evolution in quenching efficiency compared to lower redshifts.

Abstract

We introduce the largest to date survey of massive quiescent galaxies at redshift z~1.6. With these data, which cover 27.6 deg^2, we can find significant numbers of very rare objects such as ultra-massive quiescent galaxies that populate the extreme massive end of the galaxy mass function, or dense environments that are likely to become present-day massive galaxy clusters. In this paper, the first in a series, we apply our gzK adaptation of the BzK technique to select our z~1.6 galaxy catalog and then study the quiescent galaxy stellar mass function with good statistics over M_stars ~ 10^10.2 - 10^11.7 M_sun --- a factor of 30 in mass --- including 60 ultra-massive z~1.6 quiescent galaxies with M_stars > 10^11.5 M_sun. We find that the stellar mass function of quiescent galaxies at z~1.6 is well represented by the Schechter function over this large mass range. This suggests that the mass quenching mechanism observed at lower redshifts must have already been well established by this epoch, and that it is likely due to a single physical mechanism over a wide range of mass. This close adherence to the Schechter shape also suggests that neither merging nor gravitational lensing significantly affect the observed quenched population. Finally, comparing measurements of M* parameters for quiescent and star-forming populations (ours and from the literature), we find hints of an offset (M*_SF > M*_PE), that could suggest that the efficiency of the quenching process evolves with time.

Figures

Tables

Table 1. Table 3: Best fit Schechter function parameters for quiescent galaxies at 1.5 ≤ 1.5 absent 1.5\leq z ≤ 2.0 absent 2.0 \leq 2.0 . When two values are listed under area then one refers to the wide and the other to the deep component of a "wedding cake" analysis. We present two results for our sample: one with a free α 𝛼 \alpha and a fixed α = − 0.14 𝛼 0.14 \alpha=-0.14 (the value given in Tomczak et al. ( 2014 ) ). Similarly, Muzzin et al. ( 2013b ) have a free- α 𝛼 \alpha and a fixed α = − 0.4 𝛼 0.4 \alpha=-0.4 result.

Survey	Field(s)	area[deg2]	selection	log($M^{\ast }\left)\right[M_{\odot }]$	${\Phi }^{\ast }\left[{10}^{-4}Mp{c}^{-3}\right]$	$\alpha$
LARgE (this work), free $\alpha$	CFHTLS W+D	25.1+2.5	PE-$gz{K}_s$	10.59 $\pm$ 0.02	2.35 $\pm$ 0.05	0.54 $\pm$ 0.08
LARgE (this work), fixed $\alpha =-0.14$	CFHTLS W+D	25.1+2.5	PE-$gz{K}_s$	10.74 $\pm$ 0.01	2.29 $\pm$ 0.03	$-$0.14
Muzzin et al. (2013b), free $\alpha$	COSMOS	1.6	photo-z + SED	$10.67\pm 0.03$	${4.15}_{-0.08}^{+0.06}$	0.03 $\pm$ 0.11
Muzzin et al. (2013b), fixed $\alpha =-0.4$	COSMOS	1.6	photo-z + SED	$10.80\pm 0.01$	${3.61}_{-0.04}^{+0.02}$	-0.4
Ilbert et al. (2013)	COSMOS	1.6	photo-z + SED	10.73 ${}_{}{}^{+0.03}{}_{-0.04}{}^{}$	2.20 ${}_{}{}^{+0.01}{}_{-0.01}{}^{}$	0.10${}_{}{}^{+0.09}{}_{-0.09}{}^{}$
Tomczak et al. (2014)	NEWFIRM+ZFOURGE	0.4+0.1	photo-z + SED	10.76 $\pm$ 0.05	3.29 $\pm$ 0.05	-0.14 $\pm$ 0.12

Equations

$$(z-K_s)-1.27(g-z)\ge -0.022,$$

$$(z-K_s)-1.27(g-z)<-0.022\cap (z-K_s)\ge 2.55,$$

$$(z-H)>2.4(H-K_s)+1,$$

$$V_{eff}(m)=A\int \frac{dV}{dz}p(m,z)dz,$$

$$log[M_{\star }/M_{\odot }]=-0.348K_s+\mathrm{18.284.}$$

$$\Phi (M)=ln\left(10\right){\Phi }^{\ast }\left[10^{\left(M-M^{\ast }\right)\left(1+\alpha \right)}\right]exp\left[-10^{(M-M^{\ast })}\right],$$

$${\Phi }_{convolved}(M)=\underset{-\infty }{\overset{\infty }{\int }}\Phi (M^{\prime })G(M-M^{\prime },\sigma )d{M}^{\prime },$$

$$G(M-M^{\prime },\sigma )=\frac{1}{\sigma \sqrt{2\pi }}\exp \left(\frac{M-M^{\prime }}{2{\sigma }^2}\right)$$

Full text

††pubyear: 2018††pagerange: LARgE Survey. I. Dead Monsters: the Massive End of the Passive Galaxy Stellar Mass Function at Cosmic Noon–LARgE Survey. I. Dead Monsters: the Massive End of the Passive Galaxy Stellar Mass Function at Cosmic Noon

1 Introduction

The rate at which the Universe was forming stars peaked during “cosmic noon” at z$$\sim1–3 (see Madau & Dickinson, 2014, for a review), but some galaxies were already quiescent by that epoch and in some cases perhaps even earlier. In contrast to star-forming galaxies whose numbers increase rapidly with decreasing stellar mass (Sawicki, 2012; Ilbert et al., 2013; Muzzin et al., 2013b; Davidzon et al., 2017), the most common quiescent galaxies at cosmic noon were those with stellar masses111We use the Chabrier (2003) stellar Initial Mass Function (IMF) throughout. of M_{\star}$$\sim 10^{10.8}$${\rm M}_{\odot} (Arcila-Osejo & Sawicki, 2013; Ilbert et al., 2013; Muzzin et al., 2013b; Davidzon et al., 2017). While a second, distinct population of lower-mass (M_{\star}$$\la 10^{9.5}$${\rm M}_{\odot}) quiescent galaxies already appears to exist by z$$\sim2 (Tomczak et al., 2014), and could potentially extend to quiescent galaxies with very low masses (Peng et al., 2010), the observed low-mass quiescent galaxies at these redshifts are few in number compared to their more numerous massive (M_{\star}$$\sim 10^{10.8}$${\rm M}_{\odot}) cousins that dominate the quiescent population both by number and by integrated stellar mass.

A significant number of such massive (M_{\star}$$\sim 10^{10.8}$${\rm M}_{\odot}) quiescent galaxies is now spectroscopically confirmed at $1\la z\la 3$ (Dunlop et al., 1996; Cimatti et al., 2004; Glazebrook et al., 2004; Cimatti et al., 2008; McCarthy et al., 2008; Gobat et al., 2012; Onodera et al., 2012; Belli et al., 2014a, b, 2016) and their existence provides an important test of galaxy evolution models (e.g., Somerville et al., 2012; Merson et al., 2013; Wellons et al., 2015; Behroozi & Silk, 2018). However, a handful of spectroscopically-confirmed examples of even more massive (M_{\star}$$\sim 10^{11.5}$${\rm M}_{\odot}) quiescent galaxies is now known (e.g., Onodera et al., 2012; Belli et al., 2016; Kado-Fong et al., 2017) and given their extreme masses and very early inferred formation and quenching times ($z\ga 4$, Belli et al. 2018; see also Pacifici et al. 2016; Glazebrook et al. 2017) they could provide even more stringent tests of our galaxy formation models. Moreover, ultra-massive passiveley evolving galaxies (UMPEGs) at cosmic noon, which for concreteness we henceforth define to be at $1<z<3$ and to have M_{\star}$$>10^{11.5}$${\rm M}_{\odot}, could be related to other interesting galaxy populations, including massive dusty star-forming galaxies at higher redshifts and cluster brightest central galaxies (BCGs) in clusters at lower redshifts.

UMPEGs could have of course assembled from previously quenched lower-mass progenitors such as are found in some z$$\sim2 clusters (Gobat et al., 2013; Newman et al., 2014), or formed via the merger-induced quenching of merging star-forming galaxies found in, say, high-$z$ protoclusters (e.g., Steidel et al., 2000; Capak et al., 2011; Cucciati et al., 2014; Oteo et al., 2018). However, if UMPEGs formed directly through the quenching of star formation in equally ultra-massive star-forming galaxies, then these star-forming M_{\star}$$>10^{11.5}$${\rm M}_{\odot} progenitors would have to have had star formation rates (SFRs) of 200-1000 ${\rm M}_{\odot}$/yr if, just prior to their quenching, they were on (Whitaker et al., 2012; Whitaker et al., 2014) or above (Elbaz et al., 2018) the $z\ga 2$ star-forming main sequence (SFMS, Daddi et al., 2007; Elbaz et al., 2007; Noeske et al., 2007; Sawicki et al., 2007). Such high SFRs and stellar masses are typical of Submillimetre Galaxies (SMGs) at $z\ga 2$ (e.g., Michałowski et al., 2012; Michałowski et al., 2014), or super-Main Sequence massive starbursts (Elbaz et al., 2018), suggesting that massive dusty star-forming galaxies could be the direct progenitors of high-$z$ UMPEGs.

The stellar masses of UMPEGs are also consistent with the stellar masses of the brightest cluster galaxies (BCGs) of massive clusters at low redshift (M_{\star}$$\sim 10^{11.7} ${\rm M}_{\odot}$, with a few examples at M_{\star}$$>10^{12} — see compilation by Lidman et al., 2012). Central cluster galaxies of similarly high masses are also known to z$$\sim1 and, for a few cases, to z$$\sim1.5 (Stott et al., 2010; Lidman et al., 2012). Additionally, BCG formation models suggest that cluster massive central galaxies may have already been quiescent and very massive (M_{\star}$$\sim a few $\times 10^{11}$) by z$$\sim1.5 (De Lucia & Blaizot, 2007). Altogether, these comparisons suggest that UMPEGs could be the progenitors of present-day Brightest Cluster Galaxies (BCGs).

To further test these scenarios we need to compare the number densities and clustering strengths of UMPEGs with those of the other galaxy populations of interest. Additionally, the detailed shape of the stellar mass function (SMF) of quiescent galaxies can carry in it information about the quenching process that shuts down star formation (e.g., Peng et al., 2010, 2012) making the shape of the quiescent galaxy SMF a test of quenching models. However, at present even the number density of the most common quiescent galaxies at z$$\sim2, those with M_{\star}$$\sim 10^{10.8}$${\rm M}_{\odot}, is not yet satisfactorily constrained by observations given both the cosmic variance affecting the relatively small, degree-scale fields used in $z>1$ mass function studies as well as differences in how samples are selected (see § 4). The number density of the truly massive UMPEGs is even more poorly constrained because the quiescent galaxy stellar mass function appears to drop rapidly at very high masses. Our goal then is to measure the SMF of high-$z$ quiescent galaxies over a large mass range, including very high UMPEG masses.

In constraining the number density of UMPEGs we face two problems: (1) UMPEGs are rare, while (2) the selection of high-$z$ quiescent galaxies requires deep near-infra-red imaging and very deep optical imaging. Consequently, to date, at high redshift ($z\ga 1.5$) studies have been limited to a handful of degree-scale — or even smaller — fields (Hartley et al., 2008; McCracken et al., 2010; Furusawa et al., 2011; Arcila-Osejo & Sawicki, 2013; Ilbert et al., 2013; Muzzin et al., 2013b; Ishikawa et al., 2015), some of which, such as the COSMOS field, are duplicated between studies. These studies do suggest, however, that the number density of M_{\star}$$>10^{11.5}$${\rm M}_{\odot} UMPEGs at z$$\sim1.5–2.5 is $\la 10^{-5}$Mpc*-3*, albeit with large statistical uncertainties, significant but poorly constrained cosmic variance given the very biased regions that the UMPEGs can be expected to reside in, and systematic differences due to different sample selection methods. These low UMPEG number densities, one twentieth or less those of the most common $z\sim 2$ passive galaxies around M_{\star}$$\sim 10^{10.8}$${\rm M}_{\odot}, indicate that to study this rare population we need to survey areas much larger than the degree-scale surveys carried out to date.

To better understand both the evolutionary pathway of UMPEGs and their relationship to other populations at high and low redshift, we need to first identify a statistically significant sample that is embedded in a full range of environments and is free of cosmic variance. This calls for a dataset that covers an order of magnitude more area than previous studies. With this sample, it will be possible to constrain the poorly-known number density of high-$z$ UMPEGs, to relate their stellar masses to those of their dark matter halos through clustering analysis, and to examine the environments in which they reside. These are the goals of our program, which we call the Large Area Red $gzK$ Environments (LARgE) Survey, and in which we set out to find and then study the properties of a large sample of z$$\sim1.6 ultra-massive passive galaxies. In this paper, the first in the LARgE series, we describe the dataset we use and the resultant sample of colour-colour-selected high-$z$ galaxies, and present their number counts and passive galaxy stellar mass function. In subsequent papers in the series we will measure the masses of the dark matter halos of the most extreme quiescent galaxies via galaxy-galaxy clustering (Gurpreet Kaur Cheema et al., submitted to MNRAS), constrain the growth of these galaxies via galaxy-galaxy mergers (Liz Arcila-Osejo, in prep.), and identify a sample of massive galaxy cluster candidates at z$$\sim1.6 identified by the presence of massive quiescent galaxies (Liz Arcila-Osejo et al., in prep.).

Throughout this work, and unless otherwise noted, we use the AB magnitude system (Oke, 1974), the ($\Omega_{M}$, $\Omega_{\Lambda}$, $h$) = (0.3, 0.7, 0.7) cosmology, and the Chabrier (2003) stellar initial mass function (IMF).

2 Data and Catalogs

Our project makes use of archival optical and NIR imaging, primarily from the Canada France Hawaii Telescope (CFHT): The Deep and Wide fields of the CFHT Legacy Survey (CFHTLS) for the optical data, and CFHT WIRCam IR imaging from the WIRDS and VIPERS-MLS projects for the IR data; a small amount of IR imaging in the WIRDS dataset comes from UKIRT. We perform object detection in the $K_{s}$-band images over these fields, do matched-aperture photometry in other bands, and finally select high-$z$ galaxies using an adaptation of the $BzK$ technique to our filter set.

In the following sections we first describe these datasets and then detail how we used them to construct our $K$-band selected catalogs and select high-$z$ galaxies.

2.1 Data

2.1.1 Optical Data: The CFHT Legacy Survey

The Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) is a large public project composed of two main surveys: Deep and Wide. The CFHTLS Deep component consists of four independent, very deep, one square degree pointings labelled D1, D2, D3 and D4. The CFHTLS Deep surve

y was carried out to detect type Ia supernovae and study the galaxy population to faint magnitudes. CFHTLS-Wide is divided into four fields (W1–W4) and its purpose is to study matter distribution, large scale structure and clusters of galaxies. For our work, we used all four Deep fields (D1–D4) and two of the Wide fields (W1 and W4) of the CFHTLS.

CFHTLS was carried out using MegaCam, the wide-field optical imager at MegaPrime. Each MegaCam CCD pointing covers a $\sim$1 $\deg$ $\times$ $\sim$1 $\deg$ field of view with a sampling of 0.186 arc seconds per pixel. All CFHTLS fields were observed in five broad-band filters ($u{*},g^{\prime},r^{\prime},i^{\prime}$, (or its replacement, $i2^{\prime}$), and $z^{\prime}$). From brevity we will drop the superscripts and use $ugriz$ to designate this MegaCam filter set. In our work we primarily make use of the $g$ and $z$ filters, as described in Section 3.1. For the Deep fields, the 80$\%$ completeness limit for extended objects in g-band is 25.29, 25.30, 25.29 and 25.26 and for the z-band it is 23.83, 23.90, 23.71, 23.76 AB mags respectively for D1, D2, D3 and D4 (Goranova et al. 2010). Likewise, for the Wide fields, W1 and W4, in the g band the completeness limit is 24.67 and 24.71, and in the $z$ band is 22.91 and 22.90, respectively (Hudelot et al. 2012).

CFHTLS data processing was done at Terapix222http://terapix.iap.fr/. Pre-processing involved image quality checking, flat fielding, stacking from dithering, identification of bad pixels, removal of cosmic rays and saturated pixels, background estimation, and astrometric and photometric calibration. Finally, Terapix is also involved in archiving and providing data products to the scientific community. For the Deep fields we used the Terapix T0006 processed and calibrated stacks, while for the Wide fields, we use the T0007 release. One of the key differences between these two releases is a different photometric calibration and, be comparing data in the T0006 and T0007 versions of the Deep fields we noticed that there was a systematic calibration offset in the z-band, such that $z_{T0007}=z_{T0006}+0.1$ (see also Moutard et al., 2016a, who reported a shift of 0.148 $\pm$0.054). We thus applied a 0.1 mag offset to the $z$-band zeropoints in the Wide data (for which we used the T0007 release) in order to make these data consistent with the Deep fields (for which we used the T0006 dataset). We note here that we found no systematic shift in the other optical bands.

2.1.2 Infrared Data: WIRDS and VIPERS-MLS

To complement the Deep survey optical data we used the WIRCam Deep Survey (WIRDS; Bielby et al., 2012). WIRDS was a large project carried out during 2006 to 2008 to obtain near infrared, broadband imaging for the four CFHTLS Deep fields. The results for this survey included images from the J, H and Ks filters that have been re-sampled to match the pixel scale of MegaCam. Most of the WIRDS were taken with the Wide-field Infrared Camera imager (WIRCam) at the CFHT, although the $J$-band imaging of the D2 field was obtained using the Wide Field Camera (WFCAM) instrument on the United Kingdom Infrared Telescope (UKIRT). Most of the effective areas for the Deep fields are less than one square degree (with the exception of D2). We used the T0002 WIRDS data release which we obtained from Terapix. The WIRDS 50% completeness limits for point-like objects is reached at $K$=24.5 AB, except for the D2 (COSMOS) field which reaches $K$=24.0 AB (Bielby et al. 2012) .

For the Wide fields we used the $K_{s}$-band images taken and processed as part of the Visible Multi-Object Spectrograph (VIMOS) Public Extragalactic Multi-Lambda Survey or VIPERS-MLS333http://cesam.lam.fr/vipers-mls/ (Moutard et al. 2016a). These data were also obtained with WIRCam on CFHTLS and were processed into square-degree patches that match in size and pixel scale the optical pointings of the CFHTLS-Wide. These data were designed to reach a $80\%$ completeness limit of Ks $\leq$ 22.0 AB. The Ks-band calibration of VIPERS-MLS matches that of the WIRDS Survey and also matches the VIDEO and UKIDSS calibrations where overlap exists (Moutard et al., 2016a).

2.2 K-selected Catalogs

2.2.1 Effective areas

In Figure 1 we show the layout of our g, z and $K_{s}$ coverage in W1 and W2 and $g$, $z$, $H$ and $K_{s}$ coverage in D1 to D4. We masked areas where object detection and/or photometry could be affected by the presence of bright star halos or other artefacts. For the Wide optical images, external flags of bright stars and cosmic ray trails were taken from Erben et al. 2012; we performed visual inspection of each one of these masks to ensure the masking was correct. For the corresponding Wide-field Ks band images we created our own external masks by visual inspection. Similarly, for the Deep Fields, both optical and infrared bands, we constructed our own external flags by visual inspection of the different images. After masking, our total remaining area is 27.60 deg2, including 25.09 deg2 in the Wide fields and 2.51 deg2 in the Deep fields. The areas of the individual fields, after masking of bright stars and other unusable subareas, are given in Table 1.

2.2.2 Photometry

Because quiescent galaxies at high redshift are expected to be very red ($g-K_{s}\gtrsim 5$ AB), it is critical to detect sources in the reddest available band ($K_{s}$ in our case) as optical detection could miss a large number of such objects. With this in mind we used the SExtractor package (Bertin & Arnouts 1996, version 2.19.5) to detect objects in the $K_{s}$-band images and then perform forced-aperture photometry in the other bands.

A source was identified in the Ks image based on it having at least five contiguous pixels above a detection threshold. The Deep-fields images are much deeper than the Wide and suffer from more correlated noise, which means that threshold may not be the same in Deep and Wide. By trial and error and visual inspection of SExtractor detections, we determined and subsequently used a threshold requirement of five contiguous pixels each above 1.5$\sigma_{sky}$ in the Wide areas, and 1.2$\sigma_{sky}$ in the Deep. Photometry is then performed at the $K_{s}$-band positions in $g$, $z$, and $H$ (in the Deep fields), using SExtractor’s dual-image mode, as explained below.

In all the fields, total $K_{s}$-band fluxes were calculated using auto mags (total magnitudes from a Kron-like apertures). The Deep and Wide fields were then treated differently when computing colours. In the Deep fields, which we used to push into the faint limit of the galaxy population, we wanted to maximise depth. For this reason we determined colours in fairly small (10 pix diameter, or 1.86”) apertures on $g$, $z$, and $H$ images whose PSF had been homogenised to match the $K_{s}$-band PSF (Sato et al. 2014). This approach follows exactly, and on the same data, that given in Arcila-Osejo & Sawicki (2013); it gives high S/N in the colours while ensuring the fluxes are measured over the same physical areas of the objects despite the small apertures. In the Wide fields, which we use to study the bright end of the population, we did not make PSF-homogenised images, but, instead, we used SExtractor’s total magnitudes (measured over large apertures determined individually for each object in the $K_{s}$-band) to measure colours; because for our relatively bright objects these total magnitude apertures are significantly larger than the PSF, this approach gives accurate colours without unacceptably degrading the S/N.

Finally, the photometry of each object was corrected for foreground Galactic dust. This was done using the Schlegel et al. (1998) Galactic dust maps, which are based on COBE and IRAS 100 $\mu m$ to 240 $\mu m$ emission maps, but using the updated reddening curves of Schlafly & Finkbeiner (2011) that are based on stellar spectra in the Sloan Digital Sky Survey (SDSS).

The result is a catalogue that includes $\sim$700,000 $K_{s}$-band selected objects in an effective area (i.e., after masking) of 27.60 deg2. See Table 1 for the breakdown by field.

3 High-z Object Selection and Number Counts

3.1 Selection of gzK galaxies

Daddi et al. (2004) developed a simple technique that uses the $B-z$ vs. $z-K$ colour-colour plane for selecting high-redshift galaxies and classifying them as either Passively Evolving ($pBzK$) or star-forming ($sBzK$). The technique is grounded in noting how galaxy spectral models (Daddi et al., 2004 used Bruzual & Charlot, 2003 models attenuated with Calzetti et al., 2000 dust) place in the $BzK$ colour plane. Spectroscopy, along with morphological studies, support the validity of $BzK$-selected passive and star-forming galaxies (Daddi et al. 2004, 2005; Ravindranath et al. 2007; Hayashi et al. 2009; Onodera et al. 2010; Mancini et al. 2010), and the simplicity and relative frugality of the technique (only three passbands are needed) continues to make it a popular choice for selecting star-forming and quiescent samples of distant galaxies (e.g., McCracken et al., 2010; Furusawa et al., 2011; Yuma et al., 2011; Kurczynski et al., 2012; McCracken et al., 2012; Yuma et al., 2012; Lee et al., 2013; Rangel et al., 2013; Sommariva et al., 2014; Fang et al., 2015; Ishikawa et al., 2015, 2016).

In Arcila-Osejo & Sawicki (2013) we adapted the classic $BzK$ technique to the $gzK_{s}$ bandbasses available in CFHTLS fields by noting how the Daddi et al. (2004) spectral synthesis models shift when redrawn in the g$$-$$z vs. z$$-$$K_{s} plane (the shifts are small given the similarity of the filter sets). For completeness, we give here the equations from Arcila-Osejo & Sawicki (2013) for selecting high-$z$ galaxies using this approach, noting again that — for brevity — hereafter we use $g$ and $z$ to denote the CFHT Megacam $g^{\prime}$ and $z^{\prime}$ filters. These $gzK_{s}$ selection criteria are

[TABLE]

for star-forming galaies, and

[TABLE]

for passive galaxies (Arcila-Osejo & Sawicki, 2013).

To select and classify high-redshift galaxies we then apply the cuts described by equations 1 and 2 to our Wide catalogues. This is illustrated in Figure 2 for a $\sim$ 1 deg2 subarea of our Wide catalogue (i.e., approximately 3% of our dataset). In common with the classic $BzK$ selection, the upper-left region of the diagram contains star-forming high-$z$ galaxies ($SF$-$gzK_{s}$, blue points) and the upper right has quiescent high-$z$ galaxies ($PE$-$gzK_{s}$, red points). Stars follow a distinct sequence in the lower part of the diagram (green points). Yellow points mark lower-$z$ galaxies.

In the Deep fields the $gzK_{s}$ technique can identify high-redshift galaxies as well, but by itself it is not sufficient to distinguish between star-forming and passive sub-classes. This is because to separate $PE$-$gzK_{s}$ from $SF$-$gzK_{s}$ galaxies we need to be able to measure colours to (g$$- $K_{s}$)$>$4.5 and in the Deep fields the $K_{s}$-band data are so deep that some $gzK_{s}$ galaxies are undetected in the $g$-band. While clearly at high redshift (because they have red $z$-$K_{s}$ colours), we cannot tell their evolutionary state from $gzK_{s}$ data alone. Instead, following Arcila-Osejo & Sawicki (2013), we use these galaxies’ $H-K_{s}$ colours to break the degeneracy. While in the absence of a $g$-band detection it is not clear whether the red $z$-$K_{s}$ colour is due to dust or the presence of the Balmer/4000Å break complex,the $H-K_{s}$ colour can remove the age-dust degeneracy, as shown in Arcila-Osejo & Sawicki (2013, see their Fig. 3). For completeness, we give here the relation from Arcila-Osejo & Sawicki (2013, their Eq. 6),

[TABLE]

which identifies quiescent galaxies from amongst $g$-undetected $gzK_{s}$ candidates.

Applying $gzK_{s}$ selection to the Wide fields and the g$$z$$K_{s}+$H$ method to the Deep, we then construct catalogues of star-forming ($SF$-$gzK_{s}$) and quiescent ($PE$-$gzK_{s}$) galaxies. We limit our catalogue depth in the $K_{s}$-band in each of the fields, as given in Table 2 in order to avoid incompleteness and – in the Wide fields – to ensure that the vast majority of objects are detected in the $g$-band.

To check that our $gzK_{s}$-selected objects are not point sources but are resolved — and therefore likely to be galaxies — we compared their peak $K_{s}$-band surface brightnesses (PSBs) with total magnitudes. Stars are compact (high PSB for their total magnitude) as compared to galaxies and present a clear and distinct locus in a PSB-mag diagram. All of our $PE$-$gzK_{s}$ objects with K${}_{s}<18.5$ (a total of five sources) are consistent in the PSB-vs-mag diagram with point source morphologies and we excluded them from further consideration. In comparison, only $\sim 7\%$ of $PE$-$gzK_{s}$ objects detected at $18.5\leq K_{s}<19.5$ and less than $0.5\%$ of $PE$-$gzK_{s}$ objects detected at $19.5\leq K_{s}<20.5$ are point sources in our data; they, too, have been removed from further analysis.

We also visually inspected the images of all $PE$-$gzK_{s}$ objects with $18.5<K_{s}<20$ in the Deep and Wide fields ($\sim$600 objects). Only six of these $\sim$600 objects (i.e., $\sim$ 1% of the total) were found to be corrupted (edge effects, closeness to bright-star diffraction spikes, corrupted pixels, etc.) and were removed from the sample.

A closer inspection of ultra-bright galaxies ($K_{s}<19.5$; 73 objects) revealed the presence of six possible late-stage major mergers: Galaxies with double cores that were not properly segmented by SExtractor. By visually identifying a boundary between the two sources, we recovered a separate background-corrected flux for each source in the pair, confirming the fact that what appears to be an ultra-bright passive galaxy could in fact be an ongoing, major merger between two fainter ($K_{s}\sim 20$) sources. We retained these six objects in our catalogue. These major merger candidates could potentially result, with time, in an ultra-bright passive galaxies.

Finally, we matched all K${}_{s}<20$ objects with the Chandra444http://cda.harvard.edu/chaser/, XMM555http://nxsa.esac.esa.int/nxsa-web/, NED666https://ned.ipac.caltech.edu/forms/nnd.html and Simbad777http://simbad.u-strasbg.fr/simbad/sim-fbasic databases in order to discard any possible AGN sources in our sample. In the Deep fields, no X-ray sources were found in the the XMM and Chandra databases but one object in D3 was identified in the Simbad catalogue as an AGN source in the Extended Groth Strip, belonging to the AEGIS-X Deep Survey (Goulding et al., 2012). This object was removed from our sample, as were two X-ray sources in field W1.

At the end of this process we were left with a large sample of $PE$-$gzK_{s}$ and $SF$-$gzK_{s}$ objects, all likely to be high-redshift galaxies. The numbers are listed in Table 2.

3.2 Number Counts

To study the number counts of $gzK_{s}$ galaxies we used the Wide fields for the bright end and the Deep fields for the faint end. We used simulations (Arcila-Osejo & Sawicki, 2013; Sato et al., 2014) to determine incompleteness corrections in the Deep fields, while in the Wide fields we truncated our analysis below magnitudes where the Wide field number counts diverge from the Deep fields data.

The incompleteness simulations in the Deep fields are described in detail in Sato et al. (2014). Briefly, artificial objects were added at random locations into the science images with empirically motivated morphological parameters. High-z star-forming galaxies were assumed to be disk-like objects with effective radii between 1 $\leq$ $r_{e}$ $\leq$ 3 kpc (Yuma et al., 2011) while passive galaxies present a more compact morphology and effective radii between 3 $\leq$ $r_{e}$ $\leq$ 6 kpc (Mancini et al., 2010). Once artificial galaxies are added into an image, SExtractor and catalogue-making are run with the same settings as were used for creating the object catalogue (Sec. 2) and the recovered/input numbers of simulated galaxies give detection completeness as a function of magnitude. The results are shown for $PE$-$gzK_{s}$ galaxies in Fig. 3 (a similar result was obtained for SF-$gzK_{s}$ galaxies but is not shown). Incompleteness corrections are simply 1/completeness and are applied to the observed number counts in the Deep fields. The resulting, incompleteness-corrected number counts are shown in Fig. 4.

Incompleteness corrections for the Deep fields are typically $\sim$1.02 over Ks = 17 − 22 (AB) mag but become larger at fainter levels. We limited our analysis in magnitude bins where incompleteness corrections are smaller than a factor of two, namely $K_{s}$ $<$ 23.5 mag (except for D2, where our passive population is only complete to $K_{s}<23.0$). For the Wide fields, instead of carrying out simulations we simply stop the analysis at magnitudes where the observed Deep-fields number counts diverge from the Wide-fields number counts. This happens at $K_{s}=20.5$, shown in Fig. 3 with the vertical dashed line. This approach allows us to constrain the bright end of the number counts without needing to go into computationally expensive simulations.

The number counts of $PE$-$gzK_{s}$ and $SF$-$gzK_{s}$ galaxies (corrected for incompleteness at the faint end in the Deep data) are shown in Fig. 4. We show each independent field (D1–D4, W1, W4) with a different symbol. Error bars are Gaussian $\sqrt{N}$ uncertainties. The small variations between our fields are most likely due to cosmic variance. Our total number counts are also presented in Table 2.

For both star-forming and passive galaxies, our results are in good agreement with the recent $BzK$ number counts of McCracken et al. (2010) and Sommariva et al. (2014), as shown in Fig. 5. However, we have better statistics and go significantly brighter in the population.

3.2.1 Comparison with models

Figure 5 shows a comparison between the observed and predicted number counts from semi-analytical models of passive $BzK$ galaxies. Shown in red points are our observations and for comparison, we plot three semi-analytical models: Kitzbichler & White (2007), Merson et al. (2013), and Somerville et al. (2012).

Kitzbichler & White (2007) compared observations of the high-redshift population with predictions from the Millennium simulation (Springel et al. 2005), which follows the hierarchical growth of dark matter structures from z=127 to the present, including several physical processes such as star-formation, gas cooling, growth of super-massive black holes, stellar population synthesis modelling for photometry and a new treatment for radio mode feedback only for galaxies at the centre of groups or clusters (hence not applied uniformly to all massive galaxies). Somerville et al. (2012) semi-analytical models, similar to those Kitzbichler & White (2007), include hierarchical dark matter growth, star-formation, gas cooling, supernova and AGN feedback, and metal enrichment. Feedback processes are mainly driven by supernovae at the faint end of the SMF and by AGN feedback at the massive end. Finally, Merson et al. (2013) populate dark matter halos from the Millennium simulation with galaxies using the GALFORM semi-analytical prescription (Bower et al. 2006). This prescription includes gas cooling, supernova and AGN feedback (only effective in systems with quasi-hydrostatic cooling). In their semi-analytical models, Merson et al. (2013) recalculated their predicted number counts in two different ways: one representing model $BzK$ galaxies with a simple K-mag limit brighter than Ks $\leq$ 24 and a second one which also takes into account a B-band detection limit of B $\leq$ 28.

The oldest model among those we examine (Kitzbichler & White, 2007) fails to reproduce the number counts at all magnitudes, while the newer models Somerville et al. (2012); Merson et al. (2013) do better, although the Merson et al. (2013) model underpredicts counts at the bright end. None of the models extend to the brightest passive galaxies that we observe ($K_{s}$<19.5) and examining such an extension would provide an interesting test of the models, and in particular whether their AGN feedback prescription can account for the observed number counts or needs to be modified or amended by additional physics.

4 Stellar Mass Function

4.1 Procedure

We take a simple approach to estimating the stellar mass function of our $PE$-$gzK_{s}$ galaxies: we start with the number counts already detailed in Sec. 3; we then estimate the (magnitude-dependent) redshift distribution of the population using 30-band photometric redshifts for our $PE$-$gzK_{s}$ galaxies in the D2 (COSMOS) field, which allows us to estimate effective volumes of our entire survey; and we convert $K_{s}$-band magnitudes to stellar masses using an empirical relation derived, again, in our D2 (COSMOS) subsample. These ingredients are described in detail below.

4.1.1 Redshift distribution and effective volumes

Although our Wide fields contain $ugriz+K_{s}$ data from which Moutard et al. (2016b) derived photometric redshifts, the lack of $J$ and $H$ bands in this dataset means that redshift estimates at redshifts of interest to us suffer from large degeneracies due to the inability to constrain the location of the Balmer/4000Å breaks. Therefore, instead of relying on photometric redshifts of individual objects, we take a statistical approach and use the photometric redshift distribution for a subsample of our objects where the photometric redshifts are likely to be very good, namely in the D2 (COSMOS) field. We cross-matched our D2 $PE$-$gzK_{s}$ galaxies with the Muzzin et al. (2013a) 30-band COSMOS photo-$z$ catalogue, resulting in 1384 $PE$-$gzK_{s}$ objects with excellent photo-z’s. The top panel of Fig 6 shows the result and it is clear that the redshift distribution is magnitude-dependent.

We smoothed our 2D histogram at each magnitude bin using a Gaussian kernel, which extrapolates at both edges of the distribution, i.e, the brightest and faintest bins, and we normalized the result by the area under the curve to obtain a redshift probability density distribution, shown in the middle panel of Fig. 6. The redshift distribution peaks at $z\sim 1.6$ for bright $PE$-$gzK_{s}$ galaxies and at somewhat higher redshifts for fainter ones; in all cases there is a non-negligible non-Gaussian tail to higher redshifts. For comparison, in the bottom panel of Fig. 6 we also show the redshift PDFs derived in a similar way, but using the photometric redshifts of Moutard et al. (2016b). We adopt the PDFs based on the Muzzin et al. (2013a) photometric redshifts for the rest of this paper.

We note that the initial assessment of the $BzK$ galaxy population redshift distribution by Daddi et al. (2004) was based on a much smaller sample with a limited number of bright galaxies. Since most of their sample lies at fainter magnitudes, they reported a redshift distribution of $1.4<z<2.5$. With the advent of larger area surveys with good photometric redshifts we can see that the redshift distribution peaks at the lower end of this redshift range, particularly for the brighter objects.

We assume that the normalised redshift density distributions shown in the middle panel of Fig. 6 are a good proxy for our redshift probability density function, or $p(m,z)$, which we will use later in estimating the effective volume of the survey. Starting from the $p(m,z)$ we can then estimate the magnitude-dependent effective volumes $V_{eff}(m)$ for our sample. Following Steidel et al. (1999) and Sawicki & Thompson (2006), this is done using

[TABLE]

where A is the effective area of our survey (27.60 deg2 for galaxies with $K_{s}$< 20.5 and 2.51 deg2 for fainter galaxies), $dV/dz$ is the comoving volume per square degree and $p(m,z)$ is the redshift- and magnitude-dependent PDF that we estimated from the COSMOS photometric redshift analysis. Note that this approach is somewhat different from that of Steidel et al. and Sawicki & Thompson, since those authors assumed that not all the LBGs in their surveys were detected, whereas we assume that all of our $PE$-$gzK_{s}$ galaxies are.

4.1.2 Estimating masses from $K_{s}$ magnitudes

We estimated galaxy masses for our passive galaxies from their $K_{s}$ total magnitudes. At redshifts of z $\sim$ 1.6, the peak of the redshift distribution of luminous $PE$-$gzK_{s}$ galaxies, the observed $K_{s}$ band corresponds to rest-frame 9200Å and at these wavelengths the light from these galaxies is generated by the long-lived low-mass stars that are expected to contain most of the stellar mass of the galaxy.

To find a suitable conversion between magnitude and mass we used the same $\sim$1400 cross-matched $PE$-$gzK_{s}$ galaxies in D2 (COSMOS) that we used to work out the redshift distribution. To estimate galaxy stellar masses Muzzin et al. (2013a) followed the now standard SED-fitting approach first introduced by Sawicki & Yee (1998) in which multi-wavelength broadband photometry is compared with a grid of predictions derived from redshifted and dust-attenuated stellar population synthesis models. Here, Muzzin et al. (2013a) used a set of models with exponentially declining star-formation histories (SFHs), assumed solar metallicity, the Chabrier (2003) IMF, and allowed visual attenuation (AV) to vary between 0 and 4. As expected for passive galaxies, the masses scale well with observed $K_{s}$-band total magnitudes as shown in Fig. 7. The best-fit orthogonal distance correlation is shown as a solid line and is described by

[TABLE]

We use this relation to convert the observed magnitudes of our galaxies to stellar masses.

4.2 Results

Using our empirical magnitude-to-mass relation (Eq. 5), we estimated stellar masses for all our $PE$-$gzK_{s}$ galaxies from their $K_{s}$ magnitudes. We then counted $PE$-$gzK_{s}$ galaxy numbers as a function of stellar mass in logarithmically-spaced mass bins. Next, we corrected these numbers for incompleteness using the same approach as we used for the number counts in § 3 but after first using Eq.  5 to convert our magnitude-incompleteness estimates of § 3 into incompleteness as a function of mass. Finally, we divided the resulting numbers by the effective volumes from § 4.1.1 after using Eq. 5 once again to convert from magnitude-dependent effective volumes into mass-dependent effective volumes. The result is the stellar mass function shown with red points and connecting dashed red line in Fig. 8.

4.2.1 Direct comparison with previous work

Using the photometric redshift distribution along with the $K_{s}$-$M_{\star}$ relation (i.e., Eq. 5) we built a z$$\sim1.6 passive galaxy stellar mass function with our combined Deep and Wide catalogues. Our direct SMF measurement is shown with red circles in Fig. 8, along with the PE galaxy SMFs of Muzzin et al. (2013b), Ilbert et al. (2013), and Tomczak et al. (2014) at similar redshifts.

As can be seen in Fig. 8, there are significant differences among these previous studies, all of which are based on relatively small areas: 0.5 deg2 for Tomczak et al. (2014) and 1.6 deg2 for Muzzin et al. (2013b) and Ilbert et al. (2013) (both covering the same area in the COSMOS field). While cosmic variance can be expected to play a role in such relatively small fields, the observed scatter between surveys is not likely to be due to cosmic variance alone since the Muzzin et al. (2013b) and Ilbert et al. (2013) results, which differ significantly, are based on independent analyses of the same dataset in the COSMOS field. Instead, sample selection or other systematics may play a significant role in producing the discrepancies.

As Fig. 8 shows, our results are consistent with previous studies at intermediate and low masses (M_{\star}$$\la1011${\rm M}_{\odot}$), lying within the number density range spanned by these previous works. However, at larger masses, M_{\star}$$\ga1011${\rm M}_{\odot}$, we measure a lower number density than the previous studies. We note that our statistics at the massive end, based on a 27.6 deg2 dataset, are significantly better than any previous work. Our points shown in Fig. 8 have not been corrected for Eddington bias (Eddington, 1913); we will correct for this in § 4.2.2888Note that only some authors correct for Eddington bias: while (Ilbert et al., 2013) do, it appears that (Muzzin et al., 2013b) and (Tomczak et al., 2014) do not. . As noted in § 4.2.2, correcting for Eddington bias brings our number density down somewhat at the massive end, exacerbating the discrepancy with the previous studies.

The discrepancy between our measurement and previous work may be due to sample variance or cosmic variance dominating the smaller fields (we have more than an order of magnitude more area than the earlier studies and, as noted earlier, the results of Muzzin et al. (2013b) and Ilbert et al. (2013) are not independent as they are based on a single field). Alternatively, systematic differences in sample selection, assumptions about objects’ redshift distributions, or mass measurement may also play a role. Full resolution of this issue will require spectroscopy of a significant number of the massive galaxies in our sample. In the meantime, given that at intermediate and low masses our measurements are consistent with previous work, we feel it plausible that at the massive end our measurement is likely more correct than the previous work based on two relatively small fields.

4.2.2 Schechter function fitting and Eddington bias correction

The Schechter function (Schechter, 1976) is a fitting formula that, with just three adjustable parameters, is found to well describe the galaxy luminosity and stellar mass functions at many epochs (e.g., Lilly et al., 1995; Sawicki et al., 1997; Sawicki & Thompson, 2006; Peng et al., 2010; Ilbert et al., 2013; Muzzin et al., 2013b; Moutard et al., 2016b). The functional Schechter form can be physically motivated, but we leave such aspects to § 5 and in the present section focus on its purely operational role as a fundamental description of a galaxy population.

To find the Schechter function that best describes our population of passively evolving $PE$-$gzK_{s}$ galaxies, we performed a $\chi^{2}$ minimisation fit to our mass function data. However, since galaxy density decreases exponentially towards brighter magnitudes, in the presence of uncertainties in mass estimates, galaxies are more likely to scatter to higher than to lower masses. This effect, called the Eddington Bias (Eddington, 1913) can cause an apparent steepening of the stellar mass function at the massive end. To assess and correct for the effect of this bias, we also performed a $\chi^{2}$ fit on a Schechter function convolved with typical stellar mass errors.

To begin, we bootstrap resample our original sample, then perturb the masses in the bootstrapped catalogue using a normal distribution with mass-dependent width that characterises the scatter seen in Fig. 7. By performing the bootstrap and mass perturbation on our sample we account for sample variance and mass uncertainties, respectively, when calculating the best-fit model.

Next, for each bootstrap, we bin our data and perform a $\chi^{2}$ minimisation using the Schechter function

[TABLE]

Schechter (1976); Marchesini et al. (2009) that has been convolved with a Gaussian using

[TABLE]

where

[TABLE]

represents the scatter in the $K_{s}$-$M_{\star}$ relation. Because photometric uncertainties in the $K_{s}$-band are very small for the bright objects that populate the steep part of the SMF, the dominant source of Eddington Bias is the scatter in the $M_{\star}$-$K_{s}$ relation and this is taken into account by the convolution described by Eqs.  7 and 8.

The free parameters in the fit are $\alpha$, $M^{*}$ and $\Phi^{*}$, as defined by the Schechter function (Eq. 6). We performed 10,000 bootstrap resamples and carried out parameter minimisation on each realisation of the dataset. We define our best-fit model to be the peak of our output best-fit parameter space. The best-fit Schechter function, after correcting for Eddington bias, is shown as the red line in Fig. 9; the orange line shows the Schechter fit without the Eddington bias correction. We also show the original, uncorrected data (orange points) and, in red, data points corrected for Eddington bias by applying the ratio of the two (corrected and uncorrected) Schechter functions. The probability distributions of the Schecher parameters, based on the 10,000 resampled fits, are shown in Fig. 10.

In Table 3 we summarise our best-fit Schechter parameters, along with the results of three other recent studies (Ilbert et al., 2013; Muzzin et al., 2013b; Tomczak et al., 2014). Comparing the different results, we see that the values of $M^{*}$ cluster closely together, although are often formally inconsistent with one another given the small uncertainties. The values of characteristic density, $\Phi^{*}$ vary significantly between the surveys, in part, though probably not wholly, reflecting differences between PE galaxy selection criteria and — for the smaller-area surveys — cosmic variance. Finally, there is a large scatter in the values of the slope at the low-mass end, $\alpha$, along with large associated uncertainties. This is not unexpected, given that in order to constrain this low-mass end slope it is necessary to go deep below $\log(M^{*}/M_{\odot})<10.5$.

Our own Schechter fit is strongly influenced by the high-S/N points around $M^{*}$ where galaxies are most numerous intrinsically and where we also have excellent statistics from the Wide fields. Figure 8 shows that among the surveys considered only the sample of Tomczak et al. (2014) reaches deep into the low-mass end (albeit with poor statistics) and thus likely gives the best constraint on $\alpha$. For this reason we also re-fit our data fixing $\alpha$ to the $-0.14$ of Tomczak et al. Fixing $\alpha$ in this way changes our $M^{*}$ and $\Phi^{*}$ somewhat and the resulting Schechter parameters are also shown in Table 3 and illustrated with the dotted black line in Fig. 9.

5 Discussion

5.1 The growth of the quenched galaxy population

It is interesting to investigate how the quiescent galaxy population grows in time. We first examine the growth of the comoving stellar mass density, shown in Fig. 11. Given that we probe the peak of the PE galaxy SMF with excellent statistics derived from our large area, we are in an excellent position to constrain the global stellar mass density contained in passive galaxies at z$$\sim1.6.

In Fig. 11 red points show the growth measured by several studies of the stellar mass density due to passive galaxies with $8<\log$($M_{\star}$/${\rm M}_{\odot}$) $<13$ (except for Tomczak et al. (2014) who integrated down to $\log$($M_{\star}$/${\rm M}_{\odot}$)$=$8). All points have been changed, when necessary, to the Chabrier (2003) IMF. Black points show the growth of the SMD due to galaxies of all types while the red points show the contribution of the quiescent galaxy population. Our $PE$-$gzK_{s}$ galaxies’ contribution at z$$\sim1.6, shown with the filled red circle is calculated by integrating our PE galaxy Schechter function over $8<\log$($M_{\star}$/${\rm M}_{\odot}$) $<13$ and is consistent with previous studies. Overall, as is well known, the stellar mass density in passive galaxies grows more rapidly than the total stellar mass density, reflecting the growing importance of the quenched population with cosmic time.

In this project we are particularly interested in the extremely rare, ultra-massive passive galaxies and so we also calculate the stellar mass density of $PE$-$gzK_{s}$ galaxies with $11.5<\log$($M_{\star}$/${\rm M}_{\odot}$) $<13$ and show it with the red star in Fig. 11. We also show (filled red squares) the UMPEG contribution to the SMD at lower redshifts obtained by integrating the quiescent galaxy SMFs of Moutard et al. (2016b). While individually very massive, UMPEGs are also very rare, and so, as a population, do not contribute significantly to the total stellar mass density: they account for only $\sim 0.66\%$ of the mass contained in quiescent galaxies at z$$\sim1.6 and a similar fraction at lower redshifts. These are rare monsters indeed.

We next turn to the evolution of the passive galaxy stellar mass function, which is shown in Fig. 12. In the top panel of this Figure we show our z$$\sim1.6 quiescent galaxy stellar mass function (corrected for Eddington bias, § 4.2.2) along with PE SMFs at z$$\sim1 (Moutard et al., 2016b) and z$$\sim0 (Peng et al., 2010). The bottom panel shows the fractional change of the number density from z$$\sim1.6 to z$$\sim1 (orange pentagons) and from z$$\sim1.6 to z$$\sim0 (yellow diamonds). These ratios were calculated by interpolating between the points (and uncertainties) from each survey to the same log masses. Also shown, as lines, are the the ratios between the best-fit Schechter functions over the same redshift intervals.

Overall, the quiescent galaxy population can be expected to grow in numbers as galaxies quench their star formation and migrate from the star-forming population to the quiescent one (see, e.g., Peng et al., 2010). The rates at which different parts of the quiescent galaxy SMF grow will then depend on the quenching mechanism(s) and other processes that post-process quenched galaxies, such as merging. With this in mind, it is interesting to examine how the SMF evolves over time. As can be seen in the bottom panel of Fig. 12, the evolution of the PE galaxy SMF cannot be described as a uniform number density increase with time. Compared to the growth at intermediate masses, there is excess growth at the low-mass end, particularly from z$$\sim1 to z$$\sim0; there may also be excess growth at the high-mass end, albeit within large uncertainties. The excess growth at low masses can be associated with the increasing importance of the environmental quenching mechanism (a mass-independent quenching process, presumably associated with satellite quenching; Peng et al. 2010) which is expected to make an increasingly important contribution at later cosmic epochs as large scale structure develops in the Universe. In this context, it is not unexpected to see little excess growth in the low-mass end from z$$\sim1.6 to z$$\sim1, but then stronger growth down to z$$\sim0.

At the massive end it is hard to be sure of the differential growth between the high-mass and intermediate-mass parts of the SMF given the large uncertainties in the SMFs at lower redshifts. However, if the growth at the massive end is in excess of that at intermediate masses, it seems to be already happening between $z\sim 1.6$ and $z\sim 1$. Such excess growth of the SMF could mean that during these later epochs quenching galaxies are significantly more massive than those quenched at higher redshifts, or that there is significant post-quenching mass growth through, for example, ("dry") mergers.

5.2 The quenching of star formation in massive galaxies

The Schechter function is not just an empirical fitting formula that has been found to reasonably describe the galaxy SMF: its functional shape can be physically motivated. For example, Peng et al. (2010) postulated that the Shechter shape, both for star-forming and passive galaxy SMFs, is exactly produced by a “mass quenching mechanism" that is of unknown physical origin but in which probability that a star-forming galaxy will quench and join the quiescent population increases exponentially with that galaxy’s stellar mass. A separate “environmental quenching mechanism”, independent of a galaxy’s mass, is responsible for quenching star formation in lower-mass galaxies. The combination of these two mechanisms is found to be remarkably good at reproducing the observed shapes of the SMFs of both star-forming and quiescent galaxies at z$$\sim0 (Peng et al., 2010).

The effect of the environmental quenching mechanism is expected to grow with time as large scale structure develops in the Universe, and consequently it is minimised at high redshift. This means that at high redshift we can observe the mass quenching mechanism in isolation, uncontaminated by the admixture of environmentally-quenched galaxies. In this context, the mass quenching mechanism predicts that the SMF of PE galaxies should be described by a single Schechter function of the form of Eq. 6; it should be related to the SMF of SF galaxies, itself also a single Schechter function, through identical characteristic masses (i.e., $M^{*}_{PE}=M^{*}_{SF}$) and low-mass slopes that differ by $\alpha_{PE}-\alpha_{SF}\sim 1$ (Peng et al., 2010).

In § 4 (Fig. 9) we presented the SMF of z$$\sim1.6 PE galaxies. Using the statistics made possible by our wide-area survey, our work shows that the PE galaxy SMF is fit very well by the Schechter functional form all the way up to very high stellar masses (log($M_{\star}$/${\rm M}_{\odot}$)$>$11.5. The excellent Schechter-like PE galaxy SMF at z$$\sim1.6, including at the high-mass end, leads to the immediate conclusion that the mass quenching mechanism must have already been well established by z$$\sim1.6, and earlier still if one assumes a non-zero timescale for galaxies to transition from the star-forming to the quiescent populations. Furthermore, given the excellent Schechter-like PE galaxy SMF, this mass quenching mechanism is likely due to a single physical process over a wide range of galaxy mass, including at ultra-high mass (log($M_{\star}$/${\rm M}_{\odot}$)$>$11.5).

Furthermore, again given the very good Schechter-like PE galaxy SMF, gravitational lensing, which is invoked to explain the non-Schechter luminosity functions sometimes observed at higher redshifts (e.g., Hildebrandt et al. 2009; Ono et al. 2018) is unlikely to be significant for massive z$$\sim1.6 galaxies.

As noted above, the mass quenching mechanism predicts a straightforward relation between the PE and SF Schechter functions: $M^{*}_{PE}=M^{*}_{SF}$ and $\alpha=\alpha_{PE}-\alpha_{SF}\sim 1$. Given the large dispersion in the measurements of $\alpha$ in the literature (see Table 3), we are not in a position to test this prediction for the low-mass slopes. However, the PE galaxy $M^{*}_{PE}$ values are fairly similar between the different surveys (even if often formally inconsistent given the stated uncertainties), so we can test whether $M^{*}_{PE}=M^{*}_{SF}$. To this end, in Fig. 13 we compare measurements of $M^{*}_{PE}$ and $M^{*}_{SF}$ at z$$\sim1.5–2. We plot $M^{*}$ values from Ilbert et al. (2013), Muzzin et al. (2013b), and Tomczak et al. (2014) as well as our results (complemented with the $M^{*}_{SF}$ for star-forming galaxies from Ilbert et al. 2013). The diagonal line shows $M^{*}_{PE}=M^{*}_{SF}$. It is clear that in most cases $M^{*}_{SF}>M^{*}_{PE}$ (the exception is the measurement of Ilbert et al. 2013). The offset between $M^{*}_{SF}$ and $M^{*}_{PE}$ increases further by $\sim$0.1 dex (indicated by the blue arrow) when we correct for the fact that stellar masses of star-forming (but not quiescent) galaxies at high redshift are underestimated in most studies due to the effect of outshining. The blue arrow in Fig. 13 shows the outshining correction to star-forming galaxy masses inferred using Eq. 6 of Sorba & Sawicki (2018) for an estimated $SFR\sim 10^{2.1}$ for a galaxy with $M_{\star}$ $=M^{*}_{SF}=10^{11}M_{\odot}$ from Whitaker et al. (2014).

The offset between $M^{*}_{SF}$ and $M^{*}_{PE}$ seems significant, apparently in disagreement with the simple prediction of the mass quenching model of Peng et al. (2010). The offset cannot be due to the effect of galaxy-galaxy merging as merging would drive $M^{*}_{PE}$ up, while keeping $M^{*}_{SF}$ constant (the latter because a merger involving SF galaxies would increase the mass of the merger product thus making it immediately more susceptible to quenching). Since this behaviour in $M^{*}_{SF}$ and $M^{*}_{PE}$ is opposite of what we observe, we conclude that major merging does not play a significant role at high redshift.

Instead, the presence of the offset between $M^{*}_{SF}$ and $M^{*}_{PE}$ may suggest that we are seeing the efficiency of mass quenching evolve with time. This would mean that (i) the physical mechanism responsible for mass quenching is redshift-dependant (i.e., $\mu=\mu(z)$ in Eq. 23 of Peng et al. 2010), and (ii) that the mass quenching mechanism is not instantaneous but, instead, has a quenching time $\tau_{MQ}$ that is consistent with a quenching mechanism with long time-scales, which might be due to smooth feedback mechanism but also to viral-shock heating (Moutard et al., 2016b). Together, these two phenomena will cause the $M^{*}$ of both SF and PE populations to evolve with time, but with $M^{*}_{PE}$ lagging behind $M^{*}_{SF}$ by an offset related to the quenching timescale $\tau_{MQ}$.

6 Conclusions

In this paper, the first in a series on z$$\sim1.6 massive quiescent galaxies, we introduced our large sample of $gzK_{s}$-selected passive galaxies. These data have allowed us to constrain the number counts and mass function of this population with unprecedented precision, spanning the mass range $10.25\la$ log($M_{\star}$/${\rm M}_{\odot}$)$\la 11.75$. This mass range reaches extremely rare objects with ultra-high masses that were not previously included in significant numbers in smaller surveys. We find that our PE galaxy SMF can be well fit with a Schechter function over essentially this entire mass range. The fact that the SMF has the Schechter form is consistent with the mass quenching scenario of (Peng et al., 2010); interpreting our result in the context of this scenario, we draw the following conclusions:

Given the SMF of PE galaxies is well represented by the Schechter function at z$$\sim1.6, the mass quenching mechanism must have already been well established by this epoch. It may well have been established earlier than that if the time to transition from the star-forming to the quenched population is significant. 2. 2.

Given the adherence of the PE galaxy SMF to the Schechter function shape over a wide mass range, the mass quenching mechanism is likely a single physical process over a wide range of galaxy mass, including at ultra-high masses [log($M_{\star}$/${\rm M}_{\odot}$)$>$11.5]. 3. 3.

Gravitational lensing, which is invoked to explain departures from the Schechter shape at higher redshifts (e.g., Hildebrandt et al. 2009; Ono et al. 2018) is unlikely to be significant for massive z$$\sim1.6 galaxies as it would distort our SMF away from its Schechter function shape. 4. 4.

Comparing $M^{*}_{PE}$ and $M^{*}_{SF}$ measurements (from this study and from other surveys) we find $M^{*}_{SF}>M^{*}_{PE}$, in apparent conflict with the simple mass-quenching model. This offset can, however, be explained if the quenching efficiency ($\mu$ in the notation of Peng et al. 2010) evolves with time. Additionally, this explanation for the offset between $M^{*}_{PE}$ and $M^{*}_{SF}$ requires a slow quenching time-scale as the $M^{*}_{PE}$ lags behind $M^{*}_{SF}$. 5. 5.

Galaxy-galaxy merging would distort the Schechter shape and would result in $M^{*}_{PE}>M^{*}_{SF}$, opposite to what is observed. Merging is therefore not likely to be a strong mechanism for the growth of masses of ultra-massive galaxies at high redshift.

In addition to the results enumerated above, in this paper we have developed and described a very large sample of z$$\sim1.6 quiescent galaxies. In forthcoming papers in this series we will study this massive quiescent population: we will study the clustering of ultra-massive [log($M_{\star}$/${\rm M}_{\odot}$) $>$ 11.5] members of the quiescent population (G. Cheema, MNRAS, submitted) as well as their environments and growth through mergers (L. Arcila-Osejo, in prep.). And we will identify and study protoclusters of quiescent galaxies found in our catalog (L. Arcila-Osejo, in prep.).

Acknowledgements

We thank Ivana Damjanov, Laura Parker, and Rob Thacker for useful suggestions and the Natural Sciences and Engineering Research Council (NSERC) of Canada for financial support.

This work is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work uses data products from TERAPIX and the Canadian Astronomy Data Centre. It makes use of the VIPERS-MLS database, operated at CeSAM/LAM, Marseille, France. This work is based in part on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada and France. The research was carried out using computing resources from ACEnet and Compute Canada.