Calibrating VLBI Polarization Data Using GPCAL. I. Frequency-Dependent Calibration

The Very Long Baseline Interferometry (VLBI) technique is a method for obtaining very high angular resolution by using telescopes that are widely separated. The telescopes in a VLBI array can split the sky signal into two orthogonal polarization channels, which are typically right and left circularly polarized (RCP and LCP, respectively), and each can record independently. The observed polarized visibilities, which are the time average of the complex correlations between the electric fields for each polarization incident at each telescope in the array, are directly related to the distributions of the four Stokes parameters in the sky by means of the van Cittert-Zernike theorem (van Cittert, 1934; Zernike, 1938; Thompson et al., 2017).

The relationship between the measured polarimetric visibilities and the source’s intrinsic polarization can be distorted due to imperfections in the instrument’s response. This distortion is primarily caused by antenna polarization that is not perfectly circular, which results in “instrumental polarization” in the visibility data. Instrumental polarization is commonly described in terms of an antenna’s response to a wave of polarization orthogonal to the nominal antenna response (e.g., Thompson et al., 2017, their Chapter 4.7), which is often referred to as “polarization leakages” or “D-terms”.

LPCAL, a task implemented in the Astronomical Image Processing System (AIPS, Greisen, 2003) based on the linearized leakage model (Leppanen et al., 1995), has been a standard method for calibrating VLBI data for nearly three decades. This method has been successful in a variety of studies using a variety of VLBI arrays (e.g., Gómez et al., 2016; Casadio et al., 2017; Jorstad et al., 2017; Lister et al., 2018; Park et al., 2018, 2019; Event Horizon Telescope Collaboration et al., 2021). However, LPCAL has some limitations that prevent accurate calibration. It assumes that the total intensity and linear polarization structures of a calibrator are similar, and also it cannot fit the instrumental polarization model to data from multiple calibrators at the same time. The former is referred to as the “similarity approximation”¹¹1As far as we are aware, the term ”similarity approximation” was first introduced in Leppanen et al. (1995). The original idea was proposed by Cotton (1993), who also cautioned that the assumption of a direct proportionality between total intensity and linear polarization emission may not hold. (Cotton, 1993; Leppanen et al., 1995), and it is generally violated for VLBI data, especially at high frequencies at which nearly all calibrators are resolved.

The Generalized Polarization CALibration pipeline (GPCAL; Park et al. 2021b) has been developed to overcome these limitations²²2Other calibration/imaging pipelines that have recently been developed are polsolve (Martí-Vidal et al., 2021), similar to GPCAL but based on CASA (Janssen et al., 2019; THE CASA TEAM et al., 2022; van Bemmel et al., 2022), the eht-imaging software library using the regularized maximum likelihood technique (Chael et al., 2016, 2018, 2022), D-term Modeling Code (DMC, Pesce 2021), THEMIS (Broderick et al., 2020), and Comrade (Tiede, 2022) using Markov chain Monte Carlo schemes. See Park & Algaba (2022) for more details. (Park et al., 2021b). The GPCAL pipeline is written in ParselTongue (Kettenis et al., 2006) and is based on AIPS and the Caltech Difmap package (Shepherd, 1997). It enables the use of more accurate linear polarization models of calibrators for D-term estimation without being limited by the similarity approximation. It can also fit the instrumental polarization model to data from multiple calibrator sources simultaneously in order to increase the accuracy of the fitting.

GPCAL was applied to the polarization analysis of the first M87 Event Horizon Telescope (EHT) results (Event Horizon Telescope Collaboration et al., 2021) as well as other recent studies (e.g., Takamura et al., 2023). We have also demonstrated that GPCAL is capable of detecting polarizations from sources that are very weakly polarized (Park et al., 2021a). The linear polarization structure of the subparsec core of the M87 jet observed with the Very Long Baseline Array (VLBA) at 43 GHz is compact and has a very low fractional polarization ($\sim 0.2$-$0.6\%$). It was not possible to well constrain this structure in previous studies with LPCAL using the same data set due to a breakdown in the similarity approximation, resulting in less accurate calibration of the instrumental polarization (Walker et al., 2018; Kravchenko et al., 2020).

