Abstract

We present maps, source lists and derived number counts from the largest, unbiased, extragalactic submillimetre (submm) survey so far undertaken with the SCUBA camera on the James Clerk Maxwell Telescope (JCMT). Our maps are located in two regions of sky (ELAIS N2 and Lockman-Hole E) and cover 260 arcmin2, to a typical rms noise level of σ850≃2.5 mJy beam−1. We have reduced the data using both the standard JCMT surf procedures, and our own IDL-based pipeline which produces zero-footprint maps and noise images. The uncorrelated noise maps produced by the latter approach have enabled us to apply a maximum likelihood method to measure the statistical significance of each peak in our maps, leading to properly quantified errors on the flux density of all potential sources.

We detect 19 sources with signal-to-noise ratios (S/N)>4, and 38 with S/N>3.5. To assess both the completeness of this survey and the impact of source confusion as a function of flux density, we have applied our source-extraction algorithm to a series of simulated images. The result is a new estimate of the submm source counts over the flux-density range S850≃5–15 mJy, which we compare with estimates derived by other workers, and with the predictions of a number of models.

Our best estimate of the cumulative source count at S850>8 mJy is graphic per square degree. Assuming that the majority of sources lie at z>1.5, this result implies that the comoving number density of high-redshift galaxies forming stars at a rate in excess of 1000 M yr−1 is ≃10−5 Mpc−3, with only a weak dependence on the precise redshift distribution. This number density corresponds to the number density of massive ellipticals with L>3–4L* in the present-day Universe, and is also the same as the comoving number density of comparably massive, passively evolving objects in the redshift band 1<z<2 inferred from recent surveys of extremely red objects. Thus the bright submm sources uncovered by this survey can plausibly account for the formation of all present-day massive ellipticals. Improved redshift constraints, and ultimately an improved measure of submm source clustering can refine or refute this picture.

Introduction

Over the past three years deep submillimetre surveys, carried out to a variety of depths (Smail, Ivison & Blain 1997; Hughes et al. 1998; Barger, Cowie & Sanders 1999; Blain et al. 1999; Eales et al. 2000; Borys et al. 2002) have highlighted the extreme importance of dust in the determination of the global star formation history of the Universe. Optical/ultraviolet (UV) studies originally suggested a steep rise in the star formation rate (SFR) as a function of redshift between z=0 and z=1 (Lilly et al. 1996), peaking at 1.0>z>1.5, and declining to values comparable to those observed in the present day at z≃4 (Madau et al. 1996). More recent optical/UV analyses have implied a much gentler trend in SFR density with epoch, Cowie, Songaila & Barger (1999) finding a more gradual slope [∼(1+z)1.5 for an Einstein–de Sitter cosmology] at low redshift (z>1) than had previously been determined by Lilly et al. (1996), and at high redshift the Lyman break population (Steidel et al. 1999) suggesting an approximately uniform value for 1>z>4. However, in heavily dust-enshrouded star-forming regions, much of the optical/UV radiation emitted by the young stars is absorbed, leading to the possibility that a significant amount of star formation, particularly at high redshift, may have been missed in these wavebands.

In order to quantify the star formation density contributed by such highly obscured objects, we must directly observe the rest-frame far-infrared (FIR) emission from the heated dust, which at redshifts greater than z∼1 is shifted into the submm waveband. The steep spectral index of the thermal emission longwards of the peak at λ≃100 μm results in a sufficiently large negative K-correction that the dimming effect of increasing cosmological distance is effectively cancelled over a wide range in redshift. Consequently, for a galaxy with fixed intrinsic FIR luminosity, we would expect to observe approximately the same 850-μm flux density for an object at z=8 as at z=1. The strength of this effect is obviously dependent on the assumed cosmology. It is most pronounced in an Einstein–de Sitter universe, in which the predicted flux density of an object of fixed luminosity actually increases slightly beyond z≃1. If instead we adopt ΩM=0.3 and ΩΛ=0.7, a very gentle decline in flux density is predicted over this same redshift range.

Indications that a large fraction of star-forming activity at high redshift is heavily enshrouded in dust have been provided, for example, by a deep submm survey, of the northern Hubble Deep Field (HDF, Hughes et al. 1998). The discovery of five discrete sources with either very faint or presently undetected optical identifications, accounting for approximately 30–50 per cent of the submm background observed by COBE-FIRAS (Puget et al. 1996; Fixsen et al. 1998; Hauser et al. 1998), have suggested a star formation density at least a factor of ≃5 higher at z<2 than that originally deduced from optical/UV data. The results of Barger, Cowie & Richards (2000) and Eales et al. (2000) are also consistent with a high-redshift SFR density up to an order of magnitude higher than that inferred in the optical/UV. These sources are believed to be analogous to the local ultraluminous infrared galaxies (ULIRGs) and undergoing a period of massive star formation at rates of 102–103 M yr−1.

The discovery of such objects at high redshift has important implications for the formation mechanism of massive elliptical galaxies. In current cold dark matter (CDM) based models, elliptical galaxies are created at modest redshift by the merging of intermediate-mass discs (Baugh, Cole & Frenk 1996). Although some recent semi-analytic results appear supportive of the formation of massive ellipticals at modest redshift by means of hierarchical mergers (Kauffman & Charlot 1998), observations continue to imply that a massive coeval starburst of time-scales ≃1 Gyr at high redshift (z<3) is required to explain the properties of at least some ellipticals (Jimenez et al. 1999; Renzini & Cimatti 2000). Additionally, the connection between black hole and spheroid mass (Magorrian et al. 1998; Kormendy & Gebhardt 2001) suggests that galaxies hosting an active nucleus may be more representative of the general population than had previously been thought, and might indicate a link between spheroid formation and the epoch of maximum quasar activity at z=2–3.

The discovery of the Lyman-break population at z=2–4 (Steidel et al. 1999) has prompted the suggestion that the formation of present-day galactic bulges and ellipticals has already been observed at optical wavelengths. However, the most massive ellipticals contain 1012 solar masses of stars, requiring a sustained SFR of ∼1000 M yr−1 for formation over a 1-Gyr period. Even if the most optimistic corrections are made for the effects of dust, the inferred star formation rates in these objects (typically 3–30 M yr−1) fall short by 1–2 orders of magnitude. Moreover, the first direct tests of such corrections suggest they may indeed have been exaggerated. Direct measurements of the submm/millimetre emission of the strongly lensed Lyman-break galaxies MS 1512+36-cB58 (Baker et al. 2001; Sawicki 2000; van der Werf et al. 2001) and MS 1358+62-G1 (van der Werf et al. 2001) imply that the procedure used to predict the submm emission of the Lyman-break population over-predicts the observed 850-μm flux densities by a factor of up to 14. A recent X-ray detection of MS 1512+36-cB58 (Almaini, Pettini & Steidel, in preparation) also supports the view that corrections for dust have been exaggereated, the level of X-ray emission implying a much lower SFR than that derived from the ‘dust-corrected’ UV spectrum, but consistent with the measured 850-μm flux density, for a typical starburst spectral energy distribution (SED).

One possibility is that spheroid formation is too widely distributed in space and/or time to give rise to such spectacular starbursts. However, observationally we may be missing massive starbursts in the optical as a result of dust obscuration or because they are at too great a redshift. Submm detections of the z≃4 radio galaxies 4C41.17 (Dunlop et al. 1994; Hughes, Dunlop & Rawlings 1997) and 8C1435+635 (Ivison 1995; Ivison et al. 1998) confirm that the Lyman-break population does not tell the full story. These extreme radio galaxies display the properties expected of youthful massive ellipticals, and have inferred dust masses <108 M and SFRs <1000 M yr−1.

Here we present the first results from the ‘SCUBA 8 mJy Survey’, which is a wider-area, somewhat shallower survey than its earlier counterparts (HDF-N –Hughes et al. 1998; Clusters Survey –Smail et al. 1997; Blain et al. 1999; Hawaii Fields – Barger et al. 1999; CUDSS –Eales et al. 2000), undertaken with the aim of constraining the brighter end of the 850-μm source counts in the region of 8 mJy. This survey has three key advantages over previous submm surveys. First, the relatively bright flux-density limit facilitates follow-up with existing instrumentation which should eventually yield an accurate astrometric position for every source via either the VLA, or in the case of the very high-redshift sources (z<3–4) IRAM PdB. Secondly, the steep slope of the source counts makes confusion a serious issue when probing deeper than ∼2–3 mJy. The number density of sources brighter than 8 mJy is a factor of ≃3 lower than at 3 mJy and so source confusion is much less of a problem in this survey. Finally, any sources discovered in this survey must be as bright or brighter in the submm than the extreme radio galaxies mentioned above. This survey is therefore optimally placed for detecting the most luminous starbursts with SFR≃1000 M yr−1, and through an extensive follow up program already underway (Fox et al. 2002; Lutz et al. 2001) will provide clues to the formation of the most massive ellipticals.

This paper is laid out as follows. In Sections 2 and 3 we describe the observations and data reduction. In Section 4 we outline the source extraction algorithm used to create our source catalogues. This is described in full mathematical detail in Appendix A. In Section 5 we present and discuss the results of simulations carried out both in conjunction with the real data, and through the production of fully-simulated survey areas. In Section 6 we present our source lists, giving positions, along with 850- and 450-μm flux densities. In Section 7 we derive the resulting source counts (both raw and corrected using the simulation results). In Section 8 we present two-point auto-correlation functions for each of our survey fields. A full discussion of the results, including comparison with other surveys, calculated dust masses, star formation rates, estimates of comoving number density, and clustering properties is given in Section 9. Finally, our conclusions are summarized in Section 10. In subsequent papers, the cross-correlation of our SCUBA sources with deep X-ray and radio observations will be discussed by Almaini et al. (2002) and Ivison et al. (in preparation) respectively.

For the calculation of physical quantities we have assumed H0=67 km s−1 Mpc−1 throughout, and all cosmological calculations have been performed for both an Einstein–de Sitter universe and a Λ-dominated flat cosmology with ΩM=0.3, ΩΛ=0.7.

Observations

The survey is divided between two fields; the Lockman Hole East (centred at RA 10h52m8s.82, Dec. +57°21′33.8) and ELAIS N2 (centred at RA 16h36m48s.85, Dec. +41°01′48.5). These areas both have a vast quantity of multiwavelength data available for follow-up studies, and were selected to coincide with deep Infrared Space Observatory (ISO) surveys at 6.7, 15, 90 and 175 μm (Lockman Hole East –Elbaz et al. 1999, Kawara et al. 1998; ELAIS N2 –Oliver et al. 2000, Serjeant et al. 2000, Efstathiou et al. 2000).

