The SCUBA-2 Cosmology Legacy Survey: 850 μm maps, catalogues and number counts

Abstract

We present a catalogue of ∼3000 submillimetre sources detected (≥3.5σ) at 850 μm over ∼5 deg² surveyed as part of the James Clerk Maxwell Telescope (JCMT) SCUBA-2 Cosmology Legacy Survey (S2CLS). This is the largest survey of its kind at 850 μm, increasing the sample size of 850 μm selected submillimetre galaxies by an order of magnitude. The wide 850 μm survey component of S2CLS covers the extragalactic fields: UKIDSS-UDS, COSMOS, Akari-NEP, Extended Groth Strip, Lockman Hole North, SSA22 and GOODS-North. The average 1σ depth of S2CLS is 1.2 mJy beam⁻¹, approaching the SCUBA-2 850 μm confusion limit, which we determine to be σ_c ≈ 0.8 mJy beam⁻¹. We measure the 850 μm number counts, reducing the Poisson errors on the differential counts to approximately 4 per cent at S₈₅₀ ≈ 3 mJy. With several independent fields, we investigate field-to-field variance, finding that the number counts on 0.5°–1° scales are generally within 50 per cent of the S2CLS mean for S₈₅₀ > 3 mJy, with scatter consistent with the Poisson and estimated cosmic variance uncertainties, although there is a marginal (2σ) density enhancement in GOODS-North. The observed counts are in reasonable agreement with recent phenomenological and semi-analytic models, although determining the shape of the faint-end slope (S₈₅₀ < 3 mJy) remains a key test. The large solid angle of S2CLS allows us to measure the bright-end counts: at S₈₅₀ > 10 mJy there are approximately 10 sources per square degree, and we detect the distinctive up-turn in the number counts indicative of the detection of local sources of 850 μm emission, and strongly lensed high-redshift galaxies. All calibrated maps and the catalogue are made publicly available.

catalogues, surveys, galaxies: evolution, galaxies: high-redshift, cosmology: observations

1 INTRODUCTION

Nearly a quarter of a century has passed since it was predicted that submillimetre observations could provide important insights into the nature of galaxies in the early Universe beyond the reach of optical and near-infrared surveys (Blain & Longair 1993). If early star-forming galaxies contained dust, then ultraviolet photons should be reprocessed through the far-infrared (Hildebrand 1983) and redshifted into the submillimetre. Early observations certainly showed that some high-redshift sources are emitting a large fraction of their bolometric emission in the rest-frame far-infrared, detectable in the submillimetre, with integrated luminosities comparable to or exceeding local ultraluminous (ULIRG, 10¹² L_⊙) infrared galaxies (Rowan-Robinson et al. 1991; Clements et al. 1992). We now know that the far-infrared background (FIRB; Puget et al. 1996; Fixsen et al. 1998; Lagache et al. 1998) represents about half of the energy density associated with star formation integrated over the history of the Universe (Dole et al. 2006) and the peak of the volume averaged star formation rate density (SFRD) occurred at z ∼ 1–3, to which submillimetre sources are expected to contribute significantly (Devlin et al. 2009). Identifying and characterizing the galaxies contributing to the FIRB was (and remains) a major goal, and motivates blank-field submillimetre surveys.

About two decades ago, the first submillimetre maps of the high-redshift Universe were made (Smail et al. 1997; Barger et al. 1998; Hughes et al. 1998; Lilly et al. 1999), opening a new window on to early galaxies. With 20 yr of follow-up work across the electromagnetic spectrum, we now have a good grasp of the nature of ‘Submillimetre Galaxies’ (SMGs) and their cosmological significance.¹ Nevertheless, the picture is far from complete. SMGs selected at 850 μm (in single-dish surveys) with flux densities above a few mJy² lie at 〈z〉 ≈ 2–3 (e.g. Chapman et al. 2005; Pope et al. 2005; Wardlow et al. 2011; Koprowski et al. 2014; Simpson et al. 2014), are massive (Swinbank et al. 2004; Hainline et al. 2011; Michalowski et al. 2012), gas-rich (Greve et al. 2005; Tacconi et al. 2006, 2008; Carilli et al. 2010; Engel et al. 2010; Bothwell et al. 2013) and are associated with large supermassive black holes (Alexander et al. 2005, 2008; Wang et al. 2013). These properties make SMGs the obvious candidates for the progenitor population of massive elliptical galaxies today, seen at a time of rapid assembly a few billion years after the big bang (Lilly et al. 1999; Genzel et al. 2003; Swinbank et al. 2006), with star formation rates in the range 100–1000 M_⊙ yr⁻¹ derived from their integrated infrared luminosities (e.g. Magnelli et al. 2012; Swinbank et al. 2014).

The formation mechanism of SMGs remains in debate: by analogy with local ULIRGs, which are almost exclusively merging systems, it is predicted that SMGs form during major mergers of gas-dominated discs (Baugh et al. 2005; Ivison et al. 2012), triggering star formation and central black hole growth. There is certainly observational evidence to support this, perhaps most convincingly in morphology and gas kinematics (e.g. Swinbank et al. 2010; Tacconi et al. 2010; Alaghband-Zadeh et al. 2012; Chen et al. 2015). On the other hand, hydrodynamic simulations may be able to reproduce the properties of SMGs without the need for mergers; for example, if there is a prolonged (∼1 Gyr) phase of gas accretion which drives high star formation rates, where cooling is accelerated through metal enrichment at early times (e.g. Narayanan et al. 2015, see also Davé et al. 2010). In recent semi-analytic models, starbursts triggered by bar instabilities in galaxy discs are the dominant mechanism producing SMGs in model universes (Lacey et al. 2016), and indeed there is some empirical evidence that SMGs have optical/near-infrared morphologies consistent with discs (e.g. Targett et al. 2013).

Observations in the 850 μm atmospheric window offer a unique probe of the distant Universe, owing to the so-called negative k-correction (Blain & Longair 1993). For cosmological sources, the 850 μm band probes the Rayleigh–Jeans tail of the cold dust continuum emission of carbonaceous and silicate grains in thermal equilibrium in the stellar ultraviolet radiation field. As the thermal spectrum is redshifted, cosmological dimming is compensated for by increasing power as one ‘climbs’ the Rayleigh–Jeans tail as it is redshifted through the band. Thus, two sources of equal luminosity will be observed with roughly the same flux density at 850 μm at z ≈ 0.5 and z ≈ 10. As a guide, a galaxy in the ultraluminous class (with L_IR ≈ 10¹² L_⊙) is observed with a flux density of 1–2 mJy at 850 μm over most of cosmic history (Blain et al. 2002). For this reason, flux-limited surveys at 850 μm offer the opportunity to sample huge cosmic volumes, potentially probing well into the epoch of re-ionization.

Despite the large redshift depth probed by deep 850 μm surveys, the solid angle subtended by existing surveys, and their sensitivity, has been bounded by technology: until recently, submillimetre cameras have been limited in field of view and sensitivity that has made degree-scale mapping difficult. However, submillimetre imaging technology has blossomed over the past 20 yr. At first, only single channel broad-band submillimetre photometers were available in (e.g. Duncan et al. 1990), making survey work impossible. Then the first cameras came online, mounted on 10–15 m single-dish telescopes such as the Caltech Submillimeter Observatory (CSO) and the James Clerk Maxwell Telescope (JCMT): the Submillimeter High Angular Resolution Camera (SHARC; Wang et al. 1996) and the Submillimetre Common-User Bolometer Array (SCUBA; Holland et al. 1999) using small arrays of tens of bolometers covering just a few arcminutes field of view. These arrays enabled the first extragalactic submillimetre surveys (Smail et al. 1997; Hughes et al. 1998), but covering a cosmologically representative solid angle at the necessary depth was still tremendously expensive in terms of observing time.

Further cameras based on bolometer arrays were developed through the late 1990s: Bolocam (Glenn et al. 1998), MAMBO (Kreysa et al. 1998), SHARC-II (Dowell et al. 2002), LABOCA (Siringo et al. 2009) and AzTEC (Wilson et al. 2008) and the scale of extragalactic submillimetre surveys grew in tandem (e.g. Eales et al. 2000; Scott et al. 2002; Borys et al. 2003; Webb et al. 2003; Greve et al. 2004; Coppin et al. 2006; Scott, Dunlop & Serjeant 2006; Weiß et al. 2009; Austermann et al. 2010; Scott et al. 2010, 2012). Unfortunately, the semiconductor technology underlying the first and second generation of submillimetre cameras is not scalable, limiting bolometer arrays to around 100 pixels. A solution was found in superconducting transition edge sensors (TES; see Irwin & Hilton 1995) coupled with superconducting quantum interference device (SQUID) amplifiers that allowed for the construction of submillimetre sensitive bolometer arrays an order of magnitude larger than previously achieved. Clearly, this opened up the possibility of performing much larger, more efficient submillimetre surveys than had ever been possible before from the ground.

