© A. M. Kutkin et al. 2022

1 Introduction

Large astronomical surveys have always played a key role in helping characterize and understand the objects populating our Universe. Continuous technological developments have facilitated deep observations over a large field of view, further boosting the possibilities for producing new, large, and deep surveys. At radio frequencies, this growth in the number of surveys is particularly evident, triggered by the plethora of new instruments and telescopes developed to prepare for the upcoming Square Kilometer Array (SKA, e.g., Braun et al. 2019). Most importantly, in the great majority of the cases, the related source catalogs have been made available online, providing the astronomical community with invaluable tools for carrying out science that has not been possible thus far.

Among the SKA pathfinders that have started producing new surveys, there is Apertif (AperTIF – Aperture Tile In Focus), the phased array feed (PAF) receiver working at 1400 MHz and installed on 12 of the dishes of the Westerbork Synthesis Radio Telescope (WSRT). It provides forty instantaneous overlapping beams (in the case of a phased array feed, also referred to as compound beams, CB), creating a wide field of view of about 10 square degrees. These features, coupled with a wide bandwidth of 300 MHz are what make Apertif an excellent instrument for surveys. A backend with high spectral resolution (12.2 kHz) provides possibilities for simultaneous spectral line, continuum, and polarization surveys. The full system and its capabilities are described in van Cappellen et al. (2022).

Apertif not only represents an excellent proof of concept for the future generation of phased array telescopes, it also provides valuable scientific data in terms of sensitivity and angular resolution and it is a significant improvement over the two major surveys of the northern sky at 20 cm performed more than 20 yr ago: the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) and VLA Faint Images of the Radio Sky at Twenty Centimeters (FIRST; Becker et al. 1995). The former offers a good surface brightness sensitivity of 0.45 mJy beam⁻¹, but a low angular resolution of 45″; while the latter is even more sensitive (0.15 mJy beam⁻¹), but resolved the extended structures of many radio sources with its 5″ beam size. The Apertif survey provides an intermediate angular resolution of ~12″ × 12″/sin(Declination), together with a noise level reaching 20 µJy beam⁻¹ in its most sensitive single-pointing images. An overlap with the above-mentioned surveys makes the Apertif continuum data a great platform for a plethora of scientific questions; in particular, a separation in time between the Apertif and NVSS surveys of about 25 yr provides an opportunity for detecting long-term variable radio sources.

It is important to note that the resolution and sensitivity of the Apertif survey are sufficient to complement the breakthrough LOFAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2017, 2019, 2022). With the LOFAR surveys aiming to cover the entire northern sky at 60 and 150 MHz, there is a perfect synergy between these and the Apertif surveys. This has been already illustrated by the studies presented, for example, in Morganti et al. (2021a,b).

The potential that Apertif presents is further enhanced by the fact that the raw data (as well as the data products) are made publicly available to the scientific community. The data release of the first year of Apertif observations (Apertif Data Release 1, ADR1) took place in November 2020¹ (Adams et al. 2022, A22 hereafter). Processed data products for every CB were released on the basis of the continuum image quality and all released data products include validation measures. In this paper, we complement the continuum images with the extracted source catalog.

A key requirement for a source catalog that has scientific applicability is a correct flux scale. For PAFs, it is a challenging task to measure the amplification diagram, namely, the primary beam response. In Apertif, the 40 CBs are formed electronically by combining and weighting the signals received by the antenna elements installed in the focal plane of every dish (Vivaldi antennas in the case of Apertif; see van Cappellen et al. 2022 for details). As a consequence, they have different shapes, which may change in time as well as per antenna. Recovering the CB shapes in a flexible, fast, and reliable way is essential for PAF observations (see Dénes et al. 2022, for a discussion) to deliver reliable scientific products.

In this work, we propose and apply a new method to recover a compound beam shape, based on the NVSS Survey and Gaussian process regression. We apply the method to the ADR1 continuum images, mosaic them, and extract the source catalog. We also illustrate the potential for scientific application. The paper is organized as follows. In Sect. 2, we describe the ADR1 data. In Sect. 3, we introduce the “primary beam” correction. The mosaicking procedure is described in Sect. 4. The source extraction, compilation of the catalog, reliability, completeness, flux scale, and astrometry are discussed in Sect. 5. The results of cross-matching of the catalog with the LOFAR/LoTSS and the NVSS catalogs are presented in Sect. 6.

2 The data

The first imaging observations aimed at scientific verification with Apertif were begun in April 2019, followed by the launch of the full survey program on July 1, 2019. Apertif has undertaken a two-tiered imaging survey: a wide field single-pointing survey and a deeper survey with multiple pointings over a smaller area (Hess et al., in prep.).

The first-year Apertif data release (ADR1) took place in November 2020 (A22). This time period includes a total of 221 observations of 160 independent fields, covering about a thousand square degrees of sky. The released data products include continuum images, polarization images and cubes, and spectral line cubes, obtained with an automatic processing pipeline separately for every compound beam.

The WSRT array has east-west orientation and normally a field is being observed for 11-12 h using the Earth rotation to fill the uv-plane. Since the dishes have an equatorial mount, the CBs do not rotate or move with respect to the sky during an observation. In a single observation, Apertif covers forty adjacent CBs with a total area of ~10 square degrees. The raw data for every CB is processed independently with the automatic pipeline, including flagging, calibration, and imaging (Adebahr et al. 2022). For ADR1, the processed frequency range covers the 1292.5–1430 MHz interval. A dedicated validation procedure has been developed and applied based on inspection of the continuum images, which were required to have only modest calibration residuals and a noise level below 60 µJy beam⁻¹. Based on these criteria, some CB images were rejected, which is mostly dictated by the limitations of the current software to handle direction-dependent effects. More details about the first data release of the Apertif imaging surveys are given in A22. Overall, a total of 3374 continuum Stokes I images of individual CBs have been released. The sky coverage of the released continuum images is shown in Fig. 1.