In GPCAL, one assumption was made that can limit the accuracy of instrumental polarization calibration: it assumes that polarimetric leakages remain constant across the frequency bandwidth of a receiver and over time. The assumption of constant leakages over frequency may not apply to many VLBI arrays that have offered a wide bandwidth in recent years. In general, the magnitudes of leakages tend to increase as the operating frequency deviates further from the nominal frequency within the band. This is due to the design of the instruments, including polarizers, which are optimized to minimize leakages at the nominal frequency (e.g., Hull & Plambeck, 2015). Based on recent wide-bandwidth VLBA data, some previous studies have suggested that frequency-dependent D-terms can limit the dynamic range of linear polarization images (Kravchenko et al., 2020). Also, the assumption of constant leakages over time is incorrect since polarimetric leakages are not constant across telescope beams; they are direction-dependent (Smirnov, 2011). The effect is more severe for large telescopes and at high frequencies, where antenna pointing is usually inaccurate.

In this series of papers, we present GPCAL implementations of novel methods of frequency and time-dependent instrumental polarization calibration of VLBI data. In this paper (Paper I hereafter), we present a method for calibrating frequency-dependent instrumental polarization. An accompanying publication (Park et al. 2023; Paper II) describes a method of calibrating instrumental polarization that varies with time. These methods are applied to real VLBI data sets and demonstrate that, depending on the circumstances, they can significantly improve calibration accuracy.

This paper is organized as follows. In Section 2, we explain the radio interferometer measurement equation, which GPCAL is based on to model the instrumental polarization. In Section 3, we introduce the calibration procedure for frequency-dependent leakage calibration implemented in GPCAL. In Section 4, we validate the method using synthetic data. In Section 5, we present results of applying the method to real observations from the VLBA at 15 and 43 GHz. We conclude in Section 6.

The observed polarized visibilities on a baseline between stations $m$ and $n$ are given by the correlation products

where $E$ is an electric field, $R$ indicates RCP, $L$ indicates LCP, angular brackets denote a time average, and an asterisk denotes complex conjugate. One can represent the visibilities using a coherency matrix formalism:

The observed ${\bm{V}}_{mn}$ are corrupted by antenna gains and polarimetric leakages. The relationship between the true $\bar{\bm{V}}_{mn}$ and the observed ${\bm{V}}_{mn}$ can be described via a sequence of linear transformations, which is the so-called radio interferometer measurement equation (RIME; Hamaker et al. 1996, see also Sault et al. 1996; Hamaker & Bregman 1996; Hamaker 2000; Smirnov 2011),

where H is the Hermitian operator and ${\bm{J}}$ is the so-called Jones matrix (Jones, 1941):

$G$ is the complex antenna gain, $D$ is the leakage factor (“D-term”), and $\phi$ is the antenna field rotation angle. Subscripts denote antenna numbers, and superscripts denote polarization. The field rotation angle takes the form

where $\theta_{\rm el}$ is the source’s elevation, $\psi_{\rm par}$ the parallactic angle, and $\phi_{\rm off}$ is a constant offset for the rotation of antenna feed with respect to the azimuth axis. Cassegrain mounts have $f_{\rm par}=1$ and $f_{\rm el}=0$ and thus the field rotation angle is equivalent to the parallactic angle, except for the constant offset. Nasmyth mounts have $f_{\rm par}=1$ and $f_{\rm el}=+1$ for Nasmyth-Right type and $f_{\rm par}=1$ and $f_{\rm el}=-1$ for Nasmyth-Left type.

For circular feeds, the $\bar{\bm{V}}$ are related to the Fourier transforms of the Stokes parameters ($\tilde{I}$, $\tilde{Q}$, $\tilde{U}$, and $\tilde{V}$) via

