We have reduced the data taken with the Spectral and Photometric Imaging Receiver (SPIRE) photometer on board the Herschel Space Observatory in the Science Demonstration Phase (SDP) of the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS). We describe the data reduction, which poses specific challenges, both because of the large number of detectors which can have noise correlated in each array, and because only two scans are made for each region. We implement effective solutions to process the bolometric timelines into maps, and show that correlations among detectors are negligible, and that the photometer is stable on time scales up to 250 s. This is longer than the time the telescope takes to cross the observed sky region, and it allows us to use naive binning methods for an optimal reconstruction of the sky emission. The maps have equal contribution of confusion and white instrumental noise, and the former is estimated to 5.3, 6.4 and 6.7 mJy beam−1 (1σ), at 250, 350 and 500 μm, respectively. This pipeline is used to reduce other H-ATLAS observations, as they became available, and we discuss how it can be used with the optimal map maker implemented in the Herschel Interactive Processing Environment (HIPE), to improve computational efficiency and stability. The SDP data set is available from http://www.h-atlas.org/.
The Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS, Eales et al. 2010) is currently the largest area open-time key project that is being carried out by the Herschel Space Observatory (Pilbratt et al. 2010). The 600 h of observations will map ∼550 deg2 of the sky in the five far-infrared (far-IR) and submillimetric bands of the Photodetector Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRE, Griffin et al. 2010) photometers. The survey of unprecedented sky and spectral coverage, and angular resolution, will target sky regions selected in the direction of the North and South Galactic Poles, and of the Galaxy And Mass Assembly Survey, the GAMA fields (Driver et al. 2009).
The first observations of the project have been completed on a 16 deg2 GAMA target, as part of the initial Science Demonstration Phase (SDP). Here we present the maps obtained with the SPIRE photometer, and the reduction pipeline, which differs in significant aspects from the standard pipeline described by Griffin et al. (2008, 2010). The SDP data are publicly available,1 and comprise a set of three SPIRE photometer maps at 250, 350 and 500 μm (discussed in this work), and two PACS photometer maps at 100 and 160 μm (Ibar et al. 2010). Each map has its corresponding error and coverage maps as separate extensions. Sources (point-like and extended) are detected in these maps by Rigby et al. (2010), and used by Smith et al. (2010) to identify counterparts in catalogues of ultraviolet and near-IR selected sources.
SPIRE uses arrays of sparsely distributed bolometric detectors to sample the sky radiation through scanning at a specific angle with respect to the bolometer array axes, and the challenge of map making or sky reconstruction is a linear algebra problem which has been extensively investigated in the past (see section 8 of Cantalupo et al. 2010 for a list of references). The bolometric data are combined to provide optimal estimates of the sky emission, for which detailed knowledge of the detector noise properties is required. Most experiments (including the high-frequency instrument on Planck) have only a limited number of detectors, sparsely populating their focal planes, but the large format bolometric cameras, which are now available, pose new challenges, as the amount of data is increasing, with potential for correlation among detectors, complicating the analysis (Patanchon et al. 2008). This is of relevance for SPIRE, and the detector noise properties need to be investigated in order to optimally reduce the data into maps.
The choice of a map making technique is dictated by the noise properties of the data. In this paper we analyse the noise in the bolometric timelines, and we use this information to generate two sets of maps: one optimized to study emission from point sources, and the other to study the large-scale structure accessible in the observed sky region. The former have been used in the scientific analyses of the SDP H-ATLAS data published in the Herschel First Results special A&A edition (2010).
This analysis procedure provides a first investigation into the properties of this data set, and our results can be used in future studies to select appropriate map making techniques, which may improve on the one used here. However, the pipeline we discuss here will provide the necessary input to these techniques, in order to optimize their performance.
2 THE SPIRE PHOTOMETER
The SPIRE photometer images the sky in three submillimetric bands centred at 250, 350 and 500 μm. Each band is 30 per cent wide (in Δλ/λ) and the focal plane is populated with arrays of silicon nitride micromesh ‘spider web’ bolometers (Bock et al. 1996, 1998) operated at a cryogenic temperature of ≲300 mK. The detectors are coupled to the spacecraft 3.5-m mirror telescope by smooth walled feed horns, and 4.5 K re-imaging optics, which also provide shielding from stray radiation. The linear distance of adjacent feed horns is 2fλ (Rownd et al. 2003) where f is the focal ratio and λ is the wavelength. Compared to its balloon-borne predecessor, BLAST (Devlin et al. 2009; Pascale et al. 2008), SPIRE has better angular resolution, with point spread functions (PSFs) of 18.1, 25.2 and 36.6 arcsec full width at half-maximum (FWHM) at 250, 350 and 500 μm, respectively (Griffin et al. 2010). Each array has, respectively, 139, 88 and 43 light sensitive bolometers, as well as two dark bolometers and one resistor. Thermistors on each array monitor temperature variations in the substrate, due to fluctuations in the thermal bath. The bolometers are AC-biased, and the voltage RMS measured across each detector is proportional to the radiation detected from the sky. The constant of proportionality is the detector’s responsivity, which is measured by Swinyard et al. (2010) from observations of sources with known spectral energy distributions.
3 THE DATA
The SDP H-ATLAS data was obtained in parallel mode, with the two Herschel instruments, PACS and SPIRE, operating simultaneously. The H-ATLAS data set consists of two observations – one at the nominal SPIRE scan angle and the other with scans at nearly orthogonal direction (∼85°), over the same 16 deg2 square region centred at coordinates RA = 9h5m30s, Dec. = 0°5′0′. Each observation is carried in ∼8 h by scanning the telescope at the constant angular velocity of 60 arcsec s−1 over great circles. The scan lines are separated by 155 arcsec in order to achieve good coverage for large fields for both SPIRE and PACS maps. During the scanning the SPIRE bolometers are sampled2 at 10 Hz.
4 MAP MAKING METHOD
The well known Generalized Least Squares (GLS) estimator of the sky (e.g. Lupton 1993) is
When the timeline noise has a 1/f-type component, high-pass filters can be used to whiten, or suppress a drifting baseline. This results in an attenuation of the map at angular scales larger than λc=fc/ω, where fc is the filter’s cut-off frequency, and ω the scan velocity. As shown later in Section 7, the outstanding stability of the SPIRE photometer makes it possible to set λc at the map scale. Therefore essentially all Fourier modes are preserved in the telescope scan direction. Cross-scan modes, however, are suppressed, and can only be recovered by scanning at an orthogonal angle. This can be seen when studying equation (2) in Fourier space, as discussed in Patanchon et al. (2008). The two SDP H-ATLAS observations can be reduced into a map using a naive binning, and simple co-addition. Since each map preserves only one set of modes, the final, coadded map suppresses (by up to a factor 2) the large scales, and therefore it is no longer equivalent to the GLS solution of equation (2). It is, however, possible to recover the suppressed modes, since all the necessary information is present and the filtering effect can be deconvolved. This is demonstrated in Section 10, by using the map maker transfer function.
5 DATA PREPROCESSING
Each detector timeline is preprocessed using a customized version of the SPIRE standard pipeline (Griffin et al. 2008, 2010), implemented in the Herschel Interactive Processing Environment (HIPE), which is run in an interactive mode, starting from the Level-0.5 data products (DP). These are a set of FITS binary tables grouping together the uncalibrated light and dark detector timelines, and all relevant housekeeping information. Each table contains data from just one scan, so there is a large number of DP tables for each observation.
Pointing information is associated with each bolometric sample using the standard pipeline, and the timelines are converted to flux density per beam by applying the calibration of Swinyard et al. (2010, the calibration is further discussed in Section 9).
Glitches (mostly cosmic ray events or Bremsstrahlung within the instrument) in the bolometer and thermistor data are flagged, and replaced with a constrained realization of white noise. Two alternative deglitchers are used (alternative in the sense that they are available, but not used by default in the standard pipeline) to identify these transients. One localizes synchronous transients among detectors in a given array, and the other is a sigma–kappa deglitcher which looks for isolated events. This allows the detection of a larger number of low signal-to-noise ratio glitches which could create false point sources in the final maps, if not corrected. We have verified that the software implementation of these deglitchers does not erroneously flag astronomical sources, which have a significantly different spatial signature in the timelines.
Electrical DC-offset jumps of unknown origin are occasionally present in the thermistor and bolometer timelines. These jumps cause temporary decorrelation of the thermistor and bolometer timelines, so that usual attempts to correct for the temperature drift in the bolometer thermal bath (see the next section) fail. Unless corrected for, these jumps cause large scale artifacts (stripes) in the final maps, at several sigma level.
At the time of these observations, the standard pipeline did not have the capability to search for and correct the DC-offset jumps. Therefore transient events have been identified by eye and compensated for by shifting the affected timelines by an appropriate DC-offset, which is estimated from the data. We have verified that the properties of the timeline power spectrum are preserved, showing that artifacts are not introduced by applying this correction.
The electronic and bolometric instrumental transfer functions are deconvolved from the timelines (using the standard pipeline) to recover a flat spectral response, and to remove filtering phase delays.
Data products of each scan are linked together to provide continuous timelines for each observation (including spacecraft turnaround). There is a clear advantage in doing this in the analysis which follows, because signal processing Fourier methods (filtering, but also power spectral analyses) perform better, and the unavoidable artifacts in the data are confined at the map edges, corresponding to the beginning and the end of each observation. Temporal gaps in the data acquisition, corresponding to the SPIRE PCAL source flashes and other events, are filled with a constrained realization of white noise, and flagged accordingly. Flagged data are not projected into maps, but filling the gaps with noise improves the stability of filters, and reduces pollution of good data at the edges of gaps.
The final preprocessed timeline is a clean (without glitches or DC jumps) and calibrated data-set, with associated pointing information, which can be binned into a map.
6 CORRELATED NOISE
Each light detector timeline contains signal from the sky, together with noise. In a large format array the noise has two components. One is uncorrelated among detectors, the other arises from common signals:Patanchon et al. 2008; Pascale et al. 2009).
The dominant common signal in the SPIRE arrays correlates well with the substrate temperature measured by the array thermistors, as shown in the top panels of Fig. 1, where representative detector timelines in the 250, 350 and 500 μm arrays are plotted against the corresponding thermistor timelines. The Pearson correlation coefficients between bolometers and thermistors are all greater than 99 per cent. The coefficients αi are estimated from the slopes of these plots. For each bolometer timeline a temperature component is obtained from the multiplication of αi with the thermistor. This temperature component is then simply subtracted from the bolometer timeline to yield the temperature-drift-corrected bolometer timeline. To avoid introducing excess noise, each thermistor is low-pass filtered at 20 mHz first. Since gradients are not well represented in Fourier space, a fifth order polynomial is fitted to the thermistor signals and subtracted before filtering, and then added back afterwards.
The corrected timelines are shown in the bottom panels of Fig. 1. The 250 and 350 μm flat timelines indicate that most of the correlated signal has been removed. However, this is less true at 500 μm. The temperature-corrected dark 500 μm bolometer (not shown in the plot) has a residual similar to the one at 250 and 350 μm, suggesting that the excess signal in the light detectors is optical pick-up. It is possible that the 500μm array is sensitive to the blackbody light radiated by the ≲4.5 K optics box cavity, which is not traced by the thermistor. At 500 μm this is 1.5 per cent of the peak emission (in νIν of a 4.5 K blackbody), while it is negligible at the shorter wavelengths. Although it might be possible to decorrelate a second component from the 500 μm timelines, this turns out not to be necessary, as shown in the following section. We also notice that this residual is present only where the temperature monitored by the thermistor has a gradient in excess of ∼1 nV s−1, measured in the timeline units.
It should be noted that dark bolometers might also be used to decorrelate the temperature component from the timelines with similar results. This is useful when the thermometer timelines are not available. However we have not addressed in our study the possible issue of light cross talk between light and dark bolometers.
7 NOISE ANALYSIS
To understand the noise properties of the clean, temperature-corrected bolometers, it is necessary to estimate noise timelines. This is done by suppressing the astronomical signal in the timelines, because cirrus and point sources have a contribution which can be above the noise in individual samples. Naive maps with 6, 10 and 14 arcsec pixels, at 250, 350 and 500 μm, respectively, are generated by binning the timelines of the two SDP observations. These pixel sizes are the default in the standard pipeline. New timelines are obtained for each detector by rescanning the map with the telescope pointing solution, effectively simulating new observations. Since each map pixel contains ∼10 samples, the instrumental noise in the rescanned timelines is about 1/3 that of the original timelines. The rescanned timelines are subtracted from the cleaned bolometers to provide estimates of detector noise. The level of correlations introduced with this operation is negligible, and the noise estimates from the subtracted timelines can be biased low by 5 per cent, at most. The noise power spectral density (PSD) is shown in Fig. 2, for a representative 250 μm detector.
The model fit values for f0 are below 4 mHz across the arrays, corresponding to angular modes larger than (i.e. the same as the survey angular size), at the scan rate of 60 arcsec s−1. The estimates of the white noise level are, w0≃ 8.6, 9.3 and 13.2 mJy s0.5 (10 per cent estimated variability across the arrays) at 250, 350 and 500 μm, which are compatible with the results of Nguyen et al. (2010), within their quoted uncertainties.
Residual correlations among detectors are investigated by estimating the cross-power spectum between channels in an array. An example for two representative 250 μm detectors is also shown in Fig. 2. It is significant at frequencies below the noise ‘knee’, meaning that the residual 1/f noise is highly correlated among detectors. This is not a problem for map making because these temporal scales (indicated by a vertical dotted line in the figure) can be safely filtered, as they all correspond to unconstrained angular modes (next section). However, we notice that the cross-correlation is significantly different from zero at high frequency, indicating the presence of broad-band correlated noise of instrumental origin. The same level of correlation is observed in the raw timelines, indicating that this is not an artefact introduced by our analysis. The level of this correlation is estimated at 1 Hz as the ratio of the cross-power spectrum to the PSD, which is 2, 10 and 25 per cent, at 250, 350 and 500 μm, respectively, and negligible among detectors in different arrays. The origin of the effect is unclear; however, it is not a concern for map making because the correlated noise component is subdominant relative to the uncorrelated component, and because different detectors in an array point at different directions in the sky, at any given time.
A fifth order polynomial is fitted and subtracted from each temperature corrected timeline, to remove residual gradients. A line is fitted to the first and last 10 s of data of each timeline, and subtracted. Temporal frequencies longer than the 1/f‘knee’ are suppressed by high-pass filtering each timeline at the frequency of 4 mHz, which corresponds to on the sky, at the survey scan rate. This preserves all the accessible angular modes in one observation (the modes along the scan direction) while suppressing the correlated noise which would otherwise complicate the analysis or introduce artifacts. The weak filtering adopted does not distort the PSF in a significant way, as discussed in the next section. This set of timelines is the final, improved Level-1 DP, which forms the input to the map maker.
The naive map maker implemented in the SPIRE standard pipline is used to bin the improved Level-1 data into maps with pixel sizes of 5, 10 and 10 arcsec, at 250, 350 and 500 μm, with identical astrometry. This pixelization makes it trivial to rebin the 250 μm map on to the 350 and 500 μm grid, when needed, without loss of information, and Nyquist sample the PSF at each wavelength. These pixel sizes are smaller than that used by the standard pipeline (6, 10 and 14 arcsec), and caution should be taken as a few pixels will have no bolometer samples (i.e. zero coverage), which means they will be left blank. This is not a problem because blank pixels are properly flagged in a mask which is also generated by the map maker. Three maps are generated: the first two are obtained by binning data from the two individual, cross-linked observations (M1 and M2), and the other is a combined map (Mc) made using timelines from both observations. Coverage and noise maps are also generated.
The two separate SDP scan maps are similar in coverage, and number of samples in each pixel. Therefore, it is reasonable to assume similar noise properties in both. By subtracting one map from the other, we obtain a jackknife which does not contain astronomical signal, and can be used to characterize the instrumental component of the noise in the maps. A χ2 test of Gaussianity on the jackknife shows that the instrumental noise is well described by a Gaussian distribution, with reduced χ2 = 1.3, 1.4 and 1.3 and Q(χ2/ν) = 15, 7 and 14 per cent, for 250, 350 and 500 μm, respectively (ν are the degrees of freedom).
The instrumental noise RMS, σ, is estimated from the jackknife pixel histogram, and we find σ = 29.6, 31 and 36 mJy (i.e. the standard deviation of each sample in a timeline), at 250, 350 and 500 μm, respectively. These values are compatible with the noise estimated from the timelines, within ∼15 per cent, and are in good agreement with the noise estimated by the map maker, which is provided in the form of a noise map (estimated from the variance in a pixel, divided by the number of samples in the pixel).
We also estimate the confusion noise, the signal arising from faint extragalactic sources blended together by the finite telescope PSF. Smoothed maps are produced by the inverse variance weighted convolution with the telescope PSF, and a region of low cirrus emission is selected, from which a pixel histogram is formed. This histogram represents the combined contributions of instrumental noise, confusion noise and emission from detected sources within the region (Patanchon et al. 2009; Marsden et al. 2009). A Gaussian function is fitted to the negative part of the histogram only, whose width is the square root of the quadrature sum of confusion and instrumental noise, the latter estimated from the pixel histogram of the smoothed jackknife.
The 1σ instrumental noise in the PSF-convolved map is estimated to be 4.8, 4.9 and 5.7 mJy beam−1 at 250, 350 and 500 μm, respectively, with 5 per cent uncertainty. The confusion noise is, respectively, 5.3, 6.4 and 6.7 mJy beam−1, with a 7 per cent uncertainty (calibration uncertainties are not included). These can be considered as the relevant values when extracting point sources from the maps.
The astrometry has been verified using the point sources detected in the 250 μm SDP maps (Rigby et al. 2010). We produced histograms of the separations in RA and Dec. for all SDSS DR7 r-band sources within 50 arcsec of the 250 μm centroids, in bins of 1 arcsec. We then fit the resulting distribution of separations using a Gaussian model for the SPIRE positional errors and allowing for the effects of galaxy clustering in the SDSS data (see Smith et al. 2010). We find a small, ≲2 arcsec astrometric shift (≲1 arcsec uncertainty), which is compatible with pixelization effects, and is accounted for by shifting the maps. A similar analysis using stacking of a catalogue of IRAS sources, with positions determined in the optical and radio, gives compatible results at each wavelength. It is interesting to note that Pilbratt et al. (2010) report a 1σ pointing uncertainty of ∼2 arcsec, similar to what we find, although stochastic in nature.
The flux calibration is carried out in the timelines as part of the standard SPIRE data reduction (as discussed in Section 5), and the resulting maps are calibrated in Jy beam−1. We applied the flux correction factors of Swinyard et al. (2010) and Griffin et al. (2010) required by the SPIRE calibration available at the time of this work. The 250, 350 and 500 μm maps are multiplied by 1.02, 1.05 and 0.94, respectively.3 The map making process effectively convolves the sky emission with the pixel window function, introducing a gain variation for point sources of up to 5 per cent for the 250-μm array. Instead of accounting for this effect in the maps, we prefer to include it in the PSFs used to retrieve point source fluxes. Gaussian PSFs are fitted to maps of Neptune, a SPIRE primary calibrator. The PSFs are then normalized using the Neptune flux model of Moreno (2010), calculated for the day when Neptune was observed. The PSFs, available as part of the SDP release, have FWHM of 18.1, 24.8 and 35.2 arcsec, and peak amplitudes of 0.95, 0.95 and 1.02 at 250, 350 and 500 μm, respectively. Future release of H-ATLAS products will make use of PSFs estimated from Neptune maps.
10 SKY RECONSTRUCTION
In this section we investigate how well the pipeline is able to provide an estimate of the sky emission. As the noise has already been discussed, it is now necessary to study the effect that high-pass filtering the timelines has on the maps. This operation is equivalent to convolving the sky with an effective filter or transfer function
As mentioned in Section 4, Fourier modes orthogonal to the scan direction are lost in the 1/f noise, which is filtered in the timelines. The two transfer functions for each observation are shown in Fig. 3. These have been estimated from 900 realizations of white noise sky emission, reconstructed by the pipeline, and show how the largest Fourier angular modes are the most affected.
The maps preserve a unit response to point sources. To verify this, noiseless maps with Gaussian PSFs of nominal sizes for the three arrays are generated. The timelines obtained by re-scanning the maps are then processed by the pipeline into new maps. From the comparison between the input and output maps, we find that the gain variation is below 1 per cent and therefore negligible, compared to other sources of error, particularly the 15 per cent calibration uncertainty estimated by Swinyard et al. (2010).
The large angular scale Fourier modes which are lost in one observation, are constrained by the other. It is therefore possible to combine the two maps with no loss of information, when the transfer functions are known. A Maximum Likelihood (ML) estimator of the sky is in this case (assuming a constant and equal noise in each map)
To verify how effective this estimator is, tests were performed using a noiseless sky realization of Galactic cirrus emission, assuming a k−3 power spectrum (e.g. Miville-Deschênes et al. 2007), where k is the inverse angular scale. After reconstructing the emission with our pipeline, the input and reconstructed power spectra were compared (Fig. 4). The effect of attenuating the large scales is shown by the dashed line. This affects Fourier modes which are larger than 20 arcmin on the sky (and therefore irrelevant for point sources). However, the power spectrum of the ML estimator in equation (6) (shown by the black line in the figure) demonstrates that all accessible Fourier modes have been properly reconstructed.
It should be noted that the scientific analysis of H-ATLAS data taken during the SDP phase were constrained to point sources or moderately <2 arcmin extended sources only. Since these sources are not affected by this choice of filtering, the naive maps were used, instead of the ML ones described here.
Because detectors can be effectively decorrelated, it would be possible to use the well-tested Cosmic Microwave Background code madmap (Cantalupo et al. 2010) in the future. This is an optimal map maker, which does not currently handle detector–detector correlations, and does not require an estimate of the map making transfer function to compute the GLS sky estimator. The detector noise can be made stationary (at least over the time spanned by our data set) and the noise model required by madmap can be efficiently precomputed and consistently used to reduce the whole H-ATLAS ∼550 deg2 data set, when it becomes available. A madmap implementation is now being optimized in the SPIRE pipeline, and preliminary tests suggest that it takes ≲30 min to reduce the SDP data set on a single CPU machine, provided that enough memory is available for the data.
Data from the SDP H-ATLAS observations have been reduced into maps using a customized pipeline. Noise analysis of individual bolometer timelines demonstrates that the SPIRE photometer is very stable, with 1/f residual noise dominating frequencies below 4 mHz. At the scan velocity adopted for this survey, this corresponds to the angular scale of the survey region. Naive map making can therefore be used to provide the optimal reconstruction of point source sky emission. However, to constrain the large scale structure it has been necessary to estimate the map transfer function by using computationally intensive Monte Carlo techniques.
The H-ATLAS maps contain Gaussian instrumental noise, in addition to extragalactic confusion noise, and the two have roughly the same contribution. Estimates of the confusion noise are in agreement with those measured from alternative studies. Residual correlations among detectors in the same array are found to be subdominant compared to the uncorrelated (among detectors) noise component, at all relevant time scales.
This map making technique represents an intermediate solution until the implementation of madmap in the SPIRE standard pipeline is completely characterized and tested. Nevertheless, the improved Level-1 DPs, and noise models produced in this work will provide the necessary inputs in order for madmap to perform optimally. This seams likely to be the method-of-choice to reduce the ∼550 deg2 H-ATLAS Survey, as it does not require the explicit estimation of the transfer functions, and detectors can be efficiently decorrelated at the time and angular scales of interest, and moreover the noisy, correlated long temporal scales can be filtered without introducing artifacts.
The Herschel ATLAS is a project with Herschel, which is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. US participants in Herschel ATLAS acknowledge NASA support through a contract from JPL. We would like to thank George Bendo, Matt Griffin and Andreas Papageorgiou for the helpful discussions.