RESOLVING THE COSMIC FAR-INFRARED BACKGROUND AT 450 AND 850 μm WITH SCUBA-2

Chian-Chou Chen; Lennox L. Cowie; Amy J. Barger; Caitlin. M. Casey; Nicholas Lee; David B. Sanders; Wei-Hao Wang; Jonathan P. Williams

doi:10.1088/0004-637X/776/2/131

1. INTRODUCTION

Following the discovery of the far-infrared (FIR) extragalactic background light (EBL) by the COBE satellite (Puget et al. 1996; Fixsen et al. 1998), many surveys have been conducted to detect the sources producing this light. Such studies have used both ground-based telescopes (e.g., Smail et al. 1997; Barger et al. 1998; Hughes et al. 1998) and space-based satellite missions (e.g., Oliver et al. 2010; Berta et al. 2011). Given that there is a comparable amount of light absorbed by dust and re-radiated in the FIR as there is seen directly in the UV/optical (Dole et al. 2006), the dusty sources uncovered by these surveys are key in the development of a full understanding of galaxy formation. The FIR number counts provide simple yet fundamental constraints on empirical models (e.g., Valiante et al. 2009; Béthermin et al. 2011) and semi-analytical simulations (Hayward et al. 2013a, 2013b).

The construction of the FIR number counts began with 850/450 μm observations made using the SCUBA camera (Holland et al. 1999) mounted on the 15 m James Clerk Maxwell Telescope (JCMT; Smail et al. 1997; Barger et al. 1998, 1999; Hughes et al. 1998; Eales et al. 1999, 2000; Cowie et al. 2002; Scott et al. 2002; Smail et al. 2002; Borys et al. 2003; Serjeant et al. 2003; Webb et al. 2003; Wang et al. 2004; Coppin et al. 2006; Knudsen et al. 2008; Zemcov et al. 2010). Since then, many single-dish telescopes, instruments, and missions have been developed to survey the sky at FIR through millimeter wavelengths. At λ < 500 μm, the number counts have been established by the second-generation Submillimeter High Angular Resolution Camera (Dowell et al. 2003) at the Caltech Submillimeter Observatory (CSO; e.g., Khan et al. 2007), the Balloon-borne Large Aperture Submillimeter Telescope (Pascale et al. 2008; e.g., Devlin et al. 2009), the Herschel Space Observatory (hereafter Herschel; Pilbratt et al. 2010; e.g., Oliver et al. 2010; Berta et al. 2011). At λ > 500 μm, in addition to SCUBA, the number counts have been probed by the LABOCA camera (Siringo et al. 2009) on the Atacama Pathfinder Experiment (Güsten et al. 2006; e.g., Weiß et al. 2009), the AzTEC (Wilson et al. 2008) camera on both the JCMT (e.g., Perera et al. 2008; Austermann et al. 2010) and the Atacama Submillimeter Telescope Experiment (Ezawa et al. 2004; e.g., Aretxaga et al. 2011; Scott et al. 2010, 2012), the Max-Planck Bolometer array (Kreysa et al. 1998) on the IRAM 30 m telescope (e.g., Greve et al. 2004; Bertoldi et al. 2007), and Bolocam (Glenn et al. 1998) on the CSO (e.g., Laurent et al. 2005).

The biggest challenge for constructing the number counts in the FIR is poor spatial resolution (typically > 10''), due to the diffraction limits of the single-dish telescopes at longer wavelengths. Poor resolution has imposed a fundamental limitation, the confusion limit (Condon 1974), on our ability to resolve directly the faint sources that contribute the bulk of the background light. Poor resolution also prevents us from resolving close pairs within the large beam sizes. Interferometric observations (Wang et al. 2011; Barger et al. 2012; Smolčić et al. 2012; Karim et al. 2013; Hodge et al. 2013) and semi-analytical models (Hayward et al. 2013a, 2013b) have shown these to be common, and their effects must be understood in order to construct the true counts.

Techniques have been developed to work around the problem of the confusion limit. Surveys targeting massive galaxy clusters have unveiled a few faint sources with fluxes many times below the confusion limit through gravitational magnification (Smail et al. 1997, 2002; Cowie et al. 2002; Knudsen et al. 2008; Johansson et al. 2011; Chen et al. 2013), though the positional uncertainties can still cause large uncertainties in the amplifications and in the intrinsic source fluxes (Chen et al. 2011). In blank-field surveys, probability of deflection analyses, or P(D), using the number distribution of pixel values put stringent constraints on counts deeper than the confusion limit (Scheuer 1957; Weiß et al. 2009; Scott et al. 2010; Glenn et al. 2010; Lindner et al. 2011).

The new SCUBA-2 camera (Holland et al. 2013) mounted on the JCMT provides the fastest mapping capability at 450 and 850 μm with the best FIR spatial resolution (FWHM ∼ 7 farcs 5) at 450 μm among single-dish FIR telescopes. This greatly enhances our ability to resolve the 450 μm EBL, thanks to the smaller confusion limit, and makes the 450 μm number counts less affected by close pairs. In addition, the better positions provide better determinations of the lensing amplifications.

We have used SCUBA-2 to target two well-studied massive lensing clusters, A370 and A1689, and three blank fields, COSMOS, CDF-N, and CDF-S, in order to construct the 450 μm number counts over the widest possible flux range. We have previously shown some of these results on A370 in Chen et al. (2013) and on COSMOS in Casey et al. (2013). In this paper, we present the full results from all five fields, and we combine the data to measure the number counts. The details of the observations and the data reduction are presented in Section 2. In Section 3, we explain our methodology for constructing the number counts, which uses a combination of gravitational lensing and P(D) analysis, and present our results. We discuss the effects of field-to-field variance and source blending (multiplicity) in Section 4. We discuss the implications of our results in Section 5, and we provide a brief summary in Section 6.

2. OBSERVATIONS AND DATA REDUCTION

