Characterising the Apertif primary beam response

Context. Phased Array Feeds (PAFs) are multi element receivers in the focal plane of a telescope that make it possible to form simultaneously multiple beams on the sky by combining the complex gains of the individual antenna elements. Recently the Westerbork Synthesis Radio Telescope (WSRT) was upgraded with PAF receivers and carried out several observing programs including two imaging surveys and a time domain survey. The Apertif imaging surveys use a configuration, where 40 partially overlapping compound beams (CBs) are simultaneously formed on the sky and arranged in an approximately rectangular shape. Aims. This manuscript aims to characterise the response of the 40 Apertif CBs to create frequency-resolved, I, XX and YY polarization empirical beam shapes. The measured CB maps can be used for image deconvolution, primary beam correction and mosaicing of Apertif imaging data. Methods. We use drift scan measurements to measure the response of each of the 40 CBs of Apertif. We derive beam maps for all individual beams in I, XX and YY polarisation in 10 or 18 frequency bins over the same bandwidth as the Apertif imaging surveys. We sample the main lobe of the beams and the side lobes up to a radius of 0.6 degrees from the beam centres. In addition, we derive beam maps for each individual WSRT dish as well. Results. We present the frequency and time dependence of the beam shapes and sizes. We compare the compound beam shapes derived with the drift scan method to beam shapes derived with an independent method using a Gaussian Process Regression comparison between the Apertif continuum images and the NRAO VLA Sky Survey (NVSS) catalogue. We find a good agreement between the beam shapes derived with the two independent methods.


