© The Authors 2026

1 Introduction

Polarization provides valuable information about astronomical objects and the ambient medium in astronomy, as it is one of the most direct ways to access information about astrophysical magnetic fields. In conjunction with very long baseline interferometry (VLBI), which is able to achieve nominal resolutions of ~20 μas at 1 mm (Event Horizon Telescope Collaboration 2019) and ~15 μas at 345 GHz (Raymond et al. 2024), we have the ability to probe extreme magneto-ionic environments.

Polarimetric studies with VLBI have revealed the presence of helical and toroidal magnetic fields in the jet of active galactic nuclei (AGN) (Asada et al. 2008; Hovatta et al. 2012; Gabuzda et al. 2017), the magnetic field of the supermassive black holes M87* and Sgr A* (Event Horizon Telescope Collaboration 2021; Event Horizon Telescope Collaboration 2024; Event Horizon Telescope Collaboration 2025), and an alignment between polarization and the jets of AGN, suggesting the presence of magnetizing shocks in AGN jets (Lister & Homan 2005; Pushkarev et al. 2023). However, polarization calibration of VLBI data is challenging, as polarization signals typically have S/N that are an order of magnitude lower than for total intensity.

Polarized data after precalibration can still have significant residual corruption from leakage between one polarization feed to another, known as polarization leakage or “D-terms”. Moreover, the complex source structure of calibrators in VLBI, especially at millimeter-wavelengths, hinders the estimation of this polarization leakage. The conventional polarization calibration approach in radio interferometry is observing a static point-like calibrator with known polarization properties and identifying leakage corruptions by using the information about the calibrator and disentangling source polarization from the instrumentation contributions that are being solved by the time-varying polarization of the feed angle rotation. However, due to the high angular resolution, it is unusual to find static point-like calibrators for polarization calibration in VLBI (Cotton 1993).

Leppanen et al. (1995) introduced the LPCAL software to infer polarization leakage corruptions in VLBI datasets. The LPCAL software utilizes a polarization prior model with a set of components that assume that the linear polarization of a component is proportional to the total intensity of the corresponding component, known as the similarity approximation (Cotton 1993; Leppanen et al. 1995). The simple and effective prior is advantageous to determine leakage corruption from sparse and noisy VLBI datasets due to the small number of degrees of freedom in the model.

However, the similarity approximation in LPCAL is not valid for sources with a complex polarization structure (Cotton 1993). Furthermore, multisource calibration is not directly supported (see Lister & Homan 2005, for an example of a multisource method with LPCAL.) and D-terms are considered static in time and frequency for each intermediate frequency (IF) by default in LPCAL. Recently, the new CLEAN-based polarization calibration software packages GPCAL (Park et al. 2021b) and PolSolve (Martí-Vidal et al. 2021) have been able to employ reconstructed Stokes Q and U images as priors in leakage calibration, a technique called polarization self-calibration. These software support multisource models to increase parallactic angle coverage, and are able to solve for D-terms that vary with time and frequency.

Nevertheless, GPCAL and PolSolve still rely on the conventional CLEAN deconvolution algorithm. Although CLEAN is easy to use and converges robustly, it comes with several limitations. It requires user-dependent inputs, such as CLEAN windows and weighting schemes, which may introduce biases into the final results. Furthermore, CLEAN assumes that the sky is a collection of point sources (usually of the size of the synthesized beam) or Gaussian blobs in multiscale CLEAN (Cornwell 2008), but has no explicit notion of an extended structure. This assumption may lead to severe artifacts on small scales that are corrected by a convolution with a restoring beam. This convolution results in suboptimal resolution, and the potential of radio interferometric observations to super-resolve sources (in the high S/N regime) cannot be exploited. Additionally, the CLEAN restored image does not necessarily fit the data, but rather a tapered version of it. Consequently, the CLEAN image model for leakage calibration is not suitable for radio interferometric data with a complex source structure. Moreover, conventional CLEAN-based self-calibration is performed iteratively by switching between flagging and manually choosing gain solution intervals, resulting in inconsistent calibration solutions (Martí-Vidal & Marcaide 2008; Popkov et al. 2021) and a lack of reproducibility. The conventional approach utilizes the sum of the CLEAN components without convolution as a model in self-calibration. This may imprint spurious small-scale structures in the self-calibrated data. Lastly, calibration uncertainties are not taken into account in the final results, and uncertainty estimation is not supported.

Modern forward modeling imaging algorithms provide us with a new approach to overcome the limitations of CLEAN-based polarization calibration methods (Birdi et al. 2020; Pesce 2021; Event Horizon Telescope Collaboration 2021; Event Horizon Telescope Collaboration 2024). Explicit regularizers and prior assumptions in forward modeling algorithms enable us to infer more robust leakage solutions from sparse VLBI datasets. For example, polarization constraints can be enforced in the polarization imaging model to avoid unphysical polarimetric images (Birdi et al. 2018). Recent instrumental advancement in VLBI accomplishes significantly improved S/N and more wideband observations. In addition, the latest algorithmic developments offer methods for solving highly degenerate inverse problems efficiently. A modern polarization calibration and imaging pipeline that utilizes the full potential of the data is highly desirable.

In this work, we introduce a novel Bayesian calibration and imaging method using the Bayesian imaging software resolve. The posterior distribution of antenna-based gains, D-terms, and Stokes images are explored jointly using variational inference algorithms (Knollmüller & Enßlin 2019; Frank et al. 2021). Recently, the polarization constraint ($I\ge \sqrt{Q^2+U^2+V^2}$) is encoded in the the resolve polarization imaging model (Arras et al. 2025). We utilize this polarization imaging model and incorporate polarization calibration and gain self-calibration directly into the imaging process. Furthermore, multisource and multi-IF polarization calibration are supported in a probabilistic framework.

This article is structured as follows. In Section 2, we explain the conventional CLEAN-based polarization calibration method and our Bayesian method. In Section 3 and Section 4, we validate the method with synthetic and real datasets, respectively. We summarize our results in Section 5.

2 Method

2.1 Radio interferometer measurement equation (RIME)

A radio interferometer measures Fourier components of the sky brightness distribution instead of imaging the sky directly. According to the van Cittert-Zernike theorem, the two-point correlation function of the signals recorded by two antennas, i and j, called visibility data V_ij, is given by the Fourier-transformed sky brightness distribution I for total intensity under the assumption that the field of view is small (Hamaker et al. 1996; Smirnov 2011; Thompson et al. 2017) $$V(u,v{)}_{ij}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}I(x,y)e^{-2\pi \mathbf{i}(ux+vy)}dxdy=\mathbb{F}\mathbb{T}[I],$$(1)

where (u, v) are the Fourier domain coordinates, (x, y) are the image domain coordinates, and 𝔽𝕋 is the Fourier transform operator.

Since the observed data are incomplete and corrupted by atmospheric and instrumental effects, we obtain a low-fidelity image, also known as the dirty map, from a direct Fourier transform of the visibility data. Inferring the real source structure from radio interferometric data is an ill-posed inverse problem, and therefore a unique solution does not exist. Thus, additional assumptions, prior knowledge, or other regularizers of the solution space are required to obtain high-fidelity images from radio interferometric data in image reconstruction.