The SCUBA-2 data were taken between late 2011 and early 2013. The observations on A1689, A370, and COSMOS were taken under the best weather conditions (band 1, τ_225GHz < 0.05). Most of the observations carried out on CDF-N and CDF-S were under band 2 (0.05 < τ_225GHz < 0.08).⁶ To cover a wide range of submillimeter fluxes, we used the CV Daisy scan pattern on both the massive cluster fields, where the smaller field is well matched to the strong lensing regions in the clusters, and the blank fields (except COSMOS). In order to cover larger areas uniformly to find brighter but rarer sources, we also used the PONG-900 scan pattern on the blank fields. Detailed information about the SCUBA-2 scan patterns can be found in Holland et al. (2013). We summarize the details of our observations in Table 1.

Table 1. SCUBA-2 Observations

Field	Centroid Coordinate		Weather	Scan Mode	Total Exposure	Effective Area^a	$\bar{\sigma }$ ^b
	R.A.(J2000)	Decl.(J2000)			Total Exposure	[450 μm, 850 μm]	[450 μm, 850 μm]
	(H:M:S)	(D:M:S)			(hr)	(arcmin²)	(mJy Beam⁻¹)
A1689	13:11:29.8	−01:20:35.8	1	CV Daisy	14.8	[69.4, 73.0]	[ 4.5, 0.79]
CDF-N	12:36:49.6	+62:13:53.0	1+2	CV Daisy + PONG-900	24.7	[102.6, 104.7]	[ 8.8, 0.83]
CDF-S	03:32:28.0	−27:48:30.0	2	CV Daisy + PONG-900	22.9	[97.0, 102.5]	[10.3, 0.88]
A370	02:39:53.0	−01:34:38.0	1	CV Daisy	13.7	[83.1, 83.9]	[ 5.3, 1.07]
COSMOS	10:00:24.0	+02:24:00.0	1	PONG-900	38.0	[286.0, 281.7]	[ 4.7, 0.87]

Notes. ^aTotal area to 1.5–2 times the central noise level. In cluster fields they represent the total source plane area. ^bAverage 1σ sensitivity within the effective area. Note that the quoted values for A370 and COSMOS are slightly higher than those quoted by Chen et al. (2013) and Casey et al. (2013), as they were quoting maximum sensitivities.

Download table as: ASCII Typeset image

We reduced the data using the Dynamic Iterative Map Maker (DIMM) in the SMURF package from the STARLINK software developed by the Joint Astronomy Centre (Jenness et al. 2011; Chapin et al. 2013) released after 2012 July 1. DIMM performs iterative estimations on the common mode signal, the astronomical signal, and the white noise. It also does flatfield and extinction corrections and applies a Fourier Transform filter to remove low-frequency excess signal relative to the white noise that is not able to be removed through common mode subtraction (Chapin et al. 2013). We adopted the standard configuration file dimmconfig_blank_field.lis for our science purposes. We ran DIMM on each bolometer subarray individually to avoid data splitting, and we used the MOSAIC_JCMT_IMAGES recipe in PICARD, the Pipeline for Combining and Analyzing Reduced Data (Jenness et al. 2008), to coadd the products into final maps.

To increase the detectability of point sources, as nearly all the submillimeter galaxies (SMGs) are expected to be much smaller than the ∼10'' resolution, we applied a matched-filter to our maps. This provides a maximum likelihood estimate of the source strength (e.g., Serjeant et al. 2003). Assuming S(i, j) and σ(i, j) are the signal and rms noise maps produced by DIMM, and PSF(i, j) is the signal point-spread function (PSF), the filtered signal map F(i, j) would be

$\begin{equation} F(i,j) = \frac{\sum _{i,j} [S(i,j) / \sigma (i,j)^2 \times {\rm PSF}(i,j)]}{\sum _{i,j} [1 / \sigma (i,j)^2 \times {\rm PSF}(i,j)^2]}, \end{equation} \tag{ 1 }$

and the filtered noise map N(i, j) would be

$\begin{equation} N(i,j) = \frac{1}{\sqrt{\sum _{i,j} [1 / \sigma (i,j)^2 \times {\rm PSF}(i,j)^2]}}. \end{equation} \tag{ 2 }$

Ideally, the PSF for the matched-filter algorithm is a Gaussian normalized to a peak of unity with FWHM equal to the JCMT beam size at a given wavelength (i.e., 7 farcs 5 at 450 μm and 14'' at 850 μm). However, the map produced from DIMM usually has low spatial frequency structures that need to be subtracted off before performing the source extraction. Thus, before running the matched-filter, we convolved the map with a broad Gaussian normalized to a sum of unity, and we subtracted this convolved map from the original map. Note that in Chen et al. (2013), we showed that the source fluxes and the signal-to-noise ratio (S/N) are not sensitive to the size of the FWHM for reasonable choices. Thus, we simply adopted the default values (20'' at 450 μm and 30'' at 850 μm). To optimize the S/N, we processed the PSF used for the matched-filter similarly. It becomes a Gaussian with a convolved broader Gaussian subtracted off, which gives a Mexican-hat-like wavelet. We adopted the PICARD recipe SCUBA2_MATCHED_FILTER for the tasks described above.

We then calibrated the fluxes using standard flux conversion factors (FCFs; 491 Jy pW⁻¹ for 450 μm and 537 Jy pW⁻¹ for 850 μm). The relative calibration accuracy is shown to be stable and good to 10% at 450 μm and 5% at 850 μm (Dempsey et al. 2013). Ten percent upward corrections were applied to compensate for the flux lost during filtering, which we estimated from simulations using fake sources. We also tested the accuracy of the flux calibrations by—instead of using the standard FCFs—adopting the FCFs obtained from the calibrators in each night of observations, as was done in Chen et al. (2013). We found the results agree to better than 10%, which is essentially the uncertainty of the FCFs.

3. NUMBER COUNTS

In order to measure the galaxy number counts, we need source-free maps with only pure noise to estimate how many fake sources are contaminating the counts. We followed the procedure in Chen et al. (2013). For each wavelength, we generated two data maps, each with roughly half of the total exposure time, and we subtracted them to obtain the source-free maps. We rescaled the value of each pixel following the equation $\sqrt{t1\times t2}$ /(t1+t2), with t1 and t2 representing the exposure time of each pixel from the two maps. Finally, we applied the matched-filter and FCFs, as we did on the signal maps. We produced these true noise maps, sometimes referred to as jackknife maps in the literature, for each of our fields. We show the S/N histograms of the true noise maps (green curves) and the signal maps (gray shading) for the COSMOS field in Figure 1. The black curves are the expected pure noise distributions with σ = 1, and they agree nicely with the results of our true noise maps. The positive long tails, as well as excess signals relative to pure noise, are from real astronomical sources. Because of the negative trough of the matched-filter PSF, we also see a negative tail in the distribution (Chapin et al. 2013).

