Calibration and Stokes Imaging with Full Embedded Element Primary Beam Model for the Murchison Widefield Array

Abstract

The Murchison Widefield Array (MWA), located in Western Australia, is one of the low-frequency precursors of the international Square Kilometre Array (SKA) project. In addition to pursuing its own ambitious science programme, it is also a testbed for wide range of future SKA activities ranging from hardware, software to data analysis. The key science programmes for the MWA and SKA require very high dynamic ranges, which challenges calibration and imaging systems. Correct calibration of the instrument and accurate measurements of source flux densities and polarisations require precise characterisation of the telescope’s primary beam. Recent results from the MWA GaLactic Extragalactic All-sky Murchison Widefield Array (GLEAM) survey show that the previously implemented Average Embedded Element (AEE) model still leaves residual polarisations errors of up to 10–20% in Stokes Q. We present a new simulation-based Full Embedded Element (FEE) model which is the most rigorous realisation yet of the MWA’s primary beam model. It enables efficient calculation of the MWA beam response in arbitrary directions without necessity of spatial interpolation. In the new model, every dipole in the MWA tile (4 × 4 bow-tie dipoles) is simulated separately, taking into account all mutual coupling, ground screen, and soil effects, and therefore accounts for the different properties of the individual dipoles within a tile. We have applied the FEE beam model to GLEAM observations at 200–231 MHz and used false Stokes parameter leakage as a metric to compare the models. We have determined that the FEE model reduced the magnitude and declination-dependent behaviour of false polarisation in Stokes Q and V while retaining low levels of false polarisation in Stokes U.

1 INTRODUCTION

Accurate astronomical measurements with radio interferometric telescopes require correction of instrumental effects. For fixed-antenna aperture array telescopes, the primary beam of the interferometer elements varies considerably with pointing direction. It is important to accurately model and correct for the primary beam, as differences between the actual and modelled beam result in errors during both telescope calibration and imaging. Moreover, correct calibration of the primary beam effects (i.e. beam chromaticity and polarisation leakage) in low-frequency aperture arrays is critical for detection of the Epoch of Reionisation (EoR), which is a key science goal of the Square Kilometre Array (SKA) ¹ and its precursors. Asad et al. (2016) studied beam effects and their impact on the EoR science for the Low Frequency Array (LOFAR) and analysis for the Hydrogen Epoch of Reionization Array (HERA) can be found in DeBoer et al. (2017). Here we present the new primary beam model for the Murchison Widefield Array (MWA), which is a low-frequency (~75–300 MHz) telescope and SKA precursor that commenced scientific operation in 2013 (Tingay et al. 2013). The challenges of low-frequency radio polarimetry, based on experiences from the MWA, were summarised in the recent paper by Lenc et al. (2017).

The Phase I deployment of the MWA consists of 128 ‘tiles’ with separations up to 3 km, each tile being a small electronically steerable phased array of 16 dual-polarised bow-tie antennas (Figure 1). The steering is provided by beamformer units where appropriate 5-bit time delays (in discrete steps of 435 ps) are applied to each of the 16 antennas in the tile. The signals from each tile are digitised and cross-correlated with other tiles, and the visibility measurements archived. The archived data covers 24 coarse channels of 1.28 MHz, which are divided into fine channels of 10 or 40 kHz—depending on the correlator settings (Ord et al. 2015).

Figure 1. Examples of the MWA’s aperture array antenna ‘tiles’, each comprised of a 4 × 4 grid of individual bow-tie dipoles (Image credit: MWA Project, Curtin University).

With processing pipelines now in place for key science observing with the MWA, such as detecting emission from the EoR (Jacobs et al. 2016) and the GaLactic Extragalactic All-sky MWA (GLEAM) survey (Wayth et al. 2015; Hurley-Walker et al. 2017), focus has turned to the accuracy of the primary beam models for the MWA.

An incorrect beam model manifests itself during calibration and imaging of the target field. A simple observing scenario is where the tile-based complex gain calibration solutions are determined by observing a strong, point-like source. Ideally, these calibration solutions are corrected for direction-dependent effects (DDEs, i.e., the modelled gain of the tile beam in the direction of the calibrator source), resulting in measurement of the direction-independent electronic gain for each tile. These calibration solutions can then be applied as a correction to the target visibility data collected on a field with or without a suitable calibrator source. In addition, the DDEs in the target pointing direction must again be corrected for. Errors in the beam model propagate through as errors in the DDE-corrected observation. The level of error will depend on the pointing direction of the tile beam during calibration and observation, and the location of the sources within the beam. The most sensitive probe of model inaccuracies are polarimetric measurements of the celestial sources because errors in the beam model for linearly polarised receptors, like the MWA’s tiles, manifest themselves as false Stokes in the calibrated images (Lenc et al. 2016).

In response to false Stokes Q observed in the data calibrated with a simple analytical (Hertzian dipole) tile beam model, a new, feko ² software simulation-based model was implemented using an ‘average embedded element’ (AEE) pattern (Sutinjo et al. 2015a). The AEE model showed significant improvements with respect to analytical model and reduced false Stokes Q in the calibrated data from ≈30% to typically below 10%. The AEE model was used to calibrate the GLEAM survey. However, a noticeable (5–20%) false Stokes Q is still reported in the GLEAM calibrated data, which is attributed to the AEE beam model (see Figure 4 in Hurley-Walker et al. (2017) and Figure 11 in this paper). Furthermore, Offringa et al. (2016) attribute ~1% leakage into Stokes V to inaccuracies in the beam model, but the false Stokes V leakage can be higher on individual snapshot images (not averaged over long periods of time). Therefore, further improvements in the MWA beam model are required to improve the accuracy of the calibration and polarimetric measurements.