Every image has a size of 3073 × 3073 pixels, with a pixel size of 4″ (the images are 200 × 200′, much larger than a primary beam response for the purpose of encompassing the side lobes for the cleaning procedure). In the Apertif field of view (FoV) layout, the centers of the CBs are equally separated from each other by 28′, resulting in every CB having up to six equidistant neighbors. The full width at half maximum (FWHM) of a CB radial attenuation profile is about 36′, however, its shape may be significantly non-symmetric and changing in time (see below). Many of the CB images are obtained for the same sky position in the framework of the medium-deep survey. The ADR1 data products are available through a virtual observatory interface².

As mentioned above, for most scientific purposes, the radial response of each CB needs to be characterized and corrected for. The primary beam characterization for the pre-Apertif WSRT was very well known (Popping & Braun 2008). However, the characteristics of the electronically formed beams from the Vivaldi antennas carry a number of complications (Dénes et al. 2022). The corresponding CBs shapes differ from those of a single-receiver antenna and from each other, being more distorted towards the edges of the field of view. Moreover, any malfunctioning of any of the PAF elements or a leakage between them can cause a significant deformation of the CB shape on a given dish. As a result, the primary beam correction used for the old WSRT is not applicable to the 40 CBs of Apertif. Below we introduce and implement an efficient approach for this correction. The new method has been applied to the ADR1 continuum images and the corresponding CB correction models have been published along with the processed data.

Fig. 1 Sky coverage of the released continuum images.

3 Compound beam shapes

Traditionally, holography is used to measure the primary beam response of radio telescopes; however Apertif has not supported a holographic mode. An alternative approach is to perform drift scans of a bright calibrator, typically Cyg A or Cas A. However, this approach has several limitations. First, it can be affected by a complex source structure, background sources, or radio frequency interference (RFI) that are present during scanning. Second, it is extremely time consuming. To perform one such calibration procedure for all 40 CBs requires more than 14h, which is longer than a full synthesis observation of a survey pointing. Finally, these calibrators are extended sources resolved on WSRT baselines and their substructure complicates the analysis even for the auto-correlation amplitudes of the visibilities, while compact calibrators such as 3C286 are not strong enough for drift-scanning. For this reason, such drift scans are performed only once per month or even less frequently. Moreover, due to observing time limitations, the resulting models of a CB are often cut off well before the first null (Dénes et al. 2022). While this approach has been used in the Apertif project, providing a good reference, it is important to have an independent, faster, and more flexible method to recover a CB shape for any given observation date.

The approach we propose and apply here is based on the public NVSS catalog data combined with a Gaussian process (GP) regression. The method takes advantage of the fact that the NVSS survey is obtained at a very similar frequency as the Apertif survey, allowing for a determination of the CB attenuation factor at any position in the FoV, for any given date of observation.

A GP is a stochastic process with normally distributed values characterized by their mean and covariance function or the kernel. The kernel depends on hyperparameters and determines the characteristics of the GP. Data can be approximated using GP regression, providing a probabilistic prediction of the GP value at a given coordinate (Rasmussen & Williams 2006).

Although the NVSS was performed at a similar frequency, it has a lower angular resolution compared to Apertif. Therefore, the first step is to convolve the Apertif CB images with a circular PSF of 45″ to match the NVSS resolution. In the next step, a source finder is used (as described in Sect. 5) and the sources are cross-matched with the NVSS catalog. Here, we used the same approach for source extraction and cross-identification, as described below in Sects. 5 and 5.6. We used all the cross-matched sources to obtain the spatial distribution of the relation between NVSS and Apertif total flux, є = S_apertif/S_nvss, over a given CB image. This relation, on average, represents the corresponding CB shape.

The next step is to involve a GP and train it on the existing data. We selected a squared-exponential radial basis function for the kernel as the simplest way to describe smooth data variations. We constructed the kernel as the sum of two squared-exponential kernels (SE1 + SE2) and a white noise one: SE1 to represent the general CB shape and SE2 to describe smaller-scale deformations related to, for instance, individual PAF elements malfunctioning. For the training, we set a wide range of priors on the hyperparameters, and used all the data for a given CB consisting of a few thousand measurements (ADR1 implies about 40 images per CB). The hyperparameters were found to be plausible (the SE1 scale describes a CB width well and the SE2 amplitude is less than 10% of the SE1 one) and stable with respect to varying sample size. After validation, we used the trained GP with the SE1 kernel to predict a CB shape. This allowed us to get the “average” CB shapes over the whole year period of the observations considered. The corresponding FITS images have been published together with the ADR1 data and can be used to correct the continuum images. We note that the peak values of the CB radial profiles differ from 1 indicating potential differences in flux scale between Apertif and NVSS surveys.

In Fig. 2, the implementation of Gaussian process regression for compound beam 01 is shown. This beam is located at one of the corners of the Apertif field of view and, therefore, is expected to have a distorted shape (e.g., Fig. 25 in van Cappellen et al. 2022). The top-left panel shows the scatter plot of the relation є for all the sources cross-matched in the CB 01 – 7154 measurements taking values from 0.0004 to 2.6 (the sizes of the markers in the plot are proportional to the value). These data have become a training sample for GP regression, whereby the result is plotted in the top-right panel. The asymmetry of the CB shape can be clearly seen from the figure. The bottom panels show slices of the GP surface along RA and Dec coordinates, where the displacement of the peak with respect to the pointing center can be seen. Both the shape distortion and the shift of the maximum must be accounted for during a primary beam correction procedure. Moreover, it was found that the CB shapes change over time, which might be related to the re-adjustment of the signals from focal Vivaldi antennas (beam weights tuning procedure) as well as to other system drifts. This effect will be investigated for future data releases; whereas in this work, we use the averaged profiles of the CBs based on the full year of observational data. Since there are many measurements of є for a given CB, this makes the GP training process very slow. To speed it up, we used the k-means clustering of these measurements into 1000 clusters. We additionally checked that GP regression obtained on such a clustered sample does not differ significantly from the one obtained on a full sample.

