High-dimensional inference of radio interferometer beam patterns I: parametric model of the HERA beams

ABSTRACT

Accurate modelling of the primary beam is an important but difficult task in radio astronomy. For high dynamic range problems such as 21 cm intensity mapping, small modelling errors in the sidelobes, and spectral structure of the beams can translate into significant systematic errors. Realistic beams exhibit complex spatial and spectral structure, presenting a major challenge for beam measurement and calibration methods. In this paper series, we present a Bayesian framework to infer per-element beam patterns from the interferometric visibilities for large arrays with complex beam structure, assuming a particular (but potentially uncertain) sky model and calibration solution. In this first paper, we develop a compact basis for the beam so that the Bayesian computation is tractable with high-dimensional sampling methods. We use the Hydrogen Epoch of Reionization Array (HERA) as an example, verifying that the basis is capable of describing its single-element E-field beam (i.e. without considering array effects like mutual coupling) with a relatively small number of coefficients. We find that 32 coefficients per feed, incident polarization, and frequency, are sufficient to give percent-level and |$\sim$|10 per cent errors in the mainlobe and sidelobes, respectively, for the current HERA Vivaldi feeds, improving to |$\sim 0.1{{\ \rm per\ cent}}$| and |$\sim 1{{\ \rm per\ cent}}$| for 128 coefficients.

1 INTRODUCTION

Large radio interferometer arrays are increasingly being used for high dynamic range (HDR) applications, where the signal of interest is several orders of magnitude fainter than other confounding signals (Paciga et al. 2011; Tingay et al. 2013; van Haarlem et al. 2013; DeBoer et al. 2017; Santos et al. 2017; CHIME Collaboration 2022). HDR observations place extremely stringent requirements on the precision and fidelity of array calibration; whereas relatively small (e.g. percent-level) errors in calibration may be tolerable for ‘traditional’ applications such as imaging bright sources, they can cause catastrophic leakage and misidentification of bright signals into parts of the data (e.g. particular regions of Fourier space) where we might otherwise have hoped to extract the faint target signal. HDR applications are challenging for standard calibration methods, and a great deal of effort has been put into methodological development to overcome their limitations and provide calibrations at the required level of accuracy (Dillon et al. 2020; Ewall-Wice et al. 2022; Byrne 2023; Charles et al. 2023; Cox et al. 2023).

A prominent example is the detection of brightness temperature fluctuations of the redshifted 21 cm line from neutral hydrogen, commonly known as 21 cm intensity mapping (IM). In this application, the 21 cm fluctuations are expected to be around the |$\sim$| mK level, compared with Galactic and extragalactic synchrotron emission (‘foregrounds’) ranging from tens of Kelvin in brightness temperature at higher frequencies (⁠|$\sim 1$| GHz) to thousands of Kelvin at lower frequencies (⁠|$\lesssim 200$| MHz) (see e.g. Mesinger 2019; Liu & Shaw 2020). This represents a typical dynamic range of |$\sim 10^4\!-\!10^5:1$|⁠. In total intensity (Stokes I), the foreground emission is expected to have a smooth, power-law-like frequency spectrum, as compared with the rapidly fluctuating spectrum of the 21 cm fluctuations, which in principle should allow the signals to be separated using standard tools such as Fourier filtering, principal component analysis, and the like (Morales & Hewitt 2004; Morales et al. 2012). Calibration errors and artefacts greatly complicate this picture however; even relatively small errors can couple bright, foreground-dominated Fourier modes into the 21 cm signal-dominated modes, swamping the signal (Barry et al. 2016; Byrne et al. 2019; Orosz et al. 2019). We require instrumental calibrations at a level better than the dynamic range, of the order of |$10^{-5}$|⁠, to keep the leaked foregrounds below the 21 cm signal (Barry et al. 2016; Thyagarajan et al. 2016). This is likely to require knowledge of the instrument or a reference sky model at least a couple of orders of magnitude better than what is currently possible (Shaw et al. 2015; Ewall-Wice et al. 2017).

While there are many aspects of an interferometer array that need calibration, these can largely be combined into two types of complex gain parameters – direction-independent gains, which describe the frequency- and time-dependent degrees of freedom such as the bandpass, and direction-dependent gains, which incorporate the change in sensitivity with angle due to each receiving element’s antenna pattern (beam), also as a function of frequency and time. Direction-independent gain calibration has been studied extensively in the context of 21 cm IM experiments (e.g. Ewall-Wice et al. 2017, 2022; Byrne et al. 2019, 2021; Dillon et al. 2020; Byrne 2023; Charles et al. 2023). In this paper we focus on direction-dependent calibration, and in particular inference of the per-antenna complex (E-field) beam.

