A LABOCA SURVEY OF THE EXTENDED CHANDRA DEEP FIELD SOUTH—SUBMILLIMETER PROPERTIES OF NEAR-INFRARED SELECTED GALAXIES

In this paper, we aim to determine the contribution to the submm background from K_vega ⩽ 20 galaxies, as well as the sub-samples of BzKs, EROs, and DRGs and their joint contributions, and it is therefore important to determine the degree of overlap between these populations. The overlaps in terms of percentages are given in Table 1.

From Figure 2 it is seen that BzK, ERO, and DRG galaxies only start to contribute significantly (>1%) to the full sample for K_AB ≳ 20. Even so, more than half of the full sample does not fall within the BzK/ERO/DRG classifications at these faint flux levels. Of the full K_vega ⩽ 20 sample, 6269 sources (corresponding to ∼76%) do not classify as BzKs, EROs, or DRGs.

Zoom In Zoom Out Reset image size

The locations of the various samples in the BzK and RJK diagrams are shown in Figure 1. The DRGs are seen to be spread out across the BzK diagram, while the EROs lie in a much more well-defined region of the BzK diagram. EROs and DRGs make up about 30% and 36% of the sBzKs, respectively. Similarly, we find that EROs and DRGs constitute 98% and 44% of the pBzK sample, respectively. Clearly, pBzK galaxies are much more often selected as EROs than is the case for sBzKs, while the occurrence of a pBzK or sBzK galaxy being classified as a DRG is about the same. Turning to the RJK diagram we see that the sBzK galaxies are much more spread out than the pBzKs. About 18% and 37% of EROs and DRGs, respectively, are made up of sBzK, while pBzKs make up about 12% and 9% of the two populations.

The overlaps between BzK, ERO, and DRG galaxies have been discussed in detail in other studies (e.g., Reddy et al. 2005; Grazian et al. 2007; Takagi et al. 2007; Lane et al. 2007). The most statistically significant study was carried out by Lane et al. (2007) who used the UKIRT Infrared Deep Sky Survey (UKIDSS) Ultra Deep Survey Early Data Release (UDS EDR) to study large samples of BzK, ERO, and DRG galaxies. For samples selected down to K_AB = 21.2, which is close to our magnitude cutoff, they found sBzK:ERO and pBzK:ERO ratios of 32% and 95%, respectively, i.e., in excellent agreement with our values. They also find that about 30% of DRGs are sBzK, again in good agreement with our findings (see also Reddy et al. 2005).

The sample was correlated against publicly available spectroscopic redshift surveys (Szokoly et al. 2004; Vanzella et al. 2005, 2006, 2008; Popesso et al. 2009; Kriek et al. 2008). Using only the most reliable spectroscopic redshifts from these surveys, we extracted 2341 redshifts, the majority of which lie within the CDF-S region. A total of 546 galaxies from our sample were matched to a spectroscopic redshift. Of these, 28 were sBzK, 21 were ERO, and 10 were DRG galaxies. The fact that no pBzK galaxies were identified with a spectroscopic redshift is not too surprising since these are in all likelihood old, evolved galaxies with optical spectra devoid of emission features, thus making it difficult to obtain robust spectroscopic redshifts (cf. Kriek et al. 2006; Cimatti et al. 2008). Sources with spectroscopic redshifts were used as a test sample to optimize input parameters for the photometric redshift code EAZY (Brammer et al. 2008). The code works by fitting non-negative linear combinations of galaxy spectra to the observed SEDs, which in our case consisted of nine MUSYC filter bands.

The resulting photometric redshifts are compared against their spectroscopic counterparts in Figure 3. The normalized median absolute deviation of Δz = z_phot − z_spec (see Brammer et al. 2008) is ≃0.037 for z ⩽ 1.5 and ≃0.079 for z > 1.5. Significant outliers, which we define to be sources with |Δz|/(1 + z_spec) five times greater than the median, make up ∼9% of the total sample. These numbers are consistent with the typical performances of photometric redshift codes (e.g., Bolzonella et al. 2000; Quadri et al. 2007; Brammer et al. 2008).

Zoom In Zoom Out Reset image size

Adopting the parameters for the test sample, photometric redshifts were derived for the remainder of the sample without spectroscopic redshifts. The redshift distributions of the full K_vega ⩽ 20 sample as well as the BzK, ERO, and DRG samples are shown in Figure 4, where we also compare with two other photometric redshift distributions obtained from: (1) the MUSYC UBVRI and deep JHK imaging of the three 10' × 10' fields HDFS1/S2, MUSYC 1030 and 1255 (Quadri et al. 2007, hereafter Q07), and (2) deep BVRi'z'JK imaging of 1113 arcmin² in the Ultra-Deep Survey (UDS) portion of the United Kingdom Infrared Telescope Deep Sky Survey (UKIDSS; Dunne et al. 2009, hereafter D09).

Zoom In Zoom Out Reset image size