The measurement equation (Eq. (1)) can be generalized for full polarization observation, and the instrumental and atmospheric data corruption can be described by Jones matrices (Jones 1941). As a result, the radio interferometer measurement equation (RIME) for the full polarimetric visibility matrix V on the basis of circular polarization is given by (Smirnov 2011) $${\mathbf{V}}_{ij}={\mathbf{J}}_i({\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{I}(x,y)e^{-2\pi \mathbf{i}(ux+vy)}dxdy){\mathbf{J}}_j^{†}+{\mathbf{N}}_{ij},$$(2)

where V_ij is the visibility matrix consisting of four complex correlation functions by the signal from the right-hand circular polarization (RCP) R and the signal from the left-hand circular polarization (LCP) L $${\mathbf{V}}_{ij}=\left(\begin{array}{ll}{R}_i{R}_j^{\ast }& R_i{L}_j^{\ast }\\ L_i{R}_j^{\ast }& L_i{L}_j^{\ast }\end{array}\right),$$(3)

where the asterisk * denotes a complex conjugate, J_i is the Jones matrix that describes the data corruption of the antenna i by instrumental and atmospheric effects $${\mathbf{J}}_i={\mathbf{G}}_i{\mathbf{D}}_i{\mathbf{P}}_i=\left(\begin{array}{cc}{g}_i^R& 0\\ 0& g_i^L\end{array}\right)\left(\begin{array}{cc}1& D_i^R\\ D_i^L& 1\end{array}\right)\left(\begin{array}{cc}{e}^{-\mathbf{i}{\varphi }_i}& 0\\ 0& e^{\mathbf{i}{\varphi }_i}\end{array}\right),$$(4)

where G_i is the antenna-based gain matrix, D_i is the leakage (D-term) matrix describing the signal leakage between polarizers (e.g., $D_i^R$ is the signal leakage from LCP to RCP for antenna i), P_i is the field rotation angle matrix, † denotes conjugate trans-position, N_ij is the additive noise, and I is the polarimetric sky brightness distribution matrix, consisting of Stokes I, Q, U, and V $$\mathbf{I}=\left(\begin{array}{cc}I+V& Q+\mathbf{i}U\\ Q-\mathbf{i}U& I-V\end{array}\right).$$(5)

The antenna field rotation angle ϕ, representing the rotation of the receiver polarization feeds with respect to the source due to Earth’s rotation, is defined as ϕ = f_elθ_el + f_parψ_par + ϕ_off, where θ_el is the elevation angle, ψ_par is the parallactic angle, and ϕ_off is a constant offset. The field rotation angle depends on the antenna mount. Alt-azimuth (ALT-AZ) mounts with a Cassegrain focus have f_par = 1 and f_el = 0. The ALT-AZ mounts with a Nasmyth-Right-type focus have f_par = 1 and f_el = 1, and ALT-AZ mounts with a Nasmyth-Left-type focus have f_par = 1 and f_el = −1. More details regarding the antenna mount types can be found in Appendix C of Janssen et al. (2019).

We note that no real antenna feed responds to a single sense of perfectly circular polarization. The voltage induced in each nominally circularly polarized channel consists of two terms: the intended response of the feed and a smaller response to orthogonal polarization (Roberts et al. 1994). As a result, the polarization state received by each feed is slightly elliptical rather than purely circular (e.g., Thompson et al. 2017). The complex D-term coefficients parameterize this departure from ideal circular polarization: their amplitudes are related to the axial ratio of the polarization ellipse traced by each feed’s effective response, and their phases describe the orientation of the ellipse on the sky (Roberts et al. 1994; Park et al. 2023a).

At centimeter and shorter wavelengths, a typical circularly polarized feed consists of a horn followed by a polarizing device, such as a quarter-wave plate or a septum polarizer, that converts the intrinsically linearly polarized detector response to circular polarization (Roberts et al. 1994). Imperfections in this conversion process—including differential amplitude and phase errors between the two orthogonal probes and deviations of the phase shift from the ideal π/2—give rise to the D-terms. Since the polarizer performance is inherently frequency-dependent, the resulting D-terms generally vary across the observing band, which motivates per-IF calibration for observations with a wide fractional bandwidth (Park et al. 2023a; see also Section 4.1). In addition, the D-terms are direction-dependent (Smirnov 2011), such that the effective leakage of an antenna can vary with time due to changes in the antenna pointing offset caused by winds and dish deformation (Park et al. 2023b).

The polarization calibration and imaging model with the polarimetric visibility matrix data V, Jones matrices J_i, J_j for antennas, i and j, and model visibilities (ℛℛ, ℛℒ, ℒℛ, and ℒℒ) from Stokes images is $$\begin{array}{rl}{R}_i{R}_j^{\ast }=& g_i^R{g}_j^{R\ast }[e^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}+D_i^R{e}^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}\\ & +e^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}{D}_j^{R\ast }+D_i^R{e}^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}{D}_j^{R\ast }]\\ R_i{L}_j^{\ast }=& g_i^R{g}_j^{L\ast }[e^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}{D}_j^{L\ast }+D_i^R{e}^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}{D}_j^{L\ast }\\ & +e^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}+D_i^R{e}^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}]\\ L_i{R}_j^{\ast }=& g_i^L{g}_j^{R\ast }[D_i^L{e}^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}+e^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}\\ & +D_i^L{e}^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}{D}_j^{R\ast }+e^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}{D}_j^{R\ast }]\\ L_i{L}_j^{\ast }=& g_i^L{g}_j^{L\ast }[D_i^L{e}^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}{D}_j^{L\ast }+e^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{R}{e}^{\mathbf{i}{\varphi }_j}{D}_j^{L\ast }\\ & +D_i^L{e}^{-\mathbf{i}{\varphi }_i}\mathcal{R}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}+e^{\mathbf{i}{\varphi }_i}\mathcal{L}\mathcal{L}{e}^{-\mathbf{i}{\varphi }_j}],\end{array}$$(6)

where the model visibility matrix consists of Stokes images $${\mathbf{V}}_{model}=\left(\begin{array}{cc}\mathcal{R}\mathcal{R}& \mathcal{R}\mathcal{L}\\ \mathcal{L}\mathcal{R}& \mathcal{L}\mathcal{L}\end{array}\right)=\mathbb{F}\mathbb{T}\left[\right(\begin{array}{cc}I+V& Q+\mathbf{i}U\\ Q-\mathbf{i}U& I-V\end{array}\left)\right].$$(7)

Polarization calibration and imaging are equivalent to estimating the gains g, D-terms D, Stokes images I, Q, U, and V from the visibility data V. In the conventional CLEAN-based method, this polarimetric calibration and imaging problem in VLBI is solved in an iterative fashion since the problem is highly degenerate. In this work, we infer the gains g, leakages D, Stokes I, Q, U, and V jointly from precalibrated visibility matrix data V in a probabilistic approach.

2.2 CLEAN-based polarization calibration method

The CLEAN-based polarization calibration method assumes that the antenna gains have already been calibrated during the data preprocessing and imaging and/or self-calibration procedures. It also assumes that the antenna’s field rotation angles have been corrected during the data preprocessing phase (see Park et al. 2021b, 2023a; Martí-Vidal et al. 2021 for more details). With these assumptions, the cross-hand visibilities consist of both the source’s intrinsic linear polarization terms and the terms associated with the antenna’s polarimetric leakages and field rotation angles.

Determining leakages is not straightforward, as they need to be disentangled from the source polarization terms. One of the most widely used methods in the past assumes that the source’s linear polarization emission is proportional to the total intensity structure within each “submodel”, which consists of a group of neighboring total intensity CLEAN components. This method is known as the “similarity approximation” (Cotton 1993; Leppanen et al. 1995) and is implemented in LPCAL, a task within the Astronomical Image Processing System (AIPS; Greisen 2003).

While generally a good approximation, the similarity method faces challenges when dealing with very weakly polarized sources (e.g., Park et al. 2021a), particularly when utilizing global VLBI observations at millimeter wavelengths, which provide ultra-high angular resolution. In such cases, the source’s linear polarization structures are often complex, and the similarity approximation may not hold well (e.g., Event Horizon Telescope Collaboration 2021; Event Horizon Telescope Collaboration 2024; Zhao et al. 2022).

