SCUBA-2 Ultra Deep Imaging EAO Survey (STUDIES): Faint-end Counts at 450 μm

The SCUBA-2 Ultra Deep Imaging EAO Survey (STUDIES) is a three-year JCMT Large Program aiming to reach the 450 μm confusion limit in the COSMOS-CANDELS region to study a representative sample of the high-redshift far-infrared galaxy population that gives rise to the bulk of the far-infrared background. We present the first-year data from STUDIES. We reached a 450 μm noise level of 0.91 mJy for point sources at the map center, covered an area of 151 arcmin2, and detected 98 and 141 sources at 4.0σ and 3.5σ, respectively. Our derived counts are best constrained in the 3.5–25 mJy regime using directly detected sources. Below the detection limits, our fluctuation analysis further constrains the slope of the counts down to 1 mJy. The resulting counts at 1–25 mJy are consistent with a power law having a slope of −2.59 (±0.10 for 3.5–25 mJy, and − 0.7 + 0.4 for 1–3.5 mJy). There is no evidence of a faint-end termination or turnover of the counts in this flux density range. Our counts are also consistent with previous SCUBA-2 blank-field and lensing-cluster surveys. The integrated surface brightness from our counts down to 1 mJy is 90.0 ± 17.2 Jy deg−2, which can account for up to 83 − 16 + 15 % of the COBE 450 μm background. We show that Herschel counts at 350 and 500 μm are significantly higher than our 450 μm counts, likely caused by its large beam and source clustering. High angular resolution instruments like SCUBA-2 at 450 μm are therefore highly beneficial for measuring the luminosity and spatial density of high-redshift dusty galaxies.


INTRODUCTION
Since the advent of the Submillimeter Common User Bolometer Array (SCUBA, Holland et al. 1999) on the James Clerk Maxwell Telescope (JCMT) and the discovery of the first submillimeter galaxies (SMGs, Smail, Ivison, & Blain 1997;Hughes et al. 1998;Barger et al. 1998;Eales et al. 1999) two decades ago, tremendous progress has been made in understanding this important dusty galaxy population (see reviews by Blain et al. 2002;Casey, Narayanan, & Cooray 2014).Wide-field 850 µm surveys made with SCUBA and other bolometer array cameras have provided large samples of SMGs, while interferometric followup observations have enabled counterpart identification and detailed studies.
It is important to point out that attempts to understand the SMG population has been overwhelmingly focused on the bright end at 850 µm.The poor angular resolution of single-dish telescopes (e.g., FWHM = 15 ′′ for JCMT at 850 µm) produces a relatively bright limit due to source blending (S 850 ∼ 2-3 mJy).This effect is known as "confusion" (Condon 1974); it prevents the detection of faint sources and the full resolution of the extragalactic background light (EBL).In the millimeter and submillimeter (mm/submm) bands, sources detected in single-dish confusion-limited blank-field surveys typically account for only 10-40% of the EBL (Barger, Cowie, & Sanders 1999;Borys et al. 2003;Greve et al. 2004;Wang, Cowie, & Barger 2004;Coppin et al. 2006;Weiß et al. 2009;Scott et al. 2010;Hatsukade et al. 2011;Geach et al. 2017), and the bulk of the EBL remains unresolved.Surveys in strong lensing cluster fields can nearly fully resolve the mm/submm EBL (Cowie, Barger, & Kneib 2002;Smail et al. 2002;Knudsen, van der Werf, & Kneib 2008;Johansson, Sigurdarson, & Horellou 2011;Chen et al. 2013b;Hsu et al. 2016) and provide insight into the nature of the faint sources (Chen et al. 2014;Hsu et al. 2016).However, the sample sizes of such strongly lensed sources remain small, and source-plane expansion and magnification bias may make cosmic variance a stronger effect in lensing-cluster surveys.
Recently, ALMA has realized its tremendous sensitivity.Its high resolution makes it essentially confusion-free and able to detect sources that comprise the bulk of the mm/submm EBL.However, because of the small primary beam of its antennas, unbiased ALMA surveys are only able to image a few arcmin 2 to substantial depths (Umehata et al. 2015;Hatsukade et al. 2016;Aravena et al. 2016;Dunlop et al. 2017).Larger samples of fainter sources have been serendipitously detected in ALMA archival data (Hatsukade et al. 2013;Ono et al. 2014;Carniani et al. 2015;Fujimoto et al. 2016) including those in the ALMA calibration fields (Oteo et al. 2016), but biases caused by clustering and cosmic variance on these small scales are a potential concern.
The 450 µm window for deep submillimeter surveys was truly opened up by SCUBA-2 (Holland et al. 2013).At 450 µm, the nearly 2 times increase in angular resolution from 850 µm makes SCUBA-2 much less confusion limited.This enables the direct detection of fainter SMGs.This also makes multiwavelength counterpart identification less ambiguous.Previous 450 µm SCUBA-2 surveys reached sensitivities (1 σ noise) of 1 to 10 mJy in various blank fields and in lensing clusters.Among these, the deepest blank-field surveys (Geach et al. 2013;Zavala et al. 2017) resolved 20-30% of the 450 µm EBL into point sources down to 6 mJy.Geach et al. (2013) also estimated that the confusion noise of JCMT/SCUBA-2 at 450 µm is approximately 1 mJy, but a detection limit of a comparable flux has not been reached by any blank-field surveys.On the other hand, lensing-cluster surveys (Chen et al. 2013a,b;Hsu et al. 2016) had reached intrinsic fluxes of ∼ 1 mJy on a smaller sample of highly magnified sources, and nearly fully resolved the EBL.For comparison, confusion limited Herschel SPIRE imaging (FWHM ≃ 30 ′′ at 500 µm) can detect sources brighter than around 20 mJy at 250, 350, and 500 µm (Glenn et al. 2010;Clements et al. 2010;Oliver et al. 2010;Valiante et al. 2016), and these sources account for only about 15% of the FIR EBL at these wavelengths.Unlike the 850 µm and millimeter bands, the 450 µm band has weaker negative K-correction and is less sensitive to z > 4 galaxies, because the peak of the redshifted dust spectral energy distribution (SED) shifts out of the passband.However, for dust temperature of a few tens of Kelvin, the 450 µm band probes the SED peak of galaxies at z ≃ 1-2, the peak epoch of both cosmic star formation and quasar activity.There is clearly great potential for SCUBA-2 to reach much deeper than Herschel at 450 µm, and to detect the faint dusty galaxies that give rise to the bulk of the far-infrared (FIR) EBL and hence the bulk of the cosmic star formation.
To exploit the potential of SCUBA-2, we initiated a new program, the SCUBA-2 Ultra Deep Imaging EAO (East-Asian Observatory) Survey (STUDIES; program ID: M16AL006).It is one of the new JCMT Large Programs under the operation of EAO, starting from late 2015.Its goal is to map a R ∼ 8 ′ region and get close to the confusion limit at 450 µm in three years with 330 hrs of observations.The survey field is at the center of the Cosmic Evolution Survey (COSMOS, Scoville et al. 2007) within the area of the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS, Grogin et al. 2011;Koekemoer et al. 2011).The STUDIES field is also within the wider and shallower SCUBA-2 survey of Casey et al. (2013).The southern shallower region of STUDIES overlaps with the 450 µm COSMOS pointing of the SCUBA-2 Cosmology Legacy Survey (S2CLS, Geach et al. 2013).STUDIES thus takes advantage of the extremely rich multi-wavelength data in the COSMOS-CANDELS region and the heritage of previous SCUBA-2 surveys.
In its first year of observations, STUDIES had accumulated approximately 40% of the data and reached an r.m.s.noise of 0.91 mJy at the map center at 450 µm.The image has a diameter of 15 ′ .The full 3-yr STUDIES image will be 1.6× deeper over the same area.The parallel 850 µm imaging will provide an image that is confusion limited over the entire area.An early science result from STUDIES can be found in Simpson et al. (2017), where we detected a "passive" galaxy at 450 µm and demonstrated that non-detections in Herschel bands do not rule out an active star-forming system.Multi-wavelength properties of the detected galaxies will be presented in our future papers.Because the progress in 2017 is much slower, caused by the poor weather and instrument servicing work, here we present the first-year STUDIES data and our improved constraints on the 450 µm source counts.
FIR and mm/submm number counts provide sensitive tests to galaxy evolution models (e.g., Baugh et al. 2005;Valiante et al. 2009;Béthermin et al. 2012b;Hayward et al. 2013;Cowley et al. 2015;Lacey et al. 2016;Béthermin et al. 2017) and require strong evolution in the properties of FIR luminous galaxies.The counts can be also compared with satellite EBL measurements, to examine whether there remain significant galaxy populations that are unaccounted for in the imaging surveys.Ultimately, when the imaging surveys are sufficiently deep, the integrated surface brightness from the resolved source counts can put the various satellite EBL measurements to the test.We derive the 450 µm counts with both direct detections of bright sources and a fluctuation analysis for faint sources using the first-year STUDIES data.In Section 2, we describe the SCUBA-2 450 µm observations and data reduction.In Section 3, we derive the source counts at S 450 > 3.5 mJy using direct 4-σ detections and with simulations.We refer to this flux density regime as the "bright end."In Section 4, we examine the image background fluctuation below 4 σ and use a maximum likelihood method (often referred as "P (D)" analysis) to constrain the counts below 3.5 mJy.We refer to this as the "faint end."In Section 5, we compare our counts with other counts measurements and models, and estimate the contributions to the 450 µm EBL.We summarize in Section 6.We analyze the various bias effects and source blending in our observations and source extraction in Section A, and verify that our results are not biased by source clustering at the scale of our beam size in Section B.

