"Super-deblended" Dust Emission in Galaxies. I. The GOODS-North Catalog and the Cosmic Star Formation Rate Density out to Redshift 6

In our analysis, we fit source fluxes in the images using ${galfit}$ (Peng et al. 2002, 2010). We use PSF fitting (i.e., we treat the sources as unresolved). The IR data with the highest angular resolution²¹ are the Spitzer16 and $24\,\mu {\rm{m}}$ images, with FWHM 4''–6'', a scale at which the intrinsic sizes of most distant galaxies can be neglected. This might not be true for some low redshift galaxies in the images, but these are not the main focus of our efforts.

Since it is not practical to use ${galfit}$ to simultaneously fit a very large number of prior sources in a large image, our code will first divide the image under examination into small regions (boxes). When fitting prior sources in each box, we also consider prior sources from a buffer region of surrounding boxes to avoid edge effects. The buffer size is at least 2–3× the PSF FHWM of the image being analyzed. We performed tests to verify that, as long as the box size and buffer size are several times larger than the PSF FWHM, the ${galfit}$ results are unaffected by the exact choices; only the computing speed is affected. We run ${galfit}$ PSF fitting in each box and then combine all boxes to make the full source-model images and residual images. All our PSF fitting is performed at fixed R.A.–decl. positions, as determined from the IRAC catalog, after checking and correcting astrometric differences with each data set under consideration. However, in order to improve the fitting for bright sources, we perform a second-pass fit (using first-pass results as first guesses) allowing high ${\rm{S}}/{\rm{N}}$ (e.g., $\gtrsim 10$) sources to vary their positions. This second-pass exercise also allows us to evaluate the extra photometric noise (and possible flux biases) introduced by residual (uncorrected) astrometric distortions between IRAC and each IR data set. These terms are generally small but we correct for them in the analysis.

We fit simple PSFs at 24 μm rather than extended source models because the PSF FWHM at this band is about $5\buildrel{\prime\prime}\over{.} 7$, much larger than the typical sizes of $z\geqslant 0.2$ star-forming galaxies, which are ∼8 kpc or $2\buildrel{\prime\prime}\over{.} 5$ at z = 0.2 (e.g., Conselice 2014). For lower-redshift galaxies in the image, we caution that this approach will lead to underestimated flux measurements, for example, for ID 11828, which is the largest spiral galaxy at ${z}_{\mathrm{spec}}=0.106$ in GOODS-N, and for 30 more $z\lt 0.2$ FIR-to-mm detected galaxies in our final catalog. To properly remove the background²² in the 24 μm image, we run a first pass of PSF fitting, then median filter the residual image on a scale of 30 × 30 pixels, then subtract the smoothed background image from the original image. We run a second pass of PSF fitting with fixed prior source positions, and a third pass with varied prior source positions for the highest ${\rm{S}}/{\rm{N}}$ sources from the second-pass fitting. Then we obtain the final 24 μm fluxes and formal ${galfit}$ errors.

The GOODS-N field has very deep VLA ∼20 cm observations from F. Owen (2018, in preparation) and shallower observations from Morrison et al. (2010).²³ The Morrison et al. (2010) catalog covers a wider region (diameter of ∼15') than the circular area covered by the F. Owen (2018, in preparation) catalog (diameter of ∼9'). In this work, we use the radio data for two main purposes: to complete the list of prior sources at high redshifts, where galaxies might be too faint to be detected at 24 μm, and to help constrain the overall IR luminosity of each prior source based on the FIR-radio correlation.

The importance of including radio prior sources is illustrated in Figure 2, where we present a series of redshifted SEDs based on templates that we used in this work (see Section 3.1). All SEDs in the figure are normalized to a common ${S}_{350\mu {\rm{m}}}\,=6.22$ mJy, which is the $3{\sigma }_{350\mu {\rm{m}}}$ detection limit derived from our simulation (see the following sections). Normalizing to other SPIRE bands will lead to a similar conclusion: the expected ${S}_{20\mathrm{cm}}$ flux is always somewhat brighter than the empirical detection limit at 20 cm ($3{\sigma }_{20\mathrm{cm}}\approx 7.5\,\mu \mathrm{Jy}$), whereas we do not expect to detect the faintest SPIRE sources with $z\gtrsim 3$ at 24 μm. This is also illustrated in Figure 3, where we plot the predicted flux density at 24 μm and 20 cm as a function of redshift, for the same normalization at 350 μm. The radio data are effectively deeper than the MIPS 24 μm data for selecting sources at $z\gtrsim 3$, by roughly a factor of $\gtrsim 4$. The predicted radio flux densities are $\sim 2\times$ brighter than our 20 cm detection limit at any redshift, while the predicted 24 μm flux densities at $z\gtrsim 4$ are $\sim 2\times$ fainter than the detection limit. Therefore, by using the VLA 20 cm data, we can have a more complete prior source catalog especially for $z\gtrsim 3$ sources, improving our completeness for fitting PACS, SPIRE, and (sub)mm photometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the radio catalog from F. Owen (2018, in preparation) as the main source of 20 cm photometry in this study, with detections down to a 5σ level of significance. F. Owen (2018, in preparation) carefully derived best estimates for the flux, accounting for possible spatial extension of the sources by comparing fluxes in radio images of different resolution from $1^{\prime\prime}$ to $3^{\prime\prime}$. In order to match our prior catalog with IRAC-based positions to the radio we adopt the same matching scheme of F. Owen (2018, in preparation; their Equation (1)), with a tolerance depending on the radio ${\rm{S}}/{\rm{N}}$ and the sizes of the sources, and including a $0\buildrel{\prime\prime}\over{.} 2$ term in quadrature to account for typical IRAC astrometric accuracy.

In order to provide a radio measurement for sources weaker than the F. Owen (2018, in preparation) 5σ detection limit (and/or upper limit), we also performed prior-extraction photometry and Monte Carlo simulations independently on the radio imaging data of both Morrison et al. (2010) and F. Owen (2018, in preparation). Small astrometric misalignments (including nonlinear distortions) between the images under exam and the IRAC catalog have been corrected in our photometry by adjusting the prior source positions, at radio, and in all other bands discussed elsewhere in this paper. (Note that this is only important for 24, 16 μm, 20 cm and PACS image data, where the PSFs are still relatively small.) We first run ${galfit}$ to fit PSF models at fixed positions, and then we run ${galfit}$ again, allowing high ${\rm{S}}/{\rm{N}}$ source positions to vary by a small amount (maximum offset less than 1 pixel). We fit on the radio images without primary beam correction, and then calculate the primary beam correction according to each prior source position and apply it to each photometry measurement.²⁴

Note that the PSF of Owen's radio image is about ∼2'' (provided by F. Owen especially for this work), and ∼17 for Morrison's radio image.²⁵ Some sources might be resolved in the radio data, typically at low redshift. In order to verify the performance of our photometry, we compare the measurements derived with our procedure to those from the F. Owen (2018, in preparation) catalog. We find that there is a good agreement in general: the median of the ratio between the different measurements is very close to 1, while the semi-interquartile range is about 6% for bright sources (>80 μJy), and reaching 10% for fainter fluxes down to the 5σ limit. This includes effects of flux losses due to over-resolution of the sources, and is accurate enough for a statistical use of the fluxes from our procedure below the 5σ limit, despite its simplicity in assuming unresolved emission and fitting at fixed positions.

As a result of these measurements, we find 1334 IRAC sources with ${\rm{S}}/{{\rm{N}}}_{20\mathrm{cm}}\gt 3$ (554 from the Owen 2017 catalog), of which 112 are not detected with ${\rm{S}}/{\rm{N}}$ of 3 at MIPS 24 μm. The 24 μm undetected radio sources could be radio-loud Active Galactic Nuclei (AGNs), intermediate redshift (e.g., $z\sim 1.5$) galaxies with strong rest-frame $9.7\,\mu {\rm{m}}$ Silicate absorption (e.g., Magdis et al. 2011), higher redshift star-forming galaxies, although in some cases they could be spurious detections boosted by the noise in the radio data. In Figure 4, we plot the histogram of photometric redshifts for these sources. Indeed two peaks can be found, and the $z\sim 3$ peak likely represents the high-z star-forming galaxies missed by the MIPS 24 μm selection. We also note that reducing the radio ${\rm{S}}/{\rm{N}}$ threshold to 2.5 for the high-z (e.g., $z\gt 2$ in the figure) IRAC sources will almost double their number, while making the gap at $z\sim 2$ shallower (indicated by the cyan histogram in Figure 4). Thus in the next section we use 2.5 as the radio ${\rm{S}}/{\rm{N}}$ threshold for $z\gt 2$ sources.

Zoom In Zoom Out Reset image size

By combining detections at 24 μm and 20 cm, we produce our 24+radio prior source catalog, which will be used for the analysis in the next steps: SED fitting (Section 3.1) and photometry over the FIR-mm bands (Sections 3.2 –3.3). There are 2983 sources in the IRAC catalog with ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\geqslant 3$, plus 112 sources with ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\lt 3$ but with ${\rm{S}}/{{\rm{N}}}_{20\mathrm{cm}}\geqslant 3$.

In order to improve completeness for potential high-z prior sources, which might be detectable in the FIR and (sub)mm data due to their redshifted SED peaks, we also include lower ${\rm{S}}/{\rm{N}}$ sources (down to S/N = 2.5) with photometric redshift above 2.0 in the 24+radio prior source list. In this way, we selected a total of 3306 sources in the 24+radio prior source catalog, with these criteria: ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\geqslant 3$, or ${\rm{S}}/{{\rm{N}}}_{20\mathrm{cm}}\geqslant 3$, or ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\geqslant 2.5$ when ${z}_{\mathrm{phot}}\gt 2.0$, or ${\rm{S}}/{{\rm{N}}}_{20\mathrm{cm}}\geqslant 2.5$ when ${z}_{\mathrm{phot}}\gt 2.0$. A total of 47.9% of sources in this full 24+radio prior catalog have spectroscopic redshifts, but only 4.66% of sources have photometric redshifts ${z}_{\mathrm{phot}}\gt 2$.

Although the flux error measurements at 24 μm and 20 cm are fairly close to Gaussian (see Section 5.3 and the Appendices), starting from a large number of 19,437 IRAC sources would inevitably result in a number of noise-dominated spurious detections at these radio/mid-IR bands, even assuming simple Gaussian statistics. We expect a total of ∼50 spurious detections above $3\sigma$ combining the two bands (24 μm and 20 cm). Given that only ∼4000 IRAC priors have ${z}_{\mathrm{phot}}\gt 2$, they would also produce some additional ∼50 sources with spurious radio/mid-IR detection at ${z}_{\mathrm{phot}}\gt 2$ and S/N between 2.5 and 3. These $\approx 100$ objects would get included among our 3306 priors spuriously, lacking any actual 24 μm or radio detections. This number is small enough that it should not adversely affect the detectability of IR galaxies, as it only marginally increases the number density of priors fitted at all bands. In most cases, we expect these spurious priors to remain undetected in the FIR/mm bands following the super-deblending analysis.

For the remaining 16,131 IRAC sources that were not included in our prior list, we assume that their flux contributions to the PACS, SPIRE, and mm image data are negligible, and we do not consider them further in the rest of this work. Even though their fluxes will not be exactly zero, their presence will act as a background whose average level will be consistently taken into account by our procedure. Their possibly inhomogeneous distribution will act as additional confusion noise, which our procedure will also quantify.

Figures 2 and 3 suggest that our prior source list is fairly complete and will include, in principle, any FIR-mm detectable star-forming galaxy, even at the highest redshifts, thanks to our very deep 20 cm detection limit. This of course is no longer entirely true when scatter in the FIR-radio correlation and SED shape variations are considered, instead of assuming galaxies to be perfectly represented by Magdis et al. (2012) templates, as in Figures 2 and 3. Also, the presence of noise in the measurements might affect the detectability of sources close to the detection limits. Therefore, despite using the deepest 24 μm and radio catalogs, we do expect some incompleteness in our prior list for $z\gt 3$ galaxies at the faintest (but detectable) flux levels, especially over the SPIRE bands and in the (sub)mm. This motivates a posteriori addition of sources found in residual map (mainly) in these bands, as will be discussed in Section 3.4.

The combination of 24 μm and radio priors is also useful to overcome the limitations at 24 μm where we fit all IRAC priors with a density of about 1 source per beam. This will imply that the effective depth of our catalog at 24 μm will depend on local crowding (see further discussions about crowdedness in this paper). However, the beam at 20 cm is small enough that crowding is irrelevant. Using the radio data should reduce the possibility that we are missing useful priors in regions where 24 μm is most crowded.

A more sophisticated way to complete the prior sample at high redshifts would be to cull more sources from the IRAC catalog and to use the correlation between stellar mass and SFR to single out appropriate additional high redshift priors. We have carried out a preliminarily investigation of this idea in GOODS-N, and find that this does not provide an obvious improvement to what has been achieved in the current work. The stellar mass selection will be explored in more detail in future work in COSMOS (S. Jin et al. 2018, in preparation), where the prior catalogs at 24 μm and radio wavelengths are shallower than in GOODS-N.

Despite having a fairly small PSF (FWHM ≈ 45; hence a small number of IRAC sources per beam), the GOODS-N 16 μm imaging data observed by the Spitzer IRS Peak-Up Imaging (PUI) are a factor of $\simeq 2$ shallower (for the same SFR) than the MIPS 24 μm image data (∼7.5 μJy versus ∼5 μJy sensitivities; see Table 1). Therefore in this work we directly fit the 16 μm image data at all 24+radio prior positions to measure 16 μm fluxes. Only 40% of our 3306 24+radio prior sources have ${\rm{S}}/{{\rm{N}}}_{16\mu {\rm{m}}}\geqslant 3$.

