THE COSMOS-WIRCam NEAR-INFRARED IMAGING SURVEY. I. BzK-SELECTED PASSIVE AND STAR-FORMING GALAXY CANDIDATES AT z ≳ 1.4^*

WIRCam (Puget et al. 2004) is a wide-field near-infrared detector at the 3.6 m CFHT, which consists of four 2048 × 2048 cryogenically cooled HgCdTe arrays. The pixel scale is 03 giving a field of view at prime focus of 21' × 21'. The data described in this paper were taken in a series of observing runs between 2005 and 2007. A list of these observations and the total amount of on-sky integration time for each run can be found in Table 1.

Observing targets were arranged in a set of pointing centers arranged in a grid across the COSMOS field. At each pointing center, observations were shifted using the predefined "DP10" WIRCAM dithering pattern, in which each observing cube of four micro-dithered observations is offset by 1'12'' in R.A. and 18'' in decl. Our overall observing grid was selected so as to fill all gaps between detectors and provide a uniform exposure time per pixel across the entire field and to ensure a the WIRCAM focal plane was adequately sampled to aid in the computation of the astrometric solution. All observations were carried out in queue-scheduled observing mode. Our program constraints demanded a seeing better than 08 and an air mass less than 1.2. Observations were only validated by CFHT if these conditions were met. In practice, a few validated images were outside these specifications (usually due to short-term changes in observing conditions), and were rejected at later subsequent processing steps.

Observations were made K_s filter ("K-short"; Skrutskie et al. 2006), which has a bluer cutoff than the standard K filter, and unlike the K' filter (Wainscoat & Cowie 1992) has a "cut-on" wavelength close to standard K. This reduces the thermal background and for typical galaxy spectra decreases the amount of observing time needed to reach a given signal-to-noise ratio (S/N) compared to the standard K-filter. A plot of the WIRCam K_s filter is available from CFHT.¹⁵

WIRCam images were pre-reduced using IRAF¹⁶ using a two-pass method. To reduce the data volume the four micro-dithers were collapsed into a single frame at the beginning of the reduction process. This marginally reduced the image quality (by less than 01 FWHM) but made the reduction process manageable on a single computer. After stacking the sub-frames the data were bias subtracted and flat fielded using the bias and flat-field frames provided by the CFHT WIRCam queue observing team. A global bad pixel mask was generated using the flat to identify the dead pixels; the dark frames were used to identify hot pixels. A median sky was then created and subtracted from the images using all images in a given dither pattern. The images were then stacked using integer pixel offsets and the headers astrometric information using IRAF's imcombine task. These initial stacks of the science data are used to generate relatively deep object masks through SExtractor's 'CHECKIMAGE_TYPE = OBJECTS' output file. These object masks are used to explicitly mask objects when re-generating the sky in the second pass reduction. Supplementary masks for individual images are made to mask out satellites or other bad regions not included in the global bad pixel mask on a frame-by-frame basis.

In the second pass reduction, we individually sky-subtract each science frame using the object masked images and any residual variations in the sky is removed by subtracting a constant to yield a zero mean sky level. These individual sky-subtracted images of a single science frame are then averaged with both sigma-clipping to remove cosmic rays and masking using a combination of the object mask and any supplementary mask to remove the real sources and bad regions in the sky frames. The region around the on-chip guide star is also masked. These images are further cleaned of any non-constant residual gradients as needed by fitting to the fully masked (object + supplementary + global bad pixel masks) background on a line-by-line basis. Any frames with poor sky subtraction or other artifacts after this step were rejected.

Next amplifier cross talk is removed for bright Two Micron All Sky Survey (2MASS) stars which creates "donut" shaped ghosts at regularly spaced intervals. These are unfortunately variable, with their shape and level depending on the brightness of the star and the amplifier the star falls on. We proceed in three steps: we first build for each star a median image of the potential cross talk pattern it could generate. This is done by taking the median of 13 × 13 pixels sub-frames at 64 pixel intervals above and below the star position, taking into account the full bad pixel mask. Next, we check whether the cross talk is significant by comparing the level of the median image to the expected noise. If cross talk is detected, we subtract it by fitting at every position it could occur the median shape determined in the preceding steps.

After pre-reduction, the TERAPIX tool QualityFITS was used to create weight-maps, catalogs, and quality assessment Web pages for each image. These quality assessment pages provide information on the instrument point-spread function (PSF), galaxy counts, star counts, background maps, and a host of other information. Using this information, images with focus or electronic problems were rejected. QualityFITS also produces a weight-map for each WIRCam image; this weight-map is computed from the bad pixel mask and the image flat field. Observations with seeing FWHM greater than 11 were rejected.

In the next processing step the TERAPIX software tool Scamp (Bertin 2006) was used to compute a global astrometric solution using the COSMOS i* catalog Capak et al. (2007) as an astrometric reference.

Scamp calculates a global astrometric solution. For each astrometric "context" (defined by the QRUNID image keyword supplied by CFHT), we derive a two-dimensional third-order polynomial which minimizes the weighted quadratic sum of differences in positions between the overlapping sources and the COSMOS reference astrometric catalog derived from interferometric radio observations. Each observation with the same context has the same astrometric solution. Note that we use the STABILITY_MODE "instrument" parameter setting which assumes that the derived polynomial terms are identical exposure to exposure within a given context but allows for anamorphosis induced by atmospheric refraction. This provides reliable, robust astrometric solutions for large numbers of images. The relative positions of each WIRCam detector are supplied by precomputed header file which allows us to use the initial, approximate astrometric solution supplied by CFHT as a first guess.