This paper is organised as follows. In Section 2, we present implementation of the new Full Embedded Element (FEE) primary beam model of MWA. We describe the physical representation of the tiles, the mathematical description in terms of spherical harmonics, and summarise differences with respect to the previous AEE model. In Section 3, we present a beam correction procedure that we implemented to test the new beam model. In Section 4, we present the results of this procedure applied to MWA data and compare performances of the new FEE model against the previous AEE model.

2 MODEL DESCRIPTION

Sutinjo et al. (2015a) proposed three tiers of beam modelling sophistication: (1) an analytic model using array theory and pattern multiplication, (2) using the AEE and incorporating mutual coupling, as detailed in that paper, and (3) using an FEE model. The improvements from the AEE model were incorporated into the mainstream MWA data processing pipelines. Here we improve the model again by using the FEE patterns for each of the bow-tie antennas to model the tile beam.

From the MWA user perspective, the approach is consistent with previous models where the beam pattern is generated on-the-fly for a given set of antenna delays ³ which define the pointing direction of the beam. Due to computational limitations, the embedded element pattern is calculated for the centre of every 1.28 MHz coarse channel. At arbitrary frequency channels, the beam model is calculated for the closest coarse channel which is within 0.64 MHz (i.e. no frequency interpolation was implemented at the current version of the beam model, but it can be added in future if even higher precision is required). The internals are different to previously implemented models in the following three key areas:

• • The physical model of the tile is improved (i.e. size of the ground screen now reflects the actual 5 × 5 m mesh size).

• • For the first time, the beam pattern is calculated using the FEE patterns of the 16 bow-tie dipoles in the tile dipole.

• • The beam pattern is calculated using spherical wave expansion (SWE), allowing accurate replication of full wave simulation of the beam pattern in any direction.

We describe these differences in more detail below.

2.1. Physical model of the tile

An MWA tile comprises of 16 bow-tie antennas aligned in the east–west (x) and north–south (y) directions, located on a regular 4 × 4 grid with 1.1-m spacing between centres, as shown in Figure 1. The normal operating frequency range is 75–240 MHz, but is capable of observations up to ~315 MHz. Tingay et al. (2013) and Neben et al. (2016) give a detailed description of the physical tile.

Figure 2 shows the bow-tie antenna modelled in the feko ⁴ electromagnetic simulation software. It has the same dimensions as the actual antenna, but uses a 15-mm diameter wire instead of the aluminium channels to reduce simulation complexity. This diameter was selected because it best matches the measured dimensions of the antenna elements. The ground mesh has spacing of 5 cm between adjacent wires, which at MWA frequencies we model as a perfect electrical conductor (PEC). Beyond the extent of the 5 × 5 m ground plane, we simulate the ground as soil from the Murchison Radio-astronomy Observatory (MRO) with 2% moisture, based on the permittivity and conductivity properties of soil from the MRO reported in Sutinjo et al. (2015b). Antennas are elevated by 10 cm from the mesh, but the plastic legs were not included in the simulation.

Figure 2. The bow-tie antenna modelled in feko. The marks along the arms indicate the simulation segments.

The model uses loaded ports as shown in Figure 3. The low noise amplifier (LNA) impedance is modelled using a lumped circuit, meaning the complex characteristics of the LNA impedance is predictable by measuring a simple shunt-series elements made up of resistors, inductors and capacitors (RLC). We represent the LNA input impedance with a 2 nH series element attached to a RLC (914 Ω, 450 nH, and 3.2 pF) shunt network. As per Figure 4, this model shows good agreement with the LNA impedance measured using a vector network analyser. The lumped circuit model results in a more compact, self-contained simulation file, and returns an impedance at arbitrary frequency resolution.

Figure 3. Bow-tie ports modelled in feko, where the ports are loaded with a voltage source.

Figure 4. Comparison between the measured LNA impedance and lumped circuit model (simulated in feko). We have verified that the small residual difference between the measured and model impedance has an insignificant (<1%) impact on the resulting beam model at 216 MHz.

2.2. Jones matrix beam model

The primary beam of aperture array telescopes strongly varies with pointing direction, and differs between the orthogonal antenna polarisations x and y. At a particular pointing direction θ (angle from zenith) and ϕ (azimuth angle, increasing clockwise from north through east), we can describe the instrumental effect of the MWA tile on the astronomical signal as a Jones matrix J, a 2 × 2 complex matrix. This Jones matrix maps the voltage at the beamformer output (v _x and v _y ) to the signal in orthogonal sky polarisations (e _θ and e _ϕ) as follows: v = Je, which in expanded form is

(1) $$\begin{equation} \left[\begin{array}{c}v_{x}\\ v_{y} \end{array}\right]=\left[\begin{array}{c@{\quad}c}J_{x\theta } & J_{x\phi }\\ J_{y\theta } & J_{y\phi } \end{array}\right]\left[\begin{array}{c}e_{\theta }\\ e_{\phi } \end{array}\right], \end{equation}$$

where v _x and v _y are the voltages from the x (E–W) and y (N–S) measurement bases (bow-tie antennas), and e _θ and e _ϕ are far-field unit vectors in spherical coordinate bases (Smirnov 2011a; Sutinjo et al. 2015a).