The BzK selection criteria are defined to select galaxies in the redshifts range 1.4 < z < 2.5 (Daddi et al. 2004), yet both the sBzK and pBzK redshift distributions extend below and above this range. Of the sBzK galaxies, 66% lie in the range 1.4 < z < 2.5, while 22% are at z < 1.4 and 12% at z > 2.5. For the pBzK galaxies, the corresponding percentages are 60%, 31%, and 9%, respectively. Both Q07 and D09 find similar fractions for their samples of K_vega ⩽ 20 galaxies, suggesting that the BzK criteria select galaxies across a somewhat broader redshift range (1 ≲ z ≲ 3.5). The redshift distribution of EROs is seen to peak strongly at z ≃ 1.1 with a tail extending to z ∼ 3.5. This is in line with photometric and spectroscopic surveys which have shown that the redshift distribution of EROs peaks around z ≃ 1.2 (Cimatti et al. 2003; Yan et al. 2004). The redshift distribution of K_vega ⩽ 20 EROs derived by D09 peaks at slightly higher redshifts (z ∼ 1.4), but overall appears similar to ours. The DRG distribution shows prominent peaks at z ≃ 1.2 and z ≃ 2, with the former being the most dominant. A similar bimodality is also apparent in the K_vega ⩽ 20 DRG sample by Q07, although the dominant peak in their distribution lies at z ≃ 2. The distribution by D09 broadly resembles ours, with a prominent peak at z ≃ 1.2 followed by a significant high-z tail. Overall, therefore, our redshift distributions are consistent with those of Q07 and D09, given the uncertainties associated with photometric redshift derivation, and the effects of field-to-field variations.

First, however, we need to identify any galaxies that are associated with robust LABOCA sources, which we take to be mean sources detected at ⩾3.7σ significance (126 in total; see Weiß et al. 2009 for details). To this end, we adopt a search radius of 128 around each LABOCA source, which corresponds to the 95% confidence search radius given the FWHM = 19'' beam. For a galaxy to qualify as a near-IR counterpart to a LABOCA source we furthermore require that z > 0.8. This cut is motivated by the observed redshift distribution of submm galaxies which shows that ≲3% are at z < 0.8 (Chapman et al. 2005). If more than one galaxy meets these criteria for a given LABOCA source, we adopt the one closest to the submm source. In this manner, we find 57 submm-near-IR associations from the K_vega ⩽ 20 sample. Of these, 19/2 sources are classified as sBzK/pBzK galaxies and 25 as EROs (of which 13 are also DRGs). Ten of the nineteen sBzK galaxies are also EROs and of those, nine are DRGs. Both pBzK galaxies are EROs and also DRGs. As a safeguard against contamination, these submm-bright near-IR sources were removed from the stacking analysis, although their contribution was included (in a variance-weighted fashion) in the final tally of average submm flux.

Next, we proceeded to perform a stacking analysis of the remaining submm-undetected near-IR selected galaxies. Due to the slightly varying noise across the map, the average 870 μm flux (〈S_870 μm〉) and noise (〈σ_870 μm〉)) values were calculated as the variance-weighted mean, i.e.,

$$\begin{equation} \langle S_{{\rm {\,870\,\mu {\rm m}}}}\rangle = \frac{\sum _i S_i/\sigma _i^2}{\sum _i 1/\sigma _i^2} \end{equation} \tag{ 1 }$$

and

$$\begin{equation} \langle \sigma _{{\rm {\,870\,\mu {\rm m}}}}\rangle = \frac{1}{\sqrt{\sum _i 1/\sigma _i^2}}, \end{equation} \tag{ 2 }$$

where S_i and σ_i are the 870 μm flux and rms noise pixel values at the near-IR position of the ith source in the stack, respectively. To avoid the stack being contaminated by robust 870 μm sources, the stacking was performed on a "residual" version of the LABOCA map in which all of the 126 870-μm sources uncovered at ⩾3.7σ (Weiß et al. 2009) had been subtracted using a scaled beam profile.

An important aspect of any stacking analysis performed on maps with coarse angular resolution is the issue deblending sources that lie within a single resolution element. For example, if a BzK galaxy has a neighbor, A, within a LABOCA beam, we have to calculate the 870 μm flux contribution from A to the position of the BzK galaxy. Now, if A also has a neighbor, B, within a LABOCA beam (which is not necessarily a neighbor to the BzK galaxy), then the contribution from B to A will affect A's contribution to the BzK galaxy and will have to be included in the deblending. If B has no other neighbors other than A, then the process stops there, but if B has another neighbor, C, its contribution will also have to be factored into the deblending, and so on. In this way, "chains" of neighboring sources are constructed for every source, and the entire "chain" of neighbors must be included in the deblending procedure.