Our internal and external astrometric accuracies are ∼02 and ∼01, respectively. Scamp produces an XML summary file containing all details of the astrometric solution, and any image showing a large reduced χ² is flagged and rejected.

We do not use Scamp to compute our photometric solutions (the version we used (1.2.11MP) assumes that the relative gains between each WIRCam detector is fixed). Instead, we first use the astrometric solution computed by scamp to match 2MASS stars with objects in each WIRCam image. We then compute the zero point of each WIRCam detector by comparing the fluxes of 2MASS sources with the ones measured by SExtractor in a two-pass process. First, saturated objects brighter than K_AB = 13.84 magnitudes or objects where the combined photometric error between SExtractor and 2MASS is greater than 0.2 mag are rejected. An initial estimate of the zero point is produced by computing the median of the difference between the two cleaned catalogs. Any object where this initial estimation differ by more than 3σ from the median is rejected, and the final difference is computed using an error-weighted mean. The difference between 2MASS and COSMOS catalogs in the stacked image is shown in Figure 1. Note that the main source of scatter at these magnitudes (13.84 < K_AB < 17) comes from uncertainties in 2MASS photometry. The magnitude range with good sources between objects in 2MASS and WIRCam catalogs is quite narrow (around two magnitudes), so setting these parameters is quite important.

Zoom In Zoom Out Reset image size

Finally, all images and weight-maps were combined using Swarp (Bertin et al. 2002). The tangent point used in this paper was 10^h00^m15^s, +02°17'346 (J2000) with a pixel scale of 015 pixel⁻¹ to match COSMOS observations in other filters. The final image has an effective exposure time of one second and a zero point of 31.40 AB magnitudes. This data will be publicly available from the IRSA Web site. The seeing on the final stack is excellent, around 07 FWHM. Thanks to rigorous seeing constraints imposed during queue-scheduled observations, seeing variation over the final stack is small, less than ∼5%.

We also add Subaru-Suprime B_J, i⁺, and z⁺ imaging data (following the notation in Capak et al. 2007). We downloaded the image tiles from IRSA¹⁷ and recombined them with Swarp to produce a single large image astrometrically matched to the K_s-band WIRCam image (the astrometric solutions for the images at IRSA were calculated using the same astrometric reference catalog as our current K_s-image, and they share the same tangent point). Catalogs were extracted using SExtractor (Bertin & Arnouts 1996) in dual-image mode, using the K_s-band image as a detection image. An additional complication arises from the fact that the B-Subaru images saturate at B ∼ 19. To account for this, we use the TERAPIX tool Weightwatcher to create "flag-map" images in which all saturated pixels are indicated (the saturation limit was determined interactively by examining bright stars in the images). During the subsequent scientific analysis, all objects which have flagged pixels are discarded. We also manually masked all bright objects by defining polygon region files. In addition, we also automatically masked regions at fixed intervals from each bright star to remove positive cross talk in the K_s-image. The final catalog covers a total area of 1.9 deg² after masking.

We used SExtractor in dual-image mode with the K_s band as a reference image to extract our catalogs. For this K_s-selected catalog, we measured K_s-band total magnitudes using SExtractor's MAG_AUTO measurement and aperture colors. Our aperture magnitudes are measured in a diameter of 2'' and we compute a correction to "total" magnitudes by comparing the flux of point-like sources in this small aperture with measurements in a larger 6'' diameter aperture. We verified that for the B_J and z⁺ Subaru-Suprime images the difference between these apertures varies less than 0.05 mag, indicating that seeing variations are small across the images which was confirmed by an analysis of the variation of the best-fit Gaussian FWHM for B, z, and K images over the full 2 deg² field.

Obviously for extended bright objects, this color measurement will be dominated by the object nucleus as the majority of the z ∼ 2 objects studied in this paper are unresolved, distant galaxies and we will neglect this effect. We verified that for these objects, variable-aperture colors computed using MAG_AUTO gave results very similar (within 0.1 mag) to these corrected aperture colors.

Based on these considerations, we apply the following corrections to our aperture measurements to compute colors:

$$\begin{equation} i^+_{\rm tot} = i^+-0.1375, \end{equation} \tag{ 1 }$$

$$\begin{equation} K_{\rm tot} = K^{\prime }-0.1568, \end{equation} \tag{ 2 }$$

$$\begin{equation} B_{\rm tot} = B_J - 0.1093. \end{equation} \tag{ 3 }$$

For the z⁺ our corrections are more involved.

As noted in Capak et al. (2007) the Subaru z-band images were taken over several nights with variable seeing. To mitigate the effects of seeing variation on the stacked image PSF individual exposures were smoothed to the same (worst) FWHM with a Gaussian before image combination. This works well at faint magnitudes where many exposures were taken so the non-Gaussian wings of the PSF average out. However, at bright magnitudes (z⁺ ∼ 20) the majority of longer exposures are saturated, so the non-Gaussian wings of the PSF in the few remaining exposures can bias aperture photometry. To correct for this non-linear effect, we apply a magnitude dependent aperture correction in the transition magnitudes between 19 < z < 20. After first applying the correction to total magnitude,

