2. Pipeline overview

The AusEoRPipe is designed to take raw data from the MWA correlator, calibrate for instrumental and ionospheric effects, subtract foregrounds, and then estimate a PS. Here, we list the different software packages that make up the AusEoRPipe and suggest simulations needed to test them. This section is intended to give the reader a broad overview of the AusEoRPipe and the overall goals of this series of papers, to give this paper context. At the end of this Section, we detail which parts of the work are covered in this paper.

Birli Footnote a is designed to take raw MWA data and preprocess it, including geometric correction, averaging, and radio frequency interference (RFI) flagging (via AOFlagger Offringa Reference Offringa2010 or SINSS Wilensky et al. Reference Wilensky, Morales, Hazelton, Barry, Byrne and Roy2019). Birli has been written to replace Cotter (Offringa et al. Reference Offringa2015), to have one unifying preprocessing package that works with both legacy and new MWA data (see Morrison et al. Reference Morrison2023, for details of the old and new MWA correlators). Given Birli has been tested to reproduce results out of Cotter, which has been tested and used within a number of MWA pipelines, we choose not to test Birli via simulation. To do so however, one would need to simulate raw MWA data, which for the original correlator would mean producing visibilities with no phase tracking, adding a frequency-dependent bandpass from the MWA polyphase band filter, and add in realistic RFI.

hyperdrive Footnote b (Jordan et al. submitted), takes either raw or preprocessed data, calibrates it, and then subtracts foregrounds. Built to replace the RTS (Mitchell et al. Reference Mitchell, Greenhill, Wayth, Sault, Lonsdale, Cappallo, Morales and Ord2008), hyperdrive also creates model visibilities to calibrate against via an image-based sky model. It initially performs a direction-independent (DI) calibration step, that derives a gain for each MWA receiving element (tile). All frequencies are calibrated independently. After DI calibration, foregrounds can be subtracted. The Fourier transform along frequency to create the PS is expected to separate the 21-cm signal from the foregrounds, as foreground emission is supposed to be spectrally smooth, whereas the 21-cm signal has spectral structure. Incomplete sampling of the u,v plane however causes mode mixing between the two, as well as other frequency-related instrumental effects. A detection therefore relies on the foregrounds being subtracted without injecting false spectral structure, as to preserve the underlying 21-cm signal. To test this functionality, a simulation needs to contain both a 21-cm signal and foregrounds, as well as frequency-dependent instrumental gain errors. One can also add frequency-dependent effects like cable reflections (e.g. Ewall-Wice et al. Reference Ewall-Wice2016) and the bandpass to investigate how calibrating each frequency channel independently performs. These effects are direction-independent and so can be added to visibilities post simulation (functionality which is in development in WODEN). Further simulations including RFI and ionospheric refraction would test how robust the calibration is to environmental influences. These effects are direction-dependent and so must be adding during the calculations of the visibilities, which WODEN is currently incapable of.

CHIPS (Trott et al. Reference Trott2016) takes the calibrated and subtracted visibilities and grids them into a spectral cube using an optimised Blackman-harris kernel. It then Fourier-transforms along frequency, and then either cylindrically averages to create a two-dimensional PS (2D PS) or spherically averages to create a one-dimensional PS (1D PS). Integrations of tens to hundreds of hours of data into a 1D PS have set the limits that have been released by the AusEoRPipe to date. The minimal test for CHIPS is to simulate a set of visibilities with a known 21-cm signal and see if CHIPS can recover it. Beyond that, taking in calibrated visibilities out of hyperdrive, derived from simulations with realistic instrumental effects and foregrounds, allows us to probe the effects of calibration on the resultant PS.

WODEN,Footnote c as previously mentioned, simulates MWA visibilities from a sky based model. It can use the same primary beam model as hyperdrive, the FEE beam (Sokolowski et al. Reference Sokolowski2017) (via the mwa_hyperbeam Footnote d package) and creates visibilities that can be fed directly into either hyperdrive or CHIPS. WODEN calculates the interferometric measurement equation for every component in the sky model. Each component can either be a point source (a dirac-delta function upon the sky), an elliptical Gaussian, or a Shapelet model (see Line et al. Reference Line2020, for more detail on Shapelets). The measurement equation encodes baseline and sky projection effects inherently, making the simulator accurate across the large field of view of the MWA. The caveat being that the sky model must be broken into the small components listed here, meaning representations of large-scale structure such as the galactic diffuse emission need millions of components. Aside from the primary beam model, the WODEN code base is entirely independent to hyperdrive. Although they use the same methodology to generate visibilities, this redundancy in code is by design; any internal bug that may cancel out when generating visibilities and calibrating them all with the same code base is avoided.