Thus, the source’s brightness distribution of each Stokes parameter can be derived once ${\bm{G}}$ and ${\bm{D}}$ matrices are estimated and removed from the data. ${\bm{P}}$ consists of purely geometric terms and can be determined straightforwardly.

Nearly all calibrators are resolved with VLBI, and their full polarization structures change rapidly (typically over several weeks or months; for example, Jorstad et al. 2017; Lister et al. 2018; Weaver et al. 2022). This means that $\bar{\bm{V}}$ are unknown and cannot be easily modeled using a simple assumption such as $\tilde{I}$, $\tilde{Q}$, $\tilde{U}$, and $\tilde{V}$ are proportional to each other³³3This assumption holds when calibrators are point sources, which is typical for connected interferometers.. As a result, a standard calibration strategy for VLBI data has been to estimate ${\bm{G}}$ first using the parallel-hand visibilities and then ${\bm{D}}$ using the cross-hand visibilities after removing ${\bm{G}}$. It is because parallel-hand visibilities are sensitive to source’s total intensity emission and antenna gains, but insensitive to linear polarization emission and D-terms, whereas cross-hand visibilities are the opposite. GPCAL follows this strategy as well (Park et al., 2021b). Some recently developed calibration/imaging pipelines, such as polsolve (Martí-Vidal et al., 2021), THEMIS (Broderick et al., 2020), DMC (Pesce, 2021), and Comrade (Tiede, 2022), utilize the full RIME for fitting.

The antenna field-rotation angles are usually corrected at an upstream calibration stage (before performing global fringe fitting). Thus, Equation 3 becomes:

which can be rewritten element-wise as (e.g., Roberts et al., 1994):

The last two equations and their approximations are used in GPCAL (Park et al., 2021b) and LPCAL (Leppanen et al., 1995), respectively. In the original GPCAL pipeline, the D-terms are assumed to remain constant over frequency within the bandwidth and over time. We present new methods for calibrating frequency- (the present paper) and time-dependent polarimetric leakages (Paper II) that have been implemented in GPCAL.

We present the method of correcting for frequency-dependent polarimetric leakages in VLBI data. The method assumes that antenna gains are unity, same as the original GPCAL pipeline, assuming that they are accurately constrained during the upstream calibration procedure and the imaging and self-calibration procedure. There is still a single phase offset ($e^{j\phi_{RL,\rm ref}}$) remaining in the cross-hand visibilities of all baselines, which originates from the instrumental phase offset between polarizations at the reference antenna.

The equations of our interest are:

$\mathscr{RL}$ and $\mathscr{LR}$ are replaced by Equation 6 and $\mathscr{RR}$ and $\mathscr{LL}$ by the final calibrated parallel-hand visibilities $r^{RR}_{mn,{\rm cal}}$ and $r^{LL}_{mn,{\rm cal}}$, respectively. It is assumed that $e^{j\phi_{RL,{\rm ref}}}$ has already been corrected, which can be done by comparing the integrated electric vector position angles (EVPAs) of calibrators to known values. This can usually be accomplished by using sources with stable integrated EVPAs over time or from contemporaneous observations of single dishes or connected arrays. $\tilde{Q}_{mn}(\nu)$ and $\tilde{U}_{mn}(\nu)$ are obtained from the Fourier Transforms of the source’s linear polarization images. In order to obtain these images, CLEAN is applied to the polarimetric leakage corrected data acquired by running the GPCAL pipeline on the frequency-averaged data. It is assumed that the linear polarization images are constant over frequency for each baseband channel (often referred to as “IF”). We make the assumption that the fluctuations of polarimetric leakages across different frequencies are more significant than those of the source’s intrinsic linear polarization emission. To obtain the values of the former, we consider the latter to be a constant and solve for it accordingly⁴⁴4In spite of this, $\tilde{Q}_{mn}(\nu)$ and $\tilde{U}_{mn}(\nu)$ are not exactly constant over frequency since the visibility at different $\nu$ has a different $(u,v)$, although the variation due to the $(u,v)$ coverage effect is generally expected to be small.. Field-rotation angles are independent of the frequency.