Of the full K_vega ⩽ 20 sample, 5985 sources (i.e., 72% of the sample) are involved in such a "chain" of two or more galaxies, and thus must be deblended. We correct for the blending of sources by assuming Gaussian sources with FWHMs equal to the LABOCA beam. This method is similar to the one adopted by Webb et al. (2004), although they (and subsequent submm stacking studies) did not take into account the effects from neighbors' neighbors (and so on) as discussed above. Our study is the first stacking analysis that has adopted this "global" deblending scheme. A simple illustrative example of the latter where only four sources are involved is given Figure 5, where A itself has a neighbor, B, within a LABOCA beam, which in turn has a neighbor C. Neither B nor C is within the LABOCA beam as measured from the BzK galaxy's position. In this case, we have to calculate C's contribution to B, and B's contribution to A, in order to correctly calculate A's contribution to the BzK galaxy, and the system of equations to be solved is therefore

$$\begin{equation} f_{{\rm {BzK}}} = I_{{\rm {BzK}}} + I_{{\rm {A}}} e^{-r_{{\rm {BzK,A}}}^2/(2\sigma ^2)}, \end{equation} \tag{ 3 }$$

$$\begin{equation} f_{{\rm {A}}} = I_{{\rm {A}}} + I_{{\rm {BzK}}} e^{-r_{{\rm {BzK,A}}}^2/(2\sigma ^2)} +I_{{\rm {B}}} e^{-r_{{\rm {B,A}}}^2/(2\sigma ^2)}, \end{equation} \tag{ 4 }$$

$$\begin{equation} f_{{\rm {B}}} = I_{{\rm {B}}} + I_{{\rm {A}}} e^{-r_{{\rm {B,A}}}^2/(2\sigma ^2)} + I_{{\rm {C}}} e^{-r_{{\rm {C,A}}}^2/(2\sigma ^2)}, \end{equation} \tag{ 5 }$$

$$\begin{equation} f_{{\rm {C}}} = I_{{\rm {C}}} + I_{{\rm {B}}} e^{-r_{{\rm {B,C}}}^2/(2\sigma ^2)}, \end{equation} \tag{ 6 }$$

where I and f are the measured and deblended fluxes at the relevant positions, respectively, and r are the distances between the sources. In order to estimate the error one makes by only deblending the fluxes from neighbors within a beam and not taking into account neighbors' neighbors etc., we ran the stacking analysis under both scenarios. We find that the deblending scheme by Webb et al. (2004) can overestimate the fluxes by ∼10% compared to the scheme described in this paper, but more typically the error is at the ∼5% level. We emphasize that one has to not only deblend neighbors within a certain galaxy population, but also across populations. Therefore, the deblending analysis was carried out on the full K_vega ⩽ 20 sample. In doing so we have ignored sources with K_vega > 20, and they have therefore not been included in the deblending analysis. In the following section, however, we show that these fainter sources do indeed make a contribution to the submm signal, which must be subtracted from the stacked fluxes. For an extensive comparative analysis of different stacking techniques (including the one presented here) and Monte Carlo simulations that demonstrate the convergence of our deblending method with those of other deblending/stacking methods, we refer to Kurczynski & Gawiser (2010).

Zoom In Zoom Out Reset image size

For each source in the K_vega ⩽ 20 sample the (deblended) signal and noise values at its pixel position in the LABOCA map were recorded, and from those the stacked 870 μm flux density of the full sample was determined according to Equations (1) and (2). The 870 μm signal and noise values corresponding to the BzK, ERO, and DRG samples were extracted and their stacked 870 μm flux densities were derived in a similar manner. Also, postage stamp images around each source were extracted and combined in a weighted fashion resulting in stacked submm images of the K-selected samples (Figure 6). From the azimuthally averaged radial profiles of the submm signal, it is clear that the baseline level is not zero, but in fact there is a residual signal amounting to 0.065 mJy. While this baseline level probably stems from the population of the K_vega > 20 galaxies lying below the submm detection limit of the LABOCA map, we stress that the residual signal will be affected by the subtraction of large-scale structure in the map described in Section 2. In order to account for the residual signal, the stacked 870 μm flux densities had a constant signal of 0.065 mJy subtracted from them. The final stacked 870 μm flux densities (with the 0.065 mJy baseline level subtracted) are given in Table 2. As a comparison, we also derived the median flux densities from the stacks and found agreement (to within 15%) with the weighted averages.

Zoom In Zoom Out Reset image size

We find that all of our K-selected samples have very significantly (≳8σ) detected stacked 870 μm fluxes, except for the pBzKs, which are detected only at the ∼3σ level. To gauge the significance of our results, we ran a series of Monte Carlo simulations in which stacking analyses were carried out on 1000 versions of the K_vega ⩽ 20 catalog, each with randomized positions with respect to the original catalog. Each source was assigned a random position by (randomly) choosing a radius (60'' ⩽ r ⩽ 200'') and an angle from its original position. By confining the new positions to within a certain distance of the original, we ensured that the noise properties were similar to those in the original stacking analysis. As expected, the distributions of stacked signals obtained from these simulations were Gaussians centered on zero. We found that the measured 870 μm signals occurred in <0.05% of the simulation runs. Roughly, the same percentage is found if we restrict our Monte Carlo analysis to the sBzK, pBzK, ERO, and DRG samples, and clearly testify to the significance of the measured signals for these subsets.

