"Super-deblended" Dust Emission in Galaxies. II. Far-IR to (Sub)millimeter Photometry and High-redshift Galaxy Candidates in the Full COSMOS Field

In this section, we describe the creation of a K_s catalog that we wish to use to perform prior-based fitting in the MIPS 24 μm, VLA 1.4, and 3 GHz images. While in L18 we used IRAC as a parent catalog of massive galaxies from which to build prior samples, in the COSMOS field, we prefer to start from a K_s catalog because substantially larger amounts of effort were put in by the community to create value-added catalogs for K_s-selected samples, as opposed to IRAC-selected samples.¹⁷

The COSMOS2015 catalog (Laigle et al. 2016) from the UltraVISTA survey (McCracken et al. 2012) contains 1,182,108 YJHK_s sources in total, and over half of them have well-measured photometric redshifts and stellar masses. We included 528,889 sources with photometric redshift and stellar mass in the COSMOS2015 catalog in our initial prior sample. The spectroscopic redshifts are taken from the new COSMOS master spectroscopic catalog (M. Salvato et al. 2018, in preparation) that is based on a variety of spectroscopic surveys published in the literature.¹⁸ Sources with a reliable spectroscopic redshift (redshift quality >3 and $| {z}_{\mathrm{spec}}-{z}_{\mathrm{phot}}| \lt 0.1\,\times (1+{z}_{\mathrm{phot}})$) are fitted at their fixed redshift in SED fitting. Sources in the COSMOS2015 catalog that are flagged to be X-ray detected by Chandra are also included in our prior catalog, but their photometric redshifts are not used (hence all redshift ranges explored) as they might be problematic. While objects located in the regions surrounding saturated optical stars are removed from this Laigle et al. catalog, we fill up these blank regions by adding 57,862 K_s sources from the UltraVISTA catalog of Muzzin et al. (2013), ensuring good coverage and evenness of prior distribution in the full UltraVISTA area. By matching the 3 GHz catalog of Smolčić et al. (2017), we find 2962 radio sources that are not present in the K_s catalog. These radio sources lacking a near-IR counterpart (separation >1'' from any K_s source) are also included in our initial prior catalog (mainly because some might be genuine high-redshift galaxies).

In total, we obtained 589,713 priors in the K_s(+radio) prior catalog. As can be seen from the red squares marked in Figure 1, this prior catalog has a source density of $\langle {\rho }_{\mathrm{beam}}\rangle =1$ source per beam in the MIPS 24 μm image, $\langle {\rho }_{\mathrm{beam}}\rangle =0.2$ source per beam in the VLA 1.4 GHz image, and $\langle {\rho }_{\mathrm{beam}}\rangle \,=0.06$ source per beam in the VLA 3 GHz image, which is appropriate for prior fitting in all these bands. Comparing to the initial IRAC catalog used in L18, the surface density of the K_s+radio priors is compatible to the one in GOODS-N, with both showing 1 source per beam at MIPS 24 μm (i.e., 35 galaxies arcmin⁻²).

Zoom In Zoom Out Reset image size

Note that the UltraVISTA imaging in COSMOS is not homogeneously deep: more sources are detected in the ultra-deep strips. However, there is only a 2.5% difference in prior density in our K_s(+radio) catalog. This difference gets to 0.9% in the following 24 μm+radio+mass-selected prior catalog (see Section 3.5). This is negligible and does not significantly impact our results.

The 589,713 K_s(+radio) priors are used to fit MIPS 24 μm, VLA 3, and 1.4 GHz images. Given that the K_s(+radio) prior catalog over COSMOS contains 56× more sources than the 19,437 IRAC priors in GOODS-N from L18, performing parallelized computations over dozens of CPUs was essential to bring the efficiency of the prior fitting to a manageable level. The detailed image fitting in MIPS 24 μm and VLA images is described in Sections 3.2 and 3.3, respectively.

We obtain PSF-fitting photometry with galfit (Peng et al. 2002, 2010), assuming that the intrinsic size of distant sources is negligible with respect to the size of the PSF, which is a good hypothesis for 24 μm except, perhaps, with rare cases of very low-redshift galaxies that are not the main focus of this paper.

