© ESO 2021

1. Introduction

The technique of very long baseline interferometry (VLBI) provides the highest angular resolutions in the study of the Universe. With this technique, it is possible to synthesize an aperture similar to the diameter of the Earth by combining the signals from radio telescopes distributed across the planet surface (e.g., Thompson et al. 2017). This technique is based on the Van Cittert-Zernike theorem, which is a direct Fourier relationship between the sky brightness distribution and the spatial coherence of the electromagnetic field that reaches the Earth.

The heterodyne receivers of the VLBI radio telescopes are able to split the sky signal into two orthogonal polarization channels, which are usually right-handed and left-handed circular polarizations (RCP and LCP, respectively); these channels are recorded independently. This splitting allows us to recover the full-polarization information of the observed sources by extending the van Cittert-Zernike theorem to the four Stokes parameters via the so-called radio interferometer measurement equation (RIME; Hamaker et al. 1996).

The RIME relates the complex visibilities (i.e., a measurement of the spatial coherence of the radio waves), which are computed between the polarization channels of two radio telescopes, to the Fourier transforms of the four Stokes parameters of the observed source. The RIME also accounts for different instrumental effects such as electronic gains, atmospheric effects, and various kinds of instrumental polarization.

The problem of the correct calibration of the instrumental polarization in astronomical interferometry has already been studied for a long time (e.g., Cotton 1993; Leppanen et al. 1995; Martí-Vidal et al. 2012). However, the most used polarimetry algorithms to date suffer from several limitations based on several assumptions that, even though well justified for the mainstream radioastronomical hardware that was used a decade ago, may not be suitable for the new-generation instrumentation. Assumptions, such as a constant instrumental polarization across the frequency bandwidth of a receiver or the use of a linearized model to describe the leakage between antenna polarizers, are too simplistic for the calibration of (ultra-)wideband feeds and high dynamic-range observations.

In this paper, we discuss the main limitations of the classical algorithms used for the calibration of the instrumental polarization in VLBI and present new software, based on the RIME, which overcomes such limitations in different ways. In the next section, we introduce the basic aspects of VLBI polarimetry. In Sect. 2.1, we describe the polarimetry algorithm used by the majority of the VLBI community and emphasize its limitations. In Sects. 2.3 and 3, we describe our software and show realistic simulations of VLBI observations, where the advantages of our algorithm can be quantitatively assessed. In Sects. 4.1 and 4.2, we apply our algorithm to real Global mm-VLBI Array (GMVA) and Very Long Baseline Array (VLBA) full-polarization observations of sources 3C 345 and 3C 279, previously reported in Martí-Vidal et al. (2012, hereafter MV12) and Kravchenko et al. (2020, hereafter K20). In Sect. 5, we summarize our conclusions.

2. VLBI polarimetry modeling

In the RIME framework (Hamaker et al. 1996; Smirnov 2011), the complex visibilities between two VLBI antennas, a and b, are formulated as a visibility matrix, V^ab, of the form

$$\begin{array}{c}\hfill {\mathbf{V}}^{\mathit{ab}}=\left[\begin{array}{cc}{V}_{\mathit{rr}}^{\mathit{ab}}& V_{\mathit{rl}}^{\mathit{ab}}\\ V_{\mathit{lr}}^{\mathit{ab}}& V_{\mathit{ll}}^{\mathit{ab}}\end{array}\right],\end{array}$$(1)

where r and l refer to the RCP and LCP polarization channels. The measured visibilities are a combination of two independent contributions: on the one hand, the instrumental plus atmospheric effects and, on the other hand, the brightness distribution of the four Stokes parameters of the observed source. Regarding the instrumental contribution to the signal polarization, there are several quantities involved. These range from the relative gains between polarizers (i.e., differences in phases, delays, and amplitudes between RCP and LCP) to the signal leakage between polarizers; that is, the amount of RCP signal that is being injected into the LCP channel, and vice versa, with a given amplitude, delay and phase. All these instrumental and atmospheric quantities can be arranged in the form of one single Jones matrix. For antenna a, this matrix is equal to (Smirnov 2011)

$$\begin{array}{c}\hfill {\mathbf{J}}^a={\mathbf{G}}^a{\mathbf{D}}^a{\mathbf{P}}^a=\left[\begin{array}{cc}{G}_r^a& 0\\ 0& G_l^a\end{array}\right]\left[\begin{array}{cc}1& D_r^a\\ D_l^a& 1\end{array}\right]\left[\begin{array}{cc}{e}^{j{\varphi }_a}& 0\\ 0& e^{-j{\varphi }_a}\end{array}\right],\end{array}$$(2)

where j is the imaginary unit, the quantities $D_r^a$ and $D_l^a$ (called D-terms) are used to model the signal leakage between the polarization channels; $G_r^a$ and $G_l^a$ account for the electronic gains (plus the atmospheric delay/opacity) at each polarizer, and ϕ_a is the angle of the feed orientation of antenna a when pointing to the observed source. The matrix P^a accounts for the effect of the feed angle in the source frame, G^a is the gain matrix, and D^a is the D-term matrix. The D-terms may depend on frequency, but are assumed to be stable in time during the extent of the observations.

Regarding the contribution of the observed source to the RIME, the equation relates the complex visibilities between antennas a and b to the so-called brightness matrix, B^ab, which contains the brightness distributions of the four Stokes parameters (namely, ℐ, 𝒬, 𝒰, and 𝒱) in the Fourier domain, via the matrix equation

$$\begin{array}{c}\hfill {({\mathbf{J}}^a)}^{-1}{\mathbf{V}}^{\mathit{ab}}{({({\mathbf{J}}^b)}^H)}^{-1}={\mathbf{B}}^{\mathit{ab}}=\left[\begin{array}{cc}\stackrel{\sim }{\mathcal{I}}+\stackrel{\sim }{\mathcal{V}}& \stackrel{\sim }{\mathcal{Q}}+j\stackrel{\sim }{\mathcal{U}}\\ \stackrel{\sim }{\mathcal{Q}}-j\stackrel{\sim }{\mathcal{U}}& \stackrel{\sim }{\mathcal{I}}-\stackrel{\sim }{\mathcal{V}}\end{array}\right].\end{array}$$(3)