We can separate the MWA tile Jones matrix into two principal components:

(2) $$\begin{equation} \mathbf {J}=\mathbf {GE}, \end{equation}$$

where G is the direction-independent effect (DIE) due to complex electronic gain and E is the DDE due to the tile beam pattern. The latter varies as a function of the tile pointing (amplitude and phase weights of each antenna), and, for a given pointing, a Jones matrix applies for each (θ, ϕ) point in the hemisphere. In feko, E(θ, ϕ) can be calculated for a given pointing by applying an appropriate phase slope across the tile. However, such an approach would require a new simulation run for every desired pointing. Instead we model E(θ, ϕ) for each embedded element, being the in-situ radiation pattern of each bow-tie dipole in the tile (Kelley & Stutzman 1993). Each antenna is simulated in turn by setting its amplitude to a constant voltage, and the other 15 to zero. The resulting Jones matrix is normalised to the zenith (θ = 0) of a zenith-pointed beam (zero delay, thus no phase slope across the tile), therefore the absolute value of the excitation voltage is not important, as long as it is equal between antennas. The tile beam pattern for a given pointing can be determined post-simulation by weighting each element pattern with the appropriate phase and amplitude.

2.3. Computational representation via spherical wave expansion

Different methods are possible to represent the full wave simulation results from feko in the beam model. The method used previously for the AEE model was to output E for regular (θ, ϕ) intervals on the hemisphere. The disadvantage of this approach is that interpolation is required between adjacent (θ, ϕ) points. In this paper, we use SWE to calculate the tile beam pattern from the electric far-field according to the following formula (FEKO 2014, chapter ‘AS card’):

(3) $$\begin{eqnarray} \overrightarrow{E}^{\mathrm{ff}}(\theta ,\phi ) &=& \beta \sqrt{\frac{Z_{0}}{2\pi }}\frac{e^{-j\beta }}{\beta } \nonumber\\ &&\times\, \Bigg [\sum _{n=1}^{\infty }\sum _{m=-n}^{n}\frac{e^{jm\phi }C_{mn}}{\sqrt{n(n+1)}} \Big (\frac{-m}{\left|m\right|}\Big )^{m} \big ( e_{mn}^{\theta }\hat{\theta }+e_{mn}^{\phi }\hat{\phi }\big ) \Bigg ], \nonumber\\ \end{eqnarray}$$

where β is the wavenumber, Z ₀ is the intrinsic impedance of free space, C _mn = ((2n + 1)(n − |m|)!)/(2(n + |m|)!)^1/2 is the normalisation factor for the associated Legendre function, P ^|m| _n (cos θ), of order n and rank |m| (see Chapter 6 in Harrington 2001). The coefficients e ^θ _mn and e ^ϕ _mn can be calculated according to the following equations:

(4) $$\begin{eqnarray} e_{mn}^{\theta }& =& \Bigg [j^{n}\frac{P_{n}^{|m|}(\cos \theta )}{\sin \theta }\Big(|m|Q_{2mn}^{tile}\cos \theta -mQ_{1mn}^{tile}\Big ) \nonumber\\ &&+\, j^{n}Q_{2mn}^{tile}P_{n}^{|m|+1}(\cos \theta )\Bigg ], \end{eqnarray}$$

(5) $$\begin{eqnarray} e_{mn}^{\phi }& =& \Bigg [j^{n+1}\frac{P_{n}^{|m|}(\cos \theta )}{\sin \theta }\Big (mQ_{2mn}^{tile}-|m|Q_{1mn}^{tile}\cos \theta \Big ) \nonumber\\ &&-\,j^{n+1}Q_{1mn}^{tile}P_{n}^{|m|+1}(\cos \theta )\Bigg ], \end{eqnarray}$$

where s = 1 and s = 2 in $\mathbf {Q}_{\text{smn}}^{\text{tile}}$ refer to transverse electric (TE) and transverse magnetic (TM) modes, respectively. $\mathbf {Q}_{\text{smn}}^{\text{tile}}$ are vectors of coefficients formed as linear combination of 16 embedded elements beam patterns for specific pointing (θ, ϕ) according to beamformer time delays t _i and $\mathbf {Q}_{\text{smn}}^i$ are feko-generated coefficients for every ith antenna in the MWA tile. Hence, vectors $\mathbf {Q}_{\text{smn}}^{\text{tile}}$ can be calculated as

(6) $$\begin{equation} \mathbf {Q}_{\text{smn}}^{\text{tile}} = \sum _{i=1}^{16} e^{-2 \pi j\nu t_{i}} \mathbf {Q}_{\text{smn}}^{i}, \end{equation}$$

where ν is the observing frequency. The spherical harmonics approach allows rapid and accurate computation of the beam at one or more desired (θ, ϕ) coordinates.

We have implemented the new FEE beam model in python for MWAtools (an internal software package for MWA data processing) and in c/c++. The c/c++ implementation is included in the Real-Time System (rts) (Mitchell et al. 2008) and its off-line implementation by Offringa et al. (2016) known as calibrate. The c/c++ implementation can calculate the beam model for 1 million spatial points (1000 × 1000 pixels) in a few minutes. However, the python implementation requires spatial interpolation in the beam calculation to generate similarly sized beam models within a similar time frame.

3 FULL-STOKES BEAM CORRECTION PROCEDURE

The visibility matrix for the cross-correlation of tiles i and j is measured as