Table 1.  GOODS-N "Super-deblended" Photometry Results

Band	Instrument	Beam FWHMa	${\rho }_{\mathrm{fit}}$ b	${N}_{\mathrm{fit}}$ c	${N}_{\mathrm{excl}.}$ d	${N}_{{\rm{S}}/{\rm{N}}\geqslant 3}$ e	${N}_{\mathrm{add}.}$ f	$1\,\bar{\sigma }$ g
		arcsec	beam−1					mJy
24 μm	Spitzer/MIPS	5.7	1.205	19437	0	3056	0	$5.165\times {10}^{-3}$
20 cm	VLA	1.7/2.0h	0.107	19437	0	1328	0	$2.744\times {10}^{-3}$
16 μm	Spitzer/IRS/PUI	3.6	0.082	3306	0	1335	0	$7.681\times {10}^{-3}$
100 μm	Herschel/PACS	7.2	0.326	3294	12	1178	0	0.315
160 μm	Herschel/PACS	12.0	0.862	3137	169	1153	18	0.681
250 μm	Herschel/SPIRE	18.2	0.973	1540	1766	668	13	1.571
350 μm	Herschel/SPIRE	24.9	1.142	966	2340	292	10	2.072
500 μm	Herschel/SPIRE	36.3	1.181	470	2836	125	17	2.570
850 μm	JCMT/SCUBA2	11.0i	0.154	668	2638	32	15	1.249
1160 μm	AzTEC+MAMBO	19.5	0.374	515	2791	43	11	0.661

Notes.

^aBeam FWHM is the full width at half maximum of the circular-Gaussian-approximation point-spread function of each image data. ^b ${\rho }_{\mathrm{fit}}$ is the number density of prior sources fitted at each band, normalized by the Gaussian-approximation beam area. ^c ${N}_{\mathrm{fit}}$ is the number of prior sources fitted at each band. ^d ${N}_{\mathrm{excl}.}$ is the number of prior sources excluded from fitting at each band. These sources are subtracted from original image with their SED-predicted flux density at each band. ^e ${N}_{{\rm{S}}/{\rm{N}}\geqslant 3}$ is the number of prior sources with ${\rm{S}}/{\rm{N}}\geqslant 3$ (i.e., detected) at each band. ^f ${N}_{\mathrm{add}.}$ is the number of ${\rm{S}}/{\rm{N}}\geqslant 3$ additional sources that are not in the prior source catalog but blindly extracted from the intermediate residual image product at each band (see Section 3.4). ^g $1\bar{\sigma }$ is the median of the flux density uncertainties of all sources detected with S/N $\geqslant 3$ in each band. ^hThe shallower VLA image data from Morrison et al. (2010) have a beam FWHM of $1\buildrel{\prime\prime}\over{.} 7$, and the deeper VLA image data from F. Owen (2018, in preparation) are especially produced with a beam FWHM of $2\buildrel{\prime\prime}\over{.} 0$. ⁱWe use the unfiltered image, which has a narrower PSF FWHM than that of the matched-filter convolved image. We measured FWHM $\approx 11^{\prime\prime}$ in our PSF reconstruction for the unfiltered image, which is roughly consistent with Figure 14 of Chapin et al. (2013), and compares with $14\buildrel{\prime\prime}\over{.} 8$ for the matched-filter image.

Download table as:  ASCIITypeset image

We consider four distinct SED components in fitting procedure. From shorter to longer wavelengths, these are: (1) a stellar component from Bruzual and Charlot (2003, hereafter BC03) with a Small Magellanic Cloud (SMC) attenuation law; (2) a mid-infrared AGN torus component from Mullaney et al. (2011); (3) dust continuum emission as modeled by the Magdis et al. (2012) library; and (4) a power-law radio component (see the text below).