The Fourier transforms ($\stackrel{\sim }{\mathcal{I}}$, $\stackrel{\sim }{\mathcal{Q}}$, $\stackrel{\sim }{\mathcal{U}}$, $\stackrel{\sim }{\mathcal{V}}$) of the Stokes parameters are evaluated at the point in Fourier space given by the baseline between antennas a and b, projected onto the plane orthogonal to the source direction (Thompson et al. 2017). All the instrumental effects are encoded in the Jones matrices, J^a and (J^b)^H, where H is the Hermitian operator. We assume that direction-dependent instrumental effects are negligible, which holds for the small fields of view achievable with VLBI.

If the polarized structure of the observed source cannot be resolved by the synthesized resolution of the interferometer, the solution to Eq. (3) is relatively simple, given that the Fourier transforms of all Stokes parameters are constant at all the points in Fourier space being sampled by the instrument. This special case allows us to model the source contribution to Eq. (3) using one real number per Stokes parameter, and the calibration problem becomes trivial. However, the high spatial resolutions achieved in VLBI, especially at high frequencies, imply that almost all the observable sources (even the most distant and compact calibrators) have some resolvable structures with nonconstant and complex-valued Fourier transforms.

Designing an accurate way to model the structure contribution of a spatially resolved polarization calibrator in Eq. (3) is the key to any VLBI polarization calibration algorithm. The different approaches to model the contribution from the source polarization can be divided into three main families. The first group is the simple parametrization of the polarization source structure, (e.g., describing the brightness distributions via parametric models, as in Martí-Vidal et al. 2014). The second example is the use of (nonparametric) deconvolution techniques in full polarization with an iterative fitting of the D-terms, as in the new versions of the EHT imaging library (Chael et al. 2018). Reconstructing the source polarization iteratively, using the incrementally calibrated visibilities, and then fixing the source to estimate the instrumental polarization, is another approach to these calibration algorithms. This strategy was first proposed by Cotton (1993).

The third family of algorithms refer to a complete exploration of the whole parameter space using methods such as Markov chain Monte Carlo (MCMC), as in the new versions of the EHT Themis library (Broderick et al. 2020) and the Dterm modeling code, DMC (Pesce, in prep.). These MCMC methods are able to recover the posterior distributions of the D-terms and their impact on each part of the reconstructed images, but their main disadvantage is the computational costs, especially for cases in which the dimension of the parameter space makes any Monte Carlo exploration impractical (for example, the case of multi-calibrator geodesy observations, as discussed in Sect. 3.1).

2.1. LPCAL

The LPCAL program (Leppanen et al. 1995), provided as a task in the Astronomical Image Processing System (AIPS) of the National Radio Astronomy Observatory (NRAO)¹ is the best known polarimetry algorithm and that most used by the VLBI community. It has been tested and exploited for decades, to the point that the procedures followed by this algorithm can be considered as the standard for VLBI interferometric polarimetry. The LPCAL algorithm models the source polarization using a parametric approach. The total-intensity brightness distribution, ℐ, is divided in a set of disjoint regions (or polarization subcomponents), ℐ_i, such that

$$\begin{array}{c}\hfill \mathcal{I}={\displaystyle \sum _i^N{\mathcal{I}}_i.}\end{array}$$(4)

Each of these subcomponents is assumed to have a constant fractional polarization. Hence, the brightness distributions of Q and U are written as

$$\begin{array}{c}\hfill \mathcal{Q}={\displaystyle \sum _i^N{q}_i{\mathcal{I}}_i;\mathcal{U}=\sum _i^N{u}_i{\mathcal{I}}_i,}\end{array}$$(5)

where all q_i and u_i are real-valued parameters defined in the interval $q_i^2+u_i^2\in [0,1]$. We note, however, that the LPCAL algorithm imposes no bounds to these parameters. The quantity 𝒱 is given by the difference between the calibrated parallel-hand correlations (i.e., the diagonal elements of the brightness matrix, Eq. (3)), whereas the D-terms may dominate the cross-hand correlations. Therefore, 𝒱 can be neglected in the calibration of the D-terms.

The LPCAL approach allows us to model relatively complex brightness distributions using a minimum set of parameters for the source polarization. In the LPCAL procedure, the complex source structure in total intensity is first deconvolved using, for example, the CLEAN algorithm (Högbom 1974), and such a structure is then divided into regions of assumed constant polarization; that is, the algorithm assumes that the fractional polarization and electric vector position angle (EVPA) of the source is constant across the extent of each polarization subcomponent (ℐ_i in Eq. (4)). This is an enhanced version of the so-called similarity assumption (Cotton 1993). Depending on how the source division is done, the LPCAL source model may deviate from reality in cases of sources with complex polarization structures. For instance, sources with a strongly polarized emission that may be shifted with respect to the peak of ℐ, as regularly happens in AGN jets (e.g., O’Sullivan et al. 2011; Martí-Vidal et al. 2012; Molina et al. 2014; Gómez et al. 2016).

The LPCAL algorithm estimates the parameters that define the source polarization (i.e., q_i and u_i in Eq. (5)) and the D-terms of all antennas simultaneously by means of a least-squares minimization of the error function given by (Leppanen et al. 1995)

$$\begin{array}{c}\hfill {\chi }^2={\chi }_{\mathit{rl}}^2+{\chi }_{\mathit{lr}}^2,\end{array}$$(6)

with

