Characterization of the AARTFAAC-12 aperture array: radio source counts at 42 and 61 MHz Open Access

ABSTRACT

1 INTRODUCTION

Recent developments in the area of signal processing in astronomy as well as the declining cost of high-performance (computing) hardware has enabled the advent of (dense) aperture array instruments in radio astronomy (Zarka et al. 2012; van Haarlem et al. 2013; Taylor & LWA Collaboration 2014; Hallinan et al. 2015; DeBoer et al. 2017; Beardsley et al. 2019, LOFAR, MWA, HERA, LWA, NenuFar, OVRO-LWA). This is especially true for low frequencies (up to a few GHz). At the lowest frequencies (tens to hundreds of MHz), the arrays use simple dipoles as receiving elements which synthesize an aperture (hence aperture arrays) to various degrees, by sampling the electric field of the incoming electromagnetic radiation and employ simple receivers for initial signal conditioning in the field. These systems usually have very large fields of view (FoV) up to |$2\pi \, \mathrm{sr}$|⁠, depending on the antenna pattern of the individual receiving element; it can be modified by combining elements together (beam forming) if so desired, diminishing the FoV, but gaining in sensitivity.

The LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) telescope uses dipoles sensitive to frequency ranges of 10–90 MHz (low band, LBA) and 110-240 MHz (high band, HBA) grouped into stations which are then correlated with each other or their signals added (depending on the science case). In another incarnation of (part of) the same instrument, demonstrating its modular design, we have taken the LBA dipoles of the central twelve stations, extracted the dipole signals out of the LOFAR signal path before station beam forming and correlated them with each other. Thus, we have obtained a dense array, with a very large FoV, the Amsterdam-Astron Transient Facility And Analysis Center (AARTFAAC-12, henceforth A12).

A12 is able to observe commensally with LOFAR, whenever it is observing in LBA mode (we usually use the LBA_OUTER mode in which a ring of 48 outermost LBA dipoles in each station is active). We use a separate correlator to produce visibilities on baselines between each of the dipoles comprising the 12 core stations, and can form up to 16 sub-bands (SBs), each with a bandwidth of 195.3 kHz across the LBA band. These SBs can be placed at any desired frequency in various configurations (continuous or dispersed). Each SB can be subdivided into up to 64 channels, but the correlator performance starts to deteriorate if a larger number of channels per SB are recorded. For more details on the system configuration and signal path, we refer the reader to Prasad et al. (2016).

The scientific use cases of the A12 system are varied; its sensitivity to diffuse extended emission renders it very suitable for Galactic science. Its large FoV and the low frequencies at which it operates make it a powerful tool for space weather studies (Mevius et al. 2016) as well as measuring the power spectrum of Cosmic Dawn (z∼18; Gehlot et al. 2020). Since the instrument images the whole sky at high temporal cadence (one second), it is also suitable for radio transient detection, enabling an order of magnitude improvement in sensitivity compared to its predecessor, AARTFAAC-6 (Kuiack et al. 2019).

To utilize the instrument properly in imaging mode, which is its main mode of operation, we need to characterize its systematics, and assess its performance as well as that of the associated data processing pipeline. In this work, we will describe the A12 calibration and imaging results as well as image noise properties using standard LOFAR processing tools (Section 2) and go on to derive source counts using a typical A12 data set (Section 3). We discuss our results in Section 4 and give concluding remarks in Section 5.

Even though our source counts are limited to brighter sources, the number of such studies at extremely low frequencies are limited (eg. Lane et al. 2014). Also, since we derive our source counts from a single observation covering a large area of the visible sky at that moment, our results have fewer observational systematics compared with surveys resulting from multiple pointings to cover the same sky area.

2 DATA REDUCTION