The second-generation SCUBA camera, SCUBA-2, on the JCMT is the first of such large format instruments using TES technology (Holland et al. 2013). SCUBA-2 comprises two arrays (for the 450 μm and 850 μm bands) of 5120 bolometers each, covering an 8 arcmin field of view. With mapping speeds (to equivalent depth) over an order of magnitude faster than its predecessor, SCUBA-2 has enabled a huge leap in submillimetre survey science. TES focal plane arrays have also formed the basis of other recent submillimetre instrumentation, such as the South Pole Telescope (Carlstrom et al. 2011) and Atacama Cosmology Telescope (Swetz et al. 2011). Future large format submillimetre cameras are likely to make increasing use of Kinetic Inductance Detectors (KIDS; Day et al. 2003): the New Instrument of KIDS Arrays (NIKA2) on the 30 m Institut de Radioastronomie Millimétrique (IRAM) telescope (Monfardini et al. 2010) uses this new detector technology.

Soon after commissioning of SCUBA-2, five JCMT ‘Legacy Surveys’ (JLS) commenced. The largest of these is the JCMT SCUBA-2 Cosmology Legacy Survey (S2CLS). In this paper, we present the wide 850 μm survey component of the S2CLS, presenting maps and a source catalogue for public use. This paper is organized as follows: in Section 2, we define the survey and describe data reduction and cataloguing procedures; in Section 3, we present the maps and catalogues, and in Section 4, we use these data to measure the number counts of 850 μm selected sources with the best statistical precision to date, including an analysis of the impact of cosmic variance on scales of ∼1°. We summarize the paper in Section 5. Where relevant, we adopt a fiducial ΛCDM cosmology with Ω_m = 0.3, |$\Omega _\Lambda =0.7$| and H₀ = 70 km s⁻¹ Mpc⁻¹.

2 THE SCUBA-2 COSMOLOGY LEGACY SURVEY

The S2CLS survey has two tiers: wide and deep. The wide tier covers several well-explored extragalactic survey fields: Akari-Northern Ecliptic Pole, COSMOS, Extended Groth Strip, GOODS-North, Lockman Hole North, SSA22 and UKIDSS-UDS (Fig. 1, Table 1), mapping at 850 μm during conditions where the zenith optical depth at 225 GHz was 0.05 < τ₂₂₅ ≤ 0.1 and field elevations exceeded 30°. In the deep tier, several deep ‘keyhole’ regions within the wide fields were mapped when τ₂₂₅ ≤ 0.05, conditions suitable for obtaining 450 μm maps which require the lowest opacities (Geach et al. 2013). Note that SCUBA-2 simultaneously records 450 μm and 850 μm photons, and while the complementary 450 μm data exist for the wide 850 μm maps we present here, they have not been processed, since they are not expected to be of sufficient quality given the observing conditions. In this paper, we present the maps (Fig. 1) and catalogue from the wide tier only.

Figure 1.

The JCMT SCUBA-2 Cosmology Legacy Survey: montage of signal-to-noise ratio maps indicating relative coverage in the seven extragalactic fields (see also Table 1). This survey has detected approximately 3000 submillimetre sources over approximately 5 deg². The two bright sources identified are ‘Orochi’, an extremely bright SMG first reported by Ikarashi et al. (2011) in UKIDSS-UDS, and NCG 6543 in Akari-NEP. For scale comparison, we show the 850 μm map of the UKIDSS-UDS from the SCUBA HAlf DEgree Survey (SHADES; Coppin et al. 2006) and the footprint of the Hubble Space Telescope WFPC2, corresponding to the size of the SCUBA map of the Hubble Deep Field from Hughes et al. (1998) – one of the first deep extragalactic maps at 850 μm. Note that the size of the primary beam of the Atacama Large Millimeter/submillimeter Array (ALMA) at 850 μm is comparable to the size of the JCMT beam: the full S2CLS survey subtends a solid angle over 100 000 times the ALMA primary beam at 850 μm. The angular scale of 30 arcmin subtends approximately 5 comoving Mpc at the typical redshift of the SMG population, z ≈ 2.

Table 1.

S2CLS survey fields (see also Fig. 1). Right Ascension and Declination refer to the central pointing (J2000). The area corresponds to map regions where the root mean squared instrumental noise is below 2 mJy. Note that at the end of the survey, the COSMOS field was only 50 per cent completed; remainder is now being observed to equivalent depth in a new survey (S2-COSMOS, PI: Smail; Simpson et al., in preparation).

Field name	R.A.	Dec.	Area	1σ 850 μm depth	Scan recipe	Astrometric reference
			(deg2)	(mJy beam−1)
Akari-North Ecliptic Pole	17 55 53	+66 35 58	0.60	1.2	45 arcmin pong	Takagi et al. (2012) 24 μm
COSMOS	10 00 30	+02 15 02	2.22	1.6	2×2 45 arcmin pong	Sanders et al. (2007) 3.6 μm
Extended Groth Strip	14 17 41	+52 32 15	0.32	1.2	6×1 15 arcmin pong	Barmby et al. (2008) 3.6 μm
GOODS-N	12 36 51	+62 12 52	0.07	1.1	15 arcmin pong	Spitzer-GOODS-N MIPS 24 μm cataloguea
Lockman Hole North	10 46 07	+59 01 17	0.28	1.1	30 arcmin pong	Surace et al. (2005) 3.6 μm
SSA22	22 17 36	+00 19 23	0.28	1.2	30 arcmin pong	Lehmer et al. (2009) 3.6 μm
UKIDSS-Ultra Deep Survey	02 17 49	−05 05 55	0.96	0.9	60 arcmin pong	UKIDSS-UDS Data Release 8 3.6 μmb

^airsa.ipac.caltech.edu/data/SPITZER/docs/spitzermission/observingprograms/legacy/goods

^bwww.nottingham.ac.uk/astronomy/UDS

2.1 Observations

The S2CLS was conducted for just over 3 yr, from 2011 December to 2015 February; Fig. 2 shows the time distribution of observations during the survey. The wide tier used the pong mapping strategy for large fields, whereby the array is slewed around the target (map centre) in a path which ‘bounces’ off the rectangular edge of the defined map area in a manner reminiscent of the classic arcade game (Thomas & Currie 2014). The pong pattern ensures that the array makes multiple passes back and forth between the map extremes, filling the square mapping area. To ensure uniform coverage the field is rotated 10–15 times (depending on map size) during an observation, resulting in a circular field with uniform sensitivity over the nominal mapping area (but with science-usable area beyond this, see Section 2.4.1). Scanning speeds were 280 arcsec s⁻¹ for maps of size 900 arcsec up to 600 arcsec s⁻¹ for the largest single map of 3300 arcsec. Observations were limited to 30–40 min each to monitor variations in observing conditions, with regular pointing calibrations performed throughout the night. Typical pointing corrections are of order ∼1 arcsec between observations. In addition to the zenithal opacity constraints described above, elevation constraints were also imposed: to ensure sufficiently low airmass, targets were only observed when above 30°, and a maximum elevation constraint of 70° was also imposed (only relevant for the COSMOS field). This high elevation constraint was set because it was found that the telescope could not keep pace with the alt-az demands of the scanning pattern, resulting in detrimental artefacts in the maps. Since the Lockman Hole North field is observable during COSMOS transit, the strategy was simply to switch targets as COSMOS rose above 70°.

Figure 2.

Time distribution of 850 μm observations. In total CLS conducted 2041 wide-field observations on 320 nights from 2011 November to 2015 February. The increase in frequency of observations towards the end of the survey reflects the effect of ‘extended observing’ into the post-sunrise morning hours when the opacity and conditions were still suitable for observations. Note that one observation is equivalent to 30–40 min of integration time.