(7) $$\begin{equation} \mathbf{V}_{ij}=2\left[\begin{array}{c@{\quad}c}\langle v_{i,x}v_{j,x}^{*}\rangle & \langle v_{i,x}v_{j,y}^{*}\rangle \\ \langle v_{i,y}v_{j,x}^{*}\rangle & \langle v_{j,x}v_{j,x}^{*}\rangle \end{array} \right]=2 {\left[\begin{array}{c@{\quad}c}XX & {XY}\\ {YX} & {YY} \end{array}\right]}, \end{equation}$$

where v is the voltage from the beamformer of a given tile (i or j) and antenna polarisation (x or y). In the calibration observation scenario, a single bright source dominating the visibilities is observed and hence the brightness matrix B represents this single point source (delta function). Therefore, in the calibration scenario, for an ideal cross-correlator, V _ij is related to the brightness matrix B of the actual ‘single strong source sky’ by the response of each tile, represented as a Jones matrix J [equation (2)] ⁵ :

(8) $$\begin{equation} \mathbf {V}_{ij}=\mathbf {J}_{i}\mathbf {B}\mathbf {J}_{j}^{\mathbf {H}}, \end{equation}$$

where B =⟨ee ^H⟩ and the H superscript denotes the Hermitian transpose. Using equation (2), we can separate J as two Jones matrices:

(9) $$\begin{equation} \mathbf {V}_{ij}=\mathbf {G}_{i}\mathbf {E}_{i}\mathbf {B}\mathbf {E}_{j}^{\mathbf {H}}\mathbf {G}_{j}^{\mathbf {H}}, \end{equation}$$

where G describes the DIEs due to complex electronic gain and E the DDEs due to the tile beam pattern. This layered description of effects on the signal path is known as the ‘radio interferometer measurement equation’ (Hamaker, Bregman, & Sault 1996; Smirnov 2011a).

If the DIEs [G in equation (2)] are correctly calibrated for, we observe what Smirnov (2011b) calls the ‘apparent sky’, being the true sky attenuated by the tile beam E:

(10) $$\begin{equation} \mathbf {B_{{\rm app}}}=\mathbf {E}_{i}\mathbf {B}\mathbf {E}_{j}^{\mathbf {H}}. \end{equation}$$

If we also assume identical beam patterns for tiles i and j, the sky brightness matrix can be estimated using a model Jones matrix $\mathbf {\widetilde{E}}$ representing the MWA tile beam:

(11) $$\begin{equation} \widetilde{\mathbf {B}}=\mathbf {\widetilde{E}^{-1}B_{{\rm app}}}\left(\mathbf {\widetilde{E}}^{H}\right)^{-1}, \end{equation}$$

where the tilde designates a modelled matrix or a result estimated from models. We will use this convention through the remainder of the paper.

We note that the assumption of identical beam patterns for all tiles significantly simplifies the data processing, since corrections to images can all be applied in image space, as a linear combination of images made in instrumental polarisation coordinates.

From $\widetilde{\mathbf {B}}$ , we can calculate Stokes parameters (following the convention of Smi11-I):

(12) $$\begin{eqnarray} \widetilde{I} & =&(\widetilde{B}_{1,1}+\widetilde{B}_{2,2})/2,\nonumber\\ \widetilde{Q} & =&(\widetilde{B}_{1,1}-\widetilde{B}_{2,2})/2, \nonumber \\ \widetilde{U} & =&(\widetilde{B}_{1,2}+\widetilde{B}_{2,1})/2, \nonumber \\ \widetilde{V} & =&i(\widetilde{B}_{2,1}-\widetilde{B}_{1,2})/2. \end{eqnarray}$$

For the randomly polarised sky, we expect $\widetilde{Q}=\widetilde{U}=\break \widetilde{V}=0$ .

3.1. Calibration and beam correction

A standard calibration procedure is to observe a bright (dominating the visibilities), unresolved, and unpolarised source and solve for the complex gains of each tile via a least-square method. For an unpolarised calibrator source of intensity I, the sky brightness B is given by

(13) $$\begin{equation} \mathbf{B=} {\left[\begin{array}{c@{\quad}c}1 & 0\\ 0 & 1 \end{array} \right]}I, \end{equation}$$

where we follow the ‘Convention-1’ definition of Stokes I (Smirnov 2011a) which was implemented in the Common Astronomy Software Applications (CASA) (McMullin et al. 2007). Assuming identical tile beams, equation (8) becomes

(14) $$\begin{eqnarray} \hspace*{-6pt}&&\mathbf{V}_{ij} = I \nonumber\\ \hspace*{-6pt}&&\quad\times\,{\left[\begin{array}{cc}g_{i,x}(\vert E_{x\theta }\vert ^{2}+\vert E_{x\phi }\vert ^{2})g_{j,x}^{*} & g_{i,x}(E_{x\theta }E_{y\theta }^{*}+E_{x\phi }E_{y\phi }^{*})g_{j,y}^{*}\\[6pt] g_{i,y}(E_{y\theta }E_{x\theta }^{*}+E_{y\phi }E_{x\phi }^{*})g_{j,x}^{*} & g_{i,y}(\vert E_{y\theta }\vert ^{2}+\vert E_{y\phi }\vert ^{2})g_{j,y}^{*} \end{array}\right]}, \quad\qquad \end{eqnarray}$$