To overcome this limitation, two methods have been developed: GPCAL (Park et al. 2021b), based on AIPS and Difmap, and PolSolve (Martí-Vidal et al. 2021), based on CASA (CASA Team 2022). These methods conduct calibration as follows:

1. They derive the leakages using the similarity approximation, such as LPCAL, and remove them from the data.

2. They perform imaging with CLEAN using the leakage-corrected data to obtain the source’s Stokes Q and U CLEAN models.

3. They use the original data to derive the leakage solutions again, this time using the source’s linear polarization models obtained in the previous step, and remove the leakages using the updated solutions.

4. They iteratively update the source’s linear polarization models and leakage solutions by repeating steps 2 and 3 until the solutions converge.

These methods can also simultaneously utilize data from multiple calibrator sources, which typically leads to improved leakage calibration accuracy, as leakage solutions are not expected to vary between sources.

2.3 Bayesian imaging software resolve

The open-source Bayesian imaging software resolve¹ treats radio interferometric imaging and calibration (Eq. (2)) as an inverse problem and computes a probabilistic solution for it. Thus, resolve computes the probability distribution of the sky brightness I given the measured visibilities. In Arras et al. (2019b), a joint calibration and imaging approach was introduced, computing a posterior probability distribution for the sky brightness I and the antenna gains G. Roth et al. (2023) extended this approach, incorporating direction-dependent effects in the antenna gains. A further extension to resolve was made in Arras et al. (2025), enabling full Stokes imaging. In Roth et al. (2024), the idea of major and minor cycles used in CELAN-based algorithms was adapted in a Bayesian version for the resolve framework.

Building on the full Stokes imaging capabilities, this work introduces polarization calibration and imaging to resolve. Thus, we compute the posterior distribution of the full Stokes sky brightness matrix I jointly with the posterior distributions of the antenna-based gain matrix G, and the leakage matrix D. The posterior distribution can be expressed via Bayes’ theorem $$\mathcal{P}(\mathbf{G},\mathbf{D},\mathbf{I}\mid \mathbf{V})=\frac{\mathcal{P}(\mathbf{V}\mid \mathbf{G},\mathbf{D},\mathbf{I})\mathcal{P}(\mathbf{G},\mathbf{D},\mathbf{I})}{\mathcal{P}(\mathbf{V})},$$(8)

in terms of the likelihood 𝒫(V|G, D, I) and the prior 𝒫(G, D, I). In the following Section 2.4, we discuss the likelihood model. In Section 2.5, we outline the prior models, and in Section 2.6, we describe the algorithm for approximating the posterior distribution given the likelihood and prior.

2.4 Likelihood distribution

As outlined in Section 2.3, we approach the polarization calibration and imaging problem as a Bayesian inference task. This entails defining the prior and likelihood to compute the posterior distribution, which represents the result of the Bayesian inference process.

Following the framework established by Knollmüller & Enßlin (2018), we employ a coordinate transformation, or reparameterization trick (also known as inverse transform sampling), to introduce new parameters ξ. This transformation ensures that the prior follows a standard normal distribution 𝒫(ξ) = 𝒢(ξ, 1), integrating all prior knowledge into the likelihood component 𝒫(V|ξ).

The likelihood 𝒫(V|ξ) is composed of two components: 𝒫(V|G, D, I) and a function of mapping ξ to specific values of 𝒢, 𝒟, and ℐ. The former component describes the data aspect of the likelihood, while the latter incorporates the prior information for the gains g, the leakages D, and the Stokes images I, Q, U, and V. This section elaborates on the first component, while the following section will deal with the second.

Our data model integrates gain matrices G, leakage matrices D, field rotation angle matrices P, and the polarization imaging model I. Upon setting G, D, P, and I, the visibility data V can be formulated as follows $${\mathbf{V}}_{ij}={\mathbf{G}}_i{\mathbf{D}}_i{\mathbf{P}}_i\mathbb{F}\mathbb{T}[\mathbf{I}]{\mathbf{P}}_j^{†}{\mathbf{D}}_j^{†}{\mathbf{G}}_j^{†}.$$(9)

By encapsulating Jones matrices into the response function, the model can be simplified $${\mathbf{V}}_{ij}={\mathbf{J}}_i\mathbb{F}\mathbb{T}[\mathbf{I}]{\mathbf{J}}_j^{†}={\mathbf{R}}^{(g,D)}[\mathbf{I}],$$(10)

where R^{(g, D)} is the response function consisting of the Fourier operator and Jones matrices.

It is important to note that field rotation angle matrices P are often precorrected during precalibration to stabilize the phase calibration in VLBI. However, the field rotation angle matrix P and the D-term matrix D are not commutative. Consequently, our model for precorrected field rotation angle data is $${\mathbf{V}}_{ij}^{\varphi \text{\hspace{0.17em}cor}.}={\mathbf{P}}_i^{-1}{\mathbf{G}}_i{\mathbf{D}}_i{\mathbf{P}}_i\mathbb{F}\mathbb{T}[\mathbf{I}]{\mathbf{P}}_j^{†}{\mathbf{D}}_j^{†}{\mathbf{G}}_j^{†}({\mathbf{P}}_j^{†}{)}^{-1}.$$(11)

This means we must reverse the precorrection of the field rotation angle in precalibrated data before adequately solving the Jones matrices, given the noncommutative properties of the matrices G, D, and P.

Assuming additive Gaussian noise on the radio interferometric data, justified by the central limit theorem, the likelihood given the noise covariance ℕ is expressed as $$\mathcal{P}(\mathbf{V}\mid \mathbf{I},\mathbf{G},\mathbf{D})=\mathcal{G}(\mathbf{V}-{\mathbf{R}}^{(g,D)}[\mathbf{I}],\mathbb{N}).$$(12)

In this work, the noise covariance is assumed to have a diagonal covariance; thus, the noise is uncorrelated. This likelihood formulation underscores how our method, unlike CLEAN-based polarization calibration methods, facilitates the simultaneous inference of calibration solutions and Stokes images. Consequently, the uncertainty estimation of Stokes images inherently reflects uncertainties from the calibration, data noise, and incomplete UV-coverage.

2.5 Polarization calibration and imaging prior model

For total intensity imaging, we utilize the resolve lognormal sky model, which builds on a Gaussian process in NIFTy software² $$I(\overrightarrow{x})=\exp (s(\overrightarrow{x})),s\curvearrowleft \mathcal{G}(s,S),$$(13)

where S is the covariance matrix of the Gaussian process.

A priori, we assume homogeneous and isotropic statistics for the sky brightness. Due to the Wiener-Khinchin theorem (Wiener 1949; Khinchin 1934), this assumption allows us to represent the covariance matrix S of the Gaussian process by a one-dimensional power spectrum P_s in the Fourier domain. As the degrees of freedom of the power spectrum scale linearly with the number of pixels in the image, this assumption makes it numerically affordable to infer spatial correlations even for high-dimensional (N > 10⁶) radio images. We transform the lognormal model so that the formal prior distribution is a standard normal distribution $$I(\overrightarrow{x})=\exp (s(\overrightarrow{x}))=\exp (\mathbb{F}[\sqrt{P_s({\xi }_s)}{\xi }_k]),\xi \curvearrowleft \mathcal{G}(\xi ,1),$$(14)

with P_s being our model for the power spectrum and ξ_s/k the standard normally distributed parameters. As in previous resolve applications, we use a stochastic process-based nonparametric model for the power spectrum P_s. Importantly, this nonparametric model allows coverage of a wide range of possible correlation patterns, making it robustly applicable to a wide range of sources.

In this work, this Gaussian process prior is utilized for polarimetric imaging and gain models. A more detailed description of our generative Gaussian process model can be found in Appendix B.