The proposed approach has multiple advantages. First, this method provides a ready-to-use correction image for a given CB image in the sense that it does not require knowledge of a frequency or antenna dependence of a CB shape. Next, the approach is free from approximation errors of analytical models. Indeed, the evolution of a CB slope or a shift of the peak of a CB with respect to image center may be hard to model analytically. Third, since a CB shape changes over time, we can select the data within a limited time period and, therefore, we can obtain a CB correction for a given observation period. Finally, all possible sources of bias contributing to an individual CB image are merged into just one – that of the NVSS catalog itself, which is assumed to be negligible throughout this work.

For the ADR1, there are two versions of these beams released: the “original” primary beam images, which provide a correction that matches the NVSS flux scale by design, and the “normalized” primary beam images that are normalized to a peak response of one. The data release documentation contains information on how to best apply these primary beam images to the Apertif data, including details on how to scale the primary beam images with frequency. The accuracy of the flux scale with the GPR beams is discussed in Sect. 5.6. The source code used to produce the CB models is available on GitHub³.

Fig. 2 Implementation of Gaussian process regression for compound beam 01. Top left: relation of Apertif/NVSS integrated flux density for all the sources observed in the compound beam 01 (the size of the markers is proportional to the value). Top right: GP regression built on these data. Bottom row: slices of the predicted CB shape (±1σ) along the RA and Dec coordinates.

4 Mosaicking of the images

Here, we describe the procedure for the creation of the corrected images taking advantage of mosaicking with neighboring beams in order to achieve nearly uniform noise across the field, which is important for the source extraction and the catalog.

As mentioned above and described by van Cappellen et al. (2022) and Hess et al. (in prep.), the adopted layout of CBs locations leads to an overlap of the neighboring CBs from Apertif. Thus, the sensitivity for each field can be increased by making a mosaic of these images. For every CB image, we take those closest surrounding images with their centers separated by less than 0.7 degrees from the current CB center. Since it is only the images released in ADR1 that are used in this work, some CB images have fewer than six surrounding ones. Moreover, some appear to be complete “orphan,” that is, there are 38 CB images without any neighbors. This circumstance leads to a patchy sky coverage and an increase of noise at the edges of such an orphan CB image. The number of good quality CB images will significantly increase in future Apertif data releases, in particular, due to a new pipeline capable of performing direction-dependent calibration. After selecting a central CB image and its neighbor beams, the smallest Gaussian restoring beam that encloses all individual restoring beams of the selected image is calculated. This one is typically a few arcseconds larger than the individual ones. All the images are then convolved to this common restoring beam. Next, each image is corrected for the corresponding CB shape clipped at the CB attenuation level of 0.1, which was chosen to avoid very high noise levels at the mosaic edges. After this, the images are re-projected onto a common projection center and stacked together into a mosaic image. Every mosaic image therefore covers up to 3.7 square degrees and represents a stack of individual CB images. The resulting mosaic images are used for source finding. We also generated a multi-order coverage map (MOC) of the survey, indicating an overall sky coverage for the mosaics totaling 970 square degrees. The source code for mosaicking is available on GitHub⁴.

An example of the mosaic image made of 57 individual CB images of a region of the medium-deep survey (all released CB images of several observations for these coordinates) is shown in Fig. 3. In the right panel of the figure, the residual map after source extraction is shown (see the next section). The dashed circles and the numbers indicate the positions and the number of CB images used for making the mosaic. The image has one of the lowest local noise levels in its central region reaching 17 µJy beam⁻¹. Imaging artifacts around some brighter sources clearly visible in the residual image on the right are caused by direction-dependent effects, which will be solved for in upcoming data releases with an improved pipeline.

Fig. 3 Mosaic example of 57 individual CB images with the center at RA = 29.26, Dec = 34.55 degrees (left panel) and the residual noise map (right panel). Seven dashed circles show the positions of compound beams, and the numbers inside indicate the quantity of the corresponding CB images used for mosaicking. The color bar scale is given in mJy beam−1.

5 Source extraction and catalog

In each mosaic, the source extraction is done using the Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015), which is known to be one of the best performing among continuum source-finding software (Hopkins et al. 2015). It produces a source catalog from an image by grouping together emission components modeled by Gaussians. Such a group is represented by an “island” namely, a region of pixels having flux above a certain threshold. With this approach, PyBDSF efficiently groups emission blobs into sources. We set the default peak-to-noise source detection threshold to 5 and the island boundary threshold to 3. The background noise level was determined inside a sliding box of 128 pixels across the image with a step of 16 pixels (rms_box=(128,16)). The box size was adjusted dynamically, decreasing to 16 pixels in the regions of high background noise (adaptive_rms_box=True, rms_box_bright=(16,4)). The imaging artifacts near bright sources occur mostly due to direction-dependent effects and appear as concentric circles around a source, affecting the estimated noise level. The size of these rings ranges from a few to several tens of the size of the restoring beam. These characteristic scales were used for the size of the sliding box to properly detect the background noise variations. We tested a few other combinations of the rms_box and rms_box_bright parameters and we used one resulting in a lower number of false detections. We note, that a background rms estimate within a small box might experience additional variations due to confusion, however, this effect is small for Apertif which has a confusion limit of ~5 µJy beam⁻¹. A more justified and automated approach for the parameters setup should be developed in future.

