2.1. The coherency

We describe the proprieties of the electric field signal from the sky with a coherency vector, given by

(1) \begin{equation} \boldsymbol{{S}}(\boldsymbol{\unicode{x03B8}}) = \left[\begin{array}{c} \left\langle |E_{\text{R}}(\boldsymbol{\unicode{x03B8}})|^2 \right\rangle \\[2pt] \left\langle |E_{\text{D}}(\boldsymbol{\unicode{x03B8}})|^2 \right\rangle \\[2pt] \left\langle E_{\text{R}}(\boldsymbol{\unicode{x03B8}}) E_{\text{D}}^*(\boldsymbol{\unicode{x03B8}}) \right\rangle \\[2pt] \left\langle E_{\text{R}}^*(\boldsymbol{\unicode{x03B8}}) E_{\text{D}}(\boldsymbol{\unicode{x03B8}}) \right\rangle \end{array}\right]\end{equation}

(Hamaker, Bregman, & Sault Reference Hamaker, Bregman and Sault1996). Here $E_{\text{R}}(\boldsymbol{\unicode{x03B8}})$ and $E_{\text{D}}(\boldsymbol{\unicode{x03B8}})$ are the components of the electric field in two orthogonal directions. In our convention, we define these directions to align with the RA and Dec. directions on the sky, respectively. $\boldsymbol{\unicode{x03B8}}$ is a two-element vector defining the position on the sky. The angled brackets $\langle \rangle$ denote the time average and the asterisk $^*$ represents the complex conjugate.

The coherency can equivalently be described as a $2\times2$ matrix (Hamaker et al. Reference Hamaker, Bregman and Sault1996; Smirnov Reference Smirnov2011). Likewise, the vector ordering and orthogonal basis are arbitrary. Equation (1) presents the convention used in the fhd analysis.

2.2. The Stokes parameters

Polarised emission is often described with respect to the Stokes parameters I, Q, U, and V (Stokes Reference Stokes1851). Stokes I corresponds to the total intensity, Stokes Q and U correspond to linear polarisation, and Stokes V corresponds to circular polarisation. The Stokes parameters are related to the coherency vector $\boldsymbol{{S}}(\boldsymbol{\unicode{x03B8}})$ via the relationship

(2) \begin{equation} \left[\begin{array}{c} I(\boldsymbol{\unicode{x03B8}}) \\ Q(\boldsymbol{\unicode{x03B8}}) \\ U(\boldsymbol{\unicode{x03B8}}) \\ V(\boldsymbol{\unicode{x03B8}}) \end{array}\right] = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & i & -i \\ \end{array}\right] \boldsymbol{{S}}(\boldsymbol{\unicode{x03B8}}).\end{equation}

Note that, because Equation (1) defines the coherency vector with respect to the RA/Dec. coordinate system, the Stokes parameters are also referenced to that coordinate system. We emphasise that the Stokes parameters are basis-dependent, and it is critical that any polarised measurements specify the basis used (Hamaker & Bregman Reference Hamaker and Bregman1996). Under the RA/Dec. coordinate system, the Stokes parameters are undefined at the North and South Celestial Poles. Consequently, analyses of fields at or near the poles may benefit from choosing a different orthogonal basis. See Ludwig (Reference Ludwig1973) for examples of alternative polarisation bases that may be used to define the Stokes parameters.

2.3. The Jones matrix

The Jones matrix is a $2\times2$ complex matrix that represents the polarised antenna response to the sky (Jones Reference Jones1941). The matrix transforms between the true electric field on the sky and the electric field measured by the instrument (Hamaker et al. Reference Hamaker, Bregman and Sault1996; Sault, Hamaker, & Bregman Reference Sault, Hamaker and Bregman1996; Hamaker & Bregman Reference Hamaker and Bregman1996; Hamaker Reference Hamaker2000, Reference Hamaker2006; Ord et al. Reference Ord2010; Smirnov Reference Smirnov2011):

(3) \begin{equation} \left[\begin{array}{c} \epsilon_{j \text{p}}(\boldsymbol{\unicode{x03B8}}) \\ \epsilon_{j \text{q}}(\boldsymbol{\unicode{x03B8}}) \end{array}\right] = \boldsymbol{\mathsf{J}}^{\text{ZA}}_j(\boldsymbol{\unicode{x03B8}}) \left[\begin{array}{c} E_{\text{Z}}(\boldsymbol{\unicode{x03B8}}) \\ E_{\text{A}}(\boldsymbol{\unicode{x03B8}}) \end{array}\right].\end{equation}

Here $E_{\text{Z}}(\boldsymbol{\unicode{x03B8}})$ and $E_{\text{A}}(\boldsymbol{\unicode{x03B8}})$ are components of the electric field aligned with the zenith angle and azimuth directions, respectively. $\epsilon_{j \text{p}}(\boldsymbol{\unicode{x03B8}})$ and $\epsilon_{j \text{q}}(\boldsymbol{\unicode{x03B8}})$ represent the contribution of emission from sky direction $\boldsymbol{\unicode{x03B8}}$ to the measurements made by antenna j, and the subscripts p and q correspond to the two instrumental polarisation modes. For example, for an antenna with dipole elements aligned with the cardinal directions, p could refer to the east-west aligned dipole and q could refer to the north-south aligned dipole. $\boldsymbol{\mathsf{J}}^{\text{ZA}}_j(\boldsymbol{\unicode{x03B8}})$ is the Jones matrix for antenna j. The superscript ZA indicates that it corresponds to the zenith angle/azimuth basis on the sky. This is the usual basis convention for reporting the Jones matrix.