Generative model for polarization imaging prior I. The generative model for the Stokes sky emission I closely follows the model presented in Arras et al. (2025). Regarded as a complex 2 × 2 matrix, I is defined by $$\mathbf{I}=\left(\begin{array}{cc}⟨E_{i,r}{E}_{j,r}^{\ast }⟩& ⟨E_{i,r}{E}_{j,l}^{\ast }⟩\\ ⟨E_{i,l}{E}_{j,r}^{\ast }⟩& ⟨E_{i,l}{E}_{j,l}^{\ast }⟩\end{array}\right)=\left(\begin{array}{cc}I+V& Q+\mathbf{i}U\\ Q-\mathbf{i}U& I-V\end{array}\right),$$(15)

within the circular basis (Smirnov 2011). Here, E is the electric field and the indices i, j denote antenna labels and r, l refer to the respective circular feeds. The matrix must meet specific constraints: (1) I is positive definite and Hermitian, (2) strictly positive total flux I > 0, and (3) an upper bound on polarized emission I² ≥ Q² + U² + V² (Hamaker et al. 1996; Smirnov 2011).

Arras et al. (2025) identify the matrix exponential as a fitting parameterization for the polarized sky brightness distribution $$\mathbf{I}=e^x:=\exp \left(\begin{array}{cc}s+v& q+\mathbf{i}u\\ q-\mathbf{i}u& s-v\end{array}\right),$$(16)

where s, q, u, and v are real numbers for each pixel and, in particular, can be both positive and negative.

Stokes I, Q, U, and V can be obtained from 2D fields s, q, u, and v $$\begin{array}{rlrl}I& =e^s\cosh p,& Q& =\frac{q}{p}{e}^s\sinh p,\\ U& =\frac{u}{p}{e}^s\sinh p,& V& =\frac{v}{p}{e}^s\sinh p,\end{array}$$(17)

where $p:=\sqrt{q^2+u^2+v^2}$.

The total flux I is strictly positive in Equation (17) and I is Hermitian since the matrix exponential and the Hermitian conjugate commute. Furthermore, I is positive definite because the eigenvalues of a Hermitian matrix are real and the eigenvalues of the matrix exponential are the exponential of the eigenvalues. The determinant of I is the product of the eigenvalues, which is positive $$0<\det \mathbf{I}=I^2-Q^2-U^2-V^2.$$(18)

Thus, it proves that this parameterization satisfies all three conditions of I.

To complete the generative model for I, we need to define the models that generate the 2D fields s, q, u, and v. We model each of those fields using a Gaussian process (𝒢𝒫) with a nonparametric correlation kernel in Equation (14), captured through generative models that convert standard-normal distributed latent parameters ξ_s, ξ_q, ξ_u, ξ_v into the Gaussian process values s(ξ_s), q(ξ_q), u(ξ_u), v(ξ_v) (following the methods in Arras et al. (2021)). In our polarization imaging prior model, the spatial correlation between pixels and the correlation between Stokes images are taken into account since the correlation kernels of s, q, u, and v are inferred from the data.

We note that our polarization imaging model is not well suited for data showing larger amplitudes in RL, LR visibilities than RR, LL visibilities. However, those cases are limited, and our model enforces the polarization constraint only in the image domain, not in the visibility domain.

Generative model for gain prior g. The antenna-based gain prior consists of two Gaussian process models $$g(t)=\exp (\lambda (t)+\mathbf{i}\varphi (t)),$$(19)

where lognormal amplitude gain λ and phase gain ϕ are 1D Gaussian process priors.

For the 1D gain prior, we use the same generative Gaussian process model in the imaging prior (see Equation (14)). Thus, the lognormal amplitude gain λ(ξ_λ) and phase gain ϕ(ξ_ϕ) can be described by the parameters of the standard-normal distributed model ξ_λ and ξ_ϕ.

In radio interferometry the gain phase is inherently degenerate since we do not measure the absolute phase. To address the degeneracy, a standard approach is to choose a reference antenna with fixed gain phases in the polarization calibration. In our method, we can choose one reference antenna with high sensitivity or a short-baseline for each observation. Then we fix the RCP and LCP phase gains for the corresponding antenna to be zero Leppanen et al. (1995). We note that the phase gain prior with a range is analogous in part to setting up a reference antenna.

For millimeter-VLBI (mm-VLBI) data with short phase coherence time, we can use an uncorrelated normal-distribution ϕ ↶ 𝒩(0, ${\sigma }_{\varphi }^2$) for the phase gain prior. In resolve, the temporal correlation structures in gain solutions are inferred from the data. In other words, we determine the gain solution interval from the data without manual steering and we are even able to infer the solution interval per antenna, which can be advantageous for data from a heterogeneous VLBI array (Kim et al. 2025). More details on antenna-based gain calibration in resolve can be found in (Arras et al. 2019b; Kim et al. 2024).

Generative model for D-term prior D. The D-term (leakage) prior model is given by $$D=\exp (a+\mathbf{i}b),$$(20)

where a ↶ 𝒩(m_a, ${\sigma }_a^2$) and b ↶ 𝒩(0, ${\sigma }_b^2$) are normal distributed priors, m_a is the mean of the lognormal amplitude D-term, and σ_b and σ_b are the standard deviations of the lognormal amplitude and phase D-term, respectively.

Therefore, the D-term amplitude is lognormal distributed and the D-term phase is normal distributed. A summary table in Table 1 describes our calibration and imaging prior models.

2.6 Posterior distribution

The posterior distribution of all unknowns given the visibility data in Equation (8) is a very high dimensional object including the Stokes images, gains, and D-terms. Due to the high number of dimensions and the complicated relations between all the involved quantities, any representation of this probability function needs an approximation. One way of representing a posterior distribution is via a set of samples drawn from it. If s = (I, G, D) denotes the quantities to be inferred and d = V denotes the data, then the samples s_i ↩ 𝒫(s|d) with i ∈ {1, ... N} are an approximative representation of the posterior 𝒫(s|d), as any expectation with respect to some function f(s) can be calculated from those approximately $$⟨f(\mathbf{s}){⟩}_{(\mathbf{s}\mid \mathbf{d})}:=\int \mathcal{D}\mathbf{s}f(\mathbf{s})\mathcal{P}(\mathbf{s}\mid \mathbf{d})\approx \frac{1}{N}\sum _{i=1}^{N}f\left({\mathbf{s}}_i\right),$$(21)

where ∫ 𝒟s indicates a path integral.

In this work, we utilized two variational inference algorithms (MGVI, geoVI) (Knollmüller & Enßlin 2019; Frank et al. 2021) as implemented in the NIFTy software package (Selig et al. 2013; Arras et al. 2019a) to explore the posterior distribution of the Stokes images and calibration solutions. In the variational inference method, the posterior distribution is approximated as a parameterized distribution by minimizing the Kullback-Leibler (KL) divergence as a cost function. The KL divergence measures the information gain between two probability distributions.

These variational inference methods scale quasi-linearly in computational complexity with the problem size; therefore, they enable the solution of high-dimensional calibration and imaging problems. Specifically, geoVI, an extension of the MGVI algorithm, can describe non-Gaussian posteriors by constructing a coordinate transformation between the latent space, in which the prior was Gaussian, but the posterior is not, and another latent space, in which the posterior becomes approximate Gaussian. More details about MGVI and geoVI can be found in (Knollmüller & Enßlin 2019; Frank et al. 2021).

From the posterior samples, the posterior mean of the fractional linear polarization is $${⟨P_{\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}}(\mathbf{s})⟩}_{(\mathrm{s}\mid \mathbf{d})}:=\int \mathcal{D}\mathbf{s}\left(\frac{Q(\mathbf{s}{)}^2+U(\mathbf{s}{)}^2}{I(\mathbf{s}{)}^2}\right)\mathcal{P}(\mathbf{s}\mid \mathbf{d})=⟨\frac{Q^2+U^2}{I^2}⟩.$$(22)

We note that the posterior mean of the fractional linear polarization is not equal to (⟨Q⟩² + ⟨U⟩²)/⟨I⟩², where ⟨I⟩, ⟨Q⟩, and ⟨U⟩ are the posterior mean Stokes images, because of the non-linear dependence of P_frac on the Stokes parameters.

As another example, the posterior mean of the electric vector position angle (EVPA) is $$⟨\mathrm{\Phi }{⟩}_{(\mathbf{s}\mid \mathbf{d})}:=\int \mathcal{D}\mathbf{s}\frac{1}{2}\arctan \left(\frac{U(\mathbf{s})}{Q(\mathbf{s})}\right)\mathcal{P}(\mathbf{s}\mid \mathbf{d})=⟨\frac{1}{2}\arctan \left(\frac{U}{Q}\right)⟩.$$(23)

Fig. 1 Comparison between the ground truth image (left panel) convolved with the nominal CLEAN beam with the uniform weighting and the posterior mean resolve polarization reconstruction (right panel) with EVPAs in all images with colors corresponding to the fractional linear polarization Pfrac. The contours represent the total intensity of corresponding images.

3 Application to synthetic data

In VLBI, polarization calibration using multiple calibrators at different declinations helps to reconstruct more robust leakage solutions, which break the degeneracy between the field rotation angle matrix P and the leakage matrix D (Park et al. 2021b; Martí-Vidal et al. 2021). In order to validate the Bayesian polarization calibration method with multiple calibrators, we tested our method with three of the synthetic VLBA datasets presented in Park et al. (2021b). Synthetic datasets were produced using PolSimulate in the CASA software (McMullin et al. 2007) with OJ287, 3C273, BLLac UV-coverage at 15GHz by ten VLBA antennas. More details on synthetic data can be found in Park et al. (2021b). The ground truth image consists of three point sources. The ground truth image convolved with the nominal CLEAN beam with uniform weighting is in Fig. 1. We assume that there is no gain corruption in the data.

Polarization calibration with synthetic data was performed as follows. First, initial D-term estimates were obtained using the maximum a posteriori (MAP) method using all three calibrator datasets. Then, the posterior distribution of Stokes images and D-term was reconstructed with 3C273 synthetic data using the geoVI method (Frank et al. 2021) starting from the estimated MAP D-terms as the initial condition. Fig. 1 shows the resolve total intensity image and EVPAs with colors corresponding to the fractional linear polarization (right panel). We chose a spatial domain of 256 × 256 pixels and a field of view of 10 mas × 10mas. The resolve image was able to recover three linearly polarized components (with fractional linear polarization 2%, 5%, and 11% respectively).

Figure 2 represents a comparison between the resolve D-term posterior distribution and the ground truth D-terms. The resolve posterior distribution is obtained from 100 posterior samples using Gaussian kernel density distribution in Scipy Python library (Virtanen et al. 2020). Overall, the reconstructed D-term posterior distributions using resolve are consistent with the ground truth D-terms within the 2σ errors. This result demonstrates that the resolve polarization calibration and imaging method is able to reconstruct reliable D-term solutions and Stokes images using multiple calibrator datasets.

Fig. 2 Comparison between the D-term posterior using resolve and ground truth D-terms from the synthetic data. Contours show 1σ and 2σ cumulative regions of resolve posterior D-terms using Gaussian kernel density estimation. The plus signs correspond to the ground truth D-terms.

4 Application to real data

4.1 3C273 VLBA observation at 15 GHz

We applied our Bayesian polarization calibration and imaging method to precalibrated (without self-calibration and D-term calibration) 3C273 VLBA MOJAVE survey data at 15 GHz on January 28, 2017 (Lister et al. 2018) to demonstrate that resolve is able to infer D-terms per intermediate frequency (IF) from a source with complex structure. The data have eight IFs (32 MHz each) and the total bandwidth is 256 MHz.

We reconstructed the resolve Stokes images with a spatial domain of 256 × 256 and a field of view of 50 mas × 50mas. For the resolve reconstruction, we added a 10% systematic error budget in the data and selected Los Alamos (LA) as a reference antenna for the polarization calibration. The reduced χ² of the resolve reconstruction was 1.45, and the number of posterior samples was 20. The wall-clock time for the resolve reconstruction was around 7 hours on a laptop with five message passing interface (MPI) tasks.

We assumed that Stokes images do not vary over frequencies due to the narrow bandwidth (256 MHz), and performed antenna gain self-calibration and leakage calibration per IF. We assumed that each antenna has four different correlation kernels for the log-amplitude gain and phase gain and two polarization modes. The same correlation kernel is inferred over multi-IFs. The model parameters for the polarization imaging prior, gain prior, and D-term prior for the 3C273 VLBA data can be found in Tables B.1 and B.2.

Figure 3 shows a comparison of linear polarization reconstructions from CLEAN and the resolve posterior mean. The EVPA colors correspond to the fractional linear polarization, and contours represent the total intensity of the corresponding images. The CLEAN images are taken from the MOJAVE archive³. For the CLEAN reconstruction, the MOJAVE team performed self-calibration and image reconstruction using the DIFMAP software (Shepherd 1997) and the LPCAL task (Leppanen et al. 1995) in the AIPS software (Greisen 2003) for polarization calibration. The D-term solutions using LPCAL are the median solution values of all sources in the epoch after removing obvious outliers (Lister & Homan 2005).

In the core region, the resolve image (right panel) has a higher resolution compared with the core in the CLEAN image (left panel). We note that it is also possible to use an over-resolved CLEAN beam to obtain a higher-resolution restored image. The extended total intensity emission in the resolve image looks thinner than the CLEAN image. The jet emissions in the resolve image are overall consistent with the CLEAN reconstruction with bright jet emissions (see the contour from 1.2% of the peak total intensity).

The linear polarization comparison between CLEAN and resolve shows noticeable differences. The fractional linear polarization on the edges of linearly polarized emission in the CLEAN reconstruction is relatively higher than that in the resolve reconstruction. The discrepancy may result from the lack of a polarization constraint in the CLEAN imaging prior or biases from the similarity approximation in CLEAN-based polarization calibration. In contrast to CLEAN, resolve can describe complex source structure spanning a range of spatial scales using the Gaussian process sky prior we described, which is aware of spatial correlations of the polarized flux while ensuring the polarization constraints at each image pixel individually. Encoding the polarization constraint facilitates physically sensible reconstruction consistent with the theoretical fractional linear polarization limit (up to 75%) for optically thin synchrotron radiation from non-thermal electrons with a power-law distribution of energies (Rybicki & Lightman 1979).

The EVPA pattern is generally consistent except for the core region, due to the improved resolution in the resolve image. Both images show the rotational EVPA structure at the core, but the resolve image exhibits thinner emission structures that show more abrupt changes. Figure D.1 shows the 3C273 resolve EVPA standard deviation. We note that a robust rotation measure can be estimated using the Bayesian approach (Vogt & Enßlin 2005). A detailed investigation of Bayesian rotation measure analysis using our method is left for future work. In the resolve reconstruction, Stokes V emission is negligible (less than 0.2% of the total intensity emission). The Stokes V reconstruction requires additional calibration steps. This aspect will be addressed in future work.

For the leakage calibration, a total of 160 D-terms (two polarization modes × ten antennas × eight IFs) are inferred since the D-terms in MOJAVE observations with VLBA at 15 GHz tend to be different per IF. This approach of reconstructing D-terms per IF is highly desirable, especially for wideband observations, as the signal path for each IF is different, which can result in changes to the D-terms for each IF separately. We also note that baseband boundaries can cause jumps in leakage that will also differ for individual IFs (Martí-Vidal et al. 2021).