From Table 2 it is seen that while the significance of the pBzK stacked signal is only ∼3σ, the actual signal (0.28  ±   0.10 mJy) is in fact identical to that of the EROs (0.29 ± 0.03 mJy). It is possible, therefore, that the low significance of the stacked pBzK signal is partly due to the smaller sample size (and hence higher noise in the stack) rather than pBzKs being intrinsically devoid of submm emission. In order to test this we extracted random subsets of 147 EROs and stacked their fluxes. Doing so 1000 times, we found that in 10% of the cases a flux equal or smaller than that of the pBzKs was obtained. The average flux over the 1000 runs was 0.29 ± 0.10 mJy (S/N = 2.9), which is virtually identical to the stacking signal obtained for the pBzK galaxies. The fact that there is a 10% chance of obtaining a signal corresponding to that of the pBzK by randomly stacking 147 sBzKs suggest that the low-significance of the pBzK stack is at least in part due to the small sample size, and that the pBzK galaxies in fact have a non-zero submm signal.

In addition to the above stacking analysis, we compared the distributions of S/N values at the positions of the submm-faint, near-IR sources with the overall S/N distribution of the residual LABOCA map. From Figure 7, it is seen that the distributions for the K_vega ⩽ 20, sBzK, pBzK, ERO, and DRG samples all appear to be biased toward positive S/N values, although a formal Kolmogorov–Smirnov (K-S) statistic suggests that in the case of the full K_vega ⩽ 20 sample, this bias is not statistically significant. Interestingly, the pBzK distribution shows an excess toward positive S/N values, which seem to be statistically significant (P_KS = 0.003). This is consistent with the above finding of a significant (at the ∼3σ level) stacked 870 μm signal from pBzKs.

Zoom In Zoom Out Reset image size

In Section 4, we found that all of our K-selected samples have highly significant (⩾8σ) stacked 870 μm fluxes, except for the pBzKs which were only marginally detected (∼3σ). How do the stacked fluxes derived in Section 4 compare with previous submm stacking analyses of K-selected samples?

We find an average 870 μm flux density of 0.37 ± 0.04 mJy for our sample of submm-faint sBzK galaxies. In comparison, Daddi et al. (2005) estimated an average 850 μm signal of ∼0.82 mJy (corresponding to a 870 μm signal of ∼0.75 mJy¹⁷) for a sample of ∼100 K_vega ⩽ 20 submm-faint sBzK galaxies within the SCUBA map of GOODS-N. Takagi et al. (2007) reported an average 850 μm flux of 0.52 ± 0.19 mJy (corresponding to a 870 μm flux of 0.48 ± 0.18 mJy) from a sample of 112 K_vega ⩽ 20 sBzK galaxies selected within the part of the Subaru/XMM-Newton deep field (SXDF) covered by the SCUBA HAlf Degree Extragalactic Survey (SHADES). Finally, D09 stacked the 850 μm signal from 1421 K_vega ≲ 21.7 sBzK galaxies and found 0.53 ± 0.06 mJy (corresponding to 0.49 ± 0.06 mJy at 870 μm).

For the 147 submm-faint pBzK galaxies we derived an average 870 μm flux density of 0.28 ± 0.10 mJy. As discussed in Section 4.3, there is statistical evidence of a positive submm signal from the pBzK sample, and the lower significance of the signal can be ascribed to a smaller intrinsic signal and sample size. We note that D09 measured an average 850 μm flux of 0.22 ± 0.18 mJy (or 0.20 ± 0.17 mJy at 870 μm) for a sample of 147 K_vega ≲ 21.7 pBzK galaxies, which, although not statistically significant, is in agreement with our result.

For the EROs and DRGs, we find stacked 870 μm fluxes of 0.29 ± 0.03 mJy and 0.32 ± 0.04 mJy, respectively. Webb et al. (2004) used the Canada–UK Deep Submillimeter Survey 03 hr and 14 hr fields (CUDSS03 and CUDSS14, respectively) to perform an 850 μm stacking analysis of 164 K_vega ≲ 20.7 EROs in the two fields, and found a stacked signal for the entire ERO sample of 〈S_870 μm〉 = 0.52 ± 0.09 mJy. Similarly, Takagi et al. (2007) measured a stacked 870 μm signal of 0.49 ± 0.16 mJy from a sample of 201 K_vega ⩽ 20 EROs selected within SXDF/SHADES. Turning to DRGs, Knudsen et al. (2005) obtained a stacked 850 μm signal of 0.74 ± 0.24 mJy from a sample of 25 K_vega ⩽ 22.5, submm-faint DRGs (uncorrected for an average gravitational lens amplification of 20%). Converting to a 870 μm flux density and correcting for the lens amplification yields 0.54 ± 0.18 mJy. Takagi et al. (2007), using significantly shallower SCUBA maps, failed at detecting a significant 850 μm signal from an average of 67 K_vega ⩽ 20, submm-faint DRGs (〈S_870μm〉 = 0.39 ± 0.23 mJy).

