The Herschel–ATLAS Data Release 2, Paper I. Submillimeter and Far-infrared Images of the South and North Galactic Poles: The Largest Herschel Survey of the Extragalactic Sky

The SPIRE instrument consists of three imaging arrays for observations at 250, 350, and 500 μm with 139, 88, and 43 bolometers, respectively. Each array has two associated thermistors to monitor the temperature of the array, although after launch only one of the 350 μm thermistors worked, and two dark bolometers that receive no light. In Herschel parlance, the "level-0" data are the raw instrumental/telescope data and the "level-1" data consist of the calibrated flux-density versus time measurements for the individual bolometers (calibrated timelines), which can then be used to create an image of the sky. In this section we explain the way we produced the calibrated level-1 data.

We converted the level-0 data to the level-1 data, the calibrated timelines, using the Herschel Interactive Processing Environment (HIPE, Ott 2010), version 11.0 1200 (development build). Unless described otherwise, we used the standard components of the data-reduction pipeline. Forty two of 51 observations for the NGP and SGP fell in observing days 320–761, during which there were positional offsets caused by a change in the operating temperature of the star-tracker camera for which the camera was not re-calibrated. We used an updated pointing model released by the Herschel Science Centre in 2012 to correct the pointing for these observations.

We corrected glitches in the bolometer and thermistor data using a different technique from the standard module in the pipeline. Instead of the default wavelet deglitcher, we used the sigma–kappa deglitcher, since tests on both our H-ATLAS data and on data from the Herschel Virgo Infrared Cluster Survey (HeViCS, Davies et al. 2010) showed it performed better for parallel-mode observations with high scan-speeds and a reduced (10 Hz) sampling rate (the non-default settings for the sigma–kappa deglitcher task were: kappa = 3.5, gamma = 0.1, boxFilterCascade = 5, largeGlitchRemovalTimeConstant = 7, iterationNumber = 3).

At this point, we had calibrated timeline data (level-1) but the data were still affected by artefacts, including "jumps" and gradual changes in the signal caused by changes in the temperature of the bolometers. The thermistor timeline data contained the necessary information to correct for the effect of temperature. Since the 350 μm array only had one working thermistor, we used dark-bolometer 1 as a replacement. However, before we corrected the bolometer timelines, it was necessary to correct both the bolometer and thermistor timelines for the jumps.

A jump is an instantaneous change in the voltage of an individual bolometer or thermistor (a typical example is shown in Figure 6). In rare cases, rather than a step change in voltage, there is a sudden large change followed by a gradual decay back toward the original value. Jumps appear to be more common in thermistors than bolometers. Our jump correction method (see below) does not work well in correcting these rare jumps, and if one of these occurred in a thermistor, we used one of the dark bolometers as a replacement when correcting the bolometer timelines for the effect of temperature. We looked for jumps in bolometers and thermistors in different ways.

Zoom In Zoom Out Reset image size

For the thermistors, rather than using the automatic jump detector, we inspected both thermistor timelines for each array by eye to spot jumps in the timelines, using the Kst visualization tool.⁸ In the case of the 350 μm array, we carried out a similar inspection of the timelines of dark-bolometer 1.

It was not practical to search for jumps in the bolometer timelines in the same way because there were too many of them. Instead we made initial maps of each individual observation from the timelines and visually searched for the light and dark thin streaks caused by jumps; since a single H-ATLAS observation consists of scans in a single direction and one map pixel usually only includes data from a few bolometers, the effects of a jump are easy to see.

Before correcting the jumps, we combined all scan legs from an observation, including the "turn-around" data, into a single timeline. We then corrected all the jumps in the timeline by fitting a linear relationship to portions of the timeline immediately before and after the jump, and then adding the difference in these relationships to the timeline after the jump (see Smith 2012 for more details). We replaced the samples immediately around the jump (Figure 6) with random noise, and these samples were then masked and not used to make the final maps.

The advantage of combining the data from all scan legs into a single timeline is that makes it possible to remove more accurately the drift in the bolometer signals caused by temperature changes. In the standard pipeline, this correction is done separately for each scan leg and the information in the "turn-around" data is not used at all. Before we made the correction for the effect of temperature, we masked any samples in the timelines that had been flagged as bad (e.g., samples effected by glitches, samples in which the signal is saturated) and any places in the timelines where there were obvious bright sources.

We fitted the following relationship between the signal measured by a single thermistor (S_T) and the signal measured by the ith bolometer (${S}_{\mathrm{bolom},i}$): ${S}_{\mathrm{bolom},i}=a\times {S}_{T}+c$. We then subtracted this relationship from the bolometer signal, effectively removing the effect of the temperature change. A difference from our procedure in the GAMA fields (V16) is that we carried out this fit for both thermistors, and then used the thermistor that produced the best fit to the data to correct the bolometer timelines. For the 350 μm timelines, we did the fits for the one working thermistor and for dark-bolometer 1. For any observation for which one of the thermistors was saturated or affected by a jump that could not be accurately corrected (see above), we used one of the dark bolometers rather than the thermistor. There were parts of some timelines where the linear relationship given above did not provide a good fit to the data. These almost always occurred six hours after a cycle of the SPIRE cooling system and became known as "cooler burps." In cooler burp regions the timelines vary far more more rapidly than the average SPIRE timeline. For these timelines we fitted a fifth-order polynomial rather than a linear relationship.

Once the thermal drift correction had been applied, we applied a high-pass filter to remove any residual drifts. Before applying the filter, we removed the brightest sources from the timelines, and then restored these after the filtering. The high-pass filter corresponds to a scale on the sky of 42, which was chosen to minimize the 1/f noise on the images (Pascale et al. 2011). Our images will therefore not contain any structures on scales larger than this. In practice, however, one of the effects of our scanning method is that any structure on scales >20' is attenuated (Waskett et al. 2007; see Section 5.4 for more details).

The next step in the data reduction was to make maps from each individual observation. The two purposes of these initial maps were to check the astrometry of each observation and to remove more bad data, since some low-level artefacts and samples containing bad data are easier to find on the maps than in the timelines.

