Calibrating VLBI Polarization Data Using GPCAL. II. Time-dependent Calibration

We present a new method of time-dependent instrumental polarization calibration for very long baseline interferometry (VLBI). This method has been implemented in the recently developed polarization calibration pipeline GPCAL. Instrumental polarization, also known as polarimetric leakage, is a direction-dependent effect, and it is not constant across the beam of a telescope. Antenna pointing model accuracy is usually dependent on time, resulting in off-axis polarimetric leakages that can vary with time. The method is designed to correct for the off-axis leakages with large amplitudes that can severely degrade linear polarization images. Using synthetic data generated based on real Very Long Baseline Array (VLBA) data observed at 43 GHz, we evaluate the performance of the method. It was able to reproduce the off-axis leakages assumed in the synthetic data, particularly those with large amplitudes. The method has been applied to two sets of real VLBA data, and the derived off-axis leakages show very similar trends over time for pairs of nearby sources. Furthermore, the amplitudes of the off-axis leakages are strongly correlated with the antenna gain correction factors. The results demonstrate that the method is capable of correcting for the off-axis leakages present in VLBI data. By calibrating time-dependent instrumental polarization, the rms noise levels of the updated linear polarization images have been significantly reduced. The method is expected to substantially enhance the quality of linear polarization images obtained from existing and future VLBI observations.