Thus, for every mosaic image, we obtain a source list which includes the coordinates, size, signal-to-noise ratio (S/N), and other fitting parameters. We used the source position and size from the PyBDSF output to identify duplicate sources in different mosaic images. Finally, the source lists for each mosaic image were combined together and duplicate sources having lower S/N were removed.

5.1 Extra filtering

In Apertif images, a “ghost” source is present at the image center. This is the result of some wide band radio frequency interference (RFI) within the backend. This RFI signal has a random phase that averages 0 and produces a fake source during imaging. Upon noting this, we removed those ghost sources from the catalog by matching their coordinates with the CB centers.

We also note that in spite of the strong validation criteria applied for ADR1 and the carefully adjusted PyBDSF parameters, some artifacts were still detected as extended sources. They tend to have an order of magnitude higher coordinates uncertainties and can thus be easily identified. We inspected visually the sources with either RA or Dec errors higher than 10″ and dropped the detected artifacts.

5.2 Catalog format

The resulting catalog contains 249672 sources. Almost 90% of the sources are modeled with a single Gaussian (‘S’-flag), and around 10% are classified as extended (‘M’). An example of the catalog structure is shown in Table 1. The columns designations are: (1): Apertif source name; (2,4): RA and Dec; (3,5): RA and Dec errors (see Appendix A); (6,8): total and peak flux density; (7,9): integrated and peak flux density uncertainties, calculated as described in Appendix A; (10,12,14): Deconvolved major and minor source size and position angle. A source size value 0.0 means that the PSF cannot be deconvolved from the fitted source along the given axis, and the corresponding uncertainty, represents an upper size limit estimate. The position angle is given in degrees from −90 to 90 measured relative to the north celestial pole, turning positive into the direction of the right ascension. When both major and minor sizes are 0.0, the position angle is omitted; (11,13,15): uncertainties of the major and minor source size and position angle calculated as described in Appendix A; (16): local background noise rms; (17): Source type as classified by PyBDSF, with ‘S’ as an isolated source fitted with a single Gaussian; ‘C’ as sources that were fit by a single Gaussian but are within an island of emission that also contains other sources, and ‘M’ as sources fitted with multiple Gaussians.

5.3 Correction for spectral index

The central frequency of the NVSS survey observing band differs from the Apertif one by 45 MHz. Within the aforementioned approach to the primary beam correction, the source fluxes would undergo a systematic offset depending on their spectral index⁵ (see also Sect. 6.2). A flux correction for a source with a spectral index of −0.7 would be ξ(−0.7) = 1.023, or 2.3%. We note that the dependence of the correction factor on the spectral index is weak: ξ(−0.5) = 1.016, = 1.033. In order to eliminate this bias we apply the correction factor of 1.023 to all sources in the catalog.

5.4 Noise properties

Analysis of the noise properties is important for the discussion that will follow and to quantify the completeness of the catalog. The background noise varies strongly from one mosaic image to another as well as within a single mosaic image. This can be seen in the residual image in the right panel of Fig. 3. The noise in this image reaches 17 µJy beam⁻¹ in the central region and goes to above 100 µJy beam⁻¹ at the edges. For the catalog, we effectively selected sources from a central part of every individual mosaic because the sources at the edges are typically duplicated or substituted by ones from central region of another adjacent mosaic having higher S/N values. This means that a more informative parameter is the noise level around the cataloged sources, not in an entire mosaic map. The relevant noise measure is the median background RMS of the sources provided by PyBDSF (column 16 in the catalog), which is 44 µJy beam⁻¹. We also note that because of the ADR1 validation criteria, many CB images have fewer than six neighbors, which implies a higher noise level at part of the edges.

Fig. 4 Relative number of false detections as a function of flux density.

5.5 Reliability and completeness

To estimate the reliability of the catalog based on the number of false detections, we ran the same source-finding procedure on the inverted (multiplied by −1) mosaic images. Assuming a symmetric noise distribution, this offers a rough estimate of the number of fake sources in the original images. We compiled the “inverted” catalog in the same way as the original one. It contains 6108 sources or 2.4% of the number of sources in the main catalog. In Fig. 4, the false detection rate (FDR) is plotted against flux density, where the FDR is the relative number of sources from inverted catalog to the number of the sources from the original catalog taken in various flux density intervals. The peak of the FDR around a few mJy is likely related to the typical flux of imaging artifacts around bright sources, while artifacts around fainter sources are below the catalog detection limit. We also found that 64% of the false detections reside within 2′ of sources from the original catalog that are brighter than 50 mJy. Overall, among the sources weaker than 30 mJy, a few percent may be false detections, while brighter sources are shown to be more reliable detections.

We estimated the completeness of the catalog by comparing the differential source counts corrected for the false positive detection rate to the one by Matthews et al. (2021). The comparison is shown in Fig. 5. The ratio (dN/dS)_Apertif/(dN/dS)_Matthews is shown in the inset plot. Slightly lower counts of bright sources in the Apertif catalog are dictated by the ADR1 validation criteria, namely, by the lack of the released images with stronger sources due to them having imaging artifacts (see Sect. 2). Another effect might be caused by an underestimation of the flux for bright sources due to some of them having been partially resolved with Apertif. On the weaker end, we find that at the 1 mJy level, the catalog completeness is about 96% and it drops to ~75% at the 400 µJy level. We emphasize that this level was calculated for the full catalog. In reality, the number of weak sources and, hence, the “local” completeness in a given sky direction depends on the local noise level in an image. To illustrate this, in Fig. 5, we plot the source counts calculated for the low-noise mosaic image (shown in Fig. 3), which covers approximately 3 square degrees. The corresponding “local” catalog is then complete down to the 0.1 mJy level.