Pilot observations of the Lockman Hole East and ELAIS N2 fields began in 1998 March and July, respectively, using the Submillimetre Common User Bolometer Array (SCUBA; Holland et al. 1999) on the James Clerk Maxwell Telescope (JCMT) on Mauna Kea, Hawaii. In brief, SCUBA consists of two feedhorn-coupled bolometer arrays, both of which are arranged in a close-packed hexagon and have approximately a 2.3-arcmin field of view (slightly smaller on the short-wavelength array). The diffraction-limited beams delivered by the JCMT have a full width at half maximum (FWHM) of 14.5 arcsec at 850 μm and 7.5 arcsec at 450 μm, and the bolometer feedhorns on the arrays are sized for optimal coupling to the respective beams. As a result, the long-wavelength (LW) array consists of 37 bolometers, and the short-wavelength (SW) array consists of 91 bolometers. Observations may be made simultaneously with both arrays by means of a dichroic beamsplitter. A 3He/4He dilution refrigerator cools SCUBA down to ∼90 mK, making observations with the array sky-background limited.

The most efficient mapping technique for sources smaller than the field of view is ‘jiggle-mapping’, and is the method employed in this survey. As explained above, the bolometers are arranged such that the sky is under-sampled at any instant and so the secondary is ‘jiggled’ in a hexagonal pattern to fill in the gaps. In order to fully sample the sky at both wavelengths a 64-point jiggle-map is required. Temporal variations and linear spatial gradients in the sky were removed by chopping and nodding. A 30-arcsec chop throw was chosen in order to optimize the sky subtraction. The chop direction is fixed in celestial coordinates, lying directly north–south in the Lockman Hole East, and 48° east of north in ELAIS N2, parallel to the long axis of the survey area (to ensure the chop remains within the boundaries).

The preliminary observing strategy was to build up a linear strip in each field by overlapping SCUBA jiggle-maps, staggered by half the width of the array. The data collected from the pilot Lockman Hole East run in 1998 March were composed of a series of jiggle-maps, forming a deep strip at the centre of the survey area with σ850≃1.6 mJy beam−1. In subsequent observing runs, the strategy was modified to that of an overlapping ‘tripod’ positioning scheme, with three offset grids of hexagonally close-packed pointings, and the overall integration time reduced by a factor of ∼2 to make more efficient use of the allocated observing time. Each tripod position received 3.0 hr of integration time, yielding a typical noise level of σ850≃2.5 mJy beam−1. In order to give near homogeneous sky noise across the final coadded map, high and low airmasses were evenly distributed throughout. We have mapped a total area of 260 arcmin2, up to and including observations made in 2001 March.

Prior to 1999 December 15, the observations were made using the narrow-band 850-/450-μm filters. The majority of the observations carried out after this date made use of the more sensitive wide-band 850-/450-μm filters. Because the central wavelength of the 850-μm filter did not change by more than 1 μm, the data taken with the wide- and narrow-band filters were directly combined in producing the final images.

Skydips and pointings were carried out, on average, every 1.5–2 hr. The observations were carried out in dry weather conditions, under the constraint τ225 GHz>0.08. Primary/secondary calibration observations were taken at the start and end of every shift, with additional beam maps of the pointing source (0923+392 in the case of the Lockman Hole East, and 3C345 in the case of ELAIS N2) being made every 3–4 hr to provide tertiary flux-density calibration.

Data Reduction and Map Production

The data have been reduced using two independent software packages; surf (Jenness & Lightfoot 1997) and an IDL-based reduction package (Serjeant et al. 2002). The two mechanisms follow a very similar core reduction process, the main difference being in the production of the final maps. In addition to the output signal maps, the IDL routines automatically create corresponding uncorrelated noise maps by considering the signal variance. The final images produced by the surf and IDL reductions were found to be consistent to within a few per cent of each other. However, given the improved knowledge of the noise local to each source, the IDL mechanism was adopted as the primary reduction method and it is this process which we describe here.

After combining the positive and negative beams, the data were flat-fielded and then corrected for atmospheric extinction. On the nights when the Caltech Submillimetre Observatory (CSO) 225-GHz opacity meter made sufficiently good measurements around the time of a map, each jiggle pointing of the data was extinction-corrected by means of a polynomial fit to the 225-GHz opacity readings, interpolated to 850 and 450 μm by means of the relations in Archibald, Wagg & Jenness (2000b). If the CSO 225-GHz opacity meter was broken, or if the 225-GHz opacity measurements around the time of a jiggle map were poor, the extinction correction was assessed instead by a linear interpolation between consecutive skydips.

An iterative procedure was employed to deglitch the data sets and subtract any residual sky emission not removed by chopping and nodding. Each iteration made a time-dependent estimate of the noise for each bolometer, removing any spikes by means of a 3σ clip. A temporal modal sky level was determined by a fit to all bolometers in the array and subtracted from the data. With each consecutive iteration, the deglitching process makes a harder cut. Noisy bolometers were assigned a low inverse variance weight in this way.

To construct the final images the signal data were binned into 1-arcsec pixels, creating ‘zero-footprint’ maps with a corresponding noise value determined from the signal variance. The term ‘zero-footprint’ is an analogy with the drizzling algorithm (Fruchter & Hook 2002). A standard shift-and-add technique takes the flux in a given detector pixel and places its flux into the final map, over an area equivalent to one detector pixel projected on the sky. Drizzling, on the other hand, takes the flux and places it into a smaller area in the final map. Simulations have shown that this helps preserve information on small angular scales, provided that there are enough observations to fill in the resulting gaps. The area in the coadded map receiving the flux from one detector pixel is termed the footprint. Our method is an extreme example of drizzling; we take the data from each 14.5-arcsec bolometer and put the flux into a very small footprint (the ‘zero-footprint’), 1-arcsec square. Unlike in the standard surf reduction, there is no intrinsic smoothing or interpolation between neighbouring pixels in this re-bin procedure. Although there is some degree of correlation between pixels in the output zero-footprint signal maps in terms of the beam pattern, the corresponding pixel noise values represent individual measurements of the temporally varying sky noise averaged over the data set integration time, at a specific point on the sky, and are hence statistically independent from their neighbours. In essence we have a very oversampled image with statistically independent pixels. Statistical non-independence would refer to pixel-to-pixel crosstalk, which is not the case here.

A final 4σ clip on the signal-to-noise ratio was carried out to remove any remaining ‘hot pixels’. A noise-weighted convolution with a beam-sized Gaussian point spread function (PSF) produces realistic smoothed maps of the survey areas and is able to account for a variable signal-to-noise ratio between individual pixels.

As a result of the overlap of individual jiggle maps, and the fact that the observations were spread over 41 observing shifts, data missing as a result of bad bolometers are generally widely distributed, and only one region of the final 450-μm ELAIS N2 map has been noticeably affected. This is in the region of RA 16h36m48s.8 Dec.+41°02′44″, but is in fact a region which contains no significant 850-μm sources. Noise arising from undersampling by bolometers at the edge of the array has been evened out as far as possible by means of the tripod positioning scheme and only significantly affects the perimeter of the final mosaiced maps. Once these edges have been trimmed off, the resulting images have almost uniform noise across the central regions of the survey areas (σ850=2.5±0.7 mJy beam−1). The nature of the tripod positioning scheme does, however, lead to a noisier border of thickness ∼1/2 an array width surrounding the uniform noise region due to incomplete observations of one or more overlap positions, and hence less total integration time.

Calibration

Dunne & Eales (2001) have examined the issue of calibration errors in detail using the data from the ‘SCUBA Local Universe Galaxy Survey’ taken with the narrow-band 850-/450-μm filters. They find that an aperture calibration factor (ACF) is more accurate in determining source flux densities than the method of determining the flux conversion factor (FCF) by fitting to the calibrator peaks. At 850 μm the standard deviation on the ACF is ∼6 per cent, compared with ∼8 per cent when applying an FCF fitted to the peak. The improvement is particularly significant at 450 μm, where an aperture calibration method yields a typical standard deviation of ∼10 per cent compared with ∼18 per cent applying an FCF fitted to the peak. The greater variation in gain when fitting to the calibrator peaks arises because of the sensitivity of the peak flux density to changes in beam shape, due to sky noise, pointing drifts, chop throw, and the dish shape (particularly important at 450 μm).

Our IDL pipeline does, however, have the considerable advantage of producing uncorrelated noise maps, which have allowed us to apply a maximum-likelihood method to measure the statistical significance of each peak in our images, (see Section 4, and Appendix A). Our source-extraction algorithm provides a simultaneous fit to all peaks in the image, and is therefore able to decouple any partially confused sources whilst still recovering the additional signal-to-noise ratio yielded by the negative sidelobes. The resulting best-fitting model therefore allows us to measure the flux densities of our significant detections (<3.00σ) directly, as well as yielding properly quantified errors on this value. At 850 μm the accuracy gained by this approach more than compensates for the extra ∼2 per cent error in the calibration and so we have chosen to employ a calibration system based on fitting to the peak values. For consistency, we have adopted the same approach when considering the 450-μm data in this paper. The short-wavelength data will, however, be considered in more detail by Fox et al. (2002).

In order to calibrate the data, each shift was divided in two and all available information was utilized to calculate the mean flux-conversion factors applicable to each half shift by fitting to the calibrator peaks. Uranus and Mars were observed as primary calibrators when available, with the planetary flux densities taken from the JCMT fluxes program. Additionally, the secondary calibrators CRL618, OH231.8, CRL2688 and IRC10216 were used. As IRC10216 is variable with a period of 635 d, it was assumed to have constant flux density over an 8/9 night observing run, calculated using all other calibration data for that run. Beam maps of the pointing sources were used for tertiary calibration. A 30-arcsec chop throw was used for observing both the survey areas and the calibrators, in order to cancel out the effects of beam distortion as far as possible. Taking into account the errors on the calibrator flux densities, the variation in the flux-conversion factor over half a shift, and errors in the measured atmospheric opacity, the average calibration error over the survey maps was found to be 9 per cent at 850 μm, and 20 per cent at 450 μm.

Pointing problems

During the course of the survey, three different pointing problems are now known to have afflicted the JCMT.

First, prior to 1998 July 28, a bug in the chopping software of SCUBA meant that the chop direction was not updated once an observation had begun, and consequently the negative sidelobes accompanying a source in this data set will appear smeared and therefore shallower than expected. However, the effect is small when compared to the sky noise of ≃2.5 mJy beam−1 and affects only a small portion of our data (51/502 maps).

Secondly, some of the jiggle maps taken since 1999 June may have been affected by an elevation drive pointing glitch (Coulson 2000). Observations appear to suffer a 4-arcsec jump in the elevation when passing through transit. Additionally, pointing on the opposite side of the meridian to the survey field is likely to introduce a similar effect, of order 1.5–2 arcsec in magnitude. The fact that each area of sky is covered by 12 individual jiggle maps significantly reduces the likelihood of a large increase in positional uncertainty and source smearing. Under the most pessimistic assumptions, a maximum of 83/502 maps could be affected, of which a maximum of 3/12 cover any one area of sky.

Finally, data taken after 1999 July 15 were affected by a non-synchronization of the GPS and SCUBA data acquisition computer clocks. This problem has now been corrected for in the data reduction (Jenness 2000).

Source Extraction Technique