INTRODUCTION
Polarization observations using Very Long Baseline Interferometry (VLBI) are ideally suited for studying the processes of mass accretion and jet formation (see, e.g., Park & Algaba 2022 and references therein).These processes occur at small physical scales and can only be observed with instruments with high angular resolution.Polarized emission from plasma in mass accretion flows and jets is directly related to their magnetic fields.Recent observations of the nearby elliptical galaxy M87 with the Event Horizon Telescope (EHT) demonstrate the power of VLBI polarimetry (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f).The EHT total intensity image of M87 reveals the presence of a prominent ring-like structure.This structure is interpreted as the result of gravitational light bending of synchrotron emission from a hot plasma surrounding the black hole, along with photon capture occurring at the event horizon.In the corresponding linear polarization images of the ring, the polarization position angles are arranged in a nearly azimuthal pattern (Event Horizon Telescope Collaboration et al. 2021a).The comparison between these images and the General Relativistic Magnetohydrodynamic (GRMHD) simulation results indicates that a strong ordered, poloidal-dominated mag-Park et al.
netic field may exist around the black hole, which is capable of generating the powerful jets seen in this source.(Event Horizon Telescope Collaboration et al. 2021b).In addition, polarization observations of AGN jets using VLBI have provided important information about the collimation, acceleration, and particle acceleration processes in AGN jets (e.g., Asada et al. 2002;Jorstad et al. 2017;Gabuzda 2018;Lister et al. 2018;Park et al. 2019aPark et al. , 2021a;;Lisakov et al. 2021).
The linear polarization of AGN jets typically ranges from a few to a few tens of percent (e.g., Jorstad et al. 2017;Lister et al. 2018).In the past, the quality of VLBI linear polarization images of AGN jets was determined primarily by their sensitivity.However, recent VLBI arrays have increased their sensitivity considerably by enlarging the recording bandwidth and including very high sensitivity stations such as the phased Atacama Large Millimeter/submillimeter Array (ALMA; Matthews et al. 2018) and the phased Very Large Array (VLA).Consequently, systematic errors are becoming a more critical factor in determining the quality of VLBI linear polarization images.The most dominant systematic error affecting VLBI polarization data is "instrumental polarization".The source of this polarization signal is due to the antenna polarization not being exactly circular or linear, causing an unpolarized source to appear polarized.This polarization is commonly referred to as "polarization leakage" or "D-terms" and it must be corrected in the visibility data before producing an image.
LPCAL is a task incorporated into Astronomical Image Processing System (AIPS, Greisen 2003) and based on the linearized leakage model (Leppanen et al. 1995).It has been a standard program for calibrating instrumental polarization in VLBI data for a long time.In spite of its success for many studies utilizing various VLBI arrays (e.g., Casadio et al. 2017;Jorstad et al. 2017;Lister et al. 2018;Park et al. 2018Park et al. , 2019a;;Gómez et al. 2022;Zhao et al. 2022), there are some limitations which prevent accurate calibration.An important limitation is that the "similarity approximation"1 , which assumes similar total intensity and linear polarization structures for the calibrators (Cotton 1993;Leppanen et al. 1995), is often violated for VLBI.The problem becomes more severe at high frequencies, where nearly all calibrators are resolved.Recently, a number of new calibration and imaging pipelines have been developed to improve calibration accuracy and linear polarization images, including the Generalized Polarization Calibration pipeline (GPCAL; Park et al. 2021b), polsolve (Martí-Vidal et al. 2021), the eht-imaging software library (Chael et al. 2016(Chael et al. , 2018)), D-term Modeling Code (DMC, Pesce 2021), THEMIS (Broderick et al. 2020), and Comrade (Tiede 2022).Some of the pipelines were applied to the first linear polarization imaging of the M87 black hole (Event Horizon Telescope Collaboration et al. 2021a).
Nevertheless, most of the existing pipelines rely on the fundamental assumptions that instrumental polarization remains constant over a wide range of frequencies of a receiver and over a period of time during observation.VLBI arrays of recent generations can violate these assumptions, which results in systematic errors in the visibility data.This series of papers presents new methods for modeling the frequency-and time-dependent instrumental polarization in VLBI data, which have been implemented in GPCAL.The method for correcting frequency-dependent instrumental polarization is presented in a companion paper (Park et al. 2023;henceforth Paper I).The purpose of this paper is to introduce the method for correcting for time-dependent instrumental polarization.
Instrumental polarization itself is believed to not change significantly during the observation.However, it is a direction-dependent effect (e.g., Smirnov 2011), so the instrumental polarization varies depending on the direction of the antenna beam (due to the cross polarized sidelobes; see, e.g., Napier 1989;Thum et al. 2008).More specifically, instrumental polarization of an antenna consists of two main components (see Section 4.3 in Napier 1989).There is one at the center of the antenna beam that can be considered constant across the beam (D on−axis ).The other component varies across the beam (D off−axis ).The accuracy of antenna pointing during observation can be subject to variations over time, arising from stochastic winds and the deformation of antennas caused by sunlight, among other factors.Thus, effective D-terms can vary over time during an observation due to the direction-dependent effect.Only if the antenna pointing is sufficiently accurate throughout the observing period does instrumental polarization remain constant.Time-dependent leakage amplitudes are expected to be directly related to antenna pointing accuracy in this case.In general, high-frequency observations and antennas with large diameters usually have less accurate antenna pointing due to small beam sizes and stronger dish deformation at low elevations.Moreover, the instrumental polarization of phased arrays (such as the phased ALMA and VLA) is variable (e.g., Event Horizon Telescope Collaboration et al. 2019c), as is the phasing efficiency.As demonstrated in this paper, even the D-terms for the Very Long Baseline Array (VLBA) may vary greatly due to inaccurate antenna pointing and changing weather conditions.GP-CAL can increase the dynamic range of the VLBA linear polarization images by a factor of several by correcting the time-dependent instrumental polarization.
The paper is organized as follows.In Section 2, we describe the radio interferometer measurement equation, which forms the basis of the GPCAL instrumental polarization model.We discuss in Section 3 the strategy for calibrating time-dependent polarization leakages.We validate the performance of the method with a synthetic data set in Section 4. In Section 5, we apply the method to real VLBA data sets and verify that time-dependent leakage correction using GPCAL can significantly enhance the quality of VLBI linear polarization images.We summarize and conclude in Section 6.

THE RADIO INTERFEROMETER MEASUREMENT EQUATION
We use the radio interferometer measurement equation (RIME; Hamaker et al. 1996;Sault et al. 1996;Hamaker & Bregman 1996;Hamaker 2000;Smirnov 2011), as described in Section 2 in Paper I. In this paper, we provide a brief explanation of the equation for the convenience of readers.
The observed complex visibilities between two VLBI antennas, m and n, are expressed in a visibility matrix, V mn , of the form where R and L refer to the right-and left-handed circular polarizations (RCP and LCP), respectively.The observed V mn is corrupted by antenna gain and polarization leakage.It is convenient to arrange all the corruptions into a single Jones matrix (Jones 1941): where G is the complex antenna gain, D is the leakage factor (D-term), and ϕ is the antenna field rotation angle.Subscripts and superscripts denote antenna numbers and polarization, respectively.The field rotation angle is a combination of the source's elevation (θ el ), parallactic angle (ψ par ), and a constant offset for the rotation of antenna feed with respect to the azimuth axis (ϕ off ) via: Cassegrain mounts have f par = 1 and f el = 0. Nasmyth-Right type mounts have f par = 1 and f el = +1 and Nasmyth-Left type mounts have f par = 1 and f el = −1.
The observed V mn are modifications of the true Vmn through: where H is the Hermitian operator.
(5) As explained in Paper I, GPCAL assumes that the antenna field-rotation angles are already corrected at an upstream calibration stage (before performing global fringe fitting).Thus, Equation 4 becomes: 3. CALIBRATION PROCEDURE