The method fits Equation 3 to the data for each channel to determine $D^{R}(\nu)$ and $D^{L}(\nu)$ using the Scipy curve_fit package. As with the GPCAL pipeline, it is possible to fit the model to data from multiple calibrator sources simultaneously to enhance the fitting accuracy.

The method uses a multi-source UVFITS dataset (AIPS data catalog file) provided by users. To reduce calculation time and memory requirements, it bins the data in frequency, if specified by users. This approach is usually recommended because recent VLBI data often contain a large number of channels, while the D-terms are unlikely to vary significantly across adjacent channels. Determining the optimal bin size may be data-dependent and therefore, it is recommended that users determine this parameter themselves. If the phase offset ($e^{j\phi_{RL,{\rm ref}}}$) still exists in the provided data, the method can remove the offset using the offset value provided by users⁵⁵5As explained above, the offset value can be derived by performing the EVPA calibration of the frequency-averaged data., which may vary depending on the IF used. The multi-source UVFITS data is usually not self-calibrated. The method performs amplitude and phase self-calibration with the task CALIB in AIPS, if requested, using the total intensity CLEAN models (saved in Difmap fits image files) provided by users.

The method runs the GPCAL pipeline to derive the D-terms for each channel. It then calibrates the multi-source UVFITS data using the derived D-term spectra to obtain $\bar{V}$ in Equation 3. This results in the final output being a multi-source UVFITS dataset saved in an AIPS catalog after correction of the frequency-dependent D-terms.

The delay offset between polarizations at the reference antenna is corrected during data pre-processing in AIPS after global fringe fitting. It should be noted, however, that the frequency-dependent instrumental polarization may prevent accurate delay offset corrections. Thus, it is recommended to perform additional delay offset corrections using the output data after the frequency-dependent D-terms have been removed. It is recommended that users perform this additional correction after checking the output file, as it is not included in the calibration procedure. The method is implemented in GPCAL, which is publicly available at https://github.com/jhparkastro/gpcal.

We validate our method by employing synthetic data generated from actual VLBA data. Our objective is to evaluate the effectiveness of our method on synthetic data that closely resembles real data in terms of properties. To achieve this, we acquired the ground-truth models of total intensity and linear polarization for the sources, as well as the noise characteristics and frequency-dependent leakages from real VLBA 15 GHz data. These real data components served as the basis for generating the synthetic data. The data we selected is part of the Monitoring of Jets in Active Galactic Nuclei with VLBA Experiments⁶⁶6https://www.physics.purdue.edu/MOJAVE/ (MOJAVE; Lister et al. 2018) data sets, which have been monitoring a number of AGN jets with the VLBA at 15 GHz for decades. The data were collected from the observation of 25 sources on 2020 June 14 (project code: BL273AK). The data were recorded in both RCP and LCP, with two-bit quantization in four IFs, using the digital downconverter (DDC) observing system, at a recording rate of 4 Gbps, yielding a total bandwidth of 128 MHz for each IF and polarization.

We performed a standard data postcorrelation process with AIPS according to Park et al. (2021c) and hybrid imaging with CLEAN and self-calibration in Difmap. We ran the GPCAL pipeline on the self-calibrated data. We used several calibrators with linear polarization structures dominated by cores for the initial D-term estimation using the similarity approximation. Further instrumental polarization self-calibration was performed with 10 iterations using several calibrators with high signal-to-noise ratios in the cross-hand visibilities. The phase offset between polarizations at the reference antenna was corrected by comparing the integrated EVPAs of several calibrators using the D-term corrected data provided by GPCAL with the data already calibrated in the monitoring program’s public database. We produced Stokes $Q$ and $U$ images of the sources with CLEAN in Difmap.