Introduction
Phased-array feeds (PAFs) are an important advancement in radio astronomy. These instruments are multi-element receivers in the focal plane of a telescope that make it possible to simultaneously form multiple beams on the sky by combining the complex gains of the individual antenna elements. This significantly helgadenes@gmail.com increases the instantaneous field of view (FoV) and the survey speed of the telescope.
Over the last few years, several telescopes have been equipped with PAF receivers, from single dishes such as the Robert C. Byrd Green Bank Telescope (GBT; e.g. Landon et al. 2010;Roshi et al. 2018;Pingel et al. 2021) and the Effelsberg telescope (e.g. Chippendale et al. 2016;Reynolds et al. 2017;Deng et al. 2017) to interferometers such as the Australian SKA Pathfinder Telescope (ASKAP; e.g. Hotan et al. 2014Hotan et al. , 2021 Article number, page 1 of 14 Article published by EDP Sciences, to be cited as https://doi.org/10.1051/0004-6361/202244045 A&A proofs: manuscript no. 44045corr McConnell et al. 2016) and the Westerbork Synthesis Radio Telescope (WSRT), equipped with a PAF system called Apertif (APERture Tile In Focus). Out of the 14 WSRT dishes, 12 of them were upgraded with the new PAF receivers. Each of the PAFs has 121 Vivaldi antenna elements that can be combined to form multiple beams simultaneously on the sky. In the standard observing mode for imaging, there are 40 partially overlapping compound beams (CBs) formed in a rectangular configuration (Figure 1). The Apertif system is tunable to observe between 1130-1750 MHz in a 300 MHz band. The observing band is subdivided into 384 subbands of 781.25 kHz each, and each subband is further subdivided into 64 channels of 12.2 kHz each. For a more complete description of the Apertif system, we refer to van Cappellen et al. (2022). Apertif was designed as a survey instrument and was operated in survey mode between 1 July 2019 and 28 February 2022. Apertif carried out a two tiered imaging survey: 1) with a shallow, wide-area and 2) a medium-deep, small-area component (see Hess et al. in prep), as well as a time domain survey (see Maan & van Leeuwen 2017;van Leeuwen et al. 2022). All of the Apertif surveys are legacy surveys providing publicly available data products. The first data release of Apertif (DR1) was in November 2020  1 .
Accurate knowledge of the primary beam of a telescope is important to achieve the correct deconvolution for primary beam correction and for mosaicking of radio images. It is furthermore essential for the accurate instantaneous localisation of transients such as fast radio bursts (Connor et al. 2020;van Leeuwen et al. 2022). An incorrect primary beam can lead to severe imaging artefacts and an incorrect flux scale. This can prove to be especially problematic for telescopes with PAF receivers, where the numerous simultaneously formed primary beams can all have significantly different shapes; this is in contrast to telescopes with traditional receivers with only one primary beam shape. An example for imaging artefacts due to an inaccurate primary beam model is concentric ring-shaped direction-dependent errors (DDEs) around bright continuum sources. Figure 3 in Adebahr et al. (2022) shows an example of such DDEs in the Apertif data. These DDEs can be mitigated with rigorous flagging of CBs that have low-quality data (e.g. if CBs on one antenna have strong distortions compared to the other antennas; for more, see Figure 2) and advanced imaging techniques that include direction-dependent calibration or custom antenna-based primary beam models. This makes it essential to map and understand the shape of the primary beam as a function of frequency and time. As part of the Apertif imaging survey's data products, we provide primary beam maps for each of the 40 digitally formed CBs.
In this work, we present the characterisation of the Apertif CBs using two independent methods: drift scan measurements and a Gaussian process (GP) regression method. In Section 2, we briefly describe how the CBs are formed. In Section 3, we describe the drift scan observations and the method used to create the beam maps. In Section 4, we present the properties of the CBs, such as beam sizes, frequency dependence, and time dependence, as well as a comparison between the drift scan beam maps and the GP beam maps. Finally, in Section 5, we summarise our findings.

Apertif compound beams
Apertif CBs are formed by applying beam weights to the data from the individual antenna elements. This can be done with different beam forming optimisation methods (e.g. Jeffs et al. 2008;Ivashina et al. 2011). Beam forming can be optimised for a high signal-to-noise ratio (S/N) or the shape of the CBs. When maximising the S/N, the shape of the CBs gets gradually distorted with distance from the optical axis of the focal plane. Apertif uses maximum S/N beam weights, which means that CBs away from the optical axis suffer from coma distortions. Similar distortions are also observed for the CBs of ASKAP (e.g. McConnell et al. 2016).
In the case of Apertif, CBs are formed so that most of the signal is contributed by only a few elements (∼ nine elements) and, hence, the sensitivity and the shape of the CBs strongly depend on the properties and health of the main contributing elements. The health of the Apertif system is assessed through monitoring all hardware and software components of the telescope (see detailed description in van Cappellen et al. (2022) i.e. Sections 9 and 11). The components of the system health that are most critical for the quality of the science observations are the state of the individual Vivaldi antenna elements and their signal path to the correlator. To monitor this aspect, the power against frequency is recorded from all antenna elements. Occasionally, the power drops from the nominal 94 dB to 70 dB for an individual element. When this happens, the element is considered malfunctioning and the signal path from it is disabled. This means that broken elements are excluded from the beam weights measurements and the scientific observations. Depending on how many and which elements fail on a certain PAF, some CB shapes can be significantly distorted. Figure 2 illustrates this scenario via a comparison of the beam weights and beam maps of a PAF with broken elements to a PAF with no broken elements. The top row of Figure 2 shows the beam weights and the drift scan map of CB18 on antenna RT8 and the bottom row shows the same data on antenna RTB, where there are several broken elements. On the healthy PAF, the main contribution for a beam comes from approximately nine elements and the resulting CB shape is regular, while the broken elements on RTB distort the CB shape. To mitigate these distortions, it is possible to repair the broken elements on the PAFs, however, due to the complicated nature of the PAF maintenance, elements are only repaired once a critical number (from five to seven elements) fail on a PAF. On the upside, because of the small number of elements contributing to Broken elements are white pixels in the image. Right: Corresponding CB map from a drift scan measurement. White contour levels are: 0.2, 0.4, 0.6, 0.8. Red contour levels are: 0.1 and 0.5. The colour scale is the same for all four images. Most of the signal for a given beam is contributed by nine antenna elements. If one of these antenna elements is broken, the CB for the given antenna will be distorted. In this case, two of the highest contributing elements on antenna RTB were broken. each individual CB, the correlated noise between CBs is relatively low (i.e. on the level of a few percent).
The CBs are formed for each 781.25 kHz wide subband and the beam characteristics can change across the subbands. One problem in particular is strong intermittent radio frequency interference (RFI), which may be present in some subbands during the beam weight observations. This issue has been improved since December 2019; subbands strongly affected by RFI have since been flagged more efficiently and beam weight values are interpolated from the nearby RFI free subbands. Beam weights for Apertif are measured at the start of every observing run, which is typically once in every two weeks. With every new beam weight measurement, the properties of the CBs can change depending on the health of the antenna elements, the health of the system, and the RFI environment during the beam weight measurement.

Deriving CB maps
To measure the shape of the response for each of the CBs for the purposes of this study, we used two independent methods. The first is based on performing drift scans on a bright continuum source and reconstructing the beam shapes from the measured autocorrelations. The second method is based on comparing the flux of continuum point sources in the Apertif survey observations with the measured flux in an already existing, well-calibrated catalogue, such as the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) and constructing the beam shape with a Gaussian Processes regression method (see Section 4.5).
During the drift scan measurement, the PAF is at a fixed position on the sky and the observed source drifts through the field of view in a straight line. With the equatorial mount of the WSRT drift, scans can only be performed in one direction along the RA axis. The separation between the individual drifts is 0.1 deg in declination. To scan the whole field of view of the 40 Apertif CBs (up to the 10 % level of all the edge beams) 33 individual drift scans are needed that altogether form a full set of drift scans 2 . Figure 1 illustrates this process, where the numbered grey circles represent the CBs with a diameter of 0.6 deg, and the blue lines represent the individual drift scans. The start of each drift scan is marked by a grey triangle.
We performed drift scan measurements periodically on a compact, bright radio source: Cygnus A (CygA) or Cassiopeia A (CasA). CygA has an integrated flux of 1589 Jy (Bîrzan et al. 2004) at 1.4 GHz, Cas A has an integrated flux of 2204 Jy (DeLaney et al. 2014) at 1.4 GHz. Both of these sources are marginally resolved by a single WSRT dish, however they are the only sufficiently bright and compact sources available to perform drift scan measurements. Figure 3 shows the continuum images of Cyg A and Cas A from the Dwingeloo 820 MHz survey 3 (Berkhuijsen 1972). The Dwingeloo telescope has the same diameter (25 m) as the WSRT dishes; hence, the angular resolution is comparable to the Apertif autocorrelation data. When constructing the beam maps, we used the autocorrelation data from the telescope. The cross-correlation data can not be used to measure the CB shapes because the fringe rotation of Apertif is faster for the long baselines than our time sampling (10 second). This means that we only detect signal on the short baselines but not on the long baselines during the drift scan measurements. Decreasing the time sampling for the observations is not possible in the imaging mode of Apertif. In addition, a finer time sampling would be impractical, considering the already large data volume of the drift scans. A full set of drift scans on Cyg A is ∼15.5 TB and a full set on Cas A is ∼22.5 TB. Since we only used the autocorrelations for measuring the CB shapes, we reduced the data size to 3.1 TB and 4.5 TB for the Apertif Long Term Archive (ALTA) by deleting the cross-correlations from the drift scan measurement sets. The archived drift scan autocorrelation data are available to the community upon request.
We performed drift scan measurements approximately once per month because of their time-consuming nature. One set of drift scans on Cyg A takes ∼14 hours and one set on Cas A takes ∼ 20 hours. The difference in observing time and data size is due to the different declination of the sources (δ CygA = 40.73 • , δ CasA = 58.81 • ). Drift scans take longer on sources at higher declination. The relatively long observing time also means that a full set of drift scans needs to be observed in at least two separate observing windows, since the equatorial mount of the telescopes can only track sources for 12 hours in a day.
A&A proofs: manuscript no. 44045corr  (Berkhuijsen 1972). Both images are 10 × 10 • with a spatial resolution of 1.2 • . Contour levels are 10, 20, 30, 40, 50, 60, 70, 80, and 90 K. MHz to better avoid the RFI. Accordingly, for those observations we updated the CB map range to 1300-1520 MHz, the first 80 MHz of the bandwidth are continued to be the flagged due to RFI. We constructed the CB maps from the autocorrelations with the following steps: First we extracted the autocorrelation data for each CB from each drift scan measurement set. We extracted all measured polarisation data (XX, YY, XY, YX) in 10 frequency bins (1050 channels -12.8 MHz -per bin), for all individual antennas and averaged for all antennas as well. For data taken after January 2021, at the new observing frequency, we used 18 frequency bins (1000 channels -12.2 MHz -per bin). Next, from the frequency-binned autocorrelation data, we constructed FITS images for all 40 CBs. This step is done by combining the individual drift scans and gridding them onto a comoving RA -Dec grid with the observed continuum source. For this purpose, we used the scipy package: interpolate.griddata with cubic interpolation. Maps were generated for XX, YY, I, and XX-YY polarisation and cover 3.36 × 2.3 deg ( Figure 4). Finaly, spline smoothing is applied to each CB map (to minimise the effect of RFI) and a smaller 1.1 × 1.1 deg (40 x 40 pixel) beam map, centred on each CB, is written into the FITS files ( Figure 5). These smaller CB maps are designed to use for primary beam correction and mosaicking the Apertif images and image cubes. They approximately cover the CBs to the 10% level. The spline smoothing is done with the interpolate.RectBivariateSpline module of scipy, which performs a bivariate spline approximation over a rectangular mesh to smooth the data. We only produced these smaller beam maps for 9 or 17 frequency bins, since the first frequency bin (at 1.3 GHz) is almost always strongly affected by RFI.
The Python code to produce the CB maps is publicly available on GitHub 4 . To refer to specific drift scan data sets, we used CB map IDs, which are based on the observing date of the drift scans in the format YYMMDD (e.g. 190821, 201028). In this work, we use CB maps 190821, 190912, 200130, and 201028 for illustration. The 190912, 200130, and 201028 data sets are good overall representatives for most of the measured CB maps and, in addition, the 190821 data set is chosen to illustrate the effects of broken PAF elements. Figure 4 shows an example of CB maps for the data set 201028 CB 0. We show CB 0 here, since it is the central CB of the Apertif footprint, directly along the optical axis of the dish (van Cappellen et al. 2022), and is the most symmetric shaped CB. From top to bottom, the figure shows I, XX, YY, and XX-YY polarisation CB maps for the ten frequency bins. The spatial extent of the individual maps is 3.36 × 2.3 deg. The image is scaled to highlight the pattern of the side-lobes. The four-fold symmetry of the CBs, particularly visible in the side-lobes, can be attributed to the support structure of the receiver. A somewhat similar beam shape was also measured for the previous MFFE 1.4 GHz receiver on the WSRT by Popping & Braun (2008). Figure 5 shows the 1.1 × 1.1 deg smoothed CB maps for the same data set (201028) for all beams at 1.36 GHz. The beams in the centre of the footprint are relatively symmetric, while the edge beams are quite elongated due to coma distortion and fewer contributing antenna elements towards the edges of the PAF.
We wish to emphasise that the use of the classic WSRT primary beam correction is not appropriate for Apertif. Figure 6 shows one set of measured CB maps (201028) divided by the previous analytic WSRT primary beam model (cos 6 (cνr), where c=68, ν is the observing frequency in GHz and r is the distance from the pointing centre in degrees). In addition to the elongated shapes (and offsets) visible in the outer CBs, the Apertif CB value is generally smaller than the classic WSRT primary beam value. The Apertif PAF illuminates the reflector dish more uniformly than the old MFFE frontends, which leads to an increased aperture size and a smaller CB shape (see Figure 36 in van Cappellen et al. 2022).
Drift scan maps can be used to evaluate some of the polarisation properties of the Apertif system. Figure 7 and the bottom panel of Figure 4 show the leakage properties of the system, where |XX-YY| corresponds directly to the Stokes Q leakage. We find leakages on the order of 0.001 at the beam centres and 0.02 at at the 0.5 response level showing four-fold symmetry. For some of the edge beams, the leakage can be higher, up to 0.1 at the 0.5 response level of CBs 01 and 39, which have generally the most distorted CB shape. Analysing the polarisation measurements in the Apertif imaging observations, we find that the fractional polarisation of polarised continuum sources stays constant at an approximate level of 1% up to a CB response of 0.3. A more detailed analysis of the polarisation characteristics of the Apertif system is presented in Adebahr et al. 2022.  There is a systematic variation in the position of the peak intensity of the CBs compared to the nominal centre of the beams. This offset is on average 1.5 , but can be up to 5 for the edge beams. Figure 8 shows the size and the direction of the peak intensity offsets from the nominal beam centre for all drift scan data sets in red and for the GP CBs in black. These offsets are intrinsic to the CBs and are present in the maps measured with both -GP and drift scan -methods. The offsets correspond to the elongated shapes of the CBs, that is, the more elongated a CB, the larger the offset in the opposite direction of the elongation. This is also visible in Figure 9, where the nominal centre of the beam, indicated with a black dot, is offset from the most sensitive part of the beam, indicated with a red dot. The amplitude of the offset is stable for most CBs, however the direction of the offset does vary for some CBs. The CBs most affected by these variations are the ones in the left centre part of the footprint, which are also the CBs that have been the most affected by broken antenna elements during the Apertif survey observations. This also indicates that the CB shape variation over time is largely due to broken and repaired antenna elements.