For all but the EGS and COSMOS field, the targets were mapped with single pong scans with diameters ranging from 900 to 3300 arcsec (Table 1). The EGS was mapped using a chain of six 900 arcsec pong maps (each slightly overlapping) to optimize coverage of the multiwavelength data along the multiwavelength strip. In COSMOS, the mapping strategy was a mosaic consisting of a central 900 arcsec pong and four 2700 arcsec pong maps offset by 1147 arcsec in RA and Dec. from the central map, forming a 2 × 2 grid of ‘petals’ around the central pong, with some overlap. This was deemed preferable to obtaining a single very large pong map encompassing the full field, allowing depth to be built up in each tile sequentially. Only ∼50 per cent of the COSMOS area was completed to full depth, due to the end of JCMT operations by the original partners. The full 2° × 2° field is now being completed as part of a follow-on project ‘S2-COSMOS’ (PI: Smail and Simpson et al., in preparation). Fig. 1 shows a montage of the S2CLS fields to scale, and Fig. 3 shows an example of the sensitivity variation across a single pong map (the UKIDSS-UDS field), illustrating the homogeneity of the noise coverage across the bulk of the scan region, with instrumental noise varying by just ∼5 per cent across degree scales. We describe the process to create the S2CLS 850 μm maps in the following section.

Figure 3.

An example of the sensitivity coverage in a single S2CLS field. This map shows the instrumental noise map of the UKIDSS-UDS (a single pong), scaled between σ_instr = 0.8 and 1.2 mJy. Contours are at steps of 0.05 mJy starting at 0.8 mJy. This demonstrates the uniform nature of the pong map over the majority of the mapping region, radially rising beyond the nominal extent of the area scanned to uniform depth (effectively overscan regions receiving shorter integration time).

2.2 Data reduction

Each SCUBA-2 bolometer records a time-stream, where the signal is a contribution of background (mainly sky and ambient emission), astronomical signal and noise. The basic principle of the data reduction is to extract astronomical signal from these time-streams and map them on to a two-dimensional celestial projection. We have used the Dynamical Iterative Map-Maker (dimm) within the Sub-Millimetre Common User Reduction Facility (smurf; Chapin et al. 2013). We refer readers to Chapin et al. (2013) for a detailed overview of smurf, but describe the main steps, including specific parameters we have chosen for the reduction of the blank-field maps, here (see also Geach et al. 2013).

First, time-streams are downsampled to a rate matching the pixel scale of the final map, based on the scanning speed (Section 2.1). All S2CLS maps are projected on a tangential coordinate system with 2 arcsec pixels. Flat-fields are then applied to the time-streams using flat scans that bracket each observation, and a polynomial baseline fit is subtracted from each bolometer's time-stream (we actually use a linear – i.e. order 1 – fit). Then each time-stream is cleaned for spikes (using a 5σ threshold in a box size of 50 samples), DC steps are removed and gaps filled. After cleaning, the dimm enters an iterative process that aims to fit the data with a model comprising a common-mode fluctuating atmospheric signal, positive astronomical signal and instrumental and fine-scale atmospheric noise. The common mode modelling is performed independently for each SCUBA-2 sub-array, deriving a template for the average signal seen by all the bolometers. The common mode is then removed, and an extinction correction is applied (Dempsey et al. 2013). Next, a filtering step is performed in the Fourier domain, which rejects power at frequencies corresponding to angular scales θ > 150 arcsec and θ < 4 arcsec. The next step is to estimate the astronomical signal. This is done by gridding the time-streams on to the celestial projection; since each pixel will be sampled many times by independent bolometers (slewing over the sky in the pong scanning pattern), then the positive signal in a given pixel can be taken to be an accurate estimate of the astronomical signal (assuming the previous steps have eliminated all other sources of emission or spikes, etc.). This model of the astronomical signal is then projected back to a time-stream and subtracted from the data. Finally, a noise model is estimated for each bolometer by measuring the residual, which is then used to weight the data during the mapping process in additional steps. The iterative process above runs until convergence is met. In this case, we execute a maximum of 20 iterations, or terminate the process when the map tolerance Δχ² reaches 0.05.

S2CLS obtained many individual scans of each field. The dimm allows for all the scans to be simultaneously reduced in the manner described above. However, we adopt an approach where the dimm is only given individual observations, producing a set of maps for each target field which can then be co-added into a final stack. For this, we use the picard recipe mosaic_jcmt_images which uses the wcsmosaic task within the starlinkkappa package, weighting each input image by the inverse variance per pixel. With a set of individual observations for each field, we can also construct maps of sub-sets of the data and produce jackknife maps where a random 50 per cent of the images are inverted, thus removing astronomical signal in the final stack, and generating source-free noise realizations of each field (Fig. 4); useful for certain statistical tests.

Figure 4.

Distribution of pixel values in the UKIDSS-UDS flux density map, showing the characteristic tail representing astronomical emission. The shaded region shows the equivalent distribution in a jackknife map, constructed by inverting a random half of the data before co-addition. The dashed line is simply a normal distribution with zero mean and scale set to the standard deviation of pixel values in the jackknife map, illustrating that the noise in the map is approximately Gaussian.

The last processing step is to apply a matched filter to the maps, convolving with the instrumental PSF to optimize the detection of point sources. We use the picard recipe scuba2_matched_filter which first smooths the map (and the PSF) with a 30 arcsec Gaussian kernel, then subtracts this from both to remove any large-scale structure not eliminated in the filtering steps that occurred during the dimm reduction. The choice of a 30 arcsec kernel has not been optimized; however, this scale proved to be effective at eliminating any remaining large-scale structure from examination of the maps before and after the match-filtering step. In addition, since we are concerned solely with the detection of point sources, rejecting emission structure on scales will have a negligible impact on the detection of sources, whilst ensuring a uniform background across the map. After this subtraction step, the map is then convolved with the smoothed beam; a step that optimizes the detection of emission features matching the beam (i.e. point sources). A flux conversion factor of 591 Jy beam⁻¹ pW⁻¹ is applied to give the maps units of flux density. This canonical calibration is the average value derived from observations of hundreds of standard submillimetre calibrators observed during the S2CLS campaign (Dempsey et al. 2013). The filtering steps employed in the data reduction, including the match-filtering step, introduce a slight (10 per cent) loss of response to point sources. We have measured this loss by injecting a model source of known (bright) flux density into the data and recovering its flux after filtering; we correct for this in the flux calibration. The absolute flux calibration is expected to be accurate to within 15 per cent.

2.3 Astrometric refinement and registration

The JCMT pointing is regularly checked against standard calibrators during observations, with typical pointing drift corrections typically of order 1–2 arcsec; similar to the pixel scale at which the maps are gridded. To improve the astrometric refinement of the final co-added maps, we adopt a maximal signal-to-noise stacking technique: for each field, we use a mid-infrared selected catalogue and stack the submillimetre maps at the positions of reference sources to measure a high-significance statistical detection. We repeat the process many times, updating the world coordinate system reference pixel coordinates at each step with small Δα and Δδ increments. The goal is to find the (Δα, Δδ) that maximize the signal to noise of the stack in the central pixel. We iterate over several levels of refinement until no further change in (Δα, Δδ) is required. The average changes to the astrometric solution are of order 1–2 arcsec, comparable to the pixel scale and similar to the source positional uncertainty (see Section 2.5.2). Table 1 lists the reference catalogues used for each field.

2.4 Statistics

2.4.1 Area coverage

The pong scanning strategy results in maps that are uniformly deep over the nominal scanning area; however, the usable area in each map is larger than this because of overscan, with radially increasing noise due to the lower effective exposure time in these regions. Although shallower than the map centres, these annular regions around the perimeters of the fields are deep enough to detect sources. Fig. 5 shows the cumulative area of the survey as a function of (instrumental) noise. The total survey area is approximately 5 deg², with >90 per cent of the survey area reaching a sensitivity of under 2 mJy beam⁻¹.

Figure 5.

Cumulative area of the SCUBA-2 Cosmology Legacy Survey as a function of sensitivity, compared to the largest previous 850 μm surveys SHADES (Coppin et al. 2006) and (at 870 μm) LESS (Weiß et al. 2009). The majority of S2CLS reaches a sensitivity of below 2 mJy beam⁻¹, a dramatic step forward compared to previous surveys in the same waveband.

2.4.2 Modelling the PSF

The matched-filtering step described in Section 2.2 modifies the shape of the instrumental PSF, effectively slightly broadening it and increasing the depth of bowling. We derive an empirical PSF by stacking 322 >5σ significance point sources in the UKIDSS-UDS map and fit an analytic surface function to the average profile. The profile is shown in Fig. 6 in comparison to the instrumental PSF, and has an fwhm of 14.8 arcsec. Two-dimensional fitting of the stack reveals that the beam profile P(θ) is circular to within 1 per cent and can be fit with the superposition of two Gaussian functions:

\begin{equation} P(\theta ) = A\exp \left(\frac{\theta ^2}{2\sigma ^2}\right) - 0.98A\exp \left(\frac{\theta ^2}{2.04\sigma ^2}\right) \end{equation}

(1)

Figure 6.

Model of the SCUBA-2 PSF. The dashed line shows the instrumental PSF (Dempsey et al. 2013), and the points show the shape of the average point source in the UKIDSS-UDS field, derived by stacking all sources detected at 5σ significance or greater. The maps are match-filtered, which includes a smoothing step that slightly broadens the instrumental PSF and deepens ‘ringing’. The empirical PSF is well modelled with the superposition of two Gaussians (Section 2.3.2), is circular, and has an fwhm of 14.8 arcsec.

with A = 41.4 and σ = 9.6 arcsec.

2.4.3 The confusion limit

The confusion limit (Jauncey 1968) σ_c is the flux level at which the pixel-to-pixel variance σ² no longer reduces with exposure time due to crowding of the beam by faint sources. The total variance is a combination of the instrumental noise σ_i (in units of mJy beam⁻¹√s) and the confusion noise (in units of mJy beam⁻¹):

\begin{equation} \sigma ^2 = \sigma ^2_it^{-1} + \sigma ^2_c. \end{equation}

(2)

We can evaluate the confusion limit by measuring σ² directly from the pixel data in a progression of maps as we sequentially co-add new scans. Fig. 7 shows how the variance evolves as a function of inverse pixel integration time for the central 15 arcmin of the UKIDSS-UDS, which reaches an instrumental noise of 0.8 mJy beam⁻¹. The best-fitting σ_c is 0.8 mJy beam⁻¹; this confusion noise should be added in quadrature to instrumental and deboosting (Section 2.5.1) uncertainties when considering the flux density of sources. In Section 3, we revisit the estimate of the confusion limit with knowledge of the source counts which allows us to analytically assess the contribution to the noise from rms fluctuations in the flux density due to faint sources below a given limit.

Figure 7.

Measurement of the 850 μm confusion limit for SCUBA-2: we progressively co-add single exposures of the UKIDSS-UDS field, measuring the pixel-to-pixel root mean square value in the uniform central 15 arcmin of the beam-convolved flux map, whilst also tracking the fall off in the pure instrumental noise estimate. At infinite exposure, the instrumental noise is projected to reach zero, whereas the non-zero intercept of the observed flux rms is the confusion limit (equation 2). We measure this to be σ_c ≈ 0.8 mJy beam⁻¹ averaged over the field. Note that the exposure time is the average per 2 arcsec pixel.

2.5 Source extraction

The matched-filtering step optimizes the maps for the detection of point sources – i.e. emission features identical to the PSF. To extract and catalogue sources, we employ a simple top–down peak-finding algorithm: starting from the most significant peak in the signal-to-noise ratio map, the peak flux, noise and position of a source is catalogued before the source is removed from the flux (and signal-to-noise) map by subtracting a scaled version of the model PSF. The highest peak in the source-subtracted map is then catalogued and subtracted and so-on until a floor threshold significance is reached, below which ‘detections’ are no longer trusted. Note that this procedure can potentially deblend sources with markedly different fluxes. The floor detection limit is set to 3σ which allows us to explore the properties of the lowest-significance detections, noting that further cutting can be performed directly on the catalogue. In the following, we assume a cut of 3.5σ as the formal detection limit of S2CLS, where we estimate that the false detection rate is approximately 20 per cent (see Section 2.5.3).

2.5.1 Completeness and flux boosting

To evaluate source detection completeness, we insert fake sources matching a realistic number count distribution into the jackknife noise maps of each field and then try to recover them using the source detection algorithm described above. We adopt the differential number counts fit of Casey et al. (2013) as a fiducial model, which has the Schechter form:

\begin{equation} \frac{\text{d}N}{\text{d}S} = \left(\frac{N_0}{S_0}\right)\left(\frac{S}{S_0}\right)^{-\gamma }\exp \left(-\frac{S}{S_0}\right) \end{equation}

(3)

with N₀ = 3300 deg⁻², S₀ = 3.7 mJy and γ = 1.4. We insert sources down to a flux density limit of 1 mJy and each source is placed at a random position into each map (we do not encode any clustering of the injected sources). An injected source is recovered if a point source is found above the detection threshold within 1.5 × FWHM of the input position. This is a somewhat arbitrary, but generous, threshold, and if there are multiple injected sources within this radius, then we take the closest match. Note that this is a blind approach – no prior is given for the estimated position of injected sources. This procedure is repeated 5000 times for each map, generating a set of mock catalogues containing millions of sources with a realistic flux distribution, allowing us to assess the completeness and flux-boosting statistics.

The ratio of recovered sources to total number of input sources is evaluated in bins of input flux density and local (instrumental) noise. When applying completeness corrections, we use the binned values as a look-up table, using two-dimensional spline interpolation to estimate the completeness rate for a given source. Fig. 8 compares the average completeness of each field (i.e. at the average depth of each map) as a function of intrinsic flux density. Table 2 lists the average 50 per cent and 80 per cent completeness limits for each field and the number of sources above each limit.

Figure 8.

Completeness of the different S2CLS fields, derived from the recovery rate of fake sources injected into jackknife maps as a function of input flux, where a successful recovery at a detection significance of 3.5σ. Note that the completeness falls to zero at 1 mJy as this corresponds to the limit of the injected source model; in practice, it is possible that sub-mJy sources could be boosted above the detection limit. The 50 per cent and 80 per cent limits of each field are listed in Table 2.

Table 2.

50 per cent and 80 per cent completeness limits for the S2CLS fields, quoted at the median map depth (Table 1). We also present the number of sources brighter than the 50 per cent and 80 per cent limits in each field (N_{50, 80}). Note that these flux densities refer to the deboosted – i.e. intrinsic – flux densities. At the 5σ level, observed flux densities are typically overestimated by 20 per cent (Section 2.5.1).

Field	50 per cent	80 per cent	N50	N80
	(mJy)	(mJy)
Akari-NEP	4.1	5.2	132	59
UKIDSS-UDS	3.0	3.8	543	302
COSMOS	4.9	6.2	302	181
Lockman Hole North	3.6	4.6	96	49
GOODS-N	3.9	4.7	32	21
Extended Groth Strip	3.9	5.0	99	51
SSA22	3.9	4.9	78	38

We can simultaneously evaluate flux boosting as a function of local noise and observed flux density simply by comparing the recovered flux to the input flux density of each source. Flux boosting is the overestimation of source flux when measurements are made in the presence of noise and is related to both Eddington and Malmquist bias. Due to the statistical nature of boosting, a source with some observed flux density S_obs is actually drawn from a distribution of true flux density, p(S_true). Our recovery procedure allows us to estimate p(S_true), since we can simply measure the histogram of the injected flux density of sources in bins of (S_obs, σ). This method can be compared to the traditional Bayesian technique to estimate boosting (e.g. Jauncey 1968; Coppin et al. 2005), such that the posterior probability distribution for an observed flux density can be expressed:

\begin{equation} p(S_{\rm true}|S_{\rm obs},\sigma )=\frac{p(S_{\rm true})p(S_{\rm obs},\sigma |S_{\rm true})}{p(S_{\rm obs},\sigma )}. \end{equation}

(4)

The likelihood of the data is given by assuming a Gaussian photometric error on the observed flux density, and the prior is simply the same assumed number counts model used in the simulations described above. Fig. 9 compares the empirically estimated p(S_true) and the posterior probability distribution for S_true from equation (4). The empirical distributions are truncated at 1 mJy because this is the faint limit of the injected source distribution; clearly, we cannot track individual sources fainter than this. An identical counts model is used as a prior in the Bayesian approach, but note that the posterior flux density distribution does extend below 1 mJy; this is because it is effectively the product of a Gaussian (the observed flux density and instrumental uncertainty) and the histogram of pixel values in a map of sources drawn from the model number counts, convolved with the beam. The two methods return similar results, although the empirical method systematically predicts a slightly smaller boosting factor |$\mathcal {B} = S_{\rm obs}/S_{\rm true}$| than the Bayesian approach, with the two methods converging as S_obs increases. Note that neither method assumes any clustering of sources, which could well be important (Hodge et al. 2013; Simpson et al. 2015a).