**Figure 1.** 450 μm and 850 μm S/N histograms for the pixels located within the regions of the COSMOS signal map where the noise level is less than 1.5 times the central noise level (gray shading) and for the pixels in the corresponding regions of the true noise map (green). The black curves are the expected pure noise distributions with σ = 1. Dashed vertical lines are 4 σ cuts.
Download figure:
Standard image High-resolution image

3.1. Methodology

In previous work, we extracted sources down to ∼4 σ and used these robust catalogs with low contamination rates (⩽5%; e.g., Chen et al. 2013; Casey et al. 2013) to calculate the number counts. However, as shown in Figure 1, excess positive signal can be seen to ∼2 σ. For studies on individual sources, low contamination rates are essential to justify the robustness of the sample. However, for number counts analyses where positional information is no longer important, we can lower our detection threshold to a level where there are still excess counts that are statistically significant. We therefore adjusted our S/N thresholds to ∼2 σ. Because the exact thresholds vary for each field, we experimented with different binning and S/N cutoffs in order to exploit fully the statistically significant signals to obtain the deepest counts possible.

Following Chen et al. (2013), we generated our source catalogs by identifying the peak S/N pixel, subtracting this peak pixel and its surrounding areas using the PSF scaled and centered on the value and position of that pixel, and then searching for the next S/N peak. We iterated this process until we hit the S/N threshold.

We generated the PSFs used for creating the catalogs by making a weighted average of all primary calibrators taken before and after the science data. These are mostly Uranus, CRL618, and CRL2688. As an example, in Figure 2 we show how the normalized averaged PSFs of A1689 agree well with the PSFs of the individual strong sources detected in A1689. The FWHM of the PSF is slightly larger than the ideal size, which could mean the observations were slightly out of focus.

**Figure 2.** Normalized PSFs as a function of the distance, expressed in arcseconds relative to the center. Averaged PSFs for the A1689 maps are shown in red, and PSFs for the individual strong sources (S/N > 10 at 850 μm and S/N > 7 at 450 μm) detected in the A1689 maps are shown with dashed curves.
Download figure:
Standard image High-resolution image

We ran the extraction on both the signal maps and the true noise maps. We computed the number density for each extracted source by inverting the detectable area, which is the area over which the source can be detected above the S/N threshold given the noise level. We then calculated the number counts by summing up the number densities of the sources selected in each flux bin. Finally, we subtracted the counts obtained from the true noise maps, if any, from the counts obtained from the signal maps, to produce the pure source number counts.

We plot the differential number counts of all five fields in Figure 3. The black (blue) symbols are the counts from the signal (true noise) maps, and the red symbols are the pure source counts. The vertical dashed lines represent the mean 4σ depth of each field. The counts are dominated by noise at the faint end and by real sources at the bright end. The last few bins start to show signs of incompleteness. We de-lensed all the counts in the cluster fields using LENSTOOL (Kneib et al. 1996) and the latest lensing models from Richard et al. (2010) (A370) and Limousin et al. (2007) (A1689). Thanks to the gravitational lensing in the cluster fields, we can see statistically significant pure source counts much deeper than can be probed with the 4σ limits. Our statistical analysis also allows us to detect deeper counts in the blank fields.

We note that fair estimations of the pure noise counts are critical to compute the legitimate pure source counts. In SCUBA-2 maps, the noise in neighboring pixels is correlated due to the fact that as the bolometer array scans through the sky, the pixels covered by the same bolometer record the noise pattern from that bolometer. Since the noise patterns are correlated in the time domain, the noise of the pixels close to one another are correlated. Also, the performance of each subarray is different, which leads to some parts of the map being nosier than others (Holland et al. 2013), and this property of the noise is retained in the true noise maps. Thus, the pure noise counts obtained from the true noise maps should be the most representative noise counts.

However, in non-uniform maps like ours, pure noise counts based on the positions of single extractions from the true noise maps could be biased, given that occasionally some signals would happen to be located in small noisier regions, or in highly amplified regions in lensing fields. We tested for possible biased noise counts by iterating the following process: we took the source catalogs generated from the true noise maps and randomized the positions of each source. In the catalogs we have the source information on S/N, fluxes, rms values, and positions. In the process of randomization, we retained the S/N information. With the new assigned position, each source has a new rms value obtained by matching the new position to the original rms value map. In Figure 4, we show an example rms value map for A1689 at 850 μm. The rms value maps were generated by computing the variance of the data that lands in each pixel. Thus the rms value maps naturally keep the information of the performance of the bolometers, as well as the correlations among the nearby pixels. Note also that we prevented two random positions from coming closer than half of the beam FWHM, which is the case in the real catalogs, since we removed the detected signals by subtracting a PSF scaled to the peak of the detection during our process of signal extraction. Once the new rms values were assigned, the fluxes were calculated based on the new rms values and the retained S/Ns. We then computed the noise counts according to the new fluxes and positions of all the sources.

**Figure 4.** rms value map for A1689 at 850 μm in units of mJy/beam. The blue dashed contour shows two times the central noise.
Download figure:
Standard image High-resolution image

We iterated this process 50 times and calculated the average noise counts for each flux bin. We subtracted the average counts from the simulations from the counts from the true noise maps and plotted them in Figure 5. The noise counts obtained from these two methods agree with each other to within the uncertainties for all five fields. Thus, we conclude that the pure noise counts obtained from the true noise maps are robust.

**Figure 5.** Differenced differential counts per bin (i.e., average counts from the simulations minus the counts from the true noise maps) on A1689 (orange triangles), CDF-N (red stars), CDF-S (brown upside down triangles), A370 (green squares), and COSMOS (blue circles) with 1σ error bars. There is no statistical difference between these two sets of counts.
Download figure:
Standard image High-resolution image

Following Chen et al. (2013), we ran Monte Carlo simulations to estimate the underlying counts models. We first randomly populated the true noise map with simulated sources, drawn from an assumed model and convolved with the PSFs, to form a simulated image. The counts model is in the form of a broken power law