Observations
We carried out the 450 µm SCUBA-2 observations between December 30, 2015 and May 11, 2016, under the best submillimeter weather conditions on Maunakea (Band 1, τ 225 GHz < 0.05, or precipitable water vapor, PWV < 0.83 mm).SCUBA-2 also takes 850 µm data simultaneously, and the data will be presented elsewhere.The opacity was actively monitored with a water vapor radiometer using the 183 GHz water line at the telescope pointing direction (Dempsey et al. 2012), to ensure that all observations were carried out under Band-1 conditions.The median PWV is 0.615 mm and the 90th-percentile range is 0.440-0.803mm, corresponding to 450 µm opacities of τ 450µm = 0.768, and 0.590-0.958,respectively.Hourly telescope pointing checks and less frequent focusing were conducted on the nearby infrared source IRC+10216.Typical pointing offsets were within 1 ′′ , much smaller than the 450 µm beam.Nightly flux calibration was obtained by imaging Uranus, Mars, and the infrared sources CRL618, CRL2688, and Arp 220.We conducted the imaging using the "CV Daisy" scan pattern, 1 which creates a circular map of R 6 ′ with a R ≃ 1. ′ 5 deep core at the center and rapidly increasing noise at R > 5 ′ .Each scan spanned 30 minutes in time.The scans were slightly offset from each other to even out the effects of noisy bolometers.This also slightly expands the map to R ≃ 7. ′ 5.The total on-source time was approximately 120 hrs, accounting for 40% of the total allocated integration of STUDIES.

Data Reduction
We performed the data reduction using the Sub-Millimeter Common User Reduction Facility (SMURF, Chapin et al. 2013) and the PIpeline for Combining and Analyzing Reduced Data (PI-CARD, Jenness et al. 2008).The time-stream data from the SCUBA-2 bolometers contain noise and signal from the background (atmosphere and ambient thermal emission), as well as astronomical objects.To extract the astronomical signal from the time streams and to map the results onto a celestial projection, we adopted the Dynamic Iterative Map-Maker (DIMM, Jenness et al. 2011;Chapin et al. 2013) routine of SMURF.We used the standard "blank field" configuration file, which aims to detect extremely low signal-to-noise point sources from deep observations.First, flat-fields were applied to the time streams by using the flat scans that bracket each science observation.The flat-field procedure subtracts a polynomial baseline fit from each time stream and scales the data to units of pW.Next, DIMM enters an iterative stage that attempts to fit the data with a model comprising a common-mode background signal, astronomical signal, and noise.The iterations were repeated only four times, as we do not expect a single science scan to have any significant signal to well constrain the model.We verified that all the final scans do not change significantly with further iterations, meaning that the residual between the model and data has converged.
To obtain the absolute flux scale, we measured the flux conversion factors (FCFs) from a subset of 140 standard submillimeter calibrators that were observed during the STUDIES campaign, after excluding extreme values (usually from early evening and morning observations).We obtained a mean FCF of 490 Jy beam −1 pW −1 and an r.m.s.scatter of ±136 Jy beam −1 pW −1 .The nominal value for SCUBA-2 at 450 µm is 491 ± 67 Jy pW −1 (Dempsey et al. 2013).Our scatter is twice larger than this and those in earlier SCUBA-2 450 µm studies in the literature.An about twice larger yearly scatter is also seen in the SCUBA-2 Calibration Database2 for 2016.The large scatter in our FCF therefore indicates a real night-to-night variation, rather than problems in our data or our reduction.However, other than the larger scatter, we do not observe a long-term trend in the mean FCF, and our mean FCF is consistent with the nominal value.We therefore applied the value of 490 Jy beam −1 pW −1 to our data.In general, we expect a 12% uncertainty in flux calibration (Dempsey et al. 2013), which is a combination of a 5% uncertainty in the absolute calibration of planetary models and another 10% uncertainty in the determination of the FCF.
After each scan was reduced and flux calibrated, we adopted the MOSAIC JCMT IMAGES recipe from PICARD to combine all the products into a final map.To optimize the detection of point sources, we convolved the map with a Gaussian kernel that is matched to the instrumental point-spread function (PSF).For this, we adopted the PICARD recipe scuba2 matched filter, which first smooths the map by convolving it with a 20 ′′ Gaussian kernel and subtracts this smoothed map from the image to remove any large-scale structure.Then a normalized Gaussian kernel with a FWHM set to the diffraction limit of the telescope was used to convolve the map to obtain the optimal pointsource S/N for each pixel (Stetson 1987;Serjeant et al. 2003).The same process was also applied to the calibrators before the FCF was measured.Also, although the above procedure builds in a flux adjustment to compensate the flux loss caused by the subtraction of the 20 ′′ -smoothed image, additional adjustment may be required.In order to see the effect on source fluxes of the blank-field configuration of DIMM and the above convolution procedures, we inserted artificial point sources with the instrumental FWHM into the data streams.For source flux densities from sub-Jy to a few Jy, we found that the above reduction procedures attenuated the flux density by 6.2%, on average.We therefore made a 6.2% flux adjustment to the map.This is slightly less than the 10% adjustment reported by Geach et al. (2013) and Chen et al. (2013b), but the difference is well within the generally accepted 10% calibration uncertainty.Our final match-filtered map has a noise at the map center of 0.91 mJy for a point source, and the noise increases to around 10 mJy toward the map edge at a radius of 7. ′ 5.The 0.91 mJy sensitivity is about 30% deeper than that in the deepest 450 µm map in the literature (Zavala et al. 2017), and is comparable to the deepest 450 µm pointing in S2CLS (Geach et al., in prep.).We present our S/N map in Fig. 1, and the histograms of the pixel flux in Fig. 2. In Fig. 1, we can see negative rings around the brightest sources, produced by the above match-filtering process.The negative ring in the PSF is responsible for the negative excess in Fig. 2 comparing to pure noise.We verify the astrometry of our image by comparing the source positions in our image and in the Very Large Array 3 GHz image (Smolčić et al. 2017).Among the 98 4-σ sources extracted from our image (Section 3.1), 59 have 3 GHz counterparts with 4 ′′ search radii.We expect 1.3 of them to be chance alignments, given the spacial density of the radio sources.The mean positional offset between 450 µm and 3 GHz is 0. ′′ 3 along R.A. and 0. ′′ 2 along Dec., and the dispersions are ∼ 1. ′′ 5.These are all much smaller than the 450 µm beam.There is thus not an apparent astrometric offset in our image.

Source Extraction
We first constructed a PSF for source extraction.We generated a synthetic PSF by averaging the 10 highest signal-to-noise ratio (S/N) sources in our final match-filtered map.These 10 sources contribute equally to the averaged PSF.The resulting PSF has a FWHM of 10. ′′ 1, which is 9% broader than an idealized PSF (i.e., instrumental PSF with the same match-filtering), likely caused by the minor smearing effects from pointing errors and focus changes.In the steps described below, we tried using the synthetic PSF and the idealized PSF.The results are in agreement with each other, except that the idealized PSF has a slight tendency to "detect" more faint sources very close to bright sources.This is likely caused by its narrower profile and under-subtraction of the outer part of bright sources.Since the synthetic PSF contains the smearing effect from observations and should be more realistic, we adopt the synthetic PSF for the subsequent analyses.
To extract sources, we employ an iterative procedure that is similar to "CLEAN" deconvolution in radio interferometric imaging (Högbom 1974).We identify the peak pixel in the S/N map and subtract 5% of a peak-scaled synthetic PSF at its position.This 5% subtraction fraction is often called CLEAN "gain," and is typically set to a few to 10-percent in order to achieve stable convergence.We record its coordinates and subtracted flux.We identify the next peak in the image and repeat the process until we meet a floor S/N threshold, which we set to be 3.5 σ.During the process, if a subsequent S/N peak is located within 4 ′′ (approximately half the beam FWHM) from a previously identified peak, we consider them as the same source.In such a case, we only perform the iterative subtraction at the original position.Otherwise, we consider it as a new source and perform the PSF subtraction at the new position.Since the cleaning of source flux stops at 3.5 σ, we sum up the cleaned fluxes and the remaining 3.5 σ flux to be the final flux for each individual source.The very outer part of the map receives much shorter integration and contains arc and stripe patterns that are clearly not random noise.Visual inspection shows that sources extracted in such areas are not always convincing.We therefore only include sources extracted within 7. ′ 5 (r.m.s.noise < 10 mJy) from the map center in our final catalog.We detect 141 sources at 450 µm at > 3.5 σ, among which 98 are > 4.0 σ (circles in Fig. 1).
The above iterative algorithm is inspired by the radio CLEAN deconvolution, and is similar to that in Wang, Cowie, & Barger (2004).It has the capability of separating mildly blended sources whose separation is larger than roughly one beam FWHM.For example, a pair of blended SCUBA sources in Wang, Cowie, & Barger (2004) was subsequently confirmed by the interferometric imaging of Wang et al. (2011) to be multiple sources.In our 3.5 and 4-σ catalogs, there are seven and five pairs, respectively, whose separations are less than 12 ′′ (see Fig. 1 for the 4-σ cases).They will be targets of our future followup observations, to test the accuracy of our source extraction in such limiting cases.We also note that when the CLEAN gain is set to 1.0, our source extraction reduces to the standard source identification adopted by other teams for SCUBA-2 images.The few close pairs would have been identified as single sources if the gain were 1.0.In our simulations (Section A), this would increase the fraction of sources that are blends of multiple faint sources (hereafter the "multiple fraction") of the catalog.Also, in such a case, fluxes can be over-estimated for blended sources, since their fluxes are contaminated by the neighbors.However, the previous SCUBA-2 450 µm surveys are shallower and further from the confusion limit.Source blending and therefore the flux contamination may not be an issue there.The noise distribution in our image is roughly axisymmetric, so area can also be approximately mapped to radius, which is shown with the right-hand y-axis.