$$\begin{array}{c}\hfill {\chi }_{\mathit{rl}}^2={\displaystyle \sum _k{|\frac{V_{\mathit{rl}}^k}{G_r^a{(G^b)}_l^{\ast }}-\frac{D_r^a{V}_{\mathit{ll}}^k}{G_l^a{(G^b)}_l^{\ast }}-\frac{{(D^b)}_l^{\ast }{V}_{\mathit{rr}}^k}{G_r^a{(G^b)}_r^{\ast }}-{\mathcal{P}}_{1,k}|}^2{w}_k}\end{array}$$(7)

and

$$\begin{array}{c}\hfill {\chi }_{\mathit{lr}}^2={\displaystyle \sum _k{|\frac{V_{\mathit{lr}}^k}{G_l^a{(G^b)}_r^{\ast }}-\frac{D_l^a{V}_{\mathit{rr}}^k}{G_r^a{(G^b)}_r^{\ast }}-\frac{{(D^b)}_r^{\ast }{V}_{\mathit{ll}}^k}{G_l^a{(G^b)}_l^{\ast }}-{\mathcal{P}}_{2,k}|}^2{w}_k,}\end{array}$$(8)

where ω_k is the weight of the kth visibility. In these equations, ${\mathcal{P}}_{1,k}=({\stackrel{\sim }{\mathcal{Q}}}_k+j{\stackrel{\sim }{\mathcal{U}}}_k)e^{j({\varphi }_a+{\varphi }_b)}$ and ${\mathcal{P}}_{2,k}=({\stackrel{\sim }{\mathcal{Q}}}_k-j{\stackrel{\sim }{\mathcal{U}}}_k)e^{-j({\varphi }_a+{\varphi }_b)}$, where the Fourier transforms of the Stokes parameters, ${\stackrel{\sim }{\mathcal{Q}}}_k$ and ${\stackrel{\sim }{\mathcal{U}}}_k$, are evaluated at the point of uv-space corresponding to the kth visibility. We note that the antenna indices, a and b, depend implicitly on the visibility index, k, as well as the complex gains, G^a and G^b (for both polarizers, r and l), which must be known for all the observing times of the experiment. The antenna gains can be estimated by means of self-calibration to Stokes ℐ (Pearson & Readhead 1984). In these equations, we omitted k in the gains and the antenna indices for clarity.

Thanks to the parallactic-angle coverage in a set of VLBI observations (i.e., the time changes in ϕ_a and ϕ_b), it is possible to decouple the D-terms from the source polarization parameters in Eqs. (7) and (8). Essentially, the more complete the ϕ coverage in a VLBI experiment, the more accurate the D-term calibration.

2.2. Limitations of LPCAL

Equations (6)–(8) are an approximation to the exact equation

$$\begin{array}{c}\hfill {\chi }_{\mathrm{RIME}}^2={\displaystyle \sum _{k,m}{w}_k\left({|{\left[{({\mathbf{J}}^a)}^{-1}{\mathbf{V}}^{\mathit{ab}}{({({\mathbf{J}}^b)}^{-1})}^H\right]}_m-{\mathbf{B}}_m|}^2\right),}\end{array}$$(9)

where m runs over the rl and lr matrix elements of Eq. (3). This exact RIME version of the polarimetry error function contains several elements that are proportional to products among the different antenna D-terms. However, the LPCAL model (Eq. (6)) only contains elements that depend linearly on the D-terms. As long as the D-terms (the amount of signal leakage between antenna polarizers) are small enough and the signal is limited by sensitivity (i.e., the dynamic range of the polarization is low), the linear approximation of LPCAL should be good enough for the polarization calibration. However, modern VLBI interferometers may require the exact RIME equation, as we discuss in Sect. 3.2.3.

The LPCAL algorithm also suffers from other limitations that, in special situations, can be critical. For example, the D-terms of all the antennas appearing in the data always have to be solved, with no restrictions nor use of any a priori information about their polarization leakage. If, for some reason, the D-terms of a given antenna are known (e.g., they are provided by the corresponding observatory), that information would not be usable by LPCAL. Either all the data related to that antenna would not be used by LPCAL or the D-terms of that antenna would have to be fitted, together with those of the remaining antennas, with no restrictions on their values. Furthermore, if there is an antenna with a missing polarization channel, that is, only the data of one circular polarizer have been properly recorded, LPCAL cannot estimate the D-terms for that antenna. The reason for this limitation can be seen in Eqs. (7) and (8), where the observed parallel-hand visibilities of both polarizations, V_rr and V_ll, are explicitly needed for the construction of the model of both, V_rl and V_lr.

In addition to these issues, LPCAL can only use one calibrator at a time, which may have important limitations in some VLBI experiments in which some of the antennas may have a very limited parallactic-angle coverage for the source being used as calibrator. Last but not least, there is no way to parameterize any frequency dependence of either the instrumental or the source-intrinsic polarization, which is a clear limiting factor for the analysis of wideband VLBI observations.

2.3. PolSolve

We developed a new code, polsolve, which overcomes the limitations of LPCAL described in the previous section. Our code is written in C++ and compiled as a Python module for use in the NRAO’s Common Astronomy Software Applications (CASA) software infrastructure², to which VLBI capabilities have recently been added (van Bemmel et al. 2018; Janssen et al. 2019). Our algorithm is part of a VLBI polarimetry toolbox that we developed for CASA and is publicly available with a GPLv3 license³. The polsolve algorithm estimates the D-terms and the source polarization parameters from the least-squares minimization of Eq. (9), which is the error function computed with the exact RIME. Hence, polsolve is not limited to the case of small polarization leakage at the antennas and/or low fractional source polarizations. Our algorithm is also able to use many calibrator sources (each one with its own polarization brightness modeling) in a simultaneous fit, which may help optimize the parallactic-angle coverage limitations in a VLBI array.