We conclude that, within the errors, the previous stacking studies agree well with our results. We also note that our study provides the first robust (∼8σ) detection of submm-faint DRGs.

From the stacked submm fluxes we are able to estimate average IR luminosities and SFRs (Table 3). IR luminosities are derived by adopting the IR-to-submm SED of Arp 220 and scaling it to the stacked submm fluxes (at the median redshifts derived from the redshift distributions in Section 3.2) and integrating it from 8 to 1000 μm. For comparison, we also derive IR luminosities assuming that the SEDs are described by modified blackbody law with a dust temperature of T_d = 30 K and β = 1.5, and integrating from 8 to 1000 μm. SFRs are derived following Kennicutt (1998): SFR[M_☉ yr⁻¹] = 1.73 × 10⁻¹⁰L_IR[L_☉]. This conversion assumes a Salpeter initial mass function (Salpeter 1955). Of course, we stress that considerable uncertainty is associated with the derived IR luminosities and SFRs since they depend on the assumed SED and IMF.

The average IR luminosities and SFRs estimated here for the K_vega ⩽ 20 sBzK, ERO, and DRG galaxies on the basis of their stacked submm fluxes lie in the ranges ∼(1–6) × 10¹¹ L_☉  and ∼20–110 M_☉  yr⁻¹, i.e., comparable to those of luminous infrared galaxies (LIRGs; L_IR ∼ 10¹¹ L_☉ —Sanders & Mirabel 1996) in the local universe. For the sBzK galaxies, the average IR luminosity and SFR derived here are fully consistent with UV, 24 μm, and radio studies of these galaxies (Daddi et al. 2007). We find that the ERO and DRG populations have significantly lower IR luminosities (by ∼40%) than the sBzK galaxies. This is in part due to the fact that we have made no attempt to weed out passive EROs/DRGs in the stacking analysis, and the stacked submm flux from dusty, star-forming EROs/DRGs is likely to be significantly higher.

Using the photometric redshifts obtained in Section 3.2 we can stack our samples into separate redshift bins, thereby allowing us to determine which redshifts are contributing the most to the stacked submm signals. Redshift bins were chosen so that they were larger than the typical redshift uncertainty, and provided roughly the same number of sources in each bin. The latter ensured that the same sensitivity was reached in each bin, thus allowing for a direct comparison. Figure 8 shows the average flux densities of the different samples as a function of redshift. We stress that due to the essentially flat selection function of submm surveys over the redshift range 1 ≲ z ≲ 8 (Blain & Longair 1993), the comparison of stacked submm fluxes at different redshifts directly translates into a comparison between far-IR luminosities (and thus SFRs—Kennicutt 1998) between these redshift bins (to the extent that submm flux is a measure of far-IR luminosity).

Zoom In Zoom Out Reset image size

The average submm signal from K_vega ⩽ 20 galaxies is found to be roughly constant (∼0.1–0.2 mJy) out to z ∼ 1.4, and consistent with the average flux density of the full sample (Table 2). At z ∼ 1.7, however, the submm signal has increased to ∼0.4 mJy. The stacked submm fluxes of z < 1.4 and z > 1.4 K_vega ⩽ 20 sources are 0.18 ± 0.01 mJy and 0.45 ± 0.03 mJy, respectively (Table 2). A similar increase in the average submm flux at z > 1.4 is seen for the ERO and DRG populations. The average IR luminosities and SFR for z > 1.4 EROs (L_IR ≃ 5.5 × 10¹¹ L_☉  and SFR ≃ 96 M_☉  yr⁻¹) are about 2× higher than for z < 1.4 EROs (L_IR ≃ 2.5 × 10¹¹ L_☉  and SFR ≃ 43 M_☉  yr⁻¹), where we have adopted the IR luminosities and SFRs derived from the Arp 220 SED. A similar trend is seen for DRGs: z > 1.4 DRGs have on average L_IR ≃ 5.5 × 10¹¹ L_☉  and SFR ≃ 94 M_☉  yr⁻¹, while z < 1.4 DRGs have L_IR ≃ 2.7 × 10¹¹ L_☉  and SFR ≃ 46 M_☉  yr⁻¹. These findings also fit with the ERO and DRG redshift distributions (Figure 4), where we found evidence for two sub-populations separated at z ∼ 1.6. Due to our magnitude cutoff at K_vega ⩽ 20, our samples are biased toward increasingly more massive galaxies at higher redshifts. The strong submm signal indicates that a substantial fraction of massive EROs/DRGs at z > 1.4 are actively star-forming galaxies and not old, red galaxies. The gradual drop in the average submm signals at z > 2.5 could be indicative of a decline in the star formation activity in EROs and DRGs at these higher redshifts or, alternatively, it could simply reflect a bias toward less dusty (and thus star-forming) galaxies since at redshifts >2.5 K-band observations start sampling the optical (rest-frame) emission.