We carried out our astrometric calibration of the observations in the NGP using the technique that is described more fully in Smith et al. (2011). Briefly, we produced initial source catalogs for each map using our source-detection method (V16). We then produced histograms of the differences in R.A. and decl. between the positions of the sources and the positions of all objects on the SDSS DR7 r-band images (Abazajian et al. 2009) within 50'' of each source. We then fitted these distributions using a Gaussian model for the SPIRE positional errors, allowing for the effect of clustering in the SDSS data (see Smith et al. 2011 for details). This procedure allowed us to measure the average difference in positions in both R.A. and decl. 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 ranged from less than an arcsec to a 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 astromometry for each Herschel observation, so that the effective calibration of the NGP maps and catalogs should be the same as the SDSS.

The SGP field is not covered by the SDSS. To calibrate the astrometry of the SGP observations, we used the same method as above but we replaced the SDSS catalogue with the catalogue from the VLT Survey Telescope ATLAS (Shanks et al. 2015). We found astrometric shifts of a similar size to those for the NGP except a few tiles with corrections of ∼4 arcsec.

We also used these initial single-observation maps to look for any residual artefacts. In particular, the standard deglitcher modules (Section 3.1) do not completely mask very large glitches, due to a parameter which limits the maximum number of samples that can be masked. These unmasked samples produce "glitch tails" on the images, linear features in the scan direction. We looked for these by eye on the map and then masked the appropriate parts of the timelines.

We also used a new iterative technique to look for glitches that were too faint to be detected by algorithms that work on the timelines. This technique was not applied to the data for the GAMA fields, since we developed it after our first data release. In this technique, we look for bolometer samples that are discrepant by at least 5σ from the value predicted from the statistics of the map. The first step in the procedure is to use all the observations to make a low-resolution map of each field with three times the default pixel size (Section 3.3). The flux in each pixel, ${F}_{\mathrm{MapPixel}}$, is the mean of the fluxes of the ${N}_{\mathrm{Pixel}}$ timeline samples that fall within that pixel. If ${F}_{\mathrm{Sample}}$ is the flux density of a single bolometer sample and ${\sigma }_{\mathrm{MapPixel}}$ is the error in ${F}_{\mathrm{MapPixel}}$, we treat a sample as bad, and thus mask it, if

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{\mathrm{Sample}}-{F}_{\mathrm{Map}\mathrm{Pixel}}}{{\sigma }_{\mathrm{Map}\mathrm{Pixel}}\sqrt{{N}_{\mathrm{Pixel}}}}\gt 5.\end{eqnarray} \tag{ 1 }$$

After masking all the discrepant samples in the timelines, we remade the maps and looked for additional samples that met the criterion in Equation (1). We carried out four iterations of this procedure, masking in total 295,496 and 539,126 bolometer samples for the NGP and SGP, respectively, which equates to ∼0.02% of all bolometer samples.

As the final step in the processing of the timelines, we masked the "turn-around" data (i.e., regions where the satellite was not scanning at a constant speed) at the end of each scan leg.

We created the final maps by combining all the corrected and masked level-1 data for each field. We used the simple ("naive") map-maker in which the flux density in each pixel is taken to be the mean flux density of all the bolometer samples that contribute to that pixel. Despite the concern before launch (Section 2; Waskett et al. 2007) that sophisticated map-making algorithms would be necessary to remove large-scale artefacts, the use of the thermistors to correct for the thermal drift of the arrays worked well enough that this simplest of all map-making techniques was sufficient. As for the GAMA fields (V16), we used a pixel size of 6'', 8'', 12'' for the 250, 350, and 500 μm bands, respectively, which are different from the default pixel sizes of Herschel images; we chose them because they correspond to roughly one third of the size of the point-spread function (PSF) (full-width half maximum; FWHM) in each band (see below) and they are big enough that the chance of a map pixel containing no bolometer samples, thus producing a Not-a-Number pixel, is low. The 250 μm maps for both fields are shown in Figures 2 and 4.

The standard pipeline produces an estimate of the uncertainty in the flux density measured in each pixel by calculating the variance of all the timeline samples that contribute to that pixel. However, this method does not produce an accurate estimate of the flux uncertainty for two reasons: (a) the small number of samples in each pixel means that the error on the uncertainty estimate is quite large; (b) the variance will be too high if the pixel coincides with a bright object. Instead we have produced our own uncertainty maps, using ${\sigma }_{\mathrm{inst}}/\sqrt{{N}_{\mathrm{sample}}}$ as our estimate of the uncertainty in flux density, in which ${\sigma }_{\mathrm{inst}}$ is the instrumental noise for one timeline sample, calculated using the method described in Section 5.1, and ${N}_{\mathrm{sample}}$ is the number of timeline samples contributing to that pixel. This, of course, is an estimate of the uncertainty in flux density arising from instrumental noise and does not include the effect of source confusion (Section 5.2).

For our PSF, we use the same PSF that was determined by V16 from images of Neptune with the same pixel size as the H-ATLAS images (see that paper for more details). The FWHM of the azimuthally averaged PSF is 17.8, 24.0, and 35.2 arcsec at 250, 350, and 500 μm, respectively.

As part of the data release, we also produced images optimized for the detection of point sources. The first step in producing these images was to remove any large-scale structure from the images, which is mostly emission from Galactic dust ("cirrus" emission). We removed the cirrus emission using the Nebuliser⁹ algorithm developed by the Cambridge Astronomical Survey Unit and we refer the interested reader to V16 for the details of how we did this. The result of the application of Nebuliser is that the images should not contain any emission, whether from Galactic dust or from outside the Galaxy, with an angular scale $\succcurlyeq 3$ arcmin.

For an image containing only one point source and instrumental noise, the maximum signal-to-noise for the source is obtained by convolving the image with the PSF (North 1943). However, the noise in the H-ATLAS images is a combination of instrumental noise and "confusion noise," the result of the large number of submillimeter sources that are too faint to detect individually but which merge together to form an undulating background to the images. Chapin et al. (2011) have shown how to calculate a convolving function or "matched filter" that will produce the maximum signal-to-noise for an unresolved source for any ratio of confusion to instrumental noise. As part of this data release, we have produced images optimized for finding point sources by convolving the raw images with the matched filters we used for the GAMA fields (V16). The instrumental noise and confusion noise for the SGP and NGP are actually slightly different than for the GAMA fields (see Sections 5.1 and 5.2), which means that the matched filters we have used are not precisely optimized, but this small disadvantage is outweighed by our being able to use the results of the extensive simulations we carried out with the GAMA matched filters (V16). We tested the effect of using a matched filter optimized with our updated noise values and found that the difference in the number and fluxes of sources is negligible. For readers interested in measuring the flux density of a point source, the matched-filtered images are the ones to use.