Number Counts and Simulations
The map center has a noise of 0.91 mJy beam −1 .The noise slowly increases to ∼ 2 mJy beam −1 , where the map area reaches 100 arcmin 2 .After that, the noise increases more rapidly toward the outer region.For our number counts, we do not use image area that is more than five times shallower than at the map center, indicated by the vertical dashed line.
We derived 450 µm source counts using our extracted sources and the noise map.We only use 4 σ sources for this calculation because they are more secure detections.We will discuss the reliability of the 4 σ sources in Section A. Furthermore, to ensure reliable source counts on faint sources, we confine our calculations to areas where the noise is less than five times the noise of the map center.This further discards the most noisy 9% of the area used for the source extraction.The cumulative area as a function of noise is shown in Fig. 3, and the area involved in our estimates is 151 arcmin 2 , indicated by the vertical dashed line.97 of the 98 4-σ sources fall in this area.
In the number count calculation, each source contributes 1/(A e (S)dS) to the counts, where A e is the effective area where the source can be detected at > 4 σ given its measured flux density S, and dS is the flux density interval.Essentially, A e (S) is the function in Fig. 3 when the x-axis is multiplied by 4. The error is assumed to be Poissonian.The faintest flux density bin is discarded, because it only contains one source.The resulting raw differential counts are presented in Fig. 4 (open squares) and Table 1 (column 3).
The raw counts suffer from various observational biases: flux boosting produced by the Eddington bias and faint confusing sources; detection incompleteness; spurious sources; and source blending.These are analyzed in detail in Section A. To overcome these, we performed Monte Carlo simulations   A).The solid black curve is the input Schechter counts in the simulations.The dashed black curve is the measured output counts in the simulations, and can also be described by a Schechter function.The two dotted black curves show the 68th-percentile range of the outputs of the 200 realizations.The 68th-percentile range reasonably matches the observed error bars derived based on the assumption of Poisson errors.Two parameterizations to the corrected counts are shown with colored curves.The red solid curve is a Schechter fit to the corrected counts.The reduced χ 2 of this fit is 0.45 (N df = 6).This curve is almost indistinguishable from the input Schechter function (black solid curve), showing that our iterative procedure converges.The green dashed line is a power-law fit to the corrected counts after excluding the brightest flux density bin.The reduced χ 2 of this fit is 0.38 (N df = 6).
to derive the intrinsic counts.We first created a "true noise" map using the jackknife method, described by Cowie, Barger, & Kneib (2002) for SCUBA imaging.We divided the half-hour SCUBA-2 scans into two interlacing halves and made two maps of nearly identical area coverage and sensitivity.The two half-maps underwent the same beam-matched convolution as the full map.They were subtracted from one another, and each pixel of the resultant map was scaled by √ t 1 t 2 /(t 1 + t 2 ) where t 1 and t 2 are the noise-weighted integration times of that pixel in the two half-maps.(t 1 and t 2 are not identical because a fixed sky position can be swept by different bolometers in different scans.)This effectively removes any faint sources and provides a map whose noise distribution is identical to the real image.We measured rms noise locally on the true noise map and found it consistent with the r.m.s.map generated by SMURF.The brightness distribution in the true noise map is Note-S obs is the observed flux density.S corr is the flux density corrected for boosting (Section A).For the differential counts, flux densities are the center of the bins.For the cumulative counts, flux densities are the lower ends of the bins.
presented in Fig. 2 (red histogram).It is symmetric about zero.It is Gaussian-like, but not exactly a Gaussian.Instead, it is a sum of Gaussians of various widths, because of the non-uniform sensitivity distribution, We then created simulated images using the synthetic PSF.We randomly placed scaled PSFs in the true noise map with assumed source counts (see below) and flux densities between 1 and 50 mJy.We found that this 1-50 mJy flux density range is sufficient, and the results do not change if we expand the range for the input sources (see Section 4).In each fine flux density bin, the number of the simulated sources is determined by the assumed counts plus a Poissonian fluctuation.We created 100 simulated images using each of the positive and negative true noise maps.Because all of the effects of flux boosting, completeness, spurious sources, and source blending depend on the intrinsic source counts, the selection of the intrinsic counts in the simulations may be crucial.We employed an iterative procedure to approach the intrinsic counts from the observed raw counts.We started with assuming that the intrinsic counts are the observed raw counts; we fit the raw counts with a Schechter function, and used the fitted function as input to create the first set of 200 realizations.We then ran source extraction on the simulated images, derived simulated output counts, and compared them with the observed raw counts.We calculated the ratio between the simulated and observed counts, used that to adjust the input counts, and repeated the simulations.In the first two iterations, the simulated output counts fluctuated around the observed counts, and then converged after the third iteration.In the third iteration, the required adjustment was much smaller than the error bars in the counts.Our final corrected counts (C corr ) are calculated using the raw counts (C raw ) and the ratios between the input and output counts in the simulations (C sim,in and C sim,out , respectively): where S obs is the observed flux density of sources and S corr is the flux density corrected for boosting (see Section A).The corrected counts are presented in Table 1 (column 5) and in Fig. 4 (solid squares).
For readers' reference, we also present the corrected cumulative counts in Table 1 (column 7), derived with the same manner as the differential counts.
Our calculation in Eq. 1 is fundamentally different from what is often adopted in the literature for mm/submm single-dish source counts, where the counts are corrected with completeness and spurious fractions.For such corrections, one generally has to assume one-to-one relations between the input and output source lists.Such an assumption is not accurate for single-dish observations, because multiple sources in the input list can be blended in the image and detected as a single source of very different flux densities in the output list.In contrast, Eq. 1 does not rely on such an assumption.The only part in Eq. 1 that contains a one-to-one relation is the correction of observed flux density (S obs ) to intrinsic flux density (S corr ), and this is presented in Section A. However, this only slightly affects the interpretation of where the counts are measured in terms of intrinsic flux density, and does not change the measured amplitude and slope of the counts.
The corrected counts in Fig. 4 can be reasonably well fit with a Schechter function: with the best-fitting parameters listed in Table 2.The reduced χ 2 of the fit is 0.45, for a degree of freedom of N df = 6.The fitted function is the red solid curve in Fig. 4.This fitting is considered as the final (fourth) iteration to approach the intrinsic counts.It is almost indistinguishable from the input Schechter function used to create the simulations (black solid curve in Fig. 4), showing that the iteration converges very nicely.In the studies of Chen et al. (2013b) and Hsu et al. (2016), similar iterative procedures were also adopted and the fitting of the intrinsic counts was made as part of the analysis procedure.
In the literature, mm/submm counts are also often fitted with a broken double-power law.The turn-over point is typically between 20 and 30 mJy for 450 µm counts (Casey et al. 2013;Chen et al. 2013a,b;Hsu et al. 2016;Zavala et al. 2017), and 8.4 mJy for 850 µm counts constrained by ALMA (Simpson et al. 2015b).These two are consistent given that the S 850 /S 450 ratio is typically 2.5 with large scatter (Casey et al. 2013;Hsu et al. 2016).Our data seem to also suggest a turn-over point at 20-30 mJy, but we do not have sufficient data to constrain the slope in the brighter portion.If we fit our data with a broken power law, we get a turn-over flux density between our first and second brightest flux density bin, and thus no constraint on the bright-end slope.This is essentially fitting the counts with a single power law after excluding the brightest flux density bin.We therefore fit our counts at S < 25 mJy with a power-law form: with the best-fitting parameters listed in Table 2.The reduced χ 2 of the fit is 0.38 (N df = 6).This is the green dashed line in Fig. 4 and is a slightly better fit to the data at < 25 mJy, compared to the Note-The fitting is conducted on the bright-end counts above 3 mJy.For the power-law fit, the brightest flux-density bin is not included.
Schechter fit.If we exclude the brightest flux density bin and then perform a Schechter fit, we obtain a large S 0 (> 60 mJy) and strong degeneracies among the fitted parameters.This again shows that the data at < 25 mJy are consistent with a single power law.