The chopping-nodding mechanism of the telescope provides a valuable method of discriminating between real detections and spurious noise spikes in the data. The 30-arcsec chop throw is small enough to fall onto the SCUBA array and so negative sidelobes, half the depth of the peak flux, are produced on either side of a real source, 30 arcsec away from the position of the central peak. This side lobe signal can be recovered to boost the overall signal-to-noise ratio of a detection.

For well-separated sources, convolving the images with the beam is formally the best method of source extraction (Eales et al. 1999, 2000; Serjeant et al. 2002). However, following a careful examination of the reduced data it became clear that some of the sources were partially confused, particularly in the map of ELAIS N2 where the negative sidelobes of individual sources had overlapped and were therefore somewhat deepened relative to both source peaks. Consequently, in order to decouple any confused sidelobes, a source-extraction algorithm was devised based on a simultaneous maximum-likelihood fit to the flux densities of all potentially significant peaks in the beam-convolved maps.

Using a peak-normalized beam map as a source template (CRL618 in the case of the Lockman Hole East, or Uranus in the case of ELAIS N2), a basic model was constructed by centring a beam map at the positions of all peaks <3 mJy as found in the Gaussian-convolved maps. The heights of each of the positioned beam maps were then calculated such that the final multisource model provided the best description of the submm sky, as judged by a minimum χ2 calculation made feasible by the independent data points and errors yielded by the zero-footprint IDL-reduced maps. A full mathematical description of this technique is provided in Appendix A.

The actual maps and source analysis are presented in Section 6. First, however, it is necessary to understand the selection effects that arise in an analysis of this sort.

Simulations

In order to assess the effects of confusion and noise on the reliability of our source-extraction algorithm, Monte Carlo simulations were carried out using 220 out of the 260 arcmin2 of data taken by the end of 2000 August. The dependences of positional error, completeness and error in reclaimed flux density on source flux density and noise in the maps were determined by planting individual sources of known flux density into the real ELAIS N2 and Lockman Hole East SCUBA maps. This has the advantage of allowing us to test the source-reclamation process against the real noise properties of the images. However, these simulations do not allow assessment of the level to which false or confused sources can contaminate an extracted source list. We have therefore also created a number of fully simulated images of the survey areas by assuming a reasonable 850-μm source-counts model. The results of analysing these two sets of simulations are discussed in the following two subsections.

Simulations building on the real survey data

A normalized beam map was used as a source template. At 0.5-mJy flux intervals, from 5 to 15 mJy, individual sources were added into the unconvolved zero-footprint signal maps. This was done one fake source at a time, so as not to enhance significantly any existing real confusion noise within the image. The source-extraction algorithm was then re-run. This exercise was repeated for 100 different randomly selected positions on each image, at each flux level, so that source reclamation could be monitored as a function of input flux and position/noise level within the maps. Source reclamation was deemed to have been successful if the source-extraction algorithm returned a source with a signal-to-noise ratio (S/N)<3.50 within less than half a beamwidth of the input position. When analysing the output data, ≃10 per cent of the input sources were found to lie within half a beam width of another peak <5 mJy already present in the image. Any such source was excluded from statistical analyses of positional error, completeness and error in reclaimed flux density, as the proportion of blends produced by two sources, both with comparably bright fluxes in the region of 8 mJy, should be negligible within our sample.

Figs 1–4 show the dependences of mean positional error, percentage of sources returned, systematic offset in returned flux density and predicted output against input flux density for both regions of uniform and non-uniform noise. The error bars show the standard error on the mean. The solid lines in Fig. 1, and the solid curves in Figs 2 and 3, are best-fitting models to the plotted data points. From Fig. 3, we see that the effect of noise and confusion is to produce systematic ‘flux boosting’, the retrieved flux density always being greater than the input value. The solid curves in Fig. 4 show the predicted output versus input flux density based on the fitted flux-density boosting models of Fig. 3. We find that there is a strong correlation with input flux density for each of these quantities, the scatter being somewhat larger in the areas of non-uniform noise lying towards the edge, as would be expected. The positional error and flux-density boosting factors decrease with increasing input flux density, whilst the percentage of sources returned at <3.50σ increases. The effect of additional noise in the non-uniform regions considerably increases the uncertainty in positional error, as well as the mean flux-density boosting factor. It also serves to typically half the percentage of sources reclaimed. A summary of the results at 5, 8, 11 and 14 mJy is given in Table 1. The integral completeness down to each flux-density limit has been calculated by a counts slope-weighted sum of the fraction of sources returned at each specific flux-density level, based on an adopted 850-μm source counts model consistent with the measured source counts, given in Section 5.2. As detailed in Table 1, the typical completeness at 8 mJy is ≃80 per cent within the uniform noise regions of the images, with flux-density boosting factors of ≃15 per cent. By 11 mJy, completeness has risen to ≃90 per cent, while flux-density boosting has dropped to >10 per cent.

Figure 1.

The dependence of positional error on input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and then retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid lines are a best fit to the plotted data points.

Figure 1.

The dependence of positional error on input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and then retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid lines are a best fit to the plotted data points.

Figure 2.

The percentage of sources returned against input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid curves show a best-fitting exponential-based model to the plotted data points.

Figure 2.

The percentage of sources returned against input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid curves show a best-fitting exponential-based model to the plotted data points.

Figure 3.

The factor by which the input source flux is boosted (when retrieved from the submm map) plotted against input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid curves show a best-fit exponential-based model to the plotted data points.

Figure 3.

The factor by which the input source flux is boosted (when retrieved from the submm map) plotted against input source flux density in the range 5–15 mJy, based on a single template source repeatedly being added into and retrieved (<3.50σ) from the real unconvolved signal maps. The error bars give the standard error on the mean. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The solid curves show a best-fit exponential-based model to the plotted data points.

Figure 4.

Based on the fitted flux-density boosting curves shown in the preceding figure, these plots show the predicted output flux density against input flux density in the range 5–15 mJy. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The dashed line in both cases delineates the ideal case of output flux density equal to input flux density.

Figure 4.

Based on the fitted flux-density boosting curves shown in the preceding figure, these plots show the predicted output flux density against input flux density in the range 5–15 mJy. Plot (a) gives the results from the uniform noise region, whereas plot (b) gives the results for the non-uniform edge regions which have not received the full integration time. The dashed line in both cases delineates the ideal case of output flux density equal to input flux density.

Table 1.

Results of the completeness simulations at the 3.50σ significance level, in which individual sources of the specified flux density given in column 2 were added into the unconvolved survey data taken up to and including 2000 August 16, and retrieved by means of the source extraction algorithm outlined in Section 4.

Table 1.

Results of the completeness simulations at the 3.50σ significance level, in which individual sources of the specified flux density given in column 2 were added into the unconvolved survey data taken up to and including 2000 August 16, and retrieved by means of the source extraction algorithm outlined in Section 4.

Completely simulated maps

Adopting an 850-μm source counts model which is consistent with the measured point source counts, and which does not over-predict the FIR background, we have created fully simulated versions of our two survey fields. In order to produce a realistic model of the background counts, the number of sources expected at 0.1-mJy intervals was Gaussian randomized and then the appropriate sources were planted into the simulated data sets at random positions. The whole image was then convolved with the beam. At all flux-density levels the sources were assumed to be unclustered (although see Section 9.4 and Almaini et al. 2002).

We then constructed a pure noise image for each of our two fields by scrambling the astrometry of the individual SCUBA observations used to produce the maps. These may be seen in Figs 5 and 7. After creating the scrambled noise maps we first ran the source-extraction algorithm on these to assess the frequency with which spurious sources would artificially appear on the basis of random statistics. The noise images yielded only a single artificial source at <3.50σ, and a total of 12 spurious sources at <3.00σ, which is not unexpected given the ≃1200 beams across the areas of these maps. Adding the simulated background counts to the scrambled maps produces the final simulated images. The source-extraction algorithm was applied as previously.

Figure 5.

A pure noise image of the ELAIS N2 region, created from the real data by corrupting the astrometry of the individual jiggle maps to produce a ‘scrambled’ map. This image includes the jiggle maps taken up to and including 2000 August 16. There are no source detections at <3.50σ.

Figure 5.

A pure noise image of the ELAIS N2 region, created from the real data by corrupting the astrometry of the individual jiggle maps to produce a ‘scrambled’ map. This image includes the jiggle maps taken up to and including 2000 August 16. There are no source detections at <3.50σ.

Figure 7.

A pure noise image of the Lockman Hole East region, created from the real data by corrupting the astrometry of the individual jiggle maps to produce a ‘scrambled’ map. This image includes the jiggle maps taken up to and including 2000 August 16. There is one spurious source detection at <3.50σ, which is circled on the image.

Figure 7.

A pure noise image of the Lockman Hole East region, created from the real data by corrupting the astrometry of the individual jiggle maps to produce a ‘scrambled’ map. This image includes the jiggle maps taken up to and including 2000 August 16. There is one spurious source detection at <3.50σ, which is circled on the image.

Four realisations of the sky in each field were generated for the source counts model  

(1)
formula
where N0=1.1×104 deg−2 mJy−1, S0=1.79 mJy, α=0.0 and β=3.13. Examples of the simulated ELAIS N2 and Lockman Hole East maps for our adopted counts model may be seen in Figs 6 and 8.

Figure 6.

Fully simulated ELAIS N2 SCUBA 850 μm map, using the adopted counts model given in Section 5.2. The image noise was simulated using the scrambled map shown in Fig. 5.

Figure 6.

Fully simulated ELAIS N2 SCUBA 850 μm map, using the adopted counts model given in Section 5.2. The image noise was simulated using the scrambled map shown in Fig. 5.

Figure 8.

Fully simulated Lockman Hole East SCUBA 850-μm map, using the adopted counts model given in Section 5.2. The image noise was simulated using the scrambled map shown in Fig. 7.

Figure 8.

Fully simulated Lockman Hole East SCUBA 850-μm map, using the adopted counts model given in Section 5.2. The image noise was simulated using the scrambled map shown in Fig. 7.

Tables 2–4 summarize the results at 4.00, 3.50 and 3.00σ, respectively, for the simulated fields at 8 mJy. In these tables and the discussion that follows, we have adopted the practical criterion that any source retrieved at the 8-mJy level which ultimately could not be identified with a single source brighter than 5 mJy (and therefore could not be confirmed through interferometric follow-up, for example), should be regarded as spurious or the result of confusion. In the uniform noise regions, contamination from spurious sources or confused sources >5 mJy is ≃40 per cent at S/N<3.00, dropping to ≃20 per cent at S/N<3.50 and >10 per cent at S/N<4.00. The completeness at a measured flux density of 8 mJy is ≃85 per cent within the regions of uniform noise, consistent with the completeness determined from the simulations using the real survey data. We also find that although the real source counts (i.e. total number of sources out minus the number of spurious sources) at <3.50σ agree well with the number planted into the simulations at 8 mJy, this is partly by chance; the number of sources boosted from a lower flux-density level, between 5 and 8 mJy, is approximately the same as the number of real sources <8 mJy which are not retrieved by the source-extraction algorithm. By lowering our S/N threshold from 3.50 to 3.00, there is very little improvement (if any) in completeness, however the proportion of spurious sources detected is much more serious. These results provide a clear motivation for setting a S/N threshold of <3.50 when determining the source number counts, as opposed to the <3.00σ level applied in some other surveys. Our results are consistent with those of Hughes & Gaztañaga (2000), who have simulated submm galaxy surveys for a range of wavelengths, spatial resolutions and flux densities.