The flux calibration we applied to the images was publicly released as calTree v11. There is no difference between this and the most recent (at the time of writing) flux calibration (calTree v14). Because of the change in the SPIRE calibration, the flux densities for the NGP and SGP are not quite on the same scale as those in the GAMA fields (V16). To create the maps described by V16, we used the SPIRE v5 calibration tree to create the level-1 data, but applied a 1.0067 correction factor to the 350 μm data to make the effective calibration the same as calTree v8. Between v8 and v14 the average (multiplicative) change in flux density is 1.0253 ± 0.0012, 1.0182 ± 0.0045, and 1.0125 ± 0.0006 at 250, 350, and 500 μm, respectively. Therefore, to put the H-ATLAS GAMA flux densities on the same scale as those for the NGP and SGP, we need to multiply the flux densities from the DR1 release by these factors. Note, however, that because each bolometer is calibrated individually, the actual correction factor for an individual source depends on which individual bolometers crossed that position.

In determining the instrumental noise, the first step is to remove any real astronomical signal (e.g., galaxies, cirrus, confusion noise) by creating a jackknife map from subtracting two images of the same part of the sky made from individual SPIRE observations. The instrumental noise can then be measured from the jackknife map. For all pixels in both the NGP and SGP that are covered by at least two individual observations, we calculated the instrumental noise per single bolometer sample from:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{Sample}}=\sqrt{\displaystyle \frac{\sum _{i}^{\mathrm{Npix}}{({M}_{\mathrm{ortho},i}-{M}_{\mathrm{nom},i})}^{2}}{\sum _{i}^{\mathrm{Npix}}\left(\tfrac{1}{{C}_{\mathrm{ortho},i}}+\tfrac{1}{{C}_{\mathrm{nom},i}}\right)}}\end{eqnarray} \tag{ 2 }$$

in which ${M}_{\mathrm{ortho},i}$ and ${M}_{\mathrm{nom},i}$ are the flux densities in the ith pixel in the two maps out of which the jackknife is made (the flux in the jackknife map is ${M}_{\mathrm{ortho},i}\mbox{--}{M}_{\mathrm{nom},i}$) and ${C}_{\mathrm{ortho},i}$ and ${C}_{\mathrm{nom},i}$ are the numbers of timeline samples contributing to the ith pixel in the two maps. We measured the uncertainty on ${\sigma }_{\mathrm{Sample}}$ by randomly assigning the pixels to four groups and calculating ${\sigma }_{\mathrm{Sample}}$ separately for each group. We repeated this five times, which was enough to give a reasonable estimate of the uncertainty in ${\sigma }_{\mathrm{Sample}}$.

We measured the noise per bolometer sample separately for the SGP and NGP fields and the values are given in Table 1. There is generally good agreement between the two fields: a difference of $\simeq 0.1$ mJy at 250 and 350 μm and $\simeq 0.6$ mJy at 500 μm. These differences are much greater than the measured uncertainty; we do not know the reason for this but it does not seem likely to have any practical consequences. The averages of the noise values for the two fields are 31.38, 32.08, and 36.21 mJy beam⁻¹ at 250, 350, and 500 μm, respectively. After adjusting for the small difference in average calibration (Section 3.4), these values are higher at 250 and 350 μm than those reported in V16 for the GAMA fields and lower at 500 μm, but the difference is for all bands is <2%.

Table 1.  SPIRE Instrumental Noise

		Raw Maps			Matched-filtered Maps
	Field	250 μm	350 μm	500 μm	250 μm	350 μm	500 μm
		(mJy beam−1)	(mJy beam−1)	(mJy beam−1)	(mJy beam−1)	(mJy beam−1)	(mJy beam−1)
Noise per	NGP	31.327 ± 0.005	32.001 ± 0.006	35.922 ± 0.012	⋯	⋯	⋯
Sample	SGP	31.426 ± 0.005	32.149 ± 0.006	36.506 ± 0.006	⋯	⋯	⋯
${N}_{\mathrm{scan}}=1$	NGP	19.3715 ± 0.0031	19.2201 ± 0.0036	19.9354 ± 0.0066	10.1834 ± 0.0018	10.0136 ± 0.0021	15.8384 ± 0.0060
	SGP	17.9236 ± 0.0030	17.8210 ± 0.0032	19.4017 ± 0.0034	18.2855 ± 0.0029	18.5682 ± 0.0031	35.3218 ± 0.0093
${N}_{\mathrm{scan}}=2$	NGP	10.3446 ± 0.0017	10.1512 ± 0.0019	11.1793 ± 0.0037	5.0488 ± 0.0008	4.9319 ± 0.0010	5.6003 ± 0.0021
	SGP	10.7022 ± 0.0018	10.5098 ± 0.0019	11.7208 ± 0.0020	5.2081 ± 0.0008	5.1019 ± 0.0008	5.8404 ± 0.0015
${N}_{\mathrm{scan}}=3$	NGP	8.6163 ± 0.0014	8.4736 ± 0.0016	9.3358 ± 0.0031	4.2170 ± 0.0007	4.1447 ± 0.0009	4.7175 ± 0.0018
	SGP	8.9101 ± 0.0015	8.7633 ± 0.0016	9.7512 ± 0.0017	4.3557 ± 0.0007	4.2836 ± 0.0007	4.9041 ± 0.0013
${N}_{\mathrm{scan}}=4$	NGP	7.4953 ± 0.0012	7.3730 ± 0.0014	8.1258 ± 0.0027	3.6817 ± 0.0006	3.6268 ± 0.0007	4.1126 ± 0.0016
	SGP	7.7019 ± 0.0013	7.5780 ± 0.0014	8.4179 ± 0.0015	3.7694 ± 0.0006	3.7243 ± 0.0006	4.2594 ± 0.0011
${N}_{\mathrm{scan}}=5$	NGP	6.9083 ± 0.0011	6.8135 ± 0.0013	7.5298 ± 0.0025	3.4204 ± 0.0006	3.3723 ± 0.0007	3.8342 ± 0.0015
	SGP	6.9362 ± 0.0012	6.8602 ± 0.0012	7.6284 ± 0.0013	3.3952 ± 0.0005	3.3831 ± 0.0006	3.8603 ± 0.0010