We have recorded visibilities for fifteen minutes starting at 01:29:27.5 UT on 18 February 2019. Each SB has 3 channels and we have placed one contiguous group of eight SBs around 41.7 MHz and another 8 SB group around 61 MHz, at the maximum sensitivity of the LBA dipoles. The integration time was one second. A12 observes in drift scan mode, pointing at zenith with a FoV extending over the entire visible sky.

The A12 correlator produces visibilities in near-real time; to facilitate studies of radio transients it was designed to minimize latency using nominal observing parameters. The visibilities are written on disk in a way compatible with the AARTFAAC-6 reduction and imaging pipeline which produces near-real time (up to three seconds lag time) images for transient detection (Prasad et al. 2016). The recording format is a streaming format; to facilitate calibration using LOFAR data reduction packages, we convert the visibilities to measurement set (MS) format using the aartfaac2ms tool (Offringa et al. 2015). During the conversion, we flag the visibilities to mitigate RF interference using the aoflagger tool (Offringa, van de Gronde & Roerdink 2012), flag for bad dipoles (based on auto correlation values) and phase shift the visibilities to the coordinates of zenith at the mid-point of the observing run (in this case α: 11^h55^m01.2^s, δ: +52^○50′17.43″).

Before initiating the calibration, we concatenate each group of eight individual single-SB measurement sets into a single measurement set comprising 24 frequency channels. The calibration procedure consists of the following steps:

• Solving for the direction independent (DI) complex gains of each receiving element per integration time. The beam pattern of the elements is taken into account. We use CasA and CygA as calibrators to set the flux scale, using source models composed of Gaussians as well as points, set on the Scaife and Heald flux scale (Scaife & Heald 2012). We find that using the brightest (visible) 3C sources as calibrators does not result in usable calibration solutions; they have too low S/N, even after using solution intervals of a few minutes.

  In the DI case, the problem (per baseline) is cast as (Smirnov 2011):

  $$\begin{eqnarray*} V_{\mathrm{pq}} = G_{\mathrm{p}}B_{\mathrm{p}}M_{\mathrm{pq}}(B_{\mathrm{q}})^{H}(G_{\mathrm{q}})^H , \end{eqnarray*}$$

  (1)

  where V_pq are the recorded visibilities on the baseline between dipoles P and Q, G are the dipole (complex) gain Jones matrices, B is the dipole beam Jones matrix and M are the model visibilities on that baseline. H denotes the Hermitian (complex) transpose operator.
  We solve for the amplitudes and phases of all four polarization products using the full-jones mode in dppp, a package which includes various solvers and other tools for MS manipulation, one of the main components of the LOFAR data processing pipeline. The solution interval is set to equal the integration time (one second) and we solve per channel. In the full-jones solve mode, the gain Jones matrix per dipole is:

  $$\begin{eqnarray*} G_{\mathrm{p}} = \begin{pmatrix}A_{\mathrm{xx}}e^{\phi _{\mathrm{xx}}} A_{\mathrm{xy}}e^{\phi _{\mathrm{xy}}} \\ A_{\mathrm{yx}}e^{\phi _{\mathrm{yx}}} A_{\mathrm{yy}}e^{\phi _{\mathrm{yy}}} \end{pmatrix}, \end{eqnarray*}$$

  (2)

  here XX, XY, YX, and YY are the four instrumental (linear) polarizations of the recorded visibilities. A and ϕ stand for the amplitude and phase of the gains.
• Direction dependent (DDE) solve (again using the full-jones solve mode) is performed using CasA, CygA, TauA, and VirA source models (depending on source visibility and elevation for the duration of the observing run) on the Scaife and Heald flux scale. The dipole beam pattern is taken into account, and the solve is done per one second and channel. The DDE solutions are used to subtract the sources solved for from the visibilities. In this case:

  $$\begin{eqnarray*} V_{\mathrm{pq}} =\sum _{\mathrm{dir}}G_{\mathrm{p}}^{\mathrm{dir}}B_{\mathrm{p}}M_{\mathrm{pq}}(B_{\mathrm{q}})^{H}(G_{\mathrm{q}}^{\mathrm{dir}})^H, \end{eqnarray*}$$

  (3)

  i.e. we have gain solutions per each direction towards a particular calibration source.

