The Herschel (cid:2) -ATLAS data release 1 – I. Maps, catalogues and number counts

We present the ﬁrst major data release of the largest single key-project in area carried out in open time with the Herschel Space Observatory . The Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS) is a survey of 600 deg 2 in ﬁve photometric bands – 100, 160, 250, 350 and 500 μ m – with the Photoconductor Array Camera and Spectrometer and Spectral and Photometric Imaging Receiver (SPIRE) cameras. In this paper and the companion Paper II, we present the survey of three ﬁelds on the celestial equator, covering a total area of 161.6 deg 2 and previously observed in the Galaxy and Mass Assembly (GAMA) spectroscopic survey. This paper describes the Herschel images and catalogues of the sources detected on the SPIRE 250 μ m images. The 1 σ noise for source detection, including both confusion and instrumental noise, is 7.4, 9.4 and 10.2 mJy at 250, 350 and 500 μ m. Our catalogue includes 120 230 sources in total, with 113 995, 46 209 and 11 011 sources detected at > 4 σ at 250, 350 and 500 μ m. The catalogue contains detections at > 3 σ at 100 and 160 μ m for 4650 and 5685 sources, and the typical noise at these wavelengths is 44 and 49 mJy. We include estimates of the completeness of the survey and of the effects of ﬂux bias and also describe a novel method for determining the true source counts. The H-ATLAS source counts are very similar to the source counts from the deeper HerMES survey at 250 and 350 μ m, with a small difference at 500 μ m. Appendix A provides a quick start in using the released data sets, including instructions and cautions on how to use them.

. However, the exceptional sensitivity of Herschel, aided by the large and negative k-correction at submillimetre wavelengths (Blain & Longair 1993), has meant that a significant fraction of the sources in H-ATLAS actually lie at high redshift (Amblard et al. 2010;Lapi et al. 2011;González-Nuevo et al. 2012;Pearson et al. 2013). The H-ATLAS survey is therefore useful both to astronomers studying the nearby Universe and to those studying the early (z > 1) Universe. The large area of the survey means that there are also potential uses for it in Galactic astronomy (Eales et al. 2010), with one practical example being a search for debris discs around stars (Thompson et al. 2010).
We selected the H-ATLAS fields to avoid bright continuum emission from dust in the Galaxy, as seen on the IRAS 100 μm image, which shows well the cool low surface brightness dust in the interstellar medium. We also chose the fields to provide the maximum amount of data in other wavebands. In the latter respect, a high priority was to choose fields that had already been surveyed in one of the major optical spectroscopic surveys of the nearby Universe. We selected five fields: a large (∼170 deg 2 ) field close to the North Galactic Pole, which was observed in the SDSS; a second large field (∼285 deg 2 ) close to the South Galactic Pole, which was observed in the 2-Degree-Field Galaxy Redshift Survey (2dF; Colless et al. 2001); and three fields on the celestial equator, which are each about 54 deg 2 in size, which were observed in the SDSS, 2dF and GAMA redshift surveys. The existence of the GAMA data is particularly important because this survey is much deeper than the SDSS (more than two magnitudes in the r band), yielding a surface-density of galaxies with redshifts approximately four times greater. The positions of these latter equatorial fields are given in Table 1. The fields are at right ascensions of approximately 9, 12 and 15 h and hereafter we will call these fields GAMA9, GAMA12 and GAMA15.
Apart from the photometric data in the five H-ATLAS fields, there is a large amount of imaging data in many other wavebands. In particular, these fields have been imaged with the Galaxy Evolution Explorer in the FUV and NUV filters (GALEX; Martin et al. 2005) and with the Wide-field Infrared Survey Explorer in the W1, W2, W3, W4 filters (WISE; Wright et al. 2010) and as part of the UK Infrared Deep Sky Survey Large Area Survey (UKIDSS-LAS; Lawrence et al. 2007) and the VISTA Kilo-Degree Infrared Galaxy Survey in the Z, Y, J, H, K S filters (VIKING; Edge et al. 2013) and the VST Kilo-Degree Survey in the u, g, r, i filters (KIDS; de Jong et al. 2013). This is the first of two papers describing the public release of the H-ATLAS data for the entire GAMA fields. We intend to release the data for the other fields over the following few months. The second paper describing the data release for the GAMA fields is Bourne et al. (2016, hereafter Paper II) which describes the properties of the optical counterparts of the H-ATLAS sources. We released part of the data for the GAMA9 field at the end of the Herschel Science Demonstration Phase (hereafter SDP). This data set was described by Pascale et al. (2011), Ibar et al. (2010), Rigby et al. (2011) and Smith et al. (2011). The data we release now supersedes this earlier data set. A comparison between the two data sets is described in Section 9.1. This paper describes the Herschel images of the fields and the catalogues of sources that we have constructed based on the images. We have tried to describe the data set in a way that is easily comprehensible to astronomers not expert in submillimetre data, referring those who are interested in some of the more technical details to other publications. Nevertheless, there are some unavoidable technical issues that do need to be understood if the data are to be used in a reliable way. The most important of these is the importance of the noise in the images produced by the blending together of faint sources -so-called 'source confusion'. One result of this is that there are two kinds of noise on submillimetre images: instrumental noise, which is usually uncorrelated between image pixels, and confusion noise which is highly correlated between pixels. One practical result of this is that the optimum filter for detecting unresolved sources differs from the point spread function (PSF) of the telescope (Section 3.5), which is the optimum filter for maximizing the signal to noise of unresolved sources if the noise is uncorrelated between pixels (Chapin et al. 2011).
A second technical issue is the statistical bias created by the effect of random errors if the distribution of the property being measured is not uniform. Eddington (1913) recognized this effect for the first time, showing that there was a significant bias in the measured number counts of stars even when the errors on the flux densities of the stars have the usual Gaussian distribution. The cause of this bias is that there are more faint stars to be scattered to bright flux densities by errors than bright stars to be scattered to faint flux densities by errors, and thus any sample of stars in a narrow range of flux density will contain more stars whose true flux densities are lower than this range than stars whose true flux densities are brighter than this range. This effect is present in optical catalogues but is much more severe in submillimetre catalogues because of the steepness of the submillimetre source counts (Clements et al. 2010;Oliver et al. 2010;Lapi et al. 2011). This effect creates systematic errors in both the number counts, the phenomenon spotted by Eddington, and also the flux densities of individual sources. In this paper, we will refer to the first phenomenon as 'Eddington bias' and the second phenomenon as 'flux bias'. 2 The second effect is often called 'flux boosting', but we prefer our name because, in principle, the bias could be negative if there were fewer faint sources than bright sources.
We have included in this paper extensive modelling of the effect of Eddington and flux bias, as well as other simulations and modelling that will allow the astronomical community to use the H-ATLAS images and catalogues in a reliable way. Paper II describes the catalogue of the galaxies from the SDSS r-band images that we have identified as the counterparts to the Herschel sources and the multiwavelength data for each of these galaxies. All of the data described in these two papers can be found on http://www.h-atlas.org.
The layout of this paper is as follows. Section 2 contains a description of the Herschel observations. Section 3 describes the reduction and properties of the images made with one of the two cameras on Herschel, the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010). Section 4 describes the data reduction and properties of the images made with Herschel's other camera, the Photoconductor Array Camera and Spectrometer (PACS; Poglitsch et al. 2010). Section 5 describes the detection and photometry of the sources. Section 6 describes the simulations we have carried out to investigate the effects of source confusion and instrumental noise on the catalogues. Section 7 describes our determination of the 250 μm source counts, corrected for the effect of Eddington Bias using the results of the simulations. Section 8 describes our models of the completeness of the survey and of the flux bias of individual sources at 250, 350 and 500 μm. Section 9 describes a comparison of our flux-density measurements with previous measurements at these wavelengths. Section 10 describes the catalogue of sources. The main results from the paper are summarized in Section 11.
For those who wish to use the data as soon as possible without wading through a lengthy technical paper, we have provided a quick start in the form of Appendix A. In this appendix, we describe the basic data sets we have released, with instructions and cautions on how to use them and references back to the full paper. Appendix B provides a list of the asteroids detected in our maps, since these are not included in the catalogues of sources provided in the data release.

THE HERSCHEL O B S E RVAT I O N S
The GAMA fields were observed with Herschel in parallel mode, with the two Herschel instruments, PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010), operating simultaneously. The Herschel survey of each GAMA field consists of four overlapping rectangular regions, strung along the celestial equator. With one exception, 3 each quadrant was observed twice by the telescope, each observation lasting about 9 h. The telescope scanned along great circles on the sky at a constant angular speed of 60 arcsec s −1 , with the scan directions of the two observations being roughly orthogonal (about 85 • apart). By making maps out of two data sets made with orthogonal scan directions, it is possible, with the correct map-making algorithm, to correct for any large-scale artefacts that would result from the gradual changes in the response of the camera, which was a major concern before launch ( Waskett et al. 2007). In operation, SPIRE, although not PACS, proved sufficiently stable that only a very simple map-making algorithm was required and maps made from observations made with a single scan direction were usually free of artefacts from responsivity changes.
The time delay between the two observations made it possible to look for variations in the submillimetre sky by comparing maps made from the two observations of each field. We used the differences of these maps to look for variable sources but only found asteroids (see Appendix B). We also used these maps for the astrometric registration of the data (Section 3.1). We have not included these maps in this data release (maps made from the individual H-ATLAS observations can be found in the Herschel archive).
During the observations, the telescope scans at a constant velocity along a great circle across the field. At the end of each 'scan leg', the telescope decelerates to rest, moves a constant distance orthogonal to the scan leg, and then scans backwards across the field. A single observation is then made up of a large number of scan legs during which the telescope is moving at a constant speed. Data are still being taken during the sections between scan legs when the telescope is accelerating, but this 'turn-around' data is not included in the final maps. In parallel mode, the scan lines were separated by 155 arcsec in order to achieve a good coverage with both PACS and SPIRE. More details can be found in the SPIRE and PACS Observers' Manuals, which are available at http://herschel.esac.esa.int.

THE SPIRE MAPS
H-ATLAS imaged the sky with the SPIRE camera simultaneously in three submillimetre bands centred at 250, 350 and 500 μm. Each band is approximately 30 per cent wide in λ/λ. More technical detail of the camera is given in Pascale et al. (2011) and a full description of the camera is in Griffin et al. (2010).

Map making
Several map-making software packages are capable of working with time-ordered data pre-processed by the standard SPIRE pipeline. Some of these have been developed directly for SPIRE and others were adapted from different instruments. The description and a detailed comparison of the performance of some of the most used map-makers is presented in Xu et al. (2014).
We processed the data with a very similar method to one we used for the data we released at the end of the SDP (Pascale et al. 2011). In this paper we give a basic description of the method, emphasizing any differences from the method we used during SDP. Some additional technical details can be found in Pascale et al. (2011).
We processed the data within the Herschel Interactive Pipeline Environment (HIPE; Ott 2010). We used HIPE version 8. Our fundamental flux calibration is based on the value of Neptune assumed in SPIRE Calibration Tree Version 8. Note that our SDP results (Pascale et al. 2011;Rigby et al. 2011) were based on Version 5 of the SPIRE Calibration Tree. Between Version 5 and Version 8 there was no change in the calibration at 250 and 500 μm and a change of only 1.0067 at 350 μm, in the sense that the flux densities in the current data release are lower by this factor than they would have been if we had assumed the calibration we used in the SDP.
The Herschel Level-1 data consists of fully calibrated time-line data: files of flux density versus time for each bolometer. We produced the Level-1 data using the standard SPIRE data-reduction pipeline with a few exceptions. The first exception was the way we corrected the time-line data for glitches produced by cosmic rays. We used the SIGMAKAPPADEGLITCHER module instead of the default WAVELETDEGLITCHER module as we found the former performed better in masking glitches for our fast-scan parallel-mode observations.
The second exception is the way we corrected for gradual changes in the bolometer signals created by gradual changes in the temperature of the instrument. SPIRE contains thermistors which monitored this temperature change, and we designed our own procedure to use the thermistor results to correct the bolometer signals. We first visually inspected the thermistor signals for sudden jumps, which are spurious and do not represent temperature changes, and we removed these jumps by adding an appropriate constant to the thermistor timeline immediately after each jump. In the standard pipeline, each scan leg is corrected separately for the drift in the bolometer signal caused by the temperature changes, and the temperature information in the turn-around sections (Section 2) is not used at all. However, to make full use of all available information, we used the thermistor data for the entire observation at once, including the data from the turn-around sections. We fitted the following function to the relationship between the temperature measured by a suitable thermistor (T) and the signal for the ith bolometer (S bolom, i ) for the entire 9 h of data: S bolom, i = a × T + c. We then subtracted this relationship from the bolometer signals, effectively removing the effect of the gradual temperature changes. There were some parts of the time-lines where a linear relation clearly did not fit. These almost always occurred 6 h after a cycle of the SPIRE cooling system and became known in the trade as 'cooler burps'. We fitted the bolometer and thermistor data in these regions with a fifth-order polynomial rather than a linear relation, and used this function to correct the data. This polynomial fitting stage is the only difference from the method we used for our SDP data (Pascale et al. 2011).
To remove any residual gradual drift in signal caused by thermal or other effects, we then applied a high-pass filter to the time-line data after first masking bright astronomical sources. The high-pass filter corresponds to a scale on the sky of 4. • 2, and thus our images will not contain structure larger than this scale. This value is the same as the angular size of the image made from a single 9 h observation and was chosen in order to minimize the 1/f noise of the maps (Pascale et al. 2011).
As described above (Section 2), each field was mapped twice, with the scan direction of the telescope in roughly orthogonal directions. After processing the time-line data as above, we made initial maps using the data from the separate observations. We visually inspected the maps to spot any missed jumps in bolometer or thermistor signals, which we then fixed in the time-line data. We could also see linear features on the maps in the scan direction which were caused by the residual effect of bright cosmic rays. We masked these 'glitch tails' in the time-line data.
We also used these initial maps to perform an astrometric calibration of the data, as done by Smith et al. (2011). We first produced initial source catalogues for each map using our source-detection method (Section 5.1). We then produced histograms of the separations in RA and Dec. between the sources and all objects on the SDSS DR7 r-band images (Abazajian et al. 2009) within 50 arcsec of each source. We then fitted these distributions using a Gaussian model for the SPIRE positional errors and allowing for the effects of galaxy clustering in the SDSS data. More details of this procedure are given in Smith et al. (2011). This procedure allowed us to measure the average difference in positions in both RA and Dec. for each data set between the Herschel positions and the SDSS positions with a precision of ∼0.05 arcsec in each direction. The shifts we found range from less than an arcsec to few arcsec, in agreement with the 1σ pointing uncertainty of ∼2 arcsec given for Herschel (Pilbratt et al. 2010). We used these shifts to correct the astrometry for each Herschel observation, so that the effective astrometric calibration of our maps and catalogues should be the same as that of the SDSS.
We corrected the astrometry of the time-lines and reprocessed them, removing the remaining thermistor or bolometer jumps and masking the glitch tails. We then made the final maps using all the time-line data sets for each GAMA field. We made these images using the standard 'naive' map-maker in HIPE. In this map maker, the flux density in each pixel of the final image is calculated by taking the mean of the values for the samples of the time-line data that have positions that fall within this image pixel. We did not use any turn-around data to make the final images or data from regions for which there were data in only one-scan direction. We did not use the default pixel size of the standard pipeline (6, 10 and 14 arcsec for 250, 350 and 500 μm, respectively), but instead 6, 8 and 12 arcsec at 250, 350 and 500 μm, respectively. We made this choice because this corresponds approximately to 1/3 of the full width at half-maximum (FWHM) of the PSF (see Section 3.2) in all bands and provides a good compromise between having a sufficient number of samples in each image pixel to avoid high shot noise and producing a good sampling of the PSF. Note that this choice differs from the one we used for our SDP data release, in which the maps have pixel scales of 5, 10 and 10 arcsec at 250, 350 and 500 μm, respectively (Pascale et al. 2011).
We provide these images in the data release (Appendix A1). Table 1 gives the precise areas covered by each map. Fig. 1 (top) shows the basic 250 μm images of GAMA9 and Fig. 2 shows the coverage map for the same field. The field we observed during the SDP is the second quadrant from the left. The large-scale features visible in the GAMA9 field are emission from dust in the Galaxy ('cirrus emission'). This was the one field where we had to accept some significant cirrus emission as the price of some high-quality multiwavelength data.
In an H-ATLAS SPIRE image made from an observation made with a single scan-direction, it should be possible to reconstruct the submillimetre structure of the sky in the scan direction on all scales up to ∼4. • 2, the angular scale of the high-pass filter. A single observation, however, contains little information about the structure of the sky in a direction orthogonal to the scan direction. Nevertheless, by using the observations from both orthogonal scan directions and a Maximum Likelihood estimator (Pascale et al. 2011), it is possible to produce an accurate image of the sky up to the maximum scale set by the high-pass filter. However, in this data release our maps have been made with the 'naive' map-maker, and as a result the Fourier modes with scales greater than 20 arcmin may be attenuated. Therefore, our images will be adequate for measuring the submillimetre flux densities of even large nearby galaxies but should not be used to investigate structure on larger scales.
All Herschel images contain emission from dust in the Galaxy. Before trying to detect extragalactic sources, it is preferable to try and remove this emission. We used the Nebuliser algorithm, 4 developed by the Cambridge Astronomical Survey Unit, to attempt to remove the emission from Galactic dust from the SPIRE images in all the three wavebands. The algorithm is fairly simple. In each pixel, the median of the intensity values in a square of N × N pixels centred on the pixel is used to estimate the background at that position. The background map is then smoothed using a box-car mean filter with box size chosen to be N/2 × N/2. The user can also set positive and negative signal-to-noise thresholds relative to the filtered map, outside of which the pixels are not used to compute the median. We used lower and upper thresholds of −10σ and +3σ , and carried out three iterations of pixel rejection. We tested different values of N and concluded that with N = 30 (corresponding to 3 arcmin), there was no significant effect on the flux densities of point sources. As discussed in Section 4.1, some flux is lost from sources whose size is comparable to the median box size, but the effect is negligible for sources with radius 1 arcmin. We have provided these background-subtracted images as part of the data release (see Appendix A1). These images do not contain extended emission from the Galaxy and any kind of extragalactic signal with spatial scale larger than ∼3 arcmin. Fig. 1 (bottom) shows the difference between the GAMA9 map at 250 μm and the same map after Nebuliser has been applied; the application of Nebuliser has clearly efficiently removed the extended emission from the Galactic Figure 1. Top: raw image of the whole GAMA9 field at 250 µm, with strong contamination by Galactic dust emission. The field we observed during the SDP is the second rectangular region from the left. Bottom: images of a small part of the GAMA9 field. The left-hand panel shows the raw 250 µm image and the middle panel shows the same part of the 250 µm image after both Nebuliser (Section 3.1) and the matched filter (Section 3.5) have been applied. The right-hand panel shows a pseudo-colour image of this region made by combining the 250, 350 and 500 µm images. On this image, red represents a source that is particularly bright at 500 µm and blue one that is particularly bright at 250 µm. Red sources are generally ones which contain cold dust or at high redshifts (Amblard et al. 2010). Images realized using CUBEHELIX (Green 2011). dust. Nevertheless, very compact Galactic dust regions and cirrus knots are still present in the map, and our final catalogues will contain compact Galactic sources, such as debris discs (Thompson et al. 2010).