Turning to the BzK galaxies, we find that sBzKs exhibit a positive and, within the error bars, constant 870 μm signal (∼0.4 mJy) over the redshift range 1 ≲ z ≲ 3. This supports the notion that the sBzK criterion selects star-forming galaxies across this redshift range, and that for our K-band magnitude limit of K_vega ⩽ 20, the distribution of SFRs of sBzK galaxies is roughly constant within this redshift range. In contrast, the pBzK galaxies show no evidence of significant submm signal in any redshift bin, which is consistent with these galaxies being devoid of significant star formation. D09 reported a 850 μm signal of 0.89 ±  0.34 mJy (∼2.6σ) for pBzKs at z < 1.4 but no significant signal for pBzKs at z > 1.4. They argued that the submm signal for z < 1.4 pBzK galaxies was due to contamination by star-forming galaxies at z < 1.4. In comparison, we find no significant stacked 870 μm signal (0.04 ± 0.17 mJy) for z < 1.4 pBzKs and only a marginal signal for z > 1.4 pBzKs (0.38 ± 0.12 mJy).

Using the 24 μm source catalog from the FIDEL survey (which includes all sources >27 μJy, 5σ point-source sensitivity; Dickinson et al. 2007), a total of 466/53, 511, and 400 sources in the sBzK/pBzK, ERO, and DRG samples, respectively, were identified at 24 μm. Their average submm fluxes are given in Table 2 along with those of the 24 μm faint (<27 μJy) subsets. The 24 μm detected subsets have ≳5× higher average submm flux densities than the 24 μm faint sources. This is not surprising since the mid-IR is known to trace the thermal dust emission, and one might even expect a correlation between the submm and mid-IR emission.

In order to investigate the latter, we measured the stacked 870 μm signal as a function of 24 μm flux density bins (Figure 9). The sBzK, ERO, and DRG galaxies all exhibit a similar behavior, namely a significant linear correlation between the stacked 870 μm and 24 μm flux density up to S_24 μm ≃ 350 μJy, given by 〈S_870 μm〉 = 4.5 × 10⁻³〈S_24 μm〉, followed by a flattening of the relation for S_24 μm > 350 μJy. Not surprisingly, the pBzKs are not significantly detected at 870 μm in any of the 24 μm bins, and no submm-mid-IR correlation is seen.

Zoom In Zoom Out Reset image size

Although the 24 μm selection function (Daddi et al. 2007) depends strongly on redshift (unlike the selection function at submm wavelengths), the observed S_24 μm–S_870 μm correlation strongly suggests that 24 μm measurements (with S_24 μmJy ≳ 350 μJy) trace systems dominated by star formation, and can be used to derive IR luminosities and SFRs. The constant S_870 μm ≃ 1.3 mJy ratio for S_24 μm ≃ 350–1000 μJy either implies that these bright 24 μm sources are significantly contaminated by active galactic nucleus (AGN), which would boost the 24 μm signal but leave the 870 μm flux relatively unchanged, or that they are nearby (z < 0.5) sources in which case we would not expect them to exhibit a strong submm signal (due to the lack of a significant k-correction). The latter possibility can be ruled out, however, as a negligible fraction of the sources in our sample with S_24 μm ⩾ 250 μJy lie at z < 0.8. In fact, omitting these sources from the analysis does not change the S_24 μm–S_870 μm correlation in any significant way.

It therefore seems most likely that the flattening of the correlation is due to the AGN contributing more prominently to the 24 μm flux in the brightest sources. Since the mid-IR will be more sensitive to warm dust (T_d ∼ 80–200 K) than the submm, which largely traces cold dust (T_d ∼ 20–60 K), the ability of the observed 24 μm emission to reliably trace IR luminosity and star formation may be compromised by AGN-heated warm dust. The turnover from a linear to a flat S_870 μm–S_24 μm relation occurs at S_24 μm ≃ 350 μJy, which corresponds to L_IR < 1.5 × 10¹² L_☉  at z ∼ 2.

Our findings are in line with those of Papovich et al. (2007), who found that IR luminosities derived from 24, 70, and 160 μm Spitzer data (of a sample of K-selected galaxies at 1.5 ≲ z ≲ 2.5 with S_24 μmJy ≃ 50–250 μJy) were in good agreement with those derived from the 24 μm data alone. For sources with S_24 μmJy > 250 μJy, however, the latter would be ∼2–10× higher, suggesting that the AGN may contribute significantly to the high 24 μm emission. Similarly, Daddi et al. (2004, 2005) found that a large fraction (>30%) of sBzK galaxies with L_IR ≳ 1.5 × 10¹² L_☉  (which corresponds to the IR luminosity where we find a turnover in the S_870 μm–S_24 μm relation) show an excess of emission in the near-IR (rest frame) and are statistically detected in hard X-rays—evidence of powerful AGNs in these very IR-luminous systems.