Beam estimation is difficult because the beam itself has a large dynamic range, a complicated dependence on direction and frequency, and multiplies a sky that is also subject to substantial modelling uncertainties. Even antennas that are designed to have large directivity, such as the parabolic dish reflectors of the Hydrogen Epoch of Reionization Array (HERA), exhibit non-negligible sidelobes far from the centre of the beam pattern. The (peak-normalized) maximum sidelobe level of the HERA Vivaldi feed design is in the region of |$-15$| to |$-25$| dB depending on frequency, for instance (Fagnoni et al. 2021a), with substantial variation with frequency. Modelling these sidelobes accurately is a necessity if one hopes to reach a dynamic range of |$10^{-4}$| (i.e. |$-40$| dB) or better. The beam pattern is also expected to vary somewhat between different elements of the array, as slight variations in feed position and alignment, reflector geometry etc. will always occur (Orosz et al. 2019; Choudhuri, Bull & Garsden 2021; Kim et al. 2022, 2023). These variations may be exacerbated by environmental conditions, for example different wind loadings or air temperatures at different positions within the array. Electromagnetic models and lab-based beam measurements (which usually cannot be performed in the far field) are typically unable to account for these complicated variations in a robust manner, and so measuring the beam patterns of the array elements in situ is necessary.

A variety of approaches to in situ measurements have been tried. One is to keep track of the autocorrelation (zero-spacing) signal as the sky rotates through the beam as it points in a fixed direction. Using a sky model, one can try to reconstruct a model for the beam that accounts for the antenna temperature variations as different parts of the sky (e.g. bright sources) drift through different parts of the beam. This is generally difficult, as sky models are neither complete nor calibrated at the required level of accuracy. A similar approach uses the visibilities from pairs of correlated antennas, again with a sky model being used as a reference (Pober et al. 2012; Nunhokee et al. 2020). The two antenna beams modulate one another however, and other aspects of the interferometer (such as direction-independent gains) compound the measurements, in addition to the familiar issue of the sky model being incomplete. To try and get around the incompleteness problem, attempts have also been made to use bright artificial sources such as drone-mounted radio transmitters to map the beams (e.g. Virone et al. 2014; Pupillo et al. 2015; Jacobs et al. 2017), or bright emission from satellites (Line et al. 2018; Chokshi et al., in preparation). Other methods such as holography and a combined photogrammetry plus detailed EM modelling approach have also been attempted (Berger et al. 2016; Iheanetu et al. 2019).

In this paper series, we describe and demonstrate a parametric Bayesian method for inferring the per-element beam from observed visibilities. Bayesian instrumental characterization methods have been demonstrated before in the context of radio astronomy (e.g. Lochner et al. 2015; Yatawatta 2018; Sims, Pober & Sievers 2022a, b; Anstey 2023; Cumner et al. 2023; Sims et al. 2023). In this paper (Paper I), we focus on the model-building aspect of the inference, and in the second paper (Paper II) we demonstrate the feasibility of the inference on a toy problem. In particular, we build a flexible linear model based on a set of analytic basis functions that describe the simulated E-field beams of the HERA receiving elements with relatively few free parameters. An accurate and efficient basis that respects the symmetries of the problem is essential for avoiding significant model errors that could introduce spurious spectral and spatial structure. Achieving this in a compact way (i.e. with the fewest basis functions/coefficients possible) improves the prospect of directly inferring the structure of the beams from observations – particularly if this is to be done for each antenna within the array (rather than as an array average), out to the far sidelobes, and/or as a function of frequency. Conversely, a more ‘cautious’, overly flexible basis would be harder to constrain with typical data, essentially only allowing us to marginalize over a very broad range of possible beam structures. The basis we propose in this paper has many of the properties required for a successful inference scheme in which individual antenna patterns can be measured versus frequency, and is tuneable in terms of its precision, at least with respect to the complex EM simulations that we use as a reference. Note that the methods we use in Paper I to determine this ‘sparse basis’ are not themselves Bayesian; the basis is in support of a Bayesian Gibbs sampling framework expounded in Paper II.

Analytic parametrization of the beam is an active area of research in the HERA collaboration (Choudhuri et al. 2021).¹|$^{,}$|  ² In general, analytic parametrizations have demonstrated significant promise in solving problems in HDR applications, e.g. more rapid beam evaluation (Asad et al. 2021), better physical beam characterization (Nasirudin et al. 2022; Sekhar et al. 2022), and removal of calibration biases that stem from gridding errors (Barry & Chokshi 2022). This last reference also highlights the importance of accurate beam models for the sake of accurate gridding kernels in imaging power spectrum analyses (Morales et al. 2019).

We introduce a basis set that, to our knowledge, has not been used in the beam modelling literature. This basis solves several problems for us where other similar solutions from the literature did not.

• It is fully analytic (unlike principle component analysis of holographic measurements for example). This makes it highly portable between different pixellization schemes while avoiding problems associated with discretization issues (e.g. spurious fluctuations in apparent source brightness at the horizon due to the setting of discrete pixel centres, or issues associated with sampling on sub-pixel scales for arrays with long baselines relative to the beam model resolution).