The point spread function
To characterize the PSF in our images, we created images of Neptune (Griffin et al. 2013), using the same pixel scale as for the H-ATLAS images. From these images, we made azimuthally averaged and normalized PSF profiles to use in the source extraction (see Section 5.1). Although the actual SPIRE PSF is slightly elongated (see SPIRE Handbook 5 ), our images are made from multiple observations with different scan directions, so the approximation that 5 http://herschel.esac.esa.int/Docs/SPIRE/spire_handbook.pdf the PSF is circularly symmetric should be a good one. Besides it is the approximation made in the calibration of the SPIRE flux densities (Griffin, private communication). We measured the FWHM of the azimuthally averaged PSF to be 17.8, 24.0 and 35.2 arcsec at 250, 350 and 500 μm, respectively. We provide this azimuthally averaged normalized PSF profile as part of the data release (see Appendix A1).

The instrumental noise
The uncertainty maps produced by HIPE, which come as extensions to the standard images, are obtained by calculating the variance of all the time-line samples that contribute to each pixel. However, this method does not produce an accurate estimate of the instrumental noise in each pixel for two reasons. First, the small number of samples in each pixel means the estimate of the variance itself has a Figure 2. Coverage map of the whole GAMA9 field at 250 µm. The maps show the number of samples per pixel. The coverage, and thus the instrumental noise, is quite uniform through the whole map, with the exception of the parts where the rectangular regions overlap. Image realized using CUBEHELIX (Green 2011). significant uncertainty, leading to two pixels with the same number of samples from equally sensitive bolometers having in general different uncertainties even though the instrumental noise should be the same. Secondly, if the pixel coincides with a bright source the variance will be higher because each time-line sample will have been measured at a slightly different position on the source.
We therefore adopted a different technique for estimating the instrumental noise in each pixel. We used the statistics of the difference between the two maps, m 1 and m 2 , of the same region made with the different scan directions. This jackknife map only contains instrumental noise because any real astronomical structure will be removed by subtracting the images (apart from a small number of asteroids -Appendix B). Pascale et al. (2011) have shown that the histogram of pixel intensities from the jackknife map closely follows a Gaussian distribution. If σ 1 is the instrumental noise from one sample, the variance in each pixel of the jackknife map will be var(m 1 − m 2 ) = σ 2 1 (1/c 1 + 1/c 2 ), in which c 1 and c 2 are the numbers of samples in that pixel in each of the two maps. Our estimate of the instrumental noise per sample, σ 1 , is then given by in which i refers to the ith pixel and the sum is carried out over all the N pix pixels in the image. We have estimated σ 1 for each GAMA field and then calculated the average for the three fields as 30.1, 30.9 and 36.1 mJy beam −1 per time-line sample at 250, 350 and 500 μm, respectively. These values are quite similar to the values obtained by Pascale et al. (2011) from the data of the SDP and also broadly consistent with the values measured by Nguyen et al. (2010) from data taken as part of the HerMES survey (Oliver et al. 2012), after allowing for the fact that the bolometer sampling during the H-ATLAS parallel-mode observations (10 Hz) was less frequent than during the HerMES observations (18.2 Hz). The instrumental noise in each pixel of our final images is given by σ 1 / √ c i , where c i is the number of data samples in that pixel.
The instrumental noise in our fields is quite uniform, because most of the survey area has been observed by the telescope for the same exposure time (see Section 2). The exceptions are the overlapping regions between two adjacent quadrants, which had twice the exposure time of the rest of the survey (see Fig. 2). We have provided maps of this instrumental noise as part of the data release (Appendix A1) and also estimated the typical instrumental noise in each image from the following equation: in which the sum is carried out over all the pixels in the image. We give these estimates in Table 2 (top two panels). We also tried an alternative to equation (3) to estimate the average instrumental noise on a map. 6 These values are given in Table 3. The good agreement between the estimates of the instrumental noise in Tables 2 and 3 gives us confidence that our estimates of the instrumental noise are reliable.

The confusion noise
The best way to detect sources on submillimetre images depends on the relative proportions of instrumental noise and noise due to the confusion of faint sources. For example, in an image with instrumental noise and containing only one point source (and therefore 6 In this alternative method, we started with the noise map produced by the jackknife technique, in which the noise in the ith pixel is given by in which c i is the number of time-line pixels contributing to that pixel. We then used this noise map to create a Monte Carlo realization of an image with no real astronomical sources, on the assumption that the noise has a Gaussian distribution and with the standard deviation of the Gaussian function used to generate the intensity values in each pixel given by σ i . We then plotted a histogram of the intensities in all the pixels of the Monte Carlo image, and then fitted a Gaussian function to this distribution, using the standard deviation of the Gaussian as our estimate of the noise on the image. Table 2. Noise from the variance analysis: σ tot Var is calculated using equation (5) and F lim = 200 mJy, σ inst is calculated using equation (3) and σ conf Var is calculated from the first two using equation (4). These values of the confusion noise and total noise are effectively upper limits on the actual confusion noise at the position of a source and on the 1σ uncertainty in the flux density of the source (see Sections 3.4 and 6.4 for more details). σ inst (mJy beam −1 ) σ conf Var (mJy beam −1 ) σ tot Var (mJy beam −1 ) 250 µm 350 µm 500 µm 250 µm 350 µm 500 µm 250 µm 350 µm 500 µm Table 3. Noise analysis from Gaussian fits: σ tot Gauss is calculated by fitting a Gaussian to the negative side of the histograms of pixel intensities (see Fig. 4), σ inst is calculated from Monte Carlo realizations of the noise (see Section 3.3) and σ conf Gauss is calculated from the first two using equation (4). These values of the confusion noise and the total noise are lower limits on the actual confusion noise at the position of a source and on the 1σ uncertainty in the flux density of the source (see Sections 3.4 and 6.4 for more details). no confusion noise), the maximum signal to noise is obtained by convolving the image with the PSF of the telescope (North 1943).
On the other hand, in an image of confused point sources with no instrumental noise, the optimal way to find the sources would be to take the Fourier Transform of the image, divide this by the Fourier Transform of the PSF, and then take the inverse Fourier Transform of the result, a procedure which results in a perfect deconvolution of the original image. Our H-ATLAS images, with their mixture of confusion noise and instrumental noise, are somewhere between these two extremes. Chapin et al. (2011) have shown that for images in this midway regime, it is possible to calculate a convolving function or 'matched filter' that will produce the maximum signal to noise for unresolved sources for any ratio of confusion to instrumental noise. It is therefore important to measure this ratio, and in this section we describe a detailed analysis of the properties of the noise on the H-ATLAS images.
There is not a unique definition of confusion noise and this makes this matter, ironically, quite confusing. One definition of confusion is that it is the root-mean-square value of the fluctuations in an image due to the emission of the astrophysical components of the sky background. In our maps, the extended Galactic emission has been removed (or at least most of it, see Section 3.1). The only source of confusion noise will then be the galaxies we are trying to detect, in particular the ones fainter than the detection limit. We tried two different methods to estimate the confusion noise on our images. In both methods, we first estimated the total noise on the maps and then calculated the confusion noise from in which σ inst is the instrumental noise measured using the jackknife analysis described in Section 3.3.  Table 2. Note the monotonic increase of σ conf Var : at 500 µm this effect is due to gravitationally lensed sources, which are particularly numerous in the GAMA9 field (see main figure), while at shorter wavelengths is due to bright nearby galaxies, which populate the GAMA12 field (see inset figure).
In the first method, we calculated the total noise from the variance of the map: in which F i is the flux density in the ith pixel, N pix is the number of pixels in the image, F map is the mean value of the pixels in the map and F lim is an upper limit to the flux densities that are used in the calculation. We used an upper limit to the flux density because the effect of confusion on a source is only produced by signals fainter than the flux density of this source. In other words, in a pixel with flux density equal to F lim , the confusion noise is calculated from the variance of all pixels with flux density fainter than F lim . If a brighter pixel contributed to the noise, we would measure a brighter flux in that pixel. The solid lines in Fig. 3 shows the relationship between σ conf Var and F lim for the GAMA9 and GAMA12 fields. The figure shows that the confusion value does not asymptotically approach some fixed value as F lim increases. It is easy enough to see why this should be. If the differential source counts are given by dN/dS ∝ S α , a value of α > −3 leads to a monotonic increase of the confusion noise with increasing F lim . This is likely to be the case for the number counts of submillimetre sources, which at bright flux densities are dominated by a mixture of nearby galaxies, with α −2.5, and gravitationally lensed sources (Blain 1998;). In fact, as shown by Fig. 3, this effect is more evident in the GAMA12 and GAMA15 fields (the latter is not shown, but it is very similar to GAMA12), where we find big local bright galaxies. It is noticeable only at 500 μm in the GAMA9 (see inset of Fig. 3), where no particularly big galaxies are present, but where we detected a conspicuous number of gravitational lensed objects (Negrello et al. 2016). The maximum value of the confusion noise is set by the brightest source in the survey, because in that case we would calculate the variance of all pixels in the map. This limit will tend to be higher for surveys covering larger areas of sky.
We therefore have to be careful when comparing our estimates of the confusion noise with those from other surveys. In the HerMES survey (Oliver et al. 2012), which was deeper but covered a smaller area of sky, Nguyen et al. (2010) used a method very similar to ours to measure the confusion noise, finding values for σ conf Var of 5.8, 6.3 and 6.8 mJy beam −1 at 250, 350 and 500 μm, respectively. Fig. 5 in Nguyen et al. (2010) shows that the confusion values measured in the HerMES fields are changing very slowly with F lim at F lim = 200 mJy, and so to make a fair comparison we calculated σ tot Var and σ conf Var with the same limit. The upper panel of Table 2 gives our measurements for all three GAMA fields. The range of values for the three GAMA fields is ∼0.5 mJy at all wavelengths, and the mean values are slightly higher than those from HerMES. However, the second panel in Table 2 shows our measurements from the images that have had the Galactic dust emission removed with Nebuliser (the dashed lines in Fig. 3 shows how σ conf Var depends on F lim for these images). The values for the confusion noise in Table 2 are now very similar for the three GAMA fields, with a range of 0.2 mJy, with mean values of 6.4, 6.9 and 6.5 mJy beam −1 , fairly similar to the HerMES values. This comparison shows that the emission from Galactic dust makes a significant contribution to the confusion noise, and so in comparing confusion estimates one must take into account both the areas covered by the different surveys and also the emission from Galactic dust.
In the second method, we estimated the noise by fitting a Gaussian to the negative side of the histograms of pixel intensities. Fig. 4 shows that the histograms roughly follow a Gaussian at negative and small positive flux densities, but deviate strongly from a Gaussian at bright flux densities. The natural interpretation of this distribution is that the non-Gaussian tail is produced by individual bright sources, whereas the Gaussian distribution is the result of the combination of the instrumental noise and the population of very faint sources, which have such a high surface-density that the confusion noise they produce follows a Gaussian distribution, as a result of the central limit theorem. Table 3 gives the results of this analysis. As expected, the values for the total noise and the confusion noise given in this table are lower than the values in Table 2 because we are now not including the contribution of the strongly non-Gaussian tail that is seen at positive flux densities (Fig. 4).
From a practical point of view, the crucial thing we need to know is the uncertainty in the flux density of each source in the survey. If a source is close to the detection limit of the survey (∼30-40 mJy), the noise measurements in Table 2 will be larger than the true error in the source's flux density, because these measurements contain a confusion contribution from pixels brighter than this limit. On the other hand, the values in Table 3 will be too low because they do not include a contribution from the positive tail of bright sources. We return to this question in Section 6.4, in which we use the results of a simulation to estimate the true errors in the flux densities. Chapin et al. (2011) show how to calculate a 'matched-filter' that maximizes the signal to noise for unresolved sources from the ratio of the confusion noise and instrumental noise on the map. We have used this method to calculate matched-filters for the SPIRE images, using the values for the instrumental and confusion noise estimated from the Gaussian fits to the histograms of the . Histograms of flux density for the pixels in, from top to bottom, the GAMA12 250, 350 and 500 µm images. The same region, in which the noise is approximately uniform, has been used in all the images. The histograms are fitted well by a Gaussian with an excess at positive flux densities due to bright sources on the images. The continuous line is the Gaussian that fits the negative side of the histogram and represents the total noise in the map. The dotted line shows the distribution for instrumental noise predicted from the jackknife results (Section 3.3). The difference between the total and instrumental noise is the result of the noise caused by the confusion of faint sources. Our estimates of instrumental noise from the jackknife analysis and of the total noise and confusion noise from these fits are given in Table 3. background-subtracted maps (top panel of Table 3). This is because we want to maximize the detection sensitivity of the survey, so we optimize the filter using the noise appropriate for sources near the detection limit. Fig. 5 shows the PSFs and the matched-filters. The matched-filters have the negative lobes that are generated by this method (see fig. A1 of Chapin et al. 2011), which is the feature which makes it possible to reduce the effect of confusing sources by resolving sources that would otherwise be blended together.