This formalism is fully general in that it assumes dual-polarisation antenna but makes no further assumptions about the antenna type. Each antenna could consist of dish with a dual-polarisation feed (e.g., GMRT or HERA), a simple crossed-dipole element (as with PAPER or the OVRO-LWA), or a beamformed array of dipole elements (such as LOFAR’s stations or the MWA’s tiles). The two antenna polarisations could be orthogonal but need not be. They could measure linearly or circularly polarised modes. Note that here we have defined a per-antenna Jones matrix that depends on the antenna index j: antennas need not have homogeneous polarised responses. The Jones matrix is direction-dependent, and we require a $2\times2$ Jones matrix for each location on the sky $\boldsymbol{\unicode{x03B8}}$ . It is also complex-valued: the direction-dependent complex phase of the Jones matrix elements captures timing delays in the instrumental response to an incoming wavefront.

In Figure 2 we plot the four elements of the Jones matrix of an MWA tile, as modelled by Sutinjo et al. (Reference Sutinjo2015), as a function of position on the sky. Each MWA tile, pictured in Figure 1, consists of 16 dual-polarisation beamformed elements. The Jones matrix model in Figure 2 corresponds to the full tile response in its zenith-pointed mode, averaged across a frequency range of 167–198 MHz. This Jones matrix defines the instrumental response for fhd analyses including the polarised mapping in Byrne et al. (Reference Byrne2022) and the EoR analyses in Barry et al. (Reference Barry2019b) and Li et al. (Reference Li2019). While the Jones matrix is complex-valued, the complex phase is near-zero, and for simplicity we plot the real part only. Note that elements of the Jones matrix exhibit a discontinuity at zenith. This is a feature of the basis in which it is defined and not a characteristic of the physical antenna response.

Figure 1. A photo of an MWA tile. Each tile consists of 16 dual-polarisation beamformed elements, and the full array consists of 256 such tiles (Wayth et al. Reference Wayth2018). The tile’s response to incident radiation is estimated by a beam model. Three equivalent representations of the beam model are shown in each Figures 2, 3, and 4 and 5. Photo credit: Natasha Hurley-Walker, the MWA Collaboration, and Curtin University.

Figure 2. The elements of the Jones matrix for a zenith-pointed MWA tile (pictured in Figure 1) at 167–198 MHz, as modelled by Sutinjo et al. (Reference Sutinjo2015). The Jones matrix defines and instrument’s polarised response. It is a $2\times2$ complex matrix defined at each point on the sky; here we plot the real part only. The top row depicts the response of the east-west aligned, or p, tile polarisation while the bottom row depicts the response of the north-south aligned, or q, tile polarisation. Here the Jones matrix is normalised such that the peak response amplitude of each tile polarisation is one. This Jones matrix is defined with respect to the zenith angle/azimuth basis (see Equation (3)). The left and right columns corresponds to the tile’s response to electric field emission polarised in the zenith angle and azimuth directions, respectively. The Jones matrix elements exhibit a discontinuity at zenith as a result of the pole of the coordinate system.

Equation (3) and Figure 2 present the Jones matrix in the zenith angle/azimuth coordinate system, but in Equation (1) we define the coherency with respect to the RA/Dec. basis. Our analysis is therefore simplified if we define a Jones matrix $\boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_j(\boldsymbol{\unicode{x03B8}})$ with respect to RA/Dec. The transformation between bases involves rotating by the direction-dependent parallactic angle $\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})$ :

(4) \begin{equation} \left[\begin{array}{c} \boldsymbol{{e}}_Z(\boldsymbol{\unicode{x03B8}}) \\ \boldsymbol{{e}}_A(\boldsymbol{\unicode{x03B8}}) \end{array}\right] = \left[\begin{array}{c@{\quad}c} \sin [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] & -\cos [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] \\ -\cos [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] & -\sin [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] \end{array}\right] \left[\begin{array}{c} \boldsymbol{{e}}_{\text{R}}(\boldsymbol{\unicode{x03B8}}) \\ \boldsymbol{{e}}_{\text{D}}(\boldsymbol{\unicode{x03B8}}) \end{array}\right].\end{equation}

Here $\boldsymbol{{e}}_{\text{R}}(\boldsymbol{\unicode{x03B8}})$ , $\boldsymbol{{e}}_{\text{D}}(\boldsymbol{\unicode{x03B8}})$ , $\boldsymbol{{e}}_Z(\boldsymbol{\unicode{x03B8}})$ , and $\boldsymbol{{e}}_A(\boldsymbol{\unicode{x03B8}})$ are unit vectors in the RA, Dec., zenith angle, and azimuth directions, respectively. The parallactic angle is given by

(5) \begin{equation} \unicode{x03D5}(\boldsymbol{\unicode{x03B8}}) = \tan^{-1} \left( \frac{-\sin \alpha}{\cos \delta \tan \delta_{\text{zen}} - \sin \delta \cos \alpha} \right),\end{equation}