We average the calibrated data by a factor of four in frequency (sufficient to avoid bandwidth smearing out to the edge of the FoV) and scale the visibilities with the dipole beam pattern in the direction of the phase center to get intrinsic flux values. Next, we image the averaged data set using wsclean (Offringa et al. 2014). The image scale is 6 and 3 arcmin per pixel at 41.7 and 61 MHz respectively, and we image the full FoV using robust 0 weights while correcting for the dipole beam during the imaging process. Multi-scale cleaning is used with scales of 0, 10, 30, and 60 pixels. We have specified 10⁶ clean iterations, letting the clean continue until the stopping criteria determined by the wsclean parameters auto-mask and auto-threshold were reached. The σ values were set to 3 and 0.3 respectively for these parameters. We produce Stokes I and V images without and with a uv taper of 10λ (corresponding to angular scales larger than six degrees) and removing the diffuse Galactic synchrotron emission to determine the influence it has on our analysis. The images have PSF sizes of 24 × 18 arcmin and 17 × 11 arcmin at 41.7 MHz and 61 MHz respectively. The dipole sensitivity drops sharply near the horizon, which defines our total field of view (FoV). One of the images (non-tapered, Stokes I at 61 MHz), masked starting at the end of the FoV (below an elevation of 15^○) is shown in Fig. 1 for reference.

Figure 1.

The FoV of the A12 array at 61 MHz (Stokes I, no uv taper). The image extends to 15^○ elevation, centered on the zenith J2000 coordinates at mid observation. The Galactic plane is visible along the image edge, with the Galactic spur at bottom left. The PSF is 11.4 × 16.8 arcmin A dashed circle marks the 52^○ elevation limit of the central image cutout used in the source count analysis, and the thin filled line marks the horizon.

Open in new tabDownload slide

2.1 Flux scale

To check whether the flux scaling in the image was set correctly in the initial calibration step, we have used pybdsf (Mohan & Rafferty 2015) and extracted all the sources detected in a circular region with a radius of 38^○ centered on zenith. The source detection settings have been set to the default values (3σ and 5σ for the island and peak detection thresholds respectively) and used a 2D rms map for the background (rms_map=True) with an rms_box of 30 pixels and a step of 10 pixels used to calculate the rms map. Since A12 has a large synthesized PSF (see Table 1), and is sensitive to extended emission, we have matched our source catalogs to the 38 MHz 8C survey (4.5 arcmin PSF; Hales et al. 1995) and the 60 MHz AARTFAAC-6 survey catalog (60 arcmin PSF; Kuiack et al. 2019). We have used a match radius of 5 and 17 arcmin at 41.7 MHz and 61 MHz respectively and did not scale the comparison catalog fluxes due to the proximity in frequency. From the matched catalogs, we have selected the point sources as follows. At 41.7 MHz, we have selected the sources for which their A12 peak to total flux density ratio was greater than 0.5 and their size parameter in the 8C catalog was between 0.95 and 1.2. At 61 MHz, we have selected the sources for which their A12 peak to total flux density ratio was greater than 0.85. Then, to determine the flux scale correction for each image, for the matched and selected sources we have computed the peak flux density ratio between the 8C catalog and the 41.7 MHz A12 catalog as well as the total flux density ratio between the AARTFAAC-6 survey catalog and the 61 MHz A12 catalog. The flux scale correction is taken to be the mean of the computed ratios The correction factors are 0.72 ± 0.03 and 0.98 ± 0.07 at 41.7 MHz and 61 MHz respectively (as shown in Fig. 2) and we use them to scale the A12 catalog flux densities at the corresponding frequencies in the subsequent analysis.

Figure 2.