• It describes the complex, realistic HERA beam simulations with relatively few basis functions compared to alternatives.

• It behaves well at a coordinate singularity where e.g. polynomial-based bases did not (see Section 2).

• The special functions involved have numerically stable scipy implementations up to high order which allows for thorough exploration of the parameter space.

We implement the sparse basis construction code of Paper I as well as the Gibbs sampling code of Paper II as a part of the hydra Gibbs sampling framework (Kennedy et al. 2023).³ The sparse fitting code makes use of the UVBeam class in the pyuvdata  ⁴ package (Hazelton et al. 2017). This class provides a general python interface for a variety of beam file formats, meaning the code we provide can be readily applied to other experiments.

In what follows, we specialize to interferometric arrays operating in a drift-scan configuration (i.e. always pointing at the zenith), similar to HERA. Using previous EM simulations of the HERA Phase I and Vivaldi feeds, we show how an efficient basis can be constructed that allows accurate modelling of uncertain E-field beams with a relatively small number of parameters (Section 2-3). We then explore whether any spurious spectral structure arises as a complication of choosing such a sparse spatial basis (Section 4). Lastly, we conclude (Section 5).

2 MATHEMATICAL FORMALISM

In order to incorporate the beam into our Bayesian inference, we will need to posit a parametrized model for the beam that relates to the data from which the inference is drawn. In the general landscape of Bayesian inference, this is a rather arbitrary decision moderated only by the fact that some parametrizations may describe the data more efficiently than others. Without extreme prior knowledge about the beam, and in particular, uncertainties in the beam due to physical variations from the ideal scenario, we are left to choose some arbitrary basis for the beam with the hopes that our existing prior knowledge can at least constrain which basis functions are the most dominant. As a concrete example, we can take a detailed electromagnetic simulation or highly precise holographic measurement and fit a linear basis to the beam pattern for a high number of modes. Depending on the beam pattern and choice of basis functions, we may find that most of the fit is dominated by a small subset of basis functions. This has been shown to be relatively effective in de Lera Acedo et al. (2013), Young et al. (2013), Bui-Van, Craeye & de Lera Acedo (2017), Iheanetu et al. (2019), and Asad et al. (2021), though the choice of basis is often idiosyncratic to the telescope. We can then restrict ourselves to these dominant basis functions and consider the coefficients as the free parameters in our model.

An appealing alternative, though one that seems to us extremely computationally demanding for a full Bayesian computation, would be to have not only detailed knowledge of the ideal receiving element, but also have a parametrization for physical perturbations to this ideal along with detailed simulations or calculations of the response to them. Then a few parameters describing the perturbations only need to be constrained. See, for example, the characteristic basis function pattern approach in Maaskant et al. (2012) and Young et al. (2013). The primary issue is that perturbations tend to belong to a multidimensional continuum and developing simulation outputs for finely graded perturbations to understand the resulting effect on the beam pattern is computationally prohibitive, though not infeasible depending on the study design (Kim et al. 2022, 2023). If a suitable emulator were constructed (i.e. if we could apply machine learning to produce fast but necessarily incomplete outputs of the simulator), the computational overhead could be reduced in the long term with an initial investment. However this might violate intellectual property law and is a somewhat ill-defined problem for a general electromagnetic simulator. Due to these seemingly intense challenges, and since our goal is mainly to establish the tractability and usefulness of our Gibbs sampling method, we opt for a more arbitrary linear basis expansion.

We aim to infer the beam parameters from the visibilities via a forward model. We write the visibility equation in terms of the polarized beam response for antenna j (Jones matrix, ignoring other propagation effects such as ionospheric refraction), |$\boldsymbol{\mathcal {J}}_j(\hat{s}, \nu )$|⁠, sky coherency matrix, |$\boldsymbol{\mathcal {C}}(\hat{s}, \nu )$|⁠, and physical antenna separation for antennas j and k, |$\Delta \vec{x}_{jk}$|⁠, as,

where c is the speed of light, |$\nu$| is the frequency of observation, |$\hat{s}$| is a unit vector pointing to different positions on the sky, |$i^2=-1$|⁠, and depending on the coordinate system, some of these quantities are time-dependent. Different choices of coordinate system present different trade-offs. The Bayesian inference for the beam is drastically simpler if we operate in a coordinate system where the beam pattern is fixed. For the sky coordinates, we choose azimuth, |$\phi$| and zenith angle, |$\theta$|⁠. We choose sky polarization directions that align with this coordinate system. This puts all of the time dependence in the sky coherency matrix at the cost of producing a discontinuity in the beam pattern at zenith. This discontinuity arises because in this coordinate system, there is no well-defined direction for the polarization unit vectors at zenith. In math, if we think of the polarization unit vectors, |$\hat{\theta }$| and |$\hat{\phi }$| as vector fields, then