CB maps for individual dishes
We derived CB maps for individual dishes in the same way as for the whole array. Figure 9 shows an example of CB maps at 1.36 GHz derived for individual antennas for the 190821 data set. The plots show that some antennas have distorted CB shapes compared to the same CBs on other antennas, for instance, beam 33 on antenna RT6, beam 18 on antenna RTB, or beam 15 on antenna RTC. This is generally due to faulty elements in the PAF on some of the antennas (see also Figure 2). Faulty elements are periodically repaired, which improves the beam shapes. Distorted beams can also be caused by strong RFI during the beam  weights measurement. This can especially be a problem for the low frequency part of the observing band. In addition to distortions the maximum sensitivity point of some beams can also be shifted. The degree and direction of this offset varies slightly between the different dishes. This is due to the independent beam weights for the PAFs on the individual dishes. The per antenna CB maps can be used with advanced imaging algorithms (e.g. WSClean, Offringa et al. 2014) to improve the image quality of the data products. This can be particularly helpful with DDEs caused by broken elements on the PAFs.

Beam size
To characterise the properties of the CBs, we fit a 2D Gaussian to each CB map. Figure 10 shows the FWHM of the fits as a function of beam number. Blue and orange lines show the FWHM measured along the vertical (Dec) and horizontal (RA) axis and the black line shows the average of the two. There is a line representing each frequency bin. The plot shows that CBs in the centre of the field of view (beams 0, 24, 17, 18) are relatively symmetric with similar sizes in the RA and Dec direction while the corner CBs (beam 1, 7, 33, 39) are the most asymmetric and diagonally elongated, where the FWHM in different cross-sections can differ by ∼ 12 . The elongation of the beams is due to the coma effect and that beams at the edge have fewer contributing elements from the PAF compared to the beams in the centre.