$$\begin{equation} z^+{_{\rm tot}} = z^+ - 0.1093. \end{equation} \tag{ 4 }$$

We apply a further correction, for z⁺ < 19.0, z⁺_tot = z⁺_tot − 0.023; and for z⁺>20.0, z⁺_tot = z⁺_tot + 0.1. For 19 < z⁺ < 20,

$$\begin{equation} z^+_{\rm tot} = z^+_{\rm tot}+(z^+_{\rm tot}-19.0) \times 0.077+0.023. \end{equation} \tag{ 5 }$$

Flux errors in SExtractor are underestimated (in part due to correlated noise in the stacked images) and must be corrected. For a given 4000 ×  4000 image section, we compute a correction factor from the ratio of the 1σ error between a series of 2'' apertures on empty regions of the sky and the median SExtractor errors for all objects. The mean correction factor is derived from several such regions. For B, i, and z images, we multiply our flux errors by 1.5; for the K_s image, we apply a correction factor of 2.0.

We conducted an extensive set of realistic simulations to determine limiting magnitude as a function of object magnitude and profile. In our simulations, we first created a noiseless image containing a realistic mix of stars, disk-dominated and bulge-dominated galaxies using the TERAPIX software stuff and skymaker. Type-dependent luminosity functions are their evolution with redshift are taken from the VVDS survey (Zucca et al. 2006; Ilbert et al. 2006). The spectral energy distribution (SED) of each galaxy type was modeled using empirical templates of Coleman et al. (1980). The disk size distribution was modeled using the fitting formula and parameters presented in de Jong & Lacey (2000).

Next, we used SExtractor to detect all objects on the stacked image (using the same configuration used to detect objects for the real catalogs) and to produce a CHECKIMAGE in which all these objects were removed (keyword CHECKIMAGE_TYPE -OBJECTS). In the next step, this empty background was added to the simulated image, and SExtractor run again in "assoc-mode" in which a match is attempted between each detected galaxy and the output simulated galaxy catalog produced by stuff. In the last step, the magnitude histograms of the number output galaxies is compared to the magnitude histograms in the input catalog; this ratio gives the completeness function of each type of object. In total, around 30,000 objects (galaxy and stars) in one image were used in these simulations. The results from these simulations are shown in Figure 2. The solid line shows the completeness curve for stars and the dotted line for disks. The completeness fraction is 70% for disks and 90% for stars and bulges at K_s ∼ 23.

Zoom In Zoom Out Reset image size

In addition to these simulations, we compute upper limits for each filter based on simple noise statistics in apertures of 2'' (after applying the noise correction factors listed above). The 1σ, 2'' limiting magnitude for our data at the center of the field are 29.1, 27.0, and 25.4 AB magnitudes for B, z, and K_s magnitudes, respectively.