where $\alpha$ and $\delta$ are the RA and Dec. coordinates of $\boldsymbol{\unicode{x03B8}}$ , respectively, and $\delta_{\text{zen}}$ is the declination at zenith. For a full derivation of this expression see Byrne (Reference Byrne2021), Appendix D.

Here Equation (5) assumes that the Earth’s axis of rotation aligns with the poles of the RA/Dec. coordinate system. Under this assumption the parallactic angle is time-independent: for an array at a given declination, the transformation between the zenith angle/azimuth and RA/Dec. polarisation bases does not depend on time. In practice, the Earth’s precession and nutation introduce deviations in the alignment of the rotational axis. While fhd fully accounts for precession and nutation when calculating source positions on the sky, where small positional errors can substantially degrade analysis accuracy, it does not account for them when transforming between polarisation bases, instead using the parallactic angle calculation presented in Equation (5). This is a good approximation as the errors from the neglected precession and nutation effects produce only small errors in the reconstructed polarisation angle. Future extensions to fhd’s polarised imaging pipeline could add these higher-order effects to the calculation of the polarisation basis transformation. This would introduce time dependence to the parallactic angle calculation in Equation (5) and by extension the basis transformation in Equation (4).

Applying the transformation in Equation (4) to the Jones matrix, we get that

(6) \begin{equation} \boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_j(\boldsymbol{\unicode{x03B8}}) = \boldsymbol{\mathsf{J}}^{\text{ZA}}_j(\boldsymbol{\unicode{x03B8}}) \left[\begin{array}{c@{\quad}c} \sin [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] & -\cos [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] \\ -\cos [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] & -\sin [\unicode{x03D5}(\boldsymbol{\unicode{x03B8}})] \end{array}\right],\end{equation}

where

(7) \begin{equation} \left[\begin{array}{c} \epsilon_{j \text{p}}(\boldsymbol{\unicode{x03B8}}) \\ \epsilon_{j \text{q}}(\boldsymbol{\unicode{x03B8}}) \end{array}\right] = \boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_j(\boldsymbol{\unicode{x03B8}}) \left[\begin{array}{c} E_{\text{R}}(\boldsymbol{\unicode{x03B8}}) \\ E_{\text{D}}(\boldsymbol{\unicode{x03B8}}) \end{array}\right].\end{equation}

$\boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_j(\boldsymbol{\unicode{x03B8}})$ , the Jones matrix in the RA/Dec. basis, is used throughout the fhd analysis because it corresponds with the coherency vector defined in Equation (1). However, it is more typical to report models of the polarised antenna response in terms of $\boldsymbol{\mathsf{J}}^{\text{ZA}}_j(\boldsymbol{\unicode{x03B8}})$ because the RA/Dec. basis is dependent on the antenna’s location on the Earth and varies with latitude.

In Figure 3 we plot the elements of the MWA Jones matrix in RA/Dec. coordinates by applying the transformation in Equation (6) to the Jones matrix presented in Figure 2. Once again, we plot the real part of the Jones matrix elements only. Note that we no longer see a discontinuity at zenith. Instead, the Jones matrix elements exhibit a discontinuity near the lower edge of the plots, corresponding to the position of the South Celestial Pole.

Figure 3. The Jones matrix of an MWA tile, plotted in Figure 2, now recast in the RA/Dec. coordinate system (see Equation (6)). Once again, the top and bottom rows correspond to the polarised responses of the east-west and north-south aligned antenna polarisation, respectively. However, the left column now corresponds to the tile’s response to emission polarised in the RA direction while the right column corresponds to the response to emission polarised in the Dec. direction. As in Figure 2, we plot the real part of the Jones matrix elements only. Since the RA/Dec. coordinate system has poles at the North and South Poles, we see a discontinuity at the bottom edge of the plots, corresponding to the position of the South Pole relative to the MWA’s $-27^\circ$ latitude.

2.4. The Mueller matrix

The Mueller matrix is a $4\times4$ complex matrix formed by the Kronecker product of two Jones matrices:

(8) \begin{equation} \boldsymbol{\mathsf{M}}_{jk}(\boldsymbol{\unicode{x03B8}}) = \boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_j(\boldsymbol{\unicode{x03B8}}) \otimes {\boldsymbol{\mathsf{J}}^{\text{R} \text{D}}_k}^*(\boldsymbol{\unicode{x03B8}})\end{equation}

(Hamaker et al. Reference Hamaker, Bregman and Sault1996). Here we define the Mueller matrix with respect to the RA/Dec. basis to once again align with the convention in Equation (1).

The Mueller matrix defines the mapping between the sky signal and visibilities. We can write the visibilities formed from correlating antennas j and k in terms of the Mueller matrix:

(9) \begin{equation} \boldsymbol{{v}}_{jk} = \int d^2\boldsymbol{\unicode{x03B8}} \, \boldsymbol{\mathsf{M}}_{jk}(\boldsymbol{\unicode{x03B8}}) \boldsymbol{{S}}(\boldsymbol{\unicode{x03B8}}) \, e^{2 \pi i \boldsymbol{\unicode{x03B8}} \cdot \boldsymbol{{u}}_{jk}}.\end{equation}