Table 2.

Results of the completeness simulations at the 4.00σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximum-likelihood estimator outlined in Section 4.

Table 2.

Results of the completeness simulations at the 4.00σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximum-likelihood estimator outlined in Section 4.

Table 3.

Results of the completeness simulations at the 3.50σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximumlikelihood estimator outlined in Section 4.

Table 3.

Results of the completeness simulations at the 3.50σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximumlikelihood estimator outlined in Section 4.

Table 4.

Results of the completeness simulations at the 3.00σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximumlikelihood estimator outlined in Section 4.

Table 4.

Results of the completeness simulations at the 3.00σ significance level, in which each field was fully simulated four times using an adopted counts model consistent with the measured source counts. The sources were then retrieved by means of the maximumlikelihood estimator outlined in Section 4.

Sources

Sources down to a significance level of <3.00σ at 850 μm are listed in Table 5 for ELAIS N2 and Table 6 for the Lockman Hole East, ranked in order of formal significance as determined from our maximum-likelihood extraction technique. The top section of each table, containing 7 and 12 sources in ELAIS N2 and Lockman Hole East, respectively, is confined to sources which have S/N<4.00, and the middle section to sources which have S/N<3.50. As our simulations indicate that 70–80 per cent of these sources are expected to be real individual sources with flux densities within ≃25 per cent of the extracted flux density, we recommend that follow-up observations are confined to these abbreviated S/N<3.50 source lists. The 850-μm flux densities in column 4 are the best-fitting scaling factors determined by our extraction algorithm, with the quoted errors also incorporating the 9 per cent calibration error.

Table 5.

Sources retrieved in the long-wavelength ELAIS N2 image above the 3.00σ significance level using the source-extraction algorithm outlined in Section 4. We used a peak-normalized calibrator beam-map as a source template, and constructed a basic model by centring a beam map at the positions of all peaks >3 mJy as found in the 14.5-arcsec FWHM Gaussian convolved maps. The heights of each of the positioned beam maps were allowed to vary simultaneously such that the final model satisfied a minimized χ2 fit.We suggest that follow-up observations be confined to those detections >3.50σ given the serious problem of contamination by spurious sources at the 3.00σ level. Calibration errors of 9 per cent and 20 per cent at 850μm and 450μm, respectively, are included in the error on the flux densities.

Table 5.

Sources retrieved in the long-wavelength ELAIS N2 image above the 3.00σ significance level using the source-extraction algorithm outlined in Section 4. We used a peak-normalized calibrator beam-map as a source template, and constructed a basic model by centring a beam map at the positions of all peaks >3 mJy as found in the 14.5-arcsec FWHM Gaussian convolved maps. The heights of each of the positioned beam maps were allowed to vary simultaneously such that the final model satisfied a minimized χ2 fit.We suggest that follow-up observations be confined to those detections >3.50σ given the serious problem of contamination by spurious sources at the 3.00σ level. Calibration errors of 9 per cent and 20 per cent at 850μm and 450μm, respectively, are included in the error on the flux densities.

Table 6.

Sources retrieved in the Lockman Hole East field above the 3.00σ significance level using the source extraction algorithm outlined in Section 4.We used a peak-normalized calibrator beam-map as a source template, and constructed a basic model by centring a beam-map at the positions of all peaks >3 mJy as found in the 14.5-arcsec FWHM Gaussian convolved maps. The heights of each of the positioned beam maps were allowed to vary simultaneously such that the final model satisfied a minimized χ2 fit. We suggest that follow-up observations are confined to those detections >3.50σ given the serious problem of contamination by spurious sources at the 3.00σ level. Calibration errors of 9 per cent and 20 per cent at 850μm and 450μm, respectively, are included in the error on the flux densities.

Table 6.

Sources retrieved in the Lockman Hole East field above the 3.00σ significance level using the source extraction algorithm outlined in Section 4.We used a peak-normalized calibrator beam-map as a source template, and constructed a basic model by centring a beam-map at the positions of all peaks >3 mJy as found in the 14.5-arcsec FWHM Gaussian convolved maps. The heights of each of the positioned beam maps were allowed to vary simultaneously such that the final model satisfied a minimized χ2 fit. We suggest that follow-up observations are confined to those detections >3.50σ given the serious problem of contamination by spurious sources at the 3.00σ level. Calibration errors of 9 per cent and 20 per cent at 850μm and 450μm, respectively, are included in the error on the flux densities.

Figs 9 and 10 show the beam-sized Gaussian convolved survey maps at 850 μm for ELAIS N2 and the Lockman Hole East, respectively, with the sources <3.50σ circled and numbered. We have also applied our maximum-likelihood estimator to the 450-μm data (Fox et al. 2002), centring normalized 450-μm beam maps at the same positions as in the 850-μm images. Column 7 gives the 450-μm flux densities if there is a detection in the short-wavelength data with S/N<3.00σ, or, more commonly, a 3σ upper limit. We also carried out a search for peaks in the short-wavelength data lying within a search radius of 8 arcsec from the 850-μm source peaks (the combined result of adding in quadrature the 1/2 beamwidth positional uncertainties of the long- and short-wavelength data) and again employed our source-extraction algorithm. The flux densities of any detections with S/N better than 3.00σ are given in column 9, with their significances and positional offset from the long-wavelength position in columns 10 and 11. The errors quoted on the 450-μm flux densities also include the 20 per cent calibration error.

Figure 9.

The 850-μm image of the ELAIS N2 field, smoothed with a beam-size Gaussian (14.5-arcsec FWHM). The numbered circles highlight those sources found at a significance of <3.50, the labels corresponding to the numbers in Table 5.

Figure 9.

The 850-μm image of the ELAIS N2 field, smoothed with a beam-size Gaussian (14.5-arcsec FWHM). The numbered circles highlight those sources found at a significance of <3.50, the labels corresponding to the numbers in Table 5.

Figure 10.

The 850-μm image of the Lockman Hole East field, smoothed with a beam-size Gaussian (14.5-arcsec FWHM). The numbered circles highlight those sources found at a significance of <3.50, the labels corresponding to the numbers in Table 6.

Figure 10.

The 850-μm image of the Lockman Hole East field, smoothed with a beam-size Gaussian (14.5-arcsec FWHM). The numbered circles highlight those sources found at a significance of <3.50, the labels corresponding to the numbers in Table 6.

Given that the FWHM of the JCMT beam at 850 μm is 14.5 arcsec, the source positions are only likely to be accurate to within ≃3–4 arcsec. Consistent with this, our simulations suggest that an 8-mJy source, detected at <3.50σ, should have accurate astrometry to ∼3.5 arcsec in uniform noise regions and ∼4 arcsec in the non-uniform noise areas, the positional uncertainty decreasing with increasing flux. However, these simulations do not account for the elevation drive pointing glitch mentioned in Section 3.2, and so it is possible that the astrometric accuracy of a source affected by this problem could be somewhat worse than this.

Analysis of our fully-simulated survey areas suggests that of the 19 sources with S/N<4.00, only one is likely to be spurious. At the 3.50σ threshold we have 38 sources, of which ≃7 may be spurious. The situation is markedly worse at 3.00σ, where ≃30 of the supposedly significant peaks may in fact be due to noise or confusion of low flux-density sources.

In the short-wavelength data, we can find only three detections at <4.00σ within an 8-arcsec radius of the 850-μm positions. These 450-μm detections appear to be associated with the most significant Lockman Hole East 850-μm source and the two most significant ELAIS N2 850-μm sources. There are additionally another five 450-μm detections at better than 3.00σ, of which all but one is associated with a <3.50σ 850-μm source. We must be cautious when considering the credibility of the lower significance (<3.00σ) 450-μm detections. The atmospheric opacity at 450 μm is much more sensitive to changes in the amount of water vapour present than at 850 μm, and so the resulting short-wavelength maps are noisier than their long-wavelength counterparts. Furthermore, the smaller beam size (7.5-arcsec FWHM at 450 μm as opposed to 14.5-arcsec FWHM at 850 μm) means that the short-wavelength maps will contain a greater number of spurious detections. For S/N<3.00 we would expect 53 spurious sources across the full 260 arcmin2 of short-wavelength data (and 14 spurious detections in the long-wavelength data) based on Gaussian statistics alone.

A detailed discussion of individual sources may be found in the second paper of the SCUBA 8-mJy survey series (Fox et al. 2002), along with a multiwavelength analysis considering possible identifications and SED-based redshift constraints. The 450-μm images of the survey fields are also presented by Fox et al. (2002).

Number Counts

In order to obtain the best estimate of the number counts at 850 μm, we must account for the effects of flux-density boosting, incompleteness and contamination from spurious sources. Using the results of our simulations carried out on the real survey data, the observed flux densities of the <3.50σ sources were corrected statistically using the output flux density versus input flux density plots shown in Fig. 4, for the appropriate image noise region. The number of boost-corrected sources in bins 0.5-mJy wide, centred on 5.0, 5.5, 6.0 and so on up to 15 mJy, were then amended for completeness using the best-fitting results of the real-data simulations shown in Fig. 2. The integrated source counts were obtained by summing the number of sources in the bins down to the required threshold flux. As the scatter in the simulation results is substantially higher in the non-uniform noise regions (reflected by reduced χ2 values of the plotted best-fitting curves being a factor of ∼5 higher in the edge regions), we have based our estimation of the 850-μm source counts purely on the sources detected in the areas of uniform noise, covering an area of 190 arcmin2. In practice this means that 5/38 of the sources above the 3.50σ significance threshold have been excluded. This process is summarized in Table 7, with the raw and amended counts per square degree shown in Table 8, along with Poisson errors.

Table 7.

A breakdown of the raw and corrected source counts at each 0.5-mJy interval bin for the sources with S/N > 3:50 from both survey maps, excluding those uncovered in the non-uniform noise regions. Column 1 gives the flux density, and columns 2 and 3 the raw number counts in the 0.5-mJy wide bin and cumulative to that flux-density level. Column 4 gives the number in each bin once flux-density boosting has been accounted for using the fitted curve in Fig 4(a), and similarly column 5 gives the number in each bin once completeness has also been taken into account using the fitted curve in Fig 2(a). Column 6 gives the cumulative corrected counts down to each fluxdensitylevel.

Table 7.