Source flux density ratio (8C survey and AARTFAAC-6 catalog over A12 measured flux densities at a given frequency) across the analysis field. Plotted values are for point sources as described in the text. The black dashed line marks the mean value of the plotted data points, while the red area denotes the standard error for the sample.

Open in new tabDownload slide

3 RESULTS

We limit our analysis to the regions of the FoV free from Galactic emission with the highest intensity. One can model this emission and subtract it from the data, however due to the complexity of the model, we have deferred this approach for the time being and decided to use a UV taper as mentioned previously.

The local sidereal time of the observation limits the brightest of the Galactic emission to lower elevation; we have masked the relevant region of the image and will use an area centered on zenith with a radius of 38^○ in our further analysis. This area is shown in the top row of Fig. 3 and the corresponding r.m.s. images (with and without UV taper) after performing source extraction with pybdsf are given in the middle row.

Figure 3.

Image cutout of the central regions, where the analysis is performed. Top row: 61 MHz images without and with (right) UV taper. Mid row: the corresponding r.m.s. images produced by pybdsf after source extraction is performed on the image shown in the top panel. The contours are on a linear scale going from 0.6 to 2 Jy/PSF. Bottom: UV tapered image rms at 41.7 MHz. Contours are drawn on a linear scale between 0.1 and 6 Jy/PSF.

Open in new tabDownload slide

3.1 Image noise

For the theoretical thermal noise under these assumptions, the obtained value is: σ_th = 0.06 and 0.05 Jy/PSF at 41.7 and 61 MHz respectively. The Stokes V image noise is a good estimator of the theoretical thermal noise; we have produced Stokes V images on various time-scales and confirmed that the Stokes V image noise scales with the observing time as expected (⁠|$\frac{1}{\sqrt{\Delta t}}$|⁠). For the observing time we use in this analysis (t = 15^m), the measured Stokes V image noise is 0.69Jy/PSF and 0.33Jy/PSF at 41.7 and 61 MHz respectively. The excess over the theoretical values may be the result of flux leakage from Stokes I to Stokes V due to element beam model systematics.

We use the ‘probability of deflection’ analysis method (Scheuer 1957) to better understand the noise properties of the images. Detailed description of the P(D) method can be found in Vernstrom et al. (2014) and Franzen et al. (2019), and here we present an overview of the procedure we have implemented. We convert the images to Zenithal equal area (ZEA) projection using the montage package before doing the analysis, to ensure that distortions in the image plane are handled properly. For a given source count model, we compute its probability density function, i.e. the P_source(D), taking the dirty beam into account. The source count distribution is taken to be a polynomial fit to the 154 MHz LOFAR and MWA source counts scaled to our observing frequencies (Franzen et al. 2019) using a spectral index of −0.8. Then, we assume additive Gaussian (thermal) noise contribution which we estimate by measuring the standard deviation of the Stokes V images. Ideally, they should be only affected by thermal noise, assuming that the sky signal is unpolarized at the scales we sample and that there is no polarization leakage. Since the random variables describing the thermal noise and the signal produced by the source counts mentioned above are assumed to be independent, the combined probability distribution is a convolution between the probability distributions of the noise and the source model: P(D) = P_n(D)*P_source(D). We compare the convolution result with the measured probability density function obtained from the pixel values of the Stokes I image, P_obs(D). The resulting plots are shown in Fig. 4.

Figure 4.

Left: P(D) distributions at 41 MHz. The filled line Gaussian centered at zero represents the image noise probability distribution function, measured in a region of the UV tapered Stokes V image. The red dashed line represents the P(D) function of the convolution between the noise and source model (orange) probability distributions. The remaining function, marked in blue, represents the smoothed probability distribution function of the pixels in the Stokes I image. Right: The same quantities calculated for the corresponding 61 MHz images. Note that the measured probability distribution matches the theoretically derived one.

Open in new tabDownload slide