$\begin{eqnarray} \frac{dN}{dS} = \left\lbrace \begin{array}{@{}ll}N_0\left(\frac{S}{S_0}\right)^{-\alpha } & \quad {\rm if}\; S \le S_0\\ {N_0}\left(\frac{S}{S_0}\right)^{-\beta } & \quad {\rm if} S {>} S_0 \\ \end{array} \right.. \end{eqnarray} \tag{ 3 }$

The faintest fluxes we adopted for any of our models are the fluxes at which the integrated flux density agrees with the EBL measurements within errors (Puget et al. 1996; Fixsen et al. 1998). We then extracted the signal and computed the recovered number counts in exactly the same way as we did with the real data maps. We measured the ratio between the recovered counts and the input counts, which reflects the Eddington bias (Eddington 1913), and then applied this ratio to the statistically significant observed counts to correct for that bias. We did a χ² fit to the corrected observed counts using a broken power law to obtain the normalization and power law indices. We used this fit as the next iteration of the model counts in the procedure and repeated the process. We continued until the input model agreed with the corrected counts at the 1σ level throughout the statistically significant range.

There are only three (four) statistically significant points in A370 450 μm (CDF-S 850 μm) counts, so we only fitted the normalization to avoid overfitting. We also note that the fitting results were not affected by specifying the break positions of the broken power laws into the fit. We chose to fix the breaks to again avoid overfitting and to obtain better statistical constraints on the fit. We excluded the CDF-N and CDF-S from the 450 μm analysis, because we have no statistically significant 450 μm counts for these fields.

For the cluster fields, we populated the simulated sources in the source plane and imaged them onto the image plane using LENSTOOL. At 850 μm, we located the source planes at z = 3.0 based on the latest observational results (Barger et al. 2012; Vieira et al. 2013) and theoretical models (Hayward et al. 2013b). At 450 μm, two recent SCUBA-2 results on the COSMOS field have shown that the majority of the 450 μm sources are at z < 3 (Geach et al. 2013; Casey et al. 2013). We located the 450 μm source planes at z = 1.3 based on Geach et al. (2013), since their flux range is closer that of our observations. In any case, we stress that our statistical approach to estimating the counts is not sensitive to the adopted source plane redshifts (Blain et al. 1999).

On the other hand, however, the de-lensing process could be significantly affected by the uncertainties of the source positions (Chen et al. 2011). We tested this bias by again running the Monte Carlo simulations on the extracted signals by randomizing their positions. We showed in Chen et al. (2013) that the positional uncertainty is a function of S/N; thus, we randomized the position of each source by setting the offsets smaller than the 90% confidence boundary given its S/N. We then calculated the de-lensed counts based on the new positions. We iterated this process 50 times and obtained the average source counts. We then subtracted the average counts from the counts obtained from the real signal maps. We show our results in Figure 6. The counts obtained from these two methods agree with each other to within the uncertainties for both A1689 and A370. This illustrates the robustness of our methodology. Given the small effects of the positional uncertainties on the de-lensing process, we conclude that using our statistical method, the uncertainty on the overall counts caused by the de-lensing process is negligible.

**Figure 6.** Differenced differential counts per bin (i.e., average counts from the simulations minus the counts from the real signal maps) for A1689 (orange triangles) and A370 (green squares) with 1σ error bars. There is no statistical difference between these two sets of counts.
Download figure:
Standard image High-resolution image

3.2. Results

We show the corrected differential number counts for each field individually in Figure 7 (colored symbols). This figure clearly illustrates the advantage of observing both blank and cluster fields in order to probe a wide flux range. The black solid curves show the input models. The dashed curves show the minimum χ² fits. The gray shading denotes the 1σ error regions based on the fits. We also show the SCUBA-2 counts constructed using a 4σ detection threshold for A370 (Chen et al. 2013) and a 3.8σ detection threshold for COSMOS (Casey et al. 2013; Geach et al. 2013; gray symbols), in which they both extrapolated their counts to fainter end by doing source flux deboosting. These agree nicely with our results. In Tables 2 and 3, we summarize, respectively, the model parameters and the corrected counts for each field.

Table 2. The Resulting Parameters of the Broken Power Law Model Curves

Field	Wavelengths	N₀	S₀	α	β
Field	(μm)	(mJy⁻¹ deg⁻²)	(μm)	α	β
A1689	450	30	20.9	2.5	6.0
	850	160	5.45	2.25	5.0

CDF-N	850	340	4.5	1.4	3.5

CDF-S	850	270	4.5	2.1	3.0

A370	450	30	20.9	1.9	4.0
	850	200	5.45	2.15	3.0

COSMOS	450	25	20.9	2.5	6.0
	850	160	5.45	2.25	3.5

Download table as: ASCII Typeset image

Table 3. Corrected Differential Number Counts on Each Individual Field

A1689				CDF-N		CDF-S
S₄₅₀	dN/dS	S₈₅₀	dN/dS	S₈₅₀	dN/dS	S₈₅₀	dN/dS
(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)
1.26	$99640^{+67202}_{-60525}$	0.10	$5436690^{+4016800}_{-2727500}$	2.64	$609.0^{+428.4}_{-428.4}$	2.64	$813.4^{+212.2}_{-212.2}$
4.24	$884.4^{+469.3}_{-469.3}$	0.23	$133324^{+103520}_{-72100}$	4.05	$410.2^{+95.61}_{-87.84}$	4.05	$270.2^{+91.11}_{-81.75}$
11.23	$171.7^{+99.21}_{-99.21}$	0.53	$36446^{+26232}_{-23530}$	6.20	$96.39^{+44.02}_{-31.52}$	6.20	$129.0^{+44.52}_{-34.11}$
16.31	$48.85^{+30.42}_{-30.42}$	1.05	$4729^{+1698}_{-1698}$	9.51	$29.04^{+19.64}_{-12.55}$	9.51	$23.03^{+22.40}_{-12.54}$
22.65	$27.01^{+21.38}_{-13.75}$	1.72	$2649^{+745.5}_{-745.5}$	13.94	$5.21^{+11.97}_{-4.31}$	...	...
30.74	$5.41^{+12.43}_{-4.47}$	2.64	$886.0^{+219.1}_{-219.1}$	...	...	...	...
...	...	4.05	$313.4^{+68.39}_{-68.39}$	...	...	...	...
...	...	6.20	$127.7^{+76.30}_{-50.67}$	...	...	...	...
...	...	9.51	$8.02^{+18.45}_{-6.63}$	...	...	...	...