Model Equation
GPCAL assumes that antenna gains are already corrected during the upstream calibration and imaging/self-calibration procedures.The original GP-CAL pipeline (Park et al. 2021b) fits Equation 6with the assumption of G = I to the observed cross-hand visibilities averaged over the frequency bandwidth.The pipeline assumes that polarimetric leakages are constant during the observation.However, if this assumption is violated, there are residual leakages in the data.Following the Appendix in Pesce (2021), we can write the true leakage matrix D as: and the estimated leakage matrix D with the assumption of constant leakage terms during the observation, which can be different from D, as: The visibility matrix after running the GPCAL pipeline would become: where R is a residual leakage matrix: The approximation holds when dropping out secondorder terms.With this approximation, the off-diagonal terms in R are the "residual" leakage terms, i.e., [R i,j = D i,j − Di,j ] ∈i̸ =j .
First, we use the GPCAL pipeline to remove polarimetric leakages that are assumed to remain constant over time ( D).We then fit Equation 9 to the corrected data in order to derive the "residual" time-dependent leakages.According to Pesce (2021), it may be more appropriate to redo calibration entirely with raw data (after data pre-processing) rather than incrementally calibrating a partially calibrated data set.In spite of this, we use this two-step procedure because, as we will demonstrate below, the signal-to-noise ratio of data plays a crucial role in the accurate estimation of timedependent leakages.Therefore, before performing timedependent instrumental polarization calibration, it is preferable to average the data over frequency.Nevertheless, as demonstrated in Paper I, modern VLBI arrays provide a large recording bandwidth, which can lead to significant variations in the D-terms over frequency.The frequency-dependent D-terms may introduce nonnegligible non-closing errors in the data if they are not removed before averaging the data over frequency.It is possible to accurately correct frequency-dependent Dterms using our method presented in Paper I or even the GPCAL pipeline, which uses averaged data within each IF, which can remove the gross variations of D-terms across the entire frequency band.Thus, we use the data after correcting for D, where the gross variations of Dterms over frequency have already been removed, and correct for the residual time-dependent leakages.
The cross-hand visibilities in Equation 9 can be rewritten as: where we replaced RL and L R using Equation 5and RR and L L by the final calibrated parallel-hand visibilities r RR mn,cal and r LL mn,cal , respectively.We assume that Q and Ũ are constant over time during the observation and depend only on the baseline coordinate (u, v).We fit Equation 11to the data for each scan to derive ∆ R (t) and ∆ L (t) using the Scipy curve fit package.