FLUCTUATION ANALYSES FOR THE FAINT END
The derivation of the source counts in the previous section makes use of sources detected at 4 σ at above 3.5 mJy (intrinsic).The flux distribution of the residual map after 4-σ sources are removed contains additional information about fainter sources.In Fig. 5, we show the flux distribution of the pixels within 5 ′ of the map center in the residual map and in the true noise map.The noise level in this area is 0.91-1.5 mJy.Comparing to pure noise fluctuations, the residual map has excesses on both the positive and negative sides, caused by < 4 σ faint sources and the negative rings around their PSF.This shape can be used to further constrain the faint-end counts.To do this, we employ a parameter estimation method similar to that described in Patanchon et al. (2009, see discussion therein).In the literature, Maloney et al. (2005), Coppin et al. (2006), Weiß et al. (2009), andScott et al. (2010) also used similar methods to derive counts from their mm/submm images.The idea is to minimize the difference in the flux density distributions of the model-predicted and the measured residual maps, by maximizing the likelihood.In the literature, this is often referred as the P (D) method.
For each number count model, the likelihood of the data is the probability that all of the flux density values in the residual map had occurred, in a logarithmic form: .Distributions of pixel brightness in the residual maps from the observations (black), the true noise map (red), and three representative models from our fluctuation analysis that have high likelihoods (orange, green, and blue, corresponding to the three crosses in Fig. 6).Here we only show the pixels within 5 ′ of the map center, i.e., the pixels included in our fluctuation analysis.
where n i is the number of pixels in the i-th flux density bin, and p i (θ) is the integral of the PDF in the i-th bin.Eq. 5 is what we adopt for calculating the log-likelihood for a number count model, and the PDF is equivalently the flux density histogram of the pixels in an image.
For simplicity, we adopt a single power law for the faint-end counts, with only two parameters.The first one is the termination point of the counts, S min .The counts are assumed to be zero at S < S min and S min can be chosen to be much fainter than the detected sources.The second parameter is the faint-end power-law slope, denoted as α f .The power law is normalized to the corrected bright-end count at 3.94 mJy.This is the median flux density in the faintest flux density bin of our corrected bright-end counts.In principle, we could add a third parameter so that the junction point of the bright and faint ends does not have to be 3.94 mJy.Furthermore, we should allow a faint-end turnover at perhaps below 1 mJy in order to prevent the sources from over-producing the EBL, and this would require two additional parameters (the turn-over flux density and the slope in the extreme faint end).However, as we will show, the results of the two-parameter model are already highly degenerate, and our data do not allow meaningful constraints on additional parameters.
To determine the PDF of each number count model, we performed Monte Carlo simulations for the 450 µm images in a way identical to that described in Section 3.2.Because the clustering of faint 450 µm sources on the scales of our beam size is unknown (although likely to be weak), the spatial distribution of our simulated sources is random (cf.Vernstrom et al. 2014, and see discussion therein).For each faint-end count model, 200 simulated images were created using the true noise map, and another 200 using the negative true noise map.The 400 images underwent the same source extraction procedure as the observed image, and 4-σ sources were extracted and removed.Since our goal is to constrain the faint-end counts, we only use the central R < 5 ′ deep region, where the noise is less than approximately twice the noise at the map center.We also only used the flux density range of within ±4.2 mJy for the likelihood calculation, so the results are not skewed by the small number of brighter pixels.The resultant PDF averaged from the 400 model images was used as p(θ) in Eq. 5 for the calculation of the log-likelihood.
We calculated the log-likelihood over the parameter space of α f = −0.6 to −6.8 and S min from 0.0 to 2.75 mJy, with intervals of 0.2 in both.We neglected areas in this parameter space if there we observed a trend of monotonically decreasing log-likelihood toward the extreme parameter values, or if the model overproduces the COBE 450 µm EBL measured by Gispert, Lagache, & Puget (2000) by factors of more than four.The resultant relative log-likelihood distribution is presented in Fig. 6.There is a curved "ridge" of high likelihood extending from the upper-left to lower-right, representing the parameters that can reproduce the observed residual map.
To estimate the confidence levels, we used a bootstrap method.Instead of using the PDF combined from the 400 simulated images for each model, we calculated the PDF and log-likelihood of each simulated image and then determined the dispersion in the log-likelihood.Because we perturbed the number of sources in each flux density bin according to the Poisson distribution in the simulations and adopted both the positive and negative true noise maps, this dispersion represents the variance in the ensemble of images with the same set of model parameters and the underlying noise distribution.It therefore offers an estimate of the confidence range for the likelihood.This is the same method adopted by Scott et al. (2010).Around the maximum likelihood "ridge" in Fig. 6, the mean dispersion is 61.We therefore consider the 1, 2, and 3-σ confidence ranges as 61, 122, and 183, respectively, below the maximum log-likelihood.These are the white contours in Fig. 6.
To obtain a rough idea of how the above likelihood method can reproduce the observed residual map, we pick three sets of parameters (white crosses in Fig. 6) that are well within the 1-σ confidence range, and compare their PDFs with the observed one in Fig. 5.We find that all three models closely reproduce the PDF of the observed residual map.All of them are clearly different from pure noise.We thus conclude that sources fainter than our 4-σ detection limit are needed to explain the background fluctuation in the observed image.
Unfortunately, both Fig. 6 and Fig. 5 show that a broad range of model parameters can meet the requirement of reproducing the observed low-level fluctuations.The 1-σ confidence region in Fig. 6 is not localized.Instead, it shows a strong degeneracy between α f and S min .We can reproduce the observations with a relatively shallower faint-end power-law slope α f > −3 if the counts extend to S min < 1 mJy.We can also reproduce the observations with extremely steep power-law slopes (α f = −4 to −6) if the counts terminate at S min ≃ 1.5-2.5 mJy.One can argue that the latter is highly unlikely for two reasons: (1) α f < −4 is much steeper than the power-law slope we observed at the bright-end (α = −2.59± 0.10); and (2) observations in the lensing cluster fields clearly show that the nearly power-law counts should extend to at least 1 mJy (Chen et   2016, also see Fig. 7 and Section 5.1).Below we assume that the counts extend to S min < 1 mJy, and discuss the implication of Fig. 6.
At S min < 1 mJy, the contours in Fig. 6 become vertical.This means that adding sources fainter than 1 mJy to the simulated maps does not improve the results for a fixed power-law slope of α f = −2.6 +0.4 −0.7 .In other words, our data are not sensitive to either a termination or a turn-over of the counts if it occurs at < 1 mJy.At 1 mJy, the 1-σ range of the inferred counts is 2.5 × 10 4 to 1.2 × 10 5 deg −2 mJy −1 (Fig. 7).This is consistent with the lensing counts in Chen et al. (2013b) and Hsu et al. (2016), and also consistent with the power-law extrapolation of our bright-end counts.These counts also translate to 0.3-1.4sources per SCUBA-2 450 µm beam, i.e., sources are highly confused.This is consistent with the confusion limit of approximately 1 mJy estimated by Geach et al. (2013), and explains why our fluctuation analysis is no longer sensitive to sources of < 1 mJy.Scott et al. (2010) also found in their 1.1 mm analyses that adding sources fainter than the limit of roughly one source per beam does not alter the distribution.
Wang et al. Figure 7.Comparison of 450 µm counts (top panels), and counts renormalized with S 2.5 to better show the differences and changes in shape (bottom panels).Panel (a) compares results from this paper and other SCUBA-2 450 µm counts in the literature (Geach et al. 2013;Casey et al. 2013;Chen et al. 2013b;Hsu et al. 2016;Zavala et al. 2017).The solid squares are our counts derived from 4-σ sources, while the open square is our count derived from a fluctuation analysis assuming that the faint-end counts extend to ≃ 1 mJy (Section 4).The counts derived from various SCUBA-2 surveys are in excellent agreement across the entire flux density range, but the data points are not completely independent (see text).Panel (b) compares our counts with Herschel SPIRE 350 µm and 500 µm counts.The Herschel counts were derived from direct source detection at the 20 mJy bright end (Glenn et al. 2010;Clements et al. 2010;Oliver et al. 2010;Béthermin et al. 2012a;Valiante et al. 2016), and stacking analyses (Béthermin et al. 2012a) and P (D) analyses (Glenn et al. 2010) at the faint end.The x-axis shows the flux density measured at the corresponding wavebands, and we do not convert flux densities to a common wavelength.The insert shows the filter profiles of the SPIRE 350 µm (cyan) and 500 µm (orange) wavebands, and the SCUBA-2 450 µm wavebands (black).Although it is reasonable to expect the 450 µm counts to lie between the 500 µm and 350 µm counts, it is not the case here.At > 5 mJy, nearly all Herschel counts are consistently above the 450 µm counts.Béthermin et al. (2017) explained this with source clustering at the scales of Herschel beam sizes.Panel (c) compares our counts with 450 µm count models in Béthermin et al. (2012bBéthermin et al. ( , 2017) ) and Lacey et al. (2016).
To sum up, we do not find evidence of a turn-over of the counts between 1 and 3.5 mJy.A single power law with a slope of ≃ −2.6 can explain the observations between 1 and 25 mJy.Our fluctuation analysis is insensitive to the counts at below 1 mJy as long as the counts maintain a power-law slope between −1.9 and −3.1.A turn-over or a sudden termination of the counts at below 1 mJy can still be consistent with our data.Indeed, the counts must turn over at some point below 1 mJy, or the integrated surface brightness from the 450 µm sources will exceed the EBL.We will discuss this in Section 5.2.