A breakdown of the raw and corrected source counts at each 0.5-mJy interval bin for the sources with S/N > 3:50 from both survey maps, excluding those uncovered in the non-uniform noise regions. Column 1 gives the flux density, and columns 2 and 3 the raw number counts in the 0.5-mJy wide bin and cumulative to that flux-density level. Column 4 gives the number in each bin once flux-density boosting has been accounted for using the fitted curve in Fig 4(a), and similarly column 5 gives the number in each bin once completeness has also been taken into account using the fitted curve in Fig 2(a). Column 6 gives the cumulative corrected counts down to each fluxdensitylevel.

Table 8.

The 850-μm source counts per square degree based on sources with S/N<3.50 in both survey maps, and excluding those detected in the non-uniform noise regions. Column 1 gives the flux density and column 2 the cumulative raw counts per square degree with the Poisson error. Column 3 gives the cumulative corrected counts per square degree, the upper error corresponding to the Poisson error, and the lower error accounting for both the Poisson error and the presence of spurious sources based on the simulation data.

Table 8.

The 850-μm source counts per square degree based on sources with S/N<3.50 in both survey maps, and excluding those detected in the non-uniform noise regions. Column 1 gives the flux density and column 2 the cumulative raw counts per square degree with the Poisson error. Column 3 gives the cumulative corrected counts per square degree, the upper error corresponding to the Poisson error, and the lower error accounting for both the Poisson error and the presence of spurious sources based on the simulation data.

The applied correction factors have been determined from simulations carried out on the real images and so are independent of source-count models. We cannot, however, attempt to correct for the presence of random detections without making some assumptions about the background counts. Consequently, the estimated number of spurious sources is included in the lower error bars combined in quadrature with the Poisson error, rather than being applied directly to the actual tabulated source-count values.

Fig. 11 shows our corrected source counts, at 1-mJy intervals, from 5 to 12 mJy, along with data points from other SCUBA surveys (Smail et al. 1997; Hughes et al 1998; Barger et al. 1999; Blain et al. 1999; Eales et al. 2000). The upper solid and dashed lines in Fig. 11 are source-count models, derived by interpolating the 60-μm luminosity function (Saunders et al. 1990) to 850 μm, with an assumed dust temperature of 40 K and dust emissivity index β=1.2 (consistent with measurements of local ULIRGs such as Arp 220, for example Dunne et al. 2000), and assuming pure luminosity evolution of the form L(z)=L(0)(1+z)3 to z=2.0 (beyond which the luminosity function is simply frozen). The solid line assumes an Einstein–de Sitter cosmology, whereas the dashed line assumes the cosmology ΩM=0.3 and ΩΛ=0.7 (both adopting H0=67 km s−1 Mpc−1). The lower solid and dashed lines represent the same interpolated 60-μm luminosity function with no luminosity evolution included, again for cosmologies ΩM=1.0, ΩΛ=0.0, and ΩM=0.3, ΩΛ=0.7, respectively. The dotted line shows the adopted source counts model, which was assumed when creating the fully simulated images. The dot-dashed and dot-dot-dashed lines represent a more physical form of luminosity evolution (Jameson 2000; Smail et al. 2001), which is compatible with models of cosmic chemical evolution and includes a peak in the evolution function:  

(2)
formula
where b=2.2±0.1 and c=1.84±0.1. The dot-dashed line assumes a cosmology ΩM=1.0, ΩΛ=0.0, and the dot-dot-dashed line ΩM=0.3, ΩΛ=0.7. A dust emissivity index of β=1.2 and a dust temperature of 37 K were also assumed (Saunders et al. (1990)).

Figure 11.

Cumulative 850-μm source counts. The solid diamonds show the new results from the ‘SCUBA 8-mJy Survey’ derived in this paper. The upper solid and dashed curves are predicted number counts, based on the 60-μm luminosity function of Jameson 2000; Smail et al. 2001, assuming a dust temperature of 40 K, a dust emissivity index of β=1.2, and pure luminosity evolution of the form L(z)=L(0)(1+z)3 out to z=2 (beyond which the luminosity function is simply frozen), for cosmologies (ΩM=1.0,ΩΛ=0.0) and (ΩM=0.3,ΩΛ=0.7) respectively. The lower solid and dashed curves show the number counts predicted by the same models if no luminosity evolution is included. The dotted line shows the counts model assumed when creating simulated images. The dot-dash and dot-dot dash assume a more complicated, but arguably more realistic luminosity evolution of the form L(z)=L(0)(1+z)3/2sech2[b ln(1+z)-c]cosh2c, for cosmologies (ΩM=1.0,ΩΛ=0.0) and (ΩM=0.3,ΩΛ=0.7), respectively, and assuming a dust temperature of 37 K and dust emissivity index of β=1.2 (Smail et al. 2001).

Figure 11.

Cumulative 850-μm source counts. The solid diamonds show the new results from the ‘SCUBA 8-mJy Survey’ derived in this paper. The upper solid and dashed curves are predicted number counts, based on the 60-μm luminosity function of Jameson 2000; Smail et al. 2001, assuming a dust temperature of 40 K, a dust emissivity index of β=1.2, and pure luminosity evolution of the form L(z)=L(0)(1+z)3 out to z=2 (beyond which the luminosity function is simply frozen), for cosmologies (ΩM=1.0,ΩΛ=0.0) and (ΩM=0.3,ΩΛ=0.7) respectively. The lower solid and dashed curves show the number counts predicted by the same models if no luminosity evolution is included. The dotted line shows the counts model assumed when creating simulated images. The dot-dash and dot-dot dash assume a more complicated, but arguably more realistic luminosity evolution of the form L(z)=L(0)(1+z)3/2sech2[b ln(1+z)-c]cosh2c, for cosmologies (ΩM=1.0,ΩΛ=0.0) and (ΩM=0.3,ΩΛ=0.7), respectively, and assuming a dust temperature of 37 K and dust emissivity index of β=1.2 (Smail et al. 2001).

Two-Point Correlation Functions

In order to obtain a first measure of the clustering properties of each of our two fields, we have used our <3.50σ sources to calculate two-point angular correlation functions for both ELAIS N2 and the Lockman Hole East. We created a catalogue for each survey area, containing ∼40 000 fake sources placed according to a Poisson distribution within the boundaries of the real data. Although the positions were generated at random, we used our uncorrelated noise maps to weight the number density of sources across the image, as we would expect to uncover a larger density of sources above a specified signal-to-noise ratio threshold in areas of lower noise. In practise, we divided the full image into a series of sub-images, 20 by 20 arcsec in size, and calculated the mean noise level in each of these grid sections. For sources brighter than ≃5 mJy, the number of sources expected above a constant signal-to-noise ratio threshold will increase approximately as N(<S)∝S−1.5, and consequently we would expect the relative number of sources in each sub-image to scale as 〈noise−1.5.

The angular correlation function w(θ) describes the excess probability of finding a source at an angular separation θ from a source selected at random, over a random distribution. There are a variety of possible estimators for w(θ) as a function of pair-count ratios. Here we have used:  

(3)
formula
where there are n and m objects in the real and mock catalogues respectively, DD is the number of distinct data pairs in the real image within a bin covering a specified range of θ, and DR is the number of cross-pairs between the real and fake catalogues within the same range of θ.

The angular correlation functions may be seen with Poisson errors in Figs 12 and 13, where a bin size of twice the beam (29.0 arcsec) has been used.

Figure 12.

Two-point angular correlation function for the ELAIS N2 survey field. 1+w(θ)=2mDD〉/[(n-1)〈DR〉] where there are n and m real and fake <3.50σ sources, respectively, and DD and DR are the number of distinct data pairs and cross data random pairs within 29-arcsec wide bins as a function of angular separation θ. The power-law line indicates the correlation function found by Daddi et al. (2000) for extremely red objects (EROs) with R-K<5 and K>18.5 as discussed in Sections 9.3 and 9.4.

Figure 12.

Two-point angular correlation function for the ELAIS N2 survey field. 1+w(θ)=2mDD〉/[(n-1)〈DR〉] where there are n and m real and fake <3.50σ sources, respectively, and DD and DR are the number of distinct data pairs and cross data random pairs within 29-arcsec wide bins as a function of angular separation θ. The power-law line indicates the correlation function found by Daddi et al. (2000) for extremely red objects (EROs) with R-K<5 and K>18.5 as discussed in Sections 9.3 and 9.4.

Figure 13.

Two-point angular correlation function for the Lockman Hole East survey field. 1+w(θ)=2mDD〉/[(n-1)〈DR〉] where there are n and m real and fake <3.50σ sources, respectively, and DD and DR are the number of distinct data pairs and cross data random pairs within 29-arcsec wide bins as a function of angular separation θ. The power-law line indicates the correlation function found by Daddi et al. (2000) for EROs with R-K<5 and K>18.5 as discussed in Sections 9.3 and 9.4. Note that the vertical dashed line demonstrates the cut-off at which there are no more source pairings in the real data.

Figure 13.

Two-point angular correlation function for the Lockman Hole East survey field. 1+w(θ)=2mDD〉/[(n-1)〈DR〉] where there are n and m real and fake <3.50σ sources, respectively, and DD and DR are the number of distinct data pairs and cross data random pairs within 29-arcsec wide bins as a function of angular separation θ. The power-law line indicates the correlation function found by Daddi et al. (2000) for EROs with R-K<5 and K>18.5 as discussed in Sections 9.3 and 9.4. Note that the vertical dashed line demonstrates the cut-off at which there are no more source pairings in the real data.

Discussion

Comparative source counts and models

Although our survey was undertaken with the aim of constraining the brighter end of the 850-μm number counts, we have in fact achieved the most accurate determination to date of the submm source counts down to S850≃5 mJy. From Fig. 11 we see that at 5 mJy there is good agreement between our counts graphic and the values of Eales et al. (2000)(500±200 deg−2), and Barger et al. (1999) graphic, to well within the quoted errors. However, at larger flux densities, the number of sources per square degree in both the Eales and Barger survey areas is somewhat lower than our counts. For example, at 6 mJy the number count of Eales et al. is a factor of ∼3 lower than our value, and that of Barger et al. (1999) is lower by a factor of ∼2.5. The most likely explanation for this discrepancy is that it is due to small-number statistics rather than a real steepening of the source counts at <6 mJy. We would in fact expect this to be the case if SCUBA sources were generally clustered on scales of a few arcmin (as indicated by our two-point autocorrelation functions, albeit with large errors, shown in Section 8 and discussed in Section 9.4). The ‘SCUBA 8-mJy Survey’ covers 260 arcmin2 of sky, five times the area of the ‘Canada UK Deep Submillimetre Survey (CUDSS) 14h Field’ of Eales et al., and ∼2.5 times the area of Barger et al., making us less susceptible to arcmin-scale clustering and thus leaving us much better able to constrain the number counts at brighter flux densities. At 8 mJy, Blain et al. (1999) suggest a rather higher count of 800±600 deg−2, compared to our value of graphic, however the large error on the Blain et al. number count renders their result almost meaningless. At flux densities <10 mJy, our survey also inevitably becomes limited by the lower surface density of brighter sources. In a similar manner to the source counts of Eales et al. and Barger et al. beyond 6 mJy, our counts at 11 and 12 mJy suggest a steepening of the counts slope towards brighter fluxes. Borys et al. (2002) report a number count of graphic at 12 mJy, suggesting that this is again a statistical rather than a real effect.