${N}_{\mathrm{scan}}=6$	NGP	6.1648 ± 0.0010	6.0691 ± 0.0011	6.6935 ± 0.0022	3.0391 ± 0.0005	2.9969 ± 0.0006	3.4182 ± 0.0013
	SGP	6.4155 ± 0.0011	6.3283 ± 0.0011	7.0319 ± 0.0012	3.1369 ± 0.0005	3.1289 ± 0.0005	3.5077 ± 0.0009
${N}_{\mathrm{scan}}=7$	NGP	5.6386 ± 0.0009	5.5154 ± 0.0010	6.0978 ± 0.0020	2.7739 ± 0.0005	2.7573 ± 0.0006	3.2051 ± 0.0012

Note. The instrumental noise properties of the SPIRE maps. The top two rows show the instrumental noise per bolometer sample in the two fields. The other rows correspond to the average instrumental noise per map pixel for pixels with the same number of scans (${N}_{\mathrm{scan}}$) for the raw (left) and matched-filtered maps (right).

Download table as:  ASCIITypeset image

A common assumption is that the instrumental noise in Herschel images is Gaussian. We have tested this in Figure 7. We divided the pixels in the jackknife maps by the number of timeline samples each contains, only keeping pixels that have the same number of samples in both of the maps used to make the jackknife. The figure shows the noise per pixel plotted against the number of timeline samples contributing to the pixel (C_i). If the noise is Gaussian, we would expect the noise per pixel to decrease as ${C}_{i}^{-1/2}$. The dashed lines in the figure show the predictions of this Gaussian model, using our noise-per-sample estimates. The model agrees almost exactly with the observations, confirming that the instrumental noise does have Gaussian statistical properties.

Zoom In Zoom Out Reset image size

As part of the data release, we produced maps showing the noise per pixel in the two fields. As our estimate of the instrumental noise in each pixel, we have used ${\sigma }_{\mathrm{Sample}}/\sqrt{{C}_{i}}$, in which ${\sigma }_{\mathrm{Sample}}$ is given in the top panel of Table 1 and C_i is the number of timeline samples in each pixel. As a useful guide to the instrumental noise in parts of the raw and matched-filtered images made with different numbers of scans, in Table 1 (bottom panel) we have given the average instrumental noise for map pixels produced from data from ${N}_{\mathrm{scan}}$ individual observations, with values of ${N}_{\mathrm{scan}}$ from 1 to 7.

By far a more difficult quantity to define and measure is the confusion noise (as V16 says, "confusion is confusing"), and different scientific objectives require different methods for measuring it. Source confusion has several different effects on observations. Two of the most important are: (a) confusion increases the difficulty of detecting sources by increasing the overall noise on images; (b) confusion increases the error on the flux measurements. V16 used two different definitions of confusion noise: one suitable for measuring the signal-to-noise with which sources are detected on an image, and one suitable for estimating the errors in flux measurements. We have used the same two definitions of confusion, but improved the methods for measuring them described in V16.

The first method was designed to produce a confusion estimate suitable for measuring the noise in signal-to-noise estimates. V16 estimated this using the histogram of pixel values in the SPIRE maps (Figure 8). The shape of this distribution is produced by the instrumental and confusion noise in the maps and the significant individual sources, which produce the right-hand tail in the figure. On the assumption that the tail of individual sources is not relevant for estimating the noise in signal-to-noise estimates, V16 measured the confusion noise by first fitting a Gaussian to the negative part of the pixel histogram, thus avoiding the positive tail, and then calculating $\sqrt{{\sigma }_{\mathrm{tot}}^{2}-{\sigma }_{\mathrm{inst}}^{2}}$, in which ${\sigma }_{\mathrm{tot}}$ is the standard deviation of the best-fitting Gaussian and ${\sigma }_{\mathrm{inst}}$ is the standard deviation of the Gaussian distribution for the instrumental noise.

Zoom In Zoom Out Reset image size

Here we have used the same definition of confusion noise but a slightly different approach for measuring it. We used the coverage map and noise-per-sample measurements to generate an artificial image containing only instrumental noise. We then fit for the confusion noise (using the Python lmfit package) by adding the confusion noise to our artificial map assuming the confusion noise has a Gaussian distribution. A ${\chi }^{2}$ statistic is calculated from the difference between the real pixel histogram and that from the model image; the iteration with the lowest ${\chi }^{2}$ value gives our best estimate of the confusion level. To avoid the biasing effect of the positive tail produced by the significant sources, we generally only calculated ${\chi }^{2}$ for the negative side of the pixel histogram, with the bins to the right of the peak contributing only if the value for the model lay above the real distribution. We repeated this whole process 96 times to produce an estimate of the error on our measurement of the confusion noise.¹³ The biggest advantage of this method over that of V16 is that we are not assuming a single instrumental noise for the whole map.

Figure 8 shows the artificial histograms that produce the best fit to the histograms for the raw 250 μm image of the NGP, the 250 μm image from which the background has been subtracted using Nebuliser, and the 250 μm image from which the background has been subtracted and which has then been convolved with the matched filter. The values of the confusion noise that give the best fit to the data for both the SGP and NGP are given in Table 2 for all three wavebands. Our confusion noise values are the most accurate values produced at the SPIRE wavelengths, with both fields agreeing to within 0.1 mJy beam⁻¹ for the nebulised and matched-filtered maps in all bands (the exception is SPIRE 250 μm where the difference is 0.14 mJy beam⁻¹). This is a significant improvement over V16 whose estimates varied by up to 0.9 mJy beam⁻¹ between fields. Our confusion values tend to be slightly lower than those of V16 on the nebulised images and slightly higher on the matched-filtered images, which is probably due to the improved method we have used in this paper.

Although the confusion estimates for the NGP and SGP agree well, our uncertainty estimates are so small that the differences between the NGP and SGP are formally significant. Although these fields are very large, it is possible that these differences are due to large-scale interstellar cirrus or large-scale extragalactic structure (Negrello et al. 2017). Part of the explanation may be that the instrumental noise is such a large part of the total noise, especially in the raw maps, that small errors in the estimate of instrumental noise may lead to large errors in the estimate of confusion noise. Another possible problem may be that our assumption that the source population can be divided into a population of faint confusing sources, which produces a Gaussian pixel distribution, and a population of sources which are detected individually and which produce a non-Gaussian tail to the pixel distribution, may be too simplistic.