Filtering the maps
To investigate whether this method gave an improvement in the signal to noise over using the PSF, we produced two sets of convolved H-ATLAS images, one set convolved with the matched-filter and one set with the PSF. We then applied the same two methods of measuring the noise statistics that we used for the raw images (see Section 3.4). Table 2 gives the results from the variance method and Table 3 gives the results from the Gaussian-fitting method.
First, consider Table 2, which lists the values from the variance analysis. The reduction of the total noise due to the convolution with the PSF results in an increase in signal to noise of 1.09, 1.04 and 1.12 at 250, 350 and 500 μm, respectively. The gain in the signal to noise is only very small because of the large component of the variance that comes from bright sources that are not confused, which does not change when the images are convolved. The reduction of the total noise due to the convolution with the matched-filter results in an increase in the signal to noise of 1.33, 1.30 and 1.39 at 250, 350 and 500 μm, respectively. The gain in signal to noise from using the matched-filter rather than the PSF is therefore 1.22, 1.25 and 1.24 at the three wavelengths. Now consider Table 3, which lists the values from the Gaussian fitting. The reduction of the total noise due to the convolution with the PSF results in an increase in signal to noise of 1.87, 1.62 and 1.52 at 250, 350 and 500 μm, respectively. The reduction of the total noise due to the convolution with the matched-filter results in an increase in the signal to noise of 2.01, 1.87 and 1.79 at 250, 350 and 500 μm, respectively. The gain in signal to noise from using the matched-filter rather than the PSF is therefore 1.07, 1.15 and 1.18 at the three wavelengths.
Therefore, whichever way we assess the noise, the use of the matched-filter rather than the PSF when convolving the images results in a small but useful gain in signal to noise. Chapin et al. (2011) found a similar useful but modest improvement in signal to noise when they applied their technique to data from the Balloonborne Large Aperture Submillimeter Telescope (see their table 1).

Making the images
We observed the sky simultaneously in two bands with the PACS camera, at 100 and 160 μm. The two bands were designed to cover the wavelength ranges 85-130 μm and 130-210 μm (Poglitsch et al. 2010), although in practice the bands are slightly different than this. The detailed spectral response of PACS in the two bands is shown in the PACS Observer's Manual. 7 In parallel mode, the SPIRE and PACS observations were made simultaneously, with the two cameras offset on the sky by 21 arcmin. The final SPIRE and PACS images of the GAMA fields therefore have a small offset.
Processing the PACS data, which was taken in fast-scan parallel mode, presented a number of challenges: (1) the huge data sets (the data for each GAMA field consist of eight or nine observations of 8-9 h each) limited the data processing to relatively powerful machines with at least 64 GB of random access memory; (2) the noise power on PACS images has a weak dependence on spatial frequency (∝ 1/f α with α 0.5), which makes it difficult to reduce the noise by spatial filtering without affecting the properties of extended sources; (3) the length of the observations meant that the observations were affected by rare strong 'glitches' that are missed by the standard data-reduction software; (4) the on-board averaging of the PACS data that was required in parallel mode to make it possible to transmit the data to the Earth made it difficult to disentangle glitches from real sources; (5) the on-board averaging significantly affects the PSF, making it more elongated along the scan direction (especially at 100 μm).
The standard data-reduction pipeline provided by the PACS Instrument Control Centre (ICC) working within the Herschel Interactive Processing Environment (HIPE v10.2747; Ott 2010) was not optimized for the H-ATLAS data and so we made some modifica-7 http://herschel.esac.esa.int/Docs/PACS/html/pacs_om.html tions to it. The data reduction procedure we adopted is largely the one described by Ibar et al. (2010), but we made two major changes to the procedure described in that paper.
We processed the Level 1 data (the 'timelines' -signal versus time for each bolometer) largely using the method described in Ibar et al. (2010) with one major modification, which is the way we dealt with spurious jumps in timelines ('glitches'). Our very long observations meant that our data included rare types of glitch. In particular, we found that there were several types of strong glitch that occurred every 1-2 h that were missed by the de-glitching methods used in the standard pipeline (Fig. 6). We wrote a program specifically designed for the H-ATLAS data to find and correct them. We corrected these glitches by fitting them with either a Heaviside or an exponential function and subtracting this function from the timeline. If none of these functions provided a good fit to the glitch, we simply masked the part of the timeline containing it. Fig. 6 shows the effect on the final images of detecting and removing these glitches.
The second major modification we made to the old procedure was how we dealt with the problem that the noise on PACS images only depends very weakly on spatial frequency (see above). In producing our SDP images, Ibar et al. (2010) tried to reduce the noise on the images by first filtering of the time-line data with a high-pass filter. We then made 'naive' images from the timeline data, in which the sky brightness in each image pixel is the average of the values for the time-line samples that fall within this pixel. The filtering did reduce the noise, but at the price of removing extended structure and also reducing the flux densities of the point sources, an effect we attempted to correct for (Ibar et al. 2010).
The naive images do not use the information that there are at least two observations of each field, with roughly orthogonal scan directions. A number of algorithms, such as the Scanamorphos algorithm (Roussel 2013), use this information to separate genuine extended emission on the sky from noise fluctuations in the time-line data with low spatial frequency.
Our approach for the Data Release 1 was therefore not to apply high-pass filtering to the Level-1 data but instead to use one of these algorithms to distinguish the noise fluctuations with low spatial frequency from genuine large-scale emission.
As the first step, to calibrate the astrometry, we made naive images from the data for each individual 9 h observation. In principle, we should have been able to apply the astrometric corrections for SPIRE maps (see Section 3.1) to the PACS data. However, because the new PACS images were created using a different model of where Herschel was pointing (the 'pointing product' provided within HIPE v10.2747) to that used for making the SPIRE images, we could not simply assume that the sources detected by PACS were at exactly the same positions as those detected by SPIRE.
To measure any offsets between the SPIRE and PACS positions, we selected sources detected on the 160 μm images with peak flux densities greater than 65 mJy pixel −1 and measured their positions by fitting each source with a 2D-Gaussian function. All of these sources were detected on the SPIRE 250 μm images and, for each individual PACS image, we calculated the mean differences between the positions of the sources measured from the PACS and SPIRE images. These mean differences were usually 1 arcsec in both RA and Dec. We made our final PACS images after applying an astrometric correction to each individual PACS data set, which was the sum of the astrometric correction obtained from making the comparison between SPIRE and SDSS (see Section 3.1) and the offset obtained from comparing the PACS and SPIRE positions.
We tested three different algorithms for making the final PACS images: the standard implementation of the Scanamorphos algorithm (Roussel 2013), the JScanamorphos version of the algorithm provided as part of HIPE (Graciá-Carpio, Wetzstein & Roussel 2015), and Unimap (Piazzo et al. 2015). We tested the images made using the different algorithms in two different ways. First, we measured the noise on the images by measuring flux densities through a large number of apertures placed at random positions on the images (Section 5.3). Secondly, we measured the PSF for each set of images using a stacking technique that is described in the next section. We found that the noise on the Scanamorphos images was lower than on the other two sets of images, but that the PSF on these images had negative sidelobes extending ∼1 arcmin from the peak of the PSF in each scan direction. We are not sure of the reason for the sidelobes, but their existence suggests that some spatial filtering was occurring in the version of the algorithm we tested, which would also explain the lower noise.
There was no strong reason from our tests to prefer one of the other two algorithms but for practical reasons we decided on JScanamorphos.
We made the final images from the Level-1 data, after making the astrometric corrections, using JScanamorphos. We made the images using the 'parallel' and 'non-thermal' flags. Since the PACS data set for each GAMA field is too large for it to be practical to run JScanamorphos on the entire data set, we made images separately for each of the four rectangular regions (see Fig. 1 to visualize this). We used a pixel scale of 3 arcsec at 100 μm and 4 arcsec at 160 μm, which is approximately 1/3 of the FWHM at each wavelength, which is how we selected the pixel scale for the SPIRE maps (Section 3.1).
We then removed any residual large-scale artefacts from the images by applying the Nebuliser algorithm, which removes largescale emission by using a combination of mean and median filters.
There is obviously a danger in this of removing genuine emission from extended sources, such as nearby galaxies. We tested this effect by injecting artificial galaxies on to the images, applying Nebuliser and then measuring the flux densities of the galaxies. The galaxies were given exponential discs truncated at five scalelengths, and we tested a range of diameters from 12 to 96 arcsec and fluxes from 20 mJy to 1 Jy. This showed that a significant fraction of flux is lost from faint extended sources when the source radius is >1/3 of the filter scale. We chose to use a median filter scale of 300 arcsec, and mean filter scale 150 arcsec, which is means that galaxies smaller than 100 arcsec in radius will not be affected by the filter.
Finally, we masked parts of each image for which there was data from only one-scan direction, and then combined the images for each GAMA field using the SWARP algorithm 8 (Bertin et al. 2002).
The layout of the PACS mosaics is the same as that of the SPIRE mosaics shown in Fig. 2: four overlapping rectangular regions. Each quadrant has been made from two 9 h observations with orthogonal scan directions. In the overlapping regions the images have been made from four observations, and so the noise level in these regions is approximately √ 2 less than in the other areas. In contrast to the SPIRE images, the PACS images are entirely dominated by instrumental noise, which is much greater than the confusion noise expected at these wavelengths (Berta et al. 2011).

The point spread function
The PACS team has carried out a detailed study of the PACS PSF (Lutz 2015). Although the PACS team did not have available observations of a point source in our observing mode -parallel mode with a scan speed of 60 arcsec s −1 -it was possible to estimate the PSF for this observing mode in the following way. The crucial feature of the PSF in parallel mode is that it is elongated in the scan direction, especially at 70 and 100 μm, because of the on-board averaging of the PACS data necessary to transmit both the PACS and SPIRE data to the Earth. Using observations with PACS alone of the asteroid Vesta and Mars made with a scan speed of 60 arcsec s −1 and with the same scan angles used in parallel mode, the PACS team simulated in software the on-board averaging used in parallel mode to produce an estimate of the PSF out to a radius of 1000 arcsec. However, although this should provide a reasonable estimate of our PSF, the real PSF also depends on the pixel size of the map, the spectral energy distribution (SED) of the source, and the algorithm used to make the image (Lutz 2015). For these reasons, we attempted to estimate the real PSF for our observations from the observations themselves.
Since, not surprisingly, there is no very bright point source on our images, we used a statistical 'stacking analysis' to estimate the PSF. To do this, we started with the optical counterparts to the H-ATLAS sources presented by Paper II. We used only the optical counterparts with spectroscopic or photometric redshifts <1, since it is unlikely a galaxy at higher redshift will have any discernible emission in the PACS bands, and with optical sizes (SDSS parameter ISOA) less than 5 arcsec, in order to ensure that the counterpart is likely to be unresolved in the PACS wavebands. There are ∼10 4 of these counterparts that satisfy these criteria in each field. We then extracted a 80 × 80 pixel image centred on each counterpart from the overall PACS image, and then added these together, which should provide a good estimate of the real PSF of our observations. . The dashed black line shows the PACS team's estimate of the EEF for our observing mode (Lutz 2015). The coloured lines show our estimates for the three GAMA fields from the stacking procedure described in the text (see Section 4.2). These curves have been normalized so that they have the same value as the PACS team's EEF at a radius of 30 arcsec.
Our estimates of the beams show no obvious artefacts from our processing method, in particular from the Jscanamorphos imaging algorithm. The empirical PSFs we derived are not recommended to be used as a filter to smooth the map, because of the uncertainties in the peak. We therefore do not release them. For filtering the map, we recommend instead the use of our Gaussian fit. We fitted a two-dimensional Gaussian to the empirical PSF, finding that the FWHM is 11.8 × 11.0 arcsec at 100 μm and 14.6 × 12.9 arcsec at 160 μm. These are slightly larger than the measurements from the PACS team's simulated PSF: 10.7 × 9.7 arcsec at 100 μm and 13.7 × 11.5 arcsec at 160 μm (Table 6; Lutz 2015). We also fitted an azimuthally symmetric Gaussian to our empirical PSFs, which is the one we recommend for smoothing the maps, obtaining a value for the FWHM of 11.4 and 13.7 arcsec at 100 and 160 μm, respectively. It is these values we used when carrying out photometry on the PACS images (Section 5.3). Fig. 7 shows the encircled energy fraction (EEF) plotted against radius for the PSF derived from each of the GAMA fields. We have also plotted in the figure the EEF derived from the PACS team's simulated PSF (Lutz 2015). We have normalized all the curves to the value of this EEF at a radius of 30 arcsec. The PACS team's EEF increases slightly faster with radius than the estimates from the three GAMA fields, as expected given the difference between the values of the FWHM.
The EEF is the crucial function for correcting aperture photometry (Section 5.3) to a standard radius. We cannot derive an empirical EEF from our images much beyond a radius of 30 arcsec because the signal to noise produced by our stacking method is too low. We therefore constructed an EEF that can be used for aperture photometry by using the PACS team's EEF for radii between 30 and 1000 arcsec and the average value of our empirical EEFs for the three GAMA fields at smaller radii, since these should best reflect our actual observations, normalizing our empirical curve so that it has the same value as the PACS team's EEF at a radius of 30 arcsec. We supply these EEFs in the data release.