Figure 9.

Comparison of deboosted flux density distributions for a Bayesian and empirical ‘recovery’ method (Section 2.5.1), using the UKIDSS-UDS field as an example. Both deboosting methods involve considering a model source distribution (down to a flux density of 1 mJy in this case). Each panel shows an observed flux probability distribution, assuming Gaussian uncertainties, for increasing observed flux. The solid and hatched distributions show the predicted intrinsic flux distribution for the Bayesian and direct methods, respectively. In general, the average boosting measured by the two methods agree well, converging as observed flux density increases; however, the ‘direct’ method systematically predicts less boosting compared to the Bayesian approach; we discuss this in the main text.

There are two important differences in the deboosting methods that may explain this: (i) the Bayesian approach does not consider noise (aside from the confusion noise arising from convolving the fake map with the beam), and, related, (ii) the posterior flux distribution derived in equation (4) is not necessarily measured ‘at peak’, i.e. does not consider that the recovered position of a source can shift due to the presence of noise; in the empirical method, we account for such shifts. This relates to the ‘bias-to-peak’ discussed by Austermann et al. (2010). We adopt the ‘empirical’ approach in this work to deboost observed fluxes: we draw samples from the distribution of S_true for a given (S_obs, σ) and calculate the mean and variance of these true fluxes, with the latter providing the uncertainty on the deboosted flux density (provided in the source catalogue). We summarize the empirically derived completeness and boosting for each field, visualized in the plane of flux density and local instrumental noise, in Fig. 10. In Fig. 11, we show the average flux boosting as a function of signal-to-noise ratio in each field, indicating that at fixed detection significance, the level of flux boosting is consistent across the survey, with observed flux densities approximately 20 per cent higher on average than the intrinsic flux density at the 5σ level. The average boosting is well described by a power law:

\begin{equation} \mathcal {B} = 1 + 0.2\left( \frac{{\rm SNR}}{5}\right)^{-2.3}. \end{equation}

(5)

Figure 10.

Two-dimensional visualizations of the results of the recovery simulation in each field. The first column shows the number of artificial sources injected per bin of input flux density and local instrumental noise (labels are log₁₀(N)). The prominent horizontal ridges clearly show the typical depth of the map. The middle column shows the completeness as a function of true flux density and local instrumental noise and the last column shows the average flux boosting as a function of observed flux density and local instrumental noise. The dashed line shows the 3.5σ detection limit.

Figure 10

– continued

Figure 11.

Average flux boosting as a function of signal-to-noise ratio, showing consistency at a fixed signal-to-noise level across different fields. The boosting can be described by a power law; however, in practice, we deboost sources individually based on their observed flux density and local instrumental noise, and drawing on the full probability distribution of true flux densities derived from our recovery simulation.

2.5.2 Positional uncertainty

The simulations described above allow us to investigate the scatter in the difference between input position and recovered position. Like the completeness and boosting, we evaluate the average δθ between input and recovered position in bins of input flux density and local instrumental noise. Following Condon (1997) and Ivison et al. (2007), for a given (Gaussian-like) beam, the positional accuracy is expected to scale with signal to noise. Fig. 12 shows the mean difference between input and recovered source position as a function of signal-to-noise ratio for each field. We find that the positional uncertainty of S2CLS sources is well described by a simple power law, reminiscent of equation B22 of Ivison et al. (2007):

\begin{equation} \delta \theta = 1.2{arcsec}\times \left(\frac{{\rm SNR}}{5}\right)^{-1.6}. \end{equation}

(6)

Figure 12.

Average positional error based on the difference between input and recovered (peak) position from our recovery simulation. All fields follow a similar trend, with the positional uncertainty decreasing with increasing source significance. We fit the uncertainties with a simple power law to estimate the 1σ positional uncertainty as a function of observed signal-to-noise ratio (Section 2.5.2).

2.5.3 False detection rate

To measure the false detection rate, we compare the number of ‘detections’ in the jackknife maps to those in the real maps as a function of signal-to-noise ratio. By construction, the jackknife maps contain no astronomical signal and have Gaussian noise properties (Fig. 4); therefore, any detections are due to statistical fluctuations expected from Gaussian noise at the ≥3.5σ level. Fig. 13 shows the false detection rate as a function signal-to-noise ratio; at our 3.5σ limit the false detection (or contamination) rate is 20 per cent, falling to 6 per cent at 4σ and falls below 1 per cent for a ≥5σ cut. The false detection rate is as follows:

\begin{equation} \log _{10}(\mathcal {F}) = 2.67 - 0.97\times {\rm SNR}. \end{equation}

(7)

An alternative approach to estimating the false detection rate that takes into account the presence of real sources in the map uses the Bayesian estimate of the posterior probability distribution of the flux density per source; the integral of equation (4) at S ≤ 0 mJy can be taken as the probability that a source is a false detection (e.g. Coppin et al. 2006). We confirm that estimating the false detection rate in this manner gives results consistent with the ‘pure noise’ estimate captured by equation (7), indicating that false positives are dominated by Gaussian statistics.

Figure 13.

False detection rate averaged over the survey defined as the ratio of ‘detections’ in jackknife maps to real detections for sources at a fixed signal-to-noise limit. At our 3.5σ limit, the false detection (or contamination) rate is 20 per cent, falling to 6 per cent at 4σ and is negligible for a ≥5σ cut. The implication is that, although the final S2CLS catalogue is cut at 3.5σ, caution should be taken in the consideration of individual sources below a significance of 5σ.

Equation (7) implies that caution should be taken when considering individual sources in the S2CLS catalogue at detection significance of less than 5σ; follow-up confirmation and/or robust counterpart identification will be important for assessing the reality of sources detected close to the survey limit, and this work has already begun (e.g. Chen et al. 2016).

3 NUMBER COUNTS OF THE 850 μM POPULATION

In Table 3, we present a sample of the S2CLS catalogue. The full catalogue contains 2851 sources at a detection significance of ≥3.5σ. The catalogue contains observed and deboosted flux densities, instrumental and deboosted flux density uncertainties, and individual completeness and false detection rates. The full catalogue and maps (match-filtered and non-match-filtered) are available at the DOI: http://dx.doi.org/10.5281/zenodo.57792. Appendix 1 gives a complete description of the catalogue columns.

Table 3.

Sample of the full S2CLS catalogue, listing the highest and the lowest significance detections in each field. Coordinates are J2000, with the individual map astrometric solutions tied to the reference catalogue listed in Table 1. The |$S^{\rm obs}_{850}\pm \sigma _{\rm inst}$| column gives the observed flux density and instrumental noise, S/N gives the detection signal-to-noise ratio, and S₈₅₀ ± σ_tot gives the estimated true flux density and combined total (instrumental, deboosting, confusion) uncertainty. The final column |$\log _{10}(\mathcal {F})$| is the logarithm of the false detection rate for the detection signal-to-noise ratio (equa-tion 7), negligible for bright sources, but important to consider for sources at the detection limit.

S2CLS ID	Short name	R.A.	Dec.	|$S^{\rm obs}_{\rm 850}\pm \sigma _{\rm inst}$|	S/N	S850 ± σtot	|$\langle \mathcal {C} \rangle$|	|$\log _{10}(\mathcal {F})$|
S2CLSJ175833+663757	NEP.0001	17 58 33.60	+66 37 57.7	195.4 ± 1.2	158.4	195.4 ± 1.5	1.00	−147.78
S2CLSJ175416+665117	NEP.0329	17 54 16.57	+66 51 17.0	4.3 ± 1.2	3.5	2.9 ± 1.9	0.25	−0.66
S2CLSJ100015+021548	COS.0001	10 00 15.72	+02 15 48.6	12.9 ± 0.8	15.2	12.9 ± 1.2	1.00	−11.78
S2CLSJ095936+022506	COS.0733	09 59 36.09	+02 25 06.5	5.4 ± 1.5	3.5	3.6 ± 2.4	0.25	−0.65
S2CLSJ141951+530044	EGS.0001	14 19 51.56	+53 00 44.8	16.3 ± 1.2	14.1	16.3 ± 1.4	1.00	−10.69
S2CLSJ141612+521316	EGS.0227	14 16 12.05	+52 13 16.8	3.7 ± 1.0	3.5	2.6 ± 1.7	0.28	−0.66
S2CLSJ123730+621258	GDN.0001	12 37 30.73	+62 12 58.5	12.8 ± 1.0	13.2	11.9 ± 1.6	1.00	−9.82
S2CLSJ123734+620736	GDN.0068	12 37 34.33	+62 07 36.5	5.0 ± 1.4	3.5	3.3 ± 2.2	0.23	−0.65
S2CLSJ104635+590748	LHO.0001	10 46 35.78	+59 07 48.0	12.0 ± 1.0	11.6	11.5 ± 1.8	0.99	−8.31
S2CLSJ104541+584640	LHO.0219	10 45 41.47	+58 46 40.0	4.1 ± 1.2	3.5	3.0 ± 1.8	0.29	−0.67
S2CLSJ221732+001740	SSA.0001	22 17 32.50	+00 17 40.4	14.5 ± 1.1	13.0	14.5 ± 1.4	0.96	−9.72
S2CLSJ221720+002024	SSA.0198	22 17 20.23	+00 20 24.4	3.9 ± 1.1	3.5	2.8 ± 1.8	0.28	−0.65
S2CLSJ021830-053130	UDS.0001	02 18 30.77	−05 31 30.8	52.7 ± 0.9	56.7	52.7 ± 1.2	1.00	−51.18
S2CLSJ021823-051508	UDS.1080	02 18 23.12	−05 15 08.9	3.1 ± 0.9	3.5	2.4 ± 1.5	0.30	−0.65