The second definition of confusion noise used by V16 was one designed to produce an estimate suitable for estimating the errors on flux measurements. Errors on flux measurements are produced by all the other sources on the image, not just the faint ones contributing to the Gaussian distribution in Figure 8 but also the tail of significant sources. The only sources on an image that cannot contribute to the flux error for a source are the pixels in the map that are brighter than that source. To produce an estimate of the confusion noise appropriate for a source with flux density F_s, V16 measured the variance on an image but only included pixels with flux densities $\lt {F}_{s}$. They then took the confusion noise as $\sqrt{{\sigma }_{\mathrm{var}}^{2}-{\sigma }_{\mathrm{inst}}^{2}}$, in which ${\sigma }_{\mathrm{var}}^{2}$ is the variance and ${\sigma }_{\mathrm{inst}}$ is the standard deviation of the Gaussian distribution for the instrumental noise (see V16 for additional details).

The only difference between our method and that of V16 is that we have taken account of the variation in depth over the image arising from different coverage levels, C_i. We have allowed for this by estimating the confusion noise separately for pixels with different numbers of bolometer samples and then averaging the different estimates of the confusion noise.

Figure 9 shows the results for the raw maps and for the maps from which the background has been subtracted using Nebuliser. As expected, the confusion noise for the images from which the background (mostly cirrus emission) has not been subtracted is higher than for the images from which the background has been subtracted. The results for the NGP and SGP are almost the same. By definition, the confusion noise depends on F_s. As in V16, we use a value of F_s of 200 mJy to estimate a confusion noise that is easy to compare with the estimates of others. For the background-subtracted maps and with F_s = 200 mJy, we find the values of the confusion noise are 6.62, 7.16, and 6.69 mJy beam⁻¹ at 250, 350 and 500 μm, respectively, for the NGP field, and 6.66, 7.15, and 6.73 mJy beam⁻¹ at 250, 350, and 500 μm, respectively, for the SGP field. V16 found mean values of 6.53, 7.03, and 6.58 mJy beam⁻¹ for the three GAMA fields at 250, 350, and 500 μm, respectively (corrected for the calibration differences), using a very similar method. As for the previous method, the differences between the NGP and the SGP are much smaller than those between the GAMA fields found by V16. Nguyen et al. (2010) estimated the confusion noise in the HERMES survey, using a fairly similar technique to this second method, and derived estimates of 5.8, 6.3, and 6.8 mJy beam⁻¹.¹⁴ The estimates in our work using this technique are broadly similar to those of V16 but systematically higher than those of Nguyen et al. (2010). It is unknown whether these differences are genuinely the result of differences in the source populations in the different fields, or due to the differences in the methods used in the different fields.

Zoom In Zoom Out Reset image size

By combining our histogram confusion estimates with those of the instrumental noise from Section 5.1, we find the total pixel noise in our nebulised images of 11.0, 11.1, and 12.3 mJy beam⁻¹ at 250, 350, and 500 μm, respectively. If the matched filter is used to extract point sources is applied to the map, our 1σ map sensitivity estimate is 5.7, 6.0, and 7.3 mJy beam⁻¹, respectively.

The correct method to use for photometry depends on whether the object is expected to be extended or unresolved by the SPIRE PSF. If the object is unresolved, the best method is to use the flux value at the object's position on the SPIRE image that has been convolved with the matched filter. As part of the data release, we produced a map of the instrumental noise on this image. However, the error on the flux density will also include a component from the confusion noise. The correct value to use for the confusion noise depends on the purpose of the investigation (see the previous section), but for detection experiments we suggest using the values obtained from fitting the Gaussian part of the pixel histograms, which are given in Table 2. For an estimate of the error of the flux density of an individual source, the correct confusion noise value to use would be that obtained from measuring the variance on the image, which depends on the flux density of the source, as shown in Figure 9, and given in Table 6 of Appendix B. Whichever version of the confusion noise is chosen, the confusion noise and the instrumental noise, taken from the map of instrumental noise, should be added in quadrature.

Astronomers interested in carrying out a statistical "stacking analysis," in order to measure the mean submillimeter flux density of some class of object, should use the images that have been nebulised to remove cirrus emission. They should be aware that the means of the maps are not necessarily zero and so they should subtract the mean from the map before carrying out the analysis.¹⁵ We recommend using the SIMSTACK algorithm (Viero et al. 2013) or similar, that allows the user to correct for the effects of clustering. We would recommend carrying out a Monte Carlo simulation in which they measure the mean flux density at a large number of random positions. This procedure will (a) give an estimate of the mean level on the map and (b) produce an empirical estimate of the error in the stacking measurements.

Photometry of extended sources should be carried out using aperture photometry. The images supplied in the data release have units of Jy beam⁻¹. These can be converted into images with units of Jy pixel⁻¹, suitable for aperture photometry, by dividing the flux value in each pixel by a factor ${C}_{\mathrm{conv}}$, which is given by the area of the telescope's beam divided by the area of a pixel. The current values in the SPIRE Data Reduction Guide¹⁶ are 469.4, 831.2, and 1804.3 arcsec² at 250, 350, and 500 μm, respectively. Note that it is possible to produce SPIRE maps that have been optimized for aperture photometry, using the SPIRE "relative gain" method. However, for simplicity, we decided to produce only a single set of maps, optimized for point-source detection.

We recommend carrying out aperture photometry on the images from which the background has been subtracted with Nebuliser. The application of Nebuliser does mean that the flux density of any sources with a size greater than 3' might be underestimated,but our tests on the GAMA fields found no evidence of this effect (V16). We also found that the photometric errors were smaller if we used the nebulised images. Note that in this case there is no need to subtract the mean map value, since the application of Nebuliser should already have subtracted the best estimate of the sky level at that position. Some of the object's emission will fall outside the aperture because of the extended profile of the PSF (Griffin et al. 2013). As part of the data release, we have supplied a table of corrections factors for this effect.

We carried out Monte Carlo simulations to estimate the errors in the flux densities measured with aperture photometry. We placed apertures randomly on the NGP and SGP maps in areas which are made from two individual observations (${N}_{\mathrm{scan}}=2$). The aperture radii were varied in size from approximately the beam size up to 100'', in 2'' intervals, and for each radius we used 3000 random positions. Figure 10 shows the results of the Monte Carlo simulation, and very consistent results between the two fields can be seen.