We derive a Stokes I image (confusion) noise of 2.6 Jy/PSF and 0.9 Jy/PSF at 41.7 and 61 MHz respectively by dividing the interquartile range with 1.349, i.e. the rms of a Gaussian distribution. At the higher frequency this value matches nicely the value estimated by convolving the SKADS source model confusion noise with the thermal noise (the width of the red dashed Gaussian in the right-hand panel of Fig. 4). In the lower frequency band, there is an excess measured image noise, most likely coming from residual sidelobes of bright sources. The theoretical and measured noise values as well as detected source numbers (Cohen 2006) for both observing frequencies are summarized in Table 1.

3.2 Source counts

As an initial step in deriving the source counts, we have extracted sources from the (ZEA converted) image (covering 3958 square degrees, 1.2 sr) using the pybdsf source finder. We used the default parameters of the process_image task, hence the island and pixel detection thresholds were set to 3σ and 5σ respectively. We have adopted non-default values only for the rms_box which we have set to a size of 30 pixels. Also, we have selected to use a 2D background rms map.

We correct the total flux density entries in the catalog by the factors derived in Section 2.1; we do not correct for ionospheric smearing of the peak flux densities due to the size of our PSF.

Next, we bin the sources in twelve logarithmic spaced flux density bins and compute the A12 source counts using the procedure outlined in Lane et al. (2014) and Williams et al. (2016). The raw counts per bin are weighted by the image area (expressed in steradians) associated to the sources in the corresponding bin. The computation of the areas proceeded as follows. For each source in a given bin, we have summed the number of pixels in the RMS image (output by pybdsf after the source extraction) which have flux densities five times smaller than the source peak flux density. We then sum over all the sources in the bin, and knowing the pixel scale, obtain the area in steradians. We further Euclidean normalize the counts by multiplying the area corrected counts per bin by S^2.5 where S is the flux density of the bin center.

The completeness (probability that all of the sources above a given flux density are detected) of the catalog was obtained by performing a Monte Carlo simulation; ten images were generated by inserting 2000 sources per image at 41.7 MHz (and 4000 at 61 MHz) at random locations in the residual map output by pybdsf during the source extraction mentioned above. These sources were point sources simulated as 2D Gaussians using the restoring PSF size of the original image. Their brightness was randomly drawn from a power law distribution with S_min = 0.2Jy, S_max = 50Jy and a slope of α = −0.6. Out of all injected sources, 80 per cent were point sources. A pybdsf source extraction was performed on each of the simulated images using the same parameters as the ones used during the master catalog creation. The resulting source catalogs (from the simulated images) were matched with the input source catalog per image using topcat. These matched source catalogs were used to compute the source detection fraction, the fraction of sources actually detected per flux bin, accounting for the areas over which the source finding was performed. The completeness for a given flux density was calculated by integrating the detection fraction upwards of that flux density. The detection fraction and completeness are shown in Fig. 5. Our source catalog is 97.5 per cent complete above 64.8 Jy at 41.7 MHz and above 50.6 Jy at 61 MHz.

Figure 5.

Top row: source detection fraction (left) and catalog completeness (right) at 61 MHz. The red line shows the area-corrected detection fraction. Dotted lines show its 1σ uncertainty. The vertical dashed line show the 5σ rms noise value of the Stokes I image and the dotted vertical lines give its 2σ bounds. Bottom row: source detection fraction and completeness at 41.7 MHz.

Open in new tabDownload slide

Since our PSF is large, source blending is an issue. We have derived a blend correction factor as follows. We randomly draw 2000 source fluxes between 5 and 150 Jy at 41.7 MHz (4000 source fluxes between 1 and 150 Jy at 61 MHz) from a power law distribution (as was done for the completeness correction) and inject them into the Stokes I image. We perform source extraction on these images with the same parameters as before, and obtain simulated catalogs. Source matching these catalogs to the input catalog, we get to the blend correction, by looking at the fraction of simulated sources detected close to real sources per flux density bin (Franzen et al. 2019). We have also computed the reliability of the catalog and its associated false detection rate and found that the associated correction is negligible, so we have not corrected the counts for it.