where g _x and g _y are the DIEs for the respective tile polarisations and E(θ, ϕ) is the beam pattern at the (θ, ϕ) direction of the calibrator source (we drop the θ, ϕ notation in the following equations). Calibration solves for the g _x values from the XX visibilities and likewise the g _y values from the YY visibilities. The DDEs are taken into account by either correcting the calibrator source model for the beam pattern E prior to solving for complex gains (as implemented in calibrate) or by dividing the resulting complex gains by amplitudes of the electric field of X and Y dipoles [ $(\vert \tilde{E}_{x\theta }\vert ^{2}+\vert \tilde{E}_{x\phi }\vert ^{2})^{1/2}$ and $(\vert \tilde{E}_{y\theta }\vert ^{2}+\vert \tilde{E}_{y\phi }\vert ^{2})^{1/2}$ , respectively]. Both approaches lead to DIE complex gains which can be represented as diagonal matrix:

(15) $$\begin{equation} \mathbf {\widetilde{G}} ={\left[\begin{array}{c@{\quad}c}\widetilde{g_{x}} & 0\\ 0 & \widetilde{g_{y}} \end{array}\right]}. \end{equation}$$

Note that we assume the complex gain matrix [equation (15)] to be diagonal, because off-diagonal terms are very small (negligible in comparison with mutual coupling of x and y dipoles) due to high isolation between the x and y analogue chains in MWA.

We introduce matrix D _ij as a specific instance of visibility data (dataset) obtained from some specific sky observation. Subsequently, calibration of the observed visibility data D _ij corrects for the DIEs on the XX, YY, XY, and YX polarisations:

(16) $$\begin{eqnarray} \mathbf {\widetilde{D}}_{ij} & =&\mathbf {\mathbf {\widetilde{G}}}_{i}^{-1}\mathbf {D}_{ij}(\mathbf {\widetilde{G}}_{j}^{H})^{-1} \nonumber\\ & =& {\left[\begin{array}{cc}D_{ij}(0,0)/(\widetilde{g}_{i,x}\widetilde{g}_{j,x}^{*}) & D_{ij}(0,1)/(\widetilde{g}_{i,x}\widetilde{g}_{j,y}^{*})\\ D_{ij}(1,0)/(\widetilde{g}_{i,y}\widetilde{g}_{j,x}^{*}) & D_{ij}(1,1)/(\widetilde{g}_{i,y}\widetilde{g}_{j,y}^{*}) \end{array}\right]}. \end{eqnarray}$$

In the image space, this is our measurement of the apparent sky [equation (10)], and from equation (11), we can calculate the sky brightness matrix as

(17) $$\begin{equation} \mathbf {\widetilde{B}}=\mathbf {\widetilde{E}}_{{\rm obs}}{}^{-1}\widetilde{\mathbf {D}}\left(\mathbf {\widetilde{E}}_{{\rm obs}}^{H}\right)^{-1}, \end{equation}$$

where $\tilde{\mathbf {E}}_{{\rm obs}}$ is our pointed beam model in the direction of the observed target source. It is then trivial to use equation (12) to calculate Stokes parameters.

3.2. Implementation of beam calibration pipeline

We have the following three beam models to test:

• • Analytic model of an array of Hertzian dipoles above a metallic ground plane (Sutinjo et al. 2015a).

• • AEE model reported in Sutinjo et al. (2015a).

• • FEE model described in this paper.

All three models have been implemented in MWA reduction software MWAtools, rts, and calibrate. We process the same observations independently for each beam model. The steps are as follows:

1. 1. Observe a calibrator source and use calibrate (Offringa et al. 2016) to solve for the tile-based direction-independent complex gains, a diagonal matrix $\mathbf {\widetilde{G}}$ [equation (15)]. The calibrate software incorporates the beam model under test into the calibration procedure.

2. 2. Apply direction-independent calibration solutions to visibilities from target field.

3. 3. Create sky images in all instrumental polarisations (XX, XY, YX, and YY) with the wsclean software (Offringa et al. 2014).

4. 4. Use the full-Stokes beam $\mathbf {\widetilde{E}}$ , modelled for all (θ, ϕ) directions of the observation tile beam pointing [equation (17)] to calculate the sky brightness matrix $\mathbf {\widetilde{B}}$ from the instrumental polarisations.

5. 5. Calculate images in Stokes polarisation according to equation (12).

6. 6. Use the Aegean source finder (Hancock et al. 2012) to identify sources in Stokes I images, measure their flux densities in Stokes I, Q, U, and V images and measure false Stokes Q, U, and V relative to Stokes I.

4 FULL-STOKES DEMONSTRATION ON MWA DATA

In the following sections, we will apply the above primary beam correction procedure to MWA data and compare performance of the three primary beam models. In the first section, we present comparison of the three models applied to the GLEAM data using our procedure as described in the previous section. In Section 4.2, we summarise the original GLEAM calibration procedure and how the beam model was applied, then in Section 4.3, we present a three-night sample of GLEAM data re-calibrated with the FEE model and compare the resulting false Stokes Q between the AEE and FEE model. Finally, in Section 4.4, we explain the false polarisation effect observed in the original GLEAM data calibrated with the AEE model.

4.1. Comparison of models performance on MWA data