A370				COSMOS
S₄₅₀	dN/dS	S₈₅₀	dN/dS	S₄₅₀	dN/dS	S₈₅₀	dN/dS
(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)

5.20	$573.1^{+561.1}_{-561.1}$	0.35	$145886^{+203320}_{-142801}$	11.23	$55.85^{+52.68}_{-52.68}$	2.64	$705.1^{+134.0}_{-134.0}$
16.97	$106.7^{+83.48}_{-83.48}$	1.31	$8384^{+5740}_{-5740}$	16.31	$55.01^{+18.56}_{-18.56}$	4.05	$248.4^{+41.82}_{-39.62}$
35.78	$2.03^{+4.67}_{-1.68}$	2.64	$1485^{+434.4}_{-434.4}$	22.65	$11.55^{+4.55}_{-4.55}$	6.20	$118.4^{+22.78}_{-22.78}$
...	...	4.05	$386.7^{+122.6}_{-113.6}$	30.74	$2.18^{+2.12}_{-1.19}$	9.51	$10.97^{+7.42}_{-4.74}$
...	...	6.20	$239.7^{+73.32}_{-57.53}$	39.69	$0.47^{+1.09}_{-0.39}$	13.94	$9.84^{+7.78}_{-4.71}$
...	...	9.51	$37.94^{+25.66}_{-16.39}$	...	...	20.31	$1.50^{+3.46}_{-1.24}$

Download table as: ASCII Typeset image

In Figure 8, we show all the counts together for the two wavelengths. The black solid curves represent the best fit models, which we present in Table 4. At 850 μm, our results are almost indistinguishable from the SCUBA results in Knudsen et al. (2008) (dot–dot–dot–dashed curve) covering a similar flux range, and they agree within the errors with the Zemcov et al. (2010) results (long–dashed), which come from an analysis of all the SCUBA data taken on cluster fields. We do not observe a significant under-abundance in the CDF-S 850 μm counts above 3 mJy, as was found in the LABOCA wider area but shallower sensitivity survey of the Extended CDF-S (ECDF-S; Weiß et al. 2009, dot-dashed). We discuss potential causes for this discrepancy in Section 4.1. The SCUBA-2 counts for A370 (Chen et al. 2013; short-dashed) are slightly higher relative to the other fields at ∼3–6 mJy (850 μm) and ∼10–25 mJy (450 μm), which could be caused by sample variance from fields which are small in size compared to the large-scale structure. We show the hybrid Herschel 500 μm counts (Béthermin et al. 2012b) obtained from a mix of directly resolved counts above 20 mJy and 24 μm stacking counts below 20 mJy with asterisks. The noticeable overabundance of those counts relative to ours could indicate that the Herschel counts are biased upward due to source blending caused by poorer resolution (∼35'' beam FWHM compared to 7 farcs 5 beam FWHM). We also show various other results at 450 μm from the literature, and they all agree well with our results.

Table 4. Best χ² Fits on the Number Counts from All Five Fields at 450 and 850 μm

Wavelengths	N₀	S₀	α	β
(μm)	(mJy⁻¹ deg⁻²)	(mJy)	α	β
450	$24^{+30}_{-12}$	20.4	$2.53^{+0.67}_{-0.67}$	$5.57^{+4.5}_{-4.0}$
850	$120^{+65}_{-45}$	6.21	$2.27^{+0.5}_{-0.5}$	$3.71^{+2.5}_{-2.0}$

Download table as: ASCII Typeset image

4. POSSIBLE ISSUES

4.1. Field-to-field Variance: Is the CDF-S Underdense?

Several studies of the CDF-S have shown that at redshifts between 2 and 3, the massive red galaxies selected through rest-frame optical colors (DRGs and pBzKs) that constitute the bulk of the mass during that epoch are under-abundant relative to the mean density from other deep fields (van Dokkum et al. 2006; Marchesini et al. 2007; Blanc et al. 2008). However, in the same studies, there is no sign of an underdensity of non-DRGs (Marchesini et al. 2007) and less massive star-forming BzKs (Blanc et al. 2008). Together with the fact that recent 4 Ms Chandra observations also show no sign of an underdensity (Lehmer et al. 2012) at the lower flux end where high-redshift late-type galaxies start to dominate the X-ray number counts, this could imply that at 2 < z < 3, massive passive galaxies are underdense, while less massive star-forming galaxies are not.

On the other hand, an underdensity of the 870 μm sources was reported by the LABOCA survey of the ECDF-S (Weiß et al. 2009). However, we do not confirm this result. As shown in Figure 8, our SCUBA-2 observations at 850 μm of the central region of the LABOCA field show no sign of an underdensity. Recently, Scott et al. (2010) also found no apparent underdensity at 1.1 mm from an AzTEC survey toward a similar region as our coverage. Contrary to the AzTEC observations, our observed waveband is close to that of LABOCA, so both should see a similar population. Although the coverage of the LABOCA observations is wider (30' × 30') with slightly shallower but uniform sensitivity (∼1.2 mJy beam⁻¹), an underdensity in the central region of the LABOCA map should be apparent in our data.

In Figure 9, we show our 850 μm S/N map with the LABOCA detections plotted with blue circles and our 4σ detections with brown squares. While our SCUBA-2 observations recover all the LABOCA sources, many LABOCA detectable SCUBA-2 sources (sources with fluxes greater than the LABOCA limit; green circles) are missed by LABOCA. There are 19 SCUBA-2 sources with fluxes greater than 4.4 mJy (the LABOCA limit), only 8 of which are detected by LABOCA. If we correct for this factor of 2.375 in the cumulative counts, then the LABOCA counts would be in good agreement with the measurements from other fields.