Calibration Strategy
We fit Equation 11to the data for each scan.There is a limited number of data points used for fitting in this case.In a scan with N antennas, there are 4N free parameters (the real and imaginary parts of the Dterms for RCP and LCP for each antenna).If one tries to fit a model with 4N degrees of freedom to a limited number of data points within a scan as compared to the entire dataset, there is likely to be a significant correlation among the parameters.We will demonstrate below that the best-fit D-terms in this case have large amplitudes because of the high correlation between parameters, while the assumed D-terms have small amplitudes in the synthetic data (Section 4).This is similar to performing an amplitude self-calibration at a very short solution interval on the total intensity data of weak sources, which results in antenna gain solutions with very high amplitudes.
Due to the limited signal-to-noise ratio of the data (primarily as a result of the limited number of total data points) compared to the number of free parameters, it is challenging to accurately constrain the polarimetric leakages of all stations for each scan.Additionally, antenna pointing should be reasonably accurate for most stations and most scans.Nevertheless, some stations may have more inaccurate antenna pointing than others due to poor weather conditions or large diameters.For stations with inaccurate pointing during particular scans, the linear polarization models produced after correcting for on-axis instrumental polarization may differ substantially from the cross-hand visibilities of all baselines associated with those stations.If this is the case, one can assume that there are significant residual leakages (∆ in Equation 11) only for those stations for the scans, and fix the residual D-terms for the other stations to be zero for fitting.By doing so, we will be able to avoid the strong correlation between fitting parameters and improve the accuracy of the fitting2 .
The first polarization imaging of M87 from the 2017 EHT observations (Event Horizon Telescope Collaboration et al. 2021a) was also carried out using a similar approach.In some calibration pipelines, short baselines are used to calibrate the D-terms of the stations comprising the short baselines (ALMA and the Atacama Pathfinder Experiment (APEX) telescope in Chile and the James Clerk Maxwell Telescope (JCMT) and Submillimeter Array (SMA) on Maunakea in Hawaii).Those pipelines then assume that the D-terms of those stations have already been corrected and perform fitting for only the D-terms of the other stations based on the long baseline data.Through this two-step calibration strategy, which takes advantage of the fact that short baselines are less sensitive to complex source structures and have a high signal-to-noise ratio, it was possible to achieve good calibration accuracy and avoid strong correlations between the D-terms of many stations.
Accordingly, we adopt the following calibration strategy for time-dependent calibration of instrumental polarization.
1.By using CLEAN with Difmap on the data corrected for the on-axis D-terms, Stokes Q and U images are produced.
2. For each scan and for each station m, we compute a norm (12) where w mn,k is the weight of the kth visibility matrix V k in the scan for the baseline between stations m and n, rRL mn,k ≡ Qmn,k + j Ũmn,k and rLR mn,k ≡ Qmn,k − j Ũmn,k are the model cross-hand visibilities corresponding to the kth visibility for the baseline derived from the Fourier Transforms of the Stokes Q and U CLEAN images.
3. Identify l number of stations having the largest norm values (i.e., the stations giving the largest and the second largest χ 2 m for l = 2).l is a free parameter controlled by the variable timecal freepar in GPCAL, which is provided by users.The determination of the optimal value for timecal freepar may vary depending on the data and the SNR of the source, as indicated in Table 1 and Section 4. In the case of VLBA data, it is recommended to adopt a conservative approach and set a small value, such as timecal freepar = 1-2.4. We fit Equation 11 to the cross-hand visibilities in the scan.Let the residual D-terms (∆ R and ∆ L ) for the stations identified in Step 3 be free parameters and those for the other stations be fixed to be zero during the fitting.
6. Correct for the derived residual D-terms by inverting the Jones matrices in Equation 9.
7. Iterate the above steps as many times as specified by users.Utilize the results obtained in Step 6 to generate Stokes Q and U images during Step 1, which are subsequently employed in the remaining calibration steps.The method yields the reduced χ 2 of the fit at every iteration, which can serve as a valuable tool for selecting the appropriate number of iterations.
While this procedure is similar to the instrumental polarization self-calibration procedure implemented in the original GPCAL pipeline (Park et al. 2021b), it differs in that fitting is performed only for leakages of stations that exhibit large residuals between the model and cross-hand visibilities.Using synthetic data sets, we will demonstrate the effectiveness of this strategy in correcting for time-dependent leakages of large amplitudes (Section 4).The method is implemented in GP-CAL, which is publicly available at https://github.com/jhparkastro/gpcal.