The surface density of sources per observed flux density interval dN/dS – of a cosmological population – is a simple measure of source abundance and a powerful tool for model comparisons. To measure the counts, for each catalogued source we first deboost the observed flux density using the empirical approach described in Section 2.5.1, and then apply the corresponding completeness correction for the deboosted (i.e. ‘true’) flux density. When deboosting, we consider the full intrinsic flux distribution as estimated by our simulation, accounting for the fact that a range of intrinsic flux densities can map on to an observed flux density. Therefore, we evaluate dN/dS 1000 times; in each calculation, every source is deboosted by randomly sampling the intrinsic flux distribution and completeness correcting each deboosted source accordingly. We take the mean of these 1000 realizations as the final number counts, with the standard deviation of dN/dS in each bin as an additional uncertainty (to the Poisson error). We make a correction for each source based on the probability it is a false positive, using the empirical determination described in Section 2.5.3.

While the various corrections are intended to recover the ‘true’ underlying source distribution, it is important to confirm if any systematic biases remain, since the procedure for actually identifying sources is imperfect, as is the ‘recovery’ of injected model sources used to estimate flux boosting and completeness. To examine this, we inject three different source count models into a jackknife noise map (of the UKIDSS-UDS field). One model is identical to the Schechter form used in Section 2.6.1 (equation 3); in the other two models, we simply adjust the faint-end slope to γ = 0.4 and γ = 2.4, keeping the other parameters fixed. With knowledge of the exact model counts injected into the map, we can compare to the recovered counts before and after corrections have been applied. Fig. 14 shows ([dN/dS]_rec − [dN/dS]_true)/[dN/dS]_true for the three models before and after corrections. In the absence of correction, flux boosting tends to result in the systematic overestimation of the number counts in all but the faintest flux bin, where incompleteness dominates, and the overestimation increases with increasing γ, as expected. After the corrections have been applied, there remains a slight underestimation in the counts in the faintest bin (3–4 mJy) at the 10 per cent level, but in general the corrected ‘observed’ counts are in excellent agreement with the input model. The origin for the slight discrepancy is not clear, but it is likely that it simply stems from subtle effects not modelled well by our recovery simulation, and in particular what constitutes a ‘recovered’ source. One can observe a systematic effect that the γ = 2.4 and γ = 0.4 models are over- and under-estimated (respectively) at approximately the 10 per cent level for the full observed flux range, but this is not a significant systematic uncertainty compared to shot noise expected from Poisson statistics. Given that the fiducial model we use in the completeness simulation is based on observed 850 μm number counts, and the γ = 2.4 and γ = 0.4 models are rather extreme compared to empirical constraints, we consider this test as an adequate demonstration that our measured number counts are robust. Nevertheless, we apply a simple correction to the observed corrected counts by fitting a spline to the residual model counts in Fig. 14 and apply this as a ‘tweak’ factor to the number counts on a bin-by-bin basis.

Figure 14.

A comparison of recovered number counts to an ideal input model (equation 3). Three input models are considered, differing only by the faint-end slope γ: a series of fake catalogues are generated for each model by injecting sources into a jackknife map and then recovering them in a manner identical to the real data. In the left-hand panel, no deboosting, completeness or false positive correction has been applied to the recovered counts, showing the trend that steeper number counts are generally overpredicted (due to flux boosting) in all but the faint bin where incompleteness dominates. In the right-hand panel, the various corrections have been applied, illustrating that we can robustly recover the ‘true’ number counts, although there is still a slight (10 per cent) underestimation of the counts in the faintest bin. The error bars in both panels reflect the Poisson uncertainties expected in a single field. The dashed line is a spline interpolation of the mean of the three models which we use as an additional tweak factor in measuring the number counts of the population. Interestingly, the two extreme count models we consider are, in general, systematically over and under predicted for the steeper and shallower faint-end slope, respectively; we discuss this in the main text.

The S2CLS differential (and cumulative) number counts are presented in Table 4 and Fig. 15. Tables of the number counts of individual fields are available in the electronic version of the paper. S2CLS covers a solid angle large enough to detect reasonable numbers of the rarer, bright sources at S₈₅₀ > 10 mJy, allowing us to robustly measure the bright end of the observed 850 μm number counts. As a guide, there are about 10 sources with flux densities greater than 10 mJy deg⁻². The 850 μm source counts above 10 mJy clearly show an upturn in source density that is due to a mixture of local emitters and gravitationally lensed sources.³ The wide-area counts of Herschel demonstrated the same (predicted) phenomenon in the SPIRE bands (Negrello et al. 2010), and it has since been demonstrated that a simple bright submillimetre flux cut is highly effective at identifying strongly lensed sources once local galaxies have been rejected. The effect has already been observed in the millimetre regime: Vieiran et al. (2010) detect the upturn in the 1.4 mm counts at S_1.4 mm > 10 mJy from SPT over a 87 deg² survey, and Scott et al. (2012) have reported tentative evidence of an upturn in the counts at 1.1 mm at S_1.1 mm > 13 mJy with AzTEC over 1.6 deg². Much larger 850 μm surveys (exceeding 10 deg²) could utilize a similar selection to cleanly identify lensed 850 μm selected high-redshift galaxies.

Figure 15.

Number counts of 850 μm sources. The left-hand panel shows the differential number counts for individual fields and the combined survey, along with a selection of data from the literature. Two model curves show the parametric evolving luminosity function model of Béthermin et al. (2012) and the semi-analytic (galform) model of Lacey et al. (2016). The Cowley et al. (2015) line shows the same galform model but taking into account source blending due to the 15 arcsec JCMT beam. The presence of foreground sources and the effect of gravitational lensing causes an upturn in the counts at bright flux densities at a level in reasonable agreement with the models (note that the galform model does not include lensing) given the low number statistics at these bright flux densities. For clarity, we only show error bars for the S2CLS data, which are evaluated from Poisson statistics (Gehrels 1986). The right-hand panel shows the cumulative counts where, for clarity, we only plot the S2CLS data, fit, and models.

Table 4.

Number counts measured in the full S2CLS. Flux density bins ΔS are 1 mJy wide. The flux density S is the bin central and S΄ = S − 0.5ΔS. Uncertainties on the counts are written such that the first set of errors are Poissonian and the second reflect the standard deviation of each bin of dN/dS after 1000 realizations of the counts, where each source is deboosted (and completeness corrected) by sampling the deboosting probability distribution corresponding to the observed flux density and local noise. These uncertainties are of comparable magnitude to the Poisson errors.