The Magdis et al. (2012) library is based on Draine and Li (2007, hereafter DL07) templates fitted to the average SEDs of MS and starburst (SB) galaxies as a function of redshift. DL07 templates are parametrized by a number of physical properties: the minimum interstellar radiation field (ISRF) intensity ${U}_{\min }$, the maximum ISRF intensity ${U}_{\max }$, the dust mass ${M}_{\mathrm{dust}}$, the mass fraction of PAH to total dust mass ${q}_{\mathrm{PAH}}$, and the mass fraction of dust grains located in photo-dissociation regions (PDR) γ, and so on. Hence they provide quite a wide range of SED shapes. However, Magdis et al. (2012) simplified the parametrization of DL07 templates for galaxies at various redshifts with only two parameters: the IR luminosity per dust mass ${L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$ or the mean ISRF intensity $\langle U\rangle$,²⁶ and whether the source is on the MS or not. Magdis et al. (2012) find that the shapes of the dust SEDs of MS galaxies at a given redshift, as traced by ${L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$, is not expected to vary significantly with increasing ${L}_{\mathrm{IR}}$ or SFR. Meanwhile, ${L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$ changes as a function of redshift, and the variation of SED shape among MS galaxies is expected to depend only on ${L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$, or equivalently $\langle U\rangle$. Magdis et al. (2012) constructed a series of dust SED templates for MS galaxies parametrized only by $\langle U\rangle$, as shown in their Figure 16. These templates assumed an evolution of $\langle U\rangle \propto {(1+z)}^{1.15}$, no evolution in PDR fraction, and a small evolution of ${q}_{\mathrm{PAH}}$ beyond z = 2, as indicated by their data. Béthermin et al. (2015) updated the evolution of $\langle U\rangle \propto {(1+z)}^{1.8}$ based on the latest data in the COSMOS field. Therefore, in this work, we use the Magdis et al. (2012) dust SED templates (which depend on $\langle U\rangle$) with redshift evolution from Béthermin et al. (2015) (which determines $\langle U\rangle$ at each redshift) to fit galaxy SEDs and predict photometric redshift and FIR/mm fluxes. Examples of our SED templates are shown in Figure 2.

Meanwhile, at all redshifts, there is a small fraction of SB galaxies which have SFRs higher than those of MS galaxies with similar stellar masses. Their ${L}_{\mathrm{IR}}/{M}_{\mathrm{dust}}$ ratios or $\langle U\rangle$ are found to be higher than those for MS galaxies (at least at $z\lt 1.5$) and likely do not vary with redshift (Magdis et al. 2012; Tan et al. 2014; Béthermin et al. 2015). So their SEDs can be also parametrized in the same way, by a constant $\langle U\rangle$. When fitting MS galaxies at a given redshift, we allow for a range of $\langle U\rangle$ values spanning ±0.2 dex around the expected value, to allow for the observed scatter of dust temperature among MS galaxies (defined by z = 0 results; see Magdis et al. 2012).

A power-law radio SED has been added to describe the radio continuum as (${S}_{\nu }\propto {\nu }^{-0.8}$)²⁷ and tied to ${L}_{\mathrm{IR}}$ by the FIR-radio logarithmic ratio ${q}_{\mathrm{IR}}\equiv {\mathrm{log}}_{10}({S}_{\mathrm{IR}}$, 8–1000 μm/$3.75\,\times {10}^{12}\,{\rm{W}}\,{{\rm{m}}}^{-2})$ − ${\mathrm{log}}_{10}({S}_{1.4\mathrm{GHz}}$/${\rm{W}}\,{{\rm{m}}}^{-2}\,{\mathrm{Hz}}^{-1})$ (e.g., Helou et al. 1985; Condon 1992; Ivison et al. 2010; Sargent et al. 2010; Tan et al. 2014), so that the radio data point can be used in the fit together with FIR/mm to constrain ${L}_{\mathrm{IR}}$. Here we adopt the latest results that suggest a slowly evolving ${q}_{\mathrm{IR}}=2.35\,\times {(1+z)}^{-0.12}+\mathrm{log}(1.91)$, following Magnelli et al. (2015) and Delhaize et al. (2017).²⁸ We find that such a small evolution is also generally warranted by our global SEDs. We verified that adopting a constant ${q}_{\mathrm{IR}}$ would have negligible impact on the photometric analysis carried out in this paper.

3.1.1. SED Fitting Parameters

In the case of radio-excess sources or radio-loud AGNs, which are outliers of the FIR-radio correlation, the radio emission does not predominantly arise from star formation. In these cases the radio photometry should not be used in SED fitting. Also, we find that mid-IR rest-frame emission is at times hard to fully reproduce with models that fit the FIR emission. For this and other analogous situations in the following, we introduce three parameters which summarize the SED fitting approach and which are included in our final catalog: (1) a flag to distinguish between MS and SB galaxies; (2) a flag to identify radio-excess; and (3) a flag for high-quality FIR photometry (in which case radio and 24 μm fluxes are not used in the SED fitting and flux prediction; they are un-necessary and can possibly be misleading if they are affected by a radio-loud or mid-IR AGN).

The first parameter is the MS/SB classification (noted as "Type_SED"), corresponding to the use of MS/SB types of SED templates. This is determined by carrying out the SED fitting twice, as we need a first estimate of SFR to decide if the galaxy is MS or SB. In the first pass, we fit all MS and SB SED templates for each source and derive the best fitting values and uncertainties for ${z}_{\mathrm{phot}}$, $\langle U\rangle$, ${M}_{\mathrm{dust}}$, and ${L}_{\mathrm{IR}}$, based on the equations in Press et al. (1992). SFRs are computed from the integrated ${L}_{\mathrm{IR}}$, assuming a Chabrier IMF (Chabrier 2003) and the relation $\mathrm{SFR}={L}_{\mathrm{IR}}/(1\times {10}^{10}\,{L}_{\odot }){M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Daddi et al. 2010b). We then compare the $\mathrm{SFR}$ and its uncertainty ${\sigma }_{\mathrm{SFR}}$ with the MS-based SFR (SFR_MS) according to Equation (A1) of Sargent et al. (2014): ${\mathrm{SFR}}_{\mathrm{MS}}={M}_{* }\times N({M}_{* })\times \exp \left(\tfrac{A\cdot z}{1+B\cdot {z}^{C}}\right)$, where the redshift z comes from the SED fitting and the stellar mass M_* is from optical/near-IR catalogs, as described in Section 2 when available. The parameter values ${N}_{{M}_{* }}$, A, B, and C are taken from Sargent et al. (2014, appendix) and are described therein. Sources are set to pure SB (Type_SED = 1) if $\mathrm{log}(\mathrm{SFR}/{\mathrm{SFR}}_{\mathrm{MS}})\gt 0.6\,\mathrm{dex}$ and $\mathrm{SFR}/{\sigma }_{\mathrm{SFR}}\gt 3$, and to pure MS (Type_SED = −1) if $\mathrm{log}(\mathrm{SFR}/{\mathrm{SFR}}_{\mathrm{MS}})\lt 0.4\,\mathrm{dex}$ and $\mathrm{SFR}/{\sigma }_{\mathrm{SFR}}\gt 3$. For all other sources (e.g., those without M_* estimate or those for which the SFR does not meet the criteria above), we assume that we cannot conclusively decide whether they are MS or SB galaxies, and we set Type_SED = 0, using all MS and SB templates to search for the best fit to their SEDs.

The second parameter is the classification of radio-excess sources or radio-loud AGNs (noted as "Type_20 cm"). This is also determined using a two-pass SED fitting approach. In the first pass, we do not use the observed radio flux density (${S}_{\mathrm{OBS},\mathrm{radio}}$). Since the radio SED component is tied to the ${L}_{\mathrm{IR}}$ by the FIR-radio correlation, we can derive an expected radio flux density (${S}_{\mathrm{SED},\mathrm{radio}}$) and associated uncertainty from the IR SED fit. Then, the (dis)agreement between ${S}_{\mathrm{OBS},\mathrm{radio}}$ and ${S}_{\mathrm{SED},\mathrm{radio}}$ determines our choice of the parameter Type_20 cm: if $({S}_{\mathrm{OBS},\mathrm{radio}}-{S}_{\mathrm{SED},\mathrm{radio}})/\sqrt{{\sigma }_{\mathrm{OBS},\mathrm{radio}}^{2}+{\sigma }_{\mathrm{SED},\mathrm{radio}}^{2}}\gt 3$ and ${S}_{\mathrm{OBS},\mathrm{radio}}\gt 2\times {S}_{\mathrm{SED},\mathrm{radio}}$, then we set Type_20 cm = 1; otherwise Type_20 cm = 0. Then, in the second pass, if a source has Type_20 cm = 1, we do not use its radio data, while for the remaining sources we use the radio data to provide tighter constraints on the physical properties. Overall, we have flagged 10.3% of the 24+radio source as radio-excess/radio-loud candidates. Our approach is very conservative, flagging a source as a possible radio AGN when both the bolometric ${L}_{\mathrm{IR}}$ and radio are well measured and when the radio flux exceeds the value predicted by the bolometric ${L}_{\mathrm{IR}}$ by a factor of two or more, which is comparable to the scatter of the radio-IR correlation. Some of these sources will not actually be AGNs, but we refrain from using their radio fluxes for SED fitting when we have a well-determined IR SED.

The third parameter is the combined ${\rm{S}}/{\rm{N}}$ computed over the FIR/mm bands (noted as "Type_FIR"). For sources with FIR/mm combined ${\rm{S}}/{\rm{N}}\geqslant 5$,²⁹ we set Type_FIR = 1 and do not fit the 24 μm and radio data points anymore. In this way our SED fitting procedure will not be affected by the scatter of the FIR-radio relation, or by galaxy–galaxy variations of the mid-IR dust SED features (e.g., variations of the IR to rest-frame 8 μm ratio, IR8, Elbaz et al. 2011). This optimization improves the FIR/mm part SED fitting of individual sources, but does not lead to obvious difference in the statistical results in Section 8.

As has already been briefly mentioned, the SED fitting is (re)run for each FIR/mm band prior to extracting photometry for that band. For example, before we measure PACS 100 μm fluxes, we run SED fitting to the K_s, IRAC, 16 μm, 24 μm, and 20 cm (if not radio-excess/radio-loud) data points for each 24+radio prior source. Similarly, before extracting SPIRE 350 μm photometry, we run SED fitting with the aforementioned bands plus the newly measured 100, 160, and 250 μm data points. The purpose of this SED fitting is to provide the best possible prediction for choosing the final prior sources that are used for fitting. This is a crucial step toward obtaining super-deblended photometry.

For sources that have reliable spectroscopic redshifts, we fix the redshifts to those values for the SED fitting. For those with optical/near-IR photometric redshifts, we perform SED fitting, allowing the redshift values to span a range of ±10% in $(1+z)$, corresponding to an uncertainty range of about $\pm 2\sigma$ (Skelton et al. 2014; Pannella et al. 2015). We do not allow redshifts to vary beyond this range, even if the minimum ${\chi }^{2}$ were to be found at the edge of the interval, to avoid getting solutions that are too inconsistent with the optical constraints. For IRAC sources without any available optical/near-IR photometric redshift (∼10% of the total IRAC catalog, but only 2.3% of the 3306 24+radio priors³⁰ ), we allow the full range $0\lt z\lt 8$ for the SED fitting. The IR-to-radio SED redshifts are presented as ${z}_{\mathrm{IR}}$ in our catalog (e.g., Table 2).

Table 2.  GOODS-N "Super-deblended" Photometry Catalog (Example)

ID	R.A.	Decl.	${z}_{\mathrm{IR}}$	${z}_{\mathrm{type}}$	S500	${\sigma }_{500}$	${\rm{S}}/{{\rm{N}}}_{\mathrm{IR}}$	$\mathrm{log}{M}_{* }$	${\mathrm{SFR}}_{\mathrm{total}}$	${\mathrm{SFR}}_{\mathrm{IR}}$	${\sigma }_{{\mathrm{SFR}}_{\mathrm{IR}}}$	$\mathrm{goodArea}$	T20 cm	TSED
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
564	189.2994995	62.3700371	4.055	s	44.0	5.8	19.9	11.25	1812.4	1726.0	86.6	1	0	1
11499	189.2606049	62.2172966	4.50	p	10.4	1.8	18.2	11.44	983.5	977.4	53.7	1	0	−1
2592	189.3670044	62.3222733	4.74	p	14.5	2.3	16.3	11.37	1238.9	1231.0	75.4	1	0	−1
15289	188.9900818	62.1734276	3.075	s	8.7	2.3	15.4	8.06	386.1	354.9	25.5	0	0	1
16332	189.2442017	62.1585388	3.55	p	11.3	3.8	14.9	⋯	659.8	594.8	163.5	1	0	0
3532	189.3077393	62.3073006	3.99	p	9.1	2.2	14.8	11.46	689.1	687.9	190.2	1	1	−1
9053	189.0358124	62.2431564	3.48	p	14.7	1.8	14.4	11.26	579.2	576.1	88.5	1	0	−1
4500	189.4094086	62.2934532	3.190	s	14.6	2.5	14.1	10.26	450.7	317.7	22.5	1	0	1
18911	189.1133881	62.1015778	4.55	p	31.8	3.9	14.1	⋯	1323.3	1193.0	100.4	0	0	0
4990	189.1330261	62.2873993	4.16	p	12.6	2.3	13.8	10.32	698.5	687.5	49.8	1	0	1
14914	188.9597931	62.1782951	5.33	p	20.1	2.6	13.6	12.55	1760.0	1724.0	132.1	0	0	−1
15213	189.0117798	62.1743240	3.24	p	5.9	2.9	12.3	10.45	399.7	384.2	148.9	1	0	0
17624	189.0368347	62.1344681	3.36	p	12.1	4.0	11.7	11.46	512.1	511.2	88.0	0	0	−1
16810	189.0555725	62.1504669	3.05	p	10.3	1.6	11.4	10.84	338.3	327.8	71.1	0	0	−1
13107	189.1156006	62.1996002	3.49	p	9.2	1.9	11.3	11.07	433.2	427.9	128.4	1	0	−1
130	189.3960114	62.3908005	4.13	p	23.7	2.5	11.0	10.67	1109.8	1085.0	98.5	0	0	1
3827	189.5537872	62.3028259	3.90	p	29.0	3.4	10.8	⋯	1083.8	977.1	264.0	0	0	0
16121	189.1436768	62.1616745	3.74	p	12.2	2.1	10.6	11.24	454.0	446.2	42.1	1	0	−1
17381	189.2108612	62.1394119	3.10	p	3.8	9.3	10.5	10.68	491.3	445.9	160.4	1	0	0
18603	189.0650787	62.1119728	3.72	p	12.6	6.6	10.4	11.04	742.8	737.7	163.4	0	0	0

Note. Here we show a few example columns and a few example ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\geqslant 10$ sources at $z\geqslant 3$ (sorted by the ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}$). Column (1)–(3) are the IRAC catalog from GOODS-Spitzer Legacy Program (see Section 2). Column (4) is our IR-to-radio SED redshift, which is spectroscopic redshift if available and otherwise the SED fitted photometric redshift; see Section 3.1. Column (5) indicates whether ${z}_{\mathrm{IR}}$ is spectroscopic redshift ("s") or photometric redshift ("p"). The references for each spectroscopic redshifts are fully listed in the released machine-readable catalog. Column (6) and (7) are flux density and uncertainty in unit of mJy. Column (8) is the FIR+mm combined S/N in Equation (2). Column (9) is the stellar mass from 3D-HST Skelton et al. (2014) or Pannella et al. (2015), converted to Chabrier IMF when needed. Column (10) is the total SFR, adding UV-unattenuated SFR (Section 8.2) to the dust-obscured SFR. Column (11) and (12) are the dust-obscured SFR and uncertainty from FIR+mm SED fitting (see the Figure 23 caption for deriving SFR from ${L}_{\mathrm{IR}}$). Column (13) is the "goodArea" parameter, equals 1 if the source is in the inner lower rms area in 24 μm (has deeper prior catalog for surrounding sources), and otherwise is 0 (has shallower or even incomplete prior catalog for surrounding sources). Column (14) is the "Type_20 cm" parameter, equals 1 if the source is classified as radio excess, and its radio data point is not fitted in SED fitting, and otherwise is 0. Column (15) is the "Type_SED" parameter, equals 1 if the source is classified as starburst type (its SED is fixed to using starburst templates), or −1 if classified as MS type (its SED is fixed to using MS templates), and otherwise is 0 (allowing fitting all SED templates).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

For each band, the SED fitting procedure returns a most likely expected flux density (${S}_{\mathrm{SED}}$), as well as its uncertainty (${\sigma }_{\mathrm{SED}}$) based on the ${\chi }^{2}$ statistics. This procedure requires that we can associate a well-behaved, (quasi-)Gaussian uncertainty to the fluxes measured in all shorter-wavelength bands. The description of the derivation of reliable flux uncertainties will be presented in Section 5.

Figure 5 shows the number of prior sources per beam with ${S}_{\mathrm{SED}}+2\,{\sigma }_{\mathrm{SED}}$ above threshold values indicated by the X axis. We add twice the flux density uncertainty to the best-fit flux density in order to conservatively estimate a plausibly maximum flux density for each source. For each band we show a black vertical line that indicates the median $1\sigma$ flux density uncertainty that we finally obtain from our procedure (see Table 1).

Zoom In Zoom Out Reset image size

Using these figures, we can optimally determine the threshold value ${S}_{\mathrm{cut}}$ for our analysis for each band. We recall that sources with ${S}_{\mathrm{SED}}+2{\sigma }_{\mathrm{SED}}$ larger than ${S}_{\mathrm{cut}}$ are included in the list of prior sources for ${galfit}$ photometric measurements in a given band (henceforth, selected sources, for that band). Sources with ${S}_{\mathrm{SED}}+2{\sigma }_{\mathrm{SED}}$ fainter than ${S}_{\mathrm{cut}}$ in a given band are excluded from ${galfit}$ fitting in that band (henceforth, excluded sources), and their flux densities, as predicted from the SED fitting, will be subtracted from the image, which will be discussed in the next section. A larger value of ${S}_{\mathrm{cut}}$ will lead to less confusion in ${galfit}$ fitting (i.e., a smaller number density of fitted sources, ${\rho }_{\mathrm{beam}}$) but a larger number of fainter sources to be subtracted, and vice versa. Also, setting too large a value for ${S}_{\mathrm{cut}}$ might prevent us from obtaining measurements for galaxies with detectable FIR/mm emission. We therefore set ${S}_{\mathrm{cut}}$ by requiring both ${\rho }_{\mathrm{beam}}\lesssim 1$ and also satisfying $1\sigma \lesssim {S}_{\mathrm{cut}}\lesssim 3\sigma$ at each band, where σ is the empirical median noise per beam at this band calculated from simulations as described below. The lower boundary is set to $1\sigma$ to avoid fitting sources that are too faint for the data in hand. The upper boundary of $3\sigma$ is still faint enough (particularly considering that we have added $2{\sigma }_{\mathrm{SED}}$ to the SED flux densities when evaluating the cutoff threshold) that we do not run much risk of discarding priors for any sources that are individually detectable in the data.

Regardless of these criteria, we also always keep as a prior any source that was significantly detected in the previous (shorter-wavelength) band with ${\rm{S}}/{\rm{N}}\gt 3$. In this way, for sources that have high ${\rm{S}}/{\rm{N}}$ in the preceding band but that are predicted to be too faint in band currently under analysis, we can still obtain an upper limit that is useful in SED fitting.

Figure 5 shows the effect of imposing ${S}_{\mathrm{cut}}$ thresholds on each band. The number of the selected sources is reduced to only ∼30% of the total number (3306) at 350 μm and 500 μm bands (overcoming blending), and to ∼15% at 1.16 mm (which avoids fitting overly faint sources).

Sources that are excluded from ${galfit}$ PSF fitting through the previous section analysis do still contribute some weak flux in the observed images, and could thus boost the photometry for sources that are retained for fitting. Therefore we first make a model image of those excluded sources by fixing their fluxes to the predicted values, and then subtract the faint-source-model image from the observed original image and obtain a faint-source-subtracted image. The ${galfit}$ photometry is then performed on the faint-source-subtracted image for the selected prior sources in a similar way, as described in Section 2.2. In Figure 6, we present the SPIRE 250 μm images in each stage of this procedure. In the Appendix we show all of the images produced during these steps for all of the bands analyzed in this paper. It is interesting to compare the total flux subtracted from the original images (for faint excluded sources) with that actually extracted from the fitted sources. This reaches up to 30% of the total flux in the worst cases for the longer-wavelength SPIRE bands, where only 15%–20% of the original priors are actually fitted, and is generally much smaller for other bands.

Zoom In Zoom Out Reset image size

The list of prior sources that are selected for fitting extends to sufficiently faint (i.e., $\lt 3\sigma$) flux densities for each band that we do not expect all sources used for the PSF fitting to actually be detected with S/N $\gt 3$ in the extraction process. The number of S/N $\gt 3$ detections (${N}_{{\rm{S}}/{\rm{N}}\gt 3}$) is generally of order of 1/3 of all fitted sources (${N}_{\mathrm{fit}}$). ${N}_{\mathrm{fit}}$ and ${N}_{{\rm{S}}/{\rm{N}}\gt 3}$ of all bands are listed in Table 1.

Although the 24+radio combination provides a very deep prior source catalog for fitting the FIR/mm bands, some less massive, radio-quiet high-z galaxies could still be missed. For example, in Figures 2 and 3, the z ∼ 3–5 SEDs have 24 μm and 20 cm fluxes that are close to or even below the detection limits. Due to the diversity of galaxy properties (i.e., PAH fraction or FIR-radio correlation scatter), some of these sources will be missed in our 24+radio prior list, but might have detectable SPIRE or (sub)mm fluxes.

We carry out a "blind" search (i.e., unguided by priors) for additional sources by running SExtractor (Bertin & Arnouts 1996) on the residual images. We combine any new sources detected this way with the previous list of priors, and then re-run ${galfit}$ to see if the new, blindly selected sources can be detected.

We limit our blind-extraction region to fall within the 24+radio catalog coverage. This is done by measuring the perimeter of the 24+radio sources and setting pixels outside of the perimeter to not-a-number (NaN) values. Then we run SExtractor on the masked residual image with a detection threshold of about 2.5–3.0 (depending on the band, with a GLOBAL background) to detect all possible residual sources. We subject the resulting SExtractor catalogs to careful visual inspection. We run ${galfit}$ again, fitting over all previously selected prior sources in that band plus the new sources from blind detection that pass the visual inspection. The new residual source candidates are finally kept and included in our catalog if their final ${\rm{S}}/{\rm{N}}$ is above 3, typically only 1/4 of the initial SExtractor candidates.

We illustrate this blind detection in Figure 7, using SPIRE 350 μm data. Panel (1) and panel (2) show the same residual image before extracting additional sources, and the green circles in panel (2) highlight the additional sources which are extracted by SExtractor and then measured with ${galfit}$ to have ${\rm{S}}/{\rm{N}}$ $\gt 3$. Panel (3) shows the new residual image after we run the final ${galfit}$ fitting, including prior sources and the additional sources that were retained. We list the number of final detected additional sources at each band in Table 1.

Zoom In Zoom Out Reset image size

We have decided not to include these additional sources in the scientific analysis presented in the final sections of this paper for a number of reasons. Their positions are uncertain, and in most cases it is not trivial to associate them with optical/near-IR counterparts. Also, a certain fraction of these blindly detected sources are probably spurious, often a superposition of a number of fainter sources simulating a brighter object. It is unlikely that these detections are residuals from poor fits to brighter sources, as this is rarely seen in the FIR/mm residual images. Nevertheless, there are also likely genuine $z\gt 3$ sources among these objects—for example, the z = 5.3 submillimeter galaxy HDF850.1 (Walter et al. 2012) is the additional source ID 850011 in our additional source catalog.

We simulate one source at a time in the faint-source-subtracted image.³¹ The simulated source positions are generated randomly within the survey area, and do not avoid other real, bright sources. Their input flux densities, ${S}_{\mathrm{in}}$, are uniformly distributed in log space within the range from $\sim 1\sigma$ to $\sim 30\sigma$, where σ is the mean flux density uncertainty at each band (see Table 1).

We add each extra simulated source to the 24+radio prior list, which has already been filtered for that band (i.e., with all selected sources, see Section 3.2, as well as additional sources detected in the residual images, see Section 3.3), and then perform the PSF fitting photometry including all active priors and an extra prior at the position of the simulated source. We also perform two-pass ${galfit}$ fitting as described in Section 2.2 to allow variation in the positions of priors in the second pass if they have high ${\rm{S}}/{\rm{N}}$ in the results of the first pass. At the end of each pass, we run the corrections described in the following Sections 5.2 and 5.3. The ${galfit}$ output flux density of the simulated source is ${S}_{\mathrm{out}}$, and the ${galfit}$ flux density uncertainty is ${\sigma }_{{galfit}}$.

Finally, this process is repeated ∼3000 times. Therefore, for each band, we have ∼3000 values for ${S}_{\mathrm{in}}$, ${S}_{\mathrm{out}}$, and ${\sigma }_{{galfit}}$. In addition, we measured several additional properties for each simulated source: for example, the local rms noise value ${\sigma }_{\mathrm{rms}\mathrm{noise}}$ for each source, measured in the instrument noise image data; the local absolute flux density in the residual image (hereafter ${S}_{\mathrm{residual}}$), measured by computing the total absolute pixel values in the PSF aperture for each fitted prior source; the local scatter in the residual image, measured also in the PSF circle; and the crowdedness parameter, computed by summing up the Gaussian weighting of all sources at current source i position,

$$\begin{eqnarray}&&\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 current source i and ${d}_{\mathrm{PSF}}$ is the FWHM in arcsec of the PSF for the band being analyzed. Defined 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 measurable parameters are designed to provide key information about the quality of fitting and the local crowding (hence blending) of prior sources. We will make use of them with the simulation results in order to calibrate the best possible flux bias and flux uncertainty corrections.

In addition to the uncertainties estimated in this way, we add in quadrature to each object's flux uncertainty the appropriate contribution (when relevant) to account for astrometric dispersion when fitting at fixed spatial positions. This term is estimated by computing the dispersion in measured fluxes for bright sources at fitted position versus those at the positions optimized in the fit.

The statistics of the differences between input and output fluxes, ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$, can be used to verify and correct the bias of the ${galfit}$ flux measurements. In Figure 10 we present the analysis of ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$ based on simulation data for SPIRE 350 μm as an example. In panel (1), ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$ is plotted against the simulation input flux density ${S}_{\mathrm{in}}$, which covers the range 3–12× in terms of the detection limit ${\sigma }_{350\mu {\rm{m}}}$. The overall median of ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$ is generally small, indicating that the (constant) background levels used in the ${galfit}$ work are quite accurate.

Zoom In Zoom Out Reset image size

For all bands we explore the dependence of ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$ on the following three key properties: (1) the flux uncertainty ${\sigma }_{{galfit}}$, which is a direct output of ${galfit}$ fitting; (2) the residual flux ${S}_{\mathrm{residual}}$, measured in the output residual image; and (3) the crowdedness, which indicates prior source blending as described in the previous section. When testing parameters (1) and (2), we further use the rms noise ${\sigma }_{\mathrm{rms}\mathrm{noise}}$ to normalize these two parameter values in order to account for the local rms variation in case of non-uniform depth of the data over the whole GOODS-N region. For example, if sources at the map center have smaller rms noise than sources near the edges, the outer sources will likely have higher ${\sigma }_{{galfit}}$ and ${S}_{\mathrm{residual}}$ simply because of their locations.

As shown in panels (2)–(4) of Figure 10, we bin simulated sources by each of the key properties. In each bin we compute the median of ${S}_{\mathrm{in}}-{S}_{\mathrm{out}}$, which is just the flux bias S_bias. Then we fit S_bias in each bin by a three-order polynomial function of the X axis parameter. The fitted polynomials are shown as the red curves in panels (2)–(4). Note that the quantities along the X axes are measurable for each real fitted source; thus these functions are also applicable to each real source.

The Monte Carlo simulations also provide a path toward obtaining reliable estimates of flux uncertainties. The flux error returned by ${galfit}$ (${\sigma }_{{galfit}}$) only reflects the fitting uncertainty in a formal way (e.g., Peng et al. 2002). In fact, the ${galfit}$ formalism is optimized for optical/near-IR images and assumes uncorrelated noise (all pixel are independent) in the data. Thus it cannot fully account for real uncertainties, especially in the presence of correlated noise and confusion noise. Imperfect flux uncertainties will strongly affect the validity of SED fitting for individual sources and the assessment of detections or lack thereof.

We analyze the statistics of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})/{\sigma }_{{galfit}}$ with simulation data, which should have a rms dispersion of 1.0 by definition if uncertainties are well defined and Gaussian-like. Taking the 350 μm simulation as an example, in the left upper panel of Figure 11 we bin simulated objects by the measurable parameter (${\sigma }_{{galfit}}/{\sigma }_{\mathrm{rms}\mathrm{noise}}$) as in previous section, and in the middle and right upper panel we show the histograms of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})/{\sigma }_{{galfit}}$ (in blue) and $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})/{\sigma }_{1}$ (in red), where ${\sigma }_{1}$ is the first-step corrected flux uncertainty. The width of the blue histogram is clearly much broader than 1.0, indicating that the values ${\sigma }_{{galfit}}$ are substantially underestimated.