for |$\phi _0 \ne \phi _1$|⁠, and similarly for the azimuthal unit vector. See Byrne et al. (2022) for an illustration of this coordinate-based phenomenon, as well as an excellent discussion of state-of-the-art polarized imaging techniques.

For a pure drift scan observing strategy, this is relatively simple to deal with. In particular, we can choose a basis that naturally incorporates this discontinuity. We single out a Jones element and write it as

where p indexes instrumental polarization and |$p^{\prime }$| indexes incident polarization of incoming radiation. If we can find a complete, orthogonal basis such that |$f_n(0) \ne 0$| for all n, then every single basis function will have this discontinuity at zenith. This requirement makes it such that a best-fitting function for some finite number of basis modes does not fill in the discontinuity with 0, as we have seen happen when using Zernike polynomials on projections of the sky to the unit disc.⁵

Fortunately there exists such a basis that is fairly natural to use with a dipole radiator. If we project the sky to the unit disc, we can use a Fourier–Bessel series to capture the radiation pattern. From our observations, the azimuthal (Fourier) modes capture the azimuthal structure of the beam patterns in question with relatively few modes. The basis can be made complete and contain the discontinuity in every basis function by using 0th order Bessel modes and demanding that a linear combination of the basis function and its derivative go to 0 at the horizon (Jackson 1998). Generally, we write

where |$u_n$| is the nth root of the linear combination of Bessel functions specifying the boundary condition. If we just demand that either the basis function or its derivative vanishes, then |$u_n$| is the nth zero of either the 0th order Bessel function |$J_0(x)$| or its derivative, |$J_0^{\prime }(x) = -J_1(x)$| (depending on choice of boundary condition), and |$q_{n}$| is a choice of normalization for the basis so that the modulus square of each basis function integrates to 1 over the disc. Since the beam response at the horizon for azimuthally polarized radiation is non-zero, the choice of Neumann (derivative-vanishing) boundary condition seems more suitable there, however in our numerical experiments we find a better overall fit from the Dirichlet (function-vanishing) boundary condition for a finite sum of basis functions. The normalization constant is given explicitly by

For this approach, we need to pick a projection to the unit disc. While an orthographic projection is a choice with geometrically transparent properties, we instead opt for a projection defined by

where |$\alpha$| is close to 1. This projection has the useful property that

In other words the projection is area-preserving (up to an overall factor of |$2\alpha$|⁠). For |$\alpha =1$|⁠, |$\rho =0$| maps to zenith and |$\rho = 1$| maps to the horizon. In this case, the total volume contributed by a basis function to the square of a Jones element is simply given by

As a function of |$(n, m)$|⁠, we call this the ‘Fourier-Bessel Energy Spectrum’ (FBES) for the beam in question.

Using |$\alpha =1$| and Dirichlet boundary conditions produces large errors near the horizon since the electromagnetic response is non-zero there. Using Neumann boundary conditions allows us to model the horizon better at the expense of significantly worse fits throughout most of the sky. We find that the horizon problem can be better solved by noting that we can essentially extend the beam with zero-padding to zenith angles beyond the horizon since any source beneath the horizon will contribute 0 flux. We then choose |$\alpha$| to be slightly greater than 1, so that the vanishing point of the radial Bessel functions is slightly beneath the horizon. This allows for non-zero values at the horizon and thus greatly reduces errors there. This also improves errors over most of the sky since we are no longer forcing the fit to satisfy a boundary condition that is explicitly disobeyed by the simulation.

To investigate the effectiveness of this basis for the HERA use case, we fit the simulated Phase I dipole and Vivaldi beams from Fagnoni et al. (2021b, a) in a weighted least-squares sense and observe what varying the number of basis functions does to the performance. Mathematically, we write our Jones element as a vector, |$\boldsymbol {j}_{j\hat{\alpha }}^p$|⁠, where each component is the simulated beam evaluated at a pixel on the sky at 1 degree resolution. We then assume

where |$\boldsymbol {B}$| is the design matrix encoding our choice of basis, and |$\boldsymbol {a}$| are the coefficients in that basis. We then calculate a linear least-squares solution for when the number of basis functions is finite. This amounts to minimizing the loss function

This is equivalent to solving for |$\boldsymbol {a}$| in the normal equations

In this work, we do this independently for each frequency. In Section 4, we examine the spectral properties of these fits. In Paper II, we use our Bayesian prior to enforce spectral smoothness. In subsequent sections, we analyse the fit performance under different choices of allowed radial and azimuthal modes.

3 COMPRESSING THE BASIS

A primary strength of this approach is that we choose a non-pixel basis to express the beam. This means that a set of coefficients can easily map to different pixel schemes without the need for interpolation. However, the coefficients make reference to a particular coordinate system, so this does not eliminate the need for careful coordinate transformations. For the inference, this strength is only manifest if the chosen basis can compactly represent the beam, i.e. approximately reproduce the true beam with relatively few coefficients. This is important both from a computational perspective and theoretical inference perspective. We would like to constrain as few parameters as possible given a data set of fixed size, and this problem is magnified by the fact that radio telescope beams have non-trivial frequency structure.