Regarding the instrumental and source-intrinsic frequency dependence of wideband observations, polsolve offers different possibilities. On the one hand, the D-terms can be expanded as a Taylor polynomial (of any degree) in frequency space. On the other hand, we implemented the so-called multi-IF fitting mode, where the D-terms are allowed to take different values at each frequency channel, but the source polarization is forced to be consistent across the full bandwidth. By consistent, we mean that the source polarization is modeled using the same parameter values for all sub-bands with the option of adding frequency depolarization and/or Faraday rotation (see Sect. 3.1 for an example). The advantage of the multi-IF approach is that any frequency dependence of the D-terms that cannot be parameterized with a smooth polynomial of low degree (e.g., oscillatory trends and jumping close to the baseband boundaries) can still be properly described with a minimum number of fitting parameters. A detailed example of the multi-IF fitting is given in Sect. 3.

In the case of antennas with missing polarization channels, polsolve uses the model prediction of the parallel-hand visibilities (obtained from the source deconvolution), so that the polarization leakage of the single observed polarization channel can still be determined.

The structure of the source polarization can be parameterized in the same way as LPCAL (i.e., by partitioning the structure as in Eqs. (4) and (5), neglecting 𝒱), but polsolve also gives the possibility of using an iterative approach, where a full-polarization source model, obtained via such programs as CLEAN, executed on all Stokes parameters, can be fixed in the fit of the instrumental polarization. This strategy, which can be called polarimetry self-calibration, was initially proposed by Cotton (1993), although its use has been very limited in the VLBI community. This is likely because it is not directly implemented in the well-known (and most used) LPCAL program. Actually, a more recent program based on LPCAL, called GPCAL (Park 2021), has finally implemented this long-sought feature in AIPS. The polarimetry self-calibration iterations can be easily implemented in polsolve via the use of helper functions that are provided with the code. In Sect. 3.2.2, we present a case in which the performance of the LPCAL source parametrization is compared to the polarimetry self-calibration.

3. Synthetic data tests

Each of the limitations described in Sect. 2.2 may have an impact on observations performed with modern VLBI systems. The wide fractional bandwidths introduce chromatic effects in the instrumental polarization, given that the systems are designed to have an optimum polarization response at a given nominal frequency in each band; departures from that frequency thus introduce unavoidable leakage effects.

We performed a quantitative study of the effects that the limitations in LPCAL may have on the calibration of current VLBI observations. For this study, we simulated several sets of VLBI observations using our CASA polarization toolbox. We studied the effects of wide fractional bandwidths and limited parallactic-angle coverage by generating synthetic VLBI Global Observing System (VGOS) observations performed in the geodetic mode. The nonlinear effects of the instrumental polarization and the limitations related to the LPCAL source parametrization were studied by simulating a VLBI experiment at 230 GHz with the Event Horizon Telescope (EHT) array. In this case, we used high leakage factors and a simulated source with a high fractional polarization. In the next subsections, we describe these simulations in more detail. In the package provided with the polsolve code, we included all the simulation and analysis scripts needed to reproduce the results shown in the following subsections.

3.1. Ultra-wide frequency band: VGOS simulations

A worldwide new-generation interferometer, VGOS, was developed by the International VLBI Service (IVS), aimed at high-precision (up to submillimeter) geodesy observations. To achieve the signal-to-noise ratio (S/N) required for such an accuracy, the VGOS antennas are equipped with ultra-wide receivers, covering from 2 GHz to 15 GHz (e.g., Petrachenko et al. 2012). The VGOS bandwidth makes it necessary to observe in a basis of linear polarization. Such a polarization basis prevents a proper phase calibration of the resulting VLBI observations, unless the visibility matrices are converted to a circular (i.e., RCP and LCP) basis (Martí-Vidal et al. 2016). The polconvert algorithm has already been successfully applied for the conversion of VGOS data from linear to circular basis (Alef et al. 2019).

The estimates of polconvert for the cross-polarization antenna gains are based on the “cross-polarization global fringe fitting” technique (Eq. (16) of Martí-Vidal et al. 2016), which has limitations for cases of poor parallactic-angle coverage and/or calibrator observations with a low S/N. As a consequence, any small deviation in the estimate of the phase difference, Δα_a, between the linear polarizers of any antenna, a, and the polconvert estimate is converted into a polarization leakage term in the circular-basis visibility matrices (Matthews et al. 2018; Goddi et al. 2019). It can be shown (Goddi et al. 2019) that the D-terms corresponding to the polarization conversion with a Δα_a phase offset are pure imaginary and equal for RCP and LCP. Such a leakage corruption in the VGOS polarimetry can be corrected with the use of polsolve, as we discuss in this section.

We used the polsimulate CASA task, which is part of our CASA toolbox for VLBI polarimetry, to simulate a set of VGOS observations with a large fraction of the C-band frequency coverage of the system (from 4 GHz to 6 GHz). Our frequency setup consists of 64 spectral windows distributed homogeneously across the simulated bandwidth. We introduced a chromatic (frequency-dependent) phase offset between the polarizers of all the antennas to simulate the residual post-conversion leakage after a suboptimal run of polconvert. The simulated frequency dependence of the phase offsets is sinusoidal, with an amplitude of 10 degrees (which would introduce a maximum polarization leakage on the order of 10%; Goddi et al. 2019) and random phases and periods. The simulated VGOS array consists of the triplet of European sites currently operative (Alef et al. 2019). These are located at the Onsala Space Observatory (Sweden), the Wettzel Observatory (Germany), and the Yebes Observatory (Spain).

The main goal of the VGOS project is to perform a quasi-continuous monitoring of VLBI geodetic parameters via the observation of ICRF3-related active galactic nuclei (AGN; about 300 defining sources and more than 4500 supplementary sources; Fey et al. 2015). Therefore, our simulations consist of a simplified version of a geodesy-like experiment, with observations of ten ICRF3 sources, organized in snapshot scans (30 seconds long) distributed along a total of two hours. Each source is only observed in two scans during the experiment, with a typical time separation of half an hour. It is clear that the parallactic-angle coverage of each source, which is what allows us to perform a robust calibration of the D-terms, is minimal in this simulation.