S	dN/dS	N( > S΄)
(mJy)	(deg−2 mJy−1)	(deg−2)
3.5	|$451.0^{+17.1}_{-16.4}\pm 20.3$|	|$1012.3^{+19.6}_{-19.2}\pm 19.6$|
4.5	|$204.4^{+9.3}_{-8.9}\pm 8.8$|	|$508.0^{+12.3}_{-12.0}\pm 9.7$|
5.5	|$102.6^{+6.0}_{-5.7}\pm 5.1$|	|$271.9^{+8.5}_{-8.2}\pm 6.5$|
6.5	|$56.1^{+4.3}_{-4.0}\pm 3.8$|	|$151.8^{+6.2}_{-6.0}\pm 4.3$|
7.5	|$32.5^{+3.2}_{-2.9}\pm 2.5$|	|$85.3^{+4.7}_{-4.4}\pm 3.1$|
8.5	|$18.0^{+2.5}_{-2.2}\pm 2.0$|	|$47.1^{+3.6}_{-3.3}\pm 2.3$|
9.5	|$9.8^{+1.9}_{-1.6}\pm 1.4$|	|$26.4^{+2.8}_{-2.5}\pm 1.6$|
10.5	|$5.8^{+1.5}_{-1.2}\pm 1.0$|	|$14.5^{+2.2}_{-1.9}\pm 1.2$|
11.5	|$3.4^{+1.2}_{-0.9}\pm 0.8$|	|$8.7^{+1.8}_{-1.5}\pm 0.8$|
12.5	|$2.1^{+1.1}_{-0.7}\pm 0.6$|	|$5.5^{+1.5}_{-1.2}\pm 0.6$|
13.5	|$0.8^{+0.8}_{-0.4}\pm 0.4$|	|$3.2^{+1.2}_{-0.9}\pm 0.5$|
14.5	|$0.5^{+0.7}_{-0.3}\pm 0.3$|	|$2.4^{+1.1}_{-0.8}\pm 0.3$|
15.5	|$0.3^{+0.6}_{-0.2}\pm 0.1$|	|$1.8^{+1.0}_{-0.7}\pm 0.2$|

The S2CLS number counts are in reasonable agreement with previous surveys for the flux range probed (for clarity, a non-exhaustive list of previous surveys, including recent SCUBA-2 results, are shown in Fig. 15: Coppin et al 2006; Weiß et al. 2009; Casey et al. 2013; Chen et al. 2013), but with the large number of sources in S2CLS, we can dramatically reduce the Poisson errors: in the faintest bin, the Poisson uncertainty on the differential counts over the whole survey is just ∼4 per cent. We fit the combined differential counts (up to 20 mJy after which the local/lensing upturn starts to contribute significantly) with the Schechter functional form given in equation (3). We find the best-fitting parameters N₀ = 7180 ± 1220 deg⁻², S₀ = 2.5 ± 0.4 mJy beam⁻¹ and γ = 1.5 ± 0.4.

With the counts well measured, we can revisit the estimate of the confusion limit presented in Section 2.4.3, since this is driven by the rms fluctuations in the integrated flux due to faint sources below some flux limit S_c. Following Helou & Beichman (1990), we can express the confusion noise as

\begin{equation} \sigma _{\rm c}^2(S_{\rm c}) = \Omega _{\rm b}\int _0^{S_{\rm c}} S^2 \frac{\text{d}N(S)}{\text{d}S} \text{d}S, \end{equation}

(8)

where Ω_b is the SCUBA-2 beam area (242 arcsec² for the 850 μm beam). If we define the confusion limit as |$\sigma _{\rm c}^2(S_{\rm c})$| with S_c → ∞ we find σ_c = 0.86 mJy beam⁻¹, in good agreement with the value measured in Section 2.4.3. Note that this can be considered an upper limit to the confusion noise contribution since it integrates over the full population.

3.1 Field-to-field variance

Taking the full survey counts as an average measure of the abundance of submillimetre sources, with S2CLS we can now investigate field-to-field variance in the number counts in a consistent manner; this is important given that SMGs are thought to be a highly biased tracer of the matter field (Hickox et al. 2012; Chen et al. 2016). Letting ρ(S) = N( > S), for each field, we can consider the deviation of the counts compared to the mean density per flux bin: δ(S) = (ρ(S) − 〈ρ(S)〉)/〈ρ(S)〉. In Fig. 16, we show the δ(S) measured for each field as a function of flux density, where uncertainties are the combined Poisson errors (obviously dominated by the single-field counts). The field-to-field scatter on ∼0.5°–1° scales is generally within 50 per cent of the survey-averaged density and reasonably consistent with the Poisson errors. There are some hints that the GOODS-N field has a slightly elevated density compared to the mean (hints that were already apparent in the original SCUBA maps of this field, see Pope et al. 2005; Coppin et al. 2006; Walter et al. 2012), but this is marginal given the Poisson errors. However, to explore this further, and to quantify the significance of any overdensity, we can evaluate the field-to-field fluctuations on scales equivalent to the GOODS-N field taking into account cosmic variance.

Figure 16.

Field-to-field scatter in the integral number counts, relative to the mean density. The field-to-field scatter (on scales of 0.5°–1°) across S2CLS is generally within 50 per cent of the mean density, with the exception of GOODS-N, which has hints of an elevated density of SMGs compared to the mean, although this is marginal with the Poisson uncertainties. We discuss this in Section 3.1.

Field-to-field variance in the observed number counts is caused by both shot noise and cosmic variance, with the latter defined as the excess variance in addition to Poisson noise (e.g. Somerville et al. 2004). We split the S2CLS survey (barring GOODS-N) into 16 independent fields of identical size to the GOODS-N field and count the number of sources in each field with deboosted flux densities greater than the 50 per cent completeness limit in GOODS-N (S₈₅₀ ≈ 4 mJy). The mean number of sources is 18, with a standard deviation over the 16 fields of 5 sources, roughly consistent with the shot noise expected from Poisson statistics. The number of sources at the same limit in GOODS-N is 32 ± 6. It is clear from this simple analysis that the error budget on the counts is dominated by Poisson noise, but we can estimate what the expected contribution from cosmic variance is. Following the method of Trenti & Stiavelli (2008) which estimates the relative excess uncertainty in number counts due to cosmic variance in a flux limit survey, we find a contribution of 15–20 per cent to the observed counts on scales of the GOODS-N field (note, we assume a Press–Schechter approach for the halo statistics). This assumes a mean redshift of 〈z〉 = 2.2 and Δz = 1 and a wide range of halo filling factors, f = 0.1–1, corresponding to a mean bias of b = 2.7–4.3 for the SMG population. We can therefore quantify the significance of the tentative overdensity in GOODS-N as the difference in the number of sources in this field to the average over a region 16 times larger in an independent field (i.e. the rest of the S2CLS). We find ΔS(>4mJy) = 14 ± 7 taking into account Poisson noise and cosmic variance. Thus, the overdensity is significant at only the 2σ level. GOODS-N is one of the most exhaustively studied extragalactic fields, and it is worth noting that overdensities of SMGs and star-forming galaxies have previously been reported here. For example, Daddi et al. (2009) report an overdensity of star-forming galaxies at z ≈ 4, including SMGs, and Walter et al. (2012) report a z ≈ 5 structure around the source HDF850.1, which happens to be one of the first SMGs to be identified (Hughes et al. 1998).

3.2 Comparison to models

At first, it proved difficult for semi-analytic models of ΛCDM galaxy formation to reproduce the 850 μm number counts (Granato et al. 2000). The model of Baugh et al. (2005) provided a much better match to observed 850 μm (and Lyman-break Galaxy) number counts than previously achieved, but required a modification to the initial mass function (IMF) such that bursts of star formation have a more top heavy IMF than ‘quiescent’ star formation. While the motivation for this can be linked to astrophysical differences in the conditions of star formation in dense gas-rich starbursts (e.g. Padoan & Nordlund 2002), deviation from a universal IMF remains controversial. An additional problem was that the Baugh et al. (2005) model failed to predict the evolution of the K-band luminosity function. Recently, Lacey et al. (2016) presented an update to the galform model that adopts the best-fitting ΛCDM cosmological parameters available from recent experiments, implementing more sophisticated treatments for star formation in discs, distinguishing molecular and atomic hydrogen (Lagos et al. 2011, 2012); dynamical friction time-scales for mergers (Jiang et al. 2008) and stellar population synthesis models.