The source counts accumulated across all the submm surveys (as shown in Fig. 11) are consistent with a population of high-redshift (z<1), dusty, star-forming galaxies, analogous to the local ULIRGs, with dust temperatures Td∼40 K, emissivity β∼1.2, and pure luminosity evolution of the form (1+z)3 out to a redshift z=2 (and simply frozen beyond). We also find that our source counts are in very good agreement with the more physical form of luminosity evolution of Jameson (2000) and Smail et al. (2001), assuming a dust temperature of 37 K, particularly for the cosmology ΩM=0.3, ΩΛ=0.7. The lower solid and dashed lines show the predicted source counts using the same luminosity function but without any luminosity evolution. Such a model is inconsistent with the measured number counts by 2–3 orders of magnitude, which simply serves to re-emphasize the extent of the evolution in the starburst population uncovered by SCUBA.

Dust masses and star formation rates

The photometric redshift estimates of the 19 sources discussed in Fox et al. (2002) and Lutz et al. (2001) are all indicative of redshifts z<1, and in the vast majority of cases favour z<2. Even without knowing a precise redshift, this fact coupled with an almost uniform sensitivity to objects between redshifts z=1 and z=8 allows us to make an estimate of the dust-enshrouded SFRs and dust masses. Assuming a dust temperature of 40 K and emissivity index β=1.2, we have interpolated the inferred 850-μm luminosity to 60 μm using an optically thin greybody approximation. The FIR luminosity provides a measure of the current SFR of massive stars. In regions of intense star formation, dust is heated primarily by the embedded O and B stars which evolve rapidly and dispense their surrounding material on similarly short time-scales (∼107 yr, Wang 1991), hence the rate of dust production is proportional to the star formation rate. We may calculate the SFR to within a factor of a few, by means of  

(4)
formula

The value of ε is uncertain to within a factor of ∼3 (Scoville & Young 1983; Thronson & Telesco 1986; Rowan-Robinson et al. 1997; Rowan-Robinson 2000), the main sources of error arising from uncertainties in the initial mass function (IMF), the fraction of optical/UV light absorbed by the dust, and the time-scale of the burst. We have assumed a value of ε=2.1 (Thronson & Telesco 1986), which incorporates a burst of O, B and A type star formation over ∼2×106 yr and assumes a Salpeter IMF. It was also assumed that all of the optical/UV radiation was absorbed and re-radiated thermally by the surrounding dust.

In calculating the mass of dust present in each source, we assumed the submm continuum to be the result of optically-thin thermal emission from the heated dust grains, with no additional contribution from bremsstrahlung or synchrotron radiation. We may then determine the dust mass Md by means of the relation  

(5)
formula
where z is the redshift of the source, Sobs is the observed flux density, DL is the luminosity distance, graphic is the rest-frequency mass absorption coefficient, and B(νrest, Td) is the rest-frequency value of the Planck function from dust grains radiating at temperature Td. The main uncertainty in determining the dust mass Md is the uncertainty in the rest-frequency mass absorption coefficient graphic. We have adopted the same approach as Hughes et al. (1997), using their average value of kd(800 μm)=0.15±0.09 m2 kg−1, and interpolating to shorter wavelengths by assuming kdλ−1.5. A different choice of kd(800 μm) would be expected to change the dust mass estimates by a factor of up to 2.

The results are outlined in Tables 9 and 10, the unbracketed quantities assuming an Einstein–de Sitter cosmology and the values enclosed in brackets adopting ΩM=0.3, ΩΛ=0.7. SFRs range from several hundred to several thousand solar masses per year (exceeding even that of the extreme local starburst Arp 220), the star forming activity being heavily obscured by 108–109 M of dust. In fact, for ΩM=0.3, ΩΛ=0.7 and H0=67 km s−1 Mpc−1, it can be seen that an 850-μm flux density of 8 mJy corresponds closely to a star formation rate of 1000 M yr−1.

Table 9.

Inferred FIR/60-μm luminosities, star-formation rates and dust masses for the sources in ELAIS N2. The unbracketed numbers in columns 7, 8 and 9 assume an Einstein-de Sitter cosmology, whereas the numbers enclosed in brackets assume ΩM = 0.3 and ΩΛ = 0.7. A dust temperature of 40 K, and emissivity index β = 1.2 were assumed throughout, as was H0 = 67 km s−1 Mpc−1.

Table 9.

Inferred FIR/60-μm luminosities, star-formation rates and dust masses for the sources in ELAIS N2. The unbracketed numbers in columns 7, 8 and 9 assume an Einstein-de Sitter cosmology, whereas the numbers enclosed in brackets assume ΩM = 0.3 and ΩΛ = 0.7. A dust temperature of 40 K, and emissivity index β = 1.2 were assumed throughout, as was H0 = 67 km s−1 Mpc−1.

Table 10.

Inferred FIR/60-μm luminosities, star-formation rates and dust masses for the sources in the Lockman Hole East. The unbracketed numbers in columns 7, 8 and 9 assume an Einstein-de Sitter cosmology, whereas the numbers enclosed in brackets assume ΩM = 0.3 and ΩΛ = 0.7. A dust temperature of 40 K, and emissivity index β = 1.2 were assumed throughout, as was H0 = 67 km s−1 Mpc−1.

Table 10.

Inferred FIR/60-μm luminosities, star-formation rates and dust masses for the sources in the Lockman Hole East. The unbracketed numbers in columns 7, 8 and 9 assume an Einstein-de Sitter cosmology, whereas the numbers enclosed in brackets assume ΩM = 0.3 and ΩΛ = 0.7. A dust temperature of 40 K, and emissivity index β = 1.2 were assumed throughout, as was H0 = 67 km s−1 Mpc−1.

The 19 sources brighter than 8 mJy, with S/N<3.5 and located in the uniform noise regions of the survey maps, account for ≃10 per cent of the submm background observed by COBE-FIRAS (Puget et al. 1996; Fixsen et al. 1998; Hauser et al. 1998). Table 11 shows the inferred star formation rate density for a variety of redshift bands, and for both of our adopted cosmologies. These bright SCUBA sources alone imply a high-redshift star formation rate density (SFRD) in the range ≃0.01–0.07 M yr−1 Mpc−3, comparable to that observed in the optical/UV (Madau et al. 1996; Steidel et al. 1999). If we assume that the whole of the submm background may be attributed to starlight reprocessed by dust, then this would imply a high-redshift SFRD in the range ≃0.1–0.7 M yr−1 Mpc−3, varying only by a factor of 2–3 on the adopted redshift band for a given choice of cosmology. Our results agree very well with Barger et al. (2000), who considered the contribution to the SFRD of sources brighter than 6 mJy, and determined completeness corrected SFRDs ≃0.1–0.4 M yr−1 Mpc−3 in the redshift range 1>z>3, and ≃0.1–0.7 M yr−1 Mpc−3 in the redshift range 3>z>6. In contrast to the marked decline in the SFRD at z<1 originally implied by optical/UV observations (Madau et al. 1996), submm surveys suggest that the SFRD is either steady or gently increasing to perhaps as far back as z=5. Current redshift constraints (Dunlop 2001a, and references therein) suggest that only 10–15 per cent of luminous (<4 mJy) submm sources lie at z>2, and that the median redshift of this population is zmed≃3. A peak in the SFRD around this epoch would not be unexpected given the strong correlation between black hole and spheroid mass found at low redshift (Kormendy & Gebhardt 2001) and the peak in optical emission from powerful quasars at z≃2.5. However, improved redshift constraints are required to establish when (and indeed if) the comoving SFRD reached a maximum.

Table 11.

The comoving star formation rate density of sources brighter than 8 mJy, detected with S/N > 3.50, using the uniform noise regions of the two survey areas and assuming that all of these sources lie within the redshift range given in column 1. Column 2 assumes an Einstein-de Sitter cosmology, and column 3 assumes ΩM = 0.3, ΩΛ = 0.7. H0 = 67km s−1 Mpc−1 was adopted for both cosmologies.

Table 11.

The comoving star formation rate density of sources brighter than 8 mJy, detected with S/N > 3.50, using the uniform noise regions of the two survey areas and assuming that all of these sources lie within the redshift range given in column 1. Column 2 assumes an Einstein-de Sitter cosmology, and column 3 assumes ΩM = 0.3, ΩΛ = 0.7. H0 = 67km s−1 Mpc−1 was adopted for both cosmologies.

Co-moving number density of 8-mJy sources

As explained in the previous section, although we currently lack redshifts for the submm sources uncovered in this survey, existing constraints mean that we can be very confident that virtually all of them lie at z<1, and indeed can be reasonably confident that the majority lie at z<2 (Dunlop 2001a; Fox et al. 2002). It is thus still possible to make a meaningful estimate of the comoving number density of bright (S850<8 mJy) submm sources at high redshift implied by the number counts derived from this survey. The results of this calculation, for a variety of assumed redshift bands, are summarized in Table 12 (again for our two standard assumed alternative cosmologies). From this table it can be seen that the comoving number density of dust-enshrouded starburst galaxies with star formation rates <1000 M yr−1 lies in the range ≃1–10×10−5 Mpc−3 in very good agreement with Lilly et al. (1999), Barger et al. (1999) and Barger et al. (2000). In a similar manner to the comoving SFRD, this result depends only weakly (a factor of 2–3) on the precise choice of assumed redshift band for the sources, for a given choice of cosmology. Adopting the now strongly favoured flat, Λ-dominated cosmology leads us to conclude that the comoving number density of dust-enshrouded starburst galaxies with star formation rates <1000 M yr−1 is ≃1×10−5 Mpc−3.

Table 12.

The comoving number density of sources brighter than 8 mJy, detected with S/N > 3.50, using the uniform noise regions of the two survey areas and assuming that all of these sources lie within the redshift range given in column 1. Column 2 assumes an Einstein-de Sitter cosmology, and column 3 assumes ΩM = 0.3, ΩΛ = 0.7. H0 = 67 km s−1 Mpc−1 was adopted for both cosmologies.

Table 12.

The comoving number density of sources brighter than 8 mJy, detected with S/N > 3.50, using the uniform noise regions of the two survey areas and assuming that all of these sources lie within the redshift range given in column 1. Column 2 assumes an Einstein-de Sitter cosmology, and column 3 assumes ΩM = 0.3, ΩΛ = 0.7. H0 = 67 km s−1 Mpc−1 was adopted for both cosmologies.

This is an interesting number. It is over an order of magnitude smaller than the number density of galaxies brighter than L* in the present-day universe, as inferred from the K-band luminosity function (Cowie et al. 1995; Glazebrook et al. 1995; Gardner et al. 1997; Szokoly et al. 1998; Kochanek et al. 2001), but is an order of magnitude greater than the comoving number density of bright optical quasi-stellar objects (QSOs) (MV>-24) at z≃2–3 (Warren, Hewett & Osmer 1995).