There are two critical concerns when choosing a basis to represent the beam. First, assuming a priori knowledge of the beam through electromagnetic physics and knowledge of the receiving system, does the basis compactly represent the ideal beam assuming the perfect modelling of the receiving system? Secondly, assuming a class of perturbations to the receiving system, are ensuing perturbations to the beam pattern also compactly represented by the chosen basis? In this work, we provide a basis that satisfies the first question for the HERA receiving systems (both the Phase I dipole and Vivaldi feeds) using simulations of the ideal receiving system constructed in Fagnoni et al. (2021b, a), which use CST Microwave Studio (Weiland 1977; Clemens & Weiland 2001). We leave a full treatment of the second question for future work.

We show the peak-normalized simulated beams for the HERA Phase I dipole and Vivaldi feeds at 150 MHz in Figs 1 and 2. This is the middle of the operating range for both feeds. Generally, the spatial structure of the beam appears more complicated at higher frequencies. We see that the receiving elements have an obvious dipolar pattern that dominates, however more complicated azimuthal structure emerges at larger zenith angles. Since the simulated receiving elements have 90 deg rotational symmetry, the N-S dipole responses are a rotated copy of the E-W polarized responses. In order to increase the visibility of plots for the rest of the paper, we will henceforth only show quantities in terms of the E-W dipole response. In practice, this is not true, since the full embedded element pattern has its symmetries broken based on where a given receiving element is in the array (Fagnoni et al. 2021b, a). We display a proof of concept with the simpler model first, and reserve an analysis of the full embedded element pattern for future work.

We fit each beam with 80 radial modes and 91 azimuthal modes and show the residuals in Figs 3 and 4. Residuals are generally less than |$10^{-5}$|⁠, however this is far too many basis functions (7280) to be practical in a Gibbs sampling setting since the coefficients vary as a function of frequency and receiving element. The Phase I dipole residuals appear to be dominated by radial modes all the way to the horizon. The Vivaldi errors are slightly more visually noticeable, but still small. Nevertheless, we are strongly encouraged by the overall size of the errors over most of the sky.

In order to compress the basis, we need a way of determining which modes are the most important to use. Where appropriate data are available, various statistical methods may be used to perform this compression, such as principle component analysis of holographic measurements (Asad et al. 2021), or by inferring direction-dependent effects in the visibilities using a penalized likelihood (Yatawatta 2018). We are opting to solve this compression problem before setting up the visibility based inference (in Paper II), rather than doing it on line during inference with e.g. a shrinkage estimator (Gelman et al. 2013). To do this, we examine equation (8) (the FBES) for the very precise fit using thousands of coefficients, shown in Fig. 5 for both feed types and sky polarizations at 150 MHz. Almost all of the energy is at low n and m for both feeds, suggesting that we can choose a much smaller basis set with relatively little error. Since all of the plotted Jones elements have similar volume, we can interpret the number of high-energy modes as a proxy for the spatial complexity of the beam i.e. how compressible it is. For both feeds, we see that the response to azimuthally polarized radiation has longer tails to high radial modes compared to the response to zenith angle polarized radiation. Additionally, the Vivaldi feed has more higher energy modes than the Phase I dipole. For example the |$|m|=3$| modes are prominent out to fairly high radial power, and there are more significant modes at much higher azimuthal number.

To demonstrate that our intuition about mode significance holds, we perform sparse fits to the beam using only a subset of the entire basis set considered so far. In each fit, we set a certain number of basis functions to use and pick the basis functions in descending FBES order i.e. we select the M highest energy basis functions for some M and fit only to those. We then calculate the root-mean-square error (RMSE) of the fit over the whole sky, which is the same as computing |$\sqrt{L/N_\text{pix}}$| where L is defined in equation (10) and |$N_\text{pix}$| is the number of pixels in the simulated beam. We show the results in Fig. 6, where for each point we increase the number of basis functions by a power of 2. We see that the RMSE decreases roughly as a power law as we include more basis functions. The Vivaldi beam appears consistently less compressible than the Phase I dipole in the sense that we observe a higher RMSE for a given number of basis functions.

Beam errors manifest as flux and phase errors for sources observed at the location of the beam errors. It is difficult to track these flux and phase errors through the complicated analysis pipelines currently in place for power spectrum estimation. The exact error tolerance for 21-cm power spectrum estimation is thus not obviously determined. Ewall-Wice et al. (2017) suggests it should be better than 1 per cent everywhere on the sky. Better than |$\sim 100{{\ \rm per\ cent}}$| errors in the sidelobes would allow this basis to improve on recent drone-based measurements assuming the sparse selection describes the physical perturbations (Jacobs et al. 2017).