5.6 Accuracy of the flux scale

In order to assess the flux scale accuracy, we cross-matched the Apertif catalog with the NVSS one using a 20″ matching radius. This results in 44524 common sources, for which the integrated flux density differs by less than 1% (on average), with a median ratio of S_int,aptf/S_nvss = 1.027. We note that the applied correction factor for the spectral index is 1.023 (as given in Sect. 5.3). Obviously, such a good agreement is a consequence of the applied method for the primary beam correction (Sect. 3). This comparison serves as a firm additional check for the flux scale of the Apertif catalog as a whole, especially since it has a different angular resolution than the NVSS catalog. The standard deviation of the above relation is 0.345, indicating a significant scatter. This spread must be related to the variability of the sources, the varied angular resolution and sensitivity characterizing the catalogs, and errors in the flux measurements.

In order to estimate the precision of the flux measurements, we compared these measurements obtained from the CB images of the medium-deep fields, which have at least five pointings. The variability of flux measurements within individual Apertif images is shown in the top panel of Fig. 6, where the relative difference of the peak flux and the mean peak flux is plotted against S/N of the sources. We split the full S/N range into 15 bins, then for every bin that we calculated, the average relative peak flux error and 16/84 percentiles of the distribution are shown with the dashed and solid curves, respectively. We conclude that the flux measurements precision is better than 10%, within the expected errors coming from data cross-calibration.

Finally, we checked the intrinsic flux scale consistency by comparing the fluxes of the sources from different overlapping mosaics. When a source was detected in multiple images, only the corresponding record with highest S/N was stored in the catalog. Thus, comparing the fluxes of those sources provide a direct probe of the applied primary beam correction and overall correctness of the mosaicking procedure (Sects. 3 and 4). For these sources, we again considered the relative difference between the peak flux and its mean value. The distribution is plotted in the bottom panel of Fig. 6 against the S/N of the sources. The scatter indicates that the error during primary beam correction and mosaicking do not exceed 5%. We take this into account for the flux error values reported in Table 1 (described in Appendix A).

Fig. 5 Differential source counts. The source counts estimates for the full catalog area (circles) and for the low-noise mosaic image (Fig. 3, squares) plotted along with the measurements of Matthews et al. (2021; dot-dashed curve). The ratio (dN/dS)Apertif/(dN/dS)Matthews is shown in the inset plot.

5.7 Astrometry

To estimate the astrometric accuracy of the sky coordinates, it is desirable to have a reference catalog with a comparable or higher angular resolution. Fortunately, there is one in this case, namely, the first Apertif data release has a common sky area of 700 square degrees with the second LOFAR 150 MHz data release (Shimwell et al. 2022). We cross-matched the Apertif catalog with the LoTSS DR2 catalog. We note that the LoTSS images have about twice better angular resolution (6″) than those from Apertif, and some sources that are resolved by LOFAR might remain unresolved by Apertif. We used an intermediate matching distance of 10″ resulting in 152824 common sources and considered the coordinate offsets between Apertif and LoTSS DR2 sources. In Fig. 7, the offset distributions for Right Ascension and Declination are plotted against S/N of the sources. The median values are 0.0″ and 0.0″, respectively, indicating the absence of any systematic shifts between the catalogs. The 16/84 percentiles are shown with the solid curves justifying the chosen matching radius. The coordinate uncertainties for the catalog were calculated as described in Appendix A.

Fig. 6 Flux measurements precision. Relative difference of the peak flux measurements over time (top) and between mosaic images (bottom) against the S/N of the sources. The dashed curves indicate average relative peak flux error. The solid curves display the 16 and 84 percentiles.

6 Scientific applications

In this section, we discuss some scientific applications related to cross-matching the Apertif catalog with the NVSS and LOFAR/LoTSS catalogs. The cross-matching of catalogs obtained at similar frequency but at distinct epochs allows for the identification of transient sources that changed their flux density significantly between the epochs. Finding common sources between two contemporaneous catalogs at different frequencies provides information about their spectral behavior.

Fig. 7 Astrometric accuracy. Coordinate offsets of the compact sources based on Apertif – LOFAR cross-match. Top: right ascension (RA). Bottom: declination (Dec). The horizontal lines indicate median value, the solid curves show the 16 and 84 percentiles, and the dashed curves show the corresponding hyperbola fit (see details in Appendix A).

6.1 Long-term transients

As an example of transient detections, we selected the sources that are present in the NVSS but have disappeared in the Apertif images. First, we exclude the Apertif sources separated by less than 50″ from a mosaic edge, where the image quality is lower due to the noise. Then we cross-matched the remaining Apertif sources with the NVSS and obtained a list of more than 4000 NVSS sources missing in our catalog. The great majority of these sources have a total flux density below 3 mJy and represent NVSS false detections. We filtered these sources by requiring the NVSS S/N to be higher than 10. Another large group consists of multi-component emission islands which were grouped by PyBDSF into sources whose location mismatch with the corresponding NVSS source. To filter these sources, we inspected the Apertif images and required the signal-to-noise ratio at the position of each NVSS source to be less than 6. This ensures that there is no significant emission in Apertif images at the position of an NVSS source. Finally, we obtained a list of 23 NVSS sources that have a NVSS S/N higher than 10 and were not detected in the Apertif data. The list is presented in Table 2.