With reference to the present-day K-band luminosity function, a comoving number density of 1×10−5 Mpc−3 corresponds to a galaxy luminosity 1–1.5 mag brighter than L*, or equivalently to galaxies 3–4 times more massive than an evolved L* galaxy. In this regime the present-day galaxy population is completely dominated by massive ellipticals (Kochanek et al. 2001). Interestingly, the hosts of present-day FRII radio sources, and the more luminous radio-quiet quasars are also confined to this same high-mass regime (Dunlop 2001b; Dunlop et al. 2001).

Thus, if one wants to attempt to link the high-redshift population of very luminous (i.e. SFR≃1000 M yr−1) dust-enshrouded starburst galaxies to a low-redshift population, purely on the basis of number-density coincidence, then the simplest connection is that the bright SCUBA galaxies are the progenitors of present-day massive ellipticals with stellar masses ≃1012 M. This is not unreasonable, given that such objects require star formation to be sustained at ≃1000 M yr−1 for ≃1 Gyr to assemble their present-day stellar populations. It is also consistent with the discovery that bright SCUBA detections of (comparably massive) radio galaxies are largely confined to z<2 (Archibald et al. 2001).

If the luminous submm SCUBA sources uncovered by this survey are indeed the progenitors of the present-day massive elliptical population (Dunlop 2001c), and if, as current evidence suggests, the majority of bright submm sources lie at z<2 (Dunlop 2001a), then one might reasonably expect to find a comparably numerous population of passively evolving massive ellipticals at intermediate redshifts. One way to search for such a population is via surveys designed to detect extremely red objects (EROs), such as that recently undertaken by Daddi et al. (2000, 2001). One can obtain a rough estimate of the comoving number density of passively evolving massive ellipticals at intermediate redshifts by considering the surface density of EROs in the Daddi et al. survey with R-K<5.3 (setting a lower redshift boundary of z≃1) and K>18.5 (setting an approximate upper redshift boundary of z≃2 for ellipticals of comparable mass to bright radio galaxies; Jarvis et al. 2001). The surface density of such objects is ≃0.1 per square arcmin (≡350 per square degree), very similar to the surface density found here for 8-mJy SCUBA sources. For the appropriate redshift band, 1>z>2, this surface density converts into a comoving number density of 3×10−5 Mpc−3, assuming ΩM=0.3, ΩΛ=0.7, which (as can be seen from Table 12) agrees with the comoving number density of bright SCUBA sources to within a factor of 2 or 3 (depending on the choice of assumed redshift band for the SCUBA sources).

This numerical coincidence provides further circumstantial evidence for an evolutionary path for all massive ellipticals which mirrors that which has already been largely established for massive radio galaxies (Archibald et al. 2000a; Willott, Rawlings & Blundell 2001; Jimenez et al. 1999), i.e. SCUBA source at z≥2.5→ ERO at z≃1.5→3–4L* evolved elliptical at z=0.

A key test of this picture will be to establish whether or not the bright SCUBA sources display comparable or even stronger spatial clustering than the EROs, for which Daddi et al. (2000) report r0≃11 h−1 Mpc (for objects with R-K<5). As described in the next section, there is tantalyzing evidence of angular clustering in the 8-mJy survey, but it is clear that a substantially larger submm survey (approaching 1° in size) will be required to settle this issue.

Clustering

If the bright 850-μm sources are indeed the progenitors of massive elliptical galaxies then they should be strongly clustered, an inevitable result of gravitational collapse from Gaussian initial density fluctuations since the rare high-mass peaks are strongly biased with respect to the mass. There is a great deal of evidence to support the presence of this bias at high redshift. The correlations of Lyman-break galaxies at z≃3 (Steidel et al. 1999) are almost identical to those of present-day field galaxies, even though the mass must have been much more uniform at early times. Furthermore, the correlations increase with UV luminosity (Giavalisco & Dickinson 2001), reaching scalelengths of r0≃7.5 h−1 Mpc– approximately 1.5 times the present-day value. In the case of luminous proto-ellipticals we expect an even stronger bias, as we select not just massive galaxies but also those that have collapsed particularly early, in order to generate the oldest stellar populations. This is suggested by studies of the local Universe which have shown that early-type galaxies are much more clustered than late-type galaxies (e.g. Guzzo et al. 1997; Willmer, da Costa & Pellegrini 1998), and more recently by the findings of Daddi et al. (2000), who have investigated the clustering properties of extremely red objects (EROs). They detect a strong clustering signal of the EROs which is about an order of magnitude larger than the clustering of K-selected field galaxies, and also report a smooth trend of increasing clustering amplitude with increasing R-K colour, reaching r0≃11 h−1 Mpc for R-K<5. These results are probably the strongest evidence to date that the largest fraction of EROs is composed of ellipticals at z<1.

We already have some hints of the possibility of strong clustering in the bright submm population from the discovery of a strong excess of bright SCUBA sources around high-redshift active galactic nuclei (AGN; Ivison et al. 2000). Motivated by this result, and the potential connection with the intermediate-redshift EROs discussed above, we have attempted to quantify the possible strength of clustering in our two survey fields via calculation of the two-point correlation function. The results, based purely on sources with S/N<3.50, are shown in Figs 12 and 13. There is tentative evidence of clustering on scales of 1–2 arcmin in both of our survey fields, but most particularly, and somewhat stronger, in the ELAIS N2 region. Referring to Fig. 9, we can in fact see by eye that the most significant (S/N<3.50) sources in ELAIS N2 do not conform to a homogeneous distribution across the field – rather there are two apparent concentrations of 850-μm sources in the top left and bottom right of the image, approximately at RA 16h37m00, Dec. +41°05′00, and RA 16h36m30, Dec. +40°56′00 respectively, with an apparent underdensity of sources in the intervening regions. A second peak in w(θ) at ≃600 arcsec (the distance between these overdensities) and a trough at ≃400 arcsec are seen to reflect this ‘by eye’ distribution. Very interestingly, we see the same large-scale inhomogeneities in the Chandra X-ray image of the ELAIS N2 region (Almaini et al., in preparation), although the coincidence of X-ray and SCUBA sources is small (>10 per cent). Almaini et al. provide a full discussion and propose a possible explanation for this effect.

It is clear from Figs 12 and 13 that, as a result of small-number statistics, this first attempt at a direct measure of the clustering of bright submm sources has proved inconclusive. However, this should not be taken as evidence that the SCUBA sources are unclustered. In fact, as illustrated in these figures, it is worth noting that these correlation functions are certainly still consistent with the strong clustering signal detected for EROs by Daddi et al. (2000). Clearly a much larger survey (∼0.5 square degrees) will be required to obtain a meaningful measurement of the strength of clustering in the bright submm population. Near complete redshift information may also be required to quantify the extent to which any clustering signal will be partially erased by projection through a wide range in redshift. For example, suppose that the bright SCUBA population spans a relatively wide range in redshift 2>z>5, whereas, as argued above, the bright EROs are confined to the redshift range 1>z>2. In that case, if the strength of the spatial clustering in the two populations was the same, w(θ) as measured from SCUBA images would be expected to be 1.6–1.9 times smaller (depending on choice of cosmology) than that which has been measured for EROs.

Conclusions

We have used the SCUBA camera on the JCMT to conduct the largest extragalactic submm survey to date, in two regions of sky (ELAIS N2 and Lockman Hole East) covering a total area of 260 arcmin2 to a typical rms noise level of σ850≃2.5 mJy beam−1. We have reduced the resulting dataset using both the standard JCMT surf software and our own IDL-based pipeline which allows us to produce zero-footprint maps and noise images. The uncorrelated noise maps from the latter approach have enabled us to quantify the statistical significance of each peak in our maps leading to properly quantified errors on the flux density of our sources. We find 19 sources with S/N<4.00, 38 sources with S/N<3.50 and 72 sources with S/N<3.00. A series of Monte Carlo simulations have allowed us to assess the completeness, positional uncertainties and effects of confusion in this survey, leading to the most accurate estimate of the submm counts over the flux-density range S850=5–12 mJy. Our best estimate of the cumulative source count at S850<8 mJy is graphic, consistent with a range of models involving strong evolution of the dust-enshrouded starburst population.

Adopting ΩM=0.3, ΩΛ=0.7 and H0=67 km s−1 Mpc−1, the sources uncovered by this survey have inferred star-formation rates in excess of 1000 M yr−1, enshrouded in 108–109 M of dust. Assuming simply that the majority of these sources lie at z<1.5, this result implies that their comoving number density is ≃10−5 Mpc−3, coincident with the number density of massive ellipticals with L<3–4L* in the present-day Universe. This comoving number density is also the same as that inferred for comparably massive passively evolving objects in the redshift band 1>z>2 from recent surveys of EROs. Thus, while the bright submm sources uncovered by this survey contribute only ≃10 per cent of the submm background, they can plausibly account for the formation of all present-day massive ellipticals. A key test of this picture is to determine whether the bright SCUBA population is at least as strongly clustered as are the EROs at intermediate redshifts. As a first attempt at such a test we have compared the distribution of our sources with a mock catalogue of fake sources. Our results have proved inconclusive as a result of the small number of significant sources (S/N<3.50) detected in our survey areas, but there are some initial indications of clustering on scales of 1–2 arcmin, particularly in the ELAIS N2 field.

The process of multiwavelength follow-up of this survey is already well underway. SED-based redshift constraints, along with potential optical–infrared identifications of our most significant sources are presented in Paper II (Fox et al. 2002), while the results of IRAM PdB interferometer imaging of the most significant Lockman Hole East is presented by Lutz et al. (2001). Correlations with deep X-ray and radio images will also be discussed in subsequent papers by Almaini et al. (2002) and Ivison et al. (in preparation).

Acknowledgments

SES acknowledges the support of a PPARC Studentship, while JD acknowledges the enhanced research time provided by the award of a PPARC Senior Fellowship. We are very grateful to the many members of the Joint Astronomy Centre for their continued help and support of this project. We also thank Steve Eales and Amy Barger for communicating their source count values, and Omar Almaini for many useful discussions. The JCMT is operated by the Joint Astronomy Centre on behalf of the UK Particle Physics and Astronomy Research Council, the Canadian National Research Council and the Netherlands Organization for Scientific Research. We also thank our anonymous referee for useful suggestions.

References

Almaini
O.
,
2002
,
MNRAS
 , in press ()
Archibald
E. N.
Dunlop
J. S.
Hughes
D. H.
Rawlings
S.
Eales
S. A.
Ivison
R. J.
,
2000
MNRAS
 ,
323
,
417
DOI:
Archibald
E. N.
Wagg
J. W.
Jenness
T.
,
2000
Archibald
E. N.
Dunlop
J. S.
Jimenez
R.
Friaca
A. C. S.
Mclure
R. J.
Hughes
D. H.
,
2001
,
ApJ
 , submitted ()