Along with the aforementioned packages, extensive quality metrics are generated at various stages (Nunhokee et al. submitted). These metrics are used to cull any observations that are deemed unusable due to insurmountable instrumental effects, ionospheric conditions, RFI events, contamination by bright sources in the primary beam sidelobes, and other effects. The simulations detailed above could be used to check whether these metrics are able to catch poor data and check whether the limits set for each metric still allow for a detection of the 21-cm signal.

In this paper, we ignore calibration entirely, and focus on generating a realistic 21-cm sky model to simulate through WODEN and directly input into CHIPS. As visibilities are additive, once we are able to generate accurate 21-cm simulations, we can add a variety of foregrounds and instrumental effects to test the calibration and subtraction step, without the need to rerun the more expensive 21-cm simulations. We leave testing hyperdrive to the second paper in this series.

2.1 MWA EoR observing

The MWA is a low-frequency radio interferometer, with receiving elements consisting of 4 $\times$ 4 grids of dual linear-polarisation bow-tie antennas. These ‘tiles’ are electronically steered through beamforming. To estimate the 21-cm PS, the MWA EoR collaboration have identified a number of fields with lower sky temperatures, which have now been observed for hundreds of hours (e.g. Trott et al. Reference Trott2020), in an attempt to average over thermal noise. Given the electronic beamforming, this has resulted in a drift-and-shift observational campaign, where the target field drifts through a number of pointings, with scheduling keeping the field centre as close to the primary beam centre as possible (see Jacobs et al. Reference Jacobs2016, for more details).

In this paper we focus on the EoR0 field, centred at RA, DEC $=0^h, -27^\circ$ , observed between frequencies of 167–198 MHz (known as the high-band). These parameters have consistently yielded the best MWA limits (e.g. Barry et al. Reference Barry2019b; Trott et al. Reference Trott2020; Chege et al. Reference Chege, Jordan, Lynch, Trott, Line, Pindor and Yoshiura2022). We also focus on five pointings, labelled -2 through +2 (c.f. Beardsley et al. Reference Beardsley2016), where -2 means two pointings before the zenith 0 pointing, and +2 meaning two pointings after 0 pointing. There is an approximate 15 $^\circ$ shift between pointings, all along the meridian. Examples are shown in Fig. 1. MWA EoR data are typically taken in two minute chunks known as an observation.

Figure 1.

Demonstration of the MWA primary beam and its interaction with the 21-cm sky model, via the: (a) -2 pointing; (b) -1 pointing; (c) Zenith pointing. Each plot is locked to the observer in HA/Dec. The solid black lines represent the instrumental Stokes I primary beam pattern contoured at 1%, 10%, 50%, and 90% power levels. The coloured images show the apparent 21-cm sky model after attenuation by the primary beam at three different LSTs, with the corresponding colour lines marking the edges of the full 21-cm sky model. Both the beam and sky model are shown at 167 MHz, where the primary beam is largest for a high-band observation. Note that the +1, +2, pointings as described in Section 5 are simply westward reflections of the -1, -2 pointings, so aren’t reproduced here for brevity.

3. 21-cm sky-model

We use the 21-cm lightcone detailed in Greig et al. (Reference Greig, Wyithe, Murray, Mutch and Trott2022), which is derived from a simulated EoR cube with sides of length $7.5\,$ Gpc. This cube was generated using a simplified version of 21CMFACST (Mesinger&Furlanetto Reference Mesinger and Furlanetto2007; Mesinger, Furlanetto, & Cen Reference Mesinger, Furlanetto and Cen2011) and was spectrally sampled at 80 kHz to match typical MWA EoR analysis parameters. This methodology generates a volume that projects to a sky coverage of $\sim 50 \times 50$ square degrees, which is essential given the footprint of the MWA primary beam. A number of steps are necessary to translate the lightcone box into a WODEN sky model. The box is a collection of 2D projected 21 cm intensity distributions, each at a unique redshift. These cartesian projections are taken by slicing a full 3D simulation volume at regular redshift intervals, as the simulation is evolved with time. This leaves us with a number of transformations needed to allow WODEN to ingest the model:

• • Each 2D 21-cm slice has a pixel resolution $\Delta x$ in cMpc; each slice therefore has a different angular resolution, as this quantity is redshift-dependent. To run efficiently, WODEN needs a grid of pixels constantly sampled in RA/DEC, necessitating a redshift-dependent angular interpolation.

• • The box is regularly sampled in redshift ( $\Delta z$ ), which results in irregular sampling in frequency ( $\Delta \nu$ ). Interferometric data is sampled regularly in frequency, necessitating a second interpolation over frequency.

• • The box is in units of mK. WODEN ingests units of Jansky (Jy), which is an integrated flux density. The conversion from mK to Jy/sterrad is straightforward, however, the physical volume of each pixel changes as a function of frequency. Along with the angular extent, the pixel volume is determined by $\Delta z$ . Given the interpolation over angle and frequency, this volume change must be taken into account, as it effects the resultant variance of the 21-cm map. This effect is noted in the script make_flat_spectrum_eor.py Footnote e in the Python package pyradiosky.

In the following, we use the cosmology assumed by Greig et al. (Reference Greig, Wyithe, Murray, Mutch and Trott2022), a $\Lambda$ CDM cosmology with: $H_0=68;\; \Omega_{M}=0.31;\; \Omega_{\Lambda}=0.69;\; \Omega_{b}=0.048$ . We use astropy (Astropy Collaboration et al. 2013, 2018, 2022) for all cosmological calculations.

We choose to interpolate all slices to the finest angular resolution in the box to retain as much angular structure as possible. We call this the reference redshift, $z_{\mathrm{ref}} = 7.5693$ . Given $\Delta \theta = \Delta x / \mathcal{D}(z)$ , where $\mathcal{D}$ is the comoving distance, we interpolate to an angular resolution of $\sim$ 27 arcsec. We found interpolating the Cartesian slice directly into a TAN FITS projection (a gnomic projection; see Calabretta & Greisen Reference Calabretta and Greisen2002, for details) centred at RA, DEC $=0^h, -27^\circ$ returned the expected PS. We experimented with bilinear and bicubic interpolation, but found significant signal loss in the resultant PS. It’s possible some variant of Gaussian Process Regression may be more effective, however for the purposes of this work we found a simple nearest-neighbour approximation was sufficient. Similarly, we applied a nearest-neighbour interpolation along the frequency axis, as again we saw signal loss in the final PS when applying linear or cubic interpolation. Along lines of sight, the 21-cm signal rapidly fluctuates between positive and zero, and so these kinds of interpolation tend to create false signal. Once these interpolations have been applied, we then scale for the change in volume of the pixel ( $\mathrm{V}_{\mathrm{pix}}(z)$ ) and its effect on the variance. We do this for each redshift by calculating via the differential comoving volume and multiplying by $\Delta z$ . We then scale each redshift slice by a factor $C_{\mathrm{pix}}$ given by

(1) \begin{equation} C_{\mathrm{pix}} = \sqrt{\frac{\mathrm{V}_{\mathrm{pix}}(z)}{\mathrm{V}_{\mathrm{pix}}(z_{\mathrm{ref}})}},\end{equation}

with the square root ensuring the variance is scaled by the pixel volume.

Once we have interpolated, scaled, and transformed each slice in Jy/sterrad, we convert each pixel into a point source with units of Jy by multiplying by the pixel solid angle. Ideally, we would the find some way to tile these maps to give an all-sky 21-cm sky model. However, given we cannot interpolate the box without signal loss, it is computationally infeasible. At the lowest frequencies, the number of point sources necessary would approach 200 million. However, given the extreme volume of the simulated 21-cm box, without tiling, the sky model already covers the main lobe of the MWA primary beam down to the $~1$ % power level for pointings of interest (see Fig. 1).

To investigate whether the projection steps detailed in this Section induce any form of signal loss or bias, we apply them to a purely Gaussian noise simulation in Section 4. We report of the results of the 21-cm sky model in Section 5.

4. White Gaussian noise sky simulation