Here $\boldsymbol{{u}}_{jk}$ is the baseline vector for antennas j and k in units of wavelengths and $\boldsymbol{{v}}_{jk}$ is a 4-element vector with components

(10) \begin{equation} \boldsymbol{{v}}_{jk} = \left[\begin{array}{c} v_{jk\text{p} \text{p}} \\ v_{jk\text{q} \text{q}} \\ v_{jk\text{p} \text{q}} \\ v_{jk\text{q} \text{p}} \end{array}\right]\end{equation}

where, for example, $v_{jk\text{p} \text{q}}$ represents the visibility formed by correlating the p-polarised signal from antenna j with the q-polarised signal from antenna k, and p and q once again refer to the two instrumental polarisations of a dual-polarised antenna. The integral is taken across the visible sky, and the Mueller matrix $\boldsymbol{\mathsf{M}}_{jk}(\boldsymbol{\unicode{x03B8}})$ contains the appropriate visibility normalisation. Polarised imaging, discussed in Section 4, amounts to inverting Equation (9) to estimate $\boldsymbol{{S}}(\boldsymbol{\unicode{x03B8}})$ from the measured visiblities.

Equation (9) assumes a coplanar array, where each baseline can be described by the two-element vector $\boldsymbol{{u}}_{jk}$ . Non-coplanar baselines introduce an additional phase term in the integrand of Equation (9), but fhd neglects this term and assumes a reasonably coplanar array. Moderate non-coplanarity can be corrected with w-projection (Cornwell & Perley Reference Cornwell and Perley1992).

4.1. Gridding

Gridding reconstructs the uv plane from the visibilities. Here we describe gridding the $\text{p} \text{q}$ -polarised visibilities; for the $\text{p} \text{p}$ , $\text{q} \text{q}$ , and $\text{q} \text{p}$ polarisations, simply substitute for the subscripts $\text{p} \text{q}$ below.

For a given baseline formed by correlating signals from antennas j and k, we denote the $\text{p} \text{q}$ -polarised beam model $B_{jk\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}})$ . Note that this represents one element of the $4\times4$ diagonal matrix $\boldsymbol{\mathsf{B}}_{jk}(\boldsymbol{\unicode{x03B8}})$ presented in Equation (17). Under optimal mapmaking, the gridding kernel is the Fourier transform of this beam model, which we will denote $\widetilde{B}_{jk\, \text{p} \text{q}}(\boldsymbol{{u}})$ . Here $\boldsymbol{{u}}$ represents the uv plane coordinate, which is Fourier dual to the sky coordinate $\boldsymbol{\unicode{x03B8}}$ and has units of wavelengths. The tilde indicates the Fourier transformed quantity: $\widetilde{B}_{jk\, \text{p} \text{q}}(\boldsymbol{{u}}) = \mathcal{FT} \left[ B_{jk\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}}) \right]$ .

Gridding with this kernel produces the following reconstructed uv plane:

(19) \begin{equation} \hat{\widetilde{S}}_{\text{app}\, \text{p} \text{q}} (\boldsymbol{{u}}) = \frac{\sum_{jk} \widetilde{B}_{jk\, \text{p} \text{q}}(\boldsymbol{{u}}_{jk}-\boldsymbol{{u}}) \, v_{jk\, \text{p} \text{q}}}{\sum_{jk} \widetilde{B}_{jk\, \text{p} \text{q} }(\boldsymbol{{u}}_{jk}-\boldsymbol{{u}})}.\end{equation}

Here $\hat{S}_{\text{app}\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}})$ is the apparent sky in the $\text{p} \text{q}$ polarisation, defined as the estimate of the sky multiplied by the beam amplitude; $\hat{\widetilde{S}}_{\text{app}\, \text{p} \text{q}} (\boldsymbol{{u}})$ is its Fourier transform. The hat symbol $ \, \hat{\,} \,$ indicates that this is the reconstructed estimate. j and k index antennas, and $\sum_{jk}$ denotes the sum over all baselines.

The numerator in Equation (19) describes the gridding operation. If the model accurately captures the true instrumental beam, this method produces an optimal and lossless sky estimate (Tegmark Reference Tegmark1997; Bhatnagar et al. Reference Bhatnagar, Cornwell, Golap and Uson2008; Morales & Matejek Reference Morales and Matejek2009). Errors in the beam model produce errors in the reconstructed intensities, polarisation, and, in rare cases, the positions of sources. Precision beam modelling is therefore an active area of research (Newburgh et al. Reference Newburgh2014; Sutinjo et al. Reference Sutinjo2015; Berger et al. Reference Berger2016; Wayth et al. Reference Wayth, Colegate, Sokolowski, Sutinjo and Ung2016; Sokolowski et al. 2017; Line et al. Reference Line2018; Fagnoni et al. Reference Fagnoni2021; Chokshi et al. Reference Chokshi2021; etc.).

The denominator in Equation (19) is known as the uv weights and is equivalent to gridding each visibility with a value of unity. If the instrument does not have complete uv plane measurement coverage, the weights will be equal to zero in certain uv plane locations. We set the expression in Equation (19) equal to zero where the weights are zero. Other imaging approaches introduce regularisation methods to, in effect, interpolate over incomplete uv measurements (see Johnson et al. Reference Johnson2017, Eastwood et al. Reference Eastwood2018, and Eastwood et al. Reference Eastwood2019).