Turning to the contribution to the submm EBL by the different populations, we adopt the spectral approximations to the EBL at submm wavelengths from COBE/FIRAS (Puget et al. 1996; Fixsen et al. 1998), including their uncertainties. In doing so we adopt a value of the EBL at 870 μm of 44 ± 15 Jy deg⁻². From the surface densities of the samples we derive their contributions to the 870 μm EBL (Table 2). We find that the total contribution from all K_vega ⩽ 20 sources to the 870 μm EBL is 7.27 ± 0.34 Jy deg⁻² or 16.5% ± 5.7%. For the BzK galaxies we find that the passive ones contribute ≲1% to the 870 μm EBL, while the star-forming population contribute 1.49 ±  0.22 Jy deg⁻², corresponding to 3.4% ± 1.3%. The EROs contribute 1.95 ± 0.21 Jy deg⁻² (or 4.4% ± 1.6%), while the contribution from DRGs, owing to their lower surface density, is about 25% smaller (1.27 ± 0.20 Jy deg⁻² or 2.9% ± 1.1%). The combined BzK/ERO/DRG sample contribute 7.1% ±  2.5% to the EBL at 870 μm after accounting for overlap between the populations.

According to Takagi et al. (2007), sBzK galaxies down to K_vega ≲ 20 contribute 3.8 ± 1.2 Jy deg⁻² (or 8.3% ± 3.9%) to the background light at 850 μm (46 ± 16 Jy deg⁻²—Puget et al. 1996; Fixsen et al. 1998). Webb et al. (2004) found that the ERO population down to K_vega < 20.7 constitutes 7%–11% of the EBL at 850 μm, while Takagi et al. (2007) reported a 850 μm EBL contribution from K_vega ≲ 20 EROs of 5.1 ± 1.5 Jy deg⁻² (11% ± 5%) to the total background at 850 μm. Knudsen et al. (2005) found that DRG galaxies down to K_vega < 22.5 contribute ∼7.7 Jy deg⁻² (∼17%) of the EBL at 850 μm.

Next, let us look at the contribution to the EBL at 870 μm from the different samples as a function of redshift and 24 μm flux bins (right-hand panels in Figures 8 and 9). We calculate the contributions using the stacked submm flux and the surface density of sources in each redshift and 24 μm bin.

Looking at the redshift dependence first, we see that the single strongest contribution to the 870 μm EBL by K_vega ⩽ 20 galaxies in a given redshift bin is coming from sources at z ∼ 1.4. Sources at z < 1.4 contribute 3.3 ± 0.2 Jy deg⁻² to the submm EBL, which corresponds to ∼45% of the total contribution from K_vega ⩽ 20 galaxies (Table 2), and ∼11% of the total 870 μm EBL. K_vega ⩽ 20 sources at z > 1.4 contribute only 2.18 ± 0.22 Jy deg⁻² to the total EBL (or 30% of the total contribution from K_vega ⩽ 20 galaxies). Although the K_vega ⩽ 20 galaxies at z > 2 have significant submm emission (Section 5.2), their low abundance means that they contribute ≲1% to the observed submm EBL. In contrast, z < 1 K_vega ⩽ 20 galaxies are so abundant (Figure 4) that despite their low average submm fluxes they contribute significantly to the submm EBL.

The bulk (∼80%) of the contribution to the submm EBL from the 24 μm detected (i.e., S_24 μJy > 27 μJy) sBzK, ERO, and DRG galaxies comes from sources with S_24 μJy ≃ 50–350 μJy, with only a minor fraction coming from brighter, presumably, AGN-dominated sources.

The sBzK and pBzK selections are designed to locate star-forming and passive galaxies in the redshift range 1.4 ⩽ z ⩽ 2.5, while the ERO and DRG color criteria are designed to select extremely red (either due to dust extinction or old age) galaxies at z > 1. Still, the exact definitions of these color criteria are somewhat arbitrary. Given the large number of sources available to us, we are in a position to construct stacks of the 870 μm signal from statistically significant subsets of galaxies across the BzK and RJK diagrams, thus allowing us to see where in these color–color diagrams the submm signal is coming from. Ultimately, this may allow us to fine-tune the sBzK/pBzK criteria, as well as potentially identify regions of the RJK diagram containing dusty/star-forming versus old/passive EROs and DRGs.