The simulated source structures consist of polarized point sources with random fractional polarizations, p/I, (with p/I ranging between 1% and 10%), random EVPAs, χ, random spectral indices, α (with absolute values between 0 and 1), and random rotation measures, RM, introduced by external Faraday screens with RM ranging, in absolute value, between 0 and 100 rad m⁻². The flux densities of all sources are fixed to 1 Jy (at the fiducial frequency of 4 GHz). The simulated properties of the VGOS sources are given in Table 1. The system temperatures of the antennas are set to 50 K across the whole bandwidth, which results in a S/N of ∼0.6−6 for the cross-hand correlations, V_rl and V_lr, at each spectral window and scan.

Obviously, the S/N of the fringes is too low for an independent fit of the D-terms to each spectral window and source; this kind of fit would be performed with the classical LPCAL software. However, we used polsolve to fit the whole dataset (i.e., all sources and spectral windows in one simultaneous fit) by imposing consistent source polarization along the whole observing bandwidth (i.e., the same polarization parameters, including frequency de-polarization and/or Faraday rotation, across the band). This way, the effective S/N used in the fit corresponds to the whole bandwidth (i.e., S/N ∼ 50 per scan). For each source, we modeled the Stokes parameters at a fiducial frequency equal to the lowest frequency in the data, together with a RM (i.e., $RM=\frac{\partial \chi }{\partial {\lambda }^2}$) and a polarization spectral index, α (i.e., p ∝ ν^α) that accounts for frequency (de)polarization. These quantities allowed us to model the source polarization in the whole observing band using only four parameters. Regarding the D-terms, we allowed for independent values at each frequency channel; we used the multi-IF fitting approach mentioned in Sect. 2.3.

The result of the global fit with polsolve is shown in Fig. 1. We note that the model to the whole dataset consists of 768 parameters for the instrumental polarization (i.e., two complex numbers for each antenna and frequency channel) plus 40 parameters for the source polarization; these parameters include four quantities per source, namely, two Stokes parameters (Q and U), a RM, and a polarization spectral index. This number of parameters is much lower than the total of 7700 parameters that would be needed by LPCAL (i.e., one independent fit for each source and spectral channel). In the left panels of Fig. 1, we show correlation plots between the true source polarization quantities and the values fitted with polsolve. The error bars shown in Fig. 1 are derived by polsolve from the postfit covariance matrix, by assuming a reduced χ² equal to its expected value (similar to the strategy followed in Martí-Vidal et al. 2014). In the figure panels at right, we show the simulated D-terms of all three antennas, together with the values fitted with polsolve. It is remarkable that the program is able to recover the source polarization parameters and the instrumental polarization using only two scans per source (and with such a low S/N per spectral window). The use of several calibrators in one simultaneous fit allowed us to maximize the information in the fit, in such a way that, even with only two scans per source, were able to decouple the effects of the instrumental polarization from the source contribution in the cross-hand visibilities.

Fig. 1. Results of the polsolve fit to the simulated VGOS observations. Left panels: correlation plots between the true source polarization quantities (fractional polarization and EVPA at the fiducial frequency, RM, and polarization spectral index) and the polsolve estimates. The dotted lines indicate the 1:1 correlation. Right panels: simulated D-term spectra as thick lines, resulting from a frequency-dependent phase offset between the linear polarizers of the antennas, and the values fitted by polsolve.

We note that in the cases of wider bandwidths, which could even cover different radio bands, any contribution of internal Faraday rotation in the calibrators may cause their EVPAs to depart from the canonical λ² law, hence biasing the polsolve model. In addition, opacity effects in the jets may introduce position shifts of the brightness peaks as a function of frequency, which would bias the D-term estimates if not taken into account. In these cases, the use of polarization self-calibration (see Sect. 3.2.2), based either on a (full-polarization) multifrequency-synthesis deconvolution (Rau & Cornwell 2011) or on the construction of an image “cube” (i.e., different model images per band), would be an alternative way to calibrate the VGOS instrumental polarization with polsolve, which would still allow for the use of the multi-IF mode.

3.2. Complex source structures: EHT simulations

The EHT is an Earth-size millimeter VLBI array, whose main goal is the imaging of horizon-scale structures in close by supermassive black holes (Event Horizon Telescope Collaboration 2019a). In this paper, we simulated EHT observations following the array configuration of the 2017 campaign, consisting of eight observatories at six geographical locations (Event Horizon Telescope Collaboration 2019b): Atacama Large Millimeter/sub-millimeter Array (ALMA) and the Atacama Pathfinder Experiment (APEX) in the Atacama Desert in Chile; the Large Millimeter Telescope Alfonso Serrano (LM) on the Volcán Sierra Negra in Mexico; the Arizona Radio Observatory (AZ), situated on Kitt Peak, Arizona; the South Pole Telescope (SPT) at the Geographic South Pole; the James Clerk Maxwell Telescope (JCMT) and Sub-millimeter Array (SMA) on Maunakea in Hawai’i (USA); and the IRAM 30 m telescope, on Pico Veleta in Spain. Table 2 lists their respective antenna mounts, which affect the feed angles and, as a consequence, the polarimetry calibration. All these antennas, with the exception of ALMA, used circular polarization feeds. The visibilities related to ALMA were converted into a pure circular basis using the polconvert program (Martí-Vidal et al. 2016; Matthews et al. 2018), which may introduce a small pure-imaginary and symmetric instrumental polarization at ALMA (e.g., Goddi et al. 2019).

3.2.1. Generating synthetic EHT data

We used the polsimulate task to generate synthetic EHT observations of a polarized source with a complex structure to test the limitations of the LPCAL algorithm for the parametrization of the source polarization brightness distribution (i.e., the source partition approach formulated in Eqs. (4) and (5)) and compare its performance to the iterative polarimetry self-calibration approach.