VALIDATION USING SYNTHETIC DATA
To validate the performance of the method, we used synthetic data derived from real VLBA data observed on 2018 Dec 08 at 43 GHz as part of the VLBA-BU-BLAZAR monitoring program3 .Data calibration and analysis methods are described in Paper I. Stokes Q and U images were produced with CLEAN in Difmap for 10 sources; 3C 279, 3C 273, OJ 287, 3C 454.3, 1633+382, 3C 345, CTA 102, 3C 84, MKN 501, and 0235+164.The sources provide a wide range of total intensity and linear polarization structures and cover a wide range of signal-to-noise ratios.
The synthetic data for these 10 sources were generated using GPCAL as described in Park et al. (2021a).GP-CAL generates synthetic data sets using Equation 6 and the I, Q, and U CLEAN models as ground truth source structures (i.e., Ĩ, Q, Ũ in Equation 5), assuming Ṽ = 0. Based on the uncertainties of each real data point, we added thermal noise to the corresponding synthetic data point.Afterwards, the data was averaged over the entire frequency band in order to increase the signal-to-noise Park et al. ratio.We assumed unity antenna gains, i.e., G = I.We introduced on-axis D-terms, which are randomly chosen based on the D-term distribution estimated by GPCAL from the real data, and are assumed to remain constant over time.We then introduced off-axis D-terms that vary randomly from scan to scan and whose real and imaginary components follow a Gaussian distribution with a zero mean and a certain standard deviation.As explained in Section 3 and will be shown in Section 5, the off-axis D-terms for most stations are anticipated to have small amplitudes in most cases.However, some stations that experience large pointing errors may have large amplitudes in their off-axis D-terms.In order to reflect this realistic situation in the synthetic data set, we assumed a standard deviation of 0.02 and 0.01 for the VLBA PT and NL station's leakages, respectively, and a standard deviation of 0.002 for the other station's leakages 4 .The method aims to identify and remove off-axis D-terms with large amplitudes from the data.First, we ran the GPCAL pipeline on the synthetic data set to correct for the on-axis D-terms.Based on the assumption that the source 3C 84 is unpolarized, the pipeline estimates the initial D-terms using the weakly polarized source 3C 84 (with the degree of linear polarization ≲ 1% at 43 GHz; Kim et al. 2019).As a next step, it implements 10 iterations of instrumental polarization self-calibration using bright calibrators with moderate or high linear polarization: 3C 279, 3C 273, 3C 454.3, and OJ 287.The data for 3C 84 is included as well in this procedure under the assumption that it is unpolarized.For all sources, the derived on-axis Dterms are removed from the data.
Using the method, we corrected for residual timedependent polarimetric leakages for all simulated sources.We used timecal freepar = 2.This means that during each iteration of calibration, the D-terms of the two stations with the largest χ 2 m values for each scan are free parameters, while the D-terms of the other stations are fixed to be zero (see Steps 2 and 3 in Section 3).The time-dependent leakage calibration procedure was repeated ten times.
We present in Figure 1 a comparison of the real and imaginary components of the time-dependent D-terms assumed in the synthetic data generation (x-axis) with the reconstructed D-terms estimated by GPCAL (yaxis).The reconstructed D-terms are the sum of the on-axis and off-axis terms.Each panel indicates the average L 1 ≡ |D Truth − D Recon | norm.Our method can reconstruct the true D-terms accurately, except for MKN 4 Obtaining the precise standard deviation values from actual data is a non-trivial task due to the dependence of off-axis D-term amplitudes on sources and scans.To address this issue, we rely on assuming values that would result in off-axis D-term amplitudes that are reasonably comparable to those obtained from the actual data.
501 and 0235+164, which have a low signal-to-noise ratio.
Figure 2 shows the probability density function (PDF) of the difference between the reconstructed and truth D-term components for each source.As indicated in each panel, we fitted a Gaussian function to each PDF and derived its standard deviation (σ).For most sources, σ is close to 0.2%, which is the standard deviation of a Gaussian distribution used to generate random time-dependent D-term components, except for NL and PT.Consequently, the method can correct for time-dependent D-terms with large amplitudes from some stations, but cannot accurately constrain timedependent D-terms with small amplitudes.A linear polarization image's quality is primarily limited by the former D-terms, while the latter is naturally expected in most realistic scenarios.Thus, the synthetic data test validates the effectiveness of the method, which is primarily designed to correct for time-dependent leakages with large amplitudes that usually appear at certain stations and/or scans.
For the weak sources MKN 501 and 0235+164, σ was substantially greater than 0.2%.We denote the average S/N of the parallel-hand and cross-hand visibilities for each source in each panel of Figure 2. It is evident that the S/Ns of the cross-hand visibilities for these sources are notably low 5 This result demonstrates that we cannot derive accurate solutions due to the limited signal to noise ratio for these sources despite the small number of parameters available for each round of fitting.Thermal noise would likely limit the quality of linear polarization images of these sources, and systematic errors such 5 However, we note that the S/Ns of the parallel-hand visibilities are also crucial for the precise reconstruction of time-dependent leakages.In our method, we assume that RR and L L are equivalent to the self-calibrated visibilities r RR cal and r LL cal , respectively, as expressed in Equation 11.This is the underlying reason why the reconstruction for 0235+164 exhibits superior results compared to that of MKN 501, despite their cross-hand SNRs being similar to each other.

Park et al.
as time-dependent leakages would not have a significant impact.As a result, it is not recommended to perform time-dependent leakage calibration on weak sources.
Probably the most critical parameter for the method would be timecal freepar, since high values of this parameter would result in large correlations between the fitting parameters and overfitting.As a demonstration of this effect, we applied the method by changing this parameter and presented the sigma values for each source from each run in Table 1.When timecal freepar = 2, the sigma values are the lowest for most of the bright sources.For weak sources MKN 501 and 0235+164, the values tend to become smaller as timecal freepar decreases.With timecal freepar = 10, i.e., when the D-terms of all 10 stations are free parameters, the values become very high for all sources.Using small timecal freepar values is recommended for VLBA data at 43 GHz, though the optimal parameter may differ depending on the quality and array of the data.