In Figures 10 and 11, we show the signal and S/N contours of the stacked 870 μm signal obtained across the BzK and RJK diagrams, where the stacking has been carried out within equally sized grid cells. We have done this for the full sample, as well as for subsets of the sample within the redshift intervals 0 < z < 0.8, 0.8 < z < 1.5, and 1.5 < z < 3.0. These intervals were motivated by the redshift distributions in Figure 4, which showed evidence of distinct populations between 0.8 < z < 1.5 and 1.5 < z < 3.0.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering Figure 10(a) first, we see that the sBzK criterion by Daddi et al. (2004) seems to be robust in the sense that the strongest 870 μm signal is emerging from the sBzK region. Comparing with Figure 1(a) indicates that virtually the entire stacked submm signal from the sBzK galaxies is coming from the subset of sBzK galaxies (with a slight overlap into the pBzK and non-BzK regions), which have also been classified as EROs. Due to their extremely red colors and relatively strong submm signal, these galaxies are very likely to be among the most dusty, star-forming sources of the submm-faint K_vega ⩽ 20 selected galaxies, yet their brightness ensures that they are detected in the blue (thus qualifying as sBzK galaxies).

The submm signal clearly extends into the pBzK region suggesting that some contamination by star-forming galaxies occurs. This is partly due the low density of sources in the pBzK-region which implies that more sources will scatter into the region (due to photometric errors) than out of it. Based on the 870 μm signal and S/N contours in Figure 10, we propose a refinement of the sBzK and pBzK criteria, namely

$$\begin{equation} (z-K)_{{\rm {AB}}} \ge (B-z)_{{\rm {AB}}} -0.5, \end{equation} \tag{ 7 }$$

$$\begin{equation} (B-z)_{{\rm {AB}}} \le 3.4 \end{equation} \tag{ 8 }$$

for sBzKs, and

$$\begin{equation} (z-K)_{{{\rm AB}}} \ge 2.9, \end{equation} \tag{ 9 }$$

$$\begin{eqnarray} &&(B-z)_{{{\rm AB}}} > 3.4 \end{eqnarray} \tag{ 10 }$$

for pBzKs.

That these refined selection criteria might do a better job at selecting star-forming and passive BzK galaxies is confirmed if we consider how the submm signal across the BzK diagram changes with redshift (Figure 10 (c)–(h)). We find that only marginal 870 μm emission emerges from galaxies in the redshift interval 0 < z < 0.8, while in the 0.8 < z < 1.5 interval, significant submm emission (∼0.6 mJy) starts to appear from galaxies lying within the sBzK region. At 1.5 < z < 3 a strong submm signal (∼1 mJy) and nearly all of the submm emission is coming from the sBzK region. This observed migration of the stacked submm signal toward the sBzK region as a function of redshift illustrates the ability of the sBzK criterion to select star-forming galaxies at z ⩾ 1.5. Stacking the 870 μm flux of the galaxies fulfilling the new pBzK criterion (only 17 sources in total) yields −0.10 ± 0.28 mJy, i.e., consistent with zero.

Turning to the RJK diagram, the bulk of the submm signal is found to come from the overlap region between DRGs and EROs, in particular from ERO/DRG galaxies with (J − K)_AB > 1.9. The latter, together with the tentative evidence of an underdensity of sources at (J − K)_AB = 1.9 (see Figure 1), suggests that (J − K)_AB > 1.9 may be a useful criterion for selecting dusty, star-forming galaxies. This is strengthened by the fact that the submm signal becomes even stronger and further concentrated toward the ERO/DRG overlap region if we consider only the subset of sources at z = 1.5–3 (Figure 11 (g)–(h)). In the redshift range 0.8 < z < 1.5, the submm signal is much weaker (∼0.6 mJy) and comes from ERO/DRG galaxies with (J − K)_AB < 1.9.

It is well known that the K_vega ⩽ 20 ERO population is a heterogeneous population, consisting of roughly a 50-50 mix of dusty, star-forming EROs and evolved, passive EROs (Dey et al. 1999; Cimatti et al. 1999; Mannucci et al. 2002). Pozzetti & Mannucci (2000) argued that a crude separation between the two types of EROs could be made based on the criterion (J − K)_AB = 0.34(R − K)_AB − 0.22 (and (R − K_S)_AB ⩾ 3.35),¹⁸ with dusty, star-forming EROs lying above the relation and old, passive EROs below it. While our stacking analysis lends some merit to the Pozzetti & Mannucci criterion, as the bulk of the submm signal is clearly found above it, significant submm emission is also detected in the passive ERO region (in particular for sources at z > 1.5), suggesting that blindly applying the criterion does not produce clean samples of star-forming and passive EROs (see also Smail et al. 2002b).

Finally, we caution that by adopting equally sized grid cells, some cells will contain a significantly larger number of sources than others, thereby introducing a potential skewing of the measured S/N across the diagrams (which is why we also show the variation of the average submm flux density across the diagrams). As a check on our results, we therefore adopted two alternative methods for binning the sources. First, an adaptive mesh was applied by requiring that no more than 200 sources were allowed within a given cell, and, second, the 50 nearest neighbors of each source were identified and stacked. Reassuringly, the two additional binning methods gave results in good agreement with the regular grid results.