We generated synthetic data for three sources, 3C 273, 3C 279, and OJ 287, using the synthetic data generation function implemented in GPCAL, as described by Park et al. (2021a). In brief, the function generates a synthetic data set based on Equation 7, using the $I$, $Q$, and $U$ CLEAN models as ground-truth source structures (i.e., their Fourier Transforms are $\tilde{I}$, $\tilde{Q}$, $\tilde{U}$ in Equation 6), and assuming $\tilde{V}=0$. The multi-channel real VLBA 15 GHz data are replaced by the source model visibility corresponding to each $(u,v)$, and we added thermal noise based on the uncertainties of the real data for each channel. We introduce frequency-dependent D-terms, which were derived from the real data using GPCAL (Section 5). For example, the frequency-dependent leakages of the VLBA BR station shown in the left panel of Figure 4 are introduced to the synthetic data set of the same station using Equation 7. Since the purpose of the synthetic data test is to validate our method of correcting for frequency-dependent leakages, we assumed unity antenna gains, i.e., ${\bm{G}}={\bm{I}}$ (However, refer to the discussion below and Appendix B for the analysis of synthetic data, incorporating antenna gain distortions.).

From the synthetic data set, we derived frequency-dependent D-terms with GPCAL. Firstly, the synthetic data set was averaged over frequency for each IF. The total intensity CLEAN models are reconstructed using Difmap for each source from this data set. Secondly, we corrected antenna leakages by running the GPCAL pipeline on the frequency-averaged data, from which we obtained Stokes $Q$ and $U$ CLEAN images. For this leakage correction, an initial D-term estimation is performed with OJ 287 based on the similarity approximation, followed by instrumental polarization self-calibration with 10 iterations using all three sources. Thirdly, frequency-dependent leakages were derived with the method from the frequency-unaveraged data set and using the $Q$ and $U$ CLEAN images. For this estimation, all three sources were used simultaneously. Lastly, we averaged the frequency-dependent leakage corrected data set over frequency and determined Stokes $Q$ and $U$ images based on that data set.

In Figure 1, we compare the components of the ground-truth D-terms (real and imaginary parts) with those of the reconstructed D-terms. Each panel indicates an average of the $L_{1}\equiv|D_{\rm Truth}-D_{\rm Recon}|$ norm. The method was able to reproduce ground-truth D-terms with an accuracy of an $\langle L_{1}\rangle$ norm of $\lesssim 0.1\%$ for all VLBA stations.

Figure 2 showcases the impact of frequency-dependent leakage correction on the resulting linear polarization images of 3C 273, as compared to the ground-truth linear polarization image employed in the generation of synthetic data (top left, labeled as “Ground-Truth”). Aside from the linear polarization image obtained after correcting for frequency-dependent D-terms (bottom left, “channel-by-channel D-term correction”), two reference polarization images were produced. The first is an image that has been reconstructed from a synthetic data set generated in accordance with the above description, but without taking into account leakage (top right, “No D-terms”). There are only two sources of error in this image: thermal noise in the data and errors introduced by the CLEAN algorithm. As a result, this image could function as a standard point of reference for evaluating the impact of polarimetric leakages present in the data. Another image is reconstructed using the synthetic data set, however, instead of a frequency-dependent leakage correction, the leakages are corrected using the data set averaged over frequency for each IF (bottom right, “D-term correction using frequency-averaged data”). It is through this method that LPCAL and the original GPCAL pipeline would perform leakage calibration. This image can then be used as a reference to test the effect of frequency-dependent leakage calibration.

In each panel, we note the dynamic range (DR) of the linear polarization images, defined as $P_{\rm peak}/\sigma_{P}$, where $P_{\rm peak}$ represents the peak polarization intensity and $\sigma_{P}\equiv(\sigma_{Q}+\sigma_{U})/2$ represents the average of the rms noise for the off-source regions in the Stokes $Q$ and $U$ images (Hovatta et al., 2012). The DR of the image obtained with frequency-dependent D-term correction (bottom left) is slightly lower than that of the top right image, while the DR of the image obtained with leakage correction using the frequency-averaged data (bottom right) is about half that of the top right image.