Fig. 7 shows the RMSE as a function of zenith angle for both feed types and polarizations and a range of sparsity. We find that with 32 modes, there is |$\sim 1{{\ \rm per\ cent}}$| agreement in the main lobe and |$\sim 10{{\ \rm per\ cent}}$| agreement in the sidelobes. With 512 modes we achieve better than 1 per cent precision at almost all zenith angles. The residuals from some of the sparse fits are shown in Figs 8 and 9 for the Phase I dipole and Vivaldi feeds, respectively. For a given level of sparsity, errors generally range by one order of magnitude above and below the RMSE value depending on the location of interest. We explore computational complexity for our Gibbs sampling pipeline in Paper II more thoroughly, however we expect having a beam model with 64–128 parameters is computationally feasible, while 512 modes will be tractable with a more optimized implementation.

We note that translating the observed errors in this fit to the actual errors produced by this approach in practice is not necessarily straightforward. The actual beam will at best be some perturbation of the simulated beam used in this work, not even including the close-packed array effects, which are known to vary even between different simulation packages on the order of a few percent (Bolli et al. 2023). Perhaps more importantly, the approach in Paper II involves inferring the beam based on the interferometric visibilities, whereas this least-squares fit is more analogous to a problem where one infers the beam based on a holographic map. In other words, this work verifies that the spatial structure of the beam can be captured by this basis with relatively few basis functions, but the errors shown do not necessarily represent the types of errors we may observe with a Bayesian point estimate since the information content of the visibilities will be different than this simulation. Furthermore, our approach will return a full Bayesian posterior, not just a point estimate. Folding uncertainties into a forward modelling procedure should hopefully alleviate some of the difficulties that come with point estimation such as calibration errors (Byrne et al. 2021).

The Fourier–Bessel basis has the advantage that it is sparse enough for our purposes, but probably flexible enough to handle perturbations. An alternative basis can be developed based on singular value decomposition (SVD) of the Jones elements at each frequency. In some preliminary experimentation, we find that this can describe the simulations extremely sparsely, and this is also demonstrated for different instruments in other works (Asad et al. 2021; Cumner et al. 2023). To remove the dependence on a particular pixellization, we can write the SVD basis functions in terms of our Fourier–Bessel basis, thus developing an approximate SVD basis that smoothly interpolates to any desired spatial position. While this produces an extremely sparse representation of a given simulation, the modes uncovered in the SVD are organized in a non-trivial manner e.g. multiple azimuthal modes describing mixtures of dipolar and quadrupolar structures. This makes it less interpretable and therefore disadvantageous from an instrumental characterization perspective compared to the Fourier–Bessel basis. Since sparsity gains from an SVD approach are significant, we aim to explore this option more extensively in the future.

Spectral structure in the beam errors will greatly magnify their harm, particularly for foreground avoidance strategies that hope to separate the EoR and foreground signals by taking advantage of the spectral smoothness of the foregrounds. In addition, errors in the beam can produce errors in calibration, which can themselves obscure observations of the cosmic reionization signal. Indeed, part of the motivation for a Bayesian approach to beam modelling is to reduce the knock-on effect of beam errors on calibration. We explore the spectral structure of the fit beams in the next section using a suite of simulations, but leave the task of examining effects on calibration to future work.

4 SPECTRAL STRUCTURE

We investigate any spurious spectral contamination introduced by fitting the beam with the chosen basis. To quantify this contamination, we compared the delay power spectra (Parsons & Backer 2009; Parsons et al. 2012) from a suite of visibility simulations which use differing numbers of fit coefficients for the beam. All other input parameters to the simulations were kept fixed, however, to isolate the effects of changing the beam.

The delay spectrum is accessed by Fourier transforming (and then squaring) the visibilities along the frequency axis into ‘delay’ space, so-named because the visibilities in this space represent time-like correlations between the antennas’ voltage signals. The intrinsic signal of astrophysical foregrounds is generally confined to delays less than the geometric delay of the baseline, but can be broadened due to instrumental effects (such as a chromatic beam or calibration errors), as well as the choice of spectral tapering function when performing the Fourier transform. Depending on baseline length and other chromatic instrumental effects, delays below a certain scale will be inaccessible without extremely accurate foreground subtraction, while delays greater than this will likely avoid foreground-related contamination (Datta, Bowman & Carilli 2010; Morales et al. 2012; Parsons et al. 2012; Trott, Wayth & Tingay 2012). For HERA, preserving high dynamic range for delays greater than |$\sim$|200–300 ns is essential to measurements of the 21-cm power spectrum during reionization (Parsons et al. 2012; HERA Collaboration 2023). Therefore, the primary purpose of this test is to check whether there is a significant decrease in dynamic range above those delays.