Data & Analysis
Using two real VLBA data sets, we evaluate the effectiveness of the method.One is the same data analyzed in Paper I and in Section 4 (project code: BM462M) observed on 2018 Dec 08.In the period between August 2018 and October 2019, the VLBA pointing model was not optimized due to problems that occurred during the transition of the telescope control computers (Blanchard 2021).Due to this, the quality of the VLBA data observed during this period was affected, particularly at high frequencies, where the antenna beam size is small.The inaccurate antenna pointing model also affected the BM462M data we analyzed.As a consequence, off-axis D-terms will be pronounced and will significantly limit the quality of linear polarization images derived from this data.
A second set of data is from the M87 monitoring program observed using the VLBA simultaneously at 4.7, 5.1, 6.5, and 6.8 GHz in order to investigate the kinematics of the M87 jet (Park et al. 2019b;Park et al. in prep.).We have selected the data observed on 2021 July 02 at 4.7 GHz (project code: BP249H).The VLBA Los Alamos (LA) station reported rain near the end of the observation, and the Stokes Q and U CLEAN models of all sources were not able to fit the data well after the rain began.Weather conditions may affect the leakage characteristics of this station, which motivates us to test our method.
Data postcorrelation was performed with AIPS and hybrid imaging with CLEAN and self-calibration in Difmap.
On-axis D-terms as well as frequency-dependent leakages were corrected using GPCAL, as described in Paper I. Following this, we corrected for time-dependent leakages based on our method, using the parameter timecal freepar = 2 and iterating the calibration procedure 10 times.

Results
Figure 3 shows the off-axis leakages as a function of time for the BM462M data.Although we present the results for three antennas only as an example, the results for other antennas were also qualitatively similar.The primary purpose of this test is to compare the trend of the off-axis D-terms derived from sources located near each other in the sky.The accuracy of the antenna pointing model depends on the direction of the antenna and the time of day.As a result, nearby scans of nearby sources may also experience similar levels of antenna pointing errors.We can expect to observe similar trends in off-axis D-terms for these scans if the direction-dependent leakages and inaccurate antenna pointing model are the primary causes of timedependent leakages.We compare the derived off-axis D-terms between 3C 273 and 3C 279 (separated by 10.4 degrees in the sky), 3C 454.3 and CTA 102 (6.8 degrees), and 3C 345 and 1633+382 (2.2 degrees).
For all antennas, the trends of the derived off-axis Dterms from nearby sources are very similar.In each panel, we present the average L 1 norm of the D-term components between nearby scans (within ten minutes) of these source pairs.We obtained smaller average L 1 norms for closer pairs of sources, which is reasonable since antenna pointing offsets, and therefore off-axis Dterms, are direction-dependent.
The method did not attempt to derive time-dependent D-terms for the Mauna Kea (MK) station for 3C 273 (the top middle panel in Figure 3).It is due to the fact that the χ 2 m value for this station is always low during any calibration round.The 3C 273 jet displays a complex linear polarization structure with a moderate to high level of linear polarization (e.g., Attridge et al. 2005;Hada et al. 2016;Park et al. 2021b).Thus, the cross-hand visibility of the long baselines associated with MK station has low amplitudes and low signal-to-noise ratios (at a ≲ 0.1 Jy level for this data).The effects of systematic errors, such as time-dependent leakages, are not significant for these baselines, and therefore GPCAL does not attempt to model them.
Figure 4 illustrates the amplitudes of the gaincorrection factors (1/|G|) and the amplitudes of the derived off-axis D-terms as a function of time for the VLBA Pie Town station.As a result of the inaccurate antenna pointing model, this station exhibits large gain correc- In each figure, the averaged norm L1 ≡ |Di − Dj| between the neighboring scans (separated by less than 10 minutes) for the adjacent sources i and j is shown for each polarization.The symbol ∆d denotes the distance between adjacent sources in the celestial sphere.While we provide results for only three antennas as an illustration, the findings for other antennas were also qualitatively comparable.
Park et al. tion factors and off-axis D-terms.Both quantities for nearby source pairs exhibit very similar trends as expected.In addition, the trends between the gain correction factors and the off-axis D-terms are also very similar.Due to the direction-dependent nature of antenna leakage, the effects of off-axis D-terms become more prominent for large antenna pointing offsets.
Figure 5 shows the amplitudes of the derived offaxis D-terms as a function of the gain correction factors for all scans and all sources.Regardless of the antenna, there is a strong positive correlation between the two quantities.It is demonstrated that timedependent leakages are present in the data due to the direction-dependent D-terms and imperfect pointing off- sets.Their amplitudes can be very large for scans that are affected by large pointing offsets.
In Figure 6, we present the derived off-axis D-terms as a function of time for 3C 273 and M87, separated by 10.3 degrees in the sky, for the BP249H data.In line with the results obtained for BM462M data, the off-axis D-terms for these sources also exhibit similar trends.As the antenna pointing is accurate at this low frequency (4.7 GHz), the amplitudes of the off-axis D-terms are generally small.For LA station, however, large off-axis D-terms are observed after around 3 UT, when rainfall began.Both 3C 273 and M87 exhibited this trend, indicating that severe weather had an impact on the large off-axis D-terms.Observations at other frequencies have also yielded similar results, although they are not included in the present paper.