Zoom In Zoom Out Reset image size

We assumed the relationship between flux error and radius is a power law, since if the noise is dominated by instrumental random noise we should get a simple linear relationship: flux error ∝radius. We found that we needed to use two power laws to fully describe the relationship at all radii, with the change in relationship occurring at 50''. Our model is described by:

$$\begin{eqnarray}{\sigma }_{\mathrm{ap}}(\mathrm{mJy})=\left\{\begin{array}{ll}{{Ar}}^{\alpha } & \mathrm{if}\ r\leqslant 50^{\prime\prime} \\ B{(r-50)}^{\beta }+A{50}^{\alpha } & \mathrm{for}\ r\gt 50^{\prime\prime} \end{array}\right.\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{\mathrm{ap}}$ is the flux error in mJy and r is the radius in arcseconds. The best-fit values for this relationship for all bands and fields are given in Table 3. Above a radius of 50'' the relationship is similar to that expected for pure instrumental noise, with values of β between 0.98 and 1.17. Below 50'' the relationship is much steeper with values of α between 1.37 and 1.48. This may be due to small-scale cirrus emission, which would not have been removed with the filtering scale used in Nebuliser, or possibly some effect of source confusion.

Table 3.  Aperture Noise Model Best-fit Parameters

Waveband	Field	A	α	B	β	X
100 μm	NGP	0.749	1.475	6.244	0.971	⋯
	SGP	0.720	1.473	6.235	0.970
160 μm	NGP	0.642	1.444	4.193	0.995	⋯
	SGP	0.620	1.446	4.247	0.992
250 μm	NGP	0.152	1.388	0.336	1.179	$5.13\times {10}^{-4}$
	SGP	0.164	1.368	0.527	1.066
350 μm	NGP	0.117	1.389	0.539	1.016	$2.91\times {10}^{-4}$
	SGP	0.111	1.410	0.497	1.016
500 μm	NGP	0.052	1.464	0.372	1.033	$1.39\times {10}^{-4}$
	SGP	0.056	1.451	0.459	0.984

Note. The best-fit parameters for the relationship between flux error and aperture radius (Equation (3)). See Section 5.3 for details. The X column gives the constant required to correct the relationship for regions of the map with different ${N}_{\mathrm{scan}}$ values (see Equation (4)).

Download table as:  ASCIITypeset image

As the areas where ${N}_{\mathrm{scan}}\gt 2$ are limited in size, we are unable to perform a Monte Carlo simulation for the deeper regions. As we believe we understand the properties of the instrumental noise, we can account for the differences in ${N}_{\mathrm{scan}}$ by subtracting our instrumental noise for ${N}_{\mathrm{scan}}=2$ in quadrature and adding back in quadrature the appropriate noise (as measured in Section 5.1 and tabulated in Table 1). The relationship between the flux error and aperture radius, for any value of ${N}_{\mathrm{scan}}$, is then given by:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ap}}(\mathrm{mJy})\\ &&=\,\left\{\begin{array}{l}\sqrt{{({{Ar}}^{\alpha })}^{2}-X({\sigma }_{\mathrm{inst},2}^{2}-{\sigma }_{\mathrm{inst},{\rm{N}}}^{2}){r}^{2}}\\ \mathrm{if}\ r\leqslant 50^{\prime\prime} \\ \sqrt{{(B{(r-50)}^{\beta }+A{50}^{\alpha })}^{2}-X({\sigma }_{\mathrm{inst},2}^{2}-{\sigma }_{\mathrm{inst},{\rm{N}}}^{2}){r}^{2}}\\ \mathrm{for}\ r\gt 50^{\prime\prime} \end{array}\right.\end{eqnarray} \tag{ 4 }$$

where $A,B,\alpha ,\beta$ are the same as in Equation (3), X is a constant given in Table 3 (which varies between bands) to account for beam area and pixel size, and ${\sigma }_{\mathrm{inst},{\rm{N}}}$ (mJy beam⁻¹) is the instrumental noise for ${N}_{\mathrm{scan}}=N$ as given in Table 1. This equation uses the values of instrumental noise averaged over all the pixels with the same number of scans (Table 1). Purists interested in using the actual instrumental noise at the position of a source can measure this instrumental noise from the noise map, and then obtain the total flux error by modifying Equation (4) in a fairly obvious way. If users wish to use elliptical apertures they could either run their own Monte Carlo simulation on the released maps, or a reasonable estimate of the flux error can be obtained by using the estimate for a circular aperture with the same area.

The width of the SPIRE filters mean that both the size of the PSF and the power detected by SPIRE depend on the spectral energy distribution (SED) of the source. The SPIRE data-reduction pipeline is based on the assumption that the flux density of a source depends on frequency⁻¹, and all our images are ultimately based on this assumption. If the user has reason to know the SED of a source, the flux densities should be corrected using the corrections from either Table 5.7 or 5.8 from the SPIRE handbook.¹⁷ It is important to apply these corrections, since they can be quite large: for a point source with a typical dust spectrum (T = 20 K, β = 2) the multiplicative correction is 0.96, 0.94, and 0.90 at 250, 350, and 500 μm, respectively.

Finally, on top of the other flux density errors, there is an error from the uncertainty in the basic flux calibration of the instrument. At the time of writing, the error in the flux density arising from the uncertainty in the absolute flux density of Neptune is 4%, and there is an additional 1.5% error that is uncorrelated between the SPIRE bands (SPIRE Data Reduction Guide). The current recommendation is that these factors should be added linearly, and so the reader should use a calibration error of 5.5%.

The primary science goals of H-ATLAS are to investigate individual sources, and so our maps were made to optimize the detection and flux-extraction of these small-scale structures. Pascale et al. (2011) used simulations of our observing strategy and map-making techniques to show that there is attenuation of the structure in the H-ATLAS maps on scales >20'. Since Pascale et al. (2011) all-sky maps produced by the Planck observatory (Planck Collaboration et al. 2011) have been released (Planck Collaboration et al. 2016) which provide a useful "truth" map to compare with the H-ATLAS maps at 350 μm and 500 μm. In this section, we calculate the 1D angular power spectrum of our maps (using the agpy package¹⁸ ) to investigate what emission scales are preserved in our maps, and if our maps are consistent between tiles.

In Figure 11 we show the 1D power spectrum from the raw SPIRE maps, as well as for each individual tile in the mosaics. The power spectra for each field tend to be in good agreement with each other, especially at 350 and 500 μm. The differences between fields is most likely explained due to differences in the cirrus emission. Given the good agreement between tiles in a field, an average "transfer function" describing the depression of power as a function of angular scale, could be used for each of the two fields.

Zoom In Zoom Out Reset image size

To test whether the differences seen in Figure 11 are due to variations in the cirrus emission we compare our maps at 350 μm with the Planck Public Release 2 maps (Planck Collaboration et al. 2016). We first convolved both the NGP and SGP maps to the same resolution as Planck using the effective Planck beam for our field. Both the Planck and SPIRE maps were matched to the same 36'' pixel grid and converted to the same flux-density units. Figure 12 shows the ratio of the 1D power spectra of the SPIRE and Planck maps. The ratio maps show broad agreement across all individual tiles and mosaics, confirming that the differences in Figure 11 is due to cirrus emission. At small spatial scales the low ratios of SPIRE/Planck power are due to the greater sensitivity of the H-ATLAS observations. The low ratios at higher spatial scale are due to the finite size of the maps and the fact that a single SPIRE observation made from scanning the telescope in a single direction will miss large-scale power in the direction orthogonal to the scan direction (Waskett et al. 2007). Our results show the attenuation of emission begins on a scale of ∼15', broadly in agreement with the value of 20' found by Pascale et al. (2011). It is possible that our observations affected by "cooler burps" (see Section 3.1) could have greater attenuation, but, due to the complexity of isolating these regions we did not investigate this further.

Zoom In Zoom Out Reset image size

In principle, it is possible to correct for this attenuation on large scales using an alternative map-maker (Waskett et al. 2007). If users wish to create maps with alternative map-makers, the authors can be contacted to assist with data/customised timelines. For the 350 and 500 μm bands it is also possible to combine the Herschel and Planck data to create maps that have the correct power on all scales. Of course, the images on large scales will also be affected by cirrus.

The PACS maps are very different from the SPIRE maps. The higher instrumental noise means that source confusion is less important and the instrumental noise is correlated between pixels. It is more challenging to measure the confusion and instrumental noise on the PACS maps because Jscanamorphos uses the existence of multiple PACS observations to remove the effect of temporal changes in the detectors, which means that it is not possible to use jackknifes to estimate the instrumental noise.

To estimate the PACS confusion noise, we used a similar approach to that of Magnelli et al. (2013), who estimated the confusion in the GOODS-S field. We measured the total noise in regions of the map with different number of observations as seen in the ${N}_{\mathrm{scan}}$ maps shown in Figures 3 and 5. To measure the noise, we fitted a Gaussian to the negative part of the histogram of pixel values, using the positive side as an upper limit (similar to what we did for SPIRE in Section 5.2). This gave us a plot of ${\sigma }_{\mathrm{pix}}$ versus ${N}_{\mathrm{scan}}$. We then fitted a simple model to this relationship The model has an instrumental noise component, which scales with the number of observations (${N}_{\mathrm{scan}}$) contributing to each pixel, and a constant confusion term:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{pix}}(\mathrm{mJy})=\sqrt{{({\sigma }_{\mathrm{inst}}{N}_{\mathrm{scan}}^{-0.5})}^{2}+{\sigma }_{\mathrm{conf}}^{2}}\end{eqnarray} \tag{ 5 }$$