In Figure C.1, the resolve D-term posterior means and the D-terms obtained using the LPCAL software are shown to be consistent with each other, as they are strongly correlated. It is important to note that the resolve D-term solutions are from 3C273 data only, and the D-terms using LPCAL are the median values of D-term solutions from multiple sources. The results indicate that our Bayesian polarization calibration method is able to estimate reliable D-terms from data with complex source structures. A more detailed analysis with additional datasets is deferred to future work.

Fig. 3 Comparison between the 3 C 273 VLBA CLEAN and resolve posterior mean linear polarization reconstructions at 15 GHz. Colored ticks indicate EVPAs in all images, with colors corresponding to the fractional linear polarization Pfrac. The contours representing the total intensity of corresponding images increase by a factor of 4, starting from 0.3% of the peak total intensity of corresponding image.

4.2 OJ287 GMVA+ALMA observation at 86 GHz

Zhao et al. (2022) reported the GMVA+ALMA observation of the blazar OJ 287 at 86 GHz on April 2, 2017. The blazar OJ 287 is a supermassive binary black hole candidate that shows quasiperiodic optical outbursts with a period of about 12 years (Sillanpaa et al. 1988). GMVA observations jointly with the phased ALMA provide high fidelity images due to improved sensitivity and long north-south baselines (Issaoun et al. 2019; Zhao et al. 2022; Lu et al. 2023; Kim et al. 2025). However, polarization calibration and imaging for mm-VLBI observations are demanding due to the low S/N of the polarized signals, tropospheric phase corruption, and heterogeneous antenna statistics, such as in GMVA+ALMA data.

We revisit the GMVA+ALMA observation of the blazar OJ 287 at 86 GHz in Zhao et al. (2022) to validate our Bayesian polarization calibration and imaging method using precalibrated data. For polarization calibration, robust gain self-calibration is required due to the degeneracy between D-terms and gains. Kim et al. (2024, 2025) validate the resolve self-calibration methods with mm-VLBI datasets. Reconstructed antenna gain solutions from the VLBA M87 observation in 2013 at 43 GHz indicate that gain solutions from a homogeneous array with high S/N are consistent among CLEAN, ehtim, and resolve methods (Kim et al. 2024). However, CLEAN self-calibration utilizes a crude regularizer (e.g., a uniform solution interval for different antenna gain solutions) and often flags a significant fraction of the data. Therefore, these limitations may hinder robust self-calibration for mm-VLBI observations.

To mitigate this issue, we employed different log-amplitude gain correlation kernels per antenna in order to take into account heterogeneous array sensitivity. For the gain phase prior, we utilized uncorrelated normally distributed phases due to the short phase coherence time in GMVA observations at 86 GHz (Kim et al. 2025). More specifically, we did not model the temporal correlation in phase solutions, which is analogous to setting up phase solution intervals smaller than the averaged time for all antennas. Employing a different gain amplitude kernel per antenna is analogous to setting up different solution interval constraints per antenna in self-calibration. In resolve, the time correlations in the gain solutions are inferred from the data in an automated fashion, and the large uncertainties of particular data points are taken into account naturally in the Bayesian framework. Therefore, we are able to perform more robust self-calibration without flagging a significant number of data points compared to CLEAN-based self-calibration methods. Furthermore, the large gain uncertainties inherent to mm-VLBI observations are explicitly accounted for and propagated into the image domain. As a result, the reliability of the calibration and imaging is quantified through uncertainty estimation.

Polarization calibration for mm-VLBI datasets presents several additional challenges. In mm-VLBI, we often use the target data for D-term calibration due to a lack of point-like calibrators. Polarization leakages tend to be higher in mm-VLBI than in centimeter-VLBI, and the high resolution often reveals complex source polarization structures. Thus, precise D-term calibration is crucial to obtain reliable polarimetric images, and accurate polarimetric images, in turn, improve D-term calibration. However, in conventional polarization calibration methods, utilizing the similarity approximation or CLEAN images as a prior hinders the reconstruction of complex polarization structures since a few subcomponents in the similarity approximation and the collection of delta components in the CLEAN reconstruction cannot accurately represent the extended emission and complex source structures that are smaller than the CLEAN beam. On the other hand, in resolve, polarization imaging models with complex source structures, satisfying the polarization constraint and consistent with the data, are utilized for more robust polarization calibration.

The data are time-averaged with 15 seconds and frequency-averaged for ehtim and resolve software, since the antenna leakages in the data are similar over frequencies (Zhao et al. 2022). For the CLEAN image, self-calibration and image reconstruction were performed using the DIFMAP software and the LPCAL method of the AIPS software was used for leakage calibration using four individual intermediate frequencies (IFs). For the ehtim reconstruction, self-calibration and polarization calibration were performed in an iterative fashion using the ehtim software. The details regarding the CLEAN and ehtim polarization calibration and image reconstructions can be found in Zhao et al. (2022).

For the resolve reconstruction, Stokes I, Q, U, and V images with a spatial domain of 256 × 256 pixels and a field of view of 500 μas × 500 μas, as well as gain solutions with a time interval of 15 seconds and leakage solutions, are inferred simultaneously. The reduced χ² of the resolve reconstruction was 1.2, and the number of posterior samples was 100. The wall-clock time for the resolve reconstruction was around 8.5 hours on a single node of the MPIfR cluster with 25 MPI tasks. The model parameters for the polarization imaging prior, gain prior, and D-term prior for the OJ287 GMVA+ALMA data can be found in Tables B.3 and B.4.

Figure 4 depicts the total intensity with EVPAs of the blazar OJ287 at 86 GHz using three different imaging algorithms: CLEAN, ehtim from Zhao et al. (2022), and resolve. All three image reconstructions show three main components in total intensity and diverging EVPAs in the northwest component. Furthermore, the curved jet (middle component) from the core (southeast component) is recognizable only in the ehtim and resolve images. In contrast, CLEAN has limitations in reconstructing smaller scale structures than the CLEAN beam. The ehtim and resolve images show better resolution than the CLEAN image because forward modeling permits partial removal of the observational point spread function during the data inversion process. The ehtim image exhibits higher resolution than the resolve image, since the sparsity promoting regularizers in ehtim tend to produce extremely sharp structures.

For polarization imaging, polarization constraints are encoded in resolve and ehtim. We performed the absolute EVPA calibration using the reported integrated EVPA of ALMA-only OJ287 image at 86 GHz in Goddi et al. (2021). The details of the absolute EVPA calibration method can be found in Section 4.3. In Figure 4, the EVPAs in the northwest jet component diverge in all three images. The fractional linear polarization along the northeast direction in the diverging jet component increases gradually in the images from CLEAN and resolve, whereas the ehtim image shows a nearly monotonic fractional linear polarization. In the ehtim image, there are extremely highly polarized pixels ($P_{\text{frac\hspace{0.17em}}}>99\mathrm{\%}$) that may be imaging artifacts caused by the sparsity regularizers.

The southeast core component and the middle jet component show EVPAs along the jet direction. This indicates toroidal magnetic fields, given the relatively small Faraday rotation (3.05 ± 0.62 × 10³ rad m⁻²) from the ALMA observation (Goddi et al. 2021). In the ehtim and resolve images, we see bending EVPA patterns, and the pattern is more wiggling in the resolve image. The diverging EVPA pattern in the northwest jet component can be explained by oblique or recollimation shocks (Zhao et al. 2022). Figure D.2 shows the OJ287 resolve EVPA posterior standard deviation. In the resolve reconstruction, the OJ287 Stokes V emission is negligible (less than 0.1% of the total intensity emission).