The unresolved SPIRE sources
Almost all the sources detected on the H-ATLAS images are unresolved by the telescope's PSF. We have developed our own program to find unresolved sources on Herschel images: the Multi-band Algorithm for source Detection and eXtraction (MADX). 9 This will be described in detail elsewhere (Maddox, in preparation). Here, we describe the basic algorithm.
The first step in the MADX source extraction is to use Nebuliser to remove the diffuse Galactic dust emission from all the maps in the three bands (see Section 3.1).
In the second step, the images are convolved with the matchedfilter (Section 3.5, Fig. 5). In carrying out the convolution, the contribution of each pixel of the input image is weighted by the inverse of the square of the instrumental noise in that pixel (Section 3.3). The matched filter is constructed using the values of the instrumental and confusion noise given in the top one of the three panels in Table 3.
A map of the variance in this convolved map is created in the following way. First, the map of the variance of the instrumental noise (Section 3.3) is convolved with the matched-filter. Then the square of the confusion noise is added to the convolved map. For this confusion term we used the mean values of the standard deviation from the middle panel of Table 3, which are 2.0, 3.1 and 3.9 mJy at 250, 350 and 500 μm, respectively. This variance map is used by MADX to determine the signal to noise with which a source is detected. The mean values of the square root of the variance given by this map are 5.2, 5.6, 6.7 mJy beam −1 . Note the important point that these values are less than the more conservative values we have adopted for the errors in the flux densities of the sources (see Section 6.4). We supply both the convolved images used by MADX and the noise maps of these convolved images as part of this data release (Appendix A1).
In the third step, the maps at 350 and 500 μm are interpolated on to images with the same pixel scale as the 250 μm map, and the three images and their associated variance maps are then combined together to form a single signal to noise or 'detection' image. In practice, we chose to give zero weighting to the images at 350 and 500 μm, so our detection image was simply the 250 μm image.
The reason for this choice was that this image contains many more sources than the other SPIRE images, and the smaller size of the PSF leads to more accurate positions for the sources. MADX produces a list of potential sources by finding all peaks in the detection image with signal to noise >2.5.
The disadvantage of this approach is that we might miss sources that are very faint at 250 μm, but bright at the two longer wavelengths. By searching for 2.5σ detections at 250 μm, but only releasing a catalogue of sources detected at a higher significance in at least one of the three wavelengths, the catalogues at the two longer wavelengths should be fairly complete. Rigby et al. (2011) investigated the alternative 'red prior' method in which the three images are combined in a way that makes the detection image sensitive to sources that are bright at the two long wavelengths. By comparing catalogues made in the two ways, they concluded that in our 250 μm-only approach, of the sources that should appear in a 5σ catalogue, 7 per cent of the sources at 350 μm and 12 per cent of the sources at 500 μm will have been missed because they were too faint to be detected at 250 μm. We have already developed an extraction method that allows for the detection of 250 μm dropouts. The method is still under testing, but it has been used to select samples of ultrared sources (Ivison et al. 2016).
After finding the peaks detected at >2.5σ , the next step in the algorithm is to sort the peaks in order of decreasing significance. The program fits a Gaussian to each peak of the detection map to provide an estimate of the source position accurate to a fraction of a pixel. The position of the source can be significantly affected by the presence of nearby sources (Rigby et al. 2011), but the negative sidelobes of the matched-filter (see Fig. 5) help to reduce this problem.
In the next step, MADX makes a first approximation of the flux density of sources in all three wavebands simply using the pixel value nearest to the source position in each of the filtered maps. Then for each band the sources are sorted in order of decreasing brightness. The flux density of the brightest source is estimated at the precise subpixel position determined from the detection image, using a bicubic interpolation between the flux densities in the surrounding (3 × 3) pixels. Using this measured flux density and the matched-filter-convolved PSF, this source is subtracted from the map. The program then moves to the next brightest source and goes through the same set of steps. Note the sequence of image subtraction goes from brightest to faintest in each band, so the source subtraction can happen in a different order in each band. The subtraction of sources is important because this reduces the errors in flux densities of the faint sources caused by the wings of brighter sources.
The final step is to correct the flux densities for a small bias arising from the fact that sources are generally not found at the centre of a pixel, and the bicubic interpolation does not quite recover the true peak flux density. We multiplied the flux densities by 1.009, 1.013 and 1.010 at 250, 350 and 500 μm, respectively, factors we estimated from the simulations described in Section 6. For a similar reason, we need a further correction to the 350 and 500 μm flux densities: we noticed that, when sources are fainter than 40 mJy at 250 μm, their flux densities in the two longer waveband are slightly underestimated. We believe this is because the 350 and 500 μm flux densities are measured at the 250 μm positions, and when the 250 μm flux densities are very low, the positions will be inaccurate. These corrections depend on the 250 μm flux density only and are applied exclusively to sources fainter than 40 mJy at 250 μm. We divided the flux densities at 350 and 500 μm by C 350 and C 500 , respectively, factors we estimated from the simulations described in Section 6 and equals to These estimates of C 350 and C 500 are consistent with the analytical derivation of the same corrections calculated assuming a Gaussian beams and the signal-to-noise ratio (SNR) expected for sources fainter than 40 mJy at 250 μm.
The final flux densities are monochromatic flux densities at 250, 350 and 500 μm, calculated on the assumption that F ν ∝ ν −1 , in which ν is frequency.
On top of the basic errors in the flux densities, which we discuss in Section 6.4, there is an additional error due to the uncertain photometric calibration of Herschel. At the time of writing (the calibration of Herschel is still improving), the error in the flux densities in the SPIRE bands arising from the absolute uncertainty in the flux density of Neptune is 4 per cent in all three bands and there is an additional error of 1.5 per cent that is uncorrelated between the bands (SPIRE Handbook). As recommended in the SPIRE Handbook, we add these errors linearly and use a value of 5.5 per cent as the total calibration error.
A few of the sources detected by MADX are asteroids. We removed these from the catalogues, although the positions on the images are given in Appendix B.

The extended SPIRE sources
The flux density measured by MADX will be an underestimate of the true flux density if the source is resolved by the PSF of the telescope. It is often hard to be sure from the SPIRE images alone whether a source is genuinely extended or whether it is multiple unresolved sources confused together. Apart from a handful of Galactic sources, a source is only likely to be extended if it is associated with a nearby galaxy. We have therefore used our catalogue of optical counterparts to the Herschel sources, which is described in Paper II, to decide whether a source is likely to be extended.
Paper II investigated whether a potential optical counterpart on the SDSS r-band image is genuinely associated with a Herschel source using a Likelihood method, based on the positional error of the Herschel source, the angular distance between the potential counterpart and the Herschel source and the r-band magnitude of the counterpart. Using this method, they estimated the probability (the reliability -R) that the potential counterpart is genuinely associated with the Herschel source. They treated any object with R ≥ 0.8 as being the likely counterpart to the Herschel source. 10 One of the quantities available for each galaxy from the SDSS catalogues is the optical size (ISOA, in SDSS nomenclature) which is defined as the length of the semimajor axis out to an isophote corresponding to 25 mag arcsec −2 . We only used values of ISOA for galaxies with SDSS r-band magnitude (the SDSS quantity rmodelmag) < 19, because values of ISOA for galaxies fainter than this limit are unreliable. In deciding whether a SPIRE source is likely to be extended, and thus whether aperture photometry is necessary, our basic assumption is that the radius of the submillimetre emission is proportional to ISOA.
We measured flux densities using aperture photometry for all reliable optical counterparts (R > 0.8) of the sources detected by our MADX source-finding algorithm (Section 5.1) and with ISOA > 10 arcsec, centring the aperture at the position given in the SDSS for the galaxy. We converted the images from Jy beam −1 to Jy pixel −1 by dividing the values in each pixel by a factor C conv , which is given by the area of the beam divided by the area of a pixel. We used the current estimates of beam from the SPIRE handbook, which are 469, 831 and 1804 arcsec 2 at 250, 350 and 500 μm, respectively.
We measured the flux densities from the images from which the background had been removed using Nebulizer (Section 3.1) and after setting the mean of the maps to zero. The alternative would have been to measure the flux densities from the raw maps, using a annulus around each galaxy to estimate the level of the background. However, in tests we found that the two methods gave very similar results for even the biggest galaxies, whereas the errors in the flux densities when starting from the raw maps were larger.
We used a circular aperture with a radius, r ap , given by In this equation, ISOA is in units of arcsec (we converted the values in the SDSS data base, which are in units of SDSS pixels, using the SDSS pixel size of 0.3958 arcsec). To increase photometric accuracy when the aperture size is small relative to the pixel size, we divided each pixel into 16 subpixels, assigning one sixteenth of the flux density in the real pixel to each subpixel. We corrected all the aperture flux densities for the fraction of the emission that falls outside the aperture because of the extended profile of the PSF (Griffin et al. 2013). We used a version of the PSF appropriate for a source with an SED with monochromatic flux density, F ν ∝ ν −1 , which is the SED assumed in the pipeline (SPIRE Handbook). We provide the table with the SPIRE aperture corrections as part of this data release. In Appendix A3, we describe how corrections should be made to the flux densities in the catalogues to obtain flux densities for more realistic SEDs. We estimated the uncertainties by placing at random positions around each sources 2000 apertures with the same size as the one used to measure the flux of the source and at distance between 2 and 15 times the aperture radius.
In general, the flux error scales as r 3/2 ap for small apertures and is flatter for aperture radii larger than ∼45 arcsec. This relationship is different from the one expected for a Gaussian noise (∝ r ap ) because of the confusion component. The flattening at larger scales is also expected, because Nebuliser removes the large-scale power (Section 3.1). No flattening is indeed visible if the same process is applied to the raw maps.
This Monte Carlo technique takes also into account fluctuations in the background and variation of the confusion noise within the map, because it is applied to each source independently. Fig. 8 shows the ratio of the flux density measured from aperture photometry to the flux density measured by MADX (Section 5.1), in which it is assumed that the sources are unresolved, plotted against ISOA. The red points show the mean ratio of flux densities in bins of ISOA, with the error bars showing the errors on the means. At ISOA >10 arcsec, flux densities are systematically greater than the flux densities measured by the source-detection algorithm. As the result of this analysis, in first approximation we would expect the aperture flux density to be the best estimate of the flux for galaxies with ISOA >10 arcsec. However, we prefer this estimate of the flux instead of the flux densities measured by the source-detection algorithm, only when the result of the aperture photometry is significantly larger than the flux measured by MADX. This is true when where F ap and σ ap are, respectively, the flux density estimated by the aperture photometry and the uncertainty from the Monte Carlo process and F ps and σ tot Var are the flux density measured by the source-detection algorithm and the uncertainty estimated for unresolved sources (see Section 5.1 and equation 12). When this condition is not satisfied, the fluxes estimated by MADX are preferable, since the algorithm is optimized for unresolved sources and the fluxdensity errors will thus be smaller than in the aperture photometry. For galaxies with large optical sizes, we visually compared the aperture given by equation (7) with the 250 μm emission of the galaxy. In ∼100 cases, the aperture was not well matched to the 250 μm emission, either being too small, too large, with the wrong shape or including the flux from a neighbouring galaxy (see Fig. 9 for examples). In these cases, we chose a more appropriate aperture for the galaxy. The details about the customized apertures of these galaxies are given as part of the data release.

