Issue 
A&A
Volume 646, February 2021



Article Number  A52  
Number of page(s)  10  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202039527  
Published online  05 February 2021 
Polarization calibration techniques for the newgeneration VLBI
^{1}
Departament d’Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50, 46100 Burjassot, València, Spain
email: i.martividal@uv.es
^{2}
Observatori Astronòmic, Universitat de València, Parc Científic, C. Catedràtico José Beltrán 2, 46980 Paterna, València, Spain
^{3}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{4}
Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics (IMAPP), Radboud University, PO Box 9010 6500 GL Nijmegen, The Netherlands
^{5}
Centro Astronómico de Yebes (IGN), Apartado 148, 19180 Yebes, Spain
Received:
25
September
2020
Accepted:
7
December
2020
The calibration and analysis of polarization observations in very long baseline interferometry (VLBI) requires the use of specific algorithms that suffer from several limitations, closely related to assumptions in the data properties that may not hold in observations taken with newgeneration VLBI equipment. Currently, the instantaneous bandwidth achievable with VLBI backends can be as high as several gigahertz, covering several radio bands simultaneously. In addition, the sensitivity of VLBI observations with the most updated equipment may reach dynamic ranges of tens of thousands, both in total intensity and polarization. In this paper, we discuss the impact of the limitations of common VLBI polarimetry algorithms on narrowfield observations taken with modern VLBI arrays, from the VLBI Global Observing System to the Event Horizon Telescope, and present new software that overcomes these limitations. In particular, our software is able to perform a simultaneous fit of multiple calibrator sources, include nonlinear terms in the model of the instrumental polarization and use a selfcalibration approach for the estimate of the polarization leakage in the antenna receivers.
Key words: techniques: polarimetric / techniques: interferometric
© 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 CittertZernike 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 righthanded and lefthanded circular polarizations (RCP and LCP, respectively); these channels are recorded independently. This splitting allows us to recover the fullpolarization information of the observed sources by extending the van CittertZernike theorem to the four Stokes parameters via the socalled 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 newgeneration 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 dynamicrange 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 mmVLBI Array (GMVA) and Very Long Baseline Array (VLBA) fullpolarization 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
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)
where j is the imaginary unit, the quantities and (called Dterms) are used to model the signal leakage between the polarization channels; and 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 Dterm matrix. The Dterms 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 socalled brightness matrix, B^{ab}, which contains the brightness distributions of the four Stokes parameters (namely, ℐ, 𝒬, 𝒰, and 𝒱) in the Fourier domain, via the matrix equation
The Fourier transforms (, , , ) 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 directiondependent 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 complexvalued 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 Dterms, 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 Dterms 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 multicalibrator 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)^{1} 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 totalintensity brightness distribution, ℐ, is divided in a set of disjoint regions (or polarization subcomponents), ℐ_{i}, such that
Each of these subcomponents is assumed to have a constant fractional polarization. Hence, the brightness distributions of Q and U are written as
where all q_{i} and u_{i} are realvalued parameters defined in the interval . We note, however, that the LPCAL algorithm imposes no bounds to these parameters. The quantity 𝒱 is given by the difference between the calibrated parallelhand correlations (i.e., the diagonal elements of the brightness matrix, Eq. (3)), whereas the Dterms may dominate the crosshand correlations. Therefore, 𝒱 can be neglected in the calibration of the Dterms.
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 socalled 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 Dterms of all antennas simultaneously by means of a leastsquares minimization of the error function given by (Leppanen et al. 1995)
with
and
where ω_{k} is the weight of the kth visibility. In these equations, and , where the Fourier transforms of the Stokes parameters, and , are evaluated at the point of uvspace 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 selfcalibration to Stokes ℐ (Pearson & Readhead 1984). In these equations, we omitted k in the gains and the antenna indices for clarity.
Thanks to the parallacticangle coverage in a set of VLBI observations (i.e., the time changes in ϕ_{a} and ϕ_{b}), it is possible to decouple the Dterms from the source polarization parameters in Eqs. (7) and (8). Essentially, the more complete the ϕ coverage in a VLBI experiment, the more accurate the Dterm calibration.
2.2. Limitations of LPCAL
Equations (6)–(8) are an approximation to the exact equation
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 Dterms. However, the LPCAL model (Eq. (6)) only contains elements that depend linearly on the Dterms. As long as the Dterms (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 Dterms 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 Dterms 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 Dterms 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 Dterms for that antenna. The reason for this limitation can be seen in Eqs. (7) and (8), where the observed parallelhand 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 parallacticangle 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 sourceintrinsic 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^{2}, 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^{3}. The polsolve algorithm estimates the Dterms and the source polarization parameters from the leastsquares 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 parallacticangle coverage limitations in a VLBI array.
Regarding the instrumental and sourceintrinsic frequency dependence of wideband observations, polsolve offers different possibilities. On the one hand, the Dterms can be expanded as a Taylor polynomial (of any degree) in frequency space. On the other hand, we implemented the socalled multiIF fitting mode, where the Dterms 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 subbands with the option of adding frequency depolarization and/or Faraday rotation (see Sect. 3.1 for an example). The advantage of the multiIF approach is that any frequency dependence of the Dterms 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 multiIF fitting is given in Sect. 3.
In the case of antennas with missing polarization channels, polsolve uses the model prediction of the parallelhand 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 fullpolarization 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 selfcalibration, 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 wellknown (and most used) LPCAL program. Actually, a more recent program based on LPCAL, called GPCAL (Park 2021), has finally implemented this longsought feature in AIPS. The polarimetry selfcalibration 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 selfcalibration.
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 parallacticangle 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. Ultrawide frequency band: VGOS simulations
A worldwide newgeneration interferometer, VGOS, was developed by the International VLBI Service (IVS), aimed at highprecision (up to submillimeter) geodesy observations. To achieve the signaltonoise ratio (S/N) required for such an accuracy, the VGOS antennas are equipped with ultrawide 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 crosspolarization antenna gains are based on the “crosspolarization global fringe fitting” technique (Eq. (16) of MartíVidal et al. 2016), which has limitations for cases of poor parallacticangle 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 circularbasis visibility matrices (Matthews et al. 2018; Goddi et al. 2019). It can be shown (Goddi et al. 2019) that the Dterms 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 Cband 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 (frequencydependent) phase offset between the polarizers of all the antennas to simulate the residual postconversion 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 quasicontinuous monitoring of VLBI geodetic parameters via the observation of ICRF3related 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 geodesylike 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 parallacticangle coverage of each source, which is what allows us to perform a robust calibration of the Dterms, 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^{−2}. 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 crosshand correlations, V_{rl} and V_{lr}, at each spectral window and scan.
Source properties used in our VGOS simulation.
Obviously, the S/N of the fringes is too low for an independent fit of the Dterms 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 depolarization 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., ) 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 Dterms, we allowed for independent values at each frequency channel; we used the multiIF 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 χ^{2} 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 Dterms 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 crosshand 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 Dterm spectra as thick lines, resulting from a frequencydependent 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 λ^{2} 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 Dterm estimates if not taken into account. In these cases, the use of polarization selfcalibration (see Sect. 3.2.2), based either on a (fullpolarization) multifrequencysynthesis 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 multiIF mode.
3.2. Complex source structures: EHT simulations
The EHT is an Earthsize millimeter VLBI array, whose main goal is the imaging of horizonscale 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/submillimeter 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 Submillimeter 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 pureimaginary and symmetric instrumental polarization at ALMA (e.g., Goddi et al. 2019).
Names of the EHT telescopes that participated in the 2017 campaign and their respective antenna mounts.
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 selfcalibration 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 Dterms 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 Dterms resulting from this distribution are usually large, compared to the typical Dterms 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 corejet shape oriented in the eastwest direction. There are three unpolarized pointlike 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 Dterm 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 Dterm fitting with the LPCAL approach. 
3.2.2. Polarimetry selfcalibration
For the Dterm 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 totalintensity peaks of the source. The results of the Dterm fitting using this source parametrization are shown in Fig. 3 (left panel). Even though there is some correlation between the true and the estimated Dterms, the deviations in the correlation (which, in some cases, are as high as 20%) may introduce high residual instrumental effects in the crosshand 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 Dterms and the values fitted using the LPCAL similarity assumption (see text). Center: correlation for the Dterms fitted after fifty iterations of polarimetry selfcalibration. Right: L_{1} norm between true and fitted Dterms, as a function of the number of selfcalibration iterations, for the case of linear Dterm approximation (black points) and full nonlinear Dterm model (white points). 
We can quantify the quality of the Dterm fitting using the L_{1} norm of the Dterm residuals (the difference between true and estimated Dterms), that is,
Using the Dterm estimates from the LPCAL approach on our EHT simulated data (Fig. 3, left), the L_{1} norm obtained is L_{1} = 1.87.
In contrast to these results, in Fig. 3 we show the correlation between true Dterms and the polsolveDterm estimates obtained after 50 iterations of the polarimetry selfcalibration 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_{1} = 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 Dterm estimates produced by the LPCAL linear approximation of Eq. (9) may become important for cases of antennas with relatively high polarization leakage (Dterms 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 parallelhand correlations (Johnson et al. 2015).
The effects of the nonlinear contribution of the Dterms in the crosshand visibilities of a VLBI dataset can be quantified with polsolve. It is possible to disconnect the nonlinear Dterm contributions from the fit (via a special keyword in the CASA task), which allows us to compare directly the L_{1} norm of the Dterm residuals with and without nonlinear corrections. In Fig. 3, we show the L_{1} norm of the Dterm residuals (true values minus fitted values) as a function of the number of polarization selfcalibration iterations for the fitting with linear approximation of the Dterm effects and for the fitting with the full polsolveDterm model.
For a small number of iterations, both models give Dterm estimates of similar quality. However, as the number of iterations increases, it becomes clear that the limitations of the Dterm linear approximation do not produce any improvement of the calibration of the crosshand 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 Dterm linear approximation. On the other hand, if the full Dterm model is used in the fit, the fullpolarization image of the source keeps improving with the number of selfcalibration iterations with no appreciable saturation in the first 50 selfcalibration iterations; as a consequence, the estimated Dterms keep approaching the true simulated values. We note, however, that other factors (e.g., an imperfect phase selfcalibration) may contribute to increase the L_{1} norm in real observations.
A direct conclusion from these results is that the nonlinear Dterm model of polsolve allows us to reach higher dynamic ranges in the polarization images obtained from highsensitivity 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 fullpolarization 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 MaxPlanckInstitut für Radioastronomie (MPIfR)^{4}, between NRAO and several European institutions, aimed at providing an infrastructure for routine VLBI observations at millimeter wavelengths. In MV12, a fullpolarization 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 Dterms. One of the conclusions from MV12 is that, even though there were some inconsistencies in the Dterms determined from multiple sources observed within the same experiment, the final images of sources with a high linearpolarization brightness did not depend strongly on the actual combination of Dterm 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 Dterm 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 fringefitted visibilities of source 3C 345 (for which a fullpolarization image is shown in Fig. 9 of MV12) into CASA format, to perform the hybrid imaging and polarimetry calibration. We reimaged and selfcalibrated the visibilities using the CASA tasks tclean, gaincal, and applycal by combining the two parallelhand 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 tradeoff between sensitivity and spatial resolution. Then, we executed a set of 20 iterations of polarimetry selfcalibration (see Sect. 3.2.2), which resulted in a fast and robust convergence of the antenna Dterms; the estimated values did not change above numerical noise after roughly five iterations. The resulting image and Dterms 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 Dterms. 
Regarding the Dterms, 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 Dterm estimates with LPCAL. The Dterm 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 Dterms estimated with the LPCAL approach and those from the polarization selfcalibration 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 fullpolarization 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^{5}. 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 subbands 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 polsolveDterm calibration to that obtained with LPCAL (as used by K20). The amplitude calibration, global fringefitting, and crosspolarization 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 selfcalibration) also using standard procedures and Briggs weighting with a robustness parameter of 0.5. After the imaging, we executed ten iterations of polarization selfcalibration with polsolve, using the multiIF 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 Dterm 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 Dterm dispersion among subbands (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 Dterm amplitudes reported in K20 (see their Table A.2) and the amplitudes of the frequencyaveraged multiIF Dterms 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 geodesylike observations, where many weak sources are observed with a limited parallacticangle 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 millimeterwave VLBI observations of sources 3C 345 and 3C 279 and obtain fullpolarization images similar to those already published (MV12, K20). The advantages of polsolve may be critical for the optimum analysis of polarization observations taken with newgeneration VLBI arrays.
Acknowledgments
This work has been partially supported by the MICINN Research Project PID2019108995GBC22. 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
 Alef, W., Anderson, J. M., Bernhart, S., et al. 2019, Proceedings of the 24th European VLBI Group for Geodesy and Astrometry Working Meeting, 107 [Google Scholar]
 Briggs, D. S. 1995, PhD Thesis, New Mexico Institute of Mining and Technology [Google Scholar]
 Broderick, A. E., Gold, R., Karami, M., et al. 2020, ApJ, 897, 139 [Google Scholar]
 Chael, A., Bouman, K., Johnson, M., et al. 2018, https://doi.org/10.5281/zenodo.1173414 [Google Scholar]
 Cotton, W. D. 1993, AJ, 106, 1241 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) et al. 2019c, ApJ, 875, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Fey, A. L., Gordon, D., Jacobs, C. S., et al. 2015, AJ, 150, 58 [Google Scholar]
 Goddi, C., MartíVidal, I., Messias, H., et al. 2019, PASP, 131, 075003 [Google Scholar]
 Gómez, J. L., Lobanov, A. P., Bruni, G., et al. 2016, ApJ, 817, 96 [Google Scholar]
 Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, A&AS, 117, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
 Janssen, M., Goddi, C., van Bemmel, I., et al. 2019, A&A, 626, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johnson, M. D., Fish, V. L., Doeleman, S. S., et al. 2015, Science, 350, 1242 [Google Scholar]
 Kim, J.Y., Krichbaum, T. P., Broderick, A. E., et al. 2020, A&A, 640, A69 [EDP Sciences] [Google Scholar]
 Kravchenko, E., Giroletti, M., Hada, K., et al. 2020, A&A, 637, L6 [EDP Sciences] [Google Scholar]
 Leppanen, K. J., Zensus, J. A., & Diamond, P. J. 1995, AJ, 110, 2479 [NASA ADS] [CrossRef] [Google Scholar]
 MartíVidal, I., Krichbaum, T. P., Marscher, A., et al. 2012, A&A, 542, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MartíVidal, I., Vlemmings, W. H. T., Muller, S., et al. 2014, A&A, 563, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MartíVidal, I., Roy, A., Conway, J., et al. 2016, A&A, 587, A143 [CrossRef] [EDP Sciences] [Google Scholar]
 Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002 [NASA ADS] [CrossRef] [Google Scholar]
 Molina, S. N., Agudo, I., Gómez, J. L., et al. 2014, A&A, 566, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Park, J., Byun, D. Y., Asada, K., & Yun, Y., 2021, ApJ, 906, 85 [Google Scholar]
 Pearson, T. J., & Readhead, A. C. S. 1984, ARA&A, 22, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Petrachenko, W. T., Niell, A. E., Corey, B. E., et al. 2012, Proc. Geod. Planet Earth, IAGS, 136, 2012 [Google Scholar]
 Rau, U., & Cornwell, T. J. 2011, A&A, 532, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 O’Sullivan, S. P., Gabuzda, D. C., & Gurvits, L. I. 2011, MNRAS, 415, 3049 [NASA ADS] [CrossRef] [Google Scholar]
 Thompson, A. R., Moran, J. M., & Swenson, G. W. 2017, in Interferometry and Synthesis in Radio Astronomy, eds. A. Richard Thompson, J. M. Moran, & G. W. Swenson, 3rd edn. (Springer), 2017 [Google Scholar]
 van Bemmel, I., Small, D., Kettenis, M., et al. 2018, 14th European VLBI Network Symposium & Users Meeting (EVN 2018), 79 [Google Scholar]
All Tables
Names of the EHT telescopes that participated in the 2017 campaign and their respective antenna mounts.
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 Dterm spectra as thick lines, resulting from a frequencydependent 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 Dterm 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 Dterms and the values fitted using the LPCAL similarity assumption (see text). Center: correlation for the Dterms fitted after fifty iterations of polarimetry selfcalibration. Right: L_{1} norm between true and fitted Dterms, as a function of the number of selfcalibration iterations, for the case of linear Dterm approximation (black points) and full nonlinear Dterm 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 Dterms. 

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 Dterm amplitudes reported in K20 (see their Table A.2) and the amplitudes of the frequencyaveraged multiIF Dterms estimated with polsolve. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.