The derived source counts have been corrected for completeness, by multiplying them with the (area corrected) detection fraction (shown in the left-hand panel of Fig. 5) and for blending by multiplying them by the derived correction factor.

The source counts are shown in Fig. 6 along with scaled 6C survey (Hales, Baldwin & Warner 1988) and VLSSr survey derived counts as well as a SKADS model-derived source counts. Details are given in Table 2.

Figure 6.

The derived Euclidean normalized A12 source counts (orange step line) measured at 41.7 MHz (left) and 61 MHz (right), shown with scaled 6C and VLSSr survey source counts (points). A scaled LOFAR and MWA (SKADS) source count model at the observing frequencies using a spectral index of −0.8 is also plotted for comparison. The A12 count errors are Poissonian (marked with a shaded blue region), and the flux density errors correspond to the bins used for deriving the counts. Details are given in Table 2.

Open in new tabDownload slide

3.3 Spectral index

From the source catalogs at both frequencies, we have computed the two-point spectral index (⁠|$\alpha =\frac{\mathrm{log(S_{1}/S_{2})}}{\mathrm{log}(\nu _{1}/\nu _{2})}$|⁠, where S₁, S₂ and ν₁, ν₂ are the integrated flux densities at the respective frequencies), by cross-matching the sources using a match radius of 5 arcmin. The resulting spectral index distribution is given in Fig. 7; we find a mean spectral index value of |$\alpha _{41.7}^{61} = -0.78$|⁠.

Figure 7.

Spectral index distribution. The vertical dashed line represents its mean value, |$\alpha _{41.7}^{61}=-0.78$|⁠.

Open in new tabDownload slide

4 DISCUSSION

From the image noise properties listed in Table 1, we see that the measured Stokes I image noise at 61 MHz (estimated from the image during the P(D) analysis as well as arrived at by convolving the source model P(D) distribution with the thermal noise distribution) is in agreement with the theoretically derived confusion noise value. We conclude that A12 is confusion noise limited at 61 MHz at 0.9 Jy/PSF. The analysis for 41.7 MHz shows that there is an excess noise component above the theoretical confusion noise limit of around 0.6 Jy/PSF most likely due to Galactic emission residuals and/or sidelobe noise from imperfect cleaning from bright sources. At this frequency, the instrument noise floor in Stokes I is 2.6 Jy/PSF. The noise floor affects the faint end of our source counts, which start at around 8σ in the respective bands.

The A12 source counts (Fig. 6) are in general agreement with the 6C and VLSSr survey counts, but go on to much larger flux density values. At lower flux densities the source counts show a drop off due to incompleteness. We noticed that at 61 MHz the derived source counts have an excess compared to the 6C and VLSSr survey and SKADS model counts for flux densities in the range between 20 and 80 Jy. We ascribe this to the image noise morphology (Fig. 3) at the spatial scales and flux densities which affect the source count normalization. To correct for this systematic, we have convolved the noise map with a Gaussian kernel (100 pixel wide and having a support of 199 pixels) and used the corrected noise map for the computation of the counts shown in the right-hand panel in Fig. 6.

The computed spectral indices for the cataloged sources have a peaked unimodal distribution with a mean value of |$\alpha _{41.7}^{61} = -0.78$| which is to be expected from a population of extragalactic sources, consisting predominantly of active AGN.