PACS photometry
As a result of the lower sensitivity of the PACS images, all the sources detected on the PACS images are detected on the SPIRE 250 μm image. Therefore, we have measured flux densities in the two PACS bands for the sources detected with MADX on the SPIRE 250 μm image. We have measured these flux densities using aperture photometry, centring the aperture either on the position of the SPIRE source or, when possible, on the position of the optical counterpart to that source.
To determine the aperture size that would maximize the signal to noise of the photometry, we followed a similar 'stacking procedure' to that described in Section 4.2. We started with all the optical counterparts to the SPIRE sources which have z < 1 (spectroscopic or photometric) and optical size (SDSS parameter ISOA) Figure 9. Two examples of sources that, after visual inspection, required a customized aperture defined manually (green line), because the automatic aperture defined by equation (7) (blue line) was inappropriate. Top: automatic aperture having a wrong shape. Bottom: automatic aperture including flux from neighbouring sources. Images realized using CUBEHELIX (Green 2011). <5 arcsec. We measured the flux density for each galaxy from the PACS images through 20 apertures with radii equally spaced from 5 to 25 arcsec (we did not need to measure a sky value because of the application of Nebuliser to the images 11 ). In our signal-to-noise calculation, the signal for an aperture of a given radius is the sum of the measurements for all the galaxies. 11 We checked whether the application of Nebuliser led to any loss of flux for big galaxies by also carrying out aperture photometry on images to which Nebuliser had not been applied. In this photometry, we estimated the background from an annulus centred on the galaxy with radii equal to 1.5 and 2.5 times the radius of the aperture. We found that the application of Nebuliser leads to a significant loss of flux for galaxies if their diameter is larger than about 1/3 of the filter scale, confirming the results of our Monte Carlo simulation (see Section 4.1).
To estimate the noise, we used a Monte Carlo simulation in which, for each aperture size, we placed 3000 apertures randomly on the images, calculating the standard deviation of the flux densities measured in all the apertures. In the Monte Carlo simulation, we used aperture sizes radius of 5,10,15. . . , arcsec, up to a maximum radius of 100 arcsec. From this simulation, the error for aperture photometry, σ ap , on the PACS images is given by the following relationship: σ ap = C(r/10 arcsec) β 2/N scan (9) in which r is the radius of the aperture, the constant C is 0.023 and 0.020 Jy at 100 and 160 μm, respectively; β is 1.47 and 1.46 at 100 and 160 μm, respectively; and N scan is the number of PACS observations contributing to the pixels in the aperture (two in most cases, but four in the regions in which the quadrants overlap). If the noise was completely uncorrelated between pixels, the value of β should be one. The fact that β > 1 shows the limitation on the accuracy of photometry on large scales produced by the noise characteristics of PACS (Section 4.1). Fig. 10 shows signal to noise plotted against aperture radius for the three GAMA fields. At both wavelengths, the signal to noise is at a maximum at very small apertures (r < 8 arcsec). Thus, the way to produce PACS photometry with the best possible signal to noise would be to carry out photometry using this aperture, and then use the Encircled Energy Function (EEF, Section 4.2) to correct the photometry to our standard reference radius of 1000 arcsec. This argument only holds if one knows precisely the position of the object for which one wants PACS photometry; if there is significant uncertainty in the position, using a very small aperture may lead to some of the flux being missed, which is the position we face when using the positions of SPIRE sources.
In practice, because of the uncertainties in the positions of the SPIRE sources, we did not use such a small aperture. Instead, we use an aperture with a radius given by the following relationship: We used the value for the FWHM obtained from our empirical PSF (see Section 4.2): 11.4 and 13.7 arcsec at 100 and 160 μm, respectively. We deliberately made the apertures smaller for the PACS observations than for the SPIRE observations because of the poorer sensitivity of the PACS images. To improve the accuracy when the aperture sizes were small relative to the pixel sizes, we divided each pixel into four, assigning one fourth of the flux density of the real pixel to each subpixel. If a source had a reliable optical counterpart (R > 0.8, Section 5.2), we centred the aperture on the optical position and used the value of ISOA for the counterpart. If the source did not have an optical counterpart, we centred the aperture on the SPIRE position and set the value of ISOA to zero. In the latter case, the uncertainty in the SPIRE positions leads to a loss of flux density. Using a positional error for the SPIRE sources of 3.4 arcsec (Paper II), we estimated that to correct for this loss of flux we needed to increase the measured flux densities by 1.10 and 1.05 at 100 and 160 μm, respectively, and therefore if the source did not have an optical counterpart, we multiplied the measured PACS flux densities by these factors. We corrected all flux densities to our reference radius of 1000 arcsec using our empirical PACS EEF (Section 4.2). We used equation (9) to estimate the error for each flux-density measurement, scaling the error to the reference radius using the EEF. The absolute accuracy of the PACS flux scale is 5 per cent (PACS Observer's Manual). The reproducibility of measurements of the flux density of individual point sources is better than 2 per cent. To be conservative, we combine these errors linearly, giving a calibration uncertainty on our flux densities in the two PACS bands of 7 per cent. As for SPIRE (see Section 5.2), all our measurements of flux density are based on the assumption that the flux density, F ν , of a source has the spectral dependence F ν ∝ ν −1 . Most sources have very different SEDs and in Appendix A3 we describe how corrections should be made to flux densities in the catalogues to obtain flux densities for more realistic SEDs.

Overview
We need simulations to determine the completeness of the catalogues and the real errors in the flux densities, the random errors but also the systematic errors (flux bias, Section 1). Although flux bias and Eddington bias (see Section 1), the systematic errors on the number counts, affect all astronomical surveys (Hogg & Turner 1998), they have always been a particular problem in submillimetre catalogues because the submillimetre source counts are much steeper than the number counts in other wavebands (Clements et al. 2010;Oliver et al. 2010). Quantitative estimates of flux bias in ground-based surveys at 850 μm have a range of 20-40 per cent of the measured flux densities (Eales et al. 2000;Scott et al. 2002;Coppin et al. 2005).
The approach we followed for the catalogue we released at the end of the Herschel SDP (Rigby et al. 2011) was to use a theoretical model of the galaxy population to generate artificial submillimetre catalogues, which we then compared with the real catalogues. The drawback of this is that the method depends on the model being a good representation of the galaxy population. In practice, the model reproduced well the slope of the submillimetre number counts but not the overall normalization of the counts. A second disadvantage of this kind of method is that its complexity makes it very difficult for anyone else to reproduce.
Our approach in this paper is simpler and more empirical. We have broken down the problem into two parts. In this section, we describe simple simulations in which we place artificial sources on the real images and then try to find them with our source-extraction algorithm (Section 5.1). These simulations allow us to determine, for a known true flux (F t ), the conditional probability distribution of the measured flux (F m ): P(F m |F t ). From the simulations, we derive analytic distributions of P(F m |F t ) for a large range of F t in the three SPIRE bands. The simulations also allow us to estimate the completeness of the 250 μm catalogue as a function of true flux density (F t ) by measuring N d /N i , the ratio of the number of detected sources to the number of injected sources.
This simple approach does not, however, allow us to estimate the more observationally useful function which is the completeness of the survey as a function of measured flux density. It also does not allow us to estimate the flux bias in the three wavebands, the completeness in the two longer wavebands, and the effect of Eddington bias. All these questions form the second part of the problem, which we address in Sections 7 and 8.

The method
We started our simulations with an empty image (no noise or sources) with 1-arcsec pixels. We then added sources to this image at positions that lay on a grid. The grid spacing was sufficiently large that any overlap between injected sources was avoided and that when these artificial sources were added on to the real images (see below) they did not affect the statistics of the images, especially the confusion noise. We added a small random scatter to the injected positions, so that when they were added on to the real images they did not fall at the same position within a real pixel.
We then created three images from the high-resolution maps by convolving them with the PSF for the three SPIRE bands and then re-binnning these images to the actual pixel sizes of 6, 8 and 12 arcsec at the three SPIRE wavelengths. As a 'sanity test' on our software, we first injected 100 1 Jy sources and ran MADX (see Section 5.1) to find these sources on these noiseless maps, using the same settings we used for the real images. The mean of the 100 flux densities we measured was 0.991, 0.987 and 0.990 Jy at 250, 350 and 500 μm, respectively, with a scatter of less than 1 per cent. This error is caused by the fact that sources are not generally found at the centre of a pixel. The bicubic interpolation in MADX partially corrects for this effect but does not do it perfectly. We used these factors to correct the flux densities of the real sources measured by MADX (Section 5.1). Table 4. Statistics of the In-Out simulations. For each 250 µm flux bin, the table reports the ratios of the number of recovered sources to the number of injected sources, N d /N i , averaged over the whole survey and for the deepest areas only, and the parameters of the shape of the conditional probability distribution, P(F m |F t ), described as two Gaussians with different standard deviations (see Section 6.3, equation 11). We then ran simulations by adding artificial sources on to the real images. We created noiseless images with point sources covering a wide range of flux densities (see Table 4). The flux densities were chosen in order to be equally spaced on a logarithmic scale and cover a range of ∼1-50 σ tot Gauss . About 3000 sources for each flux density were injected. We then added the noiseless images to the real maps.

Double-Gaussian Fits
We ran MADX on these images, the real maps plus the artificial sources, using exactly the same procedure we used for the real data. We then compared the new MADX catalogue with the catalogue obtained before the artificial sources were added, in order to find which artificial sources were detected. We used a radius of 12 arcsec to look for matches between the sources in the new catalogue and those in the old catalogue and between those in the new catalogue and the positions of the injected sources. If there was a match between a source in the new catalogue and the position of an injected source but not with a source in the old catalogue, it was clear that we had simply found the injected source, allowing us to compare the measured flux density and position with the true (injected) flux density and position. An ambiguity occurred if there was a source in the new catalogue that matched both the position of an injected source and a source in the old catalogue. In this case, we made a decision by comparing the injected flux density with the flux density of the matched source in the old catalogue. If the injected flux density was larger than the flux density of the source in the old catalogue, we concluded that the source in the new catalogue was a match for the injected source. Otherwise, we concluded that we had not detected the injected source. Hogg (2001) used the same criterion in interpreting his simulations of the submillimetre sky. Fig. 11 shows histograms for the recovered injected sources of the measured flux density, F m , for three injected flux densities, F t , at 250 μm. About 3000 sources were injected for each true flux density, although the number of sources in each histogram is less than 3000, because not all sources were recovered.

The conditional probability distributions
Each histogram provides an estimate of the conditional probability function, P(F m |F t ), at a particular value of F t . However, the histograms are quite noisy and to provide a less noisy estimate of this function, we have fitted an analytic function to the histograms. The distributions are asymmetric and they follow the shape of the pixels distribution of the map with a suppression for faint fluxes due to the 2.5σ limit in the source detection: we found that a good fit is obtained by fitting the following function, which represents the two sides of the histogram by Gaussians with different standard deviations, σ 1 and σ 2 , centred inF and normalized to N 0 : Table 4 lists the parameters of the distribution for each value of F t . It also lists N d /N i , the ratio of the number of recovered sources to the number of injected sources.

The errors in the flux densities
In Section 3.4, we concluded that a good estimate of the error in the flux density for a source is more likely to come from the variance of the image than from fitting the negative part of the histogram of pixel values. Following equation (5), we now write the variance for a source of flux density, F s , as  Table 4.
This means that each source will have a different confusion noise, depending on its flux F s . The justification for the use of a different noise for each source is that pixels brighter than the flux density of the source cannot contribute to the calculation of the confusion for that source. If a brighter source was also in the pixel, we would have assigned all of the flux in the pixel to that source; the fainter source would contribute to the confusion noise of the brighter source, but it would not have been detected as a source in our catalogue. This approach has been suggested before (Crawford et al. 2010) and also takes into account, indirectly, the variation of confusion noise within the map: in a more confused region, the flux of the detected sources will probably have a larger contribution from undetected sources, thus it will be larger and, according with equation (12), also their confusion noise will be larger. On different premises, also Leiton et al. (2015) have recently calculated a 'customized confusion error' for the sources in the GOODS-Herschel survey, using source density arguments instead of the flux-density argument used here. We now test if equation (12) is a good estimate of the uncertainty in the flux density using the results of the In-Out simulations. These simulations are as close to reality as we can make them, since we are injecting artificial individual sources on to the real images, running our detection software on the images, and comparing the measured flux densities with the injected flux densities. As our estimate of the error on the flux density from the simulations, we combine the parameters from the double-Gaussian fits to give The circle points in Fig. 12 show this quantity plotted against flux density. The variation of σ sim with flux density does in general support the idea that the error in the flux density of a source does depend on the flux density of the source. The squares show σ tot Var calculated using equation (12). The results from the In-Out simulations and from the variance measurements agree well at a flux density of ∼30 mJy, but σ sim does not increase with flux density at brighter flux densities in the way that σ tot Var does, and it falls off faster with decreasing flux density than σ tot Var .
The flattening of σ sim at bright fluxes is probably the result of small-number statistics together with the very non-Gaussian pixel distribution of the maps. We would expect to see an increase in the average σ sim with flux density, because the noise measured for an injected source comes from all pixels fainter than itself, but not  (13). The squares show σ tot Var measured from the images using equation (12). The upturn at low flux densities is due to the fact that F map used in equation (12) is calculated on the whole map so, for faint F s , it is larger than the flux limit used to calculate σ tot Var . The continuous line shows the assumption we have made about the relationship between flux-density error and flux density to derive the flux-density errors in the catalogue. The dashed line (equation 14) shows the confusion noise that is predicted from this relationship, after using equation (4) to remove the contribution of instrumental noise. The vertical dotted line shows the 4σ catalogue limit at 250 µm. from brighter pixels (Section 6.2); as the upper threshold of pixels that contribute to the noise is increased, the non-Gaussian positive tail (Fig. 4) increases the variance. In the limit of infinite number of injected sources we would recover the error distribution expected from the pixel distribution of the map, but since there are only a very small number of bright pixels in the map, it is very likely that our simulations do not have any sources on bright pixels, and so the estimated variance is low. These rare cases of very high errors lead to a high average error which is not appropriate for a 'typical' source, so we chose to use the errors from the simulations, as these are more representative of the 'typical' error for a source.
The second difference is that the apparent noise in the simulated source fluxes falls rapidly for fluxes fainter than 20 mJy, and this is due to the fact that σ sim is calculated from only those injected sources that are detected. A measured source has to be brighter than 2.5σ (Section 5.1), so we recover only a narrow tail of a broad error distribution. The simulations show that the scatter in the flux density of a source that is just detected above the 2.5σ threshold is extremely low, but this does not imply that the uncertainty is small. When we modelled this effect by measuring the variance on an image with an upper and lower flux-density limit, we found values that followed the results of the In-Out simulations much more closely than when the lower limit is not used. A better estimate of the uncertainty for a faint source is simply the total noise as estimated from the Gaussian fit of the pixel histogram (Table 3).
Based on the comparison of the variance measurements and the results of the In-Out simulations, we adopted the following approach for estimating the error in a source's flux density. We assume that the variance is constant at flux densities 30 mJy and linear below this flux density. This relationship is shown by the solid line in Fig. 12. The justification for the constancy above 30 mJy is that this is what the simulations tell us the error is. At lower flux densities, the simulations will drastically underestimate the true flux-density errors because they only give us information about the dispersion in the flux densities of the detected sources. The relationship we assume below 30 mJy is a conservative one, since while it follows the behaviour of the simulations at brighter flux densities, it is systematically higher than it is at lower flux densities, where it follows the variance measurements. We then use our estimate of the instrumental noise (see Section 3.3) and equation (4) to obtain an estimate of the confusion noise for a source of flux density F s , which is given by in which the flux density of the source is given in Jy. The actual uncertainty in the flux density of an individual source depends on the source's position because the instrumental noise is not uniform. We assigned an error to the flux density of each source in the catalogue (see Section 10) by using equation (14) to estimate the confusion noise and then adding this in quadrature to the instrumental noise given by the map of the instrumental noise (Section 3.3).
At 350 and 500 μm, the values for the errors in the flux density estimated from the In-Out simulations are fairly constant (8.0 and 8.5 mJy at 350 and 500 μm, respectively) for input flux densities above 15 mJy. At the two longer wavelengths, we are simply measuring the brightness of a source that was detected at 250 μm, and so we expect that the relationship between flux-density error and  Fig. 13. The green dashed line show the prediction from the relationship given in text if we assume that the error on the flux density of the source is given by σ tot Gauss , the red dotted line is the predicted relationship if we assume that the flux-density error is given by σ tot Var . The vertical black dotted and black dashed lines show, respectively, the 2.5σ detection and 4σ catalogue limits. flux density should be much weaker than at 250 μm. For this reason, we have assumed that the contribution to the error in flux density from confusion is a constant. Subtracting the average instrumental noise in quadrature from the errors above (see Tables 2 and 3), we obtain a confusion noise of 6.59 and 6.62 mJy at 350 and 500 μm, respectively. As at 250 μm, we estimate the total uncertainty on the flux density of each source in the catalogue (see Section 10) by adding the confusion noise in quadrature to the instrumental noise given by the map of the instrumental noise (Section 3.3).

The positional uncertainties
The In-Out simulations can also be used to investigate the accuracy of the source positions. Fig. 13 shows histograms of the differences between the true positions of the injected sources and the positions measured with MADX for three different injected flux densities. The histograms are well fitted by a Gaussian. Fig. 14 shows the standard deviation of the Gaussian as a function of injected flux density. As  expected, the accuracy of the measured positions decreases with decreasing flux density.
Theoretically, we expect the differences to follow a Gaussian distribution with the standard deviation of the Gaussian given by (0.6×FWHM)/SNR, where FWHM is the full width at halfmaximum of the PSF and SNR is the signal-to-noise ratio (Ivison et al. 2007). The FWHM at 250 μm, the wavelength at which the positions are measured in MADX, is 18 arcsec. The green lines in Figs 13 and 14 show the theoretical prediction when σ tot Gauss , the noise calculated from the Gaussian fitting technique (Section 3.4), is used in the model. The red lines shows the theoretical prediction when σ tot Var , the variance in the image, is used. The actual distributions fall between the two predictions.
In the accompanying data-release paper (Paper II), we estimate the accuracy of the positions by comparing the positions of the H-ATLAS sources with the positions of galaxies on the SDSS. Since this technique is likely to be even closer to ground truth than the In-Out simulations, we prefer the estimates of positional accuracy given in that paper, although the estimates of the positional accuracy from the In-Out simulations are actually very similar.  Table 4) plotted against injected flux density. This figure represents the completeness of the survey but as a function of true flux density rather than measured flux density.

Completeness at 250 µm
We also carried out versions of the In-Out simulations restricted to the small deeper regions of the maps in which the quadrants overlap (see Fig. 2 to visualize this) and also to the part of the map outside these regions. Fig. 15 shows the completeness versus flux density for these regions, made from four 9 h data sets and two 9 h data sets, respectively (see Section 2). As one would expect, the results for the 'two-scan' regions are almost exactly the same as for the survey as a whole, because most of the survey consists of these regions. The completeness for the deeper 'four-scan' regions is slightly higher than for the survey as a whole, again as one would expect. Table 4 also lists the values of the N d /N i for these deeper regions. Although we do not use these values in the rest of the paper, they give enough information for anyone in the future who is interested specifically in these regions to carry out the analysis of flux bias and completeness that is described in Section 8.
The measured flux densities are likely to be systematically different from the true flux densities as the result of flux bias. Therefore, to derive the completeness of the survey as a function of measured flux density, which is sometimes what one needs from a practical point of view, we need also to be able to estimate the flux bias of the survey. This is not possible from the In-Out simulations by themselves. Although we cannot derive the flux bias and completeness of the surveys directly from the In-Out simulations, they do provide an important stepping stone. The flux bias can be represented by the conditional probability of the true flux density given the measured flux density, P(F t |F m ). The In-Out simulations, however, provide an estimate of the conditional probability of the measured flux density given the true flux density, P(F m |F t ). To go from one to the other requires knowledge of the number counts of sources down to a flux density well below the measured flux densities in the H-ATLAS catalogue. Fortunately, the deeper surveys with Herschel have provided this information. In Section 7, we use the conditional probability distributions provided by the In-Out simulations to show that the number counts of the sources from our survey are consistent, over the range of flux density covered by our catalogue, with the results of the deeper surveys in other parts of the sky. In Section 8.1, we use the results from the deeper surveys to estimate the flux bias.