Equation (19) represents just one uv weighting scheme. fhd supports alternative weightings as well, corresponding to different denominators in Equation (19). For example, the denominator can be replaced with unity (often called ‘natural weighting’) or with the number of visiblities that contribute to each uv pixel (often called ‘uniform weighting’). However, Equation (19) presents a sensible choice in which the reconstructed uv pixels are scaled by the measurement sensitivity. This weighting scheme, which we call ‘optimal weighting’, was used to produce the diffuse maps in Byrne et al. (Reference Byrne2022).

Because $\boldsymbol{\mathsf{B}}_{jk}(\boldsymbol{\unicode{x03B8}})$ is defined to be diagonal, $\text{p} \text{q}$ -polarised visiblities contribute to the $\text{p} \text{q}$ uv plane only, and likewise the other visibility polarisations contribute only to their respective uv planes. This means that each visibility is only gridded once. This is a unique feature of the instrumental basis, as any other basis choice introduces off-diagonal elements of $\boldsymbol{\mathsf{B}}_{jk}(\boldsymbol{\unicode{x03B8}})$ and requires that each visibility is gridded four times. This holds true even when the Mueller matrix has significant off-diagonal components: the instrumental basis need not be orthogonal. Because gridding is the most computationally intensive step in the imaging pipeline, the instrumental basis allows for significant computational savings.

4.2. Image reconstruction

Transforming the gridded uv plane into sky coordinates gives

(20) \begin{equation} \hat{S}_{\text{app}\, \text{p} \text{q}} (\boldsymbol{\unicode{x03B8}}) = \mathcal{FT}^{-1} \left[ \hat{\widetilde{S}}_{\text{app}\, \text{p} \text{q}} (\boldsymbol{{u}}) \right],\end{equation}

where $\mathcal{FT}^{-1}$ is the inverse Fourier transform operator. This apparent sky map can then be converted into an estimate of the true sky, in the instrumental $\text{p} \text{q}$ polarisation, by undoing the beam weighting:

(21) \begin{equation} \hat{S}_{\text{inst}\, \text{p} \text{q}} (\boldsymbol{\unicode{x03B8}}) = \frac{1}{B_{\text{avg}\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}})} \mathcal{FT}^{-1} \left[ \hat{\widetilde{S}}_{\text{app}\, \text{p} \text{q}} (\boldsymbol{{u}}) \right].\end{equation}

Here $\hat{S}_{\text{inst}\, \text{p} \text{q}} (\boldsymbol{\unicode{x03B8}})$ is an element of the instrumental coherency estimate, as defined in Equation (18). $B_{\text{avg}\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}}) = \langle B_{jk\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}}) \rangle$ is the average $\text{p} \text{q}$ beam amplitude, averaged across baselines. For a homogeneous array, $B_{jk\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}}) = B_{\text{avg}\, \text{p} \text{q}}(\boldsymbol{\unicode{x03B8}})$ for all baselines $\{j, k\}$ .

Gridding all four visibility polarisations, $\text{p} \text{p}$ , $\text{q} \text{q}$ , $\text{p} \text{q}$ and $\text{q} \text{p}$ and Fourier transforming their respective uv planes produces an estimate of the instrumental coherency $\hat{\boldsymbol{{S}}}_{\text{inst}} (\boldsymbol{\unicode{x03B8}})$ . This is a reconstruction of the polarised sky signal defined with respect to the non-orthogonal instrumental polarisation basis, plotted for the MWA in Figure 5. For science applications, we convert this instrumental coherency into a fixed, orthogonal basis on the sky or Stokes parameters.

From Equation (18), we transform the estimated instrumental coherency into the RA/Dec. basis coherency:

(22) \begin{equation} \hat{\boldsymbol{{S}}}(\boldsymbol{\unicode{x03B8}}) = \boldsymbol{\mathsf{L}}^{-1}(\boldsymbol{\unicode{x03B8}}) \hat{\boldsymbol{{S}}}_{\text{inst}} (\boldsymbol{\unicode{x03B8}}).\end{equation}

In general $\boldsymbol{\mathsf{L}}^{-1}(\boldsymbol{\unicode{x03B8}})$ is well-defined and can be calculated by numerically inverting $\boldsymbol{\mathsf{L}}(\boldsymbol{\unicode{x03B8}})$ . However, $\boldsymbol{\mathsf{L}}(\boldsymbol{\unicode{x03B8}})$ is singular when the instrumental basis vectors are parallel, as is the case for the MWA at the horizon. Above the horizon $\boldsymbol{\mathsf{L}}(\boldsymbol{\unicode{x03B8}})$ is invertible.

Note that the coherency estimate $\hat{\boldsymbol{{S}}}(\boldsymbol{\unicode{x03B8}})$ must be calculated from the instrumental coherency $\hat{\boldsymbol{{S}}}_{\text{inst}} (\boldsymbol{\unicode{x03B8}})$ , not from the apparent sky estimate $\hat{\boldsymbol{{S}}}_{\text{app}} (\boldsymbol{\unicode{x03B8}})$ . Because the instrumental beam amplitude is polarisation-dependent (see Figure 4), the apparent sky estimate is intrinsically defined in the instrumental polarisation basis.

Finally, we calculate the Stokes parameters from the coherency estimate $\hat{\boldsymbol{{S}}}(\boldsymbol{\unicode{x03B8}})$ via Equation (2):

(23) \begin{equation} \left[\begin{array}{c} \hat{I}(\boldsymbol{\unicode{x03B8}}) \\ \hat{Q}(\boldsymbol{\unicode{x03B8}}) \\ \hat{U}(\boldsymbol{\unicode{x03B8}}) \\ \hat{V}(\boldsymbol{\unicode{x03B8}}) \end{array}\right] = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & i & -i \\ \end{array}\right] \hat{\boldsymbol{{S}}}(\boldsymbol{\unicode{x03B8}}).\end{equation}

This produces an optimal, lossless estimate of the fully polarised true sky signal.

A deconvolution step can be added to the image reconstruction pipeline to reduce imaging artefacts from the array’s Point Spread Function (PSF). Deconvolution assumes that sources are compact on the sky and thereby differentiates between true source emission and the PSF. Sullivan et al. (Reference Sullivan2012) describes fhd’s deconvolution algorithm. Deconvolution may not be necessary or desirable when imaging large-scale, diffuse structure, or when processing data that nearly completely samples the uv plane, as in Byrne et al. (Reference Byrne2022). It should be noted that, without mitigation through deconvolution, the imaging approach described in this section produces images that are convolved with the array’s PSF, as in Figure 6.

Figure 6. Example of polarised image reconstruction with fhd, based on a simulation of a single Stokes Q polarised point source with a polarisation fraction of 50%. The source is located in the upper right quadrant of each image at a zenith angle of 10 $^\circ$ (RA 23h30m04s, Dec. -19.4 $^\circ$ ). Zenith is marked with a plus symbol. Visibilities were simulated with the pyuvsim simulation package (Lanman et al. Reference Lanman2019) and based on the MWA Phase I at 167–198 MHz. The visibilities were then gridded and imaged with fhd to produce instrumental polarisation images (top row; see Equation (20)) and Stokes images (bottom row; see Equation (23)). The instrumental polarised images $\text{p} \text{p}$ and $\text{q} \text{q}$ are real-valued, but the $\text{p} \text{q}$ and $\text{q} \text{p}$ images are complex-valued and complex conjugates of one another. We therefore plot the real and imaginary components of the $\text{p} \text{q}$ -polarised image and do not plot the $\text{q} \text{p}$ -polarised image. While the simulated source appears predominantly in the $\text{p} \text{p}$ and $\text{q} \text{q}$ images, a small amount of power couples into the $\text{p} \text{q}$ image as a result of non-orthogonality of the instrumental basis at the source location. The reconstructed Stokes images have a 49.52% total polarisation fraction and 49.47% Stokes Q polarisation fraction at the location of the simulated source.

Figure 6 presents an example of fhd’s image reconstruction with a simulated polarised point source. Visibilities were simulated with the pyuvsim simulation package^{Footnote b} (Lanman et al. Reference Lanman2019). Although fhd itself performs polarised source simulation, pyuvsim enables very accurate, semi-analytical interferometric simulations and serves as an independent verification of the fhd’s polarised imaging pipeline. The simulation is based on a zenith-pointed observation made with the MWA Phase I at 167–198 MHz. We simulated a 10 Jy point source located north-west of zenith at a zenith angle of $10^\circ$ and polarised in Stokes Q with a polarisation fraction of 50%. The simulated visibilities were gridded and imaged with fhd, and Figure 6 depicts the resulting images in instrumental polarisation (top row) and Stokes (bottom row). The images are not deconvolved, so the simulated source appears convolved with the array’s PSF.

The top row of Figure 6 presents images of the simulated source in instrumental polarisation, corresponding to the quantity $\hat{\boldsymbol{{S}}}_{\text{app}} (\boldsymbol{\unicode{x03B8}})$ from Equation (20). Although the simulated source appears primarily in the $\text{p} \text{p}$ and $\text{q} \text{q}$ polarisations, we see a small amount of power in the $\text{p} \text{q}$ -polarised image as well. The source’s polarisation preferentially couples power into the $\text{q} \text{q}$ image, and the $\text{q} \text{q}$ image therefore has a greater amplitude than the $\text{p} \text{p}$ image. The power in the $\text{p} \text{q}$ image is due to non-orthogonality of the instrumental polarisation basis.

The bottom row of Figure 6 depicts the reconstructed Stokes images ( $\hat{I}(\boldsymbol{\unicode{x03B8}})$ , $\hat{Q}(\boldsymbol{\unicode{x03B8}})$ , $\hat{U}(\boldsymbol{\unicode{x03B8}})$ , and $\hat{V}(\boldsymbol{\unicode{x03B8}})$ from Equation (23)). At the source location, the reconstructed images have a polarisation fraction of 49.52%, of which 49.47% appears in Stokes Q. This aligns well with the simulation input of 50% fractional polarisation in Stokes Q.

4.3. Comparison with other widefield polarised imagers

fhd’s polarised imaging pipeline joins the field of other polarised widefield imagers including wsclean (Offringa et al. Reference Offringa2014) and the rts (Mitchell et al. Reference Mitchell2008). It is perhaps most similar to the imaging approach discussed in Tasse et al. (Reference Tasse, van der Tol, van Zwieten, van Diepen and Bhatnagar2013), which, like fhd, is based on optimal mapmaking (A-projection) and grids visibilities with a kernel derived from the beam model.

Tasse et al. (Reference Tasse, van der Tol, van Zwieten, van Diepen and Bhatnagar2013) describes an analysis pipeline built for widefield imaging with LOFAR. Although LOFAR’s station configurations differ across the array, each station has identical polarisation alignment. This enables the analysis to define an implicit instrumental basis for visibility gridding, such that each visibility is gridded to just one uv plane. Tasse et al. (Reference Tasse, van der Tol, van Zwieten, van Diepen and Bhatnagar2013) transforms between the instrumental basis and an orthogonal polarisation basis in the uv plane, applying a correction factor to the gridded visibilities before deconvolving. This is in contrast to fhd, which performs this transformation in the image plane after dividing by the instrumental beam amplitude.

wsclean (Offringa et al. Reference Offringa2014) is a fully polarised widefield imager that was originally built for the MWA but has been widely utilised for analysis of data from instruments including LOFAR, GMRT, the Very Large Array (VLA), the OVRO-LWA, and the Australian Square Kilometre Array Pathfinder. It was not developed as an A-projection algorithm, although it supports a variety of gridding kernels and can grid with the instrumental beam model (Lynch et al. Reference Lynch2021). In its polarised imaging mode, wsclean produces instrumentally polarised $\text{p} \text{p}$ , $\text{q} \text{q}$ , $\text{p} \text{q}$ and $\text{q} \text{p}$ images. Much like in fhd, it then follows with a beam correction step that transforms the images into Stokes parameters. This step corrects for direction-dependent image attenuation from the beam amplitude and gridding kernel and performs the polarisation basis transformation.

Mitchell et al. (Reference Mitchell2008) describes the rts, a GPU-accelerated calibration and imaging pipeline for the MWA, also discussed in Tingay et al. (Reference Tingay2013). This pipeline also performs widefield imaging and produces fully polarised images. Much like wsclean and fhd, it produces images in the instrumental polarisation basis and then uses the Mueller matrix to transform to Stokes parameters. Like wsclean, it was not initially developed to perform A-projection, and it typically does not use the instrumental beam for gridding.

5.1. Gain parameterisation

In its most general form, $\boldsymbol{\mathsf{G}}_{jk}$ is given by

(25) \begin{equation} \boldsymbol{\mathsf{G}}_{jk} = \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} g_{j\text{p}\text{p}} \, g_{k\text{p}\text{p}}^* & g_{j\text{p}\text{q}} \, g^*_{k\text{p}\text{q}} & g_{j\text{p}\text{p}} \, g^*_{k\text{p}\text{q}} & g_{j\text{p}\text{q}} \, g^*_{k\text{p}\text{p}} \\ g_{j\text{q}\text{p}} \, g^*_{k\text{q}\text{p}} & g_{j\text{q}\text{q}} \, g^*_{k\text{q}\text{q}} & g_{j\text{q}\text{p}} \, g^*_{k\text{q}\text{q}} & g_{j\text{q}\text{q}} \, g^*_{k\text{q}\text{p}} \\ g_{j\text{p}\text{p}} \, g^*_{k\text{q}\text{p}} & g_{j\text{p}\text{q}} \, g^*_{k\text{q}\text{q}} & g_{j\text{p}\text{p}} \, g^*_{k\text{q}\text{q}} & g_{j\text{p}\text{q}} \, g^*_{k\text{q}\text{p}} \\ g_{j\text{q}\text{p}} \, g^*_{k\text{p}\text{p}} & g_{j\text{q}\text{q}} \, g^*_{k\text{p}\text{q}} & g_{j\text{q}\text{p}} \, g^*_{k\text{p}\text{q}} & g_{j\text{q}\text{q}} \, g^*_{k\text{p}\text{p}} \end{array}\right].\end{equation}

Here $g_{j\text{p}\text{p}}$ and $g_{j\text{q}\text{q}}$ are the complex gains of the p and q polarisations of antenna j, respectively. $g_{j\text{p}\text{q}}$ and $g_{j\text{q}\text{p}}$ are the cross gains. $g_{j\text{p}\text{q}}$ , for example, denotes the degree to which signal we expect to appear only in the q polarisation of antenna j also appears in the p polarisation (Sault et al. Reference Sault, Hamaker and Bregman1996).

The calibration solutions are then calculated by minimising a cost function given by

(26) \begin{equation} \chi^2(g) = \sum_{jk} \left( \boldsymbol{{v}}_{jk} - \boldsymbol{\mathsf{G}}_{jk} \boldsymbol{{m}}_{jk} \right)^{\dagger} \left( \boldsymbol{{v}}_{jk} - \boldsymbol{\mathsf{G}}_{jk} \boldsymbol{{m}}_{jk} \right),\end{equation}

where $\sum_{jk}$ denotes the sum over all baselines. Here we omit all frequency dependence and assume that the cost function is minimised independently for each frequency channel. For discussion of fhd-based precision bandpass calibration techniques, see Barry et al. (Reference Barry2019a,Reference Barryb), Li et al. (Reference Li2019).