Figure 5 shows the amplitude gain posterior distribution. We employed separate correlation kernels per antenna to infer different temporal correlation structures arising from the heterogeneous antenna sensitivities. We utilized an uncorrelated phase gain prior since the phase coherence time is comparable to the averaging time (15 seconds).

Figure 6 depicts a comparison between the resolve D-term posterior distributions and the reconstructed D-terms using ehtim. The resolve and ehtim D-terms are broadly consistent. In Figure C.2, cross-correlation plots between the resolve posterior D-terms and ehtim reconstructed D-terms are shown. They exhibit strong positive correlations, but some antennas (e.g., PT) show weaker positive correlations. This discrepancy between resolve and ehtim results from the highly corrupted phases at 3 mm and the different self-calibration and flagging routines used. For polarization calibration and imaging, ehtim used eight antennas out of 13, whereas resolve utilized all antennas (YS observed only a single polarization; therefore, we show the resolve D-terms for 12 antennas). In resolve, non-Gaussian leakage posterior distributions were obtained using the geoVI method. A few antennas (e.g., BR and ON) exhibit higher uncertainty in the D-term phase due to highly corrupted visibility phases. The ehtim D-terms are overall consistent with D-terms obtained from the LPCAL software (Zhao et al. 2022). In conclusion, the D-term solutions are in good agreement across the resolve, ehtim, and LPCAL methods.

We note that regularized maximum likelihood (RML) methods can be interpreted as MAP estimation in Bayesian language (Kim et al. 2024). However, the point estimation cannot represent calibration solutions with high uncertainties well, and the MAP estimation is prone to overfitting the data. Furthermore, ehtim does not support multiscale base functions. In contrast to resolve, it tends to produce extremely sharp features. Different correlation structures between the core and the jet cannot be adequately described due to sparsity promoting priors, such as the l₁-norm, total variation (TV), and total squared variation (TSV). Employing extremely sharp images as a model in polarization calibration may result in biased D-terms and polarization images. A more detailed discussion of the Bayesian interpretation of regularizers in the RML method can be found in Appendix B of Kim et al. (2024).

Fig. 4 Comparison among the OJ287 GMVA+ALMA CLEAN, ehtim, and resolve posterior mean linear polarization reconstructions at 86 GHz. Colored ticks indicate EVPAs in all images, with colors corresponding to the fractional linear polarization Pfrac. The contours representing the total intensity of the corresponding images increase by a factor of 2, starting from 10% of the peak resolve posterior mean total intensity. All images were processed by Gaussian interpolation.

Fig. 5 Posterior amplitude gains of the OJ287 GMVA+ALMA observation using resolve. The solid line represents the amplitude gain posterior mean with a semi-transparent standard deviation. Each row represents an individual antenna with corresponding abbreviated name in the bottom left corner of each RCP plot.

Fig. 6 Comparison between OJ287 GMVA+ALMA D-term posterior distributions using resolve and D-term solutions using ehtim. Contours show 1σ and 2σ cumulative regions of resolve posterior D-terms using Gaussian kernel density estimation. The plus signs correspond to the reconstructed D-terms using ehtim.

4.3 Absolute EVPA calibration

Most VLBI observations use circular polarization feeds. Therefore, we need to constrain the arbitrary phase difference between the RCP and LCP feeds at each frequency in the fully station-based approach to the polarization classification and calibration (Cotton 1993; Leppanen et al. 1995). This process is called absolute EVPA calibration. The simplest approach to determining this phase offset is to use a trusted external observation as an “anchor”. External anchors are typically observations with a single dish, an interferometer with linear polarization feeds, or an observation that has already calibrated the phase difference between RCP and LCP independently. This method relies on the assumption that the signal of a “target” observation is similar to that of the anchoring observation, in both time and the physical structure of the emission. As such, the best calibrated sources for absolute EVPA calibration have the following three characteristics:

1. The source has sufficient linear polarization in both the target and anchor observation.

2. The source is not highly variable in time; as such it does not vary between the two observations.

3. The source is compact; this is particularly relevant when comparing a target VLBI observation against single-dish or unresolved non-VLBI observations as we assume the emission from both observations comes from the same physical structure of the source.

To compute the RCP-LCP phase correction, we take twice the EVPA offset between the target and anchor observations. Typically, the EVPA offset is found by calculating the integrated Stokes Q and U flux density of both observations and calculating an integrated EVPA as EVPA = 0.5 arctan[U/Q]. This can be done by only considering Stokes Q and U emission that is co-spatial with Stokes I emission to ensure a direct comparison between the target and anchor.

In the MOJAVE survey, the absolute EVPA calibration relied on optically thin jet features in multiple sources that have relatively stable EVPAs. Based on the comparisons of compact AGNs with near-simultaneous single-dish observations, the MOJAVE team estimated that the VLBA EVPA measurements are accurate to ~5% (Lister et al. 2018, 2021). For our reconstruction of 3C273 in Figure 3, we applied the same correction value from the MOJAVE team for the CLEAN reconstruction. For our reconstruction of OJ287 in Figure 4, we calculated the integrated EVPA of our reconstructed resolve image of Stokes Q and U and compared it with ALMA observations with the EVPA − 69.69° at 86.3 GHz in Goddi et al. (2021). From this we found an EVPA correction of −16.07° based on the OJ287 resolve-integrated EVPA and ALMA-only EVPA.

In short, the resolve reconstructions of 3C273 at 15 GHz OJ 287 at 86 GHz demonstrate that our Bayesian polarization calibration and imaging pipeline can obtain reliable super-resolution polarimetric images, D-terms, and gain solutions from VLBI data. Reconstructed Stokes images are utilized as a polarization calibration model, enabling us to obtain more complex polarization structures than conventional CLEAN-based methods. Uncertainty estimation of reconstructed images and calibration solutions in resolve is especially beneficial in mm-VLBI due to the high calibration uncertainties from tropospheric phase corruptions and antenna leakages. In our method, calibration uncertainties are accounted for in our reconstructed images, and we are able to quantify the reliability of Stokes images and calibration solutions. The resolve reconstruction FITS files and results are available in the Zenodo archive.⁴

5 Conclusions

In this work, we have presented a novel Bayesian polarization calibration and imaging method using the resolve algorithm. Our method can simultaneously reconstruct high-resolution polarimetric images that satisfy the polarization constraint ($I\ge \sqrt{Q^2+U^2+V^2}$), antenna-based gains, and polarization leakages from the entire dataset. The polarization calibration is based on the reconstructed Stokes images instead of the similarity approximation used in traditional approaches, which does not hold for complex polarized emission patterns. Therefore, we can reconstruct reliable polarimetric images with complex structures and more robust calibration solutions, consistent with the entire dataset. Moreover, multisource and multi-IF polarization calibration is supported due to the scalability of resolve, achieved through variational inference methods (Knollmüller & Enßlin 2019; Frank et al. 2021).

We have demonstrated our method with synthetic and real observations. Examples using the quasar 3C273 MOJAVE VLBA data at 15 GHz and the blazar OJ287 GMVA+ALMA observation at 86 GHz show that resolve can reconstruct polarimetric images with super-resolution (Honma et al. 2014), revealing structures below the synthesized beam, which is the resolution limit of conventional CLEAN-based methods. The fractional linear polarization in resolve images is consistent with the theoretical estimation of synchrotron radiation (<75%) in AGN jets.

Reconstructed D-terms obtained with resolve are consistent with those from the conventional CLEAN-based method and the RML-based method ehtim. Furthermore, our RIME-based Bayesian polarization calibration and imaging method is able to reconstruct D-term solutions from target data with complex source structure without using calibrators.