To investigate the possible causes of this discrepancy, we first compared the fluxes of the sources detected in both observations. We found that both measurements agree very well, as shown in the upper panel of Figure 10. Note that we excluded two close LABOCA sources from the analyses described in this paragraph (two overlapped blue circles in Figure 9), because they are blended in the LABOCA maps, and thus the uncertainties on their flux measurements are high. We then examined the histogram of the expected S/N on the SCUBA-2 detected sources with LABOCA detectable fluxes, based on the SCUBA-2 fluxes and the LABOCA noise claim (green hatched regions in Figure 10). We found that while LABOCA succeeded in detecting all of the S/N > 5 sources (blue hatched region), LABOCA failed to detect all but one of the lower S/N sources. Ideally, if the two observations had agreed, then the two distributions would be identical.

**Figure 10.** Upper: the ratio of LABOCA 870 μm to SCUBA-2 850 μm fluxes on the sources that are detected in both observations. Lower: the green histogram shows the expected S/N for the SCUBA-2 detected sources with LABOCA detectable fluxes, based on the SCUBA-2 fluxes and the LABOCA noise claim. The blue histogram shows the same information for the fraction that are LABOCA detected.
Download figure:
Standard image High-resolution image

On the other hand, the extra detections in our maps are all at very high S/N (> 5), and many of them are detected in nightly maps. We cross-correlated them with the sources detected at other wavelengths. Of our 12 SCUBA-2 sources that are not also detected in the LABOCA survey, 10 have either 24 μmor radio counterparts, and 6 are detected at both 24 μm and radio. Moreover, six are detected by AzTEC at 1.1 mm. The fact that our extra detections are highly significant and highly correlated with sources detected at other wavebands makes it extremely unlikely that many of them are spurious.

4.2. The Effects of Multiplicity

Recently, many high-resolution (∼1'' beam FWHM) submillimeter and millimeter interferometric observations have shown that a significant percentage (20%–40%) of single-dish detected submillimeter sources are in fact composed of two or sometimes three separated sources (e.g., Wang et al. 2011; Smolčić et al. 2012; Barger et al. 2012; Hodge et al. 2013). In semi-analytical simulations it has been shown that multiplicity could dramatically impact the number counts at 850 μm obtained from single-dish observations (Hayward et al. 2013b). However, the observational constraints on the fraction of multiples are subject to small number statistics and heterogeneity in the sample selection. Above an 850 μm flux of 7 mJy, three of the eight Submillimeter Array (SMA) observed SCUBA sources in the CDF-N were found to be multiples by Barger et al. (2012), which corresponds to a multiple fraction of 37.5% with a ±1 σ range from 17%–74%. Smolčić et al. (2012) compiled a list of millimeter and submillimeter interferometric continuum follow-up observations of a sample of 36 LABOCA sources at 870 μm in the COSMOS field, and they reported that 6 of their 27 interferometric detections are multiples (22% ± 9%). With a larger sample size, Hodge et al. (2013) reported that 24 of their 69 ALMA robustly detected and LABOCA pre-selected sources (MAIN ALESS sample) are multiples, giving a multiple fraction of ∼35% ± 7%. While the FWHM of the LABOCA beam is about 37% larger than that of SCUBA-2 (19 farcs 2 versus 14''), it is unclear how much this difference in spatial resolution affects the multiple fraction. If the spatial distribution among close multiples were random, as shown by Hodge et al. (2013), then the multiple fraction for the SCUBA-2 850 μm selected sources would be ∼19% based on scaling the ALMA multiplicity of the MAIN ALESS sample.

We have obtained SMA observations (all are detections) of eight new SCUBA-2 850 μm sources in the CDF-N with arcsecond spatial resolution. Together with previous SMA observations of SCUBA detected sources in the same field (mostly from Barger et al. 2012), we have compiled a list of 24 SCUBA-2 detected (4 σ) and SMA observed sources in the CDF-N. A detailed analysis is given in A. Barger et al. (2013, in preparation). All of our new SMA observations are single source detections, and only 3 out of 24 SMA detected sources with S₈₅₀ > 3.5 mJy break into close pairs; that is, a multiple fraction of 12.5 $^{ +12.1}_{ -6.8}$ %. This is consistent with the lower end of the results discussed above.

We recomputed the multiplicity-corrected CDF-N 850 μm number counts above 3.5 mJy based on these observations using

$\begin{eqnarray} \frac{dN_{{\rm corr},i}(S)}{dS} &=& \frac{dN_{{\rm orig},i}(S)}{dS} \times (1-f_{{\rm mul},i}(S))\nonumber\\ && +f_{{\rm mul},i}(2S)\times 4\times \frac{dN_{{\rm orig},i}(2S)}{dS}, \end{eqnarray} \tag{ 4 }$

where f_{mul, i} = N_mul/N_sma represents the multiple fraction of the SMA detected SCUBA-2 sources in each flux bin i, and dN_{corr, i}/dS and dN_{orig, i}/dS are the multiplicity-corrected and the original SCUBA-2 counts, respectively. We have assumed for simplicity that the source splits into two equal components. This procedure conserves the EBL contribution of the counts, as it should.

We plot the results (black circles) in Figure 11, along with the original CDF-N SCUBA-2 counts (red circles). The multiplicity-corrected counts are essentially unchanged in the first two bright bins, since all the SMA observations on these sources correspond to single detections. In the last two faint bins, owing to the effect of the multiplicity, the corrected counts differ by a small amount from the original ones. However, the systematic changes introduced by the multiplicity are smaller than the statistical noise in our counts determination.

We compare our corrected counts with the ones obtained by the ALMA follow-up observations of the LABOCA sources (blue circles in Figure 11; Karim et al. 2013). The ALMA/LABOCA counts are systematically lower than ours, especially at the bright end above 9 mJy, where they argued that the counts drop significantly due to the fact that the all of their bright sources tend to break into fainter multiple sources. This tendency is not confirmed by our observations. In fact, all of our SMA observed single-dish detected bright sources (>9 mJy) show single SMA detections. In addition, the low counts in the ALMA/LABOCA determination are a partial consequence of the low LABOCA counts in the CDF-S, which we discussed in the previous subsection.