Baker
A. J.
Lutz
D.
Genzel
R.
Tacconi
L. J.
Lehnert
M. D.
,
2001
,
A&A
 ,
372
,
L37
Barger
A. J.
Cowie
L. L.
Sanders
D. B.
,
1999
,
ApJ
 ,
518
,
L5
Barger
A. J.
Cowie
L. L.
Richards
E. A.
,
2000
,
AJ
 ,
119
,
2092
Baugh
C. M.
Cole
S.
Frenk
C. S.
,
1996
,
MNRAS
 ,
282
,
L27
Blain
A. W.
Kneib
J.-P.
Ivison
R. J.
Smail
I.
,
1999
,
ApJ
 ,
512
,
L87
Borys
C.
Chapman
S.
Halpern
M.
Scott
D.
,
2002
, submitted ()
MNRAS
 ,
Coulson
I.
,
2000
, http://www.jach.hawaii.edu/JACpublic/JCMT/Facility_ description/Pointing/problems.html
Cowie
L. L.
Hu
E. M.
Songaila
A.
,
1995
,
Nat
 ,
377
,
603
Cowie
L. L.
Songaila
A.
Barger
A. J.
,
1999
,
AJ
 ,
118
,
603
Daddi
E.
Cimatti
A.
Pozzetti
L.
Hoekstra
H.
Röttgering
H. J. A.
Renzini
A.
Zamorani
G.
Mannucci
F.
,
2000
,
A&A
 ,
361
,
535
Daddi
E.
Cimatti
A.
Renzini
A.
,
2001
,
A&A
 , in press ()
Dunlop
J. S.
,
2001
FIRSED 2000 Conf. Proc., New Astron. Rev.
 
Elsevier
, in press ()
Dunlop
J. S.
et al
,
2001
QSO Hosts and their Environments, Granada 2001
 .
IAA
, in press ()
Dunlop
J. S.
,
2001
Deep Sub-millimetre Surveys Conf. Proc.
 
World Scientific
,
Singapore
, in press ()
Dunlop
J. S.
Hughes
D. H.
Rawlings
S.
Eales
S. A.
Ward
M. J.
,
1994
,
Nat
 ,
370
,
347
Dunlop
J. S.
McLure
R. J.
Kukula
M. J.
Baum
S. A.
O'Dea
C. P.
Hughes
D. H.
,
2001
,
MNRAS
 , submitted ()
Dunne
L.
Eales
S. A.
,
2001
,
MNRAS
 ,
327
,
697
Dunne
L.
Eales
S.
Edmunds
M.
Ivison
R.
Alexander
P.
Clements
D. L.
,
2000
,
MNRAS
 ,
315
,
115
Eales
S.
Lilly
S.
Gear
W.
Dunne
L.
Bond
J. R.
Hammer
F.
Le Fèvre
O.
Crampton
D.
,
1999
,
ApJ
 ,
515
,
518
Eales
S.
Lilly
S.
Webb
T.
Dunne
L.
Gear
W.
Clements
D.
Yun
M.
,
2000
,
AJ
 ,
120
,
2244
Elbaz
D.
et al
,
1999
,
A&A
 ,
351
,
L37
Efstathiou
A.
et al
,
2000
,
MNRAS
 ,
319
,
1169
DOI:
Fixsen
D. J.
Dwek
E.
Mather
J. C.
Bennel
C. L.
Shafer
R. A.
,
1998
,
ApJ
 ,
508
,
123
Fox
M. J.
et al
,
2002
,
MNRAS
 ,
331
,
839
(this issue)
Fruchter
A. S.
Hook
R. N.
,
2002
,
PASP
 , in press ()
Gardner
J. P.
Sharples
R. M.
Frenk
C. S.
Carrasco
B. E.
,
1997
,
ApJ
 ,
480
,
L99
Giavalisco
M.
Dickinson
M.
,
2001
,
ApJ
 ,
550
,
177
Glazebrook
K.
Ellis
R.
Colless
M.
Broadhurst
T.
Allington-Smith
J.
Tanvir
N.
,
1995
,
MNRAS
 ,
273
,
157
Guzzo
L.
Strauss
M. A.
Fisher
K. B.
Giovanelli
R.
Haynes
M. P.
,
1997
,
ApJ
 ,
489
,
37
Hauser
M. G.
et al
,
1998
,
AJ
 ,
115
,
1418
Holland
W. S.
et al
,
1999
,
MNRAS
 ,
303
,
659
DOI:
Hughes
D. H.
Dunlop
J. S.
Rawlings
S.
,
1997
,
MNRAS
 ,
289
,
766
Hughes
D. H.
et al
,
1998
,
Nat
 ,
394
,
241
Hughes
D. H.
Gaztañaga
F.
,
2000
, in
Proc. of 33rd ESLAB Symp., Star formation from the small to the large scale
 , ESA SP-445
Ivison
R. J.
,
1995
,
MNRAS
 ,
275
,
L33
Ivison
R. J.
et al
,
1998
,
ApJ
 ,
494
,
211
Ivison
R. J.
Dunlop
J. S.
Smail
I.
Dey
A.
Liu
M. C.
Graham
J. R.
,
2000
,
ApJ
 ,
542
,
27
Jameson
A.
,
2000
,
PhD thesis
 , Univ. Cambridge, Cambridge
Jarvis
M. J.
et al.  
Rawlings
S.
et al.  
Eales
S. A.
et al.  
Blundell
K. M.
et al.  
Willott
C. J.
et al
,
2001
, in
QSO hosts and their environments, Granada 2001
 .
IAA
, in press ()
Jenness
T.
Lightfoot
J. F.
,
1997
,
Jenness
T.
,
2000
, http://www.jach.hawaii.edu/JACdocs/JCMT/tr/001/84/ index.html
Jimenez
R.
Friaca
A. C. S.
Dunlop
J. S.
Terlevich
R. J.
Peacock
J. A.
Nolan
L. A.
,
1999
,
MNRAS
 ,
305
,
L16
DOI:
Jimenez
R.
Padoan
P.
Dunlop
J. S.
Bowen
D. V.
Juvela
M.
Matteucci
F.
,
2000
,
ApJ
 ,
532
,
152
Kauffman
G.
Charlot
S.
,
1998
,
MNRAS
 ,
294
,
705
DOI:
Kawara
K.
et al
,
1998
,
A&A
 ,
336
,
L9
Kochanek
C. S.
,
2001
,
ApJ
 ,
560
,
566
Kormendy
J.
Gebhardt
K.
,
2001
, in
20th Texas Symposium on Relativistic Astrophysics
 .
AIP
,
New York
, in press ()
Lilly
S. J.
Le Fèvre
O.
Hammer
F.
Crampton
D.
,
1996
,
ApJ
 ,
460
,
L1
Lutz
D.
et al
,
2001
,
A&A
 ,
378
,
70
Madau
P.
Ferguson
H. C.
Dickinson
M. E.
Giavalisco
M.
Steidel
C. C.
Fruchter
A.
,
1996
,
MNRAS
 ,
283
,
1388
Magorrian
J.
et al
,
1998
,
AJ
 ,
115
,
2285
Oliver
S.
et al
,
2000
,
MNRAS
 ,
316
,
749
DOI:
Puget
J. L.
Abergel
A.
Bernard
J. P.
Boulanger
F.
Burton
W. B.
Desert
F. X.
Hatmann
D.
,
1996
,
A&A
 ,
308
,
L5
Renzini
A.
Cimatti
A.
,
2000
, in
The Hy-Redshift Universe: Galaxy Formation and Evolution at High Redshift, Conf. Proc., Berkeley 1999
 , in press ()
Rowan-Robinson
M.
,
2000
,
MNRAS
 ,
316
,
885
DOI:
Rowan-Robinson
M.
et al
,
1997
,
MNRAS
 ,
289
,
490
Saunders
W.
Rowan-Robinson
M.
Lawrence
A.
Efstathiou
G.
Kaiser
N.
Ellis
R. S.
Frenk
C. S.
,
1990
,
MNRAS
 ,
242
,
318
Sawicki
M.
,
2000
,
A&AS
 ,
197
,
6504
Scoville
N. Z.
Young
J. S.
,
1983
,
ApJ
 ,
265
,
148
Serjeant
S.
et al
,
2000
,
MNRAS
 ,
316
,
768
DOI:
Serjeant
S.
et al
,
2002
,
MNRAS
 , in press ()
Smail
I.
Ivison
R. J.
Blain
A. W.
,
1997
,
ApJ
 ,
490
,
L5
Smail
I.
Ivison
R. J.
Blain
A. W.
Kneib
J.-P.
,
2001
, submitted ()
MNRAS
 ,
Steidel
C. C.
Adelberger
K. L.
Giavalisco
M.
Dickinson
M.
Pettini
M.
,
1999
,
ApJ
 ,
519
,
1
Szokoly
G. P.
Subbarao
M. U.
Connolly
A. J.
Mobasher
B.
,
1998
,
ApJ
 ,
492
,
452
Thronson
H.
Telesco
C.
,
1986
,
ApJ
 ,
309
,
L79
van der Werf
P. P.
Kraiberg Knudsen
K.
Labbé
I.
Franx
M.
,
2001
, in
Deep Sub-millimetre Surveys
 .
World Scientific
,
Singapore
, in press ()
Wang
B.
,
1991
,
ApJ
 ,
374
,
456
Warren
S. J.
Hewett
P. C.
Osmer
P. S.
,
1995
,
ApJ
 ,
438
,
506
Willmer
C. N. A.
Da Costa
L. N.
Pellegrini
P. S.
,
1998
,
AJ
 ,
115
,
869
Willott
C. J.
Rawlings
S.
Blundell
K. M.
,
2001
,
MNRAS
 ,
324
,
1

Appendix

Appendix A:The Source-Extraction Algorithm

Suppose we consider a normalized beam map B(x, y) as a source template and that at position (i, j) in the unconvolved image the signal is S(i, j) and the noise is N(i, j). If n peaks above a specified flux threshold are located in the Gaussian-convolved image, we may construct a model to the unconvolved image such that beam maps centred on each peak position are simultaneously scaled to give an overall best fit to the entire image. Using a minimized χ2 fit as the maximum likelihood estimator gives  

(A1)
formula
where ak is the best-fitting flux to the kth peak. Minimizing with respect to each ak;  
(A2)
formula
and m=1⋯n.

Rearranging, we obtain the matrix equation  

(A3)
formula
which may be written in the form  
(A4)
formula
where  
(A5)
formula
is an n×n matrix, and  
(A6)
formula
is a vector of length n.

To recover the best-fitting values of ak we invert matrix α such that  

(A7)
formula

The variance associated with the estimate ak is given by  

(A8)
formula

As αkm is independent of S(i,j);  

(A9)
formula

Consequently we find,  

(A10)
formula

The final term in brackets is simply the matrix [α] and so this expression reduces to  

(A11)
formula

The diagonal elements of [α]−1 are the variances of the fitted parameters ak such that the significance of the peak detection is  

(A12)
formula