Comparison with Previous Counts
We compare our counts with previous SCUBA-2 450 µm counts in the literature in Fig. 7 (a).Among the blank-field surveys, the Geach et al. (2013) and Zavala et al. (2017) surveys reached r.m.s.sensitivities of 2 mJy and 1.2 mJy, respectively, and their survey areas are comparable to ours (i.e., ≃ 150 arcmin 2 ).The Casey et al. (2013)  At the very bright end of 30 mJy, all the SCUBA-2 counts are consistent with each other and show a steeper fall-off.Here, our count has a large error, because of our smaller effective area.At the roughly 1 mJy faint end, our result derived from the fluctuation analysis and the extrapolation from our bright-end counts are both consistent with the lensing results of Chen et al. (2013b) and Hsu et al. (2016).As we mentioned in Section 4, there is no evidence for a faint-end turn-over in our data at above 1 mJy.The two lensing data points at about 1 mJy seem consistent with this.
The counts are best constrained in our data below between 3 mJy (the intrinsic flux density of our faintest detected source) and 25 mJy.Here, our results are consistent with previous ones within the error bars of each individual data point, except for those between 10 and 20 mJy from the shallower survey of Casey et al. (2013), which are higher than the other counts.Cosmic variance (i.e., the effects of clustering on small fields) does not easily explain the elevated counts in Casey et al., since their survey field fully encloses ours and that of Geach et al.For cosmic variance to be the explanation, the > 10 mJy sources would have to be strongly clustered at 5 ′ scales.Unfortunately, there do not exist deep and wide-field 450 µm dataset to test such clustering.A hint of clustering can be found in Casey et al. (2015) (also see Wang, Elbaz, et al. 2016), who detected a filamentary structure at z = 2.47 in their SCUBA-2 image.This structure, represented by their seven SCUBA-2 sources, runs through our field.Five of their sources are within our survey area and all of them are detected by us.Therefore, this structure cannot be responsible for the difference between our counts and the counts in Casey et al. (2013).Hung et al. (2016) reported another z ∼ 2.1 large-scale structure in this field.However, the only galaxy with a 450 µm detection in their sample is detected by both Casey et al. (2013) and us.So this structure does not seem to drive the higher counts in Casey et al. (2013), either.
In order to more quantitatively compare the counts in these surveys, we refitted the counts quoted by the various authors and derived the counts at 6 mJy and 10 mJy and the associated uncertainties.The results are presented in Table 3.The dispersions in these values are 127 and 18 deg −2 mJy −1 at 6 mJy and 10 mJy, respectively, or 21% and 10% relative to the mean.These are comparable to the shot noise in the counts, showing that the field-to-field variance should be even smaller.We further tested this using the simulated 2 deg 2 450 µm catalog of Béthermin et al. (2017), which contains source clustering, to conduct simulations of source extraction and number counts (see Section B for details).We found a scatter of about 10% in the counts between 6 and 10 mJy for our field size.Therefore both the existing data and the simulations of Béthermin et al. (2017) do not indicate a field-to-field variance that is much larger than 20% for 450 µm sources.However, we caution that the scatter seen in Table 3 is not exactly the field-to-field variance, since the some of the fields overlap, and the Chen et al. and Hsu et al. results are combined from multiple fields.This might underestimate the field-to-field variance, since counts from different fields are combined and averaged in these studies.A more proper evaluation of the field-to-field variance will be provided by the S2CLS 450 µm surveys (Geach et al., in prep.), along with the few 450 µm lensing fields.Nevertheless, it is still quite remarkable that the various 450 µm counts agree with each other at the roughly 20% level.Geach et al. (2017) reported a roughly 50% field-to-field scatter for 850 µm sources in S2CLS, whose field sizes are 0 • .5 to 1 • .The scatter is still roughly 50% when the S2CLS 850 µm counts are measured over 10 ′ fields at 3.5 mJy, of which roughly 60% is contributed by the Poisson errors in such smaller fields (J.Geach, personal communication).Both these results are much larger than the scatter we see in 450 µm counts.This suggests that the 450 µm sources are less clustered at scales larger than our field size ( 8 Mpc at z = 2).Two factors may contribute to this.First, 850 µm surveys are only sensitive to the most luminous dusty galaxies because of the confusion limit, while 450 µm surveys can reach sources with lower luminosities (i.e., weaker clustering).However, the source densities listed in Table 3 are not yet at the limit of 850 µm SCUBA-2 surveys.Therefore, this is unlikely to be a luminosity effect.Another reason is the strong negative K-correction at 850 µm, making the 850 µm band very sensitive to luminous dusty galaxies at z > 2. Wilkinson et al. (2017) show that such high-redshift SMGs are clustered more strongly than those at lower redshifts.Such high-redshift SMGs are less abundant in 450 µm images.This is evident in the redshift distributions of the 450 and 850 µm sources (Casey et al. 2013).It is therefore possible that the larger field-to-field variance at 850 µm is driven by the high-redshift SMGs.
In Fig. 7 (b), we show a comparison between our 450 µm counts and Herschel SPIRE 350/500 µm counts, along with the filter profiles of these three passbands.Given the passbands, is reasonable to expect the 450 µm counts to fall between the 350 and 500 µm ones, probably closer to the 500 µm ones.To our surprise, this is not the case, and the 450 µm counts fall below both the Herschel 350 and 500 µm counts.The Herschel confusion limit is approximately 20 mJy in both wavebands.The counts at > 20 mJy were derived with direct source extraction (Glenn et al. 2010;Clements et al. 2010;Oliver et al. 2010;Béthermin et al. 2012a;Valiante et al. 2016), while the counts at < 20 mJy were derived by stacking Spitzer 24 µm sources in the Herschel images (Béthermin et al. 2012a) and from P (D) analyses (Glenn et al. 2010).The Herschel counts derived by various authors are highly consistent with each other, regardless of how they were derived.All of them lie significantly above the SCUBA-2 450 µm counts over nearly the entire flux density range of interest.The 500 µm data points are individually about 1 to 2 σ away from the SCUBA-2 counts; however, the fact that collectively nearly all of them are above the SCUBA-2 points means that the over-abundance in the Herschel counts is significant.In the literature, over-abundant Herschel counts were already noted by Chen et al. (2013b), but Geach et al. (2013) and Zavala et al. (2017) considered the Herschel counts to be consistent with their SCUBA-2 counts, perhaps because of the larger error bars and sparse data points.With our better constrained SCUBA-2 450 µm counts, it should now be clear that there is a discrepancy with the Herschel -derived counts.
Béthermin et al. ( 2017) offered an explanation for the elevated Herschel counts.With simulations of dark matter halos in a 2 deg 2 lightcone, an abundance-matching technique to populate the halos with galaxies, and models of star formation and galaxy SEDs, these authors can reproduce the observed Herschel counts with much lower intrinsic counts.They therefore attribute the elevated Herschel counts to source blending in the Herschel images, caused by the clustering of FIR sources and the large Herschel beam in the SPIRE bands.We note that just blending alone does not necessarily bias the P (D) analyses and stacking analyses in the faint end, as the associated simulations would have accounted for its effects.However, the simulations typically do not include source clustering, which can amplify the observed image fluctuation and the stacking signal.Béthermin et al. (2017) also show that the intrinsic counts are much closer to the SCUBA-2 counts, because the SCUBA-2 450 µm beam is much smaller.Despite this, in Fig. 7 (c), we show that the 450 µm counts predicted by Béthermin et al. (2017), as well as by Lacey et al. (2016), are somewhat higher than the SCUBA-2 counts at flux densities above 3 mJy.The over-prediction in the counts is up to about 70% in the Béthermin et al. (2017) case.The slopes of the counts at flux densities below 3 mJy is also shallower in the models than what the observations suggest.The above 70% offset in counts can be translated to a 23% offset in flux, which is much larger than the generally expected 12% calibration uncertainty in SCUBA-2 observations (Dempsey et al. 2013).In addition, the offset between the model counts and observed counts is not a constant.It is not possible to explain this offset with a simple calibration error.The discrepancy between the models and observations thus appears to be real.Nevertheless, the work of Béthermin et al. (2017) and the comparisons in our Fig.7 (b) highlight the importance of high angular resolution in studying high-redshift galaxies.The beam size of SCUBA-2 at 450 µm is much less affected by clustering (Section B), and can therefore provide much more robust measurements of both the spatial density and flux density of sources.
Finally, we estimate the bias in the Herschel 500 µm counts, assuming that our 450 µm counts are unbiased.At around 10 mJy, the observed Herschel 500 µm counts are higher than our 450 µm counts by 1.25 times in flux density, or 1.8 times in number density.These are lower limits for the bias in the Herschel counts, since 500 µm counts should be intrinsically lower than 450 µm counts.Estimating the intrinsic offset between the 450 and 500 µm counts requires knowledge about the S 450 /S 500 flux density ratios of the sources, which are determined by their redshift distribution and SEDs.A zeroth-order estimate can be made through the model of Béthermin et al. (2017).Although Fig. 7 (c) shows that the model does not perfectly match the 450 µm counts, it nevertheless predicts 450 µm counts that lie between the 350 and 500 µm counts, meaning that the S 450 /S 500 flux density ratios in its source catalog are not too far off.If we use the intrinsic offset between the 450 and 500 µm counts in the model of Béthermin et al. (2017) to adjust the above observed offset, we then find that the observed Herschel 500 µm counts should be higher than the intrinsic 500 µm counts by roughly 1.4 times in flux, or 2.5 times in number density, at flux densities of about 10 mJy.