Apparently, the source with the highest total flux density (NVSS:221119+385955) is a diffuse lobe of giant radio galaxy that was resolved by Apertif. The next brightest source, NVSS:224624+440950, is seen in the NVSS image (with a flux density of 19.8 mJy) and is clearly absent in the Apertif one (Fig. 8). The figure also nicely illustrates the differences in resolution and sensitivity between the Apertif and NVSS surveys. The local noise value in the Apertif mosaic at this position is the 120 µJy beam⁻¹, implying that the source, if it is real and point-like, had dropped its luminosity by a factor of ~30. Such a change could, for example, have taken place if the source is a blazar whose jet shifted its angle to the line of sight, causing a significant decrease in Doppler boosting. A factor of 30, however, is unusual even for extreme blazars (e.g., Kutkin et al. 2018). Another explanation could be attributed to a supernova explosion. A detailed study of these objects (e.g., measuring redshifts, searching for counterparts in other bands, extra observations etc.) might be an interesting project for future research. We note that most sources in Table 2 are relatively faint, therefore, some might be just variable sources picked up by NVSS when they were at maximum and missed by Apertif when they were at minimum. These transients, of course, deserve a separate study.

The above-mentioned procedure illustrates a use case example for the new catalog. We note that a much larger sample of transient sources can be obtained by considering the significance of the flux density change – and not just sources that have disappeared. Such a detailed analysis, however, is beyond the scope of this paper. For the same reason, we did not consider sources that are present in the Apertif survey and absent in the NVSS.

Fig. 8 Apertif (left) and NVSS (right) images. The disappeared source is marked with the dashed circle and has total flux density of 19.8 ± 0.7 mJy.

6.2 LOFAR counterparts and spectral index

Since the observing frequency of the LoTSS survey is ten times lower than that of Apertif, the combination of the two surveys provides unique information about the spectra of the sources; in turn, this offers additional information about the emission mechanisms allowing us to distinguish among the emitting regions of the sources (e.g., flat spectrum radio cores vs. optically thin jet components with steep spectra). In addition, spectral index measurements can be used to estimate source ages (see Kardashev 1962).

Williams et al. (2016) and Mahony et al. (2016) used LOFAR images of the Böotes and Lockman Hole fields, respectively, to derive the spectral index distribution between 150 and 1400 MHz. Their analysis results in very accurate median values of the spectral index of −0.79 ± 0.01 and −0.78 ± 0.02, correspondingly. These studies were limited to a relatively small area and the effect of differences in spatial resolution between the 150 and the 1400 MHz images could be taken into account. The authors used both resolved and unresolved sources for the spectral index estimation.

Because of the different resolution between the Apertif and LOFAR catalogs, we limited our analysis to those sources unresolved in the latter, noting that these sources must be unresolved in the former as well. We selected these sources following the criteria described by Shimwell et al. (2022, in particular, Sect. 3.1; Eq. (2)), when the compactness depends on the ratio between integrated and peak flux density of a source as well as the S/N level of its detection. Consistently with this finding, we also observed that about 90% of the common sources are unresolved, which is close to the fraction of 95% obtained for all LOFAR sources. We estimate the spectral index as α = ln (S_total,A/S_total,L)/ln (1355/150), where S_total,A, S_total,L are the Apertif and LOFAR total flux density of a source, and 1355 MHz and 150 MHz are the corresponding survey frequencies.

The spectral index distribution as a function of the LOFAR total flux is shown in Fig. 9. The trend of flattening of the spectral index when going to lower fluxes occurs mostly due to the combination of increasing scatter for low flux sources and the Apertif sensitivity limit and completeness; this is illustrated by the overplotted curves. The envelope shown with the dashed line corresponds to the error due to the flux density uncertainties and is derived as follows: $${\sigma }_{\alpha ,S}=\frac{1}{\ln \left(1355\text{\hspace{0.17em}MHz/150MHz}\right)}\sqrt{{\left(\frac{{\sigma }_{S_{\text{total}}}}{S_{\text{total}}}\right)}_{\text{Apertif}}^2+{\left(\frac{{\sigma }_{S_{\text{total}}}}{S_{\text{total}}}\right)}_{\text{Lofar}}^2}$$

The left solid line shows the sensitivity limit calculated for the lowest Apertif peak flux and the right one indicates the limit for the Apertif catalog completeness level of 1 mJy. Therefore, an estimate of a median spectral index is unbiased for the sources with LOFAR flux density above 40 mJy (vertical dotted line). This estimate is α = −0.78, with the scatter of σ_α = 0.26 (horizontal dash-dotted lines). As seen in Fig. 9 most of the scatter comes from real differences in spectral index of the sources, while the uncertainty due to flux density measurements is relatively small.

In addition to the cross-matching with LoTSS DR2, we cross-matched the Apertif catalog with the value-added LoTSS DR1 catalog. The latter has, among other parameters, optical counterparts and redshifts estimates for the sources (Duncan et al. 2019), which might be useful for studying either individual objects or statistical populations. The common area of the two catalogs is about 100 square degrees resulting in 22825 matched sources.

We published the list of common Apertif-LOFAR sources in Table 3. The table contains names (Cols. 1 and 2), peak, and integrated flux density (Cols. 3–6) from the Apertif and LotSS/DR2 catalogs, spectral index estimate calculated for the peak fluxes (Col. 7), angular separation between the sources (Col. 8), and redshift estimate from LoTSS DR1 value-added catalog (Col. 9; zbest as described by Duncan et al. 2019). With this list, we significantly increase the population of radio sources with spectral indexes measurements. The results can be used for studying, for example, individual objects, while the approach itself might be relevant for deeper data releases in the future.

Fig. 9 Spectral index distribution against the LOFAR integrated flux density for compact sources. The median value α = −0.78 (horizontal line) is calculated for the sources with integrated flux density above 40 mJy at 150 MHz. The surrounding dashed curves show ±σα,S error due to flux density uncertainties. The inclined solid lines show spectral index lower limit calculated for the weakest Apertif source and for 1 mJy level corresponding to the completeness of 95%. An overall scatter, ±σα, of the spectral index is shown with dash-dotted lines.