The PS P(k) of a White Gaussian noise distribution $\mathcal{N}$ of mean $\mu=0$ and variance $\sigma^2$ is proportionate to the variance, that is, $P(k) \propto \sigma^2$ . By creating a sky model of purely White Gaussian noise, we can therefore predict an output P(k) and check the normalisation of the AusEoRPipe. Given we have outlined a method to project a 21-cm box from units of cMpc in Section 3, we can simply repeat the entire process, starting from a box of purely White Gaussian noise (from hereon in referred to as the Gaussian noise simulation). We choose to set $\sigma^2=16\,\mathrm{mK}^2$ as this yields a power comparable to those expected from the 21-cm signal. We typically report the PS in units of mK $^2$ Mpc $^3$ $h^{-3}$ , where h is the reduced Hubble constant; we use $h=0.68$ in this work. We can therefore predict our expected noise PS $P_N$ via

(2) \begin{equation} P_N(k) = \sigma^2 \mathrm{V}_{\mathrm{pix}}(z_\mathrm{ref}) h^3 = 8.923 \, \mathrm{mK}^2 \mathrm{Mpc}^3 h^{-3}\end{equation}

We use WODEN Footnote f to simulate a single zenith pointing observation with an 8 s time cadence and 80 kHz frequency resolution, using a frequency-interpolated version of the FEE beam. We use the convention that Stokes I = (XX + YY)/2, where XX and YY are the two linear polarisations. Given we are using simulated data with isotropic and unpolarised sky models, there is little difference between XX and YY, and therefore the power spectra shown are effectively Stokes I. All PS shown in this paper are from the north-south aligned polarisation. We run CHIPS on the simulated visibilities, to produce the 2D PS as shown in Fig. 2. All CHIPS plots were generated using the Python package CHIPS_wrappers.Footnote g

Figure 2.

Data from a single zenith observation, without any correction for gridding density, where: (a) shows the 2D PS; (b) the ratio to the expected value; c) the CHIPS gridding weights. $k_\perp$ modes are derived from the instrument sampling in the visibility u,v plane, averaged down to one dimension, and $k_\parallel$ from the fourier transform of the visibilities along frequency. The CHIPS gridding weights therefore show the u,v gridding density averaged into one dimension.

Figure 3.

Observations simulated in this paper, based on two real nights of MWA EoR0 observing. The blue diamonds show the simulated LST1 subset, green hexagons the LST2 subset, and orange circles the LST3 subset. Each hollow square shows a different two minute snapshot which was not simulated but exists in the real data set. Dividing lines and labels show the changes from pointings -2 through +2. Note observational constraints mean there are less +2 pointings.

Fig. 2 demonstrates a known effect of estimating the PS from gridded visibilities, where the power inferred is lower where the density of gridded visibilities is higher, as observed by Barry et al. (Reference Barry, Beardsley, Byrne, Hazelton, Morales, Pober and Sullivan2019a, see Appendix A). Note we have not corrected for this in Fig. 2, which is normally done by multiplying by a factor two (as found by Barry et al. Reference Barry, Beardsley, Byrne, Hazelton, Morales, Pober and Sullivan2019a), to explicitly show the effect.

Visibilities are gridding using a kernel which is designed to minimise aliasing and match the area of the MWA primary beam (Trott et al. Reference Trott2016). Given the array layout, and the need to match the kernel to the primary beam, the kernel has an unavoidable footprint large enough to overlap with neighbouring gridded visibilities. Depending on the covariance between neighbouring visibilities, and the density of the gridded visibilities, the power estimated will therefore vary as a function of u,v gridding location. Using simulations, Barry et al. (Reference Barry, Beardsley, Byrne, Hazelton, Morales, Pober and Sullivan2019a) found a factor of two correction was sufficient as a normalisation factor for modes of interest. This estimate was made using the FHD/ $\varepsilon$ ppsilon pipeline and a model of the MWA primary beam as a gridding kernel, whereas CHIPS uses a Blackman-Harris kernel. This estimate was also made using 2 s resolution data, rather than 8 s. Given these differences, we investigate the gridding density correction factor in Section 4.1 using CHIPS.

4.1 Gridding density correction factor