Zoom In Zoom Out Reset image size

We determine correction factors as a function of the X axis parameter, requiring that, in bins defined for quantities in the X axis, the rms dispersion of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})/\sigma$ is equal to 1.0. In each bin, we compute the rms of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})$ and the median of ${\sigma }_{{galfit}}$, and then calculate their ratio, which is exactly the factor that we need for correcting ${\sigma }_{{galfit}}$ to the real uncertainty of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})$. As in the previous section, we fit the correction factor as a three-order polynomial function against the X axis parameter. The fitted polynomial functions are shown as the red curves in the figure. In this way, we obtain the first-step corrected ${\sigma }_{1}$, whose histogram is shown in red, and is much closer to a width of 1.0.

The first-step corrected uncertainty histograms can sometimes have asymmetric shapes, which indicate that there are still some systematic biases. Thus we take the first-step corrected uncertainty ${\sigma }_{1}$ as input and perform the second and the third steps. In the second step, we first perform the flux bias correction with parameter ${S}_{\mathrm{residual}}/{\sigma }_{\mathrm{rms}\mathrm{noise}}$, and then perform the uncertainty correction with the same parameter as illustrated by the middle row panels of Figure 11. Then, similarly, in the third-step we apply the flux bias correction, followed by the flux uncertainty correction with the parameter crowdedness. Therefore, with these correction recipes, we are able to apply the correction to real fitted sources.

We have also tried multiple iterations or different combinations of the three steps. For example after the three-step processing, we use the output corrected flux and uncertainty as the input for several more iterations of the three-step processing. However, we find that the corrected values and histograms are stable and cannot be further improved in this way.

Similar figures for all examined bands are reported in the Appendix B (Figures 46–52). Generally, for objects whose ${galfit}$ uncertainties are around or below the median for that band, we find that the flux uncertainties are often quite underestimated, while for objects with the largest uncertainties (typically those affected by substantial deblending uncertainties), the values are often overestimated (i.e., the real errors are smaller as judged from our simulations). An analogous effect is seen for the crowdedness parameter. We have analyzed a number of other parameters (including flux itself, spatial positions, redshift, different kind of crowdedness and residual definitions, etc.), and we found no other significant dependency to further correct flux uncertainties.

An important step in our procedure is to predict fluxes for sources in a given band before actually performing measurements in that data set. Sources with predicted fluxes (plus twice flux prediction uncertainty) brighter than ${S}_{\mathrm{cut}}$ are retained to be fitted with ${galfit}$. For those fitted sources, we can now compare their measured fluxes to the SED predictions in order to investigate the presence of any potential bias in those predictions (and, by extension, also the predictions for the faint/excluded/subtracted sources).

These comparisons are shown in Figure 13. In the left panels we compare ${S}_{\mathrm{predict}}$ and ${S}_{{galfit}}$, and in the right panels we show the difference $({S}_{{galfit}}-{S}_{\mathrm{predict}})$ divided by the total uncertainty, plotted versus ${S}_{\mathrm{predict}}$. The error bars represent the combined uncertainties of ${\sigma }_{{galfit}}$ and ${\sigma }_{\mathrm{predict}}$.

Zoom In Zoom Out Reset image size

${S}_{\mathrm{predict}}$ and ${S}_{{galfit}}$ generally agree well with each other. There is a small fraction of sources with ${S}_{\mathrm{predict}}$ lower than ${S}_{{galfit}}$ at the faint end, as seen in the left panels (but less obvious in the right panels when actual flux uncertainties are taken into account). This small bias is driven by the requirement to have a significant ${galfit}$ detection in order to be shown in the plot, and vanishes when this selection bias is taken into account.

Another important overall check of our procedure is evaluating the goodness of fit of the global SEDs that we derived with our adopted template libraries, to verify at the same time that our templates are appropriate and that the error bars are sensible. We concentrate here on the FIR+mm range between 100 and 1100 μm, which is the range most affected by blending that we are interested to test, and also because mid-IR and radio deviations from our SEDs might be expected because galaxies can be more complicated than the models (AGN components with different shapes/slopes, strength of PAHs, etc.). As we have already shown with examples of SEDs (see, e.g., Figure 9), the reduced ${\chi }^{2}$ values are generally reasonable and close to unity. We verified this aspect more systematically for the whole sample, considering the overall IR detection ${\rm{S}}/{\rm{N}}$ ratio (see Equation (2) in Section 8). For the 518 galaxies with ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\gt 10$, we find a median ${\chi }_{r}^{2}=1.0$, while a median ${\chi }_{r}^{2}=0.8$ is found for the $5\lt {\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\lt 10$ sources. About 2% of all galaxies in either sample cannot be fit by our models, with ${\chi }_{r}^{2}\gt 4$. In most cases this appears to be due to an erroneous or problematic redshift (photometric or spectroscopic), either due to an error or because of blending/superposition between foreground and background objects, because of gravitational lensing or simply line-of-sight alignment. In some cases it is quite clear that only part of the FIR/mm emission is due to the galaxy under consideration as the IRAC prior, and that an additional component comes from a higher redshift object that contributes additional flux at longer wavelengths. Many of these could be scientifically interesting cases (e.g., lenses or the high redshift background galaxies).

Here we compare our 24 μm, PACS, and SPIRE flux measurements and uncertainties with those of Elbaz et al. (2011) and HerMES (DR3; Roseboom et al. 2010),³² and our 1.16 mm catalog with that of Penner et al. (2011).

In Figure 14, we compare our 24 μm flux and uncertainty histograms with the Elbaz et al. (2011) catalog. In the left panel, the flux histogram of ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\geqslant 3$ sources in this work is shown in blue, and that of Elbaz et al. (2011) is shown in red. The shaded area under each histogram indicates the sources with FIR+mm multiband combined ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\geqslant 5$ and ${\rm{S}}/{{\rm{N}}}_{24\mu {\rm{m}}}\geqslant 3$, where ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}$ is defined as per Equation (2) in Section 8.1. Elbaz et al. (2011) do not include (sub-)mm bands (i.e., for their sources the combined ${\rm{S}}/{\rm{N}}$ only involves 100–500 μm photometry). For the bright sources, the flux histograms from both works are consistent, while for fainter sources, the flux histogram from this work shows more prior sources. We have about 450 more 24 μm detected candidate prior sources for the longer-wavelength photometry than those used in Elbaz et al. (2011). The shaded area, or the FIR to mm and 24 μm detected sources, also has a larger number of detections in our catalog. The newly detected sources are fainter. Notice how IR-detected sources fall off very rapidly below 100 μJy, suggesting that our prior sample is quite complete and possibly redundant for Herschel sources at typical redshifts (although we might still be somewhat incomplete for the highest redshifts, as previously noted in Section 2.4). The right panel shows flux uncertainties, corrected based on simulations (see Section 5.3). They are generally higher than those in Elbaz et al. (2011), but should have less bias and be more representative of the true uncertainties in the measurements (see Sections 5.3 and B). For the FIR and 24 μm detected sources in the shaded area, our flux uncertainties are also generally higher.

Zoom In Zoom Out Reset image size

In Figure 15, we compare the 100 μm flux and uncertainty histograms from this work to those from the catalog of Elbaz et al. (2011). In the left panel, our work has about 200 more sources and extends to a slightly fainter flux range. In the right panel, after applying the simulation-based flux uncertainty correction (see Section 5.3), our flux uncertainties are lower than those in Elbaz et al. (2011). This is due both to the greater number of fainter sources in the left panel, and to our simulation-based corrections to the flux uncertainties.

Zoom In Zoom Out Reset image size

Other PACS and SPIRE bands have similar histograms to that at $100\,\mu {\rm{m}}$ and hence are not shown here. The comparisons of our 1.16 mm flux measurements and uncertainties with those of Penner et al. (2011) are shown in Figure 16. We use their deboosted fluxes when possible. The number of detected sources is small at this wavelength, but our catalog has 13 more detections. The sources span a similar range of flux in both works. After applying the simulation-based flux uncertainty corrections (see Section 5.3), the flux uncertainties in this work are slightly higher than those in Penner et al. (2011). Similar to the case at 24 μm, we believe that the difference is mainly due to the fact that our corrections yield more reasonable flux uncertainties (based on the simulation results).