7 Summary

To deliver a ready-to-use source catalog for the scientific community, we processed and analyzed 3072 continuum images of the first data release of the Apertif survey. Over the course of this work, we tackled and solved some challenging technical problems.

To obtain a primary beam model, we proposed and applied a new method based on Gaussian process regression. This machine learning approach is a very natural choice because the primary beam shape and its peak are known to be smooth (Gaussian-like) and have a value of unity, which makes the choice of a kernel function and priors on hyperparameters straightforward. This method can be easily applied to other facilities and has several advantages compared to analytical modeling, or to the drift scanning procedure. The most important is the time economy (every drift scanning takes longer than 12 h). The method provides a ready-to-apply correction model for every compound beam, effectively removing a frequency- and an antenna-dependence. Instead of using many parameters for an analytical approximation there are just a few hyperparameters entering the kernel function, which substantially simplifies further analyses. Rather importantly, it becomes easy to obtain a primary beam model for any given date by taking into account only the observations around this date. This makes tracking primary beam variations easier. Finally, all possible sources of attenuation bias are effectively removed by adjusting the flux scale to the one of the NVSS catalog. The primary beam models obtained with this method have been published, along with the data products for the first Apertif data release.

In order to increase the sensitivity of the images, we created linear mosaics of the individual images that take into account the different size and orientation of the restoring beam in each of them. The mosaic of all images of the first Apertif data release covers 970 square degrees of sky. We present a new continuum source catalog for these data, which contains 249672 radio sources, many of which have been detected for the first time at L-band (1400 MHz). The full catalog has a completeness of about 95% at the 1 mJy level. The number of false positive detections does not exceed a few percent. The table contains the coordinates, integrated and peak flux density, and angular sizes of the sources.

We cross-matched the new Apertif catalog with the NVSS and the LOFAR/DR2 catalogs. On the one hand, this procedure provides an extra check for the obtained parameters of the sources, such as their fluxes and coordinate offsets. On the other hand, it offers sizable potential for a scientific research. The first sample provides an opportunity to detect long term transient sources, which have significantly changed their flux density over the last 25 yr. The second one includes flux measurements at 150 and 1400 MHz providing information about spectral properties of more than a hundred thousand sources. The redshift estimates adopted from the value-added LoTSS catalog might be useful for studying individual objects as well as volume-limited samples.

The released images were selected based on strong validation criteria providing a non-uniform sky coverage. As our next steps, we are looking forward to applying the described methods and techniques to the future data releases, where the amount and quality of the images will be significantly improved.

Acknowledgements

We thank the anonymous referee whose comments helped to improve the paper. This work makes use of data from the Apertif system installed at the Westerbork Synthesis Radio Telescope owned by ASTRON. ASTRON, the Netherlands Institute for Radio Astronomy, is an institute of the Dutch Science Organisation (De Nederlandse Organisatie voor Wetenschappelijk Onderzoek, NWO). Apertif was partly financed by the NWO Groot projects Apertif (175.010.2005.015) and Apropos (175.010.2009.012). This research made use of Python programming language with its standard and external libraries/packages including numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2011), matplotlib (Hunter 2007), pandas (McKinney 2010) etc. This research made use of Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, Astropy Collaboration 2018). The radio_beam and reproject python packages are used for manipulations with restoring beam and reprojecting/mosaicking of the images. This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France (Bonnarel et al. 2000; Boch & Fernique 2014). E.A.K.A. is supported by the WISE research programme, which is financed by the Netherlands Organization for Scientific Research (NWO) B.A. acknowledges funding from the German Science Foundation DFG, within the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source to Signal”. T.A.O. acknowledges funding from NWO via grant TOP1EW.14.105. K.M.H. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science, Innovation and Universities through the “Center of Excellence Severo Ochoa” awarded to the Instituto de Astrofísica de Andalucía (SEV-2017-0709), from the coordination of the participation in SKA-SPAIN, funded by the Ministry of Science and Innovation (MCIN). K.M.H. and J.M.vdH. acknowledge funding from the Europeaní Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement No. 291531 (‘HIStoryNU’). L.C. and L.C.O. acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement No. 617199. J.vL. acknowledges funding from Vici research programme “ARGO” with project number 639.043.815, financed by the Dutch Research Council (NWO). D.V. acknowledges support from the Netherlands eScience Center (NLeSC) under grant ASDI.15.406.

Appendix A Error estimations

Appendix A.1 Flux scale

As shown in Figure 6, even the bright sources have about 5% uncertainty in their flux density measurements. This uncertainty likely comes from the cross-calibration procedure (A22). In order to take this into account, we included the 5% uncertainty into the final error of the flux density measurements in Table 1 by summing (in quadrature) 0.05S and the fitting error obtained from PyBDSF fit σ_*,fit: $${\sigma }_{\text{S,tot}}=\sqrt{{\left(0.05S_{\text{tot}}\right)}^2+{\sigma }_{\text{tot,fit}}^2},{\sigma }_{\text{S,peak}}=\sqrt{{\left(0.05S_{\text{peak}}\right)}^2+{\sigma }_{\text{peak,fit}}^2}.$$

Appendix A.2 Coordinates