The Lacey et al. model counts are shown in Fig. 15 , and are in reasonable agreement with the data. This model still includes a mildly top-heavy IMF (slope x = 1) for starbursts, without which it cannot reproduce the redshift distribution of 850 μm selected sources (Chapman et al. 2005; Simpson et al. 2014, see also Hayward et al. 2013). The model predicts a slightly elevated abundance of galaxies below the survey limit; however, an extrapolation of the Schechter fit to the S2CLS counts is in good agreement with the deeper observations of Chen et al. (2013). Nevertheless, the shape of the faint-end slope is still to be properly determined empirically, which will most likely be through either a P(D) analysis of confused SCUBA-2 maps (e.g. Condon 1974; Pantanchon et al. 2009; Geach et al., in preparation), with the assistance of gravitational lensing (e.g. Knudsen, van der Werf & Kneib 2008; Chen et al. 2013) or through deep, unconfused ALMA surveys that can probe to the sub-mJy level, albeit over relatively small areas (e.g. Karim et al. 2013; Ono et al. 2014; Carniani et al. 2015; Dunlop et al. 2016; Hatsukade et al. 2016; Oteo et al. 2016). An important point to consider in comparing number counts to models is the issue of source blending and confusion in low-resolution single-dish surveys. Therefore, we also show the results of Cowley et al. (2015), who take the same Lacey et al. (2016) galform model, but predict the number counts after simulating observations with a single-dish telescope with the same size beam as JCMT at 850 μm. Fig. 15 shows that, over the observed flux density range, the beam-convolved predicted counts are consistent with the ‘raw’ model counts. The issue of ‘multiplicity’ of single-dish SMG detections has already started to be examined with the advent of sensitive interferometers (e.g. Hodge et al. 2013; Simpson et al. 2015a), and it is important to stress that comparisons of source abundances (between both models and data) should adopt a consistent reference resolution.

While the semi-analytic models aim to simultaneously reproduce all the main ‘bulk’ observational tracers of the galaxy population over cosmic time (i.e. the mass function, luminosity functions, number counts, clustering, etc.) in a single framework, an alternative approach to predicting the submillimetre number counts is through phenomenological modelling. Bethermin et al. (2012) present a model that considers the evolution of the space density of so-called main sequence (i.e. normal) star-forming galaxies and luminous starbursts, fitting parametric models (with assumptions about the underlying galaxy SEDs) to observed number counts across the infrared, submillimetre and radio bands. We show the Bethermin et al. model (including the strong lensing contribution) for the SCUBA-2 850 μm band in Fig. 15. Again, this is in reasonable agreement with the observations over the flux range probed by the observations. The new 850 μm number counts presented here could be used to provide improved fits to phenomenological models such as this.

4 SUMMARY

We have presented the 850 μm maps and catalogues of the James Clerk Maxwell Telescope SCUBA-2 Cosmology Legacy Survey, the largest of the JCMT Legacy Surveys, completed in early 2015. With hundreds of hours of integration time in reasonable submillimetre observing conditions (zenith opacity τ_225 GHz = 0.05–0.1), S2CLS has mapped seven well-known extragalactic survey fields: UKIDSS-UDS, Akari-NEP, COSMOS, GOODS-N, Extended Groth Strip, Lockman Hole North and SSA22. The total scientifically useful survey area is approximately 5 deg² at a sensitivity of under 2 mJy beam⁻¹, with a median depth per field of approximately 1.2 mJy beam⁻¹, approaching the confusion limit (which we have determined is approximately σ_c ≈ 0.8 mJy beam⁻¹). This is by far the largest and deepest survey of submillimetre galaxies yet undertaken in this waveband and provides a rich legacy data source. We have detected nearly 3000 submillimetre sources at the ≥3.5σ level, an order of magnitude increase in the number of catalogued 850 μm selected sources to date.

In this work, we have used the S2CLS catalogue to accurately measure the number counts of submillimetre sources, dramatically reducing Poisson errors and allowing us to investigate field-to-field variance. The wide nature of the survey makes it possible to detect large numbers of bright (>10 mJy), but rare (∼10 deg⁻² ), submillimetre sources, and we observe the distinctive upturn in the number counts caused by strong gravitational lensing of high-redshift galaxies and a contribution from local sources of submillimetre emission. The S2CLS catalogue and maps offer a route to a tremendous range of follow-up work, both in pin-pointed multiwavelength identification and follow-up of the catalogued sources (e.g. Simpson et al. 2015a,b; Chen et al. 2016) and in statistical analyses of the catalogues and pixel data. Cross-correlation of the submillimetre maps and galaxy catalogues is already proving a treasure trove of discovery, linking UV/optical/near-infrared-selected samples to submillimetre emission (Banerji et al. 2015; Coppin et al. 2015; Smith et al., in preparation; Bourne et al., in preparation). The S2CLS survey subtends an area equivalent to over 10⁵ times the ALMA primary beam at 850 μm, and the synergy between large-area single-dish surveys such as S2CLS, and the detailed interferometric follow-up possible with ALMA (and other sensitive (sub)mm interferometers) is clear. High-resolution interferometric follow-up in the submillimetre has already proven efficient and fruitful, with ALMA and SMA imaging of the brightest (>9 mJy) sources revealing a complex morphological mix, allowing us to investigate the true nature of the SMGs identified in large-beam single-dish surveys (Simpson et al. 2015a; Chapman et al., in preparation). We release the 3.5σ-cut catalogue of all S2CLS sources as part of this publication, along with the 850 μm maps for exploitation by the community. The data are available at the DOI: http://dx.doi.org/10.5281/zenodo.57792.

Acknowledgments

The authors thank M. Bethermin for supplying the SCUBA-2 850 μm number count predictions. The James Clerk Maxwell Telescope is now operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (grant no.XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The data presented in this paper were taken as part of Program ID MJLSC02. It is a pleasure to thank the entire staff of the JCMT for their superb support throughout the S2CLS campaign. Special thanks is due to Iain Coulson, Jessica Dempsey, Jim Hoge, Harriet Parsons, Callie Matulonis, William Montgomerie and Holly Thomas. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure.

JEG is supported by a Royal Society University Research Fellowship. IRS, AMS, JMS, ALRD, DMA, ACE and CGL acknowledge support from ST/L00075X/1. IRS also acknowledges support from the ERC Advanced Grant DUSTYGAL (321334) and a Royal Society Wolfson Merit Award. JSD acknowledges the support of the European Research Council through the award of an Advanced Grant. RJI acknowledges support from ERC in the form of the Advanced Investigator Programme COSMICISM (321302). AK acknowledges support by the Collaborative Research Council 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG). KK acknowledges supports from the Swedish Research Council.

1

It is worth noting that it is now common to refer to SMGs as cosmological sources selected right across the 250–1000 μm wavelength range. With the high-altitude Balloon-borne Large Aperture Submillimeter Telescope (BLAST; Pascale et al. 2008) and then the launch of the Herschel Space Observatory in 2009 (Griffin et al. 2008), the path has been opened up to large area submillimetre surveys at λ ≤ 650 μm (e.g. Eales et al. 2010), although suffering from high confusion noise due to the limited size of dishes that can be flown in the sky and space.

2

A flux limit imposed by confusion.

3

Note that the Akari-NEP field contains the galactic object NGC 6543 (the Cat's Eye Nebula) – which is a ∼200 mJy 850 μm source.

APPENDIX 1: DESCRIPTION OF CATALOGUE

The following describes the content of the S2CLS catalogue.

S2CLS name	IAU compliant coordinate-based catalogue name
Nickname	Short source name of the format field.XXXX, where XXXX is a catalogue index
RA_pix	Right Ascension of source from the peak SNR pixel in sexigesimal format (J2000)
Dec_pix	Declination of source from the peak SNR pixel in sexigesimal format (J2000)
RA_DEG_pix	Right Ascension of source from the peak SNR pixel in decimal degrees (J2000)
Dec_DEG_pix	Declination of source from the peak SNR pixel in decimal degrees (J2000)
RA_gauss	Right Ascension of source from a 2D Gaussian profile fit to the local pixel data in sexigesimal format (J2000)
Dec_gauss	Declination of source from a 2D Gaussian profile fit to the local pixel data in sexigesimal format (J2000)
RA_DEG_gauss	Right Ascension of source from a 2D Gaussian profile fit to the local pixel data in decimal degrees (J2000)
Dec_DEG_gauss	Declination of source from a 2D Gaussian profile fit to the local pixel data in decimal degrees (J2000)
S_850_observed	Flux density measured at peak SNR pixel in mJy beam−1
delta_S_850_inst	Instrumental noise measured at peak SNR pixel in mJy beam−1
detection_SNR	Signal-to-noise ratio defined as S_850_observed / delta_S_850_inst
S_850_deboost	Deboosted flux density in mJy beam−1
delta_S_850_deboost	Uncertainty on deboosted flux density in mJy beam−1
Completeness	Completeness rate for this source based on recovery simulation
log10_false_detection_rate	logarithmic probability that source is a false positive