Evaluation
Figure 7 illustrates the linear polarization images of 3C 279 and 3C 345 derived from the BM462M data before (left) and after (right) time-dependent leakage correction.A Ricean de-biasing correction has been applied to the images (i.e., P corr = P obs 1 − (σ P /P obs )), where P corr and P obs denote the corrected and original linearly polarized intensities (Wardle & Kronberg 1974;Lister et al. 2018), respectively, and σ P ≡ (σ Q + σ U )/2 is the average of the rms noise in the off-source regions in Stokes Q and U images (Hovatta et al. 2012).In each image, we note the polarization dynamic range, which is defined as the ratio of the peak linear polarization intensity to σ P .For each image, linear polarization vectors are drawn for linearly polarized intensities exceeding three times the σ P value.Through the use of the time-dependent leakage correction method, the noise level of the polarization images has been significantly reduced, resulting in an increase in polarization dynamic range by factors of two to three.After the time-dependent leakage correction, linear polarization emission was observed in the jet regions of weak total intensity that was not evident in the original linear polarization images.In the case of 3C 279, which has a high linearly polarized flux (∼ 0.8 Jy), the noise level in the polarization images is more dominated by systematic errors than by thermal noise.Using the method, residual systematic errors in the data are removed, making it possible to detect the weak polarization emission in the jet.
A comparison of linear polarization images of 3C 273 and M87 based on the BP249H data with and without time-dependent leakage correction is presented in Figure 8.Because of the improvement in polarization dynamic range, weaker linear polarization emission was detected in both sources, similar to the results obtained from the BM462M data.

CONCLUSIONS
In this series of papers, we present new methods for calibrating frequency and time-dependent leakages in VLBI data, which have been implemented in GPCAL.A wide fractional bandwidth is provided by modern VLBI arrays.The instrumental polarization is affected by chromatic effects since telescope systems are usually designed to produce the most optimal polarization response at a given nominal frequency (e.g., Martí-Vidal et al. 2021).In Paper I, we present a method for correcting leakages that vary in frequency.It is also possible for instrumental polarization to change over time.The reason for this is that instrumental polarization is a direction-dependent effect (e.g., Smirnov 2011).Therefore, instrumental polarization can vary depending on the accuracy of the antenna pointing.Antenna pointing is subject to uncertainty, and can be adversely affected by the weather, such as strong winds and the deformation of antennas caused by sunlight, which may result in instrumental polarization that is time-dependent.In poor weather conditions and with large dishes, this effect is more pronounced.The purpose of this paper is to introduce a method for correcting time-dependent leakages.This method works on the data after correcting the "on-axis" D-terms, which are calculated assuming that the D-terms are constant during observation.In this manner, "residual" leakages are derived.It calculates how well the Stokes Q and U images fit the cross-hand visibilities for each antenna and for each scan based on the data.Next, it determines which antennas exhibit the greatest difference between the data and the model.By assuming that the D-terms of all other antennas are zero, the method determines the best-fit D-terms for Park et al.
these antennas6 .The procedure is iterated until the solutions are converged, and the number of iterations can be controlled by the user.
The method was tested using a synthetic data set containing time-dependent leakages generated from real VLBA data observed at 43 GHz.Leakages consist of two components: the on-axis, stable D-terms drawn from the distribution of the D-terms from the real data, as well as the off-axis, variable D-terms that are randomly selected from Gaussian distributions for each scan.It is assumed that two antennas have large standard deviations (0.02 and 0.01 for PT and NL stations, respectively) for the Gaussian distributions and other antennas have small standard deviations (0.002), in order to mimic some recently observed VLBA data sets at 43 GHz.
Based on the method, we derived time-dependent leakages for the synthetic data set.We found that the method was able to successfully reconstruct the timedependent leakages that were assumed in the generation of synthetic data.In most cases, the distribution of the difference between the estimated and ground-truth Dterm components (the real and imaginary parts) could be described by Gaussian distributions with standard deviations less than 0.002.This result indicates that the method is capable of capturing the time-dependent leakages of large amplitudes, which are the dominant factors in limiting the quality of linear polarization images.The purpose of developing the method is actually to correct for the largest systematic errors in the data by using prior information that the off-axis D-term amplitudes of most antennas in realistic VLBI arrays are usually small.
Our method was applied to two sets of real data obtained with the VLBA at 4.7 and 43 GHz.In the VLBA 4.7 GHz data, significant deviations are observed between the source polarization models and cross-hand visibilities for all baselines associated with the VLBA Los Alamos station after the time the station reported rain.An antenna control computer upgrade caused inaccurate antenna pointing models in the VLBA 43 GHz data.We derived the off-axis D-terms for several pairs of sources located close to one another in the sky as a function of time.The trends of the derived off-axis D-terms between nearby sources are very similar.Off-axis D-terms for nearby scans of less separated sources displayed a higher degree of consistency.Additionally, the amplitudes of the off-axis D-terms and the gain-correction factors show very similar trends.Based on these results, it can be concluded that the method is capable of capturing time-dependent leakages caused by directiondependent leakages and imperfect antenna pointing.
For all sources analyzed in this paper, the dynamic ranges of linear polarization images have been improved by factors of ≳ 2 after correcting for time-dependent leakages using the method.According to these results, the method can significantly enhance the polarization image fidelity of VLBI data.Additionally, the method will be useful for improving global VLBI arrays operating at millimeter wavelengths, such as the GMVA and EHT, which have very small antenna beams because of the high observing frequencies and large dishes.