Zoom In Zoom Out Reset image size

The one-to-one flux comparisons for the three SPIRE bands and at 1.16 mm are presented in Figures 17 and 18. We cross-match common sources (within a limiting separation of 1'') and compare the fluxes from this work to those from Elbaz et al. (2011) or the HerMES catalog. The majority of the sources show consistent fluxes, but some sources in the Elbaz et al. (2011) catalog have higher SPIRE fluxes, which is due to source blending. The two most extreme outliers are marked with their IRAC IDs in the first panel. They are ID 14896 (ID 2132 in Elbaz et al. 2011) and ID 10611 (ID 1512 in Elbaz et al. 2011). ID 14896 is one of very few sources that have higher fluxes in our catalog than in that of Elbaz et al. (2011). This is due to the complex blending situation; Elbaz et al. (2011) assign most of the flux to a nearby source but little to this object. Figure 19 shows the multiband image cutouts around this source. The three sources, ID 14896, ID 14867, and ID 15027, become blended at 160 μm and longer wavelengths. Their SEDs are shown in Figure 20. For comparison, the red data points report values from Elbaz et al. (2011) for ID 14896 and ID 15027. ID 14867 has no counterpart in Elbaz et al. (2011). In the third panel, ID 15027, with a relatively low spectroscopic redshift z = 1.371, is not likely to have SPIRE flux as bright as the red data points. In the first panel, the 20 cm and 500 μm fluxes of ID 14896 suggest that its 250 μm flux is not likely to be as low as that reported by Elbaz et al. (2011). Thus, in our opinion, for this case our deblending leads to more reasonable results. For the other extreme case mentioned above (ID 10611 and its surrounding sources), we have also checked their image cutouts and SEDs. The blending situation is very similar to that described above, and in the catalog of Elbaz et al. (2011) most of the SPIRE flux is attributed to ID 10611 (e.g., like ID 15027). Thus the detailed figures and SEDs are not shown here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Similar considerations apply to the HerMES catalogs (see Figure 18). In this case, we note that there are large differences with our measurements at 350 and 500 μm, substantially larger than the uncertainties from both catalogs. Generally, we find that many sources in the HerMES catalog have their flux densities overestimated, likely because of unresolved blending. The number of single-band detected sources is almost doubled in our catalog compared to the HerMES catalog, for all SPIRE bands. The flux density uncertainties in HerMES catalog (using confusion noise computed from the inverted covariance matrix and the PSF matrix; see Roseboom et al. 2010, 2012) tend to be larger than ours.

Our method attempts to use all available information on galaxies in the region of sky under consideration in order to define an ideal set of priors to fit FIR/(sub)mm data, using the smallest possible number of useful sources. Increasing the number of priors will increase photometric errors for all galaxies on average because of increased blending. This effect is exemplified in Figure 21, where we show the normalized flux errors for the 100 and 250 μm bands as a function of the crowdedness parameter. Despite these two bands having quite different numbers of priors per beam (0.3 and 1, respectively; see also Figure 1), the normalized behavior of the flux error versus crowdedness is similar for the two bands. To first order, the relative flux error scales proportionally to crowdedness, although a second order polynomial fit is a better description for values of crowdedness $\gt 2$. At 250 μm, as a result of the much larger beam, we are fitting a smaller absolute number of priors (see Figure 21 and Table 1), which on average have higher values of crowdedness. The median value for the crowdedness is 1.64 (1.05) at 250 μm (100 μm), and the median relative flux error is similarly larger 1.63 (versus 1.02). The increase of the typical "noise" is the price to pay for fitting all reasonable priors. We note that this effect is not so dramatic, and errors are still quite acceptable, even if we are dealing with a density of priors approaching 1 per beam (in possible contradiction with results and discussions in Karim et al. 2013; Scudder et al. 2016).

Zoom In Zoom Out Reset image size

Figure 5 shows that, given the information in hand, it is not possible to reduce much further the number density of priors without failing to fit objects that could be potentially fairly bright in the FIR/mm imaging data. For example, raising ${S}_{\mathrm{cut}}$ from 3 to 6σ would reduce ${\rho }_{\mathrm{beam}}$ only from 1 to ∼0.8. On the other hand, if we were to reduce ${S}_{\mathrm{cut}}$ to 1.5σ, we would reach ${\rho }_{\mathrm{beam}}$ of about 2. It is interesting at this point to evaluate what the impact of these choices would be. Figure 22 shows the distribution of crowdedness values for the different cases. Reducing the source density from 1.0 to 0.8 per beam has a fairly minor effect on the implied crowdedness, and hence on the noise. Increasing the source density much beyond 1 has a much stronger effect: the fraction of sources in reasonable isolation, with crowdedness $\lt 1.5$, drops from 42% (${\rho }_{\mathrm{beam}}=1.0$) to 15% (${\rho }_{\mathrm{beam}}=2.0$), while the sources with poor measurements due to blending with crowdedness $\gt 2.5$ rise from 12% (${\rho }_{\mathrm{beam}}=1.0$) to 41% (${\rho }_{\mathrm{beam}}=2.0$). Therefore, if we were to accept more sources for fitting, including some of those below ${S}_{\mathrm{cut}}=3\sigma$, this would result in a significantly increase in the noise and poorer performance for all objects. Staying within the limit ${\rho }_{\mathrm{beam}}\lesssim 1$, as we have done for SPIRE in this paper, appears to be a reasonable choice.

Zoom In Zoom Out Reset image size

We are subtracting faint-source models from the images before running the photometric analysis. This procedure is prone to errors, as we might being over- (or under-)estimating the flux values to subtract. We expect that any such error would be automatically taken into account by our simulation-based approach. For example, if we were over-subtracting the fluxes of many faint sources, this would create, to first order, a negative (false) background, which we would detect and correct with our bias measurements, which in practice define the effective background in the data. To second order, the over-subtracted sources in this example, after zero-level determination, would result in increased background fluctuations and noise, because they are not likely to be spatially homogeneous, and they will have a variety of flux values. This extra noise would be reflected in our simulation-based calibration of uncertainties and thus taken into account in our error bars.

We are aware of a number of possible remaining problems that might still affect our flux measurements and their uncertainties.

One is the presence of additional sources (most often at high redshifts), which could significantly contribute flux in the data but which are not in our prior list. Such objects could spuriously boost the inferred flux (and hence lead to underestimated errors) of other priors. Judging from our simulations, this effect does not seem to be significant, since the distribution of $({S}_{\mathrm{in}}-{S}_{\mathrm{out}})/\sigma$ closely resembles a Gaussian with width equal to 1 (see Figure 11 and Appendix B).

Also, a fraction of photometric redshifts or even spectroscopic redshifts could be incorrect. Experience from various surveys have shown that a few percent of galaxies may have substantially incorrect redshift measurements or estimates, even when using values classified as reliable, as we do here. Note that in order to accept a spectroscopic redshift from the literature, we require that it is in agreement within 10% in $(1+z)$ with the photometric redshift, to avoid possible wrong redshifts. Using a wrong redshift can affect SED predictions. To some extent, allowing flexibility in dust temperatures (or, more specifically, in $\langle U\rangle$ values) moderates the impact because of the redshift-temperature degeneracy, especially if the redshifts are not catastrophically wrong. In this case, due to SED flux prediction uncertainties, there could be a few objects that are retained for ${galfit}$ fitting instead of being dropped. This would modestly increase crowdedness, but this is probably a negligible effect. Some objects might be erroneously dropped from fitting while still being bright enough for detection in the data. The latter cases could give rise to sources that would be detectable in the residual images if they are isolated. If they fall atop or near other priors, they may corrupt their photometry, biasing recovered fluxes high. Of course, in all cases, the impact of wrong redshifts would extend beyond simple photometric accuracy to the derivation of wrong bolometric ${L}_{\mathrm{IR}}$ and SFR, thus affecting the science.

Finally, another limitation is that we are using ${galfit}$ with a constant background. Angular variations in the background would act as increased noise. Carefully accounting for background variations could improve photometric uncertainties, but is limited by our knowledge of PSF wings, and of course by the very large beams in the SPIRE and (sub)mm imaging. We attempted a correction of background variations at 24 μm and 16 μm, where the beam is still fairly small and manageable. Dealing with this for PACS and SPIRE imaging is beyond the scope of the present work.

7.6.1. Some Considerations about Simulation Methodology

In Figure 23, we show the SFR versus redshift for all 3306 24+radio sources in the full GOODS-N field. The SFR are derived from the pure dust component IR luminosity over 8–1000 μm, excluding any AGN torus components as derived from our SED fitting (see Section 3.1). Among all the detected sources (henceforth, the FIR+mm sources), 106 sources have $z\geqslant 2.5$ (FIR+mm SED based photometric redshift or spectroscopic redshift when available), 70 have $z\geqslant 3$, and 20 have $z\geqslant 4$. Within the goodArea, we find 71 sources at $z\geqslant 2.5$, 49 at $z\geqslant 3$, and 14 at $z\geqslant 4$.

Zoom In Zoom Out Reset image size

In Figure 23 we indicate the SFRs of typical MS galaxies at three fiducial stellar masses (${M}_{* }={10}^{12}\,{M}_{\odot }$, ${10}^{11}\,{M}_{\odot }$, and ${10}^{10}\,{M}_{\odot }$ respectively), computed with the Equation (A1) of Sargent et al. (2014). This FIR+mm sample has a non-uniform distribution in both SFR and stellar mass. At lower redshift (e.g., $z\lt 1$), the limits of the SFR and stellar mass both rise quickly. At $z\lt 1$, the GOODS-N FIR+mm data can detect emission from MS galaxies with ${M}_{* }\lt {10}^{10}\,{M}_{\odot }$, but at $z\gt 1$ all detected galaxies are more massive. At higher redshift (e.g., $z\gt 3$), the sample is clearly biased toward high SFRs and stellar masses (${M}_{* }\gt {10}^{11}\,{M}_{\odot }$). Both the effective limits of SFR and stellar mass have risen by one to two orders of magnitude compared to those at $z\sim 1$.

We present our FIR+mm catalog in Table 2. The full machine-readable version can be accessed via the online journal; here we only show a few important columns whose descriptions are in the caption.