The gridding density is directly affected by the u,v coverage, which in turn is dictated by the array layout and phase centre. The covariance of neighbouring u,v points will therefore be affected by LST (as the phase centre is always set to EoR0 field centre). The primary beam pointing will also have an effect as this changes the amplitude and spectral behaviour of the visibilities. To thoroughly investigate the gridding density we therefore simulate a number of observations spanning realistic ranges of LST and pointings. To do so, we take the first (subset LST1), central (subset LST2), and final observation (subset LST3) across five different pointings (−2 to +2) from three real nights of MWA phase I EoR observing (see Fig. 3). For the LST2 subset, we select the observation where the primary beam pointing is closest to the EoR0 field centre, to maximise the beam coverage over the sky model coverage. As the visibilities are normally averaged to $8\,$ s, $80\,$ kHz we simulate at this cadence to save compute time. Each observation takes a total of $\sim$ $32\,$ GPU hours. For 45 simulated observations this is a total of 60 GPU days. We use the Pawsey Garrawarla clusterFootnote h and split each simulation across 24 GPUs, meaning these simulations take less than 3 days real-time.

By running CHIPS on various combinations of these observations, we can investigate both primary beam pointing and LST effects. All integrations over multiple observations are done coherently, that is, all observations are gridded to the same u,v grid, which is then used to estimate the PS. We take the median of the ratio of the observed PS to the expected value as a function of $k_\perp$ (i.e. take the median along the y-axis of the middle panel in Fig. 2) to see the effects of gridding, as shown in Fig. 4. We take the median, rather than the mean, as the distribution of ratios along each $k_\perp$ bin display significant skew.

Figure 4.

The median ratio of the recovered noise PS to expected value, as a function of $k_\perp$ , for: (a) the five different pointings, each for a single observation; (b) the three different LST subsets integrated over 5 observations; (c) the three different LST subsets integrated over 15 observations; (d) all three LST subsets integrated together over 15, 30, and 45 observations. The horizontal dashed line shows the mean of the median ratios for $k_\perp < 0.1$ as measured from the bottom right plot. This dashed line Any dataset which is a multiple of five observations has an even split across the five pointings.

Fig. 4a shows that for a single observation, the estimated power is consistent across pointings, with lower power estimated a low $k_\perp$ (higher gridding density). This shows the beam volume correction applied internally in CHIPS returns a consistent power level across the five pointings. Fig. 4b shows that the power recovered from the three LST subsets over 5 observations is roughly consistent, but the power recovered does vary somewhat. This recovered power depends on the exact u,v coverage, and the way the primary beam interacts with the sky model. Fig. 4c shows that when integrating over 15 observations for each LST subset, the recovered power at low $k_\perp$ is consistent with the power recovered when integrating over 5 observations. However, the power recovered at higher $k_\perp$ is actually higher than at when integrating over 5 observations. This highlights the fact that the power estimated comes from a combination of not just the amount of data gridded, but the covariance between the neighbouring gridded visibilities. This is further highlighted by showing that integrating the three LST subsets together across 15 observations yields significantly lower power at high $k_\perp$ . High $k_\perp$ is derived from longer baselines, where the drift of a baseline u,v coordinate over time is greatest. The LST1 and LST3 groups are also closely separated in time, yielding close gridding locations, but overlap as the beam pointing changes. The amplitude and the spectral behaviour of the primary beam therefore sharply changes on the sky, reducing the covariance between neighbouring gridded visibilities. This results in the lower power recovered. Fig. 4d shows that when integrating more and more data with the same combination of LST coverage and pointings, there is no change in power at low $k_\perp$ .

Given the stability of the recovered power when combining all LST subsets with an integration of 15 observations or more, one could fit a functional form to derive a correction. We investigate this approach in Appendix A and find little to no difference in using a fit when compared to just applying a scalar normalisation factor. PS upper limits are derived from integrating hours of data to overcome thermal noise, and modes where $k_\perp > 0.1\,h\,\mathrm{Mpc}^{-1}$ are typically cut when using the AusEoRPipe. To obtain a normalisation factor for the gridding density we therefore simply take the mean of the median observed ratios for the integrations over all LSTs and pointings where $k_\perp < 0.1\,h\,\mathrm{Mpc}^{-1}$ . This value is plotted as the horizontal grey dashed line in Fig. 4. Inverting this gives a normalisation factor of 2.651. We apply this normalisation and create a number of 1D PS for various integrations, shown in Fig. 5, and compare the outcomes to the previous correction factor of two. This shows that the new normalisation factor causes too much power to be recovered for a single observation, and when only using the LST2 subset, but recovers the correct power to within 10% for most k modes when integrating over all LSTs and pointings. The factor two correction (coincidentally) does well for a single observation, but does not recover the correct power for higher integrations.