Figure 1 .Figure 2 .
Figure1.Comparison of the time-dependent, ground-truth D-term components (real and imaginary parts) assumed in the synthetic data sets shown on x-axis with the reconstructed D-term components derived by GPCAL shown on y-axis for each scan.Each station's result is presented using distinct colors.For each source, a correlation is shown between the true and estimated D-terms in percentage units.Over all antennas, the L1 ≡ |D Truth − DRecon| norm is averaged over the real and imaginary components of the D-terms for RCP and LCP.

Figure 3 .
Figure 3.The derived off-axis D-terms for the BM462M data as a function of time.In the diagram, the filled circles and open squares represent the D-terms for RCP and LCP, respectively.The upper and lower panels in each panel display the real and imaginary parts of the results for each antenna.Different rows of results are presented for different pairs of adjacent sources in the sky (3C 273 and 3C 279 in the top, 3C 454.3 and CTA 102 in the middle, and 3C 345 and 1633+382 in the bottom).In each figure, the averaged norm L1 ≡ |Di − Dj| between the neighboring scans (separated by less than 10 minutes) for the adjacent sources i and j is shown for each polarization.The symbol ∆d denotes the distance between adjacent sources in the celestial sphere.While we provide results for only three antennas as an illustration, the findings for other antennas were also qualitatively comparable.

Figure 4 .Figure 5 .
Figure 4. Amplitudes of gain correction factors (upper) and off-axis D-terms for the VLBA Pie Town station for different pairs of adjacent sources in the BM462M data set as a function of time.The filled circles and open squares represent the results for RCP and LCP, respectively.

Figure 6 .
Figure6.Same as Figure3but for the BP249H data.After UT 3, rain was reported at the VLBA Los Alamos station (LA; middle panel), following which the off-axis D-term amplitudes increased substantially.

Figure 7 .
Figure 7. Linear polarization images of 3C 279 (upper) and 3C 345 (lower) from the BM462M data set.The hot color indicates the total intensity distribution and the colored ticks indicate the electric vector position angles (EVPAs), with the color scales indicating Ricean de-biased linearly polarized intensity.Presented in the left and right panels are images obtained before and after time-dependent leakage correction.There is a note in each panel indicating the linear polarization dynamic range, which is defined as the ratio of the peak linear polarization intensity to the linear polarization off-source rms noise.In the bottom right corner, the shape of the synthesized beam is shown.

Table 1 .
Goodness of reconstruction of time-dependent leakages by changing the timecal freepar parameter Note-Standard deviations of Gaussian distributions fitted to PDFs of D Truth − DRecon are shown in units of % depending on the timecal freepar parameter.Different columns are used to display the results from different sources.