The coordinate uncertainties of the LoTSS DR2 catalog are much smaller (~ 0.2″) than those of Apertif, which has lower resolution. Therefore, the former catalog can be used as a reference. Additionally to the fitting errors, σ_RA,fit and σ_Dec,fit, provided by PyBDSF there are astrometric errors, σ_RA,ast, and σ_Dec,ast, coming from any systematic effects in an individual mosaic image. As shown in Figure 7, the coordinates offsets are scattered with a dispersion ranging from ~2″ to 0.5″ depending on the source’s S/N. We fit the percentile curves with hyperbola ${\sigma }_{*,\text{ast}}=1\text{/}{\left(a_{*}·SNR+b_{*}\right)}^{n_{*}}$, resulting in a_ra = 0.031, b_ra = −0.132, n_rA = 0.220 and a_Dec = 0.047, b_Dec = −0.207, and n_Dec = 0.281. The fit was used to predict astrometric errors of the coordinates depending on the source’s S/N. The final coordinate errors (reported in Table 1) were then calculated as a quadratic sum of the astrometric error and the error from PyBDSF fitting: $${\sigma }_{\text{RA}}=\sqrt{{\sigma }_{\text{RA,ast}}^2+{\sigma }_{\text{RA,fit}}^2},{\sigma }_{\text{Dec}}=\sqrt{{\sigma }_{\text{Dec,ast}}^2+{\sigma }_{\text{Dec,fit}}^2},$$

with the average coordinate errors being $\overline{{\sigma }_{\text{RA}}}=1.7^{″}$ and $\overline{{\sigma }_{\text{Dec}}}=2.1^{″}$.

Appendix A.3 Deconvolved size

Along with the FWHM of a Gaussian fitted to a source, PyBDSF also provides the deconvolved source size and position angle (PA), and the corresponding errors. The latter values, however, are not calculated correctly and remain the same as the errors of the original fit (the PyBDSF version used in this work is 1.9.2). We estimated the errors of a source size and PA after deconvolution, σ_maj,min and σ_PA, by taking the maximum of the absolute values of the differences of the best fit deconvolved value of the given parameter and the deconvolved size of the eight possible combinations of 1) FWHM major axis +/− major axis error; 2) FWHM minor axis +/− minor axis error; and 3) position angle +/− position angle error. If the source cannot be deconvolved from the beam we report the upper limits which come from the maximum major and minor axes of the deconvolved Gaussians taken from trying all eight of the aforementioned combinations. This approach was adopted from the CASA imfit⁶ task.

In the case that none of these combinations produce a deconvolved size, we estimate the upper limit using the major and minor axis FWHM of the restoring PSF and source S/N. A point source has a fitted size of the restoring PSF, and the corresponding error of the fitted size is ${\sigma }_{\text{maj},min}^{\text{fit}}\approx {\text{PSF}}_{\text{maj,min}}\text{/}\left(\text{S/N}\right)$ (e.g., Fomalont 1999). Then the error after deconvolution can be estimated as the smallest deconvolved size detectable at 1σ level: $${\sigma }_{\text{maj,min}}={\left({\left({\text{PSF}}_{\text{maj,min}}+{\text{PSF}}_{\text{maj,min}}\text{/}\left(\text{S/N}\right)\right)}^2+{\text{PSF}}_{\text{maj},\min }^2\right)}^{1/2}=\frac{\text{PSF}}{\left(\text{S/N}\right)}{\left(1+2\text{S/N}\right)}^{1/2}.$$

If both the major and minor deconvolved size are zero, the reported position angle is zero and its error is omitted. The deconvolved source size, PA, and the corresponding errors are listed in columns 10–15 of Table 1.

References

All Tables

All Figures

	Fig. 1 Sky coverage of the released continuum images.
In the text

	Fig. 2 Implementation of Gaussian process regression for compound beam 01. Top left: relation of Apertif/NVSS integrated flux density for all the sources observed in the compound beam 01 (the size of the markers is proportional to the value). Top right: GP regression built on these data. Bottom row: slices of the predicted CB shape (±1σ) along the RA and Dec coordinates.
In the text

	Fig. 3 Mosaic example of 57 individual CB images with the center at RA = 29.26, Dec = 34.55 degrees (left panel) and the residual noise map (right panel). Seven dashed circles show the positions of compound beams, and the numbers inside indicate the quantity of the corresponding CB images used for mosaicking. The color bar scale is given in mJy beam−1.
In the text

	Fig. 4 Relative number of false detections as a function of flux density.
In the text

	Fig. 5 Differential source counts. The source counts estimates for the full catalog area (circles) and for the low-noise mosaic image (Fig. 3, squares) plotted along with the measurements of Matthews et al. (2021; dot-dashed curve). The ratio (dN/dS)Apertif/(dN/dS)Matthews is shown in the inset plot.
In the text

	Fig. 6 Flux measurements precision. Relative difference of the peak flux measurements over time (top) and between mosaic images (bottom) against the S/N of the sources. The dashed curves indicate average relative peak flux error. The solid curves display the 16 and 84 percentiles.
In the text

	Fig. 7 Astrometric accuracy. Coordinate offsets of the compact sources based on Apertif – LOFAR cross-match. Top: right ascension (RA). Bottom: declination (Dec). The horizontal lines indicate median value, the solid curves show the 16 and 84 percentiles, and the dashed curves show the corresponding hyperbola fit (see details in Appendix A).
In the text

	Fig. 8 Apertif (left) and NVSS (right) images. The disappeared source is marked with the dashed circle and has total flux density of 19.8 ± 0.7 mJy.
In the text

Fig. 9 Spectral index distribution against the LOFAR integrated flux density for compact sources. The median value α = −0.78 (horizontal line) is calculated for the sources with integrated flux density above 40 mJy at 150 MHz. The surrounding dashed curves show ±σα,S error due to flux density uncertainties. The inclined solid lines show spectral index lower limit calculated for the weakest Apertif source and for 1 mJy level corresponding to the completeness of 95%. An overall scatter, ±σα, of the spectral index is shown with dash-dotted lines.

In the text