In order to estimate the total SFRs in our IR-detected galaxies, we need to add the unattenuated SFR that can be directly recovered in the rest-frame UV (${\mathrm{SFR}}_{\mathrm{UV}}$). We thus cross-matched the 3D-HST catalog for the rest-frame UV flux densities of the FIR+mm galaxies (${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\geqslant 5$). A total of 94% of them have rest-frame 1400 $\mathring{\rm A}$ flux densities from the 3D-HST catalog. The ${\mathrm{SFR}}_{\mathrm{UV}}$ is in general a small contribution (10% typically) compared to ${\mathrm{SFR}}_{\mathrm{IR}}$ for our IR selected sample. Figure 24 shows the ratio of ${\mathrm{SFR}}_{\mathrm{UV}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ versus stellar masses and redshifts (and in goodArea). The left panel shows sources with rest-frame 1400 $\mathring{\rm A}$ flux densities and stellar masses, colored by their redshifts. The right panel shows all FIR+mm sources (also in goodArea), with different colors indicating the method adopted to estimate ${\mathrm{SFR}}_{\mathrm{UV}}$, as discussed below. Method (1): if a galaxy is in the 3D-HST catalog (matched within $1\buildrel{\prime\prime}\over{.} 0$ radius) and has rest-frame 1400 $\mathring{\rm A}$ flux densities, then we use the Kennicutt (1998) calibration: ${\mathrm{SFR}}_{\mathrm{UV}}=1.4\times {10}^{-28}\times 4\pi {\left(\tfrac{{d}_{L}}{[\mathrm{cm}]}\right)}^{2}\left(\tfrac{{S}_{1400\mathring{\rm A} }}{[\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{Hz}}^{-1}]}\right)$. We also check that using instead the rest-frame $1700\,\mathring{\rm A}$ flux densities (also provided by 3D-HST) leads to no more than 0.1 dex difference in average. Method (2): if a galaxy has no rest-frame 1400 $\mathring{\rm A}$ but has a stellar mass, then we use a simple correlation between ${\mathrm{SFR}}_{\mathrm{UV}}/{\mathrm{SFR}}_{\mathrm{IR}}$ and stellar mass (as shown in the left panel of Figure 24, analogous to the stellar mass–attenuation relation that is fairly constant with redshift; Pannella et al. 2015) that we obtained with the galaxies in method (1): ${\mathrm{log}}_{10}({\mathrm{SFR}}_{\mathrm{UV}}/{\mathrm{SFR}}_{\mathrm{IR}})=13.3\mbox{--}1.4{\mathrm{log}}_{10}({M}_{* })$. This relation appears only to be valid for ${\mathrm{log}}_{10}({M}_{* })\geqslant 9.6$, and does not hold at very low stellar masses, where IR-detected objects are mainly SB sources (e.g., see the right panel of Figure 25). For the only two sources with ${\mathrm{log}}_{10}({M}_{* })\lt 9.6$ and no UV continuum, we just adopted the average ${\mathrm{SFR}}_{\mathrm{UV}}/{\mathrm{SFR}}_{\mathrm{IR}}$ value of the ${\mathrm{log}}_{10}({M}_{* })\lt 9.6$, UV-continuum-available sample (${\mathrm{SFR}}_{\mathrm{UV}}/{\mathrm{SFR}}_{\mathrm{IR}}=0.22$). Finally, for eight galaxies that do not have UV continuum nor stellar mass, as well as for a small sample of ∼28 objects where UV continuum is most likely dominated by the emission of an AGN (e.g., indicated by the mid-IR AGN SED fits and optical/near-IR images), we adopted the median ${\mathrm{SFR}}_{\mathrm{UV}}/{\mathrm{SFR}}_{\mathrm{IR}}$ of the whole sample (=0.11). In the following sections, we add ${\mathrm{SFR}}_{\mathrm{UV}}$ and ${\mathrm{SFR}}_{\mathrm{IR}}$ to obtain the total SFR (${\mathrm{SFR}}_{\mathrm{UV}+\mathrm{IR}}$) for each galaxy, and use the total SFR for all further analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 25 we compare our sample with empirical MS correlations in redshift bins up to z = 6. Only goodArea FIR+mm sources with stellar masses from 3D-HST or Pannella et al. (2015) are shown here; seven FIR+mm sources lack stellar mass estimates and are not shown. The differences between the three MSs considered here (from Sargent et al. 2014; Béthermin et al. 2015; Schreiber et al. 2015) are relatively large at both $z\lt 1.5$ and $z\gt 4$, but the different estimates are more consistent at $1.5\lt z\lt 4$. In the highest redshift bins, the Sargent et al. (2014) MS has the lowest SFR normalizations. The difference between these three MS relations (in some cases their extrapolations) can reach a factor of three at ${M}_{* }\sim {10}^{11.5}\,{M}_{\odot }$ and z ∼ 4–6, where global uncertainties and lack of direct IR data in the previous derivations still prevent a solid assessment of the MS level.

The red symbols in each panel indicate sources classified as strong starbursts during the SED fitting (Section 3.1), relative to the MS of Sargent et al. (2014). They are relatively rare, and mostly have masses ${M}_{* }\lesssim {10}^{10.5}\,{M}_{\odot }$. The fraction of SBs in this sample is 9.8% at $1\lt z\lt 4$, intermediate between the fraction of 30%–50% for submillimeter galaxy-like selection at the highest luminosities and the fraction of 2%–3% derived for a sample at $1.5\lt z\lt 2.5$ selected by stellar mass (complete to ${M}_{* }\geqslant {10}^{10}\,{M}_{\odot };$ Rodighiero et al. 2011 ³³ ).

The shaded area in each panel indicates the stellar mass range over which the sample is incomplete, evaluated by analyzing SMFs for star-forming galaxies as described in the next section. Our FIR+mm sample should have a high degree of completeness at the characteristic stellar mass ${M}_{\mathrm{star}}^{* }$ in the four lowest redshift bins, where ${M}_{\mathrm{star}}^{* }$ for star-forming galaxies is estimated to be in the range of ${10}^{10.3}$–${10}^{10.9}\,{M}_{\odot }$ at $0.2\lt z\lt 2.5$ according to Davidzon et al. (2017), Ilbert et al. (2013), and Muzzin et al. (2013).

At $z\gt 3$, estimates of ${M}_{\mathrm{star}}^{* }$ fall in the range ${10}^{10.74\mbox{--}11.56}\,{M}_{\odot }$ (Ilbert et al. 2013; Muzzin et al. 2013; Grazian et al. 2015; Davidzon et al. 2017). At these redshifts, the FIR+mm sample can only probe a small number of galaxies at the high-mass end of the SMF.

To evaluate the incompleteness of our FIR+mm sample we consider sources in each redshift bin, further dividing them into bins of stellar masses in Figure 26, and compare their distributions with various star-forming galaxy SMFs from the literature.

Zoom In Zoom Out Reset image size

The literature SMFs are convolved to a common stellar mass uncertainty (${\sigma }_{{M}_{* }}$, see labels in each panel) in order to account for measurement errors when comparing with our observed stellar mass distributions (blue and red histograms in the figure, representing the full FIR+mm sample and the SB sub-sample respectively). This is done by first using the Eddington bias-corrected (i.e., intrinsic) SMFs in these literature, then convolving them with a common ${\sigma }_{{M}_{* }}$ (which is redshift-dependent), following the approach described in Appendix A of Ilbert et al. (2013; except for Song et al. 2016 where the authors have not made that correction). We multiply by the goodArea co-moving volume to convert SMF (in units of number per volume) to the absolute number of galaxies within each stellar mass bin and each redshift bin.

The literature SMFs do not always agree with each other, especially at the low-mass end and at high redshifts, but in general they agree with our observed histograms up to the $3\lt z\lt 4$ bin without doing any fitting. Muzzin et al. (2013) presented two best-fits to the SMF at each redshift: one with fixed slope at the low-mass end, the other with non-fixed slope. Their fixed-slope SMFs agree better with the other two $z\lt 4$ SMFs in the literature, while the non-fixed-slope SMFs show large discrepancies at the low-mass end compared to other SMFs at $z\gt 2$. In order to probe the variety of SMFs, we show their non-fixed-slope SMFs in the figures. Note that the SMF of Davidzon et al. (2017) does not fully cover our highest redshift bin. We have linearly extrapolated their SMFs at $3.5\lt z\lt 4.5$ and $4.5\lt z\lt 5.5$ to our bin at $5\lt z\lt 6$ for purpose of comparison. Because their SMF is even higher at $4.5\lt z\lt 5.5$ than at $3.5\lt z\lt 4.5$, the extrapolated SMF is also higher.

To evaluate the stellar mass range (shown as the shaded area in Figure 25) over which our sample becomes incomplete in each redshift bin, we estimate at which value of stellar mass the number of observed galaxies in our sample starts to be substantially smaller than the number predicted by each SMF.

In the three bins with the highest redshifts in Figure 26, our FIR+mm sample tends to have more galaxies at the highest masses than would be predicted by the literature SMFs. In the $3.0\lt z\lt 4.0$ panel, the excess is not large, but in the $4.0\lt z\lt 5.0$ panel, our FIR+mm sample has eight more sources than is predicted by the SMFs of Grazian et al. (2015) and Song et al. (2016), or six more sources than predicted by the Davidzon et al. (2017) SMF. Note that in the GOODS-N field, there is a well-known z = 4.055 proto-cluster, which includes three FIR+mm sources detected in our catalogs: GN20, GN20.2a, and GN20.2b (Daddi et al. 2009b; Tan et al. 2014). Thus small number statistics is already an important factor at these high redshifts.

Moreover, Grazian et al. (2015) discussed the importance of cosmic variance (i.e., clustering of galaxies) between the GOODS-S and UDS fields in their Figure 3 and Section 4.2.4. The SMF derived from UDS sources is much higher than the SMF derived from GOODS-S sources at the high-mass end; the UDS SMF agrees better with the SMFs of Ilbert et al. (2013) and Muzzin et al. (2013) at similar redshifts. Hence cosmic variance may also be responsible for the excess in the highest redshift bins found in this work.

Our forthcoming work in the 2 square degree COSMOS field using the same approach (S. Jin et al. 2017, in preparation) will provide results that are less affected by cosmic variance in the high-mass and high redshift bins.

In Figure 27, we analyze the SFR contribution of our sample of galaxies in bins of stellar mass. We divide our sources into several stellar mass bins and then compute the sum of SFRs in each bin.

Zoom In Zoom Out Reset image size

For comparison, we use SMFs from the literature, first multiplying by the co-moving volume to obtain the absolute number of galaxies, and then multiplying by the SFR per unit mass predicted by the MS correlation to get the SFR contribution for each stellar mass bin (SMF-converted SFR). In this way, we are able to see how the SFR of this FIR+mm sample is distributed as a function of stellar mass, and to estimate what fraction of the co-moving SFR density this sample contributes to the whole population of star-forming galaxies at each redshift up to 6. Here we show the results computed assuming the MS correlation of Sargent et al. (2014), but the main results do not change substantially if the MSs of Béthermin et al. (2015) or Schreiber et al. (2015) were used instead.

The shape of these SMF-converted SFR histograms usually displays a peak around ${M}_{\mathrm{star}}^{* }$. In the three lowest redshift bins, the predicted SFR distribution based on the Muzzin et al. (2013) SMF has the best agreement with our observations, but implies a higher SFR contribution from the most massive galaxies (${M}_{* }\gt {10}^{11.0}\,{M}_{\odot }$), where we might be incomplete due to the limited volume of our survey. The SFR distribution predicted from the Davidzon et al. (2017) SMF matches our data well at the high-mass end, but is too low overall. The prediction based on the Ilbert et al. (2013) SMF falls in between the other two.

The FIR+mm sample is already accounting for most of the global SFR density at these redshifts. In the intermediate redshift bins ($1.5\lt z\lt 4.0$), the Ilbert et al. (2013) prediction gives the best fit to the normalization. The Muzzin et al. (2013) distribution gives a better fit to the shape at the low-mass end, which is probably due to the different stellar mass limit in Muzzin et al. (2013; e.g., ${M}_{* ,95 \% \mathrm{lim}.}={10}^{10.54}\,{M}_{\odot }$ at $2\lt z\lt 2.5$) compared to that in Ilbert et al. (2013; e.g., ${M}_{* ,90 \% \mathrm{lim}.}={10}^{10.01}\,{M}_{\odot }$ at $2\lt z\lt 2.5$), keeping in mind the current difficulties in determining the low-mass end slope of SMFs (see Appendix C of Muzzin et al. 2013). In the two highest redshift bins, the SFR distribution predicted from the Davidzon et al. (2017) SMF seems to be the closest one to the observed SFR histogram, but the SFR excess of our FIR+mm sample is still obvious.

At redshifts $z\gt 3$, our understanding of the cosmic SFR density is still largely limited to UV- and optical-based studies. Direct knowledge of dust-obscured star formation at such high redshift (i.e., from the FIR to mm data) mainly comes from a small number of rarest and most ultraluminous galaxies, leaving considerable uncertainty about true SFR densities at those high redshifts (Madau & Dickinson 2014). Recent attempts to directly measure dust-obscured SFR densities from FIR and submm data have been pushed to redshift $3\lt z\lt 4$ (e.g., by Gruppioni et al. 2013 and Burgarella et al. 2013), and to redshift 6 (e.g., by Rowan-Robinson et al. 2016 and Bourne et al. 2017).

In this and the next section, we first directly compute the inferred SFR densities from our FIR+mm goodArea sample up to redshift 6 (Figure 28) and compare with results from the literature, and then estimate the incompleteness corrections using literature SMFs (Figure 29). Finally, we repeat the comparison after correcting for incompleteness in our sample (Figure 30).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use the non-parametric 1/V_max method (Schmidt 1968) to estimate the SFR density of our FIR+mm sample in each redshift bin in order to eliminate the Malmquist bias. For each source, we compute the farthest redshift z_max at which it can still be detected in our catalog (i.e., ${\rm{S}}/{{\rm{N}}}_{\mathrm{FIR}+\mathrm{mm}}\geqslant 5$). Then for each redshift bin (lower and higher boundaries are z₁ and z₂), we compute the co-moving volume for each source (V_max) from z₁ to $\min ({z}_{\max },{z}_{2})$ (i.e., the smaller value; e.g., Pozzetti et al. 2003). Then the SFR density of this sample is the sum of all individual $\mathrm{SFR}/{V}_{\max }$ in each redshift bin. To account for more realistic uncertainties, we compute $\mathrm{SFR}/{V}_{\max }$ in a boot-strapping method: we randomly select the same number of sources in each redshift bin, allowing duplication, and then repeat 1000 times to compute the mean and rms (i.e., the uncertainty of the measured SFR density).

Figure 28 shows the SFR densities directly measured from our observed sample, and compares them to previous results from the literature (Gladders et al. 2013; Madau & Dickinson 2014). The fluctuations from bin to bin at $0\lt z\lt 3$ (where we expect a smooth redshift dependence) are significant at the $\sim 2\sigma$ level and suggest total uncertainties in our determinations of order of ∼0.2 dex, probably dominated by cosmic variance (although note that the bins do not have equal co-moving volumes).

We estimate the completeness of our measured SFR density from the FIR+mm sample to the total cosmic SFR density in each redshift bin in the following way, with the help of the SFR histograms in bins of stellar masses. We assume that SMF can predict the total cosmic SFR density based on the assumption of MS correlation (e.g., Sargent et al. 2014); thus the SMF-converted SFR histogram should be able to fit the observed SFR histogram when the FIR+mm sample is complete at certain stellar mass bin (e.g., as indicated by Figure 27). The fitting is done by adjusting the normalization of SMF-converted SFR histogram (see Figure 53 in the Appendix), so as to achieve a minimum ${\chi }^{2}$ in the stellar mass range where the FIR+mm sample is at least 50%-complete in stellar mass. The 50%-complete mass (${M}_{* ,50 \% \mathrm{lim}.}$) is determined by the renormalized version of Figure 26 (Figure 53 in the Appendix): we calculate above which stellar mass the number of FIR+mm galaxies starts to be larger than half of what is predicted by SMF (i.e., ${N}_{\mathrm{FIR}+\mathrm{mm}}\geqslant 0.5\times {N}_{\mathrm{SMF}}$).

We then integrate the renormalized SMF-converted SFR histogram down to ${10}^{8}\,{M}_{\odot }$ (as was done in Ilbert et al. 2013 and Muzzin et al. 2013) to obtain an estimate of the total cosmic SFR density. The ratio of the total observed SFR to the total cosmic SFR gives the completeness fraction (${f}_{\mathrm{comp}.}$).

Note that, as mentioned at the end of Section 8.2, the directly observed SFRs for our FIR+mm sample are SFR_UV+IR. We emphasize that the aim of the incompleteness correction is to account for the SFR_UV+IR from fainter/lower-mass galaxies that are not included in our sample. Although one might consider correcting the observed IR-only SFRs in the same fashion, we caution that this may not be appropriate if fainter/lower-mass galaxies have less dust attenuation, as is observed at z < 3 (e.g., Pannella et al. 2015).

In Figure 29, we present the cosmic SFR completeness fractions versus redshift, estimated using three different SMFs from the literature. We present the incompleteness-corrected total SFR densities in Figure 30.

At $z\lt 2$, the completeness of our sample drops linearly with $\mathrm{log}(1+z)$, and the three SMF-based results are fully consistent. However, at $z\gt 2$, the three SMF-based estimates differ by large factors, largely due to the different low-mass slopes of the literature SMFs. The SMF from Muzzin et al. (2013) has a shallower slope at low masses, and fits our observed histograms much better, leading to a completeness factor of about 80% at $z\sim 3$. The SMFs from Ilbert et al. (2013) and Davidzon et al. (2017) have steeper low-mass slopes, and hence lead to much larger incompleteness corrections (see Table 3).

Table 3.  GOODS-N FIR+mm Star Formation Rate Densities

${z}_{\mathrm{low}}$	${z}_{\mathrm{high}}$	${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs},\mathrm{UV+IR}}$	${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs},\mathrm{IR}-\mathrm{only}}$	${\sigma }_{{{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs}}}$	${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{corr}.}$	SMF${}_{\mathrm{corr}.}$	${f}_{\mathrm{comp}.}$
0.0	0.2	−1.852	−2.00	0.164	−1.852	Muzzin	100.0%
					−1.852	Ilbert	100.0%
					−1.852	Davidzon	100.0%
0.2	0.4	−1.818	−1.93	0.114	−1.805	Muzzin	97.0%
					−1.830	Ilbert	102.8%
					−1.776	Davidzon	90.7%
0.4	0.6	−1.349	−1.42	0.077	−1.315	Muzzin	92.5%
					−1.331	Ilbert	96.0%
					−1.293	Davidzon	88.0%
0.6	0.8	−1.531	−1.56	0.096	−1.476	Muzzin	88.1%
					−1.482	Ilbert	89.3%
					−1.462	Davidzon	85.4%
0.8	1.0	−1.199	−1.22	0.067	−1.137	Muzzin	86.8%
					−1.137	Ilbert	86.8%
					−1.125	Davidzon	84.2%
1.0	1.5	−1.263	−1.30	0.062	−1.199	Muzzin	86.4%
					−1.191	Ilbert	84.8%
					−1.183	Davidzon	83.2%
1.5	2.0	−1.383	−1.41	0.075	−1.239	Muzzin	71.8%
					−1.233	Ilbert	70.8%
					−1.150	Davidzon	58.5%
2.0	2.5	−1.170	−1.19	0.057	−1.098	Muzzin	84.8%
					−0.983	Ilbert	65.0%
					−0.806	Davidzon	43.3%
2.5	3.0	−1.534	−1.54	0.146	−1.509	Muzzin	94.4%
					−1.322	Ilbert	61.5%
					−1.161	Davidzon	42.3%
3.0	4.0	−1.353	−1.38	0.086	−1.263	Muzzin	81.3%
					−1.120	Ilbert	58.6%
					−0.762	Davidzon	25.7%
4.0	5.0	−1.619	−1.64	0.145	−1.313	MD14	∼49.4%
5.0	6.0	−2.040	−2.07	0.164	−1.610	MD14	∼37.1%

Note. The measured SFR density (${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs}.,\mathrm{UV}+\mathrm{IR}}$ and ${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs}.,\mathrm{IR}-\mathrm{only}}$, in units of ${M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}$) and uncertainty (${\sigma }_{{{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs}}}$) for the observed FIR+mm galaxies in this work as described in Section 8.6, as well as the incompleteness-corrected cosmic SFR density (${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{corr}.}$) and the completeness percentage in SFR (${f}_{\mathrm{comp}.}$, i.e., $={10}^{({{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{obs}.,\mathrm{UV}+\mathrm{IR}}-{{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{corr}.})}\times 100 \%$) as described in Section 8.7. We use literature SMFs and the Sargent et al. (2014) MS correlation to estimate the completeness percentage; therefore each SMF leads to a corrected value of ${{\rm{\Psi }}}_{\mathrm{SFR},\mathrm{corr}.}$ at each redshift bin. The SMF${}_{\mathrm{corr}.}$ column indicates which SMF is used (Ilbert et al. 2013; Muzzin et al. 2013; Grazian et al. 2015; Song et al. 2016; Davidzon et al. 2017), except for the last two rows at $z\gt 4$, where it becomes prohibitive to estimate any incompleteness correction. Given the evidence that our FIR+mm sample does not substantially overlap with the UV-based LBG samples at $z\gt 4$ (see Section 8.9), we correct the $z\gt 4$ cosmic SFR densities by adding the observed SFR densities from our FIR+mm sample to the UV-based, attenuation-corrected SFR densities from LBG samples (i.e., the values reported in Madau & Dickinson 2014, MD14). The percentages ∼49.4% and ∼37.1% thus represent the SFR ratio of our sample to the sum of our sample and the LBGs (Madau & Dickinson 2014, MD14).

Download table as:  ASCIITypeset image

At $z\gt 4$, the completeness fractions that we estimate from the renormalized SMFs are very low (see Appendix, Figure 53), ranging from $\sim 0.2 \%$ to $\sim 5 \%$. If we do not re-normalize the SMFs but just integrate the SMF-based SFR histogram (i.e., as shown in Figure 27), we would get a total cosmic SFR density that is fully consistent with the fit by Madau & Dickinson (2014) at $z\sim 5$, and would not reflect the excess SFR that we see in the observed SFR histograms. The main reason is that at $z\gt 4$ our FIR+mm sample is almost entirely disjoint from the UV-based samples that have traditionally been used to estimate the cosmic SFRD, as we will demonstrate in the next sections. Thus to estimate the total cosmic SFR density at this high redshift range, we need to add our directly observed FIR+mm SFR density to the UV-based SFR density. In Figure 30, we show the summed total cosmic SFR densities at $z\gt 4$ as the red open squares. Our directly observed FIR+mm SFR densities are also shown as the red arrows for comparison. We list all corrected cosmic SFR density values in Table 3.

Our corrected measurements imply total SFR densities that are consistent with those reported in Madau & Dickinson (2014) at $0.5\lt z\lt 2$, but possibly up to 0.4 dex higher at $2\lt z\lt 4$, where the results are less reliable and depend strongly on the assumed SMF. However, we also note that several recent FIR, (sub)mm and radio studies also seem to find similar trends. In Figure 30, we show the recent cosmic SFR density measurements from Rowan-Robinson et al. (2016) and Novak et al. (2017) for comparison. Based on much wider (20 ${\mathrm{degree}}^{2}$) but shallower FIR surveys, Rowan-Robinson et al. (2016) studied the FIR luminosity function with 3035 SPIRE $500\,\mu {\rm{m}}$ sources. Their cosmic SFR densities (the gray filled squares) are much higher than the UV estimates at $3.5\lt z\lt 6$ (i.e., the Madau & Dickinson 2014 curve) by factors of ∼2–6×. They assumed a fixed 500 μm luminosity function at $z\gt 3.5$ (Saunders et al. 1990, functional form with fixed parameters $\alpha =1.2$ and $\sigma =0.60$), but whether such fixed form is applicable is still controversial. Novak et al. (2017) used deep VLA imaging data in the 2 deg² COSMOS field to derive and integrate the 3 GHz luminosity function (assuming pure luminosity evolution) up to z = 5–6, and estimated the cosmic SFR density assuming a redshift-dependent FIR-radio correlation. Their results are shown with gray triangles in Figure 30, and are in good agreement with the Madau & Dickinson (2014) curve and the average of our incompleteness-corrected data points at $z\lt 4$. However, at $z\sim 5$, their cosmic SFR density is $\sim 3\times$ the value from Madau & Dickinson (2014).

Recently, Bourne et al. (2017) and Dunlop et al. (2017) have investigated dust-obscured SFR densities by stacking submm and mm data for faint galaxies selected from deep HST imaging data. Bourne et al. (2017) used SCUBA2 data at $850\,\mu {\rm{m}}$ and $450\,\mu {\rm{m}}$ covering 230 arcmin² in three of the CANDELS fields, stacking emission for galaxies selected from the 3D-HST catalog of Momcheva et al. (2016). Dunlop et al. (2017) used ALMA 1.3 mm data covering 4.5 arcmin² in the HUDF, stacking for galaxies selected from the deep HUDF HST data. They convert the stacked single (sub)mm band flux into dust-obscured SFRs using submillimeter galaxy template SEDs, and then sum this with the raw UV SFRs, uncorrected for dust attenuation. At $z\lt 2$, the total cosmic SFR densities that they derive agree well with the compilation of Madau & Dickinson (2014), similar to our own findings here. At higher redshift, their results are ∼0.1–0.2 dex larger than the Madau & Dickinson (2014) curve, but still ∼0.2–0.3 dex lower than the present work. It is possible that the most luminous IR-bright galaxies at $z\gt 2\mbox{--}3$ were not part of their stacked samples, which were selected based on near-IR (1.6 μm) HST data.

As a further step, we make a simple parametric fit to our incompleteness-corrected data points (i.e., all colored data points in Figure 30) using the Madau & Dickinson (2014) double power-law and Gladders et al. (2013; see also Abramson et al. 2016) log-normal equations. From the ${\chi }^{2}$ distribution statistics, we compute the median value of each parameter within the ${\chi }_{\min }^{2}$ to ${\chi }_{\min }^{2}+d{\chi }^{2}$ range,³⁴ and then obtain the following newly fitted equations and parameters:

$$\begin{eqnarray}{\rm{\Psi }} & = & {{\rm{\Psi }}}_{0}\displaystyle \frac{{(1+z)}^{3.0}}{1+{[(1+z)/2.9]}^{5.6}}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}\\ {{\rm{\Psi }}}_{0} & = & 0.00587\ (\mathrm{Chabrier}\ \mathrm{IMF})\end{eqnarray} \tag{ 3 }$$

(the red thick solid curve in Figure 30), and

$$\begin{eqnarray}{\rm{\Psi }} & = & \displaystyle \frac{{A}_{0}}{t\sqrt{2\pi {\tau }^{2}}}\exp \left[\displaystyle \frac{-{(\mathrm{ln}(t)-{T}_{0})}^{2}}{2{\tau }^{2}}\right]\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\,{\mathrm{Mpc}}^{-3}\\ & & t\ \mathrm{is}\ \mathrm{the}\ \mathrm{cosmic}\ \mathrm{age}\ \mathrm{in}\ \mathrm{Gyr},\ {T}_{0}=1.50,\tau =0.66\\ & & {A}_{0}=0.575\ (\mathrm{Chabrier}\ \mathrm{IMF})\end{eqnarray} \tag{ 4 }$$

(the red thick dashed curve in Figure 30).

The new curves are similar to the previous fits at $z\lt 1$. The double power-law (solid curve) is higher at $z\sim 2$, but the log-normal (dashed curve) peaks at values similar to the original fits at $z\sim 2$, where our SMF-converted SFR densities agree well with the literature data. At $z\gt 3$, both new curves can generally fit the new FIR/radio-based data points. But they diverge in the highest redshift bin, where Rowan-Robinson et al. (2016) measured a very high value, but we derive a lower value of the FIR SFR density. There might still be some less extreme dusty galaxies, with lower infrared luminosities and SFRs, that are missed in both our FIR+mm sample and in UV-based samples. Such objects might be recovered by future ALMA/NOEMA surveys.

Using our new FIR+mm photometric catalog, the completeness-corrected cosmic SFR densities are quite well constrained up to redshifts z ∼ 2–3. Variations up and down around the average trend reflect small sample statistics together with cosmic variance in the relatively small GOODS-Northfield.

At $z\gt 3$, despite our concerted effort on photometry, the uncertainties are still much larger, not only due to cosmic variance but also due to other factors that we briefly discuss here. There are remaining limitations in the photometry, even for bright sources due to residual blending. For example, we have evaluated the incompleteness in recovering IR sources as a function of their SFR/masses and redshifts. Figures 26 and 27 show that this is not likely to be a big issue at $z\lt 3$ for bright galaxies (while we correct for faint galaxies), but this clearly becomes more important at $z\gt 3$ (Figure 29). This is further exemplified by the "additional sources" that we found and neglected in the SFRD study (Section 3.4). Assuming that some of these are not spurious, some of them are likely to have redshifts $z\gt 3$, as shown by the example of HDF850.1, which is not included in the SFRD calculation. The limited sensitivity of the IR data also is a factor, permitting direct detections of sources only with SFR $\gt 500\,{M}_{\odot }\,{{yr}}^{-1}$ at $z\gt 4$. This may be far brighter than the threshold at which UV selection starts to include the large majority of galaxies, although the flux range for which this is true is currently unknown at $z\gt 3\mbox{--}4$. In any case, this may be another source of significant incompleteness, even for UV+IR samples at these redshifts. Meanwhile the accuracy of photometric redshifts can also be a problem. The photometric data used to estimate redshifts are very good for well-studied deep survey fields like GOODS-North, but catastrophic photometric redshift failures might still occur (e.g., in the case of red galaxies with apparent photometric redshifts, $z\gt 3\mbox{--}6$). It is noteworthy that some of the most distant and luminous star-forming galaxy candidates that we find, with $z\gt 5$ (see next sections), are also embarrassingly massive (see Figure 31). If their redshifts are correct, which must be tested with spectroscopy, these objects could pose interesting challenges for galaxy formation models, or may demonstrate limitations in our ability to derive accurate stellar masses for such objects (see Section 8.10 for more discussions about these sources).

Zoom In Zoom Out Reset image size

We have found a substantial population of IR-detected galaxies all the way out to redshifts z = 3–6, and have estimated that their contribution to the SFR density of the universe is significant and comparable to or larger than that derived from samples selected at UV rest-frame wavelengths. It is therefore interesting to discuss the possible overlap between galaxy samples selected at IR and UV wavelengths in order to evaluate whether they trace independent components of the cosmic SFR density (which should then be summed to derive the total), or if we are measuring the same population.

At $z\gt 3$, the 912 Å Lyman limit and 1216 Å Lyman α forest spectral breaks are redshifted through optical passbands, causing the galaxies to become very faint and to "drop out" at bluer wavelengths while still being clearly detected in redder filters. Optical color–color diagrams can then be used to select samples of LBGs in specific redshift ranges based on these "dropout" color signatures. For the HST ACS filters that were used to observe GOODS-N, most B-dropout LBGs fall into the redshift range $2.8\leqslant z\leqslant 4.4$, while V-dropout LBGs are found to have $z\approx 4.4\mbox{--}6.5$ (Giavalisco et al. 2004; Vanzella et al. 2009; Dahlen et al. 2010). However, the Lyman-break technique relies on detecting and measuring rest-frame UV light, and may therefore miss a fraction of dusty galaxies whose UV emission is strongly attenuated and reddened. To investigate how dusty galaxies are missed in LBG samples, we compare our GOODS-N FIR+mm detected sources at z = 2.8–4.4 and $z=4.4\mbox{--}5.5$ with the B- and V-dropout criteria (e.g., Dickinson et al. 2004; Giavalisco et al. 2004; Stark et al. 2009; Vanzella et al. 2009; Bouwens et al. 2012).

We use the HST optical to near-infrared photometry from Skelton et al. (2014; cross-matched in Section 2) and use the color–color criteria from Vanzella et al. (2009) to select B- and V-dropout LBGs in GOODS-N:

• 1.  
  B-dropout

  $$\begin{eqnarray}&&\left\{\begin{array}{l}\ ({B}_{435}-{V}_{606})\gt ({V}_{606}-{z}_{850})+1.1\\ \ \mathrm{and}\ ({B}_{435}-{V}_{606})\geqslant 1.1\\ \ \mathrm{and}\ ({V}_{606}-{z}_{850})\leqslant 1.6\end{array}\right.\end{eqnarray} \tag{ 5 }$$
• 2.  
  V-dropout

  $$\begin{eqnarray}&&\left\{\begin{array}{l}\left\{\begin{array}{l}\ ({V}_{606}-{i}_{775})\gt 0.9\times ({i}_{775}-{z}_{850})+1.5\\ \ \mathrm{or}\ ({V}_{606}-{i}_{775})\gt 2.0\end{array}\right.\\ \ \mathrm{and}\ ({V}_{606}-{i}_{775})\gt 1.2\\ \ \mathrm{and}\ ({i}_{775}-{z}_{850})\lt 1.3\\ \ \mathrm{and}\ {\rm{S}}/{\rm{N}}({B}_{435})\lt 2\\ \mathrm{and}\ {\rm{S}}/{\rm{N}}({i}_{775})\geqslant 3.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

In Figure 32, FIR+mm detected sources at any redshift that also have V₆₀₆ and z₈₅₀ three-sigma detections are plotted in the ${B}_{435}-{V}_{606}$ versus ${V}_{606}-{z}_{850}$ diagram. The subscripts indicate the HST observing bands (i.e., F435W for B band, F606W for V band, and F850LP for z band). We also show two modeled tracks as the colored curves, which are calculated by redshifting SED templates (BC03+DL07; e.g., Figure 2). Here we show two example model tracks with $E(B-V)=0.05$ and 0.1. Heavier dust attenuation will make the model track shift toward redder colors. Data points within the upper-left region identified by the dotted lines indicate galaxies that meet the B₄₃₅-dropout color criteria.

Zoom In Zoom Out Reset image size

Among all the FIR+mm and V₆₀₆ and z₈₅₀ detected galaxies, we highlight in Figure 32 the 32 ones having $2.8\leqslant z\leqslant 4.4$. Only 9 of them can be identified as B-dropouts (we use 1σ limits, as in Vanzella et al. 2009). In comparison, there are a total of 78 FIR+mm galaxies with photometric redshift within $2.8\,\leqslant z\leqslant 4.4$ (59 have valid optical photometry in 3D-HST catalog, often thus with near-IR WFC3 detections but lacking ACS detections in the optical). Therefore, only half of our FIR+mm sample at $2.8\lt z\lt 4.4$ has a significant HST optical counterpart in V and i to start with, and only about 10% of the parent sample (20% of those with optical counterparts) could be identified as B-dropouts at the depths of the current HST imaging. The optically detected FIR+mm galaxies with $2.8\lt z\lt 4.4$ fail to be selected as B-dropouts for a number of reasons: photometric noise affecting these faint detections, effect of photometric redshifts (some objects might be slight out of the boundary, in reality), lack of depth in the optical data (upper limits in the B band are not always stringent enough). In addition, even for the very small number of galaxies that can be selected as B-dropouts, the rest-frame UV-based SFR significantly underestimates the dust-obscured SFR. For example for ID 3276, a source that is detected in B band and is within the B-dropout region in Figure 32 (with its ID number labeled in the figure), the optical/near-IR (rest-frame UV) based SFR, corrected for dust attenuation, is $\sim 32\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (stellar mass ${10}^{11.06}\,{M}_{\odot }$) in the 3D-HST catalog. For comparison, its FIR+mm based SFR is $256\pm 37\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ in our FIR+mm SED fitting.

In Figure 33, we show the distribution of FIR+mm galaxies with i₇₇₅ and z₈₅₀ detections in the ${V}_{606}-{i}_{775}$ versus ${i}_{775}-{z}_{850}$ diagram. None of the seven galaxies in our sample with a photometric redshift within the V-dropout range of $4.4\,\leqslant z\leqslant 6.5$ is detected in i₇₇₅ and z₈₅₀. There are some galaxies within the V-dropout color selection region for which the photometric redshift estimate is $z\lt 4.4$.

Zoom In Zoom Out Reset image size

All in all, there is very little overlap between our FIR+mm-selected sample of high redshift galaxies and dropout samples. For those few FIR+mm-bright objects that can be selected as LBGs, their UV-estimated SFRs are largely underestimated compared to those derived from their far-infrared luminosities. We thus conclude that the SFR density contributions we have derived from IR data are largely independent from the ones estimated from UV-selected samples, at least at $z\gtrsim 3.5$.

In our FIR+mm catalog, we find three sources that have FIR+mm SED photometric redshifts $z\geqslant 5$, and whose SFRs and stellar masses appear to be consistent with the MS at these redshifts (at least within the current, large uncertainties on its locus). These are ID 14914 (which is not in the goodArea), ID 6406, and ID 12646 (the right-most blue data point in Figure 23, our highest redshift candidate, with $z\sim 6.8$). We discuss the properties of these three particularly interesting sources in some details below.

Zoom In Zoom Out Reset image size

ID 14914 is detected in all FIR+mm bands. Multi-wavelength SED and cutout images are presented in Figures 34 and 35. Weak extended emission can be seen in the HST ACS and WFC3 images. Two sources from the 3D-HST catalogs are very close to ID 14914. We cross-matched the closer (0585) 3D-HST source, to the north-east, as the counterpart. It is clearly detected by IRAC, as well as in the (sub)mm and at 20 cm. Blending is not significant. The 3D-HST ID number and other properties derived from our analysis are listed below:^{35 ,36}

$$\begin{eqnarray}{\mathrm{ID}}_{\mathrm{IRAC}} & = & 14914\\ {\mathrm{ID}}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 7338\\ {z}_{\mathrm{FIR}+\mathrm{mm}} & = & 5.33\pm 0.18\\ {z}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 5.57\ (5.335\mbox{--}5.90)\\ {M}_{* ,3{\rm{D}} \mbox{-} \mathrm{HST}} & = & {10}^{12.55}\ {M}_{\odot }\\ {M}_{* } & \approx & {10}^{11.68}\ {M}_{\odot }\\ {\mathrm{SFR}}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 436.5\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\\ {\mathrm{SFR}}_{\mathrm{FIR}+\mathrm{mm}} & = & 1724\pm 132.1\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\\ \langle U\rangle & \approx & 45.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

This source is also identified as AzGN04 and GN1200.12 in the 1.1–1.2 mm catalogs of Perera et al. (2008), Greve et al. (2008), and Penner et al. (2011). We show the photometry of Chapin et al. (2009) and Penner et al. (2011) in Figure 34 for comparison (as red circles and magenta triangle respectively). Our measurements are in good agreement with these previous values.

Note that the 3D-HST catalog (Skelton et al. 2014) reports a stellar mass of ${M}_{* ,3{\rm{D}} \mbox{-} \mathrm{HST}}={10}^{12.55}\,{M}_{\odot }$ for this object, at a photometric redshift z = 5.57. However, from our SED fitting, or simply by comparing IRAC fluxes with other, similar sources, such as GN20, we find a much lower stellar mass at a similar redshift, z = 5.33. (see footnote 36) This value is not far from our own SED fitting results, and it seems more reasonable than the 3D-HST value. Also note that this source is not in the goodArea; thus it was not included in our cosmic SFR density analyses.

ID 6406 is detected in the long-wavelength FIR and (sub)mm bands. Image cutouts are presented in Figure 36. It is almost invisible in ACS images, but is detected in the near-IR WFC3, and becomes clear in the IRAC images, as well as the (sub)mm and radio data. Its 3D-HST and FIR+mm properties are listed below:

$$\begin{eqnarray}{\mathrm{ID}}_{\mathrm{IRAC}} & = & 6406\\ {\mathrm{ID}}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 26337\\ {z}_{\mathrm{FIR}+\mathrm{mm}} & = & 5.24\pm 0.23\\ {z}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 5.47\ (3.725\mbox{--}5.82)\\ {M}_{* ,3{\rm{D}} \mbox{-} \mathrm{HST}} & = & {10}^{11.44}\ {M}_{\odot }\\ {\mathrm{SFR}}_{3{\rm{D}} \mbox{-} \mathrm{HST}} & = & 0.1\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\\ {\mathrm{SFR}}_{\mathrm{FIR}+\mathrm{mm}} & = & 918\pm 109\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\\ \langle U\rangle & \approx & 50.\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

This source has no millimeter counterpart in Perera et al. (2008), Greve et al. (2008), and Penner et al. (2011). The closest mm source is AzGN31 in the Perera et al. (2008) catalog, at a distance of $4\buildrel{\prime\prime}\over{.} 9$, and ID 33 in the Penner et al. (2011) catalog, at a distance of $8\buildrel{\prime\prime}\over{.} 8$. Recent NOEMA mm follow-up observations detect a source that accurately corresponds with the IRAC and VLA 20 cm position (D. Liu et al. 2018, in preparation).

Multi-wavelength cutouts of ID 12646are presented in Figure 37. There is no obvious detection in ACS or WFC3 images at the exact IRAC position. The closest 3D-HST source is 3D-HST ID 12669 at a photometric redshift z = 4.65 (4.248–5.06), but this is located $1\buildrel{\prime\prime}\over{.} 52$ away from the IRAC position. Given the radio localization, this is unlikely to be the counterpart of ID 12646. ID 12646 is also not matched to any source in the catalog of Pannella et al. (2015); hence we do not have a stellar mass for it. The emission starts to be seen in the WIRCAM K_s image, but has an irregular morphology and low ${\rm{S}}/{\rm{N}}$ ratio. It is fairly bright in IRAC images and bright at (sub)mm wavelengths, but with some offset from the emission peak. It matches AzGN10 in Perera et al. (2008), Chapin et al. (2009), and Penner et al. (2011). Its physical properties are listed below:³⁷

$$\begin{eqnarray}{\mathrm{ID}}_{\mathrm{IRAC}} & = & 12646\\ {z}_{\mathrm{FIR}+\mathrm{mm}} & = & 6.8\pm 1.0\\ {M}_{* } & \approx & {10}^{11.44}\ {M}_{\odot }\\ {\mathrm{SFR}}_{\mathrm{FIR}+\mathrm{mm}} & = & 814\pm 249\ {M}_{\odot }\,{\mathrm{yr}}^{-1}\\ \langle U\rangle & \approx & 31.\end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

The photometry of Chapin et al. (2009) and Penner et al. (2011) is also shown in the ID 12646 panel of Figure 34. The measurements generally agree with each other. To estimate the photometric redshift, Chapin et al. (2009) adopted a simple method which infers redshift as a function of IRAC and 24 μm flux according to Pope et al. (2006), and derived a photometric redshift $z\sim 3$. Note that using their equation, GN20 (ID 564 in this work) with a spectroscopic redshift of 4.055 would have a photometric redshift $z\sim 2.6$.

While these sources are interesting candidates for the most distant IR-detected galaxies in GOODS-N, we emphasize that we still only have photometric redshift estimates. In particular, the optical/near-IR SEDs are very red, faint at optical wavelengths, and rising at >3 μm. It is hard to obtain precise redshifts in these cases. However, note that the redshifts estimated from the FIR+mm SEDs agree well with the optical/near-IR photometric redshifts in all cases, while using templates with $\langle U\rangle \sim 30\mbox{--}40$, lower than what is expected for MS galaxies at $z\gt 4$ (Béthermin et al. 2015). Because of the dust temperature versus redshift degeneracies, using warmer templates would further increase the estimated FIR+mm redshifts. We cannot exclude that these sources are actually SB galaxies, and their inferred $\langle U\rangle$ values would be appropriate in that case. However, even in that case, substantially lower redshifts would require very cold dust SEDs, which seems unlikely.

Spectroscopic observations will be essential to measure redshifts for these objects, and to understand the reasons behind the different SED fitting results. Due to the dusty nature of these FIR+mm galaxies, optical/near-IR spectroscopic observations will likely not be as efficient as for LBGs. Mid-IR spectroscopic surveys with JWST, and (sub)mm spectral line scans with ground-based interferometric arrays, will be the most efficient ways to measure redshifts.