Frequency dependence
CB sizes change linearly with frequency. The expected frequency response is θ = 1.22 λ/D, where θ is the beam size, λ is the wavelength, and D is the diameter of the telescope. The measured average frequency dependence for Apertif CBs is: where θ is the FWHM in arcminutes and ν is the frequency in Hz. This is based on fitting a 2D Gaussian to each CB map at each frequency, taking the average FWHM from the 2D Gaussian fit and then fitting a first-order polynomial to the FWHM versus frequency. The results were then averaged for all beams in 16 different drift scan measurements. The uncertainty of the parameters is estimated with the standard deviation between the different data sets. The theoretical expected (for a 25m dish) slope and intercept of the relation are: -2.7 ×10 −08 and 73.75; the slight difference in the average of the measured numbers is due to the more effective illumination of the dish. Figure 11 shows the average beam size for each 40 CBs as a function of frequency for a set of drift scans (190912). The dashed black line shows the average fitted line to the data and the red dotted line shows the theoretical expectation. Some of the beams occasionally show non-smooth variations with the beam size. In most cases, the cause for this is RFI in certain frequency bins. Beam weights are calculated for each individual subband, which means that the beam shape can be slightly different for each of them, especially if strong RFI is present during the beam weights measurement. Unfortunately, it is not feasible in practice to calculate CB maps per subband because of the insufficient S/N of the data.

Time variability
The beam shapes derived from drift scans have average FWHM sizes ranging from 33.5 to 40.4 at 1.35 GHz. Over time the average FWHM typically varies by less then 1 , which is at the few percent level, for CBs in the inner part of the footprint. For CBs at the edge of the footprint, with significantly distorted beam shapes, the average FWHM can vary by up to 2 (∼ 5%). Correspondingly, the median primary beam correction for the CBs typically varies by a few percent between different drift scan measurements. Figure 12 shows the average FWHM size for 16 different drift scan measurements in 2019 and 2020 at 1.36 GHz. The beam sizes show relatively small variation for most of the beams, ∼ 1 -2 . These can be due to different system conditions, for example broken elements, different RFI conditions, different ambient temperature, or small random variations in the electronics of the telescope.

GP CB maps
The GP CB maps are produced by comparing the flux of continuum point sources in the Apertif survey observations with the measured flux of the same sources in the NVSS survey and then constructing the beam shape with a Gaussian Processes regression method. We used the scikit-learn gaussian_process python library to implement a GP kernel that is a combination of two squared-exponential functions and white noise. One of the squared-exponential functions describes the overall beam shape, while the other accounts for small scale deformations. This method requires a certain minimum number of sources for each CB to accurately capture the full CB shape. In practice, this means that a minimum of 5 to 7 survey pointings (containing about 500 cross-matched continuum sources) are needed to reconstruct the beam shape. There are typically between 6-12 individual pointings observed with the same beam weights. This means that the GP method can be used to reconstruct the CB shapes for each set of beam weights used for imaging survey observations. For simplicity, the Apertif DR1 is accompanied by a single set of GP CB maps that are produced from all the highquality observations released in DR1. These maps represent the average CB shapes for the first year of imaging data in DR1. In this work, we present results from this single set of GP CB maps. Considering that we find only a few percent variation in the size of the CBs over this time period (Section 4.4), the average GP CB maps are good representations of the Apertif CB shapes. Further details of the calculation of GP beam shapes are discussed in the Apertif DR1 documentation 5 and in a separate publication (Kutkin et al. 2022).

Comparison of drift scan CB maps to GP CB maps
The comparison of CB maps derived with two independent methods, using drift scans and a GP regression method, gives us the opportunity to better understand the CB shapes of Apertif. Figure 13 compares the CB maps derived with the two different methods. The left panel of Figure 13 shows a set of drift scan maps in grey scale overlaid with contours of the drift scan data (blue) and the GP map (orange). The two sets of contours match each other very closely, hence for a more detailed comparison the right panel of Figure 13 shows the ratio between the drift scan maps (1.36 GHz) and the GP maps. Overall the beam shapes and the average beam sizes agree well between the drift scan CB maps and the GP CB maps. Figure 14 shows the distribution of the pixel ratio between the two beam maps (drift scan map (1.36 GHz)/GP), where the drift scan maps are scaled to the same pixel size as the GP maps. The peak of the distribution is 1.03, demonstrating that the two methods produce consistent CB shapes.
The two methods are highly complementary to each other. The drift scan method has the advantage that it can measure CB shapes for I, XX, and YY polarisation of the Apertif data, for all individual dishes and at several different frequencies. The disadvantages of this method are that the time-consuming observing method makes it impractical to measure beam maps for each new set of beam weights and the bright continuum sources that are used for the measurements are not point sources for the individual WSRT dishes. These two things can introduce a bias in the CB maps. The advantages of the GP method are that it only depends on the NVSS catalogue and that CB maps can be derived for each set of beam weight measurements by using the Apertif continuum images observed with the same beam weights within a two-week time period. This way the CB maps can reflect the actual status of the system at the time of observing. The disadvantage of the GP method is that it can only be used for the broad band continuum images and cannot be derived for different frequencies within the Apertif observing band, for the different polarisation products or for individual dishes. In conclusion, the GP maps can be used for primary-beam corrections of the multi-frequency continuum images that the Apertif pipeline produces. In addition, the drift scan maps can be used in combination with the GP maps to carry out frequency-dependent primary-beam corrections for spectral line data cubes or polarisation data cubes; they can be used with advanced imaging algorithms (e.g. Offringa et al. 2014) for antenna-based deconvolution and primary-beam correction. In addition, they can help identify problems with Apertif observations, for instance, badquality CBs on certain antennas at certain times. Figure 15 shows maps of Apertif/NVSS flux ratios derived with three different sets of drift scan-based CB maps and a set of the average GP CB maps aimed at evaluating how representative the different beam maps may be for all of the Apertif DR1 data. The plotted flux ratios are derived by first primary beam correcting all continuum images in the Apertif DR1 and then crossmatching point sources in the corrected Apertif images with the NVSS catalogue. There is clearly either a bias or gradient in the flux ratios for most beams in the drift scan data sets, with an average flux ratio (Apertif/NVSS) of 1.2. The dominant direction of this gradient is different for measurements on Cyg A (higher flux ratio towards the south-west) and on Cas A (higher flux ratio towards the south-east), as well as the GP beam maps (slightly higher flux ratio towards the south-east). For the GP beam maps, this gradient is much milder compared to the drift scan-based beam maps. This bias is partially due to the fact that neither CygA nor CasA are point sources. Figure 3 shows continuum images of Cas A and Cyg A from the Dwingeloo 820 MHz continuum survey (Berkhuijsen 1972). In addition, both Cyg A and Cas A are close to the Galactic Plane, with diffuse large-scale continuum emission near both sources. In Figure 3, the map of Cyg A shows some extended continuum emission towards the south-west and the map of Cas A has its peak emission towards the south from its nominal coordinates. To try to mitigate this issue, we also performed drift scan measurements on 3C 295, which is one of the primary calibrator sources used for Apertif imaging observations; however, due to the significantly lower flux (∼ 22 Jy) of this source, we were not able to produce highquality CB maps from the observations.
Another cause for the gradient in the Apertif/NVSS Flux ratios in Figure 15 is the quality of the beam weights at the time of the drift scan observations. As discussed in Section 4.6, drift scan-based CB maps represent the beam shapes for a given beam weight measurement and they depend on the health of the system at a given time. This means that if the quality of certain CBs is not adequate at the time of the drift scan measurements, these maps would be unsuitable for PB correction for observations with good-quality CBs. Examples of low-quality beam weights are CB 12, 13, 32, and 39 shown in the top left panel of Figure 15 (beam models 190912). These CB maps show a strong Apertif/NVSS flux ratio bias compared to other CBs in the same drift scan measurements and also compared to the same CBs in other drift scan measurements. This can also account for the smaller gradient in the 200130 data set compared to the 190912 data set, since the quality of the beam weights was substantially better after December 2019 with the improved RFI mitigation strategy. The flux ratio maps also indicate that CB maps from drift scan measurements are generally only suitable for primary beam correction of observations that were taken within a few weeks of the drift scan measurement. On the other hand, science observations with the same low-quality CBs will have a larger uncertainty in their measured fluxes when using the GP CB maps for primary beam correction.