The whole set of antennas that participated in the EHT 2017 observing campaign were included in the simulations. There is only one source in these synthetic data, whose coordinates are set to those of the Galactic center (Sagittarius A^*). The exact observing times are set to those of the real EHT observations of the Galactic center taken on April 7, 2017. The observing frequency is set to 230 GHz, with a bandwidth of 2 GHz, and the system temperatures, T_sys, of all antennas are set to a single value of 50 K. The EHT is a heterogeneous array, with different ranges of T_sys at each station, and a T_sys of 50 K is a rather low value for a receiver at 230 GHz (only the phased ALMA approaches that value; see, e.g., Table 2 of Event Horizon Telescope Collaboration 2019c). However, we note that real calibrator observations with the EHT (e.g., 3C 279; Kim et al. 2020) are currently limited by dynamic range, that is, sensitivity may not be the main limitation in the calibration of the instrumental polarization of the EHT.

The visibility noise is estimated by polsimulate via the CASA simulation tool, sm, based on the effective diameters of the telescopes (Event Horizon Telescope Collaboration 2019b) and the T_sys. The real and imaginary parts of the D-terms for all the antennas are taken from a random Gaussian distribution, centered at zero with a standard deviation of 0.2. We note that the D-terms resulting from this distribution are usually large, compared to the typical D-terms of EHT antennas (e.g., Johnson et al. 2015) and are not representative of the true instrumental polarization of the EHT.

The simulated source structure has a core-jet shape oriented in the east-west direction. There are three unpolarized point-like components (a 1 Jy core plus two 0.8 Jy jet components, separated by 40 and 80 μas from the core) and two polarized components (at 10 and 30 μas from the core, with 0.4 Jy each) with high fractional polarizations (0.94 and 0.90) and different EVPAs (16 deg and −45 deg). In Fig. 2, we show the CLEAN reconstruction (using uniform visibility weighting) of the simulated data after correcting for the effects of the instrumental polarization. We made the distribution of the source components in such a way that the peaks of polarized intensity are shifted with respect to the peaks in total intensity. These shifts in polarization brightness are a problem for the LPCAL source parametrization, since the polarization structure is not proportional to that of Stokes ℐ, either across the extent of the core or on each of the jet components. Therefore, the D-term estimates resulting from the LPCAL approach may be biased by the contribution of the source structure that, owing to the shifts in the polarization intensity, cannot be reproduced by the source model (Eq. (5)).

Fig. 2. CLEAN image of the synthetic EHT observations free from instrumental polarization. The full width at half maximum (FWHM) of the restoring beam is shown at bottom left. The 16 contour levels are separated in logarithmic scale and represent total intensity (peak of 1.28 Jy beam−1; lowest contour of 1.3 mJy beam−1); the blue color scale represents polarized intensity (peak of 0.36 Jy beam−1) and the red lines represent the EVPA spatial distribution. The yellow crosses indicate the locations of the polarization intensity peaks, which are shifted with respect to those of Stokes ℐ. The green boxes indicate the extent of the polarization subcomponents used in the D-term fitting with the LPCAL approach.

3.2.2. Polarimetry self-calibration

For the D-term fitting using the LPCAL approach, we divided the source structure into three polarization subcomponents. These subcomponents are shown as green rectangles in Fig. 2. Each rectangle is centered on each of the main total-intensity peaks of the source. The results of the D-term fitting using this source parametrization are shown in Fig. 3 (left panel). Even though there is some correlation between the true and the estimated D-terms, the deviations in the correlation (which, in some cases, are as high as 20%) may introduce high residual instrumental effects in the cross-hand visibilities, which could limit the dynamic range (and fidelity) of the polarization image. The quality of the fit could be improved by increasing the number of polarization subcomponents (i.e., dividing the green rectangles of Fig. 2 into smaller pieces), but the final result would still depend on the exact location and extent of the subcomponents with respect to the location of the true polarization features of the source. In any case, the simple source subdivision used in this example reflects the limitations of the LPCAL strategy in a clear way.

Fig. 3. Results of polsolve applied to our simulated EHT observations. Left: correlation between the true antenna D-terms and the values fitted using the LPCAL similarity assumption (see text). Center: correlation for the D-terms fitted after fifty iterations of polarimetry self-calibration. Right: L1 norm between true and fitted D-terms, as a function of the number of self-calibration iterations, for the case of linear D-term approximation (black points) and full nonlinear D-term model (white points).

We can quantify the quality of the D-term fitting using the L₁ norm of the D-term residuals (the difference between true and estimated D-terms), that is,

$$\begin{array}{c}\hfill L_1={\displaystyle \sum _{\mathrm{ants}}|D_{\mathrm{fit}}-D_{\mathrm{true}}|.}\end{array}$$(10)

Using the D-term estimates from the LPCAL approach on our EHT simulated data (Fig. 3, left), the L₁ norm obtained is L₁ = 1.87.

In contrast to these results, in Fig. 3 we show the correlation between true D-terms and the polsolveD-term estimates obtained after 50 iterations of the polarimetry self-calibration approach (using the CLEAN algorithm as implemented in CASA). The correlation is, in this case, remarkably better than that obtained with the LPCAL calibration approach. We can reach values as low as L₁ = 0.2 (i.e., an order of magnitude lower than with the LPCAL approach) after 20–30 iterations.

3.2.3. Linear versus nonlinear leakage determinations

The bias in the D-term estimates produced by the LPCAL linear approximation of Eq. (9) may become important for cases of antennas with relatively high polarization leakage (D-terms higher than 10–20%) and/or sources with high fractional polarizations in Fourier space. This quantity, in principle, is unbounded, because of the possible presence of low visibility values, or even visibility nulls, in the parallel-hand correlations (Johnson et al. 2015).