In Figure 3, we present the dirty images of the residual visibilities, which are the visibilities after subtracting the CLEAN model visibilities for Stokes $Q$ and $U$. These images showcase the noise levels within the residual visibilities, which can serve as an indicator of the accuracy of the leakage calibration. It is evident that the application of frequency-dependent leakage correction has resulted in a significant reduction in image noise compared to the image obtained with leakage correction using frequency-averaged data. Based on this difference, it can be concluded that frequency-dependent leakage correction with the method improves the quality of linear polarization images significantly when antenna gains are perfectly calibrated.

Compared to the past, modern VLBI arrays offer a much wider bandwidth. The large bandwidth of the data allows for a significant reduction in thermal noise. However, systematic errors in the data do not decrease as bandwidth increases, and this can pose a significant limitation in terms of image quality. In the case of linear polarization images, this effect is more severe, since conventional calibration methods are not able to easily remove the systematic errors. Some of the systematic errors are polarimetric leakages varying over frequency within the bandwidth and over time. Leakages usually become larger at frequencies far from the nominal frequency in the band, as the instruments are designed to have small leakages at the nominal frequency. Also, polarization leakage is direction-dependent effect (see, e.g., lecture 3 in Taylor et al. 1999; Thum et al. 2008; Smirnov 2011). Accordingly, D-terms can vary depending on the antenna pointing accuracy, which varies from scan to scan. Most existing calibration tools, however, assume that leakages are constant over time and frequency.

We present in this series of papers new calibration methods for frequency- (this paper) and time-dependent (Paper II) polarimetric leakages in VLBI data, which have been implemented in GPCAL. The method presented in this paper determines leakages for each channel based on the Stokes $Q$ and $U$ calibrator models, which are derived from the frequency-averaged data set. As in the original GPCAL pipeline, data from multiple calibrator sources can be used simultaneously for fitting to improve accuracy.

We verified the method using synthetic data generated based on real VLBA 15 GHz data. The method was able to reconstruct the ground-truth D-terms with an accuracy of $\langle L_{1}\rangle\lesssim 0.1\%$. The linear polarization images are significantly improved with the frequency-dependent leakage calibration using the method in the case of no antenna gain corruption assumed in the synthetic data.

We applied the method to two sets of real data obtained with the VLBA at 15 and 43 GHz. The D-term spectra were derived using several bright and moderately to highly polarized calibrators. There is continuous and smooth variation in the spectra across all bandwidths, which indicates that the spectra are an accurate representation of the characteristics of the VLBA’s instruments. The D-terms obtained from the data averaged over frequency within each IF are a good approximation of the D-term spectra. After correction of frequency-dependent D-terms, substantial variations in the cross-hand visibility spectra over frequency disappear. This result indicates that the variation is caused by frequency-dependent polarimetric leakages, rather than by intrinsic signals of the source. This result also indicates that the method is effective in removing the frequency-dependent leakages from real VLBI data.

We evaluated the effects of frequency-dependent instrumental polarization calibration on linear polarization images. There were no significant differences between the images obtained after correcting for frequency-dependent D-terms and those obtained from frequency-averaged data, although the dynamic range of the linear polarization images with frequency-dependent leakage corrections were slightly higher.

Based on these results, the calibration of D-terms based on frequency-averaged data appears to be a reasonable approximation for the particular VLBA data sets analyzed in this paper. The reason for this is that the primary contribution of D-terms to the cross-hand visibilities is linear. When the cross-hand visibilities are averaged over frequency before correcting frequency-dependent D-terms, there are second-order terms that can cause distortion. The distortion is expected to be small if the D-term amplitudes of frequency channels are small and do not vary significantly over frequency within each IF. Furthermore, other systematic errors in the data, such as residual antenna gain errors, may dominate the errors caused by averaging data over frequency before correcting for frequency-dependent leakage. This appears to be the case for the VLBA data sets presented in this paper. Nevertheless, we observe some improvement in the quality of polarization images as a result of calibration of frequency-dependent D-terms, indicating that users should consider calibration of frequency-dependent D-terms to increase calibration precision. Additionally, it will contribute significantly to the calibration of future VLBI data, for which significant variations in D-terms are expected over frequency within large bandwidths.