Wide-field radio sky surveys yield information for large samples of Galactic and extra-galactic objects, allowing analysis of emission physics, source populations and cosmic evolution of radio sources such as active galactic nuclei, star-forming galaxies, pulsars, and supernova remnants. Surveys in new regions of parameter space, such as frequency, are particularly useful for identifying new and unusual objects and adding constraints to emission models.
While radio surveys such as the Sydney University Molonglo Sky Survey (SUMSS; Bock, Large, & Sadler Reference Bock, Large and Sadler1999; Mauch et al. Reference Mauch, Murphy, Buttery, Curran, Hunstead, Piestrzynski, Robertson and Sadler2003) at 843 MHz and the Molonglo Reference Catalogue (MRC; Large et al. Reference Large, Mills, Little, Crawford and Sutton1981) at 408 MHz cover the southern sky, there is no deep survey of the southern sky below these frequencies. The Culgoora Circular Array (Slee Reference Slee1995) measured approximately 1 800 high-frequency selected sources over a Dec range of − 48° to + 35° at 80 and 160 MHz with limiting flux density of 4 and 2 Jy, respectively. However, these were targeted observations rather than a survey. A recent low-resolution (26 arcmin) survey by the Precision Array for Probing the Epoch of Reionisation (PAPER; Jacobs et al. Reference Jacobs2011) detected ≈ 500 sources at 145 MHz over 4 800 deg2.
The new generation of aperture-array telescopes such as the Murchison Widefield Array (MWA; Tingay et al. Reference Tingay2013; Lonsdale et al. Reference Lonsdale2009) and the Low-Frequency Array (LOFAR; van Haarlem et al. Reference van Haarlem2013) offer a chance to conduct all-sky low-frequency surveys at higher angular resolution, sensitivity, and speed than ever achieved before. Observations using the MWA 32-element prototype array such as Williams et al. (Reference Williams2012) and Bernardi et al. (Reference Bernardi2013) have demonstrated the scientific potential of low-frequency observations with the MWA design on a radio-quiet southern site.
The MWA is the only low-frequency precursorFootnote 1 for the Square Kilometre Array (SKA) and the first of the three SKA precursors to be operational for science. The MWA is located at the Murchison Radio-astronomy Observatory in outback Western Australia, the planned site of the future multimillion element SKA-low array. Thus, the MWA explores the characteristics of the site at low frequencies in the pursuit of challenging science, exercising much of the physical infrastructure that will be applied to the SKA over a geographic distance of 800 km (on-site infrastructure, long haul data transport, data archive infrastructure), and provides a valuable base for SKA verification systems.
This paper presents the Murchison Widefield Array Commissioning Survey (MWACS), work undertaken during the commissioning period of the full MWA, using a subset of its capabilities. The survey aimed to verify instrumental performance, verify our understanding of the MWA primary beam, motivate development of new data processing techniques, and create an initial sky model of brighter radio sources at MWA frequencies. A comprehensive sky model, in particular, is a pre-requisite for calibrating the full MWA, due to its huge field-of-view (Mitchell et al. Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008). The observed field also includes two of the three regions chosen for deep ( ≈ 1 000 h) integrations for the MWA’s Epoch of Reionisation key science program (see Bowman et al. Reference Bowman2013 for a summary of the MWA’s key science programs). MWACS is the first large-scale survey performed by the MWA in its role as an SKA precursor and the first such survey by an SKA precursor.
These observations cover approximately 6 100 deg2 roughly centred on the south Galactic pole. This is a region of sky with low brightness temperature, populated mostly by extragalactic sources, which we detect, catalogue and characterise in unprecedented detail at these frequencies. The MWACS is the first systematic exploration of the end-to-end MWA system and the characteristics of the SKA site at these frequencies, and lays a base of understanding for the much larger and more comprehensive GaLactic and Extragalactic All-sky MWA (GLEAM) survey, which is currently underway and will survey the entire sky south of δ = +20° between 72 and 230 MHz (Wayth in preparation). In turn, this will lay a base for the massive continuum surveys that will take place at low frequencies with the SKA, from the same site.
This paper is laid out as follows. Section 2 describes the observations, and the calibration and imaging strategy used in data reduction. We demonstrate how a hybrid strategy using conventional radio astronomy tools can generate high-quality mosaics by taking advantage of the MWA’s very good instantaneous u, v coverage and near-coplanarity. Section 3 describes the source-finding and flux density calibration procedures, including correcting for the MWA primary beam and how we have dealt with resolved structures. Section 4 describes the properties of the source catalogue and discusses issues affecting reliability and sensitivity. Section 5 contains a discussion and conclusion.
2 OBSERVATIONS AND DATA REDUCTION
As detailed by Tingay et al. (Reference Tingay2013), the MWA consists of 128 32-dipole antenna ‘tiles’ distributed over an area approximately 3 km in diameter. Each tile observes two instrumental polarisations, ‘X’ (16 dipoles oriented east–west) and ‘Y’ (16 dipoles oriented northsouth). The zenith field-of-view of a beamformed tile is ≈ 30° at 150 MHz. The signals from the tiles are collected by 16 in-field receiver units, each of which services eight tiles. For engineering reasons, during the commissioning period only four receivers were active at any one time, hence the tiles and receivers were commissioned as an overlapping series of six 32-tile sub-arrays labelled alpha through zeta. The sub-arrays were chosen to have good snapshot u, v-coverage within the various technical constraints in place during commissioning. The data presented in this paper were recorded by sub-arrays beta and gamma. The antenna layout and snapshot u, v-coverage of the combination of these two arrays are shown in Figure 1. The baselines are 8–1 530 m in length, and their combined effective angular resolution at 180 MHz is of order 3 arcmin (5 arcmin at 150 MHz and 6 arcmin at 120 MHz.). For simplicity, only the ‘XX’ and ‘YY’ polarisations are used in the following analysis; the cross-polarisation terms are discarded. The effect of ignoring these terms is constant throughout the night, as the instrument gains are very stable; any small loss in flux is later fixed during the flux calibration stage (Section 3.5).
A number of observation programs were undertaken during commissioning, one of which was night-time ‘drift scans’ where the MWA tiles were pointed to a single declination (Dec) along the meridian and data were collected as the sky drifted through the tile beams. This form of observation has previously been shown to be an effective way to observe large fractions of the sky with good sensitivity for the MWA (Bernardi et al. Reference Bernardi2013) and many other instruments, past and present, that use fixed dipole arrays also use drift scans. Using drift scans, uncertainty about the system calibration is minimised because the settings of all analogue components of the system are unchanged over the entire observation. The MWA has excellent stability of both the amplitude and phase of antenna gains using this mode, especially at night, when the ambient temperature changes slowly. We found that the standard deviation of primary-beam-corrected peak flux densities of typical bright unresolved (> 50 Jy) sources was only 1.5% as they drifted through the zenith beam. Typical commissioning calibration scans show phase stability of better than than 1° over the band. Overall, we found the quality of drift scan data to be limited by the stability of the ionosphere, which generates slow astrometric changes (see Section 3.9.3); in calm conditions a single calibration solution can be applied to all the data for an entire night, in the same fashion as Bernardi et al. (Reference Bernardi2013).
The data presented in this paper are from drift scans taken at two Dec settings: δ = −26.°7 (the zenith; hereafter referred to as Dec − 27) and δ = −47.°5 (hereafter referred to as Dec − 47). The observations are broken into a series of scans that cycle between three frequency settings with centre frequencies 119.04, 149.76, and 180.48 MHz (henceforth given as 120, 150, and 180 MHz); the scans are two minutes long and are separated by eight-second gaps. After each scan, the frequency is changed and a new scan begins. The bandwidth of the MWA is 30.72 MHz, hence the total frequency range of our observations is 104 to 196 MHz. Observations at a particular frequency are thus separated into many two-minute scans that begin every six minutes. In October, the night-time drift field-of-view encompasses right ascension (RA) in the range ≈ 20–09 h; we restricted our analysis to the range where both beta and gamma data were available at identical local sidereal times: 21 h < RA < 8 h. The observations are summarised in Table 1, and in total comprised 3 TB of unaveraged visibilities (reduced to 540 GB after flagging and averaging; see Section 2.1). Figure 2 shows an animation of snapshots of one frequency produced by a typical drift scan.
a All observations used three frequency settings centred on 119.04, 149.76 and 180.48 MHz with 30.72 MHz bandwidth.
Unlike Bernardi et al. (Reference Bernardi2013), whose focus was the large-scale polarisation features of the low-frequency sky, these commissioning data are more suitable for measuring the flux densities of individual faint compact sources. Bernardi et al. (Reference Bernardi2013) used the Real Time System (Mitchell et al. Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008; Ord et al. Reference Ord2010) and a forward-modelling scheme to peel 250 sources from 2 400 square degrees of sky; given the instantaneous sensitivity of the commissioning sub-arrays, following such a procedure for our larger, more resolved sky, was infeasible. The following data reduction therefore uses standard astronomy tools such as Common Astronomy Software Applications (casa Footnote 2 ) and miriad (Sault, Teuben, & Wright Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995).
2.1 Flagging and calibration
The MWA band is comprised of 3 072 10 kHz ‘fine’ channels, derived from 24 1.28 MHz ‘coarse’ channels. The eight edge fine channels of each coarse channel suffer aliasing and are flagged in every observation. The central fine channel of each coarse channel contains the (non-zero) DC component of the polyphase filterbank, so is also flagged. Known misbehaving antennas were flagged, and the Orbcomm transmission frequency range 136–138.5 MHz was completely excised. Radio-frequency interference was excised using aoflagger (Offringa, van de Gronde, & Roerdink Reference Offringa, van de Gronde and Roerdink2012); due to the radio-quiet nature of the observatory, less than 3% of the data were affected and removed, mostly at the 136–138 MHz range of the Orbcomm satellite downlink. From a native frequency resolution of 10 kHz, the data were averaged by a factor of four in frequency, appropriately downweighting any frequency bin containing a reduced number of fine channels. From the original time-sampling of 0.5 s, gamma data were averaged by a further factor of two in time, while beta data were time-averaged by a factor of eight, the difference arising from the larger number of long baselines in gamma.
A single calibration solution was determined for each frequency, for each night, based on the calibrator listed in Table 1.
2.1.1 Dec − 27 phase calibration: 3C444
The Dec − 27 drift scans were phase-calibrated using 3C444: RA = 22h 14m 25s.752; Dec = − 17°01′ 36.′′29, which has a flux density 80 Jy and spectral index − 0.9 at 160 MHz (Slee Reference Slee1995). 3C444 dominates the visibilities when observed near the meridian, and is unresolved at the resolution of the commissioning array. A two-minute observation produces solutions of S/N ≃ 5, while integrating over the three two-minute snapshots in which 3C444 was closest to the meridian increases the S/N to ≃ 7. As the sky has drifted by a total of 18 min by the time data have been taken at all three frequencies, the instrumental response is different for each snapshot, so the expected flux of 3C444 was scaled by a primary beam model for each snapshot (see Section 3.5 for more details on the primary beam model). A least-squares fit over the full duration of each observation was performed on those visibilities with u, v-distance > 0.1 kλ using bandpass, a casa routine, to produce a single complex gain solution for each tile, for each polarisation, for each frequency interval, which was then applied to the night’s observations.
2.1.2 Dec − 47 phase calibration: Pictor A
The Dec − 47 drift scans were calibrated using Pictor A: RA = 05h 19m 49s.7; Dec = − 45°46′ 43.′′70. Pictor A was used only as a phase calibrator; its flux density of S ≃ 400 Jy (Slee Reference Slee1995) over the MWA observing band dominates the visibilities such that a single two-minute observation is sufficient to obtain per-channel solutions with S/N ≫ 20. Since Pictor A is resolved by our instrument, the Very Large Array (VLA) images at 1 400 and 333 MHz (Perley, Roser, & Meisenheimer Reference Perley, Roser and Meisenheimer1997) were used to extrapolate a spatial and spectral model for each MWA observing frequencyFootnote 3 . The same procedure was followed as with 3C444: the model was Fourier-transformed and used to generate a calibration solution via bandpass, which was then applied to the night’s observations.
After calibration, the beta and gamma data were concatenated. We note that the ionosphere was different between the two nights, but only the gamma data possesses baselines of sufficient length to experience phase offsets significant compared to the width of the synthesised beam. Therefore, it is possible to combine the beta and gamma visibility data for the same area of sky, but not to combine all of the gamma data together, primary beam effects aside.
2.2 Subtraction of bright sources in primary beam sidelobes
The first sidelobes of the MWA primary beam are estimated to have a response of approximately 3% of that of the main beam, and the second around 0.1% (see Section 3.5 for more details on the primary beam model). When an extremely bright source such as Cygnus A passes through such a sidelobe, it contributes signal to the visibilities. If not subtracted directly from the visibilities (‘peeled’) or deconvolved, this contribution is visible in the image as striping from the sidelobes of the synthesised beam at the position of the source. In the absence of computing constraints, one could attempt to image the entire sky, such that all sources are deconvolved. However, we found that bright sources sufficiently far (> 45°) from the phase centre were poorly deconvolved even in a snapshot image due to the array not being perfectly coplanar. At the time of analysis, techniques such as w-projection (Cornwell, Golap, & Bhatnagar Reference Cornwell, Golap and Bhatnagar2008) were too computationally expensive for the ~ 10 000 pixel-wide grid required to cover the sky. The sidelobes also introduce a large spectral variation to the image pixels, as they vary much more quickly with frequency than the main beam.
To reduce the effects of bright sources in the sidelobes of the primary beam, we phase-rotated the calibrated, concatenated snapshot datasets to the source positions and performed a shallow (threshold = 10 Jy) clean of a 1° × 1° around each source, with an additional Taylor term, allowing each pixel to have its own spectral index (Sault & Wieringa Reference Sault and Wieringa1994). The model pixels were then subtracted from the visibilities, which were then phase-rotated back to their original positions, ready for imaging. The sources removed in this way were Pictor A, the Crab Nebula, Hydra A, Cygnus A, PKS2153 − 69 and PKS2356 − 61. Only Pictor A was present in the primary beam main lobe, as well as the sidelobe, and we note that this strategy means we are then unable to recover its flux density in the final mosaic.
2.3 Imaging strategy
At 150 MHz, the full-width half maximum (FWHM) of the MWA primary beam at zenith is approximately 30°, corresponding to more than two hours of RA. Sources therefore appear in many of the two minute scans as the sky drifts through the beam. Within any two-minute scan (hereafter a ‘snapshot’), the hybrid beta/gamma array is sufficiently coplanar that snapshot images do not suffer from wide-field effects within the main primary beam lobe (< 15° from the phase centre) and the image coordinates can be described by a slant-orthographic (generalised SIN) projection (Perley Reference Perley, Taylor, Carilli and Perley1999; Calabretta & Greisen Reference Calabretta and Greisen2002). The data were thus divided into snapshots with the phase centre defined as the local sidereal time in the middle of the snapshot. The goal of data processing is to form an image of the entire RA range of the observation using all of the available data.
Given the drift scan observation strategy, standard synthesis imaging techniques are not suitable. In a drift scan, instead of pointing the beam at adjacent locations on the sky, the beam stays constant and the sky moves, relatively. The data are therefore best processed using a mosaicking methodology. Typically the goal of a mosaic observation is to image objects larger than the primary beam of the telescope (e.g. Holdaway Reference Holdaway, Taylor, Carilli and Perley1999). Much of the mosaicking literature and software is therefore focused on recovering the large-scale structure in a synthesis image using an interferometer with very short baselines. In the case of these data, the commissioning sub-array has only 12 baselines shorter than 20 m and thus sensitive to scales > 5°, unlike the full MWA which fully samples up to scales of 12° (at 180 MHz). As we are also observing a part of the sky mostly devoid of objects of large angular scale, our focus is on combining the data to maximise sensitivity to compact structures.
To form the mosaic images, we perform the following steps, for each frequency band, for each polarisation:
• generate an output accumulation image for the data and primary beam in an equal area coordinate system with 0.5 arcmin pixel resolution and over a sufficient region of the sky for all data;
• for each snapshot:
– invert using multi-frequency synthesis over a field-of-view of 40° × 40°, with a pixel resolution of 1 arcmin, weighting each u, v-cell equally (‘uniform weighting’);
– clean to a threshold of three times the typical snapshot root-mean-square (RMS), which itself is typically ≃ 240–80 mJy beam− 1 from 120–180 MHz; the synthesised beam sidelobes are only ≈ 10% of the peak response, and a gain factor of 0.1 is used, thus avoiding cleaning the noise, and clean bias;
– restore using an elliptical Gaussian approximation of the snapshot synthesised beam;
– regrid and add the deconvolved image multiplied by the appropriate primary beam model into the accumulation data image;
– likewise, regrid and add the square of the beam image into the accumulation beam image;
• after all additions, divide the data image by the beam-squared image to form the final mosaic.
In essence, this follows Equation (1) from Sault, Staveley-Smith, & Brouw (Reference Sault, Staveley-Smith and Brouw1996) which maximises the signal-to-noise ratio (S/N) in the output mosaic given the changing S/N over the field in each snapshot due to the primary beam. We use the miriad for the imaging, deconvolution and image arithmetic, and use our own code to generate MWA primary beam models (see Section 3.5 for more details on the primary beam model). miriad’s understanding of advanced World Coordinate Systems (WCS) and ability to regrid between arbitrary coordinate frames is especially useful.
The process of adding snapshots improves the quality of the final image by reducing the thermal noise and improving the synthesised beam, revealing many fainter sources. These fainter sources have not been deconvolved in the snapshots, so their sidelobes remain in the mosaic with ~ 1 h of effective earth rotation synthesis (the time it takes for sources to move through the utilised area of the primary beam). We estimate the (5σ) classical confusion limit (Condon Reference Condon1974) to be approximately 5 mJy beam− 1 given our effective synthesised beam size; this is consistent with the confusion measurement by Bernardi et al. (Reference Bernardi2009, Reference Bernardi2010) at 150 MHz at a similar angular resolution and frequency. While we therefore begin to approach the classical confusion limit, the dominant components of noise in our mosaics are sidelobe confusion and calibration errors. A randomly-chosen section of the 180 MHz Dec − 47 mosaic is shown in Figure 3, showing the well-behaved point spread function (PSF) and averaging down of calibration errors that results from the mosaicking process.
The zenith-pointed MWA primary beam is virtually identical for the XX and YY polarisations. We therefore averaged the XX and YY Dec − 27 mosaics to form a (pseudo) Stokes I image. For the Dec − 47 mosaic, we found that the primary beam response of the XX and YY polarisations was sufficiently different that Dec-dependent corrections were required. Hence, the XX and YY mosaics for Dec − 47 were not combined at this stage. Details of this process and absolute calibration are discussed in Section 3.5.
3 SOURCE-FINDING AND FLUX DENSITY ESTIMATION
The imaging process results in nine mosaics, consisting of three frequencies for the Dec − 27 scan, and three frequencies × two polarisations for the Dec − 47 scan. The RMS noise level of these mosaics is typically a factor of five lower than that of the individual snapshots, allowing much deeper source-finding; see Section 4.1.1 for a discussion of the origins of noise in the mosaics, and their effect on the accuracy of our source flux densities.
Maps of the RMS value of the 180 MHz mosaic pixels were generated using aegean (Hancock et al. Reference Hancock, Murphy, Gaensler, Hopkins and Curran2012) which calculates the RMS for blocks of pixels with dimensions of approximately 20 × 20 beamsFootnote 4 . These are reproduced for two representative mosaics in Figure 4 (see also Figure 16). Two effects can clearly be seen in these figures: an increase in RMS at the edges of the mosaics, and isolated patches of locally high RMS. The first is due to the effect of the primary beam and the second is due to the effect of sidelobe confusion and calibration errors around the brighter sources. The latter is particularly severe near extremely bright, resolved structures which can most easily be seen close to the Galactic plane (at the easternmost edge of the Dec − 47 scan). Fornax A produces more calibration-error noise in the Dec − 27 scan, and less in the better-calibrated Dec − 47 scan; deconvolution errors from poorly reconstructing the largest scales of its structure are present in both.
3.2 Source identification
We wish to identify and characterise the morphology and flux density of all sources within the field which are bright enough to be discerned from the background noise. Having characterised the noise properties of the mosaics, we identify those pixels of the mosaics whose flux is greater than five times the local RMS, which itself is hereafter denoted σ. blobcat (Hales et al. Reference Hales, Murphy, Curran, Middelberg, Gaensler and Norris2012) was used to identify ‘blobs’: islands of pixels containing a 5σ peak, expanded to include all contiguous pixels brighter than 3.5σ, while accounting for various biases.
This was repeated for all three (Dec − 27, Dec − 47 XX and Dec − 47 YY) of the highest frequency (180 MHz) mosaics. These were used exclusively for source-finding since the highest resolution allowed the greatest separation of neighbouring sources into discrete blobs. The objects detected in the three maps were combined into a single list (i.e. from this point on no regard was taken as to which map the source was detected in).
3.3 Characterising the point spread function
For an interferometer coplanar to the plane perpendicular to the zenith, observing at zero hour angle, the synthesised beam becomes elongated N-S relative to the zenith, as csc ζ, where ζ is the zenith angle. Regridding from a slant orthographic projection to a zenithal equal area projection (using standard miriad WCS tools) conserves the synthesised beam volume, ignoring this effect, such that a correction factor of csc ζ must be applied to rescale the N-S extent of the synthesised beam. We apply this immediately, before making further measurements.
While the synthesised beam shape for an individual snapshot is known very accurately (short-timescale ionospheric effects aside), determination of the effective synthesised beam in the mosaics is a much more difficult problem. The mosaic contains approximately an hour (depending on frequency) of effective earth rotation synthesis which significantly circularises the PSF. However, refraction from the ionosphere acts differently in each snapshot and enlarges the PSF as snapshots are combined to form each mosaic.
To measure this effect, we fitted 2D Gaussians to all detected sources, and plotted the distribution of ellipse values against various parameters (such as RA, Dec) and checked for correlations: none were seen. Plotting against S/N, as in Figure 5, shows that the ellipse sizes of unresolved sources form a line, with the resolved sources having larger values. For each Gaussian parameter, a line was fitted to the unresolved sources, and compared with the predicted value in order to calculate the correction needed to increase the size of the PSF to match the data.
Henceforth, all source-finding uses the corrected PSFs, shown in Table 2. Note that the PSFs are identical for the XX and YY polarisation Dec − 47 mosaics, as all of these effects are polarisation independent. A further projection factor of csc ζ is applied when measuring sources in the mosaics. Figure 6 shows an exaggerated illustration of the effect of these changes. After these corrections, it is possible to recover the flux densities of extended sources, which would otherwise be over- or under-estimated depending on the difference between the local PSF and average PSFFootnote 5 .
3.4 Characterisation of discrete sources
We now aim to determine flux density and morphology as accurately as possible given the three (Dec − 27), six (Dec − 47) or nine (overlapping region) measurements that we have for each source. On visual inspection it is clear that the vast majority of our objects can be characterised by a single elliptical Gaussian similar in shape to the synthesised beam. This reflects the fact that the vast majority of our sources are unresolved, as expected.
In order that they be fit with the minimum number of free parameters possible, we begin by modelling each of our sources as a single elliptical Gaussian, only attempting to characterise more complex morphology when this assumption is determined to be invalid (see Section 3.4.1). In another refinement to avoid the effects of confusion due to nearby sources, the fits were performed only over those pixels which lie within the FWHM of the a priori fit (with an extra margin of 10%). Levenburg–Macquart least-squares fits were made to these pixels, with free parameters: amplitude; major axis, a; minor axis, b; position angle, θ; RA offset, ΔRA; and Dec offset, ΔDec. The initial parameters were the synthesised beam (with frequency and position dependence calculated as described in Section 2.3) with an amplitude of the peak pixel in the region of the fit, centred on the flux weighted centroid of the object.
For error calculation purposes the RMS in the region of the blob in each map was determined by taking the mean of the RMS map over an array of 24 × 24 pixels centred on the centroid of the blob. This provides a modest amount of smoothing of the RMS shown in Figure 4 (where the RMS is determined over blocks a factor of ~ 5 × 5 larger).
3.4.1 Separating single-component and multi-component sources
Next we determine which of the blobs are indeed well-characterised by a single elliptical Gaussian and which require further characterisation. Our criterion is that any source whose centroid lies within the half-maximum of the a priori elliptical Gaussian is ‘single-component’. This simple criterion is effective due to our selection of the a priori centroid on the basis of the amplitude-weighted centroid of the blob.
When there is no peak (i.e. positive curvature in both spatial dimensions) within our stringent pixel range, the centroid of the fit is forced outside the a priori ellipse. If there is a peak, but it is not coincident with the a priori centroid, this is indicative of signal within the blob which is not associated with the fitted Gaussian.
Using this method, 12 863 (≈ 91 %) of our 14 110 sources are classified as single-component.
3.5 Flux density calibration
The southern sky at low frequencies is poorly surveyed and flux densities derived from extrapolation of various radio catalogues to the MWA frequency range disagree at the 10–20% level. Work is under way in the community to develop a unified radio flux density scale from MHz–GHz frequenciesFootnote 6 , and GLEAM will expand this scale over the entire southern sky. At the time of writing, the best course of action was to bootstrap from the fairly well-understood northern sky to the south. Conveniently these declinations pass close to the zenith for the MWA, allowing the use of a simple model of the MWA primary beam.
The MWA primary beam has previously been approximated by an analytic model incorporating a dipole over groundscreen and geometric array factor with good results (e.g. Williams et al. Reference Williams2012; McKinley et al. Reference McKinley2013; Bell et al. Reference Bell2014), and we adopt the same approach in this paper. Bernardi et al. (Reference Bernardi2013) find that the zenith primary beam model is incorrect at the ~ 2% level at higher frequencies, and we expect the model to be less correct further from zenith as the simple dipole approximation becomes less valid. Due to the uncertainty in the primary beam model, especially away from the zenith, it is impossible to correctly perform the absolute flux density calibration in a single step. In addition, phase errors from imperfect phase calibration will reduce the recovered peak flux of radio sources, necessitating a small, position-independent correction. We note that understanding and improving the primary beam model is an ongoing activity within the MWA collaboration. Recent work (Sutinjo et al. Reference Sutinjo, O’Sullivan, Lenc, Wayth, Padhi, Hall and Tingay2014) to incorporate mutual coupling effects into the model has improved our understanding of the response of the MWA beam. However, the work is still ongoing and for the purposes of this paper our empirical method to bootstrap the flux density scale, described below, was sufficient.
In order to provide an absolute flux density scale, a list of sources was drawn up which appear in our catalogue as well as all of the Culgoora (80, 160 MHz Slee Reference Slee1995), VLSS (74 MHz Cohen et al. Reference Cohen, Lane, Cotton, Kassim, Lazio, Perley, Condon and Erickson2007), WISH (325 MHz De Breuck et al. Reference De Breuck, Tang, de Bruyn, Röttgering and van Breugel2002), MRC (408 MHz Large et al. Reference Large, Mills, Little, Crawford and Sutton1981), and NVSS (1.4 GHz Condon et al. Reference Condon, Cotton, Greisen, Yin, Perley, Taylor and Broderick1998) catalogues. A good choice of flux calibrator might be our northern phase calibrator, 3C444; however, there are two reasons to avoid using this source. Firstly, it is towards the edge of our RA range, and is only sampled by a few snapshots. Secondly, it is resolved into two components by NVSS, so measurements by different instruments with different resolutions may give conflicting results.
The second brightest candidate is 3C32. This source as it appears in NVSS (see Figure 7) would be unresolved in our data; its flux density is well-modelled by a simple power law (Figure 8) and shows no evidence of variability or being resolved. In particular, the VLSS 74 MHz and Culgoora 80 MHz flux densities both lie on the power-law fit within the errors, despite their differing resolutions and epochs of observation. We performed an image-plane Gaussian fit (Figure 9) to measure the flux density of the source and confirmed it as unresolved. We therefore scaled the integrated flux densities of each of the Dec − 27 maps appropriately to fit the least-squares fit to the catalogue fluxes (see Table 3). Since the zenith analytic primary beam model had already been divided out in creating the Dec − 27 mosaic, the relative flux densities of sources were already correct and only this single scaling factor should be required to set the absolute flux density scale across the whole mosaic, assuming that the zenith primary beam model is correct.
The flux densities measured in the XX and YY mosaics of the Dec − 47 field disagree at the ~ 5% level due to unmodelled primary beam effects. A second-order polynomial fit (weighted by the two flux densities added in quadrature) was made to the Dec-dependent discrepancy between XX and YY as shown in Figure 10. The XX mosaic was then scaled to match the YY mosaic.
The Dec − 47 mosaic was corrected to be consistent with the Dec − 27 mosaic by using the ≈ 600 unresolved sources which lay in the overlapping region. The mean of the flux density ratio between the two maps of all of these sources was taken to be the additional correction factor to the Dec − 47 mosaic to produce a corrected absolute flux density calibration scale consistent with that used for Dec − 27 (Table 3). After this correction, the XX and YY mosaics of the Dec − 47 field were combined, weighted by their respective primary beam responses (which differ by around 5%), to form a pseudo-Stokes I mosaic, on which source-finding can be performed.
3.6 Combining the mosaics
Many MWACS sources are detected both in the Dec − 27 and Dec − 47 mosaics. Combining the information from both mosaics is likely to give superior results for a number of reasons. This would recover some of the signal-to-noise in the overlapping region of the two mosaics (at Dec ≃ −36°) which is otherwise attenuated by the primary beam, especially at the highest frequency. The differing primary beam between the two scans also means that the signal from sources in the distant sidelobes of the primary beam will be somewhat different, so combining the mosaics would reduce this noise.
Unfortunately, the very different state of the ionosphere between the two gamma drift scans makes a simple image-plane combination difficult; without extensive modelling of the ionosphere, there would be large position- and frequency-dependent smearing of the sources. This is evident from Figure 14; the astrometry errors are time-dependent (seen as a change in RA) and differ night-to-night. However, an improved flux density estimate can be determined for each source in the overlapping region by taking the mean of the two peak flux densities, weighted by the S/N of the fit in each of the maps (as defined in Section 3.8).
Though the size of the overlapping region varies dramatically with frequency, for simplicity and transparency the map flux densities were only combined in the region well-sampled by both the Dec − 27 and Dec − 47 scans at the highest frequency: a Dec range of − 38.°20 < δ < −35.°25. Furthermore, the positions and morphologies of the sources in the overlap region were taken as those measured in the Dec − 47 mosaic, as this scan has a somewhat more compact synthesised beam and marginally better phase calibration (see Section 4.1.1).
3.7 Fitting 180 MHz flux densities
The majority of the MWACS source-flux densities are well-characterised by a simple power-law spectral index; as described in Section 3.2, we perform source detection at the highest resolution available, at 180 MHz, so we use the fitted integrated source flux densities at each frequency to measure the spectral indices. Alongside this, we report the 180 MHz flux density for each source.
We performed weighted least-squares fits to the flux densities at the three frequencies, propagating the confusion and fitting errors described in Section 3.8.
3.8 Estimation of errors
3.8.1 Noise-like errors
Condon (Reference Condon1997) addresses the problem of determining the error on the parameters of a least-squares fit of an elliptical Gaussian G(xk, yk ) to a set of pixel values ak . Here, we adopt the same notation where A, θ M , θ m , ϕ x 0 and y 0 are the amplitude, FWHM major and minor axes, position angle and centroid offsets, respectively, of a Gaussian fit over ak pixel amplitudes, each having the same Gaussian error distribution of RMS μ. In our case, the noise, whether it be thermal or confusion noise, will be correlated on the scale of the synthesised beam. For the case of unresolved or weakly resolved sources, where the correlation of the noise matches the apparent source size, the error on A 2, μ2(A) is μ2, i.e. the error on the flux density of the source due to noise is simply the local RMS. The errors on A, θ M and θ m are therefore:
By analogy to Condon (Reference Condon1997) Equation 21,
where μ(ϕ) > 90°, μ(ϕ) is set to null.
Finally, the errors on RA (α) and Dec (δ) due to noise-like errors are given by
3.8.2 Calibration errors
In the bright source regime, the error on the flux density is dominated not by noise but by calibration errors. This is added in quadrature to the other sources of error as a fractional error of 4% on the integrated flux density. This was determined by measuring the local increase in RMS flux density around bright sources of known flux density, where calibration errors are the dominant source of noise, and dividing it by the flux density of each bright source.
3.9 Remaining systematic errors
3.9.1 Error on the absolute calibration
To test the consistency of the absolute flux density calibration across the field and across different sources, we first cross-matched our catalogue with the only other catalogue which lies within our frequency range: the Culgoora-3 list of radio source measurements at 160 MHz, at a resolution of 1.85 arcmin. These flux densities were scaled to 180 MHz using the MWACS spectral index and the scaled flux densities compared with the MWACS flux densities, shown in Figure 11. While the scatter in the ratio of flux densities is large, and considerably larger than the estimated error on the Culgoora fluxes (Slee Reference Slee1977, Tables 1 and 2), there is no evidence of a systematic offset with respect to declination, which would have indicated uncorrected primary beam errors.
In addition, we identified all 41 sources which appeared in all the five catalogues used for absolute calibration (Section 3.5) and which had an MRC (408 MHz) flux density greater than 4 Jy (hereafter the MRC4 sample). Since these catalogues only cover a restricted Dec range near the top of our field, these sources were supplemented by 73 sources from the ‘MS4’ sample, defined by Burgess & Hunstead (Reference Burgess and Hunstead2006) as those sources in the MRC catalogue with a flux density greater than 4 JyFootnote 7 .
From the resulting sample we excluded all sources which were resolved in our catalogue—defined as those with a fitted ellipse > 5% larger than synthesised beam ellipse (19/73 MRC4 sources and 9/41 MS4 sources were excluded). Two fits were made to the catalogued flux densities: one a simple power law (i.e. a straight-line fit in log–log space); another incorporating curvature (a parabolic fit in log–log space). Those sources where the 180 MHz flux density predicted by the two fits were discrepant by more than 10% were classified as ‘curved’. For the MS4 sources we used the multi-frequency flux density measurements compiled by Burgess & Hunstead (Reference Burgess and Hunstead2006). Many of these sources had not listed flux density below 300 MHz (i.e. the 180 MHz flux density is extrapolated from flux density measurements at 365 or 408 MHz and above). These sources were classified as ‘extrapolated’. The resulting flux density ratios are also shown in Figure 11, which overall illustrates the challenge of defining a low-frequency flux scale in the southern hemisphere. There are few sources with well-characterised spectra at this frequency, and there are inconsistencies between these sources.
Figure 12 replots these data as cumulative histograms of the difference of the flux density ratio from unity. Three subsets of (unresolved) sources are shown: Culgoora, MS4 & MRC4, and MS4 & MRC4 sources which do not have curved spectra and are not pure extrapolations from 408 MHz downward. Using this latter, most high-quality subset, we see that 68% of sources lie within 10% of unity. We therefore add an absolute flux calibration error of 10% in quadrature to the final flux density error after fitting.
3.9.2 Primary-beam errors
For a mosaicked drift scan, an error in the primary beam model would be expected to manifest itself as a Dec-dependent error in flux density. Overall there is little evidence of a Dec-dependent error in Figure 11.
Any error in the primary beam model would very likely have strong frequency-dependence, such that it may systematically affect the spectral index. This can be tested by binning the sources by spectral index and looking for systematic changes with Dec. The results of this analysis are shown in Figure 13.
There are clearly systematic changes in the median spectral index with Dec, though there is no obvious functional form which could be modelled to take account of this. However the change in spectral index across the band corresponds to an error across the frequency range of less than 10%. An error of 0.1 is added in quadrature to the to the error on the fitted spectral index to take account of the systematic shift in spectral index seen here.
3.9.3 Astrometry errors
Figure 14 shows the astrometric error for all unresolved sources cross-matched between our catalogue and the NVSS catalogue for the Dec − 27 mosaic, and the SUMSS catalogue for the Dec − 47 mosaic. The astrometric errors in RA and Dec are not entirely correlated, vary with RA (time, for a drift scan), and are likely to be functions of ionospheric activity and viewing angle through the ionosphere (i.e. zenith angle). As we do not attempt to solve this systematic error, we measure the RMS without weighting each measurement by the flux density of the sources, which leads to a conservative astrometric RMS estimate of 0.3 arcmin. This value is added in quadrature to the errors on the RA and Dec derived in Section 3.8.1.
3.10 Resolved sources
While the majority of sources in the catalogue are unresolved at the resolution of the commissioning array, and are therefore described well by a single flux density, some objects are more complex and must be examined more closely to determine an integrated flux density. These are automatically detected as described in Section 3.4.1.
For each of these sources, a 100 × 100 pixel subsection of the mosaic was regridded to a local SIN projection, and the appropriate PSF determined in Section 2.3 was added to the header. The source-finding software aegean was then run on the sub-image, detecting and fitting multiple components simultaneously. The results of this fitting at the highest frequency only were used to produce the integrated flux density measurement at 180 MHz.
Due to the varying resolution over the bandwidth of the array, cross-matching from the 180 MHz-detected sources to the other two frequencies is not always straightforward when sources are extended, or close to each other.
• For some sources, particularly those which are only slightly extended and not confused with other nearby sources, the components match easily, and integrated flux densities are calculated for each freqency.
• For confused sources, where different numbers of components are detected at different frequencies, an integrated flux density is calculated for each frequency. This is done by flood-filling an aegean-detected island down to a level four times the local RMS flux density, and integrating.
• For those (≈ 5σ) sources which were not easily fit by a Gaussian, or by a flood-fill, at 120 and 150 MHz, the position centroid from the source fit to the 180 MHz data was used to perform a ‘forced’ measurement on that pixel position in the lower-frequency data.
The flux densities at each frequency are then fitted in the usual way to produce a single spectral index for the 180 MHz source components.
Examples of each method are shown in Figure 15. The fitted spectral indices and 180 MHz flux densities for the components of resolved sources are reported in the catalogue in the same way as for the unresolved sources. As shown in Table 4, a column in the table indicates whether the spectral index was determined by integrating an extended source, or performing a forced measurement.
Error measurements for the flux densities of these extended sources were identical to those discussed in Section 3.8, with the fitting error combined in quadrature with the RMS noise. It should be noted, however, that the changing resolution of the instrument at different frequencies, combined with the intrinsic complexity of these sources, mean that the automatically derived spectral indices should be treated with caution. Postage stamps of the extended sources are available on request should the reader wish to perform their own measurements.
4 SOURCE CATALOGUE
The MWACS catalogue comprises 14 110 sources and is available via VizierFootnote 8 , and in the electronic edition of this journal. We report the positions of sources, the 180 MHz integrated flux density and spectral index, the Gaussian fit parameters (major axis, minor axis, and position angle) as well as errors on all of the above parameters. Flags indicate whether a source was confused or extended and thus required special treatment (see Section 3.10.) The estimated synthesised beam at that position is also provided. Postage stamps are available on request.
4.1 Catalogue reliability and completeness
Recall that all of the source-finding for this survey was performed on the highest-frequency mosaic. For point sources, after source identification, the source was fitted at all other frequencies. If any of these fits failed, the source would be treated as resolved; see Section 3.10 for details. We searched for and identified ten spurious sources in the sidelobes of bright sources by searching the areas around (≈ 100) sources of S > 10 Jy, comparing images at all frequencies, and removed these from the catalogue; the faintest source which produced a spurious sidelobe source was 14 Jy. In typical circumstances away from bright sources, we require that all sources are detectable at 5σ at 180 MHz and well-fitted by a Gaussian in the other bands, so with a Gaussian random background, the chance of any source being spurious is vanishingly small. Our noise is somewhat non-Gaussian due to the presence of undeconvolved source sidelobes, but in the general field, these sidelobes are faint due to our good u, v-coverage, which suppresses sidelobes down to a 10% level even before deconvolution; and the positions of change sidelobe alignments would vary with frequency and thus be unlikely to be detected in every channel. Therefore, we can have an extremely high confidence that all of the sources in the catalogue are real.
The varying detection threshold affects the catalogue completeness for two reasons. More obviously, the variations in the local RMS cause a varying absolute detection limit across the field. In addition, the noise in the map introduces a non-zero probability that a source whose intrinsic flux is > 5σ, has a measured flux < 5σ due to instrumental effects (and vice versa). We therefore present the noise properties of our maps in some detail.
4.1.1 Variations in RMS sensitivity
Figure 4 shows maps of the local RMS of the mosaics rescaled to take into account absolute calibration (with XX and YY combined together for the Dec − 47 scan). This information is summarised in Figure 16.
This figure shows that the extremely high RMS regions associated with bright extended sources, particularly near the Galactic plane, account for a small proportion of survey area. Overall, the sensitivity of the survey varies about the median of ~ 40 mJy by less than a factor of two over 80% of the survey area, a variation in sensitivity comparable with the VLSS (cf. Cohen et al. Reference Cohen, Lane, Cotton, Kassim, Lazio, Perley, Condon and Erickson2007, Figure 7).
The RMS in the Dec − 27 and Dec − 47 mosaics are reasonably comparable, though the Dec − 47 mosaic has a somewhat lower RMS, particularly for regions below the median. We ascribe this to a better phase calibration for the Dec − 47 data, as a brighter phase calibrator was used. The RMS becomes only marginally worse towards the edge of the primary beam in the Dec − 47 mosaic, however it becomes significantly worse towards the edge of the primary beam in the Dec − 27 mosaic. This is because the sky rotates more quickly through the primary beam closer to the celestial equator.
Figure 17 shows the histogram of flux density values across a significant proportion of each mosaic. When calculated over a wide area, the RMS for the two maps agrees to ~ 10%; 90% of the surveyed area has an RMS lower than 50 mJy, except for the very northern edge.
4.2 Fornax A
Fornax A is a nearby (z ≈ 0.006) radio galaxy, situated in the Fornax cluster of galaxies. It comprises two large lobes extending NNW and SSE 20 arcmin (≈ 200 kpc) from the lenticular galaxy NGC 1316, at RA = 03h 22m 41s.789 Dec = − 37°12′ 29.′′52. On scales ≳ 20′ , structure is not sampled by MWACS, so we cannot measure the integrated flux density for this object, only set lower limits. This under-sampling also causes deconvolution errors across our mosaics, resulting in increased noise levels around the Fornax A region.
Using blobcat to flood-fill the entire region down to a 4σ level and measure the integrated flux density, we set lower limits across the entirety of Fornax A at 786 ± 39, 668 ± 33, and 514 ± 26 Jy at 120, 150, and 180 MHz, respectively; these are consistent with the values measured by Bernardi et al. (Reference Bernardi2013) and McKinley et al. (Reference McKinley2014). An image of the radio galaxy at the highest MWACS resolution is shown in Figure 18.
5 DISCUSSION AND CONCLUSIONS
We have performed a survey of approximately 6 100 square degrees of the southern sky during the MWA’s commissioning period in 2012 October: the MWA Commissioning Survey (MWACS). The observations were made in meridian drift scan mode at two Dec settings (− 26.°7 and − 47.°5) centred on three frequencies: 119, 150, and 180 MHz, resulting in an effective frequency coverage of 104 to 196 MHz. Because the data were taken as drift scans, and because the MWA’s primary beam X and Y polarisation responses are different when off-zenith, some unconventional data processing was required. Six separate mosaic images were formed, one for each frequency and Dec covering the entire RA range of the drift scan. By taking advantage of the instantaneous (near) co-planarity of the array and good instantaneous u, v-coverage (including large fractional bandwidth), the data were processed as continuum multi-frequency synthesis images (at the three central frequencies) by forming a mosaic of weighted regridded snapshot images.
The source 3C32 was used to set the absolute flux density scale. Stokes I images at 119, 150, and 180 MHz were formed by adding the primary beam weighted XX and YY mosaics. The absolute flux density was bootstrapped from the Dec − 27 mosaic to the Dec − 47 mosaic using sources in the overlapping five degrees of declination. We found the primary beam response at Dec − 47 to be sufficiently different between the XX and YY polarisations that an independent scaling factor was used for each before being combined to Stokes I.
The data do not have sufficient short baselines to recover large-scale structure, so only compact sources were extracted from the mosaics. We first identified sources as islands of > 5σ flux density, then separated sources into single-component (unresolved) and multi-component (resolved, and/or multiple nearby components) categories, and measured their flux densities at each frequency, appropriately. We present the result of the MWACS a single catalogue comprised of 180-MHz flux densities and spectral indices calculated from the three individual measured frequencies, for 14 110 sources.
The sensitivity of MWACS, as measured by local image RMS, is mostly consistent around 40 mJy beam− 1 (1σ) across the mosaics but degrades near bright sources and near the edges of the fields. All sources found in MWACS at 180 MHz were also found in the two lower frequency mosaics and had reasonable spectra given their S/N. We therefore consider detections to be reliable subject to the usual considerations on confusion and resolution. Users are advised, however, that non-detections can only be used to set an upper limit to flux density in the context of the local RMS. Regions within a few degrees of the edge of the image or bright sources will have locally higher RMS.
This commissioning survey demonstrates the impressive survey capabilities of the MWA and the feasibility of processing MWA drift scan datasets using fairly conventional radio astronomy tools with simple mosaicking techniques. The full MWA has near-complete snapshot u, v-coverage between approximately 5 and 500 wavelengths (as multi-frequency synthesis images) and with 128 fully cross-correlated tiles, its quality of calibration is higher than our combined commissioning sub-arrays. We therefore expect the GLEAM survey to be confusion limited in Stokes I, to be able to capture large-scale structures (up to ~ 10°) and to have much reduced imaging artifacts from better calibration and deconvolution.
The progression of the MWA from 32T prototype (circa 2009–2011) through the commissioning array (late 2012 through mid-2013) and into operations (mid-2013 on) has shown a steady improvement in area and depth of sky coverage and of our understanding of the instrument and data. The GLEAM survey, currently underway, has derived great benefit from the accumulated expertise of the 32T and commissioning phases and will likewise provide increased understanding of the instrument, especially the primary beam and resulting polarisation performance of the MWA, and data processing methods. Looking forward to SKA-low, we expect the expertise gained from processing the large GLEAM dataset on world-class supercomputers to be provide equally valuable lessons when SKA-low Phase 1 is complete.
This scientific work makes use of the Murchison Radio-astronomy Observatory, operated by CSIRO. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site. Support for the MWA comes from the U.S. National Science Foundation (grants AST-0457585, PHY-0835713, CAREER-0847753, and AST-0908884), the Australian Research Council (LIEF grants LE0775621 and LE0882938), the U.S. Air Force Office of Scientific Research (grant FA9550-0510247), and the Centre for All-sky Astrophysics (an Australian Research Council Centre of Excellence funded by grant CE110001020). Support is also provided by the Smithsonian Astrophysical Observatory, the MIT School of Science, the Raman Research Institute, the Australian National University, and the Victoria University of Wellington (via grant MED-E1799 from the New Zealand Ministry of Economic Development and an IBM Shared University Research Grant). The Australian Federal government provides additional support via the Commonwealth Scientific and Industrial Research Organisation (CSIRO), National Collaborative Research Infrastructure Strategy, Education Investment Fund, and the Australia India Strategic Research Fund, and Astronomy Australia Limited, under contract to Curtin University. We acknowledge the iVEC Petabyte Data Store, the Initiative in Innovative Computing and the CUDA Center for Excellence sponsored by NVIDIA at Harvard University, and the International Centre for Radio Astronomy Research (ICRAR), a Joint Venture of Curtin University and The University of Western Australia, funded by the Western Australian State government. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.