The simulations were performed using pyuvsim  ⁶ (Lanman et al. 2019), a high-precision visibility simulator which has been rigorously tested (Aguirre et al. 2022; Lanman, Murray & Jacobs 2022). We simulated a selection of baselines from a close-packed hexagonal array layout with various orientations and lengths (⁠|$0\le b\le 90$| m). For the sky, we used the 2016 Global Sky Model (GSM; Zheng et al. 2017), a map of galactic diffuse emission. For each antenna feed type (dipole or vivaldi), we ran simulations using varied numbers of beam fit coefficients, from 32 to 512 (‘FB simulations’). For a given number of fit coefficients (⁠|$N_{\text{coeff}}$|⁠), we use the FBES to select the |$N_{\text{coeff}}$| modes with the largest FBES amplitudes. We do this independently at each frequency, and take the union of all basis functions selected in this way across the simulated band. Since the top |$N_\text{coeff}$| contributing basis functions varies slightly as a function of frequency, this means each selected |$N_\text{coeff}$| is less than the number of modes actually used in each instance of the test. Note that this is different in previous sections, where we were only examining one frequency. As a reference against which we can compare the FB simulations, we also ran a simulation for each antenna feed type using the corresponding raw electromagnetic simulations from Fagnoni et al. (2021b, a) (‘Fagnoni simulations’).

This test is greatly complicated by the fact that the reference beam simulations are only available at 1 MHz resolution, which cannot access most of the delays relevant to EoR science. To accommodate a more straight-forward comparison, we fit the Fourier–Bessel beam at 1 MHz resolution and then interpolated the coefficients (as well as the reference beam) during simulation using a cubic spline to a resolution of |$\sim 97$| kHz to match the HERA spectral resolution. The Nyquist frequency for a 1 MHz bandwidth is 500 ns. Delays in the range |$|\tau |\gt 500$| ns are therefore susceptible to interpolation artefacts regardless of which beam is used. In other words, there is no obvious ground truth of what HERA should observe with this sky model above the Nyquist delay. With this caveat in mind, we list three important observations of the delay spectra from these visibility simulations, and then discuss them in more detail below. The delay spectra are shown in Figs 10 and 11 for the dipole and Vivaldi feeds, respectively.

• The Fourier–Bessel beams generally produce slightly smoother visibilities than the simulated beams of the corresponding feed i.e. there is a loss of power at intermediate delays and beyond, meaning there is no hit to dynamic range.

• There is a significant discrepancy below the Nyquist delay, visible in the XX delay spectra, which appear to have pronounced shoulders when using the reference beam. We argue that these are largely an artefact that stems from a complex interplay between spatial interpolation of the beam and the chosen sky model. If this is the case, this structure is spurious and it is desirable that the FB fits should not reproduce it.

• There is also a significant discrepancy at |$|\tau | \sim 1500$| ns in the XX beams, which is above the Nyquist delay. We suspect these are the pronounced shoulders from the reference beams being aliased to a harmonic of the Nyquist delay as a result of the frequency interpolation. This therefore represents high-delay contamination that is happily lacking when using the Fourier–Bessel beams.

To better understand the choice of interpolation spline, we also ran a pyuvsim simulation using linear frequency interpolation of the Fagnoni beam. For linear interpolation, interpolation artefacts were highly visible as a set of peaks (duplicates of the peak at |$|\tau |\le 500$| ns) spaced uniformly every 500 ns (the Nyquist frequency for 1 MHz spectral resolution). While the amplitude of these duplicated peaks decreased with increasing |$|\tau |$|⁠, the dynamic range between |$\tau = 0$| ns and the highest delays was only |$\sim 10^8$|⁠. Whereas for the cubically interpolated Fagnoni simulations, we see a dynamic range of |$\sim 10^{11}$| in Figs 10 and 11 and little evidence of such duplicated low-delay structure.

There are, however, also differences between the Fagnoni beam and FB beam simulations inside the Nyquist delay. The raw Fagnoni beams were simulated on a rectilinear grid in altitude and azimuth with a 1 deg resolution. These beams thus also require spatial interpolation during visibility simulation to get the beam value at the location of each source (pixel) in the sky catalogue. Our FB beam fits are analytic and thus interpolate smoothly to any pixellation scheme. To determine the importance of this spatial interpolation on our results, we also performed a set of simulations using a mock, gridded point source catalogue with one point source at the centre of each Fagnoni beam pixel at a single LST. In this case, the Fagnoni beams require no spatial interpolation. For this gridded point source catalogue, we see almost no discrepancy inside the Nyquist delay. This suggests that the discrepancies we see for |$|\tau |\sim 500$| ns are related to the spatial interpolation of the Fagnoni beam during visibility simulation i.e. it does not represent intrinsic spectral structure of the beam that is subsequently imprinted on any sky model for which it is used. However, since the gridded sky model does not have the same intrinsic spectral structure as the GSM, this additional test by itself leaves open the question as to whether these shoulders are intrinsic beam structure that somehow goes missing when we use the FB beams. By running the simulations of the GSM with an achromatic beam equal to the FB beam at the lowest interpolated frequency in the band, we find that these shoulders are not present, suggesting that the shoulders are not intrinsic to the GSM. In summary, we suspect the large discrepancies between the Fagnoni beam simulations and the best-fitting FB beam (see Figs 10 and 11) are caused by interpolation issues with the simulated beams that the FB basis successfully avoids, rather than being due to a serious reconstruction error on the part of the FB basis.