where ${\sigma }_{\mathrm{pix}}$ is the total pixel noise in mJy, ${\sigma }_{\mathrm{inst}}$ is the instrumental noise in mJy for a single PACS observation (${N}_{\mathrm{scan}}=1$), and ${\sigma }_{\mathrm{conf}}$ is the confusion noise (in mJy). In principle, this procedure allowed us to estimate ${\sigma }_{\mathrm{inst}}$ and ${\sigma }_{\mathrm{conf}}$.

We initially applied this method to our final maps, but found that the noise in the regions in which the tiles overlap is significantly reduced due to the re-projecting procedure used to create the mosaics.¹⁹ We therefore decided to use only the individual tiles, which limited the range of ${N}_{\mathrm{scan}}$ to 2–5, reducing the sensitivity of the method. To regain the sensitivity, we used some observations from the HeViCS survey (Davies et al. 2010), which mapped ∼55 sq. deg. of the Virgo Cluster using the same fast-scan parallel observing mode that we used. While most of the Virgo Cluster was observed in $4^\circ \times \ 4^\circ$ tiles with eight observations per field, the northernmost Virgo tile was observed 10 times by PACS. We reduced the observations of this tile using the same Jscanamorphos method we used for H-ATLAS, starting with the level-1 data produced by the standard pipeline. We then made five independent maps from each pair of observations and applied Nebuliser to each map, which gave us five maps of the same region of sky. We then averaged various combinations of maps and estimated the total noise on each combined map using the method above, giving us estimates of the total noise from ${N}_{\mathrm{scan}}=2\mbox{--}10$. The results of pixel noise versus ${N}_{\mathrm{scan}}$ are shown in Figure 13 for both the H-ATLAS and HeViCS results. The values of ${\sigma }_{\mathrm{inst}}$ and ${\sigma }_{\mathrm{conf}}$ obtained from fitting Equation (5) to the results for the individual fields are given in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  PACS Pixel Noise Model Parameters

Waveband	Field	${\sigma }_{\mathrm{inst}}$	${\sigma }_{\mathrm{conf}}$
		(mJy)	(mJy)
100 μm	NGP	3.578 ± 0.013	0.393 ± 0.047
	SGP	3.539 ± 0.030	0.603 ± 0.059
	HeViCS	3.774 ± 0.002	0.184 ± 0.004
160 μm	NGP	3.515 ± 0.017	0.389 ± 0.055
	SGP	3.532 ± 0.017	0.380 ± 0.058
	HeViCS	3.714 ± 0.001	0.240 ± 0.003

Note. The best-fit parameters for the relationship between pixel noise and number of scans (Equation (5)). See Section 6.1 for details.

Download table as:  ASCIITypeset image

As expected, the estimates of instrumental noise in Table 4 are much higher than the estimates of confusion noise. The errors on the confusion noise estimates for the HeViCS field are much less than those for the NGP and SGP fields because of the larger range of ${N}_{\mathrm{scan}}$. The confusion noise estimates for the different fields are formally inconsistent, which we suspect arises because the instrumental noise is so much larger than the confusion noise, making any estimate for estimating the confusion noise sensitive to systematic errors (e.g., if the assumption that the noise is Gaussian is slightly wrong). Our most reliable estimates of confusion noise come from the HeViCS tile because of the larger range of ${N}_{\mathrm{scan}}$, and are 0.184 ± 0.004 and 0.240 ± 0.003 mJy, at 100 and 160 μm, respectively. These estimates are broadly similar at 100 μm, but differ at 160 μm, to those presented by Magnelli et al. (2013) of 0.15 and 0.68 mJy at 100 and 160 μm, although these values may not be directly comparable due to differences in pixel size and beam size. Assuming beam areas of 207 and 308 arcsec² at 100 and 160 μm (calculated from our measured PSFs), the confusion noise is 4.23 and 4.62 mJy beam⁻¹ at 100 and 160 μm, respectively.

The PACS PSF is not a simple Gaussian and in fast-scan parallel mode is significantly extended in the scan direction (Section 4), which means that it must vary within both fields, especially between points on the maps that are composed of different numbers of individual observations. For this reason, the technique of maximising the signal-to-noise for point sources by convolving the images with the PSF is not as accurate as for the SPIRE images. Instead, we use aperture photometry with a small aperture (see below). However, if the reader does prefer to convolve the map with the PSF, for example for detecting a faint point source, we recommend the use of our Gaussian fit to the empirical PSFs, which have a FWHM of 11.4 and 13.7 arcsec at 100 and 160 μm, respectively. Anyone carrying out a stacking analysis should be aware that the means of the PACS maps are not zero, and so the mean of the map should be subtracted before proceeding with the stacking.²⁰ Errors for a stacking analysis should be obtained from a Monte Carlo analysis in which flux densities are measured at random points in the image.

For photometry of an unresolved source, the alternative to measuring the flux density from an image that has been convolved with the PSF is aperture photometry with an aperture not much larger than the PSF. V16 found that the signal-to-noise peaks for an aperture with a radius of $\simeq 8$ arcsec at both wavelengths. We suggest that astronomers wishing to carry out photometry of point sources should use this aperture, although since such a small aperture contains only a small number of pixels, they should think carefully about pixelization effects when using this approach. The units of the PACS maps are Jy pixel⁻¹, so aperture photometry can be carried out by adding up the flux density values for all the pixels within the aperture; there is no need to estimate a sky value because we have already subtracted any residual background emission using Nebuliser. 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 (see Section 4). Both the flux densities and the flux errors (see below) should be corrected for the fraction of the PSF that is outside the aperture using this table.

Photometry of sources that are expected to be extended, for example nearby galaxies, should also be carried out by aperture photometry. There is no need to estimate a sky value because we have already subtracted any residual background emission using Nebuliser. Both the flux densities and the flux errors (see below) should be corrected for the fraction of the PSF that is outside the aperture using the EEF.

We have carried out a similar Monte Carlo simulation to estimate the errors in the flux densities measured in aperture photometry to the one we carried out for SPIRE (Section 5.3). For each aperture radius, we placed 3000 apertures at random positions on the part of each image with ${N}_{\mathrm{scan}}=2$. We used radii ranging from approximately the size of the beam up to 100 arcsec. The results are shown in Figure 14. As for SPIRE, we fit a power-law relationship (Equation (3)) to the results of the Monte Carlo simulation. The values of the best-fit parameters in this relationship are given in Table 3. As for SPIRE, we found that at large radii the flux error is approximately proportional to the aperture radius, which is the relationship expected for instrumental noise that is not correlated between pixels. At smaller radii, as for SPIRE, we found the flux error increases more rapidly with radius. We are not sure of the explanation but possibilities include cirrus emission that has not been removed because of the large filtering scale used in Nebuliser and residual 1/f noise not removed by the map-maker. We have not produced maps of the PACS instrumental noise for the data release. Instead, we have produced maps showing the number of individual data sets (${N}_{\mathrm{scan}}$) contributing to each pixel. This map and the following equation can then be used to obtain an estimate of the flux density error for any object and any aperture size:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{ap}}(\mathrm{mJy})\\ &&=\,\left\{\begin{array}{l}\sqrt{{({{Ar}}^{\alpha })}^{2}-\displaystyle \frac{{\sigma }_{\mathrm{inst}}^{2}\pi }{{P}^{2}}\left(\displaystyle \frac{1}{2}-\displaystyle \frac{1}{{N}_{\mathrm{scan}}}\right){r}^{2}}\\ \mathrm{if}\ r\leqslant 50^{\prime\prime} \\ \sqrt{{(B{(r-50)}^{\beta }+A{50}^{\alpha })}^{2}-\displaystyle \frac{{\sigma }_{\mathrm{inst}}^{2}\pi }{{P}^{2}}\left(\displaystyle \frac{1}{2}-\displaystyle \frac{1}{{N}_{\mathrm{scan}}}\right){r}^{2}}\\ \mathrm{for}\ r\gt 50^{\prime\prime} \end{array}\right.\end{eqnarray} \tag{ 6 }$$

where A and B are best fit parameters from Table 3, ${\sigma }_{\mathrm{inst}}$ is the noise in mJy given in Table 4 and P is the pixel size of the maps in arcseconds (3'' at 100 μm and 4'' at 160 μm). As with SPIRE, if the user wishes to use bespoke elliptical apertures a reasonable estimate can be obtained by using the estimate of the flux uncertainty for a circular aperture with the same area, or alternatively, they could run their own Monte Carlo simulation on the released maps.

Zoom In Zoom Out Reset image size

On top of the flux density uncertainty given by our power-law model, there is also a fundamental calibration error. As for SPIRE, the dominant uncertainty is due to the models of the calibration objects, in the case of PACS stars and asteroids, which is estimated to be 5% (PACS Calibration page²¹ ). The reproducibility of calibration sources is measured to be ∼2% (Balog et al. 2014) and so, as in V16, we add the uncertainties and thus make the conservative assumption that the calibration uncertainty is 7%. As with SPIRE, all our measurements of flux density are based on the assumption that flux density is proportional to frequency⁻¹, which introduces an error if the source does not have this SED because of the width of the PACS bandpass filters. We refer anyone wishing to make a correction for this effect to the PACS Color-Correction document.²²