To test the accuracy of the different beam models by measuring false leakages in all Stokes parameters (Q, U, and V), we developed a pipeline implementing our full calibration procedure according to steps described in Section 3.1 and applied it to typical GLEAM observations at high frequencies (200–230 MHz). These frequencies are most severely affected by inaccuracies of beam models, since some approximations of the physical model representation in feko become less accurate at shorter wavelengths (for example, arms of the dipoles are represented as 15-mm diameter wires instead of aluminium channels). At frequencies below 170 MHz, the false Stokes polarisation is below 5% for both AEE and FEE models. For calibration, we used a 116 s observation of Hydra A starting at 13:24:48 UTC on 2014-03–06 (LST ≈8.14 h) at (ϕ, θ) ≈ (52.1°, 22.0°). The following steps were performed in order to probe the differences between the models:

I. Solve for the tile-based direction-independent complex gains (Step 1 in Section 3.2)

As Hydra A is partially resolved at MWA frequencies (Sutinjo et al. 2015b; Lane et al. 2004), we only use visibilities with baselines of length 30–300~λ during calibration so that Hydra A appears as an unresolved source. The calibrate software uses the beam model to correct the calibrator model. Hence, the resulting complex gains are already direction-independent complex gains.

II. Measure false Stokes on the calibration observation (Steps 2–6 in Section 3.2)

As described earlier, the beam model is applied at two stages: to correct calibration solutions for the beam response in the direction of the calibrator source and later to beam-correct an image of another field observation. Therefore, we verify the accuracy of the calibration correction at both stages. In the first check, we apply the calibration solutions to correct the visibilities of the calibration observation, estimate the sky brightness matrix, and calculate the Stokes parameters. The Stokes Q, U, and V images of the calibrator fields are consistent with noise.

III. Transfer of calibration to another observation (Steps 2–6 in Section 3.2)

In order to test beam correction in a typical observation scenario, we applied the direction-independent complex gains to a set of drift-scan observations at (ϕ, θ) ≈ (0°, 28.3°) performed on 2014-03-06 between 11:44:47 and 13:16:40 UTC. Stokes I, Q, U, and V de-convolved images obtained with the wsclean and beam corrected with the FEE model are shown in Figure 5. The corresponding second-order polynomial surfaces fitted to the false Q, U, and V leakages of the brightest sources are shown in Figure 6. The false Stokes Q, U, and V averaged in 5° declination bins are shown in Figures 7–9, respectively.

Figure 5. Stokes I, Q, U, and V images (a,b,c,d, respectively) obtained from 2-min observation started at 13:14:48 UTC on 2014-03-06. The images were beam corrected in 1.28 MHz coarse channels and averaged in the 200–212 MHz band. Only part of the band was used to avoid radio-frequency interference that affected the upper part of the band (most likely due to digitial TV) which caused subtle artefacts in the Q, U, and V images. The images obtained with wsclean were beam corrected using the FEE model. The false Stokes leakages are within ±5% in image centres and get a bit higher closer to the edges. The second-order polynomial surfaces fitted to false Stokes Q, U, and V leakages are shown in Figure 6 and leakages averaged in declination bins are show in Figures 7–9. (a) Stokes I. (b) Stokes Q. (c) Stokes U. (d) Stokes V.

Figure 6. False Stokes Q, U, and V leakages (a,b,c, respectively) surfaces obtained from fit of second-order polynomial to leakages of the brightest sources from the Q, U, and V images in Figure 5 (the colour scale is the same for all three images). (a) Q leakage. (b) U leakage. (c) V leakage.

Figure 7. Comparison of Stokes Q leakage measured in images at 200–231 MHz calibrated with the three different models. The leakage data from individual sources were averaged in 5° bins (100–200 sources per bin). The analytic model performs the worst of all three models. The FEE models performs better at positive declinations and converges to AEE model at negative declinations (both are within measurements errors on this relatively small data sample).

The new FEE model has false Stokes Q below 5% falling down from ≈5 to 0% with increasing declination (Figures 7–9). The data calibrated with the AEE model has slightly higher false Stokes Q, but also within 5%. Both FEE and AEE models are better than the analytic model which has a false Stokes Q ≈ 30% (Figure 7). The false Stokes U is ≈5% for all three models and the false Stokes V in the data corrected with FEE model is consistent with 0%, whilst the V leakage in the data calibrated with the AEE model is ≈1–2% [of similar magnitude to that reported by Offringa et al. (2016) below 200 MHz].

However, the errors on this relatively small data sample are quite high and the results from all three models agree within the errors (except the false Stokes Q of the analytic model). These errors result from hour angle (HA) dependence of false Stokes leakage across the image (Figure 6), which is averaged in declination bins and the errors calculated as standard deviation correspond to variation in HA. Note that the errors of false Stokes U for the analytical model are significantly larger than for the FEE and AEE models (Figure 8) because the HA dependence of false Stokes U was reduced significantly for the non-analytic models. In the next section, we will show the effects of the new model on a larger (three nights) GLEAM data sample and how it improves the false Stokes Q originally observed in the calibrated GLEAM data.

Figure 8. Comparison of Stokes U leakage measured in images at 200–231 MHz calibrated with the three different models. The leakage data from individual sources were averaged in 5° bins (100–200 sources per bin). All three models have similar values of the U leakage ≈5%. Note that the errors of false Stokes U for the analytical model are significantly larger than for the FEE and AEE models (Figure 8) because the hour angle dependence of false Stokes U was reduced significantly for the non-analytic models.

Figure 9. Comparison of Stokes V leakage measured in images at 200–231 MHz calibrated with the three different models. The leakage data from individual sources were averaged in 5° bins (100–200 sources per bin). All the models have values within the errors bars, but AEE model has V leakage ≈1–2% consistent across the band, which is of similar magnitude to that reported by Offringa et al. (2016) below 200 MHz.