Contributions to the 450 µm EBL
We compare the 450 µm EBL resolved in our SCUBA-2 observations with various COBE EBL measurements.The low angular resolution satellite EBL measurements and resolved source counts contain different sets of systematics.Comparing them provides insight into the nature of the sources that give rise to the EBL, as well as testing them against each other.
The 450 µm EBL estimations based on the COBE FIRAS experiment are 109 Jy deg −2 (Puget et al. 1996), 142 Jy deg −2 (Fixsen et al. 1998), and 150 Jy deg −2 (Gispert, Lagache, & Puget 2000).The uncertainty in them is about 30%.The difference between these results arises from the subtractions of the foregrounds, and should be considered as a measure of systematic uncertainty.In the subsequent discussion, we compare our resolved EBL with the full range of 109-150 Jy deg −2 from the above COBE results.
We show the 450 µm surface brightness integrated from our source counts in Fig. 8.For this calculation, we adopted our Schechter fit for above 10 mJy and our power-law fit for 0.1-10 mJy.For uncertainty estimation, we perturbed our counts according to their errors, re-conducted the fitting, and calculated the dispersion in the results.Down to 3.5 mJy, which is the corrected flux density of our faintest 4-σ source, we have resolved 35.5 ± 4.3 Jy deg −2 of the 450 µm EBL (black solid curve in Fig. 8).This corresponds to 24 ± 3% to 33 ± 4% of the total EBL, depending on which EBL estimate we adopt.If we integrate down to 1 mJy, where we have constraints on the counts from our fluctuation analysis, we obtain a value of 90.0 ± 17.2 Jy deg −2 (black dashed curve in Fig. 8), corresponding to 60 ± 11% to 83 +15 −16 % of the EBL.If we continue the integration with the same power-law slope, we will reach 100% resolution of the EBL from Puget et al. (1996) at 0.77 mJy, while in order to fully account for EBL from Gispert, Lagache, & Puget (2000), the counts have to extend to approximately 0.48 mJy with the same faint-end slope.Roughly speaking, the flux of 0.48 mJy corresponds to infrared luminosities of 3.8×10 10 L ⊙ and 7.9×10 10 L ⊙ at z = 1 and z = 2, respectively, if we assume the luminosity-dependent dust SEDs of Chu et al. (in prep.), which incorporate the latest WISE and Herschel photometry for local infrared selected galaxies (Chu et al. 2017).The luminosity would be 30% (z = 2) to 40% (z = 1) higher if we assume the median SED of bright 870 µm-selected SMGs in Danielson et al. (2017), because such luminous SMGs have higher dust temperature and therefore more emission from 200 µm.The luminosity of 3.8-7.9× 10 10 L ⊙ corresponds to SFR of roughly 6-13 M ⊙ yr −1 .These are in the regime of the typical SFR of rest-frame UV selected high-redshift galaxies.The black curve is the result derived from our 450 µm counts, and the light shaded area is its uncertainty range.The colored curves are results from previous 450 µm surveys.The solid portions of the curves represent counts constrained by detected sources, while the dotted portions represent the faint-end extrapolations of the counts.The dashed portion of the black curve is where our counts are constrained through the fluctuation analysis.The models of Béthermin et al. (2012b) and Lacey et al. (2016) are also shown for comparison.The right-hand y-axis is the fraction of the resolved EBL, and the 100% resolution level is set to be the mid point between the minimum and maximum values among the COBE determinations.
The above assumes that the extension of the counts below 1 mJy maintains a power-law slope of α f = −2.6.A more likely scenario is that the counts become shallower below 1 mJy.Shallower (Schechter) slopes in the UV luminosity functions are typical for faint rest-frame UV selected populations, ranging from −0.8 to −1.7 at z ≃ 1.7 to 3 (Steidel et al. 1999;Sawicki & Thompson 2006;Reddy & Steidel 2009).Steeper slopes are observed at higher redshifts: −1.6 at z ≃ 4 (Bouwens et al. 2015), and −2.0 at z ≃ 6-8 (Bouwens et al. 2015(Bouwens et al. , 2017)), but are still shallower than the slope we observed on 450 µm sources.
Another hint of a shallower slope in the extreme faint end comes from ALMA imaging at longer wavelengths.Recent ALMA observations reported nearly full resolution of the 1.1 mm EBL with serendipitous detections of faint sources (Fujimoto et al. 2016) and stacking analyses (Aravena et al. 2016;Wang, Kohno, et al. 2016;Dunlop et al. 2017).Roughly speaking, it requires the detections of S 1100 ≃ 0.01-0.02mJy sources to achieve a full resolution of the 1.1 mm EBL.For z ≃ 1-2, these sources would have 450 µm fluxes of about 0.1-0.2mJy, a few times fainter than the above 0.48 mJy extrapolated based on a constant slope.This again suggests a slope shallower than −2.6.
Our STUDIES survey was allocated 330 hrs of observations (including some overheads).Roughly 40% of the observations were carried out in the first year and included in this paper.Once the survey is complete, we can expect to detect sources with intrinsic flux densities of 2 mJy in the deepest region and directly resolve around 55 Jy deg −2 of the EBL (30% to 50%).Furthermore, an extension of STUDIES was recently approved to conduct an equally deep pointing in the Subaru/XMM-Newton Deep Survey field (Furusawa et al. 2008).Both the increases in depth and survey area will lead to dramatically improved source counts in the 2-30 mJy bright end.The deeper image may enable meaningful fluctuation analysis below 2 mJy with more model parameters than our current simple two-parameter model.If this is the case, better constraints in the 1-2 mJy regime may help to narrow down the parameter space for the extremely faint end below 1 mJy (e.g., a shallower slope and/or a turn-over flux density).It will be interesting if the improved fluctuation analysis can reach > 100% resolution of the COBE EBL, and hence start to pin down the EBL better than the satellite estimates.The most robust EBL estimate will likely require extremely deep and wide interferometric imaging in the future.

SUMMARY
In our JCMT Large Program, STUDIES, we are carrying out extremely deep 450 µm imaging with SCUBA-2 in the COSMOS-CANDELS region.The 7 ′′ resolution of the 450 µm band provides the opportunity to probe beyond the conventional 850 µm confusion limit to detect fainter galaxies.Our goal is to detect faint 450 µm sources close to the confusion limit of SCUBA-2, to study a representative sample of the high-redshift FIR galaxy populations that give rise to the bulk of the FIR background and the cosmic star formation.With the first year of STUDIES data, we reached a noise of 0.91 mJy at the map center for point sources, about 30% deeper than previous deepest 450 µm map, and covered a deep area of R = 7. ′ 5. We detected 98 and 141 sources at 4.0 and 3.5 σ, respectively.Our source counts are best constrained between 3.5 and 25 mJy (4 σ), and are consistent with most of the previous counts derived from blank-field and lensing-cluster surveys.The field-tofield variance among the various surveys is about 20%, much smaller than the variance observed in 850 µm surveys, perhaps suggesting weaker clustering in the 450 µm population at scales larger than our field size.In this flux density range, our counts are consistent with a power law with a slope of α = −2.59± 0.10.We further constrain the counts at 1-3.5 mJy with a fluctuation analysis.We see evidence of a termination or turn-over of the faint-end counts between 1 and 3.5 mJy.The power-law slope at 1-3.5 mJy is α f = −2.6 +0.4  −0.7 , consistent with the slope at > 3.5 mJy.This is also consistent with the counts at around 1 mJy derived from previous lensing-cluster surveys.On the other hand, ours and all other SCUBA-2 450 µm counts appear significantly lower than Herschel counts at 350 and 500 µm.This discrepancy is likely caused by source blending under the coarse Herschel beam, amplified by source clustering at the scales of the beam.Because of its higher angular resolution, SCUBA-2 counts at 450 µm do not suffer from these effects.Our extremely deep SCUBA-2 map has resolved a substantial fraction of the 450 µm EBL estimated by COBE, 24 ± 3% to 33 ± 4% from the 4-σ sources at > 3.5 mJy, and 60 ± 11% to 83 +15 −16 % if we include the 1-3.5 mJy faint-end counts derived from the fluctuation analysis.STUDIES is an ongoing survey and we expect that our future deeper image can be used to better determine the EBL at 450 µm, as well as providing accurate counts for constraining galaxy evolution models.2015b).Later we will show that source blending is not entirely negligible in our deep image.For the other effects, previous studies estimated the amplitude of the corrections through simulations similar to ours, and applied these corrections to the observed flux densities and the counts.The correction to the counts typically looks like where F spu is the fraction of spurious sources, and F comp is the detection completeness fraction.To derive the two F terms and the conversion between S obs and S corr , one normally has to assume a one-to-one relation between an output source and an input source in the simulation.The possibility of source blending makes this assumption no longer valid.We avoid this and do not attempt to correct the observed counts using Eq.A1.Instead, we fully rely on the iterative procedure and Eq. 1 described in Section 3.2 to derive the intrinsic counts.However, we can still estimate these effects using our simulations, in order to obtain a picture of our observation and source extraction efficiency, and to compare with previous studies.
To estimate flux boosting caused by the Eddington bias and faint confusing sources, we matched sources in the input and output catalogs of our simulations.For each output source, we searched for input sources within a radius of a beam HWHM.We consider that we have a match when the flux densities of the input and output sources are within a factor of two of each other.This ensures that the input and output sources have similar flux densities, since for any output source, it would be possible to associate it with an arbitrarily faint nearby input source, given that faint sources have much higher spatial density (we will come back and discuss the choice of this factor of two at the end of this section).When there are multiple input sources meeting the above distance and flux ratio criteria, the brightest one is considered as the match.We started from the brightest output source and moved down the list.The matched input sources were not considered in the subsequent searches.We repeated this for the 200 realizations.The output-to-input flux ratio of each matched pair is the flux boosting factor.Fig 9 shows the flux boosting factor as a function of the output (observed) source flux density (panel a) and the input (intrinsic) source flux density (panel b).Overall, the required correction for measured flux densities is around 20% in panel (a).Panel (b) show that sources as faint as about 2 mJy (intrinsic) can be strongly boosted into our 4-σ detection range at 4 mJy, but with very low detection completeness (panel c).The flux boosting factor in (b) is artificially capped at 2.0 for the faintest sources, because we require the flux densities of the input and output sources to be within this factor.However, the curves do not suggest a much higher boosting factor if we remove this requirement.
To estimate the completeness and spurious fraction, we only considered input sources that have a chance of being detected, given the noise levels at their locations.This is necessary because our map has a highly non-uniform sensitivity distribution.To perform this estimate, we used the +1σ flux boosting factor in Fig. 9 (a) to boost the flux densities in the input source list.Only those sources with boosted flux densities greater than 3.5 times and 4 times the local r.m.s.noise are considered, corresponding to 3.5 σ and 4 σ detection thresholds.We then matched the sources in the input and output catalogs using the same criteria as for the flux boosting estimate.An output source without a match is considered as spurious (although its flux may still be real, contributed by multiple fainter sources).An input source without a match is considered as undetected (although it may contribute flux to a nearby source).The detection completeness calculated this way is presented in Fig. 9 (c),  d) and (e), the black and yellow curves are derived by requiring the matched input and output sources to have flux density ratios less than 2 and 4, respectively.The difference between the yellow and black curves shows that the required correction for spurious sources sensitively depends on the details of the matching between the input and output source lists.Panel (f) shows the fraction of observed sources that are blends of multiple sources.The three curves represent the cases where the flux densities of the blended sources in the input list are at least 1/2, 1/4, and 1/8 of the flux densities of the detected sources.The 1/4 curve best represents the multiple fraction to be found in sensitive interferometric followup, if the sources are not clustered at ∼ 10 ′′ scales.
for sources detected at > 4 σ and > 3.5 σ.The spurious rate is presented in Fig. 9 (d) and (e) (black curves), as functions of the output (measured) flux density and detection S/N, respectively.
In Fig. 9 (e), we see that the spurious fraction for ∼ 3.5 σ sources is approximately 50%.Even at ∼ 4 σ, the fraction is approximately 30%.Both values appear somewhat high.However, sources detected with such low S/N only comprise a small fraction in our sample.Fig. 9 (d) better reflects the overall spurious fraction, as here sources detected at different S/N levels are mixed according to the realistic counts and sensitivity distribution of the map.In Fig. 9 (d), there is a small bump in spurious fractions between 10 and 20 mJy.This is because some of these brighter sources are detected in the outer area where noise increases rapidly.For > 4 σ detections, the spurious rate is below 15% above 7 mJy and below about 20% in nearly the entire flux density range of interest.If we lower the threshold to 3.5 σ, the completeness would increase by almost a factor of two at the faint end of 3-4 mJy, but the spurious rate would also increase to 35%.Given this, it is not obvious to us whether we should favor 3.5-σ detections over 4-σ detections for number count estimates.Furthermore, because of the nature of our Eq. 1, the choice of detection threshold will not change the corrected source counts; it will only slightly change the interpretation of the intrinsic flux density range over which the counts are measured.
We compare the above results with other SCUBA-2 450 µm surveys and do not find major discrepancies.Our flux boosting factor of 20% is comparable to those in Casey et al. (2013) and Chen et al. (2013a), but seems much higher than that in Zavala et al. (2017, their figure 3).The very low flux boosting factor in Zavala et al. might be caused by the difference in how the simulations were conducted.The behavior of the completeness curves in Fig. 9 (c) is qualitatively similar to those presented by other authors (Geach et al. 2013;Casey et al. 2013;Zavala et al. 2017).The differences in sensitivities and mapping strategies prevent further quantitative comparison.Our spurious rates of 20-40% near the detection limit is also comparable to that in Casey et al. (2013).
Finally, to obtain a rough idea of source blending, we estimated the fraction of output sources that consist of multiple blended faint sources.We again conducted source matching with a search radius of half a beam FWHM.However, this time, for each output source, we consider input sources that are brighter than 1/2, 1/4, and 1/8 times the output flux density.Fig. 9 (f) shows the fraction of output sources that contain contributions from more than one input source that meet the above flux ratio thresholds.The "> 1/2" curve (solid) represents cases where sources of nearly equal brightness are blended and detected; these are rare, and the blending fraction does not exceed 10%.The "> 1/4" curve (dashed) probably better mimics the multiple fractions to be found in sensitive interferometric surveys, since the aimed detection S/N for the single-dish fluxes is typically between 10 and 20 and a four times fainter companion can be detected at > 3 σ.This suggests that we will find roughly 10-20% of our sources to be multiples if we conduct deep ALMA Band-9 imaging followup.The majority of the multiple source are pairs, and blends of more than three sources only account for 6% of them.This predicted multiple fraction is much lower than the 35%-60% multiple fractions typically found in interferometric followup of 850 µm SCUBA-2 and LABOCA sources (Barger et al. 2012;Hodge et al. 2013;Simpson et al. 2015a), This is a consequence of the smaller SCUBA-2 beam at 450 µm.Nevertheless, the fraction of 10% to 20% is comparable to the error in the amplitudes of our counts, meaning that its effect is not negligible below 10 mJy.The "> 1/8" curve (dotted) steeply raises from above 20 mJy to below 10 mJy, but artificially saturates at 6 mJy because we do not include < 1 mJy sources in our simulations.
We point out that the results presented in this section are highly sensitive to the details of the source matching criteria between the input and output lists.For example, if we allow the flux density ratio between the matched input and output sources to be greater than two (i.e., make easier to find a match), we can artificially decrease the spurious fraction, increase the flux boosting factor, and extend the derived intrinsic counts to fainter flux density limits.To demonstrate this, we show the spurious fraction derived with input and output flux ratios relaxed to 4 in Fig. 9 (d) and (e) with yellow curves.For 4-σ detections, the spurious fraction is below 5% over nearly the entire flux density range, but many of the added detections have flux densities that are dramatically boosted (more than a factor of two).We stress again that we do not base our source count correction on the corrections of flux boosting, and completeness and spurious sources, because of all the above ambiguities and the arbitrariness in the choice of source matching criteria.

B. EFFECT OF CLUSTERING
The results in the previous section are a strong indication that our counts will be biased if we do not properly take source blending into account when we correct the observed counts using Eq.A1.Our Eq. 1 does not suffer from this complication.On the other hand, our simulations (as well as similar simulations in the literature) assume a random (unclustered) spatial distribution.Our counts will still be biased if the clustering of the 450 µm sources at the 10 ′′ scale (20 comoving kpc) is strong enough to alter the blending fraction.The simulations of Cowley et al. (2015) and Béthermin et al. (2017) show that clustering is not likely to have a strong effect on the observed counts for our beam size.The ALMA observations of 870 µm single-dish selected sources of Hodge et al. (2013) seem to also support weak or no clustering at such small scales.Future ALMA observations of 450 µm sources are needed to further test the small-scale clustering in the 450 µm population.In the mean while, we use the simulations of Béthermin et al. (2017), which includes clustering effects, to test whether clustering would bias our results.
We used the 2 deg 2 lightcone provided by Béthermin et al. (2017) and the associated 450 µm fluxes of the simulated sources.The general procedure is identical to our simulations described in Section 3.2, except that we use sources in Béthermin et al. (2017) for the input.We convolved the simulated sources with our SCUBA-2 PSF, and randomly pick a location in the 2 deg 2 region to simulate an observed image.This is repeated 100 times each for the positive and negative true noise map.We performed source extraction over the 200 simulated maps and derived the raw (observed) source counts.The results represent the observed counts where sources are clustered, and are shown with red symbols and curve in Fig. 10.We then randomized the positions of the simulated sources in the 2 deg 2 region, and repeated the procedure.The results represent the observed counts where sources are not clustered, but with identical intrinsic counts, and are shown with blue symbols and curve in Fig. 10.With the fitted Schechter functions (dashed curves), we find that at 5, 10, and 20 mJy, the counts with clustering are higher than the unclustered cases by −3.2%, +0.6%, and +3.6%, respectively.These are essentially negligible, and do not suggest a systematic bias caused by source clustering.By observing the data points in Fig. 10, it is apparent that the accuracy of this test is limited by the size of the lightcone as well as the area coverage of our observations.We thus confirm that our counts are not biased by clustering, thanks to the small beam size of SCUBA-2 at 450 µm.
Next, we test whether clustering would bias the fluctuation analyses for the faint end.We again use the simulations of Béthermin et al. (2017), with the original simulated source positions for the clustered case and randomized source positions for the unclustered case.The procedure is identical to that described in Section 4. For simplicity, we tested the cases with faint-end cutoffs of S min = 0.5 and 1.0 mJy.We created 400 simulated maps for each case, and calculate the log-likelihood against the observed residual map.The pixel brightness distributions of the simulated residual maps are shown in Fig. 11.For S min = 1.0 mJy, the dispersions of the log-likelihoods of the 400 simulated maps are 104 (clustered) and 73 (unclustered), while the mean difference between the clustered and unclustered log-likelihoods is 37.For S min = 0.5 mJy, the dispersions of the log-likelihoods are 109 (clustered) and 89 (unclustered), while the mean difference between the clustered and unclustered log-likelihoods is 25.In both cases, the difference induced by source clustering is much smaller than the variance in the ensemble of simulated images.This can also be seen in Fig. 11, where the clustered and unclustered curves essentially overlap, and the difference caused by clustering is much smaller than changing the cutoff flux.Based on these, we conclude that our fluctuation analyses are not likely to be biased by source clustering at the scale of the SCUBA-2 450 µm beam.2017) and the true noise map (red).For both S min = 0.5 and 1.0, the distributions for the clustered and unclustered cases are nearly identical.The differences in the log-likelihood are also negligible.

Figure 1 .
Figure 1.STUDIES 450 µm S/N map.Only the deeper R = 7. ′ 5 is shown here.The map center is at R.A. = 10:00:30.2,Dec. = +02:26:50.The contours represent noise levels from 1 mJy with a multiplicative step of √ 2. Sources detected at 4 σ are marked with R = 12 ′′ circles; there are 98 such sources.The easternmost source is not included in our number counts because of the requirement on noise of the source flux.Several of them have nearby companions with distances of ∼ 10 ′′ , identified by our CLEAN-like source extraction algorithm (Section 3.1).There are an additional 43 sources detected at 3.5-4 σ.These fainter sources are not labeled here to avoid confusion.

Figure 2 .
Figure2.Distribution of pixel brightness in the 450 µm image (black) and in the "true noise" map (red, see Section 3.2).The excesses in the both positive and negative sides in the image are produced by celestial objects and the negative bowling in the PSF.The noise histogram is a sum of Gaussians of various widths, because of the non-uniform sensitivity distribution.

Figure 3 .
Figure3.Cumulative map area versus noise level.The noise distribution in our image is roughly axisymmetric, so area can also be approximately mapped to radius, which is shown with the right-hand y-axis.The map center has a noise of 0.91 mJy beam −1 .The noise slowly increases to ∼ 2 mJy beam −1 , where the map area reaches 100 arcmin 2 .After that, the noise increases more rapidly toward the outer region.For our number counts, we do not use image area that is more than five times shallower than at the map center, indicated by the vertical dashed line.

Figure 4 .
Figure 4. Differential 450 µm counts.Open squares are the observed raw counts and solid squares are the corrected counts.The raw and corrected counts have different fluxes because of the correction for flux boosting (SectionA).The solid black curve is the input Schechter counts in the simulations.The dashed black curve is the measured output counts in the simulations, and can also be described by a Schechter function.The two dotted black curves show the 68th-percentile range of the outputs of the 200 realizations.The 68th-percentile range reasonably matches the observed error bars derived based on the assumption of Poisson errors.Two parameterizations to the corrected counts are shown with colored curves.The red solid curve is a Schechter fit to the corrected counts.The reduced χ 2 of this fit is 0.45 (N df = 6).This curve is almost indistinguishable from the input Schechter function (black solid curve), showing that our iterative procedure converges.The green dashed line is a power-law fit to the corrected counts after excluding the brightest flux density bin.The reduced χ 2 of this fit is 0.38 (N df = 6).
Figure5.Distributions of pixel brightness in the residual maps from the observations (black), the true noise map (red), and three representative models from our fluctuation analysis that have high likelihoods (orange, green, and blue, corresponding to the three crosses in Fig.6).Here we only show the pixels within 5 ′ of the map center, i.e., the pixels included in our fluctuation analysis.
survey covers a large area of ∼ 400 arcmin 2 but with a shallower depth of 4.1 mJy.TheGeach et al. and Casey et al. surveys  and ours are in the COSMOS field.The Casey et al. survey area fully encloses ours, as well as that of Geach et al., while our survey area partially overlaps with that of Geach et al.The Zavala et al. survey field is in the Extended Groth Strip, a different line of sight.Both the surveys of Chen et al. (2013b) and Hsu et al. (2016) include lensing cluster fields of various depth and the COSMOS field for their 450 µm analyses.For COSMOS, both teams used the data taken by Casey et al.The three cluster fields in Chen et al. with 450 µm data are all included in Hsu et al., who added additional integrations, and also two further clusters.The sensitivities of the lensing-cluster surveys are determined by the amplification factors of the lensed sources, rather than the instrumental noise.Because of the overlapping of data and survey fields, the results in these papers are not completely independent of each other, except for those in Zavala et al.

Figure 8 .
Figure8.Integrated 450 µm surface brightness from various SCUBA-2 surveys.The dark shaded area shows the range of 450 µm EBL values inferred from various COBE studies.The black curve is the result derived from our 450 µm counts, and the light shaded area is its uncertainty range.The colored curves are results from previous 450 µm surveys.The solid portions of the curves represent counts constrained by detected sources, while the dotted portions represent the faint-end extrapolations of the counts.The dashed portion of the black curve is where our counts are constrained through the fluctuation analysis.The models ofBéthermin et al. (2012b) andLacey et al. (2016) are also shown for comparison.The right-hand y-axis is the fraction of the resolved EBL, and the 100% resolution level is set to be the mid point between the minimum and maximum values among the COBE determinations.

Figure 9 .
Figure 9. Bias effects in the observations and source extraction.Panel (a) shows the flux boosting factor versus output (observed) source flux density.Panel (b) shows the flux boosting factor versus input (intrinsic) source flux density.In both (a) and (b), the solid curves are mean values, while the dashed curves are the 68th (±1 σ) and 50th-percentile (median) values.The overall boosting correction for measured flux densities is around 20%.A small number of intrinsically faint sources (S input 2 mJy) can be strongly boosted into our detection range.Panel (c) shows the completeness versus input (intrinsic) source flux density for a 4-σ detection threshold (solid curve) and a 3.5-σ detection threshold (dashed curve).Panel (d) shows the spurious fraction versus output (measured) flux density.The spurious fraction for 4-σ detections are shown with the solid curve, and for 3.5-σ detections with the dashed curve.Panel (e) shows the spurious fraction versus detection S/N.In both (d) and (e), the black and yellow curves are derived by requiring the matched input and output sources to have flux density ratios less than 2 and 4, respectively.The difference between the yellow and black curves shows that the required correction for spurious sources sensitively depends on the details of the matching between the input and output source lists.Panel (f) shows the fraction of observed sources that are blends of multiple sources.The three curves represent the cases where the flux densities of the blended sources in the input list are at least 1/2, 1/4, and 1/8 of the flux densities of the detected sources.The 1/4 curve best represents the multiple fraction to be found in sensitive interferometric followup, if the sources are not clustered at ∼ 10 ′′ scales.

Figure 10 .
Figure10.Raw SCUBA-2 450 µm counts derived with the simulations ofBéthermin et al. (2017).The red symbols are the results with source clustering and the blue symbols without clustering.The dashed curves are Schechter fits to the counts.

Figure 11 .
Figure11.Distributions of pixel brightness in the residual maps from the simulations using the lightcone ofBéthermin et al. (2017) and the true noise map (red).For both S min = 0.5 and 1.0, the distributions for the clustered and unclustered cases are nearly identical.The differences in the log-likelihood are also negligible.

Table 2 .
Parameterizations for the Corrected Counts