Figs 10 and 11 also show that there is some variation between the FB simulations as a function of |$N_{\text{coeff}}$|⁠. Most notably, the dynamic range changes slightly as |$N_{\text{coeff}}$| increases. This is most obvious in the XX polarization results for both feed types. This effect, however, is quite small and only appears significant in the fractional difference plots because of the low amplitude of the delay power spectrum of the Fagnoni beam simulations. Overall, the dynamic range is large (⁠|$\sim 10^{13}$|⁠) and roughly consistent for both feed types, both polarizations, and all |$N_{\text{coeff}}$| values. This suggests there is minimal spectral structure imposed on the visibilities by the FB beam fits as a function of |$N_{\text{coeff}}$|⁠.

5 SUMMARY AND CONCLUSION

We have investigated parametric models of the HERA beams as a means of exploring the addition of primary beam inference in Bayesian 21-cm intensity mapping pipelines. Though we have not specifically investigated it in this work, we expect that allowing for variations in the beam model will help mitigate systematic inference errors in other nuisance parameters such as the direction-independent gains, which is a common problem in pipelines that assume an exact, but incorrect, beam model. Using an analytic basis also simultaneously solves a spatial interpolation problem that arises when using pixel-based beams. Namely, since we choose a spatially smooth basis, the interpolated beam values are also spatially smooth. Since spatial modes couple into spectral modes in 21-cm power spectrum estimation, this enhances the quality of such measurements.

We found that using a modified Fourier–Bessel basis in an area-preserving projection of the Azimuth–Zenith–Angle coordinate system to the unit disc provided beam models that were relatively compressible for a given desired performance at zenith. Good performance at zenith is a particularly important feature for a drift-scan telescope such as HERA, while compressibility is important to reduce the computational cost of the inference pipeline. To compress the basis, we formed a Fourier–Bessel Energy Spectrum and determined a fixed number of important modes by choosing the strongest contributors to the spectrum. To assess the effectiveness of this compression, we took the least-squares fit in the compressed basis and compared fit residuals and root-mean-square error over the sky. Without a full test of how such errors affect e.g. direction-independent calibration, it is difficult to know what constitutes a sufficiently complex beam. However, a baseline test of the basis’ effectiveness is to simulate delay spectra using the fit beams with varying number of basis functions, assuming perfect direction-independent calibration, and compare to the same simulation with the simulated beam.

When we perform this simulation, we find that there is no sign of deleterious spurious spectral structure. For both feeds, performance at high delays greater than 2000 ns is nearly identical to the simulated beams. Performance at intermediate to high delays is within the dynamic range requirements set by the reionization signal. We observed the largest differences for both feed types at intermediate delays surrounding the Nyquist delay of the CST simulation (500 ns) and around 1500 ns. The simulated beams exhibit stronger power at these delays. We suppose this is likely to be an interpolation artefact, since structure at delays beyond the Nyquist delay should not be recovered by a spline interpolation except by coincidence, and various simulations suggest the structure is not an intrinsic structure of the simulated beams or sky model.

Overall, while this study is somewhat limited in scope, we find the results encouraging. Since HERA and other instruments seeking the EoR signal are interferometers, an important extension of this framework will be to examine the compressibility or even general effectiveness of this basis in the presence of mutual coupling between receiving elements. However, results from this work should be directly applicable to global 21-cm signal experiments. Another important effect to study in future work is whether the basis is equally effective at characterizing physical perturbations to the receiving elements due to imperfect design, effects of weather conditions such as high wind, improper soil or ground plane modelling, etc. Despite that these types of errors tend to produce highly non-trivial structure in the beam pattern, the fact that the basic problem is tractable in this framework suggests a way forward for more complicated considerations. Most importantly, the application of this within a Bayesian inference framework will allow the observer to infer what types of perturbations their physical beam exhibits relative to their a priori model. This will then prove to be a powerful tool for enhancing the fidelity of precision measurements of 21-cm reionization signals.

ACKNOWLEDGEMENTS

We are grateful to Aman Chokshi for helpful discussions.

This result is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 948764; MJW, JB, PB, KAG). We acknowledge the Science and Technology Facilities Council (STFC) for their support of Eloy de Lera Acedo. We acknowledge use of the following software: matplotlib (Hunter 2007), numpy (van der Walt, Colbert & Varoquaux 2011), and scipy (Virtanen et al. 2020).

DATA AVAILABILITY

The computer code used to model the beam patterns in this paper is available from https://github.com/HydraRadio/Hydra.

Footnotes

REFERENCES