The effects of the nonlinear contribution of the D-terms in the cross-hand visibilities of a VLBI dataset can be quantified with polsolve. It is possible to disconnect the nonlinear D-term contributions from the fit (via a special keyword in the CASA task), which allows us to compare directly the L₁ norm of the D-term residuals with and without nonlinear corrections. In Fig. 3, we show the L₁ norm of the D-term residuals (true values minus fitted values) as a function of the number of polarization self-calibration iterations for the fitting with linear approximation of the D-term effects and for the fitting with the full polsolveD-term model.

For a small number of iterations, both models give D-term estimates of similar quality. However, as the number of iterations increases, it becomes clear that the limitations of the D-term linear approximation do not produce any improvement of the calibration of the cross-hand data beyond a given limit. That is to say, there is a maximum achievable dynamic range in the model image, which only appears in the case of the D-term linear approximation. On the other hand, if the full D-term model is used in the fit, the full-polarization image of the source keeps improving with the number of self-calibration iterations with no appreciable saturation in the first 50 self-calibration iterations; as a consequence, the estimated D-terms keep approaching the true simulated values. We note, however, that other factors (e.g., an imperfect phase self-calibration) may contribute to increase the L₁ norm in real observations.

A direct conclusion from these results is that the nonlinear D-term model of polsolve allows us to reach higher dynamic ranges in the polarization images obtained from high-sensitivity VLBI observations, assuming that the antenna gains are properly calibrated.

4. Real observations

In this section, we present the results of polsolve applied to real VLBI experiments, for which full-polarization images are available in the literature. Since the published images are obtained from visibilities calibrated with LPCAL, we can test the performance of polsolve by comparing them to the results obtained in CASA.

4.1. 3C 345 at 86 GHz

The GMVA is an international collaboration led by the Max-Planck-Institut für Radioastronomie (MPIfR)⁴, between NRAO and several European institutions, aimed at providing an infrastructure for routine VLBI observations at millimeter wavelengths. In MV12, a full-polarization calibration pipeline was presented for the GMVA. That pipeline was based on AIPS and made use of the LPCAL program for the determination of the antenna D-terms. One of the conclusions from MV12 is that, even though there were some inconsistencies in the D-terms determined from multiple sources observed within the same experiment, the final images of sources with a high linear-polarization brightness did not depend strongly on the actual combination of D-term values used in the calibration. A rule of thumb from the results presented in MV12 is that any polarization intensity higher than ∼50% of the polarization peak has a very weak dependence on the different GMVA D-term estimates.

We reprocessed the data reported in MV12, which corresponds to the GMVA campaign of May 2010. We refer to that publication for the details of these observations, such as participating antennas, frequency configuration, and scheduling issues. We imported the fringe-fitted visibilities of source 3C 345 (for which a full-polarization image is shown in Fig. 9 of MV12) into CASA format, to perform the hybrid imaging and polarimetry calibration. We reimaged and self-calibrated the visibilities using the CASA tasks tclean, gaincal, and applycal by combining the two parallel-hand correlations for the gain solutions (to not affect the instrumental polarization). We also used Briggs weighting (Briggs 1995) with a robustness parameter of 0.5, which offers a good trade-off between sensitivity and spatial resolution. Then, we executed a set of 20 iterations of polarimetry self-calibration (see Sect. 3.2.2), which resulted in a fast and robust convergence of the antenna D-terms; the estimated values did not change above numerical noise after roughly five iterations. The resulting image and D-terms are shown in Fig. 4. We applied a polarization cutoff of 50% of the peak, following Fig. 9 of MV12. Even though the synthesized CLEAN beam of our image is different from that of Fig. 9 in MV12, the main features of the source (in total intensity and in linear polarization) are similar between the two images. In particular, there is linear polarization at the north of the core (and a bit upstream of the jet) with a clear bending of the EVPA, from parallel to the jet to perpendicular. Then, at about 100 μas and 200 μas downstream from the core, another two polarized components appear with their EVPAs oriented parallel to the jet direction.

Fig. 4. Results of polsolve applied to real GMVA observations. Left: GMVA image of 3C 345 obtained after the polsolve calibration of the observations reported in MV12. The FWHM of convolving beam is shown at bottom left. Contours are logarithmically spaced between 35 mJy beam−1 and 1.22 Jy beam−1 (the peak intensity). Right: correlation between the LPCAL and polsolve estimated D-terms.

Regarding the D-terms, we executed LPCAL by dividing the core region of 3C 345 into two subcomponents: one that covers the northern side of the jet, where there is the maximum polarization brightness, and another that covers the southern part, where there is no detected polarization. We also set other subcomponents, which cover the two polarization features downstream from the core. We note that it is not possible to model the EVPA rotation in the core with such a source subdivision, so we expect to have a degradation in the quality of the D-term estimates with LPCAL. The D-term estimates are shown in Fig. 4 (right). Similar to the case discussed in Sect. 3.2.2 (see Fig. 3, left), the correlation between the D-terms estimated with the LPCAL approach and those from the polarization self-calibration is rather poor (Fig. 4). This is likely because of the limited accuracy of the former method to account for the source polarization structure. In any case, the resemblance between the image obtained with AIPS/LPCAL (shown in MV12) and the image obtained with CASA/polsolve (Fig. 4) indicates a satisfactory calibration of the GMVA data with polsolve.

4.2. 3C 279 at 43 GHz

K20 report on several VLBA full-polarization observations at wavelengths of 7 mm and 1.3 cm of the active nucleus in M87. We downloaded the data corresponding to their project BG251A from the NRAO archive⁵. This project was observed in May 2017 at 7 mm (i.e., 43 GHz) and included several scans of source 3C 279 used as a polarization calibrator. The frequency configuration was arranged in eight contiguous sub-bands of equal width, covering a total of 256 MHz (see K20 for more details on the observations).