Conclusions
We used drift scan measurements to map the primary beam response of each of the 40 digitally formed CBs from Apertif. We derived beam maps for all individual beams in I, XX, and YY polarisations in ten frequency bins. We derived beam maps for the whole Apertif system and also for all 12 individual dishes. We summarise our findings as follows.
1. The beams towards the edge of the Apertif field of view are more asymmetric compared to the centre beams and are diagonally elongated. This effect is the strongest for the corner beams (CB 1, 7, 33 and 39). We wish to emphasise that the use of the classic WSRT primary beam correction is not appropriate for Apertif. 2. There is a systematic trend across the Apertif field of view that the CB peak intensity is offset from the nominal beam centres. These offsets correspond to the elongated shape of the CBs and vary slightly between the different dishes. 3. We measured Stokes Q leakages on the order of 0.001 at the beam centres and 0.02 at at the 0.5 response level of the CBs. 4. We measured the frequency dependence of the Apertif CB sizes. 5. We investigated the time variability of the FWHM of the CBs and found that for the inner CBs the standard deviation is typically less than 1 and for the edge CBs it is approximately 2 . This corresponds to a 3-5% variation in the average FWHM of the CBs.
6. We compared CB maps derived with two independent methods: using drift scans and a Gaussian regression method comparing continuum fluxes measured with Apertif and in NVSS. We find that the overall shapes of the CB maps derived with the two different methods are in good agreement. 7. The CBs derived from drift scans are essential for measuring the polarisation leakage and the frequency dependence of the CBs. However, GP beams serve as a better representation of the average shape of the beams, since they do not depend on measurements on a single bright continuum source and can be derived for each individual beam weight setup.  A&A proofs: manuscript no. 44045corr Fig. 15. Apertif/NVSS Flux ratio maps with four different sets of compound beam maps (190912, 200130, 201028 (Cas A) averaged Gaussian maps). The black dashed circle has a radius of 36 , which depicts the average FWHM for the Apertif compound beams. The drift scan based maps clearly show on average higher flux ratios to the south-west side of the beams.