4.2. Original GLEAM calibration procedure

One of the goals of the GLEAM survey was to catalogue the flux density of all radio sources below +30° declination in the 72–231-MHz frequency band, but polarisation measurements were not initially a priority. However, because beam-corrected instrumental XX and YY images have been calibrated independently to the Molonglo Reference Catalogue (MRC) catalogue (Large et al. 1981), it was identified that the ratio of YY and XX fluxes deviate from unity (equivalent to non-zero false Stokes Q) away from the image centres. The ratio has noticeable structure as a function of declination [sources were grouped in bins in declination and frequency as shown in Figure 4 in Hurley-Walker et al. (2017)]. The observed false Stokes Q were attributed by the authors to deficiencies in the AEE primary beam model, as the structure of YY/XX ratio versus declination remains the same between different nights (sometimes separated by three months), which excludes the possibility of the effect being due to variations in the ionosphere.

In the next section, we will show how application of the new FEE beam model reduces the false Stokes Q observed in a sample of GLEAM data taken over three nights. Inaccuracies in the AEE model (which was used for GLEAM processing) resulted in flux-scale variations as a function of beam pointing and declination. To overcome this, a robust flux calibration procedure was developed for the GLEAM pipeline, as outlined in Hurley-Walker et al. (2017).

The GLEAM pipeline only provided calibrated XX and YY images and so beam model effectiveness can only be tested against Stokes Q. We used the new FEE beam model in the original GLEAM calibration pipeline in order to verify if it corrected the originally reported false Stokes Q.

4.3. Application to the GLEAM data

We applied the new FEE beam model to the GLEAM data from three pointings at the local meridian and declinations δ = −13°, +1.6°, and +18.6° collected during nights starting on the 5th, 7th, and 11th of November 2013, respectively. The GLEAM data were processed as described in Hurley-Walker et al. (2017), but instead of using AEE model, the new FEE beam model was applied. We analysed images in four 7.68 MHz bands (200–208, 208–216, 216–223, and 223–231 MHz). The false Stokes Q leakages were measured for many sources (nearly 85 000 in total, with ≈40 100 at δ = −13°, ≈29 000 at δ = 1.6° and ≈16 000 at δ = 18.6°) using Aegean (step 6 in Section 3.2).

The comparison of false Stokes Q leakages by applying the AEE (as originally used in the GLEAM survey) and FEE models is shown in Figure 10. The false Stokes Q leakages resulting from calibration with the AEE model (upper plot in Figure 10) have significant declination-dependent structure (reaching values around and above 10% at the edges of the images) as a function of declination for each of the analysed fields. In contrast, sources in the images obtained from calibration with the new FEE model have smaller false Stokes Q, which is consistent with zero for the two pointings (δ = −13° and δ = +1.6°) and non-zero (within ±10%) for the δ = 18.6°, which is at the lowest elevation (≈45°). Moreover, only the lowest elevation has a noticeable structure as a function of declination.

Figure 10. False Stokes Q leakages calculated from GLEAM images corrected with the originally used AEE model (upper image) and the new FEE model (lower image). The small scattered data points represent false Stokes Q calculated for all the individual sources (around 85 000 in total in both cases) identified in the images of the three fields at declinations δ = −13°, 1.6°, and +18.6° (with approximately 40 100, 29 000, and 16 000 sources, respectively). The large data points with error bars were calculated as the mean and standard deviation in 3° bins.

4.4. False Stokes Q expected in GLEAM due to inaccuracy of the beam model

The GLEAM survey identified that the ratio of YY/XX flux densities deviates from unity away from the pointing centres and has a declination-dependent structure as shown in Figure 4 in Hurley-Walker et al. (2017) and Figure 11 in the current paper (only the data >200 MHz). The effect was reduced, but still noticeable, by the mosaicking procedure which up-weighted (by square of the beam pattern E ²) the measurements taken closer to the beam centre (local meridian). Using the YY/XX ratios observed in GLEAM, we calculated the equivalent in terms of false Stokes Q, which is in a range of approximately 5–20%.

Figure 11. The data points are the ratio $\widetilde{B}_{2,2} / \widetilde{B}_{1,1}$ based on the GLEAM data from bottom row (frequency ranges 200–208, 208–216, 216–223, and 223–231 MHz) of Figure 4 in Hurley-Walker et al. (2017). The data were collected with the beam pointing at local meridian at declination ≈2°. The GLEAM data were binned in 3° bins in declination (with about 1 000–1 200 sources averaged in the central bins down to 100–200 in the bins near the image edges). The solid lines represent the same ratio predicted when the AEE beam model is used to correct unpolarised sky brightness propagated through the FEE beam model (assumed to represent the ‘true’ MWA tile beam). The simulated curves were normalised by values in the image centre in order to replicate the normalisation of flux to The Molonglo Reference Catalogue (MRC) catalogue (Large et al. 1981) performed on the GLEAM data. (a) 200–208 MHz. (b) 208–216 MHz. (c) 216–223 MHz. (d) 223–231 MHz.

Because of the way the GLEAM data were processed, i.e. images in XX and YY instrumental polarisations were beam corrected and their flux densities normalised independently to the MRC catalogue, the observed structure was attributed to inaccuracies of the AEE beam model (ionosphere was excluded by intra-night stability of the declination structure). Therefore, with the new FEE model available, we were able to verify if the observed declination structure could be reproduced by assuming a priori the FEE beam model to be a ‘true’ (or a better representation) of tile beam and applying the AEE model to calibrate the data. In this test, we only tested beam-related effects and did not take electronic DIE gain into account. We assumed that the unpolarised sky brightness B_unpol is given by equation (13) and hence we calculated the ‘calibrated’ sky brightness according to the following equation:

(18) $$\begin{equation} \widetilde{\mathbf {B}} = \mathbf {E_{AEE}}^{-1} \Big [ \mathbf {E_{FEE} B_{unpol} E_{FEE}^H} \Big ] \left( \mathbf {E_{AEE}}^H \right)^{-1}. \end{equation}$$

The comparison between the ratio YY/XX observed in the original GLEAM data [frequency bands 200–231 MHz in Figure 4 in Hurley-Walker et al. (2017)] and our prediction for the same quantity (calculated as $\widetilde{B}_{2,2} / \widetilde{B}_{1,1}$ ) is shown in Figure 11. The GLEAM authors attributed the declination-dependent structure in the ratio YY/XX as a result of an inaccurate beam model. Our predictions, shown in Figure 11, support this statement.

5 CONCLUSIONS AND DISCUSSION

An accurate primary beam model of the telescope is required to pursue many science goals of the MWA and the future SKA-low telescope. We have presented a new FEE, which is the most rigorous realisation of the beam model for the MWA, superseding the previous AEE and analytical beam models (Sutinjo et al. 2015a). The new model was generated in the feko software using an improved physical representation of MWA tile. In the simulation, every dipole in an MWA tile was simulated in transmit mode with other dipoles not transmitting (voltages set to zero). The feko simulation results were exported using spherical harmonic representation, which enables calculation of the beam pattern and Jones matrices in arbitrary pointing directions without re-running the simulations (instead of discrete directions as in the previous AEE model).

We used polarisation measurements to compare the three beam models. We have developed a beam calibration pipeline to calibrate the direction-dependent and independent effects in MWA observations and measure false Stokes Q, U, and V as a metric for accuracy of the models. We applied this procedure to a set of 12 × 2-min GLEAM observations. The analytical model shows very high (~30%) Q leakage, considerably higher than the measured leakage in both the AEE and FEE models. The new FEE model has false Stokes Q leakage 0–5%. The AEE model gives higher absolute value of leakage at higher declinations and similar at lower declinations. The false Stokes U is similar for all three models (within ±5%). The V leakage resulting from the FEE model is consistent with zero, whilst the V leakage resulting from AEE models is ≈1–2% [of similar magnitude to that reported by Offringa et al. (2016) below 200 MHz]. However, they are both within measurement errors on the relatively small data sample we used for this test. The pipeline enables further tests of false Stokes leakages in all pointing directions in order to identify further avenues to improve the model.

We have also applied the FEE model to the original GLEAM calibration pipeline (Stokes Q only) to correct three nights of GLEAM data and calculated false Stokes Q for nearly 85 000 sources in these images. The FEE model reduces the declination-dependent structure of false Stokes Q in comparison to the AEE model and the false Stokes Q is consistent with zero in two out of three pointing directions.

We used the new FEE model to understand the Q leakage (and its structure as a function of declination) observed in the original GLEAM data. Using the new model as a ‘hypothetical true’ MWA beam model and calibrating it with the previous (AEE) beam model, we were able to reproduce the declination-dependent structure in false Stokes Q observed in GLEAM survey (Wayth et al. 2015; Hurley-Walker et al. 2017).

Although the current model is the best representation of the tile beams yet, there are still possibilities for further improvements in the physical model of the MWA tile by incorporating finer details into the simulation. However, the better physical representation of the physical tile the longer the simulation takes. Therefore, the accuracy of the physical model is a trade-off between the required precision and practical constraints (for instance, simulation time).

Another possibility is that the observed false Stokes leakages result from deviation of the actual MWA telescope from the ‘ideal’ instrument represented by the model. Although most of the component variations should average to zero and not contribute significantly to differences between XX and YY polarisations causing false Stokes, some effects may not cancel out. We performed simulations and generated 128 beams taking into account information about faulty dipoles (which were disabled at beamformers), calculated the mean beam and performed a test similar to the one described in Section 4.4 (treating the mean beam as the ‘true’ MWA beam and beam-corrected with the ideal model beam). Results of these tests showed that faulty dipoles (usually in one polarisation), especially when fault distribution is non-uniform within tiles, can lead to false Stokes leakage of the order similar to that observed [as just one faulty or disabled dipole represents 1/16 (≈6%) of the MWA tile]. In future, this effect can be tested with observed data. Any residual leakage as a result of inaccuracies in the physical model can be reduced by measuring, mapping, and subtracting the residual leakage on a per snapshot basis as described by Lenc et al. (2017).

Finally, it should be noted that electro-magnetic simulations have certain limitations and the resulting models are loaded with error (due to imperfection of the physical representation, variations of the components, numerical limitations, etc.). The acceptable error of the simulation is considered to be of the order of a few percent which also corresponds to typical (~0.2–0.4 dB) spread in RF component parameters. The currently used simulation resolution is already optimal as we have verified that increasing the mesh resolution by factor of two results in ≲1% difference in the beam response, but at a significantly higher computing cost (factor of 12). The false Stokes leakages at the level of a few percent are close to current limitations of the electro-magnetic simulations, highlighting the fact that astronomical polarisation measurements offer a very sensitive and rigorous ways of testing and validating electro-magnetic simulations.