T H E S O U R C E C O U N T S
We can combine the conditional probability distributions from the In-Out simulations as a matrix, P, each element of which, P ij , is the probability that a galaxy with true flux F i is detected with a measured flux F j . This matrix contains the same information as the conditional probability distribution, and it also contains information about the incompleteness of the survey, which is given by the recovered fraction of sources in the In-Out simulations; the sums of the columns of the matrix are the recovered fractions listed in Table 4.
Let us now represent the number of measured sources in each bin of flux density as a vector,N m , and the mean number of sources predicted, from the true source counts, to fall in these bins as another vector,N t . We calculate the number of sources in each bin, given that the flux densities we used in the simulations (see Table 4) are the centre of the bins and that the bin limits are half way between two contiguous fluxes. We can then relate the two vectors by the matrix equation: This equation expresses in matrix form the combined effects of incompleteness and Eddington Bias on the intrinsic number counts. The matrix P contains all the information we need about the effects of the noise properties of our images and of the algorithm we have used to find the sources. We want to know the true source countsN t , so in principle we should simply invert equation (15). Unfortunately, the elements of the matrix are not perfectly known, and simply inverting the matrix amplifies the errors in the matrix elements, producing a very uncertain and unstable solution (this is analogous to a deconvolution, in which errors in the deconvolution kernel lead to an increase of noise in the solution).
To solve the problem, we need to use some extra prior information about the source counts. It is not necessary to assume any strong hypothesis, we can simply introduce what is called a stabilizing functional or regularizing operator. For example, one assumption that is reasonable to make is that we expect the source counts to have a smooth form.
A detailed account of the theory of linear regularization and its implementation is given in Numerical Recipes (Press et al. 1982, Section 18.5). Using the method of Lagrange multipliers, we determine our vectorû =N t by minimizing the quantity: In this equation, A and B are two positive functionals ofû and λ is a Lagrange multiplier. The first functional comes from minimizing the Chi-squared difference between the observed source counts and the model, and is given by The elements of the matrix A are A ij = P ij /σ i , with σ i being the Poisson error on the number of sources with measured flux densities F i . The vector b contains the observational data and its elements are b i = N m, i /σ i , in which N m, i is the number of sources with measured flux densities F i . The second functional in equation (16) represents our prior assumptions about the source counts and is given by The matrix H is based on our assumptions about the source counts. In our analysis we assume that the logarithm of the source counts are approximately linear, and we therefore minimize the second derivatives of the counts. With this assumption, H = B T · B, in which B is given by (see section 18.5 of Numerical Recipes) The value of λ will give different weights to the two parts of the minimization. When λ = Tr(A T · A)/Tr(H), the two functionals have comparable weights.
The minimum value of A + λB occurs when The standard uncertainties on the elements ofû are given by the diagonal terms of the covariance matrix: The covariance between the elements of the solution strongly depends on the choice of λ: the larger the weight given to the prior, the larger will be the covariance. We do not want our results being too dependent on the choice of the prior and we also want to reduce the covariance between the elements ofû. For these reasons, we have chosen the minimum value of λ which gives a solution in which all the elements ofû are positive, and thus a solution that is physically reasonable. We never need the weight of the prior to be larger than 1/100 of the weight of the data to find a solution that satisfies this condition.  (Béthermin et al. 2012). They underestimate the counts at the bright end because the HerMES area is smaller than the H-ATLAS one, so they do not have adequate statistics above ∼100 mJy. The vertical dashed line shows the 4σ limit of the H-ATLAS catalogue (see Section 10). (19) is only one of the possible priors that can be used. We also tried a different prior (minimizing the first derivative of the counts) and the results did not change significantly, due to the small weight we give to the prior.

The matrix B in equation
This approach to retrieve the number counts shows some similarities with Reduction C used for the 850 μm survey SHADES (see section 5.2 from Coppin et al. 2006), which in turn was built on the methods used by Borys et al. (2003) and Laurent (2005) for the analyses of SCUBA and Bolocam data, respectively. The new element we have introduced in this study is the regularizing operator. The use of a regularizing operator is quite common in imagine reconstruction, for example in the study of lensed objects (see for example Warren & Dye 2003), but, to our knowledge, this is the first time it has been used to derive the source number counts of a galaxy population. We postpone a more exhaustive analysis of this technique to future papers.
The source number counts obtained from the matrix inversion are shown in Fig. 16 and Table 5. The coloured symbols in the plots show the observed number counts,N m , in each of the three GAMA fields, while the coloured lines show the true number counts for each field,N t , estimated using our inversion technique. Points are slightly correlated because of the covariance between flux bins.
The uncertainties are derived by Monte Carlo simulations of synthetic data sets: we generate 10 000 times a vector of 'measured' number counts,N m , assuming a Poisson uncertainty on the real N m , and we apply the inversion procedure each time, obtaining 10 000N t . The scatters in the results are the uncertainties shown in Fig. 16 and listed Table 5 and may be up to a factor 2 larger than the nominal uncertainties calculated from the covariance matrix (see equation 21). The quoted uncertainties do not include the correlation between flux bins. Because our counts are derived from simulations only and do not assume any prior on the existing galaxy population, we can only trust them for fluxes brighter than ∼20 mJy. At fainter fluxes, the conditional probability distribution P(F m |F t ) is too uncertain. The number counts for the three GAMA fields are quite similar except at the brightest flux densities, where the difference is due to cosmic variance. The observed counts differ quite markedly from the source counts measured by the HerMES team from their much deeper observations, both the source counts inferred from a P(D) analysis (Glenn et al. 2010) and the source counts measured by a stacking analysis (Béthermin et al. 2012). However, after we have applied our inversion technique, the corrected source counts show consistency with the HerMES measurements, except at the faintest flux densities, well below the flux limit of our catalogue.
In Section 8.2, we show that it also possible to apply this method, with some modifications, to the source counts at 350 and 500 μm.

C O R R E C T I O N S F O R F L U X B I A S A N D C O M P L E T E N E S S
By injecting artificial sources on to the images with a low enough surface density that they do not affect the statistical properties of the images and by using the real data-reduction pipeline, the In-Out simulations (Section 6) come as close as possible to the ground truth of how real sources on the sky are turned into sources in the H-ATLAS catalogue. In Section 7, we showed that this information can be used directly to correct the 250 μm source counts for Eddington bias. However, we also need to know the completeness of the survey and the flux bias correction for individual sources at all three wavelengths. The flux bias problem is an inverse problem, requiring approximations and prior assumptions about either the submillimetre sky (e.g. Rigby et al. 2011) or about the submillimetre source counts (Crawford et al. 2010). This section describes our best current attempt to estimate the flux bias, but we recognize that it may be possible to improve our analysis in the future, if only because of improved knowledge of the submillimetre sky. We recommend that anyone interested in improving our estimates should start from the ground truth provided by the In-Out simulations.

Flux bias and completeness at 250 µm
There are analytic methods to correct for flux bias (e.g. Hogg & Turner 1998;Crawford et al. 2010) and also ones based on simulations (e.g. Coppin et al. 2005). In either case, it is necessary to make assumptions about the underlying source counts down to flux densities well below the detection limit of the catalogue because it is these faint sources which produce the majority of the confusion noise. In the previous section, we showed that the 250 μm source counts in the H-ATLAS fields, after a correction for Eddington bias, are in good agreement above the H-ATLAS detection limit with the much deeper source counts of Glenn et al. (2010) and Béthermin et al. (2012). In this section, we estimate the flux bias at 250 μm, using these deeper source counts to provide a prior probability distribution for the flux density of a source. We emphasize that we only use these deeper source counts as a prior probability distribution, and the method does not require the assumption that the source counts in the H-ATLAS fields below the catalogue limit are exactly the same as those measured in the deeper surveys.
We have used the method of Crawford et al. (2010). In this method, the measured flux density in a pixel is F m and the true flux density of the brightest individual source in that pixel is F t . Given a measured flux density in a pixel, the probability of the true flux density of the brightest source in that pixel is then in which P(F m |F t ) is the likelihood of measuring a flux density equal to F m in a pixel given that the brightest source in that pixel has flux density F t , and P(F t ) is the prior probability that the brightest source in that pixel has flux density F t . The prior probability, P(F t ), is the differential source counts with an exponential suppression at low flux density (see section 2.1 of Crawford et al. 2010 for more details). In calculating the prior probability distribution, we use the differential 250 μm source counts given in Béthermin et al. (2012). The likelihood, P(F m |F t ), can be written using a Gaussian likelihood approximation (Crawford et al. 2010): in which F is the mean of the image, which in general is not precisely zero (see Section 3.1). The value of σ we use for each source is the uncertainty in the flux density of the source, derived individually for each source using the method of Section 6.4. In this way, we produce an individual estimate of the flux bias for each source. Fig. 17 (top) shows the probability distributions P(F t |F m ), P(F m |F t ) and P(F t ) for a pixel with a measured flux density equal to the catalogue limit of 29 mJy. We estimate the true flux density of the brightest source in the pixel from the maximum in the posteriori probability distribution P(F t |F m ). The bottom panel in Fig. 17 show the mean flux bias factor (F m /F t ) plotted against the measured flux density F m . We have calculated the mean flux bias for each bin of F m by estimating F t for each pixel that falls in this bin and then calculating the mean value of F m /F t for all the pixels in the bin. This relationship is listed in Table 6.
We have tested the results of our flux-bias analysis using the In-Out simulations (see Section 6). 12 In the top panel of Fig. 18, we have plotted the mean value of the measured flux density against the injected flux density of the artificial sources in the simulations. Figure 17. The results of our analysis, using the method of Crawford et al. (2010), of the flux bias at 250 µm. In this analysis, the true flux density of the brightest source in a pixel is F t and the measured flux density in the pixel is F m . Top: probability distributions for a pixel with F m = 29 mJy: P(F t |F m ), P(F m |F t ) and P(F t ) (see the text for more details). Note that the probabilities in this figure have not been normalized and have been scaled so that the shapes of the three distributions can be easily compared. Bottom: flux bias plotted against measured 250 µm flux density. We have calculated the flux bias in each bin of measured flux density by calculating the mean value of F m /F t for the sources in this bin. This relationship is listed in Table 6. The vertical black dotted and black dashed lines show, respectively, the 2.5σ detection limit (Section 5.1) and the 4σ catalogue limit (see Section 10).
The red lines show a linear fit to these mean flux-density values. As expected, the measured flux densities are systematically higher than the injected flux densities, particularly at fainter flux densities. The bottom panel shows same plot after we have used the flux-bias corrections, shown in Fig. 17 and listed in Table 6, to correct the measured flux densities. The measured flux densities are now much closer to the injected flux densities. Table 6 shows that the flux bias at the flux-density limit of the catalogue is approximately 20 per cent. This is larger than the value of 6 per cent, which we estimated for our SDP catalogue (Rigby et al. 2011). In our SDP analysis, we were forced to use a theoretical model of the submillimetre sky to estimate the distribution of true flux densities, whereas now we have had the major advantage of being able to use the deeper Herschel surveys to produce a prior Table 6. Flux biases. These numbers are the ratio between the measured and the true flux density of a source (F m /F t ). In order to derive the corrected flux densities, the flux values reported in the released catalogues need to be divided by the numbers in this   Table 6. probability distribution for the true flux densities, and so our new estimates supersede the former ones. The completeness of the survey shown in Fig. 15 was derived from the In-Out simulations and is the completeness of the survey as a function of true flux density. Our derivation of the flux bias now allows us to calculate the completeness of the survey as a function of measured flux density. First, we use the flux-bias relationship shown in the bottom panel of Fig. 17 (listed in Table 6) to estimate the true flux density of each source in the catalogue. We then use this estimate of the true flux density and the completeness estimates shown in Fig. 15 (listed in Table 4) to estimate the probability that this source would have been detected in the survey. We then derive the completeness for each bin of measured flux density by calculating the mean probability value for all the sources in the catalogue in this bin. These completeness estimates are given in Table 7 and shown in Fig. 21. In Section 7, we corrected the observed H-ATLAS counts for the effect of Eddington bias using a method that did not require us to make corrections to the flux densities of individual sources. After making this correction we found good agreement between the H-ATLAS source counts and the HerMES source counts (Béthermin et al. 2012) at flux densities 30 mJy. We can now derive the source counts in a different way, using the results of our flux-bias analysis and the estimates of completeness from the In-Out simulations. First, we use the flux-bias relationship (listed in Table 6) to estimate the true flux density of each source in the catalogue, which allows us to calculate the source counts as a function of true flux density. We then correct these source counts for incompleteness, this time using the completeness values in Table 4, which give the completeness as a function of true flux density. Fig. 19 shows the corrected source counts compared with the ones from Béthermin et al. (2012). Once again there is good agreement between the HerMES and H-ATLAS counts.
The two methods to derive the number counts described in Sections 7 and 8 give quite similar results, but they have different positive and negative aspects. The matrix inversion method requires no prior and take into account all the effects included in the real maps that can influence the measure of the flux density: residual cirrus emission, large extended sources, clustering (even if this information is partially lost during the simulation). Unfortunately, running the In-Out simulation is quite computationally expensive and the accuracy of the conditional probability distribution is not good enough at the faintest flux densities to provide reliable number counts. On the other side, correcting the flux density of the single Figure 19. The number counts at 250 µm plotted against flux density. We first corrected the measured flux density of each source for flux bias using the relationship in Fig. 17 (bottom panel) and listed in Table 6. We then used this estimate of the true flux density and the completeness estimates from the In-Out simulations (Table 4) to calculate the number counts. The circles show the average source counts for the three GAMA fields after these corrections for flux bias and incompleteness have been made. The error bars take into account uncertainties in the corrections and the scatter between the fields. For comparison, the dotted line shows the numbers counts for the HerMES survey from Béthermin et al. (2012) and the crosses are the counts derived in Section 7 and reported in Table 5 averaged among the three fields. The vertical black dotted and black dashed lines show, respectively, the 2.5σ detection (Section 5.1) and 4σ catalogue limits (Section 10).
sources it is possible to derive reliable number counts down to the detection limit of the survey, but requires the assumption of a prior. We let the reader to decide which method to use, according to the scientific goals pursued.

Completeness at 350 and 500 µm
The completeness of the catalogue at 350 and 500 μm must be derived with a different method because the 350 μm and 500 μm sources were actually detected on the 250 μm images. We have used the following empirical technique to estimate the completeness of the catalogues at these wavelengths.
Whether a source is detected should only depend on its true 250 μm flux density and not on its 350 μm flux density. We take the measured 250 μm flux density of each source in the catalogue, correct it for flux bias (Table 6) and use the completeness estimated from the In-Out simulations for the true 250 μm flux densities (Table 4) to calculate the probability this source would have been detected by the survey. In each bin of 350 μm flux density, we then take the average value of these probabilities as our estimate of the completeness in this bin. We follow the same procedure at 500 μm.
We used the results of the In-Out simulations (Section 6) to perform a simple sanity test of this technique. Although the distribution of 350/250 μm flux-density ratios in the simulations was not the same as in the real sky, 13 our correction technique should work equally well. In the simulations, the distribution of this ratio for the detected sources was quite different from the distribution for the injected sources, because sources with high values for this flux ratio, and thus low values of 250 μm flux density, will tend not to be detected. However, after we applied the correction for flux bias at 250 μm and then the corrections for incompleteness from Table 4, we recovered a distribution of 350/250 μm flux-density ratios that was similar to the distribution for the injected sources. Fig. 20 shows how this technique is applied to the real data. The figure shows histograms of observed 250 μm flux densities for  Table 7, plotted against measured flux density. The vertical lines show the 4σ catalogue limit at each wavelength, 29, 37 and 40 mJy at 250, 350 and 500 µm, respectively. The vertical blue dotted line shows the 2.5σ detection limit at 250 µm (Section 10).  (Table 6). The red dashed line shows the effect of the second correction, after we have corrected for the incompleteness at 250 μm (Table 4). The completeness in each 350 μm flux-density bin is the ratio of the number of sources under the black line to the number of sources under the red line.
Our completeness estimates are listed in Table 7 and shown visually in Fig. 21. The reason why completeness does not decline so rapidly with decreasing flux density as at 250 μm is that a source may be very faint at 350 or 500 μm but still be bright enough at 250 μm to be easily detected.
In Section 7, we described an inversion technique for correcting the source counts at 250 μm for Eddington bias. We can now use the completeness estimates at 350 and 500 μm to extend this method to the two longer wavelengths. As before, we represent the results of the In-Out simulations as a matrix, P, in which each element, P ij , is the probability that a galaxy with intrinsic flux F i is detected with a measured flux F j . At 250 μm, the matrix contains all necessary information about the effect of noise on the flux density of a source and on the incompleteness of the survey, with the sum of the elements in each column of the matrix giving the incompleteness of the survey for each injected flux density.
At the two longer wavelengths, however, the sum of the columns of each matrix does not give the true completeness of the survey because the detection of the sources was performed at 250 μm. 14 Nevertheless, we can still apply the inversion method using a simple modification.
First, we normalize each column of the matrix, so that the sum of the elements in the column is unity, effectively removing the 14 The sum of the columns does not give the true incompleteness because the distribution of the 250/350 µm flux-density ratios used in the In-Out simulations was not the same as in the real sky. We could have used a more realistic distribution in the simulations but chose not to because we could not have been sure how close our assumed distribution was to the real distribution, since the distribution we observe in the survey is strongly distorted by the effects of completeness and flux bias.  Béthermin et al. (2012). They underestimate the counts at the bright end because the HerMES area is smaller than the H-ATLAS one, so they do not have adequate statistics above ∼100 mJy. The vertical dashed line shows the 4σ limit in that band. Top: counts at 350 µm. Bottom: counts at 500 µm.
incorrect completeness information. As before (see Section 7), the observed source counts are represented by a vector,N m . We now correct this vector using the completeness estimates from Table 7. We now have an equation similar to equation (15), linking the observed counts and the true countsN t : in which P is the same as P, but with sum of the columns normalized to unity, andN m corr isN m divided by the completeness at 350 or 500 μm.
The results at 350 and 500 μm are shown in Fig. 22 and listed in Tables 8 and 9. Uncertainties are estimated by Monte Carlo simulations, in the same way as explained in Section 7 and do not include the correlation between flux bins. Because our counts are derived from simulations only and do not assume any prior on the existing galaxy population, we can only trust them for fluxes brighter than ∼20 mJy. At fainter fluxes, the conditional probability distribution P(F m |F t ) is too uncertain. As at 250 μm (see Fig. 16), the corrected H-ATLAS counts at 350 μm are in good agreement with the HerMES source counts above the flux-density limit of our catalogue. At 500 μm, the corrected H-ATLAS source counts are a little higher than the HerMES source counts but consistent within the errors.

Corrections for flux bias at 350 and 500 μm
Determining the best method to correct for flux bias at 350 and 500 μm requires careful thought because of the way we have carried out the survey. Crawford et al. (2010) have presented an elegant method for correcting for flux bias when a survey is carried out at more than one wavelength. However, in this method the wavelengths are treated as equivalent, and a simple thought experiment shows that this method will not produce the right answers for our survey. Let us consider observations at 250 and 350 μm. We detect the sources at 250 μm, selecting peaks in the map, and thus the flux densities of the sources at 250 μm will be affected by the flux bias caused by the noise on the map (instrumental plus confusion). We described in Section 8.1 how we estimated the size of this bias. Now let us consider the 350 μm observations. In our survey, we do not have to detect the source above a flux-density threshold at 350 μm. What we actually do is go to the position of the 250 μm source and measure the 350 μm flux density at this position. Since the instrumental noise on the 350 μm image is uncorrelated with the instrumental noise on the 250 μm image, the instrumental noise on the 350 μm image is equally likely to produce a positive or negative fluctuation on the map. Therefore, in this thought experiment, there should be no flux bias at 350 μm due to instrumental noise. The confusion noise, on the other hand, is highly correlated between the two bands, because we expect the sources producing confusion to emit at both wavelengths, 250 and 350 μm. The flux bias at 350 μm therefore must depend on the flux bias at 250 μm and the covariance between the 250 and 350 μm images.
In conclusion, using a pure intuitive approach, we would expect the differences between measured and true fluxes at 350 and 500 μm to be proportional to the fraction of the contribution of the confusion noise to the total noise at 250 μm, to the beam size and to the average colour of the sources.
If we treat the problem in a more formal way, we reach very similar conclusions. We can quantify this argument in the following way. Let us suppose that a source has measured flux densities of F m, 250 and F m, 350 at the two wavelengths. We then use the analysis of Section 8.1 to estimate the flux-bias factor fb 250 = F m, 250 /F t, 250 at 250 μm. The source will therefore be sitting on a fluctuation on the 250 μm image, the combination of instrumental and confusion noise, given by We make the same assumption as with the analysis of the flux bias at 250 μm that the fluctuations on the two images follow a Gaussian distribution, with the simultaneous probability of the fluctuations at 250 and 350 μm following a bivariate Gaussian distribution: in which r is the vector ( 250 , 350 ) and C is the covariance matrix. If we define the magnitude of the fluctuation under the source at 250 μm using equation (25) (27) in which Covar is the covariance between the two images and Var 250 is the variance of the image at 250 μm. Note the important point that this analysis does not require any assumptions about the proportions of the variance that come from instrumental noise and confusion. For a source with measured flux densities of F m, 250 and F m, 350 , we estimate the variance and the covariance from the matched-filterconvolved image using the following relationships: for all F i,250 < F m,250 , F i,350 < F m,350 .
We only include pixels with flux densities less than the measured flux densities in the calculations for the same reason as in Section 6.4; a pixel brighter than the flux density of a source should not be used in the calculation of the variance or covariance because  Table 6. Error bars take into account differences in covariance between the fields. The vertical lines show the 4σ catalogue limit at each wavelength, 37 and 40 mJy at 350 and 500 µm, respectively (Section 10).
otherwise we would have measured a higher flux density for that source. The flux bias at 350 μm is then calculated from This analysis leads to an individual estimate of the flux bias at 350 μm for each source. To generate an average flux-bias relationship, we calculate the mean value of the flux bias in bins of 350 μm flux density. This relationship is shown in Fig. 23 and is listed in Table 6. The figure and table also include the relationships between average flux bias and measured flux density at 500 μm, calculated in the same way. As for the 250 μm source counts (see Section 8.1), we can now correct the observed 350 and 500 μm source counts for completeness and flux bias, and then compare the corrected counts with the HerMES source counts (Béthermin et al. 2012). To correct the counts, we first estimate the incompleteness for each source using its measured 350 or 500 μm flux density and the completeness values listed in Table 7. We then correct the measured flux density of each source for flux bias using the values listed in Table 6. We then add up the number of sources in each bin of corrected (true) flux density. Fig. 24 shows the results. There is very good agreement with the corrected counts and the counts of Béthermin et al. (2012) at 350 μm. The agreement is less good at 500 μm, with the HerMES counts falling below the corrected H-ATLAS counts. We found a similar discrepancy when we corrected the 500 μm counts for the effect of Eddington bias (see Section 8.2, Fig. 22) rather than correcting the individual sources for flux bias and incompleteness, although the direct inversion method produced a smaller discrepancy. One possible explanation of the disagreement, supported by the fact that we see a similar discrepancy from both methods for reconstructing the source counts, is either that the 500 μm source counts in the HerMES and H-ATLAS fields are genuinely different or that there is some systematic error in the 500 μm source counts produced by the method used to produce either the H-ATLAS catalogue or the HerMES catalogue.
A second possibility is that the discrepancy is caused by the Gaussian approximation we have made in the flux-bias analysis, since the bright sources in the non-Gaussian tail of the distribution would increase the flux-bias factor, an effect which one  Tables 8 and 9 averaged among the three fields. Figure 25. The probability that a second source is so close to a primary source that the two sources cannot be distinguished. The x-axis gives the flux density of the second source expressed as a percentage of the flux density of the primary source, and the y-axis gives the probability that there is a second source at least this bright and close enough to the primary source that the two cannot be distinguished. The solid lines show the predictions at 250 µm (blue), 350 µm (green) and 500 µm (red) if the effect of source clustering is not included in the calculation. The dashed lines show the predictions for 350 µm (green) and 500 µm (red) if the effect of source clustering is included in the calculation. See Section 8.3 for more details. We have not including a clustering model at 250 µm because Maddox et al. (2010) did not find any significant clustering in that waveband. might expect to be greatest at 500 μm because of the larger beam size.
We refer to Section 8 for a discussion about the differences between the counts derived by the matrix inversion method and the counts obtained correcting the flux densities of single sources.
In Fig. 25, we have used the counts from Béthermin et al. (2012) to predict the probability of a second source falling so close to a primary source that the two sources cannot be distinguished. We have assumed that the maximum distance between two such sources is equal to the half width half-maximum of the PSF, and we have set the flux density of the primary source equal to the catalogue limit (Section 10) at the chosen wavelength. The solid lines in Fig. 25 show the predictions if the effect of source clustering is not included.
There is very little difference between the predictions for the three wavelengths. Therefore, the larger beam size at 500 μm makes little difference.
The dashed lines, however, show the predictions if we include the effect of source clustering. We have used the results of Maddox et al. (2010), who found no evidence for the clustering of 250 μm sources but strong clustering of 350 and 500 μm sources. We have used the values of the clustering amplitude given by Maddox et al. (2010) for the samples selected at 350 and 500 μm with no colour selection and for a value for the index of the correlation function of 0.8. The dashed lines show the results. In this case, the probability of a second bright confusing source is highest for 500 μm, so it is possible that the discrepancy in the source counts at 500 μm is caused by the breakdown of our Gaussian approximation caused by the increase clustering of sources at 500 μm.
However, whatever the cause of the discrepancy between the H-ATLAS and HerMES source counts at 500 μm, there is little practical effect on the catalogues. Our estimated flux bias at 500 μm for a source at the 4σ limit of the catalogue is 4 per cent. Increasing the flux-bias factor to 10 per cent would be enough to make the H-ATLAS and HerMES counts agree. Therefore, there is a 6 per cent systematic uncertainty in the 500 μm flux density for a source at the catalogue limit, much less than the statistical uncertainty of 25 per cent.

Comparison with SDP data
We have compared our catalogue with the one we released after the SDP (Rigby et al. 2011). The two catalogues are based on almost exactly the same Herschel observations, so any differences in the flux densities represents either a change in the flux calibration or in the methods used to produce the catalogues.
We first consider the SPIRE results. The first difference the reader can notice is in the instrumental noise of the SPIRE maps: the values measured in the SDP release were 4.1, 4.7 and 5.7 mJy beam −1 at 250, 350 and 500 μm, which are different, in particular at 350 μm, with respect to the values in Tables 2 and 3. The differences are due to two effects. The first effect is caused by the use of the matched filter instead of the PSF to smooth the maps: as explained in Section 3.5, the matched filter is optimal to maximize the signal to noise in presence of both instrumental and confusion noise. This means that it might not be the optimal filter to reduce the instrumental noise alone, in fact the released filtered maps at 250 and 500 μm have a slightly larger instrumental noise than the SDP ones had. The second effect is caused by the difference in the choice of pixel sizes: the values used in the SDP maps were 5, 5 and 10 arcsec, while the released maps have pixel sizes of 6, 8 and 12 arcsec at 250, 350 and 500 μm, respectively. A larger pixel size increases the number of samples per pixel, decreasing the instrumental noise (and vice-versa). This is particular evident at 350 μm, where the pixel area is now about 2.5 larger than the SDP one.
The second difference is in the confusion noise of the SPIRE maps. SDP estimates of the confusion noise were 5.3, 6.4 and 6.7 mJy beam −1 at 250, 350 and 500 μm. In this release, we have used equation (14) to calculate a 'customized' confusion noise for each flux density and argued that the values in Tables 2 and 3 can be considered upper and lower limits, respectively (see Section 3.4).
The SDP values lie between the two limits, but of course any particular estimate of the noise for a single source will be different with respect to the previous release.
In this data release, we have chosen to include in the catalogue sources detected above a detection limit of 4σ , which corresponds in two-scan regions to flux-density limits of 29.6, 37.6 and 40.8 mJy at 250, 350 and 500 μm, respectively. Table 7 shows that the survey is 90 per cent complete at these limits. Our flux-density limits are actually very similar to the flux-density limits for our SDP catalogue: 34, 38 and 44 mJy at 250, 350 and 500 μm, respectively (Rigby et al. 2011). The SDP limits were 5σ limits, and the difference in the signal to noise is because of the more rigorous analysis of noise we have carried out for the current release.
We first compared our SDP catalogue with a new catalogue produced from the images smoothed with the PSF (not included in the data release), since this comparison allowed us to look for any differences that are not caused by the introduction of the matched filter. We found the median ratio of SDP flux density to flux density in the new catalogue is 1.02, 1.07 and 1.06 at 250, 350 and 500 μm, respectively. The changes are well within the error of 15 per cent we gave for our SDP catalogues (Rigby et al. 2011). These changes are not due to changes in the overall flux calibration, since there was no change in the flux calibration at 250 and 500 μm and change of <1 per cent at 350 μm (Section 3.1).
Apart from the introduction of the matched filter, the three significant modifications we have made to our method since the SDP are (a) a different way of removing the emission from the Galactic dust, (b) the use of a measured PSF rather than a Gaussian function (Section 3.2; Pascale et al. 2011) and (c) the new method of sequentially removing sources from an image once the flux density of the source has been measured (Section 5.1). We cannot be sure of the reason for the small difference in the flux densities, but the third modification should reduce the flux densities of sources which are close to other bright sources and thus works in the right direction.
If we now compare the matched-filter catalogues to the SDP catalogue, the median ratio of SDP flux density to flux density in the new catalogue is 1.06, 1.12, and 1.05 at 250, 350 and 500 μm, respectively. The fact that the SDP fluxes are on average slightly higher may represent the improvement given by the matched-filter technique in removing the contribution to the measured flux density of a source from a confusing nearby source. Figure 26. Flux density measured by PACS at 100 µm versus flux density measured by IRAS at 100 µm for the galaxies in the GAMA fields that are also in the IRAS Faint Source Catalogue (Wang et al. 2014). The line shows the median value of the ratio of PACS flux density to IRAS flux density for the 75 galaxies with IRAS flux density >1 Jy.
A last difference to notice about SPIRE measurements between SDP and this data release is about flux bias corrections (see Table 6). The differences are mainly due to the adopted number counts of the source populations. The model used by Rigby et al. (2011) has been superseded at faint fluxes by the deepest Herschel surveys, so we can assume our corrections are more robust and reliable than the previous ones.
We now consider the differences in the PACS flux densities. We selected the sources in the SDP catalogue detected at >5σ at one of the PACS wavelengths and looked for a source in the new catalogue within 5 arcsec of each SDP source. We found the median ratio of SDP flux density to flux density in the new catalogue is 0.96 at both 100 and 160 μm. The small change is reassuring, given the quite large changes we have made in the analysis of the data set and the improvements in the calibration of PACS since the publication of our SDP results.

Comparison between PACS 100 μm and IRAS
We have compared our measurements of the PACS 100 μm flux densities with the 100 μm flux densities for galaxies in the IRAS Faint Source Catalogue (Wang et al. 2014). There are 184 galaxies for which there are both PACS and IRAS flux densities. Fig. 26 shows a comparison of the PACS and IRAS flux densities. For galaxies with IRAS flux densities less than 1 Jy, there is almost no correlation between the IRAS and PACS flux densities, which we suspect is due to the effect of cirrus emission and flux bias on the IRAS flux densities. However, for galaxies above this limit there is, with a few exceptions, good agreement between the two sets of flux densities. Excluding the three very discrepant points, the median ratio of the PACS flux density to the IRAS flux density for the 75 galaxies with IRAS flux density greater than 1 Jy is 1.09.
This result is without making a correction for the real SEDs, which are unlikely to be the standard SED assumed for both PACS and IRAS: F ν ∝ ν −1 . The sources detected by IRAS are mostly nearby galaxies, and therefore a reasonable assumption is that the real SED of a source is likely to be a modified blackbody with a dust temperature of T d 20 K. For this SED, the correction to the PACS 100 μm flux densities is a reduction by a factor of 1.15 (Müller, Okumura & Klaas 2011). We have been unable to find any corrections listed for the IRAS 100 μm filter for modified blackbodies. No colour corrections are listed in the IRAS Explanatory Supplement for the 100 μm filter for black bodies with dust temperatures less than 40 K, presumably because the long-wavelength side of the IRAS 100 μm filter response function was not known very accurately (IRAS Explanatory Supplement). Given these uncertainties, the systematic 9 per cent difference between the PACS and IRAS flux densities seems well within the bounds of experimental uncertainty.

0 T H E C ATA L O G U E
In the data release, we include a catalogue of all sources above the detection limit of 2.5σ at 250 μm and with a measured flux density ≥4σ at least one of the three SPIRE wavelengths (250, 350 and 500 μm). The 4σ limit is approximately 29.4, 37.4 and 40.6 mJy at 250, 350 and 500 μm, respectively, and the catalogue is 90 per cent complete at all three wavelengths. These flux-density limits are very similar to the 5σ limits for the catalogue we released at the end of the SDP (Rigby et al. 2011), the reason for the difference being the more accurate analysis of noise for the present data release. The catalogue contains 113 995, 46 209 and 11 011 sources detected at >4σ at 250, 350 and 500 μm, respectively.
The catalogue contains measurements of the flux densities of each source in all five photometric bands, including aperture photometry for sources for which the diameter of the optical counterpart implies the 250 μm source is likely to be extended and the flux from aperture photometry is significantly larger than the flux measured by the source-detection algorithm (see Section 5.2). The catalogues contain detections at >3σ at 100 and 160 μm for 4650 and 5685 sources, respectively, and the typical noise at the two wavelengths is 44 and 49 mJy, respectively (see Section 5.3).
Given the area of the survey and the noise, Gaussian statistics imply that 260 of the sources in our 4σ catalogue should be spurious, which is 0.2 per cent of the total. However, this is likely to be an overestimate. The method we used to calculate the value of σ for each source was designed to provide a good estimate of the error in the flux density of the source. The more commonly used method of deriving a value of σ for a signal-to-noise calculation is to fit a Gaussian to the negative part of the pixel flux distribution (Section 3.4, Fig. 4). This method produces a much lower value of σ than we have actually used, and therefore our signal-to-noise estimates are likely to be conservative and our catalogue more reliable than the Gaussian statistics imply. In practice, we suspect a more important problem than spurious sources is likely to be single sources that are actually multiple sources, since our estimates of the fraction of single sources that are likely to actually be more than one source depends sensitively on the assumptions made about the correlation function for submillimetre sources, although this seems unlikely to be an issue at 250 μm where the sources appear to be only very weakly correlated (see Section 8.3, Fig. 25).
The catalogue also contains information about the optical counterparts to the sources, with fluxes and redshifts measured by the GAMA survey. How these counterparts were found is described in the second data-release paper (Paper II).

A P P E N D I X A : T H E DATA P RO D U C T S
In this section, we describe the basic data products we are releasing to the astronomical community. These consist of SPIRE and PACS maps of the GAMA fields and catalogues containing the sources detected at >4σ in at least one of the three SPIRE wavebands. The data products can be obtained from the website h-atlas.org. In this section, we give an overview of the different data products, describe what data products are appropriate for different scientific projects, and discuss some of the limitations of the data products. More detailed technical information can be found on the website or in the rest of the paper.

A1 The SPIRE maps
We provide three classes of image which are as follows.
(i) The raw images produced using the method of Section 3.1 but without the subtraction of any large-scale background emission, for example from interstellar dust ('cirrus emission'). We also supply maps of the instrumental noise on the images obtained using the coverage map calibrated by the jackknife method (Section 3.3). These images are our most basic data product and should be a good representation of the submillimetre sky up to an angular scale of 20 arcmin (Section 3.1), but these images should not be used without further modelling to investigate the emission on larger scales, for example the extragalactic background radiation.
(ii) The raw images with the large-scale background subtracted using Nebuliser (Section 3.1).
(iii) The images from (ii) convolved with the matched-filter (Section 3.5). These images were created to give the highest possible signal to noise for point sources, and the value in each pixel of the images is our best estimate of the flux density of a point source at that position. We also provide a noise map, which contains our best estimate of the instrumental noise at each position. equation (14) gives our estimate of the confusion noise as a function of the flux density of the source: this value should be added in quadrature to the instrumental noise (Section 5.1).
All of these images have different scientific purposes. For the reader interested in aperture photometry of a nearby galaxy, the appropriate images to use are the second set. It is possible to use the first set but we recommend the second set because any background emission has been removed with Nebuliser -Section 5.2. Photometry performed using apertures larger than 1.5 arcmin on images (ii) might be slightly underestimated due to the use of Nebuliser. The images we supply are in the units often called 'Jy beam −1 ', which means that the value in any pixel is the flux density that a point source would have if it were centred in that pixel. These images can easily be turned into a form suitable for aperture photometry by dividing the images by the SPIRE beam area. The current estimates of the beam area are 469, 831 and 1804 arcsec 2 at 250, 350 and 500 μm, respectively, but see the SPIRE Handbook for updates. We have calculated the errors on the flux densities by carrying out Monte Carlo simulations in the vicinity of the source. The reader can either use the same technique or use the following formula which describe quite well the results of our Monte Carlo simulations: in which σ inst,i is the instrumental noise in the ith pixel within the aperture, C conv is the area of the beam divided by the area of a pixel, N ap and N beam are the angular areas in pixels subtended by the aperture and the beam, and σ conf is the confusion noise given in Table 2. To calculate the number of pixels in the beam, it is recommended to use the area of the beam given by π (FWHM/2) 2 rather than the areas given in the SPIRE handbook, since the latter are estimated by integrating the PSF to very large radii and are thus likely to overestimate the effect of confusion. All flux densities measured by aperture photometry must be corrected or the fraction of the PSF that falls outside the aperture using the provided SPIRE aperture correction table.
If the reader is interested in running a different source-extraction program to MADX (Section 5.1), the second set of images is also the one to use. These images have had the large-scale emission subtracted using a method that does not affect the flux density of point sources.
For those interested in simply measuring the flux density for an object that should be unresolved in the SPIRE bands, the third set of images are the ones to use, since these have been designed to provide flux-density measurements for a point source of the highest possible signal to noise (Section 3.5).
If the reader is interested in carrying out a statistical 'stacking analysis', the set of images from (iii) is also the one to use. In this case, the reader should be aware that the mean of the images is not zero (Section 3.1), and so the mean of the map should be subtracted before performing any stacking analysis.
Note that if the reader is carrying out a stacking analysis or even simply measuring from the H-ATLAS image the submillimetre flux density of a previously known object, there is no need to make a correction for flux bias (Section 8.1), because flux bias only affects the flux densities in our catalogue.
The width of the SPIRE filters means that both the size of the PSF and the power detected by SPIRE depend on the SED of the source. The SPIRE data-reduction pipeline is based on the assumption that the flux density of the source depends on frequency −1 , and all our images are ultimately based on this assumption. The recipe for aperture photometry outlined above is also based on this assumption. This SED is very different from that of most sources detected by Herschel, and so the user must make a correction to the measured flux densities to allow for this difference. The reader should multiply the measured flux densities by the K ColE parameter, which is given in tables 5.6 and 5.7 in the SPIRE Handbook, which is a correction for all the effects produced by the difference in the SED.
The noise maps we have provided do not include a contribution from the uncertainty in the basic calibration of flux density. The calibration of Herschel is still improving. At the time of writing, the error in the flux density arising from the uncertainty in the absolute flux density of Neptune is 4 per cent and there is an additional 1.5 per cent error that is uncorrelated between the bands (SPIRE Handbook). The current recommendation (SPIRE Handbook) is that these factors should be added, and so the reader should use a calibration error of 5.5 per cent.

A2 The PACS images
We provide a single set of images. These are the images made with JScanamorphos which have had any residual large-scale emission subtracted with Nebuliser (Section 4). The units of these images are Jy pixel −1 , so it is straightforward to carry out aperture photometry. Because of the complicated nature of the noise on the PACS images (Section 4), we do not provide noise images. Instead, we provide maps showing the number of observations, N, contributing to each pixel in the final image. We have used Monte Carlo simulations (Section 5.3) to show that the error for aperture photometry, σ ap , on the PACS images through an aperture with a radius r can be represented by equation (9).
The PACS PSF is not a simple Gaussian and in fast-scan parallel mode is significantly extended in the scan direction (Lutz 2015), which means that it must vary even within a single GAMA field. For this reason, we recommend that no attempt should be made to maximize the signal to noise for point sources by convolving the images with the PSF. Instead, we recommend that all scientific projects should use aperture photometry. If the reader prefers to filter the map, we recommend the use of our Gaussian fit of the empirical PSFs, which gives an FWHM of 11.4 and 13.7 arcsec at 100 and 160 μm, respectively.
We recommend that aperture photometry of sources that are expected to be extended, for example nearby galaxies, should be carried out by adding up the flux density in a suitable aperture; there should be no need to estimate a sky value because we have already subtracted any residual background emission using Nebulizer (see Section 4.1). Errors in the flux density should be estimated using equation (9). As part of the data release, we have supplied a file listing the EEF in the two bands out to a reference radius of 1000 arcsec. Both the flux densities and the errors should be corrected using the EEF. Note that the flux densities of sources with a diameter larger than 2.5 arcmin may be underestimated because of the use of Nebuliser.
For sources that are expected to be unresolved, we have shown that the maximum signal to noise can be obtained by aperture photometry using very small apertures (Section 5.3, Fig. 10), with a diameter <8 arcsec at both 100 and 160 μm. We recommend that if the position of an object is known precisely, the reader should use a small aperture to carry out the photometry and then use the EEF to correct the flux density and the error to the reference radius. However, when using a small aperture that only contains a small number of pixels, the reader should think carefully about pixelization effects. This is also our recommended procedure for measuring the average flux density for a class of objects in a stacking analysis, as long as the positions of the objects are known very accurately.
If carrying out a statistical 'stacking analysis', the reader should be aware that the mean of the images is not zero (Section 4.1), and so the mean of the map should be subtracted before proceeding.
On top of the flux density error given in equation (9), there is also a fundamental calibration error. Our conservative estimate of this error is that it is 7 per cent (see Section 5.3).
All our measurements of flux density are based on the assumption that the flux density, F ν , of a source has the spectral dependence F ν ∝ ν −1 . Because of the width of the PACS spectral response and because most sources do not have this SED, a correction must be made to the flux densities. A table of corrections for different SEDs is given in Müller et al. (2011;Tables 1 and 2). We recommend that the user adopts a common-sense procedure here. To give one easy example, if the class of sources of interest are all nearby galaxies, we suggest that a sensible correction would be obtained by assuming that sources have SEDs that are modified black bodies with a dust temperature chosen by the reader. Note that the corrections can be quite large e.g. 1.16 for a modified blackbody with a dust temperature of 15 K (Müller at al. 2011). In the case of nearby galaxies the answer is fairly obvious. In other more tricky examples, it should be possible to use the ratios of the flux densities in the PACS and SPIRE bands or the redshifts of the optical counterparts to choose what correction to make.

A3 The catalogues
We provide catalogues of the sources detected above the detection limit of 2.5σ at 250 μm and at >4σ in at least one of the SPIRE bands. These catalogues also contain the optical counterparts to the Herschel sources. The details of the part of the data release that consists of measurements by telescopes other than Herschel are given in Paper II.
For each detected source, we provide measurements at 100, 160, 250, 350 and 500 μm bands of the flux densities and errors on the assumption that the source is unresolved by the telescope (Sections 5.1 and 5.3). For the sources that are expected to be extended, we also provide aperture photometry designed to estimate close to total flux densities for the sources (Sections 5.2 and 5.3). The errors given in the catalogues do not contain any calibration errors, which are described above (Appendix A1 and A2).
There are two significant corrections that we have not made to the catalogued flux densities that the user should be aware of. First, both the SPIRE and PACS data-reduction pipelines are based on the assumption that all sources have an SED in which F ν ∝ ν −1 . If the flux-density measurement in the catalogue was made with aperture photometry, which is the case for all PACS measurements and some SPIRE measurements, and the reader has some knowledge of the true SED of the source, it is possible to correct the flux density using the prescriptions given in Appendix A1 and A2. If the fluxdensity was measured with the MADX detection software there is no simple prescription for correcting the flux density. However, the systematic error on the flux density arising from a different SED is generally much less than the statistical error, and is thus not usually a concern. For any reader interested in precision photometry and in thus calculating a correction to the MADX flux densities for a different SED, we note that there are two effects. First, the precise PSF of a source depends on its SED, whereas in MADX we have used a PSF derived from observations of Neptune. Secondly, there is a correction that needs to be made to the flux density arising from the variation in the SED across the SPIRE filter. The second effect at least is relatively easy to calculate using the information in the SPIRE Handbook.
The second correction we have not made to the measured flux densities is for the effect of flux bias. Although flux bias does not affect Herschel flux-density measurements of objects detected by other telescopes (Appendix A1), it does affect the flux-density measurement of the sources in our catalogues. Despite the extensive modelling described in Section 8, we have not made corrections to the flux densities in the catalogue because we are aware that there may be improvements in the modelling of this effect in the future. We suggest that the user either uses our estimate of this effect, which are listed in Table 6, or in their scientific analysis make allowance for a systematic error in the SPIRE flux density approximately of the size that we estimate.

A P P E N D I X B : A S T E RO I D S
We found the asteroids by looking at the jackknife images, the images made by taking the difference of the two images of the same region made with the different scan directions (Section 2). The only sources on the jackknife images should be objects that either vary in flux density or move. In practice, the only time-varying sources we found were asteroids. The list of source positions is given in Table A1. In the table, we also list the times and dates at which