The multiple fraction could be affected by the depth of the follow-up observations. In principle, the deeper the interferometric observations, the more likely that a multiple system would be revealed. The sensitivity of our SCUBA-2 and SMA observations are comparable, with a mean depth ratio of ∼1.1. By comparing the fluxes measured by SCUBA-2 and by the SMA on the sources with a single SMA detection, we can estimate how many sources could be further resolved into multiples if the follow-up interferometric observations were deeper. We show the flux comparison in Figure 12, where most of the SMA measurements statistically agree with the those made by SCUBA-2. That a significant number of SMA measurements are larger than the SCUBA-2 measurements could be caused by the calibration uncertainties (The average flux ratio of the SMA to the SCUBA-2 is ∼1.1). We cannot rule out the possibility that sources with fluxes fainter than the sensitivity limit of our observations are contributing a small fraction of the flux of a SCUBA-2 source and the multiple fraction could be higher, but such faint sources would be unlikely to affect our counts within our flux range of interest.

**Figure 12.** SMA 860 μm flux vs. SCUBA-2 850 μm flux for the CDF-N SCUBA-2 4σ detected sources with single SMA detections. Error bars show 2σ uncertainties.
Download figure:
Standard image High-resolution image

In our counts calculation, we did assume that the multiple fraction obtained from the SMA observations is applicable to the overall sample. We caution that our counts below ∼5 mJy could change significantly if the multiple fraction were much higher, since our SMA observations primarily targeted the brighter SCUBA-2 sources.

We also compare our counts with the semi-analytical models recently presented by Hayward et al. (2013b) in Figure 11. Their single-dish (they adopted a 15'' FWHM beam size) counts predictions are plotted as the red curve, and their predictions for counts derived with arcsecond resolution interferometric observations are plotted as the blue curve. While their single-dish predictions agree very nicely with our results, their interferometric predictions are in general agreement with the ALMA counts but are significantly lower than ours. Our results and our discussion of the Karim et al. (2013) results suggest that the models proposed by Hayward et al. (2013b) may be overestimating the multiplicity correction.

5. DISCUSSION

Since the discovery of the FIR EBL, resolving the diffuse emission into discrete sources has been one of the primary goals of FIR surveys. In particular, we would like to determine the flux range of the galaxies that contribute the bulk of the submillimeter light. Thanks to gravitational lensing, our deep SCUBA-2 observations of cluster fields are able to probe to fluxes of ∼1 mJy at 450 μm and ∼0.1 mJy at 850 μm, an unprecedented depth at 450 μm. While showing the counts from each field in the same figure as we did in Figure 8 enables us to visualize the effect of cosmic variance, the amount of EBL resolved based on the model fits to the counts in a single field are poorly constrained (∼20%–200% at 450 μm, for example), due to the small number statistics on each field and also cosmic variance. To better constrain the counts models, we combine the counts from all five fields. In each flux bin, we average the corrected, statistically significant differential number counts from each field using an area weighting. We then estimate the Poisson errors using the combined source counts from the contributing fields. We plot the results in Figure 13. As the figure shows, errors are smaller on the flux bins that are covered by more than one field, and the minimum χ² fits on these counts (black curves with errors in gray shading) are much better constrained than the ones in Figure 8. The reduced χ² are 1.2 and 1.3 at 450 and 850 μm, meaning the counts are well described by the models, which are summarized in Table 5. The area-weighted combined differential counts are given in Table 6.

**Figure 13.** Area-weighted combined counts from all five fields. The best fit broken power law models and the number counts are summarized in Tables 5 and 6, respectively.
Download figure:
Standard image High-resolution image

Table 5. Best χ² Fits on the Combined Differential Number Counts at 450 and 850 μm

Wavelengths	N₀	S₀	α	β
(μm)	(mJy⁻¹ deg⁻²)	(mJy)	α	β
450	$22^{+51}_{-15}$	$21.1^{+9.0}_{-7.0}$	$2.73^{+0.0}_{-0.0}$	$5.77^{+9.5}_{-2.00}$
850	$120^{+180}_{-70}$	$6.2^{+2.0}_{-2.0}$	$2.26^{+0.2}_{-0.2}$	$3.79^{+2.00}_{-1.00}$

Download table as: ASCII Typeset image

Table 6. The Combined Differential Number Counts at 450 and 850 μm

S₄₅₀	dN/dS	S₈₅₀	dN/dS
(mJy)	(mJy⁻¹ deg⁻²)	(mJy)	(mJy⁻¹ deg⁻²)
1.26	$99640^{+67202}_{-60525}$	0.10	$5436690^{+4016800}_{-2727500}$
4.24	$1504^{+689.2}_{-689.2}$	0.23	$133324^{+103520}_{-72100}$
11.23	$78.74^{+47.92}_{-47.92}$	0.53	$36446^{+26232}_{-23530}$
16.31	$58.09^{+16.95}_{-16.95}$	1.05	$4729^{+1698}_{-1698}$
22.65	$14.60^{+4.77}_{-4.77}$	1.72	$2021^{+553.1}_{-553.1}$
30.74	$2.82^{+2.23}_{-1.35}$	2.64	$829.4^{+101.3}_{-101.3}$
39.69	$0.47^{+1.09}_{-0.39}$	4.05	$308.8^{+31.50}_{-31.50}$
...	...	6.20	$131.1^{+15.34}_{-15.34}$
...	...	9.51	$20.03^{+5.73}_{-4.55}$
...	...	13.94	$8.62^{+5.83}_{-3.72}$
...	...	20.31	$1.50^{+3.46}_{-1.24}$

Download table as: ASCII Typeset image

Based on the combined counts models, we measure $113.9^{+49.7}_{-28.4}$ ( $37.3^{+21.1}_{-12.9}$ ) Jy/Degree² at 450 (850) μm, which corresponds to a full resolution of the 450 and 850 μm EBL based on the measurements by Puget et al. (1996). However, our measurements also correspond to $80.0^{+34.9}_{-19.9}$ ( $85.7^{+48.4}_{-29.6}$ )% of the 450 (850) μm EBL if we adopt the EBL measurements from Fixsen et al. (1998). The percentage ranges expand to 48%–153% (44%–178%) if we consider the errors on the EBL measurements. We note that the slope of the faint end counts results in a weak divergence in both bands so that the counts must flatten at fainter fluxes than those measured here. The slopes must turn over before 0.5 mJy at 450 μm and 0.1 mJy at 850 μm in order not to significantly exceed the measured EBL.