We perform PSF fitting at the positions of 589,713 K_s+radio priors in Spitzer/MIPS 24 μm GO3 images from the S-COSMOS team (PI: D. Sanders; Le Floc'h et al. 2009). The PSF used for the MIPS 24 μm image is identical to the PSF profile used for the 24 μm image fitting of L18 in GOODS-N, which is consistent with the COSMOS one from Le Floc'h et al. (2009) but at much higher S/N. Based on the first-pass results with fixed R.A.–decl. positions, we run a second-pass PSF fitting allowing for up to 1 pixel (12) variations of prior source positions for those high-S/N sources (galfit S/N > 10). This improves the fitting for bright sources while returning a residual image that is cleaner than the one from the first-pass fitting.

After the galfit PSF fitting, we run Monte Carlo simulations in the MIPS 24 μm map, then correct the galfit outputs via the three-step correction recipes based on the simulations (σ_galfit, S_residual, and crowdedness); see Section 4.4 for details. The figures describing simulation performances are shown in Figure 29 in the Appendix C. The flux bias has also been calibrated, and we obtain well-behaved, quasi-Gaussian distributions of flux uncertainties in a similar way to that of L18.

In Figure 2, we compare our 24 μm photometry to the catalogs from Le Floc'h et al. (2009) and Muzzin et al. (2013). Our flux measurements are in excellent agreement with the ones from Le Floc'h et al. (2009), while our flux errors are significantly smaller (see right panel in Figure 2), suggesting that our fitting is indeed more accurate and thus reaches much deeper than this blindly extracted catalog. Note that some galaxies have significantly lower 24 μm fluxes in our catalog as compared to that of Le Floc'h et al. (2009). We believe that the lower flux measurements are due to the resolution of blending of the 24 μm image from our method: fluxes for sources with close neighbors (e.g., with distances less than the beam size, 57) are difficult to measure reliably with the blind extraction method.

Zoom In Zoom Out Reset image size

The comparison with the 24 μm photometry reported in Muzzin et al. (2013) that also adopts fitting at prior positions shows that our measured fluxes are systematically higher than those from Muzzin et al. (2013) by a factor of 1.2, which appears to be a calibration offset that also applies consistently with respect to the Le Floc'h et al. (2009) catalog. Given that Muzzin et al. (2013) did not compare to the catalog of Le Floc'h et al. (2009), the source of the observed difference is unclear. However, we notice that the scaling goes in the opposite direction of what is inferred elsewhere in our work, namely that the 24 μm photometry might need to be recalibrated higher by some factor even with respect to our measurements (and those from Le Floc'h et al. 2009). More details about the calibration of the 24 μm photometry are discussed in Appendix A. Once accounting for this calibration difference, our results suggest that the photometric uncertainties of 24 μm fluxes from Muzzin et al. (2013) are often mildly underestimated.

Radio continuum emission is also an important tracer of star formation. Yun et al. (2001) found a strong correlation between radio and FIR, expressed in terms of the logarithmic ratio q_IR between IR and radio fluxes. Recently, some evolution of q_IR has also been revealed: q_IR appears to be decreasing with increasing redshift (Magnelli et al. 2015; Delhaize et al. 2017). This helps for the detection of high-redshift star-forming galaxies in the radio, because their radio emission is expected to be brighter.

We use the 1.4 GHz Deep Project map (Schinnerer et al. 2010) and the 3 GHz Large Project image (Smolčić et al. 2017) from the VLA-COSMOS team. The 1.4 GHz Deep Project map was combined with the existing data from the VLA-COSMOS 1.4 GHz Large Project map (Schinnerer et al. 2010). It covers 1.7 deg² with an angular resolution of 25, reaching σ ≈ 12 μJy in the central 50' × 50' but shallower elsewhere. The VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017) is based on 384h of VLA observations. The final mosaic reaches a median rms of 2.3 μJy beam⁻¹ over 2 deg² at an angular resolution of 075, corresponding to the best sensitivity and highest resolution of available radio surveys in COSMOS. Smolčić et al. (2017) presented a catalog of 10,830 blindly extracted S/N > 5 radio sources. With our prior-based PSF-fitting technique, we can push to deeper radio flux levels with high fidelity and completeness and a very low spurious detection rate expected.

We use a circular Gaussian¹⁹ PSF with FWHM = 25 and 075, respectively, and run PSF fitting at the positions of 589,713 K_s+radio priors in the 1.4 and 3 GHz images, also keeping account of their rms maps. To improve the fitting, a second-pass fitting is performed in both images, allowing for variation of positions for bright sources (galfit S/N > 20) of up to 2 pixels (04 at 3 GHz, 07 at 1.4 GHz).

Given the high resolution of the radio images, we only ran two-step correction recipes in the Monte Carlo simulations. Correction for crowding (see L18) is not required, as priors in the image are basically never crowded (crowdedness ∼1 always). The simulation recipes return a typical effective rms sensitivity of 2.5–2.7 μJy at 3 GHz with well-behaved Gaussian-like uncertainties, which is close to the nominal 2.3 μJy rms noise, allowing for reliable S/N > 3 detections at 3 GHz down to ∼8 μJy. The 1.4 GHz fitting gives a median uncertainty of 10.22 μJy (albeit better in the central area), shallower than the 3 GHz photometry but useful for SED fitting.

We compare our radio photometry to the blind-extracted catalogs of Schinnerer et al. (2010) and Smolčić et al. (2017), as shown in Figure 3. At 3 GHz, our flux measurements are consistent with the unresolved fluxes in the catalog of Smolčić et al. (2017), while we underestimate the fluxes of resolved sources. In the right panel of Figure 3, we show the comparison with the catalog of Schinnerer et al. (2010): our 1.4 GHz measurements are consistent with the peak fluxes in their catalog. As a result of these measurements, we obtain an additional 7005 S/N > 3 radio sources that have a K_s counterpart but were not reported in the catalog of Smolčić et al. (2017). We consider most of these to be reliable radio detections given that our photometric uncertainties are quasi-Gaussian, and we would therefore expect only of order of 10% of these extra sources to be the result of noise fluctuations (for purely Gaussian noise). Including sources already detected in the Smolčić et al. (2017) catalog, for which we adopt their fluxes, we have a total of 15,645 reliable radio detections in the UltraVISTA area.

Zoom In Zoom Out Reset image size

To improve PSF fitting for resolved sources in the 3 GHz image, and attempting to recover the radio flux contributions that we might be resolving out, we perform prior-based fitting in Gaussian-convolved 3 GHz images whose PSFs were increased to 15 and 2''. Obviously, the convolution increases the noise in the data, raising the median flux uncertainties to σ = 5.57 and 7.66 μJy in the fitting of convolved images with 15 and 2'' PSFs, respectively.

On the other hand, this procedure might also introduce biases, e.g., for sources with radio lobes on the sizes of these beams, which could eventually contribute to the recovered radio flux. Weighting the advantages and disadvantages of using fluxes from convolved radio images, we decided to use the photometry from the original image fitting for the rest of this work because it is much deeper. However, we publicly release the photometry from the convolved images as well, as it might be useful to users from the community.

Early results from our radio photometry catalog constructed in this way were used for the work reported by Daddi et al. (2017), which found that a strong overdensity of radio sources is hosted by the z = 2.5 X-ray-detected cluster Cl J1001, demonstrating interesting prospects for future deep and wide-area radio surveys to discover large samples of the first generation of forming galaxy clusters.

In Figure 4, we plot the predicted flux at 24 μm and 3 GHz as a function of redshift for the faintest SPIRE sources that we could hope to be detected with our survey, using dust continuum SED templates from Magdis et al. (2012) and adopting the evolving FIR–radio correlation from Magnelli et al. (2015). All SEDs are normalized to a common S_{350 μm} = 8.0 mJy, which is the posterior 3σ_{350 μm} detection limit in our final catalog. The flux at 24 μm decreases with increasing redshift, getting below the detection threshold at z > 3–3.5 (lower if the actual 24 μm calibration has to be altered). Meanwhile, as can be seen from the blue solid line in Figure 4, the radio emission at 3 GHz is always above the 3σ_{3 GHz} detection limit at redshift 0–7, suggesting that galaxies within reach of detection at 350 μm are always, in principle, also detectable at 3 GHz. Therefore, including detections at VLA 3 GHz, we can obtain a more complete prior catalog, particularly for z ≳ 3 galaxies, and improve the completeness of our prior sample for fitting PACS, SPIRE, and (sub)mm photometry. However, in reality, given the scatter in the FIR–radio correlation, the effect of noise, and possible flux losses from the use of 075 data, our effective depth in the radio is so close to the actual required limit that we do expect substantial amounts of incompleteness in our priors at z > 3, even when considering radio.

Zoom In Zoom Out Reset image size

We have performed PSF fitting in the negative 3 GHz image to test our procedure. We found 795 S/N > 3 (negative) detections out of the total of 589,713 positions fitted. This is 0.13%, very consistent with the expected one-tail Gaussian probability at 3σ. Comparing to the 15,645 (positive) S/N > 3 detections, the spurious detection rate in the 3 GHz catalog is estimated to be at the level of 5.1%: pretty low. This will not spuriously alter our prior source density while ensuring that we are fitting the IR maps at truly detected positions in most cases. The same test in the 24um negative image shows an S/N > 3 spurious detection rate of 0.07% with respect to fitted positions and 0.5% with respect to actual detections, a return rate of spurious even slightly lower than expected from Gaussian statistics.

From the prior-based fitting at positions of 589,713 K_s+radio priors in 24 μm and radio images, we obtained 88,008 galaxies with 24 μm and/or radio S/N > 3. They are shown as blue dots in Figure 5. As mentioned before, the full sample of 589,713 K_s+radio sources in COSMOS is too dense to be useful in the fitting of FIR/(sub)mm images, where the source density reaches 5–50 sources per beam at FIR/(sub)mm bands. The 88,008 detections at 24 μm and/or radio bands are preferentially included in the prior catalog for FIR/(sub)mm deblending and have much more manageable sky densities, but they are lacking in completeness, as discussed in the previous sections. Thanks to well-measured stellar masses and photometric redshifts in the UltraVISTA catalogs, we can select a supplemented prior sample by stellar mass, exploiting once again the tight connection between stellar mass and SFR and thus IR luminosity.

Zoom In Zoom Out Reset image size

In Figure 5, S/N_FIR+mm > 5 sources from the GOODS-N "super-deblended" catalog (L18) are shown as red dots. We scaled down their stellar masses, multiplying their values by a factor of $\tfrac{5}{{\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}}$ to renormalize the location of all sources to the 5σ detection limit of the catalog. There is a close connection between FIR luminosity and stellar masses in star-forming galaxies (as widely discussed in the literature in terms of the so-called star-forming main sequence; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Pannella et al. 2009; Karim et al. 2011; Rodighiero et al. 2011; Schreiber et al. 2015) that can be seen by the relative small dispersion of ∼0.3 dex of the red dots along their average redshift trend.

The scaled GOODS-N sources statistically locate the positions expected for galaxies detectable, on average, at S/N_FIR = 5 in this mass–redshift diagram. Scaling their trend a factor of 5 lower, we obtained the stellar mass limit shown as a dashed line in Figure 5, which corresponds to the S/N_FIR = 1 limit or to deviations of ∼2.5σ from the average trend. Therefore, a prior catalog including galaxies down to this stellar mass limit will be able to include the vast majority of FIR/(sub)mm-detectable sources. The number of IR-detected outliers in GOODS-N that would drop below our mass limit is only 4% at all redshifts and 7% at z > 2. Notice, though, that we still have 24 μm– and radio-selected priors below this stellar mass limit.

We thus selected 106,420 K_s sources with stellar mass log M_* > 1.8 × log (1 + 4 × z) + 8 as a supplement for the 24 μm+radio priors. In total, we have 194,428 sources in this 24 μm+radio+mass-selected prior catalog. The green dashed line in Figure 5 visually shows this stellar mass selection threshold. Some 96% of all UltraVISTA sources with stellar mass ${M}_{* }\gt {10}^{9.5}\,{M}_{\odot }$ are included in this prior catalog, and 68% of sources with ${M}_{* }\gt {10}^{9.0}\,{M}_{\odot }$ in the z = 0–4 redshift range. This prior catalog has a source density $\langle {\rho }_{\mathrm{beam}}\rangle \sim 0.5$ at PACS 100 μm (i.e., surface density ∼11 sources arcmin^–2), 2× the surface density of the 24 μm+radio prior catalog in L18 (i.e., ∼5 sources arcmin^–2 in GOODS-N).

The remaining 395,285 K_s sources are not included in our prior catalog. Similar to L18, we assume that their flux contributions to the PACS, SPIRE, and (sub)mm images are negligible and do not consider them for the rest of this work. As discussed in L18, their presence will act as a background whose average level will be accounted for consistently by our procedure, while their possibly inhomogeneous distribution will also be accounted for by error bars in the finalized photometry.

We notice that, in principle, we might have supplemented stellar mass–selected priors only for z > 2, as in general we expect our 24 μm catalog to be highly complete at z < 2, at least for SPIRE bands (see Figure 4). However, in this way we might be missing silicate-dropout sources at 1 < z < 2 (Magdis et al. 2011) that are otherwise detectable in FIR, e.g., in PACS images. On the other hand, for the SPIRE bands, our "super-deblended" procedure effectively removes from the fitting pool most of the stellar mass–selected extra priors at z < 2, so their presence is not negatively impacting our results, while guaranteeing a higher completeness for priors at all redshifts and for the PACS bands.

We list below all relevant differences with respect to the L18 GOODS-N work for the process of FIR/(sub)mm deblending (details of some of the most crucial differences are further provided in the following sections).

• (1)  
  We do not subtract faint priors from the PACS 100 μm map, given that the 24 μm+radio+mass-selected prior catalog is less confused in this image ($\langle {\rho }_{\mathrm{beam}}\rangle \sim 0.5$ sources per beam). This step was also almost irrelevant in L18.
• (2)  
  We deblend the SCUBA2 850 μm image before fitting the SPIRE 350 μm map, given that the SCUBA2 850 μm image has a much higher spatial resolution (FWHM = 11'' in the non-match-filtered image) and is very deep in about one-fourth of the COSMOS field. This change of order provides better constraints on the SEDs and helps to improve flux predictions at 350 and 500 μm (particularly useful given that PACS imaging is generally shallower in COSMOS).
• (3)  
  The SCUBA2 850 μm deblending is performed in the non-match-filtered image after removing hot pixels that appear to quite substantially plague the released COSMOS maps (both match-filtered and non-match-filtered). These hot pixels were identified as strong >5σ outliers by median filtering the SCUBA images normalized by their rms maps on scales of twice the beam (20 pixels × 20 pixels) and replacing them with the actual median from 10 × 10 (PSF) filtering. The SCUBA2 PSF profile is obtained by stacking bright sources deemed to be reasonably isolated based on our modeling and fitting the result with a 2D Gaussian. We note that we display the match-filtered version in the multiwavelength cutouts (Figures 7, 22, and 23 and Appendix D) for illustrative purposes.
• (4)  
  Given that the SCUBA2 and AzTEC images display significant depth variations over the COSMOS field, we adapted the flux threshold for accepting priors to the actual local depth, removing from the fitting pool all priors with predicted flux (plus twice the flux uncertainty, as in L18) below the 1σ depth (i.e., we exclude sources with S_SED + 2σ_SED < σ_rmsnoise from the fitting, where σ_rmsnoise varies locally).
• (5)  
  Because the 24 μm, radio, and, most importantly, PACS imaging is shallower in COSMOS with respect to GOODS-N, the flux prediction process at each wavelength step displays considerably larger uncertainties. This implies that for a large fraction of priors, while we can confidently infer that they are faint (below the threshold) and should not be considered for fitting, we cannot accurately determine their actual expected fluxes. This poses a problem to appropriately subtract their best/predicted flux from the images, as they are ill-defined, and doing such might actually degrade the performances of our "super-deblended" process rather than improving them. We deal with this problem in Section 4.3.
• (6)  
  Meanwhile, in order to further cope with the limitation discussed above, we added a fourth step in our calibration recipes, analyzing results as a function of the total flux locally subtracted at each position, thus testing for flux biases and noise variations induced by this step (Section 4.4).

Finally, we notice that the COSMOS field contains a subarea observed by the ESA CANDELS-Herschel key program (PI: M. Dickinson) with much deeper PACS imaging data (Schreiber et al. 2015). We have fitted the deep PACS 100 and 160 μm images in the CANDELS field and used the results based on these deeper images for the CANDELS area. Their fluxes and flux uncertainties have been calibrated by dedicated Monte Carlo simulations executed in the deep images.

Note that PACS fluxes measured from images need to be scaled up by a factor of 1.12× because of the flux losses from the high-pass-filtering processing of PACS images (e.g., Popesso et al. 2012; Magnelli et al. 2013). We have applied this factor in the SED fitting and the released catalog, while this factor is not applied for the flux comparison in Figure 15.

The SED fitting recipes and parameters are identical to those in L18. We recall them in some detail here for clarity. Four distinct SED components are used in the fitting procedure: (1) a stellar component (Bruzual & Charlot 2003) with a Small Magellanic Cloud attenuation law (we recall that we fit only down to the K_s band); (2) a mid-infrared active galactic nuclei (AGN) torus component (Mullaney et al. 2011); (3) dust continuum emission from the Magdis et al. (2012) library with the more updated L_IR/M_dust-redshift evolution taken from Béthermin et al. (2015) to fit galaxy SEDs and predict photometric redshift and FIR/mm fluxes; and (4) a power-law radio continuum with an evolving ${q}_{\mathrm{IR}}=2.35\times {(1+z)}^{-0.12}\,+\mathrm{log}(1.91)$ (Magnelli et al. 2015; Delhaize et al. 2017). We perform SED fitting at fixed redshift for sources with reliable spectroscopic redshifts, while we do allow redshift variations within ±10% × (1 + z_phot) for sources with an optical/near-IR photometric redshift. Using the newly fitted SFR_IRs from our catalog (as updated at each wavelength step) and the stellar masses from Laigle et al. (2016), we perform main sequence (MS)/starburst (SB) classification by measuring the distance to the main sequence from Sargent et al. (2014) at the fitted redshift. Sources are considered to be pure SBs if log(SFR/SFR_MS) > 0.6 dex and SFR/σ_SFR > 3 and to be pure MS if log(SFR/SFR_MS) < 0.4 dex and SFR/σ_SFR > 3 and fitted with the appropriate templates. When a clear MS/SB classification cannot be obtained, they are fitted with all SB+MS templates. In the case of radio-excess sources classified as radio-loud AGNs, we do not include the radio photometry in the fit. Our fairly conservative criterion requires observed radio fluxes 2× higher than the prediction from the FIR–radio correlation with >3σ significance. For sources with a combined S/N ≥ 5 over the 100 μm–1.1 mm range, we do not fit the 24 μm and radio photometry so as to avoid being affected by the scatter of the FIR–radio correlation and the variation of mid-IR features. As already briefly mentioned, the SED fitting is performed before the photometric measurement work at each band from 160 to 1200 μm, predicting the FIR flux at each band in exam before actual measurements. The SED fitting is also performed a last time on the final catalog once measurements in all bands have been made, to provide the final physical quantities released with the catalog.

We show the adopted limits for retaining sources for fitting at each band in Figure 6: priors with S_SED + 2σ_SED > S_cut are maintained and fitted (hereafter selected sources), while sources fainter than this limit are excluded from the fitting (hereafter excluded sources). Here S_SED is the predicted flux based on SED fitting, and σ_SED is its uncertainty based on the χ² statistics.

L18 subtracted the fluxes S_SED of all excluded sources from the GOODS-N maps. However, this blind approach would be problematic in the photometric work in COSMOS. As mentioned in Section 3, the 24 μm, radio, and especially PACS photometry are shallower than that in GOODS-N. Although we have included the 100 μm photometry in the first run of SED fitting to better constrain the SEDs, the PACS 100 μm photometry (σ ∼ 1.44 mJy) in COSMOS is still a factor of 4.5 shallower than the one in GOODS-N (σ ∼ 0.32 mJy), while at PACS 160 μm, the sensitivity ratio is 5.2 (3.55 versus 0.681 mJy), leaving the SEDs less well-constrained in COSMOS than in the GOODS-N field. Moreover, the 106,420 mass-selected K_s priors do not have any 24 μm or radio detection (by construction, although they do have upper limits in these bands), leaving a large uncertainty on their fluxes at the FIR/(sub)mm bands. The large uncertainty on flux predictions introduces errors and possible systematics on faint-source subtraction, which would imply, in turn, flux biases and increased photometric errors on measured fluxes.

To minimize the uncertainty in faint-source subtraction, we decided to subtract SED predicted fluxes only for sources with S_SED/σ_SED > 2 (hereafter subtracted sources) from the original map. This method is applied on PACS 160 μm, SPIRE, and AzTEC images, where the numbers of subtracted sources are listed in Table 1. Using simulations, we verified that this criterion improves the overall performance and produces lower uncertainties for fitted sources. Also, we tested that a threshold around 2σ is optimal for this step, in terms of reducing the final noise with the procedure discussed here. In the left panel of Figure 7, the orange solid circles in each cutout mark the sources subtracted from the map, while the sources shown by orange dashed circles are neither subtracted nor fitted because their predicted flux is definitely faint but uncertain. Treating these sources differently makes a difference because they are obviously a lot more faint than brighter sources intrinsically, so even if their flux is small, combining large numbers of them can have an impact on the result.

Second, we add one more step in the Monte Carlo simulation analysis to identify and correct the flux bias introduced from the subtraction and to include the errors in the finalized photometry. In the bottom panel of Figure 8, we can inspect the difference between the input and measured fluxes S_in − S_out from the simulations as a function of a parameter S_subtracted that is the sum of the total flux of subtracted galaxies convolved with the beam at the position of each specific source examined in the simulation. For positions where S_subtracted is high, S_in − S_out is also large, and the opposite is also true (the median bias is zero at this fourth processing step by construction). This implies that subtracted fluxes are too large, i.e., overestimated. This is expected in a regime of low-accuracy predictions, given that fluxes are always positively defined.

Zoom In Zoom Out Reset image size

Note that we do not apply the 2σ requirement to subtract sources in the SCUBA2 and MAMBO images, where we subtract all the excluded sources, because we found from simulations that this did not make any difference in these bands. This might be owing to the smaller beam size (FWHM = 11'' for SCUBA2) with respect to the SPIRE images, as well as the fact that most excluded priors are at low redshift and have very low fluxes predicted at ≈1 mm.

In Figure 9, we present a portion of the SPIRE 350 μm image as an example of the overall procedure. Panel (2) in this figure shows the image combining all faint sources to be subtracted. It is quite apparent how their crowding in certain regions (clustering) simulates brighter individual objects that are not really individual galaxies but just a spurious superposition of many faint galaxies. In Appendix B, we show the same set of images produced during these steps at each wavelength band.

Zoom In Zoom Out Reset image size

Different from the case in the GOODS-N field, it is impractical and unnecessary to randomly distribute simulated galaxies over the whole COSMOS field. We thus chose an experimental and representative area of 10' × 10' (similar to the size of the GOODS-N field) in the center of the COSMOS field, where the depth of the imaging data was typical/average for the whole map. Monte Carlo simulations are performed in the experimental area of the original images at the 24 μm, 1.4 GHz, 3 GHz, and 100 μm bands (as no sources are subtracted at these bands) and in the same area but in the faint-source-subtracted images at other bands. We simulate ∼5000 galaxies per band, one at a time.

Following the method defined by L18, Monte Carlo simulations are performed, simulating one source at a time in the actual image, with the following steps.

First, the position of each simulated source is randomly generated within the experimental area. Its flux, S_in, is drawn from a uniform distribution in log space within a range of ∼3σ to ∼12σ, where σ is the median flux density uncertainty at each band. We model the source as a PSF and add it in the actual image (which includes all the other, real sources). Second, we include the coordinates of the simulated source in the fitted prior list (i.e., selected sources, as well as additional sources from residual images) and perform our photometric measurements, simultaneously fitting all priors together with the extra simulated source. From the galfit output, the flux S_out and flux uncertainty σ_galfit of the simulated source will be measured. Finally, we repeat this process ∼5000 times. In this way, we obtain ∼5000 values of S_in, S_out, and σ_galfit, realistically representing the measurement process in the real images. Several additional properties are measured for each of the simulated sources: (1) the instrument rms noise value σ_rmsnoise at the position of the simulated source; (2) the local flux in the absolute-valued residual image (hereafter S_residual), measured by considering the absolute values of the pixels of the residual image, within the PSF aperture; and (3) the crowdedness parameter, defined by summing up the Gaussian-weighted distances of all sources at the position of source i (and including it),

$$\begin{eqnarray}&&\mathrm{crowdedness}\equiv \displaystyle \sum _{j=1}^{N}{e}^{(-{d}_{j,i}^{2}/{d}_{\mathrm{PSF}}^{2})},\end{eqnarray} \tag{ 1 }$$

where d_j,i is the angular distance in arcsec from source j to source i and d_PSF is the FWHM in arcsec of the PSF. In this way, the crowdedness is a weighted measure of the number of sources present within the beam, including the specific source under consideration. These parameters, measurable in the same way for all fitted priors, provide key information to check the expected quality of fitting and the actual local crowding (hence blending) of prior sources. They will be used to calibrate the flux bias corrections and flux uncertainties of each source.

We have analyzed the simulations following the method detailed above, with some important changes with respect to L18. First, we use the newly defined S_subtracted variable (i.e., the sum of the subtracted fluxes from excluded sources within the PSF circle) in a fourth-step calibration in the analysis of simulations (L18 had only three steps, similar to the first three shown in Figure 8 for flux biases). We introduced this new parameter because of the lower overall depth of COSMOS in many bands, so we expect that the actual noise will be appreciably higher in regions where more flux was subtracted from faint sources (see Figure 10).

Zoom In Zoom Out Reset image size

Using the simulations, we derive the dependence of the flux bias S_in − S_out on the four adopted observational parameters. As shown in Figure 8, we fit the flux bias in each bin using a third-order polynomial function. Meanwhile, we determine a correction factor (i.e., the "σ corr. factor" in Figure 10) on the flux uncertainty at each step, defined as the multiplicative factor that is required to scale the rms dispersion of (S_in − S_out)/σ in each bin to 1. These polynomial functions and correction factors are then applied to the four measured observational parameters of each real source. As shown in the right column of Figure 10, the four-step correction significantly improves the performance of the deblended photometry by lowering the overall rms noise of the photometry and reducing (often removing) the prominent non-Gaussian tails in the flux uncertainty distributions.

Also, we have to account in COSMOS for more substantial effects of depth variations of the data across the field. This is often important toward the edges of each data set but also, in some cases, across the data sets, as for SCUBA2. We find that we can rescale simulation results for calibration of flux uncertainties to areas with different depths in the same data set by calibrating performances based on normalized observables (flux bias terms are not affected). Notably, we have normalized the parameters used in the simulations by the local rms noise; i.e., we consider quantities as σ_galfit/σ_rmsnoise, S_residual/σ_rmsnoise, and S_subtracted/σ_rmsnoise (this is not necessary for crowdedness, which does not depend on the noise in the data set).

We tested this scaling by directly performing additional simulation in data portions with different depths and comparing the direct simulation results from those rescaled from a region with different depth. For example, for the SCUBA2 image, we executed our primary simulations in the deep area; we also carried out our independent simulations in a 3.6× shallower area, then applied the scaled deep area simulation results on the galfit outputs to the shallow area. We find that the flux uncertainties scaled from the deep area simulation are slightly larger than those directly performed on the shallower data, but they are overall consistent (or at least conservative). As a further test, we performed PSF fitting in the inverted SCUBA2 image and found 97 (negative) S/N > 3 detections, after calibration of the fluxes and uncertainties as for positive sources. The spurious fraction at S/N > 3 is 0.4% with respect to the number of fitted priors. This is close to Gaussian expectations, albeit a factor of 2 higher. Compared to positive S/N > 3 detections, the SCUBA2 spurious detection rate is 9.5%.

Figures for all band simulations are reported in Appendix C.

Although the 24 μm+radio+mass-selected prior catalog is expected to have high completeness, some FIR emitters are still likely missing in this prior catalog. For example, in Figure 4, there are several GOODS-N S/N_FIR+mm > 5 galaxies below our stellar mass limit (red dots below the green dashed line), suggesting that galaxies with lower stellar masses could be detectable in the FIR imaging, e.g., if they are starbursts. On the other hand, K_s-undetected galaxies (K_s dropouts) that are not already included in the Smolčić et al. (2017) 5σ radio catalog are also missing from our prior selection. These missing priors, which are probably very high-z galaxies or extremely low-mass starbursts, might have detectable fluxes at FIR/(sub)mm bands and emerge in the residual images of our photometric products (see Figure 11). Extracting these sources will improve the photometry of the sources around them, in addition to providing potentially interesting high-redshift candidates.

Zoom In Zoom Out Reset image size

Similar to L18, we blindly extract residual sources with SExtractor (Bertin & Arnouts 1996) in the residual maps at each wavelength. The positions of the residual sources are then added to the selected prior list and fitted together in the second pass. The galfit outputs of the second-pass fitting are then corrected via simulation recipes and archived in the finalized photometry catalog. We show an example of extraction of SCUBA2 residual sources in Figure 11. We convolve the residual image (panel 1) with a Gaussian filter to improve the visibility and detectability (panel 2) and show the extracted sources with red circles in panel 3. The final residual image (panel 4) is cleaner after prior+residual source fitting. Residual sources with S/N_{850 μm} > 2.5 are kept in the prior list and fitted in subsequent fitting of SPIRE, AzTEC, and MAMBO images. The counts of residual sources at each band are listed as N_add in Table 1.

Note that we have included residual sources in the Monte Carlo simulations, although there is no obvious difference from masking residual sources in simulations. We do not extract any residual sources in PACS 100 and 160 μm images because of their shallow depths and clean residual maps. Also, the residual images of AzTEC and MAMBO are clean; no residual sources are extracted there.

At PACS 100 and 160 μm, we matched our catalog to the PEP catalog (Lutz et al. 2011) with a tolerance of 1''. In Figure 15, our flux measurements are generally consistent with the measurements in the PEP catalog but showing a tail of sources for which we suggest systematically lower fluxes than PEP. The PEP catalog is produced by fitting 47,437 priors selected at 24 μm from Le Floc'h et al. (2009), while 5× more priors are fitted in our work. We thus believe that the lower flux measurements found in our work are due to sources blending in the PEP catalog. There is also a systematic effect affecting fluxes of all sources, as can be seen by plotting the flux differences (Figure 15). There is a constant difference of 0.8 mJy at 100 μm and 1.8 mJy at 160 μm, which are both of the order of 0.2× the rms noise at each band. We attribute these constant offsets to a small background difference applied in the two works, although we advocate that we could measure the background to better than this accuracy with our simulations; hence, we tend to believe that our measurements are correct. Our uncertainties agree well with those in the PEP catalog, as shown by the right panels of Figure 15.

At the SPIRE bands, we compare our results to the XID⁺ catalog (Hurley et al. 2017), which deblends SPIRE images via an MCMC-based prior-extraction method on 52,092 priors selected as 24 μm detections. The XID⁺ catalog covers an area of 2.27 deg²; we limit the comparison to sources with goodArea = 1 and obtain 30,372 matches with a tolerance of 1''. Apart from the matched sources, there are 161,252 priors in our catalog that are not listed in the XID⁺ catalog and 8200 XID⁺ priors missing from our list. The first is due to our deeper 24 μm photometry and the use of radio and especially of mass priors. The latter is likely mis-associations of fluxes at 24 μm.²⁰ These are ∼83% and ∼4% of the 24 μm+radio+mass-selected priors, respectively. In the 161,252 priors missed by XID⁺, we have 2198 detections at 250 μm, 1175 detections at 350 μm, and 1097 detections at 500 μm. Hence, we obtain more detections benefiting from the higher completeness of our prior catalog. Also, the lack of these extra priors in XID⁺ likely exacerbates blending and flux-boosting problems, leading to flux discrepancies with sources in common among the two works, as discussed below.

In Figure 16, we show the flux comparison of matched sources that are detected with S/N > 3 either in our "super-deblended" catalog or in the XID⁺ catalog. We find that our measurements are generally consistent with the measurements in the XID⁺ catalog at high fluxes, while there are significant discrepancies toward faint fluxes. We highlight two populations with significant flux discrepancies, which are shown in blue and green in each panel. The blue dots show sources that have S/N > 3 only from XID⁺ and lower fluxes and S/Ns in our catalog. By inspecting the priors around these sources, we find that they are located in crowded regions. We are generally deblending the signal among different galaxies based on our physical rather than statistical approach. Hence, we believe that XID⁺ often overestimated their fluxes.

Meanwhile, there are some sources that are only detected in this work (green dots in each panel) and not by XID⁺. We have visually checked these sources on the original and the faint-source-subtracted images. We find that most of them have visible signals in both maps. We do not know why XID⁺ catalogs do not retain these sources that appear to be significant according to our work. Quantitatively, in the first panel of Figure 16, there are 12,947 sources in total; 1360 of them have S_{250 μm,us}/S_{250 μm,XID+} > 1.5, while 4479 sources have S_{250 μm,XID+}/S_{250 μm,us} > 1.5. In the second panel, there are 6887 sources in total; 1388 of them have S_{350 μm,us}/S_{350 μm,XID+} > 1.5, while 2,014 sources have S_{350 μm,XID+}/S_{350 μm,us} > 1.5. In the third panel, there are 2498 sources in total; 978 of them have S_{250 μm,us}/S_{250 μm,XID+} > 1.5, while 444 sources have S_{500 μm,XID+}/S_{500 μm,us} > 1.5. So it seems that at the shortest SPIRE wavelength, we are deblending many of the XID⁺ detections into smaller/weaker sources, but this is getting less imbalanced at 350 μm and reversed at 500 μm, probably because of the rapidly decreasing fraction of sources seen by XID⁺ at these longest wavelengths.

We show the flux uncertainties for each band from our and the XID⁺ catalog in the bottom panels of Figure 16. The XID⁺ catalog has lower uncertainties in each SPIRE band. Since the flux uncertainties in our catalog have been carefully calibrated by simulations and display well-behaved Gaussian-like uncertainties, we believe that XID⁺ generally underestimates the flux uncertainties in the SPIRE bands.

In Figure 17, we cross-match common sources (within a limiting separation of 6'' because of the large SCUBA2 PSF) and compare our SCUBA2 850 μm measurements to the deboosted fluxes in Geach et al. (2017). Our "super-deblended" fluxes are consistent with the deboosted ones. Flux uncertainties are also fairly consistent with the deboosted uncertainties in Geach et al. (2017). We have 1020 galaxies (536 priors and 484 additional sources) detected with S/N_{850 μm} > 3 from SCUBA2. This is a 3.3× larger sample than the 306 detections reported by Geach et al. (2017) in the COSMOS field.

We further use 1000+ public ALMA archival data at (sub)mm wavelengths in the COSMOS field to verify our deblended 850 μm photometry. These ALMA data were reduced, imaged, and analyzed in the ongoing ALMA archive mining project A³COSMOS (D. Liu et al. 2018, in preparation).²¹ The A³COSMOS team has processed all public ALMA archives in the COSMOS field and obtained accurate (sub)mm continuum photometry with both blind source extraction (mainly with the code PyBDSM) and prior source fitting. The catalogs are already available and verified with extensive Monte Carlo simulations and will be publicly released after the publication of their first paper (D. Liu et al. 2018, in preparation). Here we use their catalogs in advance (private communication), benefiting from the exquisite ≈1'' ALMA resolution that resolves all blending issues, to verify our super-deblended flux measurements.

Figure 18 shows the flux comparison for common sources between this work and the A³COSMOS prior source fitting catalog, including all prior sources within the ALMA primary beam. Both works used optical+near-infrared+radio priors; therefore, common sources can be easily cross-matched without confusion/multiplicity issue. This figure shows good agreements between the two catalogs without systematic bias. The error bars from this work are in general larger than those of the ALMA photometry, as expected (while ALMA, on the other hand, currently only covers a very small fraction of the COSMOS field).

Zoom In Zoom Out Reset image size

In Figure 19, we further show the histogram of the flux difference between the two works normalized by the total uncertainty (quadratic sum of both errors from this work and A³COSMOS). This comparison is very similar to the histograms shown for simulations for all bands in the Appendix C figures. All cross-matched sources, including low S/N and nondetections, are included in this histogram. The median of the distribution is −0.14, and the σ (scatter) is 0.9, which is very close to 1, indicating that the scatter is consistent with the photometric error. A Gaussian fit with σ fixed to 1 is overlaid on the histogram in Figure 19. If we limit this comparison to the 136 ALMA sources brighter than 4 mJy, we find a small (0.3σ) average flux overestimate by ALMA, and the (uncertainty normalized) scatter rises only to 1.05. There are two >3σ outliers out of 136 galaxies (1.5%). The ALMA preselection is slightly biasing the comparison toward brighter ALMA fluxes for this subsample, by construction. This allows us to conclude that our fluxes and flux uncertainties at 850 μm from the "super-deblending" technique are well defined and correctly derived.

Zoom In Zoom Out Reset image size

There are three main classes of candidate z > 4 galaxies: (1) galaxies in our prior sample and with an existing optically based photometric redshift (as in L18, we eventually rederive a redshift from FIR SED fitting, constrained to be within 10% of the optical (1 + z) to avoid catastrophic failures); (2) galaxies in our prior sample for which no optical photometric redshift exists, mostly because they were selected from radio and lack an obvious match to the K_s-selected catalogs (we only have an FIR redshift here); and (3) additional sources added from the residual maps after the first-pass photometry (here also we only have an FIR redshift, and some of these are possibly spurious sources or unresolved blends; see also L18).

The FIR-based photometric redshifts from SED fitting are derived by considering the χ² of the fit as a function of redshift and applying a Δχ² = 2.5 criterion (corresponding to the 90% probability confidence) to determine the uncertainty of the redshift estimate following Avni (1976). High-z candidates are fitted using the Magdis et al. (2012) SB template for the dusty SF component, which in practice is GN20. This is appropriate regardless of the nature of SB or MS galaxies, given that even MSs have fairly warm SEDs at z > 4 (Béthermin et al. 2015; Schreiber et al. 2017). Still, the intrinsic SED shape as driven by dust temperature might have substantial variations among high-z galaxies. However, our approach is likely conservative, as GN20 at z = 4 has a relatively low dust temperature (T_dust = 35 K) compared to other z > 4 dusty star-forming galaxies (e.g., Smolčić et al. 2015; Riechers et al. 2017). If using some of the latter as templates, the redshifts would have been higher by up to 20%–25% (or about 5% × (1+z), depending on the template). We were also careful of using full AGN (Mullaney et al. 2011) and SF (Magdis et al. 2012) decomposition of the SED only in the presence of sufficient S/N in the photometry for allowing this, e.g., good detections in the mid-IR as well as FIR/(sub)mm ranges. When such a condition was lacking, we considered only the dusty SF component, to avoid cases in which the hot AGN SED would spuriously produce very high-redshift solutions with low ${\chi }^{2}$ by dominating the SED fit. Further details of the dust temperature–redshift degeneracies are discussed in Section 7.2.

Figure 21 (top panel) shows the distribution of redshift errors for the three classes. It is obvious that the redshift uncertainty can grow to be fairly large for objects where only the FIR SED is used to constrain it. As a result, we decided to limit the discussion of high-redshift candidates to those objects satisfying $z-{\rm{\Delta }}z\gt 4$. Consequently, among significantly FIR-detected sources with existing ${z}_{\mathrm{phot},\mathrm{opt}}$, we obtain 31 sources with ${z}_{\mathrm{phot},\mathrm{FIR}}-{\rm{\Delta }}{z}_{\mathrm{phot},\mathrm{FIR}}\gt 4$ for class (1) above. Among 420 (mainly radio) sources lacking an optical photometric redshift, we select 20 sources with ${z}_{\mathrm{phot},\mathrm{FIR}}-{\rm{\Delta }}{z}_{\mathrm{phot},\mathrm{FIR}}\gt 4$ for class (2) above. There are 126 additional sources from the SCUBA2 residual image that reach ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\gt 5$ (see Section 4.5), and 34 of them have ${z}_{\mathrm{phot},\mathrm{FIR}}-{\rm{\Delta }}{z}_{\mathrm{phot},\mathrm{FIR}}\gt 4$ and are in class (3). In total, we have 85 sources with ${z}_{\mathrm{phot},\mathrm{FIR}}-{\rm{\Delta }}{z}_{\mathrm{phot},\mathrm{FIR}}\gt 4$ as our final sample of high-redshift candidates. In the first three panels of Figure 23, we present examples that represent the different classes of candidates: (1) prior with ${z}_{\mathrm{phot},\mathrm{NIR}}$, (2) radio prior without ${z}_{\mathrm{phot},\mathrm{NIR}}$, and (3) additional source with IRAC counterpart.

Zoom In Zoom Out Reset image size

Note that no additional candidates with ${z}_{\mathrm{phot},\mathrm{FIR}}\,-{\rm{\Delta }}{z}_{\mathrm{phot},\mathrm{FIR}}\gt 4$ are selected from the SPIRE residual image, which is typically at a lower redshift. The relative sensitivity of the SCUBA2 850 μm map is substantially higher than that of any SPIRE band for $z\gt 4$ star-forming galaxies. Also, no further residual source candidates at $z\gt 4$ are found from the AzTEC and MAMBO maps. These bands were fitted after the SCUBA2 residual sources were already added. They generally support the reality of the SCUBA2 residual sources.

Using the far-IR colors of a galaxy to determine its redshift (in the absence of z_spec) is known to suffer from the so-called T_dust−redshift degeneracy. In particular, the same colors could be fit by a colder template placed at lower redshift or by a warmer template placed at higher z. To quantify this degeneracy but also determine how the uncertainties in optical z_phot affect the derivation of the dust temperature, or similarly of the mean radiation field, <U> = L_IR/M_dust, we perform the following simulation. We build DL07 models (Draine & Li 2007) with representative γ, q_PAH, U_min, and <U> parameters following Magdis et al. (2012) and calculate flux densities at the MIPS 24 μm, PACS, SPIRE, and SCUBA2 850 μm bands by placing the template at a wide range of fixed redshifts (z_or = 0–6 with a step of 0.01). We then perform SED fitting using the full suite of DL07 models, fixing the redshift first at z_max = z_or + Dz × (1 + z_or) and then at z_min = z_or −Dz ×(1 + z_or), where Dz = 0.03, corresponding to the average photo-z uncertainty in the COSMOS field (Laigle et al. 2016). Thus, at each redshift, we derive three <U> values—${U}_{{z}_{\mathrm{or}}}$, ${U}_{{z}_{\min }}$, and ${U}_{{z}_{\max }}$—and quantify the impact of the photo-z uncertainty in <U> by considering ${DU}=({U}_{{z}_{\mathrm{or}}}-{U}_{{z}_{\min }})/{U}_{{z}_{\mathrm{or}}}$ and $({U}_{{z}_{\mathrm{or}}}-{U}_{{z}_{\max }})/{U}_{{z}_{\mathrm{or}}}$. Our analysis yields an offset between the input and the extracted <U> with the case of z_min (z_max) systematically overestimating (underestimating) the true <U> by 15%. We conclude that the uncertainty in z_phot,FIR introduces an extra small (at least for the case of Δz = 0.03 × (1 + z)) uncertainty of 15% in the determination of <U>. Repeating the process for the case of a modified blackbody (MBB) with a single T_dust and fixed β = 1.8 yields an uncertainty in the derived T_dust of ∼5%.

Furthermore, both MS and SB galaxies at any redshift are expected to exhibit a range of intrinsic <U> values (ΔU). Consequently, ΔU is expected to introduce an uncertainty in the far-IR-based z_phot,FIR (Δz_phot,FIR = z_spec − z_phot,FIR). In order to quantify Δz_phot,FIR, we first need to adopt a reasonable value for ΔU. Since <U> ∝ L_IR/M_dust ∝ L_IR/(M_gas × Z) ∝SFE/Z (e.g., Magdis et al. 2012), where Z is the gas-phase metallicity, we can estimate the intrinsic range of <U> in terms of variations of star formation efficiency (SFE) and Z within and outside the MS. Assuming the various relations reported in the literature between SFE, distance from the main sequence (ΔMS), and Z (e.g., Magdis et al. 2012; Sargent et al. 2014; Tacconi et al. 2018), we estimate an intrinsic variation of ΔU ≈ 0.1–0.2 dex for MS galaxies. This is in agreement with the empirically determined 0.2 dex variation reported by Magdis et al. (2012) for local normal galaxies. Performing simulations as outlined above, we find that ΔU = 0.1–0.2 dex corresponds to Δz_phot,FIR ≈ 0.06–0.1 × (1 + z). We note that a similar estimate is reached for SBs at any redshift, assuming a range in gas depletion timescale of 50–200 Myr. Thus, we conclude that the intrinsic variation of the shape of the far-IR SED of both MS and SB galaxies at any redshift introduces an intrinsic uncertainty of Δz_phot,FIR ≈ 0.06–0.1 × (1 + z) that should be taken into account along with the uncertainties introduced by the photometric errors and the varying photometric coverage.

In some cases, the high-redshift candidates are interestingly coincident with previously existing low-redshift priors in our sample that were discarded for fitting at some longer wavelengths, but where a strong detection was later found in our procedure (often a residual source). This is similar to the case of source ID20003080 (AzTEC/C160) shown in Figure 22. This source is coincident with an early-type galaxy with known spectroscopic redshift $z=0.36$ that is associated with a mass-selected prior (ID659416), which is fitted at PACS bands while excluded at longer wavelengths because it is obviously too faint. Thanks to its radio detection, we already knew in this case of the presence of the very close (12 separation) radio source ID20003080 (i.e., Cosbo-7 in Bertoldi et al. 2007, AzTEC/C160 in Aretxaga et al. 2011) that is afterward solidly detected at SCUBA2 850 μm (S/N_{850 μm} = 8.2) and MAMBO 1.2 mm (S/N_{1.2 mm} = 7.4). Source ID20003080 appears to be a high-redshift galaxy that is aligned with (and perhaps gravitationally lensed by) an early-type galaxy at z = 0.36, likely part of a group or poor cluster of galaxies at the same redshift (see the RBG image in Figure 22, top right panel). We obtain a photometric redshift z = 4.0 ± 0.6 by fitting 24 μm to radio photometry.

Zoom In Zoom Out Reset image size

Other very similar cases, with background sources at possibly even higher redshift, are ID85004261 and ID20003117 (see Table 3 and Appendix D). Investigating their nature in detail is left to future work.

We find 17 sources among our robust sample at $z\gt 4$ that require an AGN component for their SED fitting (see Table 2). The criterion adopted to conclude this requires that the AGN flux normalization divided by its total uncertainty range (${S}_{\mathrm{AGN}}/{\sigma }_{\mathrm{AGN}}$) is larger than 2. This is supported by visual inspection and appears to be reliably returning reasonable results. In this AGN sample, there are eight sources with ${S}_{\mathrm{AGN}}/{\sigma }_{\mathrm{AGN}}\gt 3$ that we consider to be solid candidates, and the remaining nine sources with $2\lt {S}_{\mathrm{AGN}}/{\sigma }_{\mathrm{AGN}}\lt 3$ are considered to be tentative. The bolometric luminosities of the AGNs range over a quite remarkable ${L}_{\mathrm{bol}}^{\mathrm{AGN}}\sim {10}^{46-47}$ erg s⁻¹, fully in the QSO regime. We show their multiband cutouts and SEDs in Appendix D, where the AGN component is marked by the red curve in the SED panel.

Table 2.  High-redshift Candidates that Fitted with a Significant AGN Component

ID	zphot,opt	zspec	S/NFIR+mm	zphot,FIR	LIR,SF	Lbol,AGN
					$\times {10}^{12}\,{L}_{\odot }$	$\times {10}^{12}\,{L}_{\odot }$
223720	4.19	⋯	6.9	4.66 ± 0.28	5.5 ± 0.3	7.8 ± 2.4
339785	4.92	⋯	5.2	5.45 ± 0.45	5.3 ± 0.2	54.6 ± 17.3
347052	4.34	⋯	5.6	4.82 ± 0.14	3.8 ± 0.4	21.3 ± 7.9
551174	4.21	⋯	5.3	4.68 ± 0.16	6.1 ± 0.3	6.9 ± 2.7
556890	4.49	4.63	9.6	4.63	8.9 ± 0.1	16.1 ± 1.6
578482	4.32	⋯	10.7	4.75 ± 0.22	6.1 ± 0.4	8.8 ± 3.3
599184	4.41	⋯	7.6	4.89 ± 0.49	4.7 ± 0.4	8.5 ± 4.1
613818	4.47	⋯	6.8	4.67 ± 0.49	5.9 ± 0.8	12.5 ± 4.3
632541	4.29	⋯	6.5	4.77 ± 0.38	4.4 ± 0.5	12.0 ± 5.4
695002	4.31	⋯	13.7	4.79 ± 0.02	4.1 ± 0.2	48.8 ± 4.9
739920	4.42	⋯	7.7	4.76 ± 0.49	4.7 ± 0.3	13.6 ± 6.7
786213	⋯	4.34	32.6	4.34	15.3 ± 0.1	5.1 ± 0.7
965647	4.34	⋯	6.2	4.58 ± 0.19	4.5 ± 0.4	5.6 ± 2.7
10008348	4.16	4.45	5.2	4.45	5.6 ± 0.1	8.4 ± 1.2
10015010	3.91	⋯	8.4	4.35 ± 0.18	5.5 ± 0.3	11.0 ± 2.9
10137954	⋯	4.64	8.0	4.64	0.2 ± 0.2	34.9 ± 2.4
10213589	4.66	⋯	7.1	4.15 ± 0.03	4.2 ± 0.2	23.6 ± 4.9

Download table as:  ASCIITypeset image

For the 85 z > 4 candidates, the sum of the SFR is $7.1\times {10}^{4}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The bolometric AGN luminosity contained among AGN detections adds up to ${L}_{\mathrm{bol}}^{\mathrm{AGN}}={10}^{48}$ erg s⁻¹, corresponding to an integrated black hole accretion rate of $160\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The ratio of the black hole accretion rate to the SFR is thus 2 × 10⁻³, close but higher than the universal ratio discussed in Mullaney et al. (2012) despite the IR selection that should favor star formation. This might suggest some evolution toward stronger AGN activity in these sources that should be confirmed with further studies. The number density of the 17 objects in COSMOS (assuming volume within 4 < z < 5) is ϕ ∼ 8 × 10⁻⁷ Mpc⁻³ dex⁻¹, which matches reasonably well with the X-ray luminosity function extrapolations from Vito et al. (2018), assuming a bolometric correction of L_bol = 25 × L_{X,2–10 keV}. Our AGN detections are predicted to have X-ray fluxes of (2.9–31.1) × 10⁻¹⁵ erg cm⁻² s⁻¹ at 0.5–2 keV in the case of no obscuration, which would make them all detectable with the limiting depth of available X-ray imaging data (2.2 × 10⁻¹⁶ erg cm⁻² s⁻¹ at 0.5–2 keV) in the Chandra-COSMOS legacy survey (Civano et al. 2016). However, only two sources, ID10008348 and ID10137954, are cross-matched with the point-source catalog of Civano et al. (2016). These X-ray sources show very high obscuration with L_bol/L_{X,2–10 keV} ∼ 300, even higher than the X-ray bolometric correction factors of Compton-thick AGNs in Brightman et al. (2017). The low X-ray detection rate and large L_bol/L_X ratio clearly suggest that our IR-selected AGNs are heavily obscured in the X-ray, consistent with the current understanding of high-redshift populations (Vito et al. 2018).

We list the 34 additional high-z candidates in Table 3. Given that the positions of these additional sources are blindly extracted from the SCUBA2 residual image that has a beam FWHM = 11'', we visually searched for counterparts in SPLASH, VLA 3 GHz, and ALMA 1.2 mm (M. Aravena et al. 2018, in preparation) images with a tolerance of 5'' and set the coordinates of their counterparts as their final positions. We find that 13 additional sources have a well-defined counterpart in the SPLASH and/or VLA and/or ALMA images. The counterparts are presented in cutouts in Appendix D and listed in Table 3. These associations are unlikely to happen by chance and appear robust: by cross-matching to the COSMOS2015 catalog for SPLASH and to our VLA catalog, we would expect only 0.03 chance coincidences per 5'' radius aperture, down to the flux limits reached in these probes. This suggests at most ∼1 spurious association among the 34 additional candidates. The identification of these 13 sources with counterparts and solid detections at multiple FIR/(sub)mm wavelengths suggests that the spurious fraction among the additional candidates is quite contained.

Table 3.  High-redshift Candidates Found among Additional Sources from the SCUBA2 Residual Map

ID	R.A. (J2000)	Decl. (J2000)	Counterparta	Nameb	Distancec	S/NFIR+mm	S/N850 μm	SFRd	ΔSFRd	zphot,FIRe	Δzphot,FIRe
85000019	150.24172	1.60795	SPLASH, 3 GHz	⋯	32	6.8	3.3	1194	309	5.8	1.3
85000436	149.40348	1.73843	3 GHz	⋯	19	5.7	4.2	981	205	5.4	0.9
85000496	149.49761	1.77057	⋯	⋯	⋯	5.3	4.3	702	324	6.0	1.8
85000552	150.04723	1.78372	⋯	⋯	⋯	5.5	3.3	534	104	5.5	1.4
85000762	149.53400	1.85461	⋯	⋯	⋯	6.4	3.8	528	131	4.8	0.8
85000922	150.49660	1.90440	SPLASH	⋯	22	6.2	3.8	1232	343	6.3	1.3
85001050	149.98028	1.93565	⋯	⋯	⋯	5.1	4.1	704	119	5.8	1.6
85001505	149.46728	2.09288	SPLASH	⋯	23	6.4	4.9	816	205	5.9	1.4
85001571	150.02743	2.11300	SPLASH	⋯	04	5.0	2.8	488	24	7.0	2.3
85001674	150.30121	2.14766	SPLASH, 3 GHz	⋯	16	6.1	4.2	836	128	5.1	0.9
85001756	149.87579	2.17826	⋯	⋯	⋯	5.8	4.1	621	151	7.2	2.2
85001929	150.10971	2.25753	SPLASH	⋯	08	12.9	6.8	692	151	5.2	0.5
85001969	149.97064	2.31366	ALMA	AzTEC/C71	30	5.5	4.5	606	103	7.9	2.2
85002134	150.24795	2.38802	⋯	⋯	⋯	5.7	3.3	850	107	10.0	1.1
85002215	150.10063	2.33484	SPLASH, 3 GHz	AzTEC/C114	35	9.6	6.7	721	47	5.9	0.8
85002966	149.42140	2.57688	⋯	⋯	⋯	7.6	5.2	692	112	4.7	0.7
85003151	150.02104	2.63323	⋯	AzTEC/C132	⋯	5.5	3.3	678	154	6.2	1.7
85004261	150.05635	2.57327	ALMA, 3 GHz	AzTEC/C10	26	8.3	5.8	1009	267	7.2	2.2
85005253	150.41618	1.90769	⋯	⋯	⋯	5.2	4.3	1413	223	8.6	1.7
85005277	150.57793	1.93306	⋯	⋯	⋯	5.0	3.3	1425	316	10.0	2.3
85005285	150.68263	1.95067	3 GHz	⋯	49	5.1	3.1	1053	509	7.6	2.2
85005338	149.96455	1.98236	⋯	⋯	⋯	5.2	3.3	898	144	10.0	1.9
85005422	149.98861	2.09413	⋯	⋯	⋯	5.1	3.4	601	57	7.8	2.1
85005464	149.68279	2.09387	⋯	⋯	⋯	5.5	2.7	570	152	5.4	1.4
85005517	150.64576	2.14272	⋯	⋯	⋯	5.0	2.5	934	272	6.7	1.8
85005620	150.70247	2.25888	⋯	⋯	⋯	5.1	2.6	989	312	7.1	2.0
85005670	150.60669	2.31664	⋯	⋯	⋯	5.3	2.9	828	136	5.7	1.4
85005722	150.61282	2.41777	⋯	⋯	⋯	5.1	2.9	1171	267	5.6	1.5
85005759	149.78323	2.40500.	⋯	⋯	⋯	5.1	3.2	504.	93	5.9	1.8
85005769	150.36338	2.41572	⋯	⋯	⋯	5.9	3.2	993	287	5.3	1.1
85005926	149.95379	2.56544	SPLASH, 3 GHz	⋯	20	7.0	2.6	656	114	5.5	1.0
85005933	149.75172	2.58143	⋯	⋯	⋯	5.2	2.7	582	168	7.2	2.5
85005963	150.69766	2.60405	⋯	⋯	⋯	7.5	2.7	1176	574	6.0	1.6
85006141	150.45836	2.81048	3 GHz	⋯	18	8.3	2.5	1199	473	6.0	1.4

Notes. We report here the result for our search of Spitzer IRAC and radio counterparts.

^aCounterpart image. SPLASH: P. Capak; 3 GHz: Smolčić et al. (2017); ALMA: 1.2 mm continuum image (M. Aravena et al. 2018, in preparation). ^bReference name of the counterpart (Aretxaga et al. 2011). ^cDistance of additional source to its counterpart. ^dSFR and ΔSFR: median fit and uncertainty of SFR based on FIR+(sub)mm SED fitting, Chabrier IMF (Chabrier 2003). ^ez_phot,FIR and Δz_phot,FIR: the photometric redshift and uncertainty based on FIR+(sub)mm SED fitting.

Download table as:  ASCIITypeset image

Like for the other candidates, we fit their photometry at 100–1200 μm and 3 GHz using the GN20 template and an evolving q_IR without an AGN component. The 3 GHz photometry is taken from counterparts in the VLA 3 GHz image, while for candidates without any counterpart, we do not fit the 3 GHz photometry. Note that there are three additional sources that result in a best-fit z = 10, the highest value that we allow. These three sources are only detected at (sub)mm wavelengths while their SEDs are weakly constrained by Herschel. Given the redshift errors Δz_phot,FIR = 1.1–2.3, it is quite plausible that these sources are at more reasonable, lower redshifts. On the other hand, we notice that no counterpart is found for these sources, so they might also possibly be spurious detections in the SCUBA2 residual image, in some cases.

The cosmic volume sampled by COSMOS, adopting a generous redshift range 4 < z < 8, is 7 × 10⁷ Mpc³, and the SFRD from our population in this volume is about $1\times {10}^{-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{3}$. This seems quite low with respect to the total SFRD at these redshifts (Madau & Dickinson 2014; L18), as expected due to the shallower depths reached. We are likely still sampling the high-end tail of the luminosity function, albeit now with pretty fair statistics.

Ivison et al. (2016) selected a sample of 109 dusty star-forming galaxies with Herschel colors (S_{500 um} ≥ 30 mJy, i.e., S_{500 um}/S_{250 um} ≥ 1.5 and S_{500 um}/S_{250 um} ≥ 0.85) in a 600 deg² survey, estimated 32% of that sample to have z =4–6, and reported a number density of ${\rho }_{z\gt 4}\approx 6\times {10}^{-7}$ Mpc⁻³. The sources of Ivison et al. (2016) have a median infrared luminosity (LIR) ∼ 1.3 × 10¹³ L_⊙, while our z > 4 sample reaches to fainter levels with a median LIR ∼ 8.0 × 10¹² L_⊙ and down to LIR ∼ 3.6 × 10¹² L_⊙. We detected 64 dusty star-forming galaxies at z = 4–6 (σ_{500 μm} ∼ 3 mJy), implying a much higher space density than that in Ivison et al. (2016). Given that the completeness of our z > 4 sample is not yet constrained, detailed values for the space densities will be reported in a future paper.

Among the 85 candidates, there are six sources with confirmed spectroscopic redshift z > 4. Four of them are listed in Table 2, as they are fitted with an AGN component, while five of them are shown in Appendix D with their actual "zspec" marked in the SED panels. The last source with spectroscopic redshift is the source ID20007898 (i.e., AzTEC/C1 in Smolčić et al. 2012), as shown in the fourth panel of Figure 23, whose spectroscopic redshift z_spec = 4.7 has been reported by Brisbin et al. (2017). As a test, we ignored the redshift of ID20007898 and fit its SED over z = 0–10 to derive a photometric z_phot,FIR = 4.8 ± 0.25. This agrees well with the z_spec, suggesting that our SED procedure is reasonable for giving photometric redshifts on z > 4 dusty galaxies. We have run the same test on IR-detected sources with z_spec > 3 and shown the redshift comparison in Figure 24. Our IR SED-driven photometric redshift is in good agreement with the spectroscopic redshift for galaxies with significant FIR/(sub)mm detections, while it is perhaps somewhat underestimated in the case of galaxies with important AGN torus emission. The relevance of this hint is limited, of course, by low-number statistics. We find that in order to obtain a fair comparison between photometric and spectroscopic redshifts, the redshift uncertainties must be increased by adding in quadrature an extra component of about 10% of (1+z). This is of order of the systematic effect expected from dust temperature variations, as discussed in Section 7.2. The redshift errors associated with our sources currently only reflect the accuracy of their SEDs and do not include such systematic uncertainties. Once weighted by the (total) redshift uncertainties, the comparison in Figure 24 suggests that our FIR photometric redshifts are underestimated by 6% × (1 + z) (Δz ∼ 0.4), on average, confirming the idea that the FIR SED of GN20 is somewhat colder than the high-z average. Larger samples are required to bring these results to firmer footing.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The multiband cutouts and SEDs of all candidates are presented in Appendix D. Redshift confirmation of as many as possible of these candidates is clearly needed, and we are planning observations with NOEMA, ALMA, and hopefully JWST in the future. Such further studies of these high-z candidates will be needed to finally validate their selection and understand critical issues like completeness, spurious fraction, actual redshift distribution (hence number densities), and others. Our photometric catalog in COSMOS will be the basis for such explorations.