If an instrument experiences minimal cross-polarisation signal coupling, as is the case for the MWA, we can safely set all the cross gains to zero. $\boldsymbol{\mathsf{G}}_{jk}$ is then diagonal, and we can denote each antenna gain with a single polarisation index: $g_{j\text{p}}$ and $g_{j\text{q}}$ . The calibration cost function then becomes

(27) \begin{equation} \chi^2(g) = \sum_{jk} \sum_{ab} \left| v_{jk\,ab} - g_{ja} g^*_{kb} m_{jk\,ab} \right|^2,\end{equation}

where a and b each index the two instrumental polarisation modes $\{\text{p}, \text{q}\}$ . fhd’s polarised calibration currently does not support non-zero cross gains.

5.2. Constraining the cross-polarisation phase

fhd initially did not use the cross-polarisation visibilities $\boldsymbol{{v}}_{\text{p} \text{q}}$ and $\boldsymbol{{v}}_{\text{q} \text{p}}$ in calibration. Excluding the cross-polarisation visibilities makes calibration separable in polarisation. Two independent cost functions are minimised, corresponding to the two instrumental polarisations:

(28) \begin{equation} \chi_{\text{p}}^2(g_{\text{p}}) = \sum_{jk} \left| v_{jk\,\text{p} \text{p}} - g_{j\text{p}} g^*_{k\text{p}} m_{jk\,\text{p} \text{p}} \right|^2\end{equation}

and

(29) \begin{equation} \chi_{\text{q}}^2(g_{\text{q}}) = \sum_{jk} \left| v_{jk\,\text{q} \text{q}} - g_{j\text{q}} g^*_{k\text{q}} m_{jk\,\text{q} \text{q}} \right|^2.\end{equation}

However, calibrating in this way introduces a new calibration degeneracy, corresponding to the overall phase difference between the p and q gains across all antennas. We can identify this degeneracy by noting that the transformation $g_{\text{p}} \rightarrow g_{\text{p}} \, e^{-i\Delta/2}$ and $g_{\text{q}} \rightarrow g_{\text{q}} \, e^{i\Delta/2}$ does not affect the calibration solutions. We call the parameter $\Delta$ the ‘cross-polarisation phase’.

The cross-polarisation phase is not degenerate when we calibrate with the cross-polarisation visibilities, as in Equation (26) or (27)—provided these cross visibilities are non-zero. Although it is sometimes asserted that this phase can only be constrained by calibrating to a polarised source (Sault et al. Reference Sault, Hamaker and Bregman1996; Bernardi et al. Reference Bernardi2013; Lenc et al. Reference Lenc2017), in the widefield limit the non-orthogonality of the instrumental polarisation basis couples appreciable unpolarised source power into the cross-polarisation visibilities. This allows for constraint of the cross-polarisation phase even while using a fully unpolarised sky model. (Note that this technique was independently developed in Anderson Reference Anderson2019.)

When adapting fhd to perform fully polarised calibration, we chose to largely retain the original calibration pipeline and to supplement it with an additional step to constrain the cross-polarisation phase. As a result, all calibration parameters other than the cross-polarisation phase are fit from the single polarisation visibilities $\boldsymbol{{v}}_{\text{p} \text{p}}$ and $\boldsymbol{{v}}_{\text{q} \text{q}}$ . Since the original calibration pipeline constrains the relative phase of the gains across frequency (Beardsley et al. Reference Beardsley2016; Barry et al. Reference Barry2019a,Reference Barryb; Li et al. Reference Li2019), the degenerate cross-polarisation phase amounts to a single parameter across all antennas and frequencies. If each frequency channel were calibrated independently, we would need to calculate a cross-polarisation phase at each frequency. We fit the cross-polarisation phase by plugging the solutions into the fully polarised cost function given by Equation (27). We let $g_{\text{p}} = \hat{g}_{\text{p}} \, e^{-i\Delta/2}$ and $g_{\text{q}} = \hat{g}_{\text{q}} \, e^{i\Delta/2}$ , where $\hat{g}$ are the gains calibrated up the the cross-polarisation phase. This gives

(30) \begin{align} \chi^2(\Delta) = & \sum_{jk} \left[ \left| v_{jk \, \text{p} \text{q}} - e^{-i\Delta} \hat{g}_{j \text{p}} \, \hat{g}^*_{k \text{q}} \, m_{jk \, \text{p} \text{q}} \right|^2 \right. \nonumber\\ & \left. + \left| v_{jk \, \text{q} \text{p}} - e^{i\Delta} \hat{g}_{j \text{q}} \, \hat{g}^*_{k \text{p}} \, m_{jk \, \text{q} \text{p}} \right|^2 \right].\end{align}

We can calculate the cross-polarisation phase analytically by finding the value $\hat{\Delta}$ that minimises $\chi^2(\Delta)$ . We find that

(31) \begin{equation} \hat{\Delta} = \operatorname{Arg} \left[ \sum_{jk} \left( v^*_{jk \, \text{p} \text{q}} \, \hat{g}_{j \text{p}} \, \hat{g}^*_{k \text{q}} \, m_{jk \, \text{p} \text{q}} + v_{jk \, \text{q} \text{p}} \, \hat{g}^*_{j \text{q}} \, \hat{g}_{k \text{p}} \, m^*_{jk \, \text{q} \text{p}} \right) \right],\end{equation}

where $\operatorname{Arg}$ denotes the complex phase.