We plot the cumulative EBL as a function of flux in Figure 14, where we show that 90% of the 450 μm EBL is contributed by the sources with S_450 μm < 10 mJy, and the majority (∼70%) comes from the sources with 1 mJy <S_450 μm < 10 mJy. We computed the IR luminosities (L_IR) for the dominant sources that contribute the bulk of the 450 μm EBL by assuming the average redshift (z = 1.3; Geach et al. 2013) and a single temperature modified blackbody spectral energy distribution (SED) with a fixed dust emissivity (β = 1.5). The dust temperatures of the SMGs are diverse, ranging from 25−50 K (Magnelli et al. 2012); Thus, we determined all possible L_IR values within this dust temperature range for SMGs with 450 μm fluxes from 1−10 mJy. We also determined L_IR values for the dominant 850 μm sources (these have fluxes between 0.1−1 mJy) by assuming z = 3.0 (Barger et al. 2012; Vieira et al. 2013; Hayward et al. 2013b). At both wavelengths, the majority of the dominant sources have L_IR between 10¹¹ and 10¹² L_☉, analogous to local luminous infrared galaxies (LIRGs; Sanders & Mirabel 1996). Note that this result is not sensitive to the assumed redshifts. While most SMG studies have focused on the more luminous galaxies that can be detected in single-dish observations of blank fields, it is these lower luminosity galaxies that dominate the submillimeter EBL.

In Figure 14 we also show other EBL measurements from the literature at both wavelengths. They all agree nicely with our results, except for the one from the confusion limited Herschel SPIRE 500 μm survey (Oliver et al. 2010). This is almost certainly a consequence of the fact that the results from Oliver et al. (2010) are based on a wide area survey where they are able to recover very bright but rare sources that are missed by our observations. However, the contribution of these bright sources to the EBL is small, and correcting for them does not significantly change the calculation of the EBL contribution of the full counts.

We also compare our results with model predictions. The empirical model predictions based on the IR SED templates by Béthermin et al. (2012a) are plotted as the red curve in both panels in Figure 14. The predictions match our results very well at 450 μm, while the ones at 850 μm are slightly lower at the bright end. While there is no treatment of multiplicity in the model by Béthermin et al. (2012a), the semi-analytical model by Hayward et al. (2013b) took this issue into account, and their single-dish 850 μm predictions (blue solid curve) agree with our results from ∼1–2 mJy but are low at higher fluxes.

Hayward et al. (2013b) also predicted the counts without any blending (blue dashed curve). At the bright end these are almost an order of magnitude lower than the single-dish results. However, as shown in Section 4.2, subject to the uncertainty on the multiple fraction, as well as the apparent disagreement between the interferometric counts at 850 μm, it is hard to draw a solid conclusion about the true shape of the 850 μm counts. In fact, our interferometric counts are statistically indistinguishable from our single-dish counts. Future interferometric follow-up of full submillimeter samples should resolve this issue.

Interestingly, without considering source blending and by fitting the counts from the mid-infrared to the radio, Béthermin et al. (2012a) were able to reasonably model the 450 μm counts. This could mean that the majority of the 450 μm sources are at z < 3, in agreement with recent deep SCUBA-2 COSMOS 450 μm results by Geach et al. (2013), Casey et al. (2013), and Roseboom et al. (2013).

6. SUMMARY

We have presented very deep number counts at 450 and 850 μm obtained using the SCUBA-2 submillimeter camera mounted on the JCMT. Based on a mixture of cluster lensing fields (A1689 and A370) and blank fields (CDF-N, CDF-S, and COSMOS), we have determined the number counts over a wide flux range. At 450 μm, we detected significant counts down to ∼1 mJy, an unprecedented depth at this wavelength. By integrating the number counts to our flux limits, we find that we have measured $113.9^{+49.7}_{-28.4}$ Jy/Degree² of the 450 μm EBL, which corresponds to $80.0^{+34.9}_{-19.9}$ % of the EBL measured by Fixsen et al. (1998) using the COBE satellite. The results show that the majority of the 450 μm EBL is contributed by the sources with S_450 μm between 1–10 mJy, and these sources are likely to be the ones that are analogous to the local LIRGs. At 850 μm, we resolved $37.3^{+21.1}_{-12.9}$ Jy/Degree² of the EBL, corresponding to $85.7^{+48.4}_{-29.6}$ % of the COBE measurement by Fixsen et al. (1998). If we consider the large uncertainties on the COBE measurements, the uncertainty of the percentage of the EBL resolved expands to 48%–153% (44%–178%) at 450 (850) μm. Our SCUBA-2 observations revealed statistically different number counts in the CDF-S field compared with the ones obtained by the LABOCA observations, which could be explained by the discrepancy in the number of sources discovered by the two observations. We also showed that there is little field-to-field variance and that source blending (multiplicity) has only a small effect on the shape of the counts determined with SCUBA-2.

We thank the referee for a careful review to improve the manuscript, Axel Weiß, Ian Smail, and Fabian Walter for giving helpful comments on the LABOCA counts toward the CDF-S region, and Matthieu Béthermin for the advice on the model predictions. We also like to thank Tim Jenness, Per Friberg, David Berry, Douglas Scott, and Edward Chapin for helpful discussions on data reduction; JAC support astronomers Iain Coulson, Holy Thomas, and Antonio Chrysostomou; and JCMT telescope operators Jim Hoge, William Montgomerie, Jan Wouterloot, and Callie McNew. We gratefully acknowledge support from NSF grant AST 0709356 (C.C.C., L.L.C.), the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation (A.J.B.), the David and Lucile Packard Foundation (A.J.B.), and the National Science Council of Taiwan grant 99-2112-M-001-012-MY3 (W.-H.W.). C.M.C. is generously supported by a Hubble Fellowship provided by Space Telescope Science Institute, grant HST-HF-51268.01-A. This research is also supported by a Grant-In-Aid of Research from the National Academy of Sciences, administered by Sigma Xi, The Scientific Research Society. The James Clerk Maxwell Telescope is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the United Kingdom, the National Research Council of Canada, and (until 2013 March 31) the Netherlands Organisation for Scientific Research. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

RESOLVING THE COSMIC FAR-INFRARED BACKGROUND AT 450 AND 850 μm WITH SCUBA-2

Article metrics

Permissions

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION

2. OBSERVATIONS AND DATA REDUCTION