Validation with synthetic and real datasets indicates that incorporating polarization calibration into imaging is beneficial for sparse and noisy VLBI datasets in order to estimate reliable and reproducible calibration solutions and images. Calibration uncertainties are explicitly taken into account in our final resolve results, and the reliability of the gain and leakage solutions was estimated. For future work, the EVPA posterior distribution can be utilized for rotation measure analysis, and frequency-dependent D-term calibration within the Bayesian framework should be investigated for wideband radio interferometric observations.

Data availability

The pipeline is publicly accessible at https://github.com/JongseoKim/resolve_polcal.

Acknowledgements

We thank the anonymous referee for constructive comments and suggestions, Yuri Kovalev for comments on the manuscript, Guang-Yao Zhao and Jose Gomez for providing the GMVA+ALMA OJ287 data, MOJAVE team for providing the VLBA 3C273 data. J.K. received financial support for this research from the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. This work was supported by the M2FINDERS project funded by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement No. 101018682). J.R. acknowledges financial support from the German Federal Ministry of Education and Research (BMBF) under grant 05A23WO1 (Verbundprojekt D-MeerKAT III). This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope, the Green Bank Observatory, and the Very Long Baseline Array (VLBA). The VLBA and the GBT are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The data were correlated at the correlator of the MPIfR in Bonn, Germany. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team Lister et al. (2018).

References

Appendix A Comparison of VLBA 3C273 resolve reconstructions

Figure A.1 represents the comparison between VLBA 3C273 Stokes I resolve posterior mean reconstructions using Bayesian self-calibration and imaging method (left figure) and Bayesian polarization calibration and imaging method (right figure). We note that VLBI data can have higher LR, RL amplitudes than RR, LL amplitudes in some baselines. Even if we do not enforce the polarization constraint in the visibility domain, encoding the polarization constraint in the image domain might generate spurious Stokes I emission. To validate this issue, we performed Bayesian self-calibration and imaging method (Kim et al. 2024) using VLBA 3C273 multi-IF data. Two images look consistent and the Pearson correlation coefficient between two images is 0.96. In conclusion, the polarization constraint does not produce spurious Stokes I emission for the 3C273 data. However, the Stokes I image from the polarization calibration and imaging method (right panel) has slightly more artifacts outside of the emission region. It might result from the sparse UV-coverage (The on-source time is around 1 hour), which complicates the solution of the degenerate inverse problem.

Fig. A.1 Comparison between the 3C273 VLBA resolve posterior mean Stokes I reconstructions at 15 GHz using self-calibration and polarization calibration. The colorbar starts from 0.1% of the peak total intensity of resolve image using the polarization calibration.

Appendix B Hyperparameter setup for the image and gain prior model

In the Gaussian process prior model in Equation 14, the offset mean represents the logarithmic mean flux of the field in the unit of Jy/str. The zero mode variance mean and zero mode variance std correspond to the standard deviation of the offset mean and standard deviation of the standard deviation of the offset mean respectively. The fluctuation activates the fluctuation of the field, related to the dynamic range of the image. The flexibility and asperity are related to stochastic processes, such as the integrated Wiener process and Wiener process in the model. The flexibility parameter ensures flexibility of the power spectrum shape and asperity activates certain amplitude mode in the power spectrum, which can generate periodic patterns in the field. Average slope represents the slope of the power spectrum. As an example, a steep power spectrum corresponds to a smooth image since high Fourier modes are suppressed. The exact mathematical definition of the generative Gaussian process model can be found in Arras et al. (2021, Sec. 3.4).

Table B.1 and Table B.2 show the hyperparameter setups for the VLBA 3C273 at 15 GHz resolve image prior and calibration prior respectively. The hyperparameter setups for the GMVA+ALMA OJ287 at 86 GHz resolve image prior and calibration prior are in Table B.3 and Table B.4 correspondingly.

Appendix C D-term cross correlation plots

Fig. C.1 Cross correlation plots between 3C273 resolve posterior mean D-terms and LPCAL D-terms.

Fig. C.2 Cross correlation plots between OJ287 GMVA+ALMA D-term posterior means using resolve and D-term solutions using ehtim.

Appendix D EVPA posterior standard deviation plots

Fig. D.1 EVPA posterior standard deviation map of the 3C273 observation at 15 GHz using resolve. The contours representing the total intensity resolve posterior mean image increase by a factor of 4, starting from 0.3% of the peak resolve total intensity.

Fig. D.2 EVPA posterior standard deviation map of the OJ287 GMVA+ALMA observation at 86GHz using resolve. The contours representing the total intensity resolve posterior mean image increase by a factor of 2, starting from 10% of the peak resolve posterior mean total intensity.

All Tables

All Figures

	Fig. 1 Comparison between the ground truth image (left panel) convolved with the nominal CLEAN beam with the uniform weighting and the posterior mean resolve polarization reconstruction (right panel) with EVPAs in all images with colors corresponding to the fractional linear polarization Pfrac. The contours represent the total intensity of corresponding images.
In the text

	Fig. 2 Comparison between the D-term posterior using resolve and ground truth D-terms from the synthetic data. Contours show 1σ and 2σ cumulative regions of resolve posterior D-terms using Gaussian kernel density estimation. The plus signs correspond to the ground truth D-terms.
In the text

Fig. 3 Comparison between the 3 C 273 VLBA CLEAN and resolve posterior mean linear polarization reconstructions at 15 GHz. Colored ticks indicate EVPAs in all images, with colors corresponding to the fractional linear polarization Pfrac. The contours representing the total intensity of corresponding images increase by a factor of 4, starting from 0.3% of the peak total intensity of corresponding image.

In the text

Fig. 4 Comparison among the OJ287 GMVA+ALMA CLEAN, ehtim, and resolve posterior mean linear polarization reconstructions at 86 GHz. Colored ticks indicate EVPAs in all images, with colors corresponding to the fractional linear polarization Pfrac. The contours representing the total intensity of the corresponding images increase by a factor of 2, starting from 10% of the peak resolve posterior mean total intensity. All images were processed by Gaussian interpolation.

In the text

	Fig. 5 Posterior amplitude gains of the OJ287 GMVA+ALMA observation using resolve. The solid line represents the amplitude gain posterior mean with a semi-transparent standard deviation. Each row represents an individual antenna with corresponding abbreviated name in the bottom left corner of each RCP plot.
In the text

	Fig. 6 Comparison between OJ287 GMVA+ALMA D-term posterior distributions using resolve and D-term solutions using ehtim. Contours show 1σ and 2σ cumulative regions of resolve posterior D-terms using Gaussian kernel density estimation. The plus signs correspond to the reconstructed D-terms using ehtim.
In the text

	Fig. A.1 Comparison between the 3C273 VLBA resolve posterior mean Stokes I reconstructions at 15 GHz using self-calibration and polarization calibration. The colorbar starts from 0.1% of the peak total intensity of resolve image using the polarization calibration.
In the text

	Fig. C.1 Cross correlation plots between 3C273 resolve posterior mean D-terms and LPCAL D-terms.
In the text

	Fig. C.2 Cross correlation plots between OJ287 GMVA+ALMA D-term posterior means using resolve and D-term solutions using ehtim.
In the text

	Fig. D.1 EVPA posterior standard deviation map of the 3C273 observation at 15 GHz using resolve. The contours representing the total intensity resolve posterior mean image increase by a factor of 4, starting from 0.3% of the peak resolve total intensity.
In the text

	Fig. D.2 EVPA posterior standard deviation map of the OJ287 GMVA+ALMA observation at 86GHz using resolve. The contours representing the total intensity resolve posterior mean image increase by a factor of 2, starting from 10% of the peak resolve posterior mean total intensity.
In the text