Figure 5.

1D PS for various integrations of the Gaussian noise simulation, where (a) a factor of two was used to normalise for gridding density and (b) a factor of 2.5845 was used. The vertical dashed line shows the maximum $k_\perp$ -mode.

5. Recovered 21-cm signal

To test the 21-cm model detailed in Section 3, we simulate 30 observations, matching the first two nights detailed in Fig. 3, with the same parameters as used for the Gaussian noise simulation (detailed in Section 4). We leave out simulating the third night purely to save on computational resources. To qualitatively assess the accuracy of the model through WODEN and into CHIPS, we compare a CHIPS 2D PS to one derived directly from the lightcone box, shown in Fig. 6. We produce the 2D PS directly from the lightcone using the OSIRIS packageFootnote i (Cook, Trott, & Line Reference Cook, Trott and Line2022), where no primary beam or instrument sampling was applied, meaning the entire u,v plane was sampled. We see that the CHIPS PS is broadly consistent with the OSIRIS PS. Noticeably, it seems there is a potential signal loss at high $k_\parallel$ in the CHIPS PS (the darker blue region at the top of Fig. 6b). There is also excess power seen along high $k_\perp$ , around the dashed black line, which is expected as the instrument baseline sampling causes mode mixing, moving power up from lower to higher $k_\parallel$ .

Figure 6.

2D PS from the 21-cm simulation where: (a) produced using OSIRIS directly from the lightcone box from Greig et al. (Reference Greig, Wyithe, Murray, Mutch and Trott2022); (b) a CHIPS 2D PS made from integrating over 30 simulated observations; (c) the same PS from (b) but with standard cuts made to remove foreground contamination. The solid black line indicates an estimate of the horizon, and the dashed black line the full-width half-max of the MWA primary beam.

To test the gridding normalisation, we apply it to various integrations and produce 1D PS, shown in Fig. 7. We compare them to an expected signal, again derived directly from the 21-cm lightcone box. Excellent agreement is shown, and similarly to the Gaussian noise simulation, the gridding density normalisation estimates too much power for short integrations, but renders consistent results for longer integrations with the same mix of LSTs and pointings.

In Fig. 8 we show the ratio of the final recovered 21-cm signal to the expected value, for a 60-min integration. While good agreement is shown, there does seem to be signal loss around $0.1 < k < 1$ , which was not seen for the Gaussian noise simulation. This is likely due to the frequency interpolation of the original lightcone as detailed in Section 3. For reference, Fig. 8 shows the maximum $k_\perp$ mode that went into the average from 2D to 1D. k modes above this are binned from only high $k_\parallel$ modes, which is where signals with the most spectral structure manifest. As this area is where the signal loss is found, it is likely the frequency interpolation is causing this. As the Gaussian noise sky model has no correlation over frequency, it should not suffer from this interpolation, which is why it doesn’t display the same signal loss.

Figure 7.

Recovered 21-cm signal from various integration lengths of simulated observations, when normalised with a factor 2.5845. The panel bottom right is included to illustrate the normalisation effects between $0.1 < k < 1.0\,h\, \mathrm{Mpc}^{-1}$ . The two minute data set is for a zenith pointing; the 10 minute data set is for the LST2 subset; all other data sets are an even split between the five pointings and LST subsets as described in Section 5.

Figure 8.

Ratios of the recovered 21-cm signal to the expected value, for a 60 min integration. The orange with crosses line shows the ratio when using the factor two gridding density correction, and the blue with circles when using the fitted correction. The vertical dashed line shows the maximum $k_\perp$ -mode, meaning any k-mode above this was derived purely from the Fourier transform of the visibilities along frequency.

Fig. 8 also shows the effect of performing the ‘wedge-cut’, where modes normally dominated by foreground sources are removed. While broadly consistent, the signal loss seems to worsen after the cuts. This is likely due to the reduced number of bins going into the average from 2D to 1D. As can be seen from Fig. 6c, the horizon cut (marked by the solid black line) removes more samples for smaller $k_\parallel$ modes. However, lower sampling is not the only possible reason. The manifestation of the primary beam, as well as the exact u,v distribution, could both cause a systematic bias at specific spots in the 2D PS, which could contribute here.