A related issue is the uniformity of the limiting magnitude over the full image. Figure 3 shows the limiting magnitude as a function of position for the K_s stack. This was created by converting the weight-map to an rms error map, scaling this error map from 1σ per pixel to units of 5σ in the effective area of an optimally weighted 07 aperture, and then converting this flux to units of AB magnitude. This depth map agrees well with the completeness limit for point sources shown Figure 2. Note that future COSMOS-WIRCam observations (which will be available in around one year's time) are expected to further reduce the depth variations across the survey area.

Zoom In Zoom Out Reset image size

Based on the considerations outlined in this section, we adopt K_s = 23 as the limiting magnitude for our catalogs. At this limit, our catalog is greater than 90% and 70% complete for point sources and galaxies, respectively. At this magnitude, the number of spurious sources (based on carrying out detections on an image of the K_s stack multiplied by −1) is less than 1% of the total.

One of the principal objectives of this paper is to produce a reliable catalog of objects at z ∼ 2 using a color–color selection technique. A number of different methods now exist to select galaxies in color–color space. For instance, the "dropout" technique (Steidel et al. 1996) makes use of Lyman break spectral feature and the opacity of the high-redshift universe to ultraviolet photons to select star-forming galaxies at z ∼ 3, provided they are not too heavily reddened. Similar techniques can be used at 1 < z < 3 (Adelberger et al. 2004; Erb et al. 2003) and large samples of UV-selected star-forming galaxies now exist at these redshifts. At intermediate redshifts, spectroscopy has shown that "ERO" galaxies which are galaxies selected according to red optical-infrared colors contains a mix of old passive galaxies and dusty star-forming systems in the redshift range 0.8 < z < 2 (Cimatti et al. 2002; Yan et al. 2004). The "DRG" criteria, which selects galaxies with (J − K)_Vega>2.3 is affected by similar difficulties (Papovich et al. 2006; Kriek et al. 2006). On the other hand, the "BzK" criterion introduced by Daddi et al. (2004) can reliably select galaxies in the redshift range 1.4 ≲ z ≲ 2.5 with relatively high completeness and low contamination. Based on the location of star-forming and reddened systems in a spectroscopic control sample lie in the (B − z) and (z − K) plane and considerations of galaxy evolutionary tracks, it has been adopted and tested in several subsequent studies (e.g., Kong et al. 2006; Lane et al. 2007; Hayashi et al. 2007; Blanc et al. 2008; Dunne et al. 2009; Popesso et al. 2009). BzK-selected galaxies are estimated to have masses of ∼10¹¹ M_☉ at z ∼ 2 (Daddi et al. 2004; Kong et al. 2006).

Compared to other color criteria, it offers the advantage of distinguishing between actively star-forming and passively evolving galaxies at intermediate redshifts. It also sharply separates stars from galaxies, and is especially efficient for z>1.4 galaxies. The criterion was originally designed using the redshift evolution in the BzK diagram of various star-forming and passively evolving template galaxies (i.e., synthetic stellar populations) located over a wide redshift interval. Daddi et al. carried out extensive verifications of their selection criteria using spectroscopic redshifts.

To make the comparison possible with previous studies, we wanted our photometric selection criterion to match as closely as possible as the original "BzK" selection proposed in Daddi et al. and adopted by the authors cited above. As our filter set is not the same as this work we applied small offsets (based on the tracks of synthetic stars), following a similar procedure outlined in Kong et al. (2006).

To account for the differences between our Subaru B-filter and the B-VLT filter used by Daddi et al. (2004), we use this empirically derived transformation, defining $bz=B_{J_{\rm total}}-z^+_{\rm tot}$, then for blue objects with bz < 2.5,

$$\begin{equation} bz_{\rm cosmos}=bz+0.0833\times bz+0.053, \end{equation} \tag{ 6 }$$

otherwise, for objects with bz>2.5,

$$\begin{equation} bz_{\rm cosmos}= bz + 0.27. \end{equation} \tag{ 7 }$$

This "bz_cosmos" quantity is the actual corrected (B − z)_AB color, which we use in this paper.

Finally, we divide our catalog into galaxies at z < 1.4, stars, star-forming galaxies, and passively evolving galaxies at 1.4 < z < 2.5, by first defining the BzK quantity introduced in Daddi et al. (2004):

$$\begin{equation} BzK\equiv (z-K)-(B-z). \end{equation} \tag{ 8 }$$

For galaxies expected at z>1.4, star-forming galaxies (hereafter sBzK) are selected as those objects with BzK> − 0.2. One should also note that the reddening vector in the BzK plane is approximately parallel to the sBzK selection criteria, which ensure that the selection is not biased against heavily reddened dusty galaxies.

Old, passively evolving galaxies (hereafter pBzK) can be selected as those objects which have

$$\begin{equation} BzK<-0.2, (z-K)>2.5. \end{equation} \tag{ 9 }$$

Stars are selected using this criteria:

$$\begin{equation} (z-K) < -0.5+(B-z) \times 0.3. \end{equation} \tag{ 10 }$$

Finally, the full galaxy sample consists simply of objects which do not fulfill this stellarity criterion. The result of this division is illustrated in Figure 4. The solid line represents the colors of stars in the BzK filter set of Daddi et al. using the empirically corrected spectra presented in Lejeune et al. (1997), and it agrees very with our corrected stellar locus.

Zoom In Zoom Out Reset image size

Figure 5 shows our differential galaxy number counts compared to a selection of measurements from the literature. We note that at intermediate magnitudes (20 < K_s < 22) counts from the four surveys presented here are remarkably consistent (Elston et al. 2006; Huang et al. 1997; Hartley et al. 2008). At 16 < K_s < 20, discrepancies between different groups concerning measurement of total magnitudes and star–galaxy separation leads to an increased scatter. At these magnitudes, shot noise and large-scale structure begin to dominate the number count errors.

Zoom In Zoom Out Reset image size

The COSMOS-WIRCam survey is currently the only work to provide unbroken coverage over the range 16 < K_s < 23. In addition, our color-selected star–galaxy separation provides a very robust way to reject stars from our faint galaxy sample. These stellar counts are shown by the asterisks in Figure 5. We note that at magnitudes brighter than K_s ∼ 18.0 our stellar number counts become incomplete because of saturation in the Subaru B image (our catalogs exclude any objects with saturated pixels, which preferentially affect point-like sources). Our galaxy and star number counts are reported in Table 2.

Figure 6 shows the counts of star-forming BzK galaxies compared to measurements from the literature. These counts are summarized in Table 2. We note an excellent agreement with the counts in Kong et al. (2006) and the counts presented by the MUYSC collaboration (Blanc et al. 2008). However, the counts presented by the UKIDSS-UDS group (Lane et al. 2007; Hartley et al. 2008) are significantly offset compared to our counts at bright magnitudes, and become consistent with it by K_s ∼ 22. These authors attribute the discrepancy to cosmic variance, but we find photometric offsets a more likely explanation (see below).

Zoom In Zoom Out Reset image size

Figure 7 shows, in more detail, the zone occupied by passive galaxies in Figure 4. Left of the diagonal line are objects classified as star-forming BzK galaxies. Objects not detected in B are plotted as right-pointing arrows with colors computed from the upper limit of their B-magnitudes. An object is considered undetected if the flux in a 2'' aperture is less than the corrected 1σ noise limit. For the B band this corresponds to approximately 29.1 mag. This criterion means that in addition to the galaxies already in the pBzK selection box, fainter sBzK with B-band non-detections (shown with the green arrows) may be scattered rightward into the pBzK region.

Zoom In Zoom Out Reset image size

Counts for our passive galaxy population, including these "additional" objects, are represented by the hatched region in Figure 8. The upper limit for the source counts in this figure represents the case in which all the (z − K)_s>2.5 sources undetected in B are scattered into the pBzK region. Even accounting for these additional objects, we unambiguously observe a flattening and subsequent turnover in the passive galaxy counts at around K_s ∼ 22, well above the completeness limit of either our K_s- or B data in agreement with Hartley et al. (2008).

Zoom In Zoom Out Reset image size

This upper limit, however, is a conservative estimate. We have made a better estimate of this upper limit by carrying out a stacking analysis of the objects not detected in B- in both the passive and star-forming regions of the BzK diagram. For each apparent K_s magnitude bin in table, we median-combine Subaru B-band postage stamps for objects with no B-band detection, producing separate stacks for the star-forming and passive regions of the BzK diagram. In both cases, objects below our detection limit are clearly visible (better than a three-sigma detection) in our stacked images at each magnitude bin to K_s = 23. By assuming that the mean B magnitude of the stacked source to be the average magnitude of our undetected sources, we can compute the average (B − z) color of our undetected sources, and reassign their location in the BzK diagram if necessary. This experiment shows that at most only 15% of the star-forming BzK galaxies undetected in B- move to the passive BzK region.

Our number counts are summarized in Table 3, which also indicates the upper count limits based on B-band observations. As before, our counts are in good agreement with those presented in Kong et al. (2006) and Blanc et al. (2008) but are above the counts in Hartley et al. (2008).

Table 3. Differential Number Counts for the Passive BzK Population

KAB	Passive BzK		Passive BzK (Upper Limits)a
	NpBzK	log (NpBzK) deg−2	NpBzK	$\log\, (N_{\rm {\mathit {pBzK}}}) \,{\rm deg} ^{-2}$
19.25	13	0.84	13	0.84
19.75	69	1.56	69	1.56
20.25	265	2.15	280	2.17
20.75	553	2.47	621	2.52
21.25	837	2.65	1015	2.73
21.75	963	2.71	1285	2.83
22.25	757	2.60	1229	2.81
22.75	475	2.40	984	2.72

Notes. Logarithmic counts are normalized to the effective area of our survey, 1.89 deg². ^aThe upper limit to the pBzK was computed by including all the sBzK galaxies undetected in B with (z − K)>2.5.

Download table as:  ASCIITypeset image

To investigate the origin of this discrepancy, we compared our BzK diagram with Hartley et al.'s, which should also be in the Daddi et al. filter set. We superposed our BzK diagram on that of Hartley et al. and found that the Hartley et al. stellar locus is bluer by ∼0.1 in both (B − z) and (z − K) compared to our measurements.¹⁸ We have already seen that our stellar locus agrees well with the theoretical stellar sequence computed using the Lejeune et al. (1997) synthetic spectra in the Daddi et al. filter set and also with the stellar locus presented in Daddi et al. and Kong et al. We conclude therefore that the number counts discrepancies arise from an incorrect transformation to the Daddi et al. filter set.

In Figures 5, 6, and 8, we show counts of galaxies extracted from the semi-analytical model presented in Kitzbichler & White (2007).

Semi-analytic models start from either an analytic "merger tree" of dark matter halos or, in the case of the model used here, merger trees derived from a numerical simulation, the millennium simulation (Springel et al. 2005). Galaxies are "painted" onto dark matter halos using a variety of analytical recipes which include treatments of gas cooling, star formation, supernovae feedback, and black hole growth by accretion and merging. An important recent advance has been the addition of "radio mode" AGNs feedback (Bower et al. 2006; Croton et al. 2007), which helps provide a better fit to observed galaxy luminosity functions. The Kitzbichler & White model is derived from the work presented by Croton et al. (2006) and further refined by De Lucia & Blaizot (2007). It differs only from these papers in the inclusion of a refined dust model. We refer the readers to these works for further details. An extensive review of semi-analytic modeling techniques can be found in Baugh (2006).

To derive counts of quiescent galaxies, we follow the approach of Daddi et al. (2007) and select all galaxies at z>1.4 in the star formation rate–mass plane (Figure 18 from Daddi et al.), which have star formation rates less than three times the median value for a given mass. Star-forming objects were defined as those galaxies which do not obey this criterion, in this redshift range. (Unfortunately, the publicly available data do not contain all the COSMOS bands so we cannot directly apply the BzK selection criterion to them.)

In all three plots, the models over-predict the number of faint galaxies, an effect already observed for the K-selected samples investigated in Kitzbichler & White. We also note that adding an upper redshift cut to the model catalogs to match our photometric redshift distributions (see later) does not change appreciably the number of predicted galaxies.

Considering, in more detail, the counts of quiescent galaxies, we find that at 20 < K_s < 20.5, models are below observations by a factor of 2; whereas at 22.5 < K_s < 23.0, model counts are in excess of observations by around a factor of 1.5. Given the narrow redshift range of our passive galaxy population, apparent K_s magnitude is a good proxy for absolute K_s magnitude, which can itself be directly related to underlying stellar mass Daddi et al. (2004). This implies that these models predict too many small, low-mass passively evolving galaxies and too few large high mass passively evolving galaxies at z ∼ 1.4.

It is instructive to compare our results with Figure 7 from Kitzbichler & White, which shows the stellar mass function for their models. At z ∼ 2, the models both under-predict the number of massive objects and over-predict the number of less massive objects, an effect mirroring the overabundance of luminous pBzK objects with respect to the Kitzbichler & White model seen in our data.

A similar conclusion was drawn by Fontana et al. (2006) who recently compared predictions for the galaxy stellar mass function for massive galaxies for a variety of models with observations of massive galaxies up to z ∼ 4 in the GOODS field. They also concluded that models incorporating AGN feedback similar to Kitzbichler & White under-predicted the number of high mass galaxies.

Thanks to the wide-area, deep B-band data available in the COSMOS field, we are able to make reliable measurements of the number of faint passive BzK galaxies. Reassuringly, the turnover in counts of passive galaxies observed in our data is qualitatively in agreement with the measurements of the faint end of the mass function of quiescent galaxies at 1.5 < z < 2 made in Ilbert et al. (2009b).

We have cross-correlated our catalog with photometric redshifts to derive redshift selection functions for each photometrically defined galaxy population. Note that although photometric redshifts are based on an optically selected catalog, this catalog is very deep (i' < 26.5) and contains almost all the objects present in the K_s-band selected catalog. At K_s < 23.0, 138,376 were successfully assigned photometric redshifts, representing 96% of the total galaxy population.

Figure 9 shows the redshift distribution for all K_s-selected galaxies, as well as for BzK-selected passively evolving and star-forming galaxies in the magnitude range 18.0 < K_s < 23.0. We have computed the redshift selection function in several magnitude bins and found that the effective redshift z_eff does not depend significantly on apparent magnitude for the sBzK and pBzK populations.

Zoom In Zoom Out Reset image size

By using only the blue grism of the VIMOS spectrograph at the VLT, the zCOSMOS-Deep survey is not designed to target pBzK galaxies, and so no spectroscopic redshifts were available to train the photometric redshifts of objects over the COSMOS field. At these redshifts, the main spectral features of pBzK galaxies, namely Ca ii H & K and the 4000 Å break, have moved to the near infrared. Hence, optical spectroscopy can only deliver redshifts based on identifying the so called Mg-UV feature at around 2800Å in the rest frame. All in all, spectroscopic redshifts of passive galaxies at z>1.4 are now available for only a few dozen objects (Glazebrook et al. 2004; Daddi et al. 2005; Kriek et al. 2006; McGrath et al. 2007; Cimatti et al. 2008). We note that the average spectroscopic redshift of these objects ($\bar{z}\sim 1.7$) indicates that the average photometric redshift of $\bar{z} \sim 1.4$ of our pBzK galaxies to the same K_s-band limit may be systematically underestimated. For the medium to short term, one has to unavoidably photometric redshifts when redshifts are needed for large numbers of passive galaxies at z>1.4.

From Figure 9, we can estimate the relative contribution of each classification type to the total number of galaxies to at least z ∼ 2. This is shown in Figure 10, where we show for each bin in photometric redshift the number of each selection class as a fraction of the total number of galaxies. Upper and lower confidence limits are computed using Poisson statistics and the small-number approximation of Gehrels (1986); in general, we can make reliable measurements to z ∼ 2. In the redshift range 1 < z < 3 at K_s ∼ 22, the pBzK population represents around ∼20% of the total number of galaxies, in contrast to ∼70% for sBzK-selected galaxies. The sum of both components represents at most ∼80% of the total population at z ∼ 2.

Zoom In Zoom Out Reset image size

We estimate the fraction of DRG galaxies using J-band data described in Capak et al. (2007). DRG-selected galaxies remain an important fraction of the total galaxy population, reaching around ∼50% of the total at z ∼ 2. This in contrast with Reddy et al. (2005) who found no significant overlap in the spectroscopic redshift distributions of pBzK and DRG galaxies.

Our work confirms that most bright passive-BzK galaxies lie in a narrower redshift range than either the sBzK or the DRG selection. We also see that the distribution in photometric redshifts for the DRG galaxies is quite broad. Similar conclusions were reached by Grazian et al. (2007) and Daddi et al. (2004) using a much smaller, fainter sample of galaxies.

For each object class, we measure w, the angular correlation function, using the standard Landy & Szalay (1993) estimator:

$$\begin{equation} w (\theta) ={\mbox{DD} - 2\mbox{DR} + \mbox{RR}\over \mbox{RR}}, \end{equation} \tag{ 11 }$$

where DD, DR, and RR are the number of data–data, data–random, and random–random pairs with separations between θ and θ + δθ. These pair counts are appropriately normalized; we typically generate random catalogs with 10 times higher numbers of random points than input galaxies. We compute w at a range of angular separations in logarithmically spaced bins from log (θ) = −3.2 to log (θ) = −0.2 with δlog (θ) = 0.2, where θ is in degrees. At each angular bin, we use bootstrap errors to estimate the errors in w. Although these are not in general a perfect substitute for a full estimate of cosmic variance (e.g., using an ensemble of numerical simulations), they should give the correct magnitude of the uncertainty (Mo et al. 1992).

We use a sorted linked list estimator to minimize the computation time required. The fitted amplitudes quoted in this paper assume a power-law slope for the galaxy correlation function, w(θ) = A_wθ^1−γ; however, this amplitude must be adjusted for the "integral constraint" correction, arising from the need to estimate the mean galaxy density from the sample itself. This can be estimated as (e.g., Adelberger et al. 2005),

$$\begin{equation} C = {1 \over {\Omega ^2}} \int \!\!\! \int w(\theta)\, d\Omega _1\, d\Omega _2. \end{equation} \tag{ 12 }$$

Our quoted fitted amplitudes are corrected for this integral constraint, i.e., we fit

$$\begin{equation} w(\theta)=A_w(\theta ^{1-\gamma }-C). \end{equation} \tag{ 13 }$$

For the COSMOS field, C = 1.42 for γ = 1.8. An added complication is that the integral constraint correction depends weakly on the slope, γ; in fitting simultaneously γ and A_w, we use an interpolated look-up table of values for C in our minimization procedure.

Finally, it should be mentioned that in recent years it has become increasingly clear that the power-law approximation for w(θ) is no longer appropriate (see, for example, Zehavi et al. 2004). In reality, the observed w(θ) is the sum of the contributions of galaxy pairs in separate dark matter halos and within the same halo of dark matter; it is only in a few fortuitous circumstances that this observed w(θ) is well approximated by a power law of slope γ = 1.8. We defer a detailed investigation of the shape of w(θ) in terms of these "halo occupation models" to a second paper, but these points should be borne in mind in the forthcoming analysis.

To verify the stability and homogeneity of our photometric calibration over the full 2 deg² of the COSMOS field, we first compute the correlation function for stellar sources. These stars, primarily residing in the galactic halo, should be unclustered. They are identified as objects below the diagonal line in Figure 4. This classification technique is more robust than the usual size or compactness criterion, which can include unresolved galaxies.

The amplitude of w(θ) as a function of angular scale for stars and faint galaxies is shown in Figure 11. For comparison, we have also plotted the clustering amplitude for our faintest K_s-selected galaxy sample. The inset plot shows a zoom on measurements at large scale scales where the amplitude of w is very low. At each angular bin, our stellar correlation function is consistent with zero out to degree scales down to a limiting magnitude of K_s = 23. If we fit a power-law correlation function of slope 0.8 to our stellar clustering measurements, we find A_w = (1.7  ±   1.7) × 10⁻⁴ (at 1°); in comparison, the faintest galaxy correlation function signal we measure is A_w = (9.9  ±   1.5) × 10⁻⁴, around ∼6 times larger.

Zoom In Zoom Out Reset image size

Figure 12 shows w(θ) for galaxies in three magnitude slices. It is clear that the slope of w becomes shallower at fainter magnitudes. At small separations (less than 1''), w decreases due to object blending. Our fitted correlation amplitudes and slopes for field galaxies are reported in Table 4.

Zoom In Zoom Out Reset image size

In Figure 13, we investigate further the dependence of slope γ on K_s limiting magnitude. Here, we fit for the slope and amplitude simultaneously for all slices. At bright magnitudes, the slope corresponds to the canonical value of ∼1.8; toward intermediate magnitudes, it becomes steeper and fainter magnitudes progressively flatter. It is interesting to compare this figure with the COSMOS optical correlation function presented in Figure 3 of McCracken et al. (2007), which also showed that the slope of the angular correlation function becomes progressively shallower at fainter magnitudes. One possible interpretation of this behavior is that at bright magnitudes our K_s-selected samples are dominated by bright, red galaxies which have an intrinsically steeper correlation function slope; our fainter samples are predominantly bluer, intrinsically fainter objects with shallower intrinsic correlation function slope.

Zoom In Zoom Out Reset image size

Finally, it is instructive to compare our field galaxy clustering amplitudes with literature measurements as our survey is by far the largest at these magnitude limits. Figure 14 shows the scaling of the correlation amplitude at 1° as a function to limiting K_s magnitude, compared a compilation of measurements from the literature. To make this comparison, we have assumed a fixed slope of γ = 1.8 and converted the limiting magnitude of each of our catalogs to Vega magnitudes.

Zoom In Zoom Out Reset image size

In general, our results are within the 1σ error bars of most measurements, although it does appear that the COSMOS field is slightly more clustered than other fields in the literature, as we have discussed previously (McCracken et al. 2007).

In the previous sections, we have demonstrated the reliability of our estimates of w and our general agreement with preceding literature measurements for magnitude-limited samples. We now investigate the clustering properties of passive and star-forming galaxy candidates at z ∼ 2 selected using our BzK diagram. Figure 15 shows the spatial distribution of the pBzK galaxies in our sample; a large amount of small-scale clustering is evident.

Zoom In Zoom Out Reset image size

The upper panel of Figure 16 shows the angular correlation functions for our pBzK, sBzK and for all galaxies. In each case, we apply a 18.0 < K_s < 23.0 magnitude cut. For comparison, we show the clustering amplitude of dark matter computed using the redshift selection functions presented in Section 5 and the non-linear power spectrum approximation given in Smith et al. (2003). At intermediate to large scales, the clustering amplitude of field galaxies and the sBzK population follows very well the underlying dark matter.

Zoom In Zoom Out Reset image size

The lower panel of Figure 16 shows the bias b, as a function of scale, computed simply as $b(\theta)=\sqrt{(w_{\rm gal}(\theta)/w_{\rm dm}(\theta)}$. Dashed, dotted, and solid lines show b values for pBzK, sBzK, and field galaxies retrospectively (in this case, our w measurements have been corrected for the integral constraint). The bias for the faint field galaxy population is 1.2 at 1' indicating that the faint K_s-selected galaxy population traces well the underlying dark matter. In comparison, at the same scales, the bias values for the passive BzK and star-forming BzK galaxies are 2.5 and 2.1, respectively.

Our best-fitting γ and amplitudes (quoted at 1') for pBzK and sBzK galaxies are reported in Table 5. Given that for the sBzK galaxies γ = 1.8, we may compare with previous authors who generally assume a fixed slope γ = 1.8 for all measurements. At K_VEGA < 20, corresponding to K_AB ∼ 22, Kong et al. (2006) find (4.95 ± 0.52) × 10⁻³ whereas (at 1°) we measure (2.1  ±   0.6) × 10⁻³, closer to the value of (3.14  ±   1.12) × 10⁻³ found by Blanc et al. (2008). We note that both Hayashi et al. (2007) and Adelberger et al. (2005) also investigated the luminosity dependence of galaxy clustering at z∼ although with samples considerable smaller than those presented here. It is plausible that field-to-field variation and large-scale structure are the cause of the discrepancy between these surveys.

The best-fitting slopes for our pBzK populations is γ ∼ 2.3, considerably steeper than the field galaxy population (no previous works have attempted to fit both slope and amplitude simultaneously for the pBzK populations due to small sample sizes). In the next section, we will derive the spatial clustering properties of both populations.

To de-project our measured clustering amplitudes and calculate the comoving correlation lengths at the effective redshifts of our survey slices, we use the photometric redshift distributions presented in Section 5.

Given a redshift interval z₁, z₂, and a redshift distribution dN/dz, we define the effective redshift in the usual way, namely, z_eff is defined as

$$\begin{equation} z_{\rm eff} ={\int _{z_1}^{z_2}z({\it dN}/dz)dz/{\int _{z_1}^{z_2}({\it dN}/dz)dz}}. \end{equation} \tag{ 14 }$$

Using these redshift distributions together with the fitted correlation amplitudes in presented in Sections 6.2 and 6.3, we can derive the comoving correlation lengths r₀ of each galaxy population at their effective redshifts using the usual Limber (1953) and Peebles (1980) inversion. We assume that r₀ does not change over the redshift interval probed.

It is clear that our use of photometric redshifts introduces an additional uncertainty in r₀. We attempted to estimate this uncertainty by using the probability distribution functions associated with each photometric redshift to compute an ensemble of r₀ values, each estimated with a different n(z). The resulting error in r₀ from these many realizations is actually quite small, ∼0.02 for the pBzK population. Of course, systematic errors in the photometric redshifts could well be much higher than this. Figure 9 in Ilbert et al. shows the 1σ error in the photometric redshifts as a function of magnitude and redshift. Although all galaxy types are combined here, we can see that the approximate 1σ error in the photometric redshifts between 1 < z < 2 is ∼0.1. Our estimate of the correlation length is primarily sensitive to the median redshift and the width of the correlation length. An error ∼0.1 translates into an error of ∼0.1 in r₀. We conclude that, for our pBzK and sBzK measurements, the dominant source of uncertainty in our measurements of r₀ comes from our errors on w.

We note that previous investigations of the correlation of passive galaxies always assumed a fixed γ = 1.8; from Figure 12, it is clear that our slope is much steeper. These surveys, however, fitted over a smaller range of angular scales and therefore could not make an accurate determination of the slope for the pBzK population. In all cases, we fit for both γ and A_w.

Our spatial correlation amplitudes for pBzK and sBzK galaxies are summarized in Table 6. Because of the degeneracy between r₀ and γ, we also quote clustering measurements as r^γ/1.8₀. These measurements are plotted in Figure 17. At lower redshifts, our field galaxy samples are in good agreement with measurements for optically selected redder galaxies from the CFHTLS and VVDS surveys (Meneux et al. 2006; McCracken et al. 2008). At higher redshifts, our clustering measurements for pBzK and sBzK galaxies are in approximate agreement with the measurements of Blanc et al. (2008). We note that part of the differences with the measurements of Blanc et al. arises from the their approximation of the redshift distribution of passive BzK galaxies using simple Gaussian distribution.

Zoom In Zoom Out Reset image size

Interestingly, a steep slope γ for optically selected passive galaxies has already been reported at lower redshift surveys; for example, Madgwick et al. (2003) found that passive galaxies had a much steeper slope than active galaxies in the 2dF galaxy redshift survey.

The highly biased nature of the pBzK galaxy population indicates that these objects reside in more massive dark matter halos than either the field galaxy population or the sBzK population, and we intend to present a more detailed discussion of the spatial clustering of each galaxy sample in the framework of the halo occupation models in a future paper.