We reprocessed the BG251A observations of source 3C 279 to compare the polsolveD-term calibration to that obtained with LPCAL (as used by K20). The amplitude calibration, global fringe-fitting, and cross-polarization delay corrections were done in AIPS using standard procedures. Then, the data were exported into CASA (version 5.7) for hybrid imaging (several iterations of phase and amplitude self-calibration) also using standard procedures and Briggs weighting with a robustness parameter of 0.5. After the imaging, we executed ten iterations of polarization self-calibration with polsolve, using the multi-IF mode. The resulting CLEAN image in mfs mode, shown in Fig. 5 (left), is remarkably similar to the results reported in K20 (see their Fig. A.2, bottom left panel). Regarding the D-term estimates (Fig. 5, right), the polsolve values are correlated with those reported in K20 (see their Table A.2), although the correlation is also rather poor. The worst outliers in the correlation happen to be the cases with the largest D-term dispersion among sub-bands (i.e., the cases with the largest error bars; see Table A.2 of K20), which correspond to the antennas at Hancock (RCP polarizer) and St. Croix (LCP polarizer).

Fig. 5. Results of polsolve applied to real VLBA observations. Left: VLBA image of 3C 279 obtained from the polsolve calibration of the observations reported in K20. The FWHM of the convolving beam is shown at bottom left. The contours are logarithmically spaced between 50 mJy beam−1 and 8.4 Jy beam−1 (the peak intensity). Right: correlation between the D-term amplitudes reported in K20 (see their Table A.2) and the amplitudes of the frequency-averaged multi-IF D-terms estimated with polsolve.

5. Summary and conclusions

New VLBI arrays, with high sensitivities and wide instantaneous bandwidths, as compared to the most updated arrays in recent decades, require the use of improved calibration algorithms to account for the effects of instrumental polarization. In the case of wide bandwidths and geodesy-like observations, where many weak sources are observed with a limited parallactic-angle coverage, an algorithm that is able to combine the wide bandwidths (to increase the solution S/N) and the use of several calibrators simultaneously to maximize the amount of information in the fit is required. In the case of sensitive arrays and/or interferometers with very high spatial resolutions, special algorithms that account for the source polarization structures and correct for the nonlinear effects of the instrumental polarization may also be needed.

We present a new set of algorithms, implemented for use in the CASA interface and fully available to the VLBI community, that allow the calibration of instrumental polarization in VLBI observations to be performed without the limitations of the classical calibration algorithms used so far. Our algorithms (a realistic simulator, polsimulate, flexible solver, polsolve, and several other helper functions) allow us to perform advanced polarization calibration free of the limiting assumptions of the standard calibration software. These mainly include observations with narrow fractional bandwidths, no use of a priori information, and limited sensitivity in the polarization signal. We show the results of our software applied to realistic VLBI simulations and compare its performance with that of the most used VLBI polarimetry algorithm, LPCAL. We also apply polsolve to real millimeter-wave VLBI observations of sources 3C 345 and 3C 279 and obtain full-polarization images similar to those already published (MV12, K20). The advantages of polsolve may be critical for the optimum analysis of polarization observations taken with new-generation VLBI arrays.

Acknowledgments

This work has been partially supported by the MICINN Research Project PID2019-108995GB-C22. I.M.V. and A.M. also thank the Generalitat Valenciana for funding, in the frame of the GenT Project CIDEGENT/2018/021. M.J. is supported by the ERC Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (Grant 610058).

References

All Tables

All Figures

Fig. 1. Results of the polsolve fit to the simulated VGOS observations. Left panels: correlation plots between the true source polarization quantities (fractional polarization and EVPA at the fiducial frequency, RM, and polarization spectral index) and the polsolve estimates. The dotted lines indicate the 1:1 correlation. Right panels: simulated D-term spectra as thick lines, resulting from a frequency-dependent phase offset between the linear polarizers of the antennas, and the values fitted by polsolve.

In the text

Fig. 2. CLEAN image of the synthetic EHT observations free from instrumental polarization. The full width at half maximum (FWHM) of the restoring beam is shown at bottom left. The 16 contour levels are separated in logarithmic scale and represent total intensity (peak of 1.28 Jy beam−1; lowest contour of 1.3 mJy beam−1); the blue color scale represents polarized intensity (peak of 0.36 Jy beam−1) and the red lines represent the EVPA spatial distribution. The yellow crosses indicate the locations of the polarization intensity peaks, which are shifted with respect to those of Stokes ℐ. The green boxes indicate the extent of the polarization subcomponents used in the D-term fitting with the LPCAL approach.

In the text

Fig. 3. Results of polsolve applied to our simulated EHT observations. Left: correlation between the true antenna D-terms and the values fitted using the LPCAL similarity assumption (see text). Center: correlation for the D-terms fitted after fifty iterations of polarimetry self-calibration. Right: L1 norm between true and fitted D-terms, as a function of the number of self-calibration iterations, for the case of linear D-term approximation (black points) and full nonlinear D-term model (white points).

In the text

	Fig. 4. Results of polsolve applied to real GMVA observations. Left: GMVA image of 3C 345 obtained after the polsolve calibration of the observations reported in MV12. The FWHM of convolving beam is shown at bottom left. Contours are logarithmically spaced between 35 mJy beam−1 and 1.22 Jy beam−1 (the peak intensity). Right: correlation between the LPCAL and polsolve estimated D-terms.
In the text

Fig. 5. Results of polsolve applied to real VLBA observations. Left: VLBA image of 3C 279 obtained from the polsolve calibration of the observations reported in K20. The FWHM of the convolving beam is shown at bottom left. The contours are logarithmically spaced between 50 mJy beam−1 and 8.4 Jy beam−1 (the peak intensity). Right: correlation between the D-term amplitudes reported in K20 (see their Table A.2) and the amplitudes of the frequency-averaged multi-IF D-terms estimated with polsolve.

In the text