Following the work of Hardcastle et al. (2021) we estimate the brightness of the cosmic radio background (CRB), as reported by Fixsen et al. (2011) using the A12 source counts we derived. By integrating our source counts (corrected for completeness), we get a value of 153 K at 41.7 MHz and 133 K at 61 MHz. Both of these values are below the values predicted by the power law model proposed by Fixsen et al. (2011) (1814 K and 4909 K respectively when extrapolated to our observing frequencies) and below the values obtained if we extrapolate the values found by Hardcastle et al. (2021) at 144 MHz using a T ∼ ν^β, β = −2.7 model (1249 K and 447 K respectively). Similarly, Subrahmanyan & Cowsik (2013) and Dowell & Taylor (2018) find higher (extrapolated) values compared to what we report. On the other hand, Vernstrom, Scott & Wall (2011) have used source counts at seven different frequencies spanning the range from 150 MHz to 8400 MHz to derive a power law dependence for the CRB, which gives a value of 18 K at 150 MHz using their derived power law index of β = −2.28. At our observing frequencies, their derived dependence gives the values of 333 K and 140 K at 41.7 MHz and 61 MHz respectively, roughly matching our measurement at our higher observing frequency. Our findings suggest that any extended emission on arcminute scales to which we are sensitive (for example, extended diffuse emission associated with the sources we detect but resolved out by other surveys) does not significantly contribute to the proposed background. While fainter sources detected in other (higher frequency) surveys can partially account for the difference between the CRB values we derive from our source counts and other studies, it is more likely that the cause can be ascribed to residual systematics (especially pertaining to the lowest frequencies we probe) coupled with a possible low frequency turnover in the power law fits to existing source counts.

5 CONCLUSIONS AND FUTURE PROSPECTS

The AARTFAAC-12 aperture array is a very capable instrument, delivering confusion limited images which we have analyzed in this work and shown that they are applicable in various science cases. The instrument is particularly suitable for transient science as well as imaging large scale diffuse emission, necessary for updating existing sky models at the lowest radio frequencies. In this work we have:

• Described in detail the procedure used to calibrate and image A12 data.

• Characterized and validated the images, using P(D) analysis to derive noise parameters. This represents a crucial precursor for science exploitation. A12 is confusion limited at around 0.9 Jy at 61 MHz.

• Derived source counts at 41.7 and 61 MHz, which are in agreement with previous studies and which showcase the ability of the instrument to perform large scale surveys on a short time-scale.

• Estimated the spectral index between our observing bands and used the source counts to characterize any excess radio emission at the frequencies observed.

The A12 calibration and imaging pipeline is modular and it will be modified to speed up processing, leading to a near-real time transient detection pipeline with an imaging cadence of around one second. The instrument will evolve as the LOFAR telescope is upgraded to its more capable successor, LOFAR2.0.

ACKNOWLEDGEMENTS

The authors would like to thank the anonymous referee for the constructive comments which have helped to improve this work. AARTFAAC development and construction was funded by the European Research Council (ERC) under the Advanced Investigator grant no. 247295 awarded to Prof. Ralph Wijers, University of Amsterdam. We thank The Netherlands Institute for Radio Astronomy (ASTRON) for support provided in carrying out our observations. AARTFAAC is maintained and operated jointly by ASTRON and the University of Amsterdam. WLW also acknowledges support from the CAS-NWO programme for radio astronomy with project number 629.001.024, which is financed by the Netherlands Organisation for Scientific Research (NWO). BKG acknowledges the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 884760, ‘CoDEX’) and the National Science Foundation through an award for HERA (AST-1836019). We would also like to thank the LOFAR telescope operators and science support group for their assistance in obtaining the data used in this work. We use data obtained from LOFAR, the Low Frequency Array designed and constructed by ASTRON, which has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative (Grant No. EINF-245). This research made use of montage. It is funded by the National Science Foundation under Grant Number ACI-1440620, and was previously funded by the National Aeronautics and Space Administration’s Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. This research made use of astropy¹ a community-developed core python package for Astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), as well as: numpy (van der Walt, Colbert & Varoquaux 2011), scipy (Virtanen et al. 2019) , matplotlib (Hunter 2007) and topcat (Taylor 2005). Accordingly, we would like to thank the scientific software development community, without whom this work would not be possible.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES