© ESO, 2011

1. Introduction

SiO maser emission acts as an important tracer of the astrophysical conditions in the near-circumstellar environments of late-type, evolved stars (Habing 1996). The maser emission is detected in vibrationally-excited, rotational molecular transitions of SiO, denoted ν = ν′,  J = J′ − (J′ − 1), which have rest frequencies of order ν₀ ~ 43J′ GHz. Commonly, the brightest and most frequently observed lines are the ν = 1,  J = 1−0 and ν = 1,  J = 2−1 transitions near 43 GHz and 86 GHz respectively.

In single-dish observations, the SiO maser emission from late-type, evolved stars shows an appreciable degree of linear polarization: m_l ~ 15−30%, or greater (Troland et al. 1979), but a significantly lower degree of circular polarization, with typical median values m_c ~  several percent (Barvainis et al. 1987). In a single-dish survey of SiO linear polarization variability in evolved stars, Glenn et al. (2003) report a mean sample linear polarization with a dispersion of 7%. A larger single-dish survey of late-type, evolved stars in full polarization by Herpin et al. (2006) was able to classify mean polarization by stellar type or variability class. They report and for Mira variables, and for semi-regular variables, and and for supergiant stars. Overall, the measured single-dish values are consistent with theoretical expectations, given the non-paramagnetic nature of the SiO molecule and the associated general theory of the propagation of polarized maser emission in the limit of small Zeeman splitting (Goldreich et al. 1973).

Radio-interferometric observations of SiO maser emission toward late-type, evolved stars can measure the radiation properties of the SiO maser emission at high angular resolution in the full set of Stokes parameters  {I,Q,U,V}. For Very Long Baseline Interferometry (VLBI), this angular resolution is at the sub-milliarcsecond (mas) level, and is sufficient to resolve individual SiO maser components. However, the low level of mean circular SiO maser polarization, and the intrinsic millimeter observing wavelengths required for these transitions, pose challenges in Stokes V instrumental calibration of VLBI observations of SiO masers. Measurements of SiO maser circular polarization at mas angular resolution are important however, both to constrain theoretical models describing the propagation of polarized maser emission, and for the subsequent application of these theories to infer the magnitude and orientation of the underlying magnetic fields from measurements of component-level SiO maser polarization properties. Recently, for example, periodic changes in single-dish fractional circular polarization have been detected toward two SiO stars (Wiesemeyer et al. 2009), with a proposed explanation that these changes are tracing precessing planetary magnetospheres in the circumstellar environment (Wiesemeyer 2009).

In this paper we consider the instrumental factors that constrain measurements of SiO circular polarization using VLBI arrays, evaluate the degree to which these effects can be calibrated and corrected, and provide an assessment of the net accuracy of SiO circular polarization observations at 43 GHz and 86 GHz.

The development and first application of VLBI polarimetry was reported by Cotton et al. (1984) and Wardle et al. (1986). It has since grown to become a mature technique, of broad applicability; see Cawthorne et al. (1993) and Pollack et al. (2003) for example. For the frequency bands considered in the current paper, the first continuum polarization VLBI were reported at 43 GHz by Kemball et al. (1996), and at 86 GHz by Attridge (2001).

These prior works were concerned with linear polarimetry of continuum extra-galactic radio sources. Homan & Wardle (1999) developed a method to allow circular VLBI polarimetry of continuum sources however, by assuming a source ensemble has ; this technique has been demonstrated successfully at wavelengths longer than λ ≥ 1 cm (Homan et al. 2009), and is discussed in further detail below. The first spectral-line polarimetry in full Stokes  {I,Q,U,V} at 43 GHz – of the ν = 1,J = 1−0 SiO maser emission toward the Mira variable, TX Cam – was reported by Kemball & Diamond (1997), based on technical development described earlier by Kemball et al. (1995). Critical circular spectral-line VLBI polarimetry of water masers in the adjacent 22 GHz band was succesfully performed by Vlemmings et al. (2002) following this same method.

The current paper examines the sources of instrumental and propagation errors in VLBI circular polarimetry in the λ = 7 mm and λ = 3 mm bands in more detail, and develops and tests refinements to existing calibration methods to improve the accuracy to which Stokes V and m_c can be measured in such observations. We show that it is possible to refine existing calibration methods based on autocorrelation data that avoid the need for amplitude self-calibration. We demonstrate these refinements on an observed dataset, and show that an error in m_c of  ≤ 0.5% or better is possible in practice at λ = 7 mm and m_c ≤ 0.5−1% or better at λ = 3 mm.

The paper is structured as follows: Sect. 2 describes the instrumental and propagation data model, and Sect. 3 discusses the associated data reduction methods. The analysis of a test dataset is described in Sect. 4. General discussion and conclusions are presented in Sects. 5 and 6 respectively.

2. Theory

This section describes the data model, containing propagation and instrumental effects relevant to precise amplitude calibration of millimeter-wavelength spectral-line VLBI observations.

Given the focus of the paper described in the Introduction, this discussion is confined primarily to the case of the VLBA¹ operating at observing frequencies of 43 GHz (λ = 7 mm) and 86 GHz (λ = 3 mm). However, these techniques are broadly applicable to spectral-line VLBI observations at millimeter wavelengths in general. Implicit to this discussion are the instrumental parameters applicable to the VLBA; these are summarized in Table 1, individually referenced as indicated from Romney (2010) or Napier (1995). For each of the two observing bands considered here, the table contains columns summarizing the receiver band frequency limits, adopted center frequency, typical zenith system equivalent flux density (SEFD), primary beam angular full-width half-maximum θ_b, typical root-mean-square (rms) pointing error σ_p, beam squint separation  △ _s, and the typical baseband bandwidth assigned to an individual maser transition at these observing frequencies. The beam squint results from the offset Cassegrain optics adopted in the VLBA antenna design (Napier 1995). We note also that each of these VLBA bands is dual polarization, with nominal orthogonal RCP and LCP polarization receptors.

2.1. Data model

The data model includes all propagation and instrumental calibration effects in the signal path from the frame of the astronomical source to the correlator output. Kemball et al. (1995) present a signal-path analysis of this type for spectral line polarization VLBI in general. We extend this formalism here to include a more detailed assessment of the problem of amplitude calibration for millimeter observing wavelengths in the case of low, but non-zero, Stokes V. For brevity, we use the same notation as the original paper by Kemball et al. (1995) wherever possible, and enumerate only the data model differences here. These are presented in the same approximate sequential order in which they appear in the signal path.

Kemball et al. (1995) presented a reduction heuristic for spectral line polarization VLBI observations for the case where no prior assumption can be made concerning the magnitude of Stokes V. In contrast for example, in continuum VLBI polarimetry of compact, extra-galactic sources, the approximation Stokes V ~ 0 can be applied during data reduction to a relatively high degree of accuracy. The approach described by Kemball et al. (1995) defines a data model and develops a full calibration solution in amplitude, phase, group delay, and fringe rate for a reference receptor polarization p ∈  {R,L} using standard single-polarization VLBI techniques; and then defines a method to tie that solution to the calibration solution for the orthogonal receptor polarization q by solving for and applying differential polarization offsets in amplitude, phase, and group delay. We retain that methodology here, but focus specifically on total and differential polarization amplitude calibration, and its role in precise Stokes V measurement at millimeter wavelengths. This discussion includes a deeper examination of the data model for autocorrelation data specifically, given their relevance to amplitude calibration for spectral line VLBI observations as approached here.

2.2. General

We carry over several general assumptions made in Kemball et al. (1995) to the current analysis. Given the small angular size of astrophysical maser regions observed with VLBI and the millimeter observing wavelengths considered here, the interferometric image formation problem is in the narrow-field regime. No correction is required for the non-coplanar baselines effect in image formation, and all image-plane calibration effects can be treated as visibility-plane calibration effects. All calibration terms are factorized accordingly as antenna-based. We adopt a dual-circular orthogonal polarization basis ((p,q) ∈  {R,L}), in keeping with our instrumental focus on the VLBA.

2.3. Atmosphere

Ionospheric Faraday rotation (IFR) produces a differential polarization phase offset ${\mathit{\gamma }}^{\mathit{R}\mathrm{-}\mathit{L}}\mathrm{\propto }\frac{\mathrm{1}}{{\mathit{\nu }}^{\mathrm{2}}}$, where ν is the observing frequency. Using a maximum daytime IFR estimated in Thompson et al. (2001) from Evans & Hagfors (1968) to be 15 rotations at ν = 100 MHz for a line-of-sight at a zenith angle z = 60°, we predict approximate IFR values of γ^R−L ~ 0.03° and γ^R−L ~ 0.007° at 43 GHz and 86 GHz respectively, consistent with intuitive expectations. This effect can therefore be ignored in the current analysis: γ^R−L ≪ 1°.

The troposphere is non-dispersive and produces no differential polarization calibration offsets in any of the calibration quantities. However, the time-variable water vapor and dry atmosphere constituents are major contributors to absorption in the millimeter-wavelength bands considered here. We regard the atmosphere as isothermal, with temperature T_atm. For an atmospheric line of sight of optical depth τ, the noise power contribution from the isothermal atmosphere is (1 − e^−τ)  T_atm.

2.4. Image-plane effects

As noted above, the narrow-field constraint applicable in this case (and in most standard VLBI observations) allows relevant image-plane effects to be factorized as antenna-based corrections in the visibility plane. This does not detract from the origin of these effects in the image plane however. For the observing wavelengths and specific instrumental emphasis considered here, the most important image-plane effects originate from antenna pointing errors. These pointing errors change the region of the antenna primary beam located toward the angular direction of the correlated field center. This causes variation in the calibration terms that vary with angular separation from the optical center of the antenna primary beam, such as the amplitude and polarization response.

The nominal pointing performance and beam squint of the VLBA antennas at these observing wavelengths are summarized in Table 1. In the presence of pointing errors, the VLBA primary beam squint (△ _s ~ 5% of θ) introduces a differential polarization amplitude error. This problem can be analyzed in one dimension without loss of generality. We model the one-dimensional normalized primary beam power pattern at each antenna along angular axis ζ as a Gaussian function of width θ, beam squint  △ _s, and pointing error ϵ_p, for nominal orthogonal receptor polarizations (p,q) ∈  {R,L} as: $\begin{array}{ccc}{\mathit{A}}^{\mathit{p}}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{\theta }}^{\mathrm{2}}}\mathrm{(}\mathit{\zeta }\mathrm{-}\frac{{\mathrm{△}}_{\mathrm{s}}}{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{p}}{\mathrm{)}}^{\mathrm{2}}}\\ {\mathit{A}}^{\mathit{q}}\mathrm{\left(}\mathit{\zeta }\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{\theta }}^{\mathrm{2}}}\mathrm{(}\mathit{\zeta }\mathrm{+}\frac{{\mathrm{△}}_{\mathrm{s}}}{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{p}}{\mathrm{)}}^{\mathrm{2}}}\mathit{.}\end{array}$The differential polarization amplitude error in this model takes the form: $\begin{array}{ccc}\frac{{\mathit{A}}^{\mathit{p}}}{{\mathit{A}}^{\mathit{q}}}& \mathrm{=}& {\mathrm{e}}^{\frac{\mathrm{8}\ln \mathrm{2}}{{\mathit{\theta }}^{\mathrm{2}}}\mathrm{(}\mathit{\zeta }\mathrm{+}{\mathit{ϵ}}_{\mathrm{p}}\mathrm{)}{\mathrm{△}}_{\mathrm{s}}}\\ & \mathrm{\simeq }& \mathrm{1}\mathrm{+}\frac{\mathrm{8}\ln \mathrm{2}}{{\mathit{\theta }}^{\mathrm{2}}}\mathrm{(}\mathit{\zeta }\mathrm{+}{\mathit{ϵ}}_{\mathrm{p}}\mathrm{)}{\mathrm{△}}_{\mathrm{s}}\mathit{.}\end{array}$(3)The R/L amplitude ratio at the nominal pointing position ζ = 0 is plotted as a function of pointing error ϵ_p in Fig. 1, for beam squint ${\mathrm{△}}_{\mathrm{s}}\mathrm{=}\frac{{\mathit{\theta }}_{\mathit{b}}}{\mathrm{20}}$, as applicable to the VLBA.

Fig. 1 Plot of normalized primary beam patterns in 1-d cross-section for receptor polarizations  {R,L} with a beam squint equal to 5% of the full-width at half-maximum θb; the x-ordinate for the beam pattern plot is in units of θb. The ratio of the two beam patterns is also plotted in the upper part of the diagram as a function of pointing error ϵp, with x-ordinate $\frac{{\mathit{ϵ}}_{\mathrm{p}}}{{\mathit{\theta }}_{\mathit{b}}}$).

This effect is clearly significant for the representative pointing rms values enumerated in Table 1 for the VLBA operating at these observing wavelengths. We discuss calibration methods for this effect in the second half of this paper.

The net instrumental polarization response at each antenna in each nominal circular receptor polarization (p,q) ∈  {R,L} can be parametrized in terms of the net receptor ellipticity and orientation or equivalently as a complex D-term representing the net leakage of unwanted signal from the orthogonal nominal polarization (Conway & Kronberg 1969). The net instrumental polarization response contains contributions from the antenna optics design, feed illumination function, polarizer performance, and to a lesser extent, other electronic components in the downstream signal path that have a differential polarization response. The net instrumental polarization is an image-plane (direction-dependent) effect but, in the current narrow-field scientific application, we again factorize it as an antenna-based visibility plane correction factor, in common with the other image-plane calibration effects considered in this section. Pointing errors introduce time-variability in the net instrumental polarization response, expressed as an antenna-based visibility plane correction, by sampling different regions of the primary beam polarization pattern, a contributor to the net D-term. The associated dependence on frequency is neglected here due to the small fractional bandwidth.

2.5. Total power

Total power relationships for single-dish telescopes are discussed by Kutner & Ulich (1981), Downes (1989), and Wilson et al. (2009). Adapting this references, we represent the total system temperature at an individual antenna in nominal receptor polarization p in the form: $\begin{array}{ccc}{\mathit{T}}_{\mathrm{sys}}^{\mathit{p}}& \mathrm{=}& {\mathit{\eta }}_{\mathit{b}}^{\mathit{p}}{\mathit{\eta }}_{\mathrm{l}}{\mathit{T}}_{\mathit{B}}^{\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\eta }}_{\mathrm{l}}\mathrm{)}{\mathit{T}}_{\mathrm{spill}}\mathrm{+}{\mathit{T}}_{\mathrm{rx}}^{\mathit{p}}\mathrm{+}{\mathit{\eta }}_{\mathrm{l}}{\mathit{T}}_{\mathrm{atm}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}\\ {\mathit{T}}_{\mathrm{sys}}^{\mathit{p}}& \mathrm{=}& {\mathit{\eta }}_{\mathit{b}}^{\mathit{p}}{\mathit{\eta }}_{\mathrm{l}}{\mathit{T}}_{\mathit{B}}^{\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{+}{\mathit{T}}_{\mathit{N}}^{\mathit{p}}\end{array}$where ${\mathit{\eta }}_{\mathit{b}}^{\mathit{p}}\mathrm{=}{\mathit{A}}^{\mathit{p}}\mathrm{\left(}{\mathit{ϵ}}^{\mathit{p}}\mathrm{\right)}$ is the amplitude loss for an angular pointing error ϵ^p, the term ${\mathit{T}}_{\mathit{B}}^{\mathit{p}}$ is the source Rayleigh-Jeans brightness temperature (radiation temperature) in receptor polarization p, η_l includes antenna scattering and antenna resistive (ohmic) losses, τ = τ(z) is the optical depth along the line of sight at zenith angle z through the atmosphere, T_spill is the mean ground spill-over temperature, ${\mathit{T}}_{\mathrm{rx}}^{\mathit{p}}$ is the net receiver temperature in receptor polarization p, and T_atm is the mean temperature of the isothermal atmosphere along the line of sight.

2.6. Autocorrelation and cross-power spectra

Autocorrelation and cross-power spectra are formed by interferometric correlation of all polarization pairs for the voltages recorded in each nominal receptor polarization at each antenna in the array. We assume idealized analog correlation here for clarity but without loss of generality – the effects of quantization, sampling, and other digital signal processing effects are considered in subsequent sections. Following the notation of Kemball et al. (1995), we denote the normalized autocorrelation and cross-power spectra as ${\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}$ and ${\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}$ respectively for antenna indices (m,n) and nominal receptor polarizations (p,q). The spectra are normalized by the geometric mean of the total power in the signal received at each antenna (m,n) forming part of the correlation pair.

In this section we express the data model for ${\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}$ and ${\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}$ in the more modern measurement equation (ME) matrix formalism introduced by Hamaker et al. (1996) and Sault et al. (1996).

2.6.1. Autocorrelation data model

We denote the normalized analog autocorrelation spectrum as a 4-vector $R\mathrm{=}{\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}$; this is real-valued for polarization p = q and complex for p ≠ q. $R\mathrm{=}\left(\begin{array}{c}\\ {\mathit{r}}_{\mathit{mm}}^{\mathit{RR}}\\ {\mathit{r}}_{\mathit{mm}}^{\mathit{RL}}\\ {\mathit{r}}_{\mathit{mm}}^{\mathit{LR}}\\ {\mathit{r}}_{\mathit{mm}}^{\mathit{LL}}\end{array}\right)\mathrm{(}\mathit{u}\mathrm{=}\mathrm{0}\mathit{,v}\mathrm{=}\mathrm{0}\mathit{,\omega }\mathrm{)}$(6)where (u,v) are visibility-plane coordinates and ω denotes angular frequency. We denote the true source correlation spectrum as J. In the autocorrelation case, this is the (antenna-independent) source total-power spectrum, $J\mathrm{=}\left(\begin{array}{c}\\ {\mathit{J}}^{\mathit{RR}}\\ {\mathit{J}}^{\mathit{RL}}\\ {\mathit{J}}^{\mathit{LR}}\\ {\mathit{J}}^{\mathit{LL}}\end{array}\right)\mathrm{(}\mathit{u}\mathrm{=}\mathrm{0}\mathit{,v}\mathrm{=}\mathrm{0}\mathit{,\omega }\mathrm{)}\mathit{.}$(7)We denote a total-power noise term contribution , where ${\mathrm{𝒩}}^{\mathit{p}}\mathrm{=}{\mathit{P}}_{\mathit{m}}^{\mathit{p}}{\mathit{T}}_{\mathit{N,m}}^{\mathit{p}}$. Following Kemball et al. (1995), we use ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}$ to denote the point-source sensitivity (in Jy/K) of antenna m in receptor polarization p. The total system temperature at antenna m in polarization p is denoted as ${\mathit{T}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}$ and follows Eq. (5).

The autocorrelation data model can then be represented in matrix form as, $\begin{array}{ccc}& & R\mathrm{=}{K}_{\mathit{mm}}{L}_{\mathit{mm}}\mathrm{(}J\mathrm{+}{N}_{\mathit{m}}\mathrm{)}\\ & & K_{\mathit{mm}}\mathrm{=}\mathrm{(}{G}_{\mathit{m}}\mathrm{\otimes }{G}_{\mathit{m}}^{\mathrm{\ast }}\mathrm{)}\mathrm{(}{B}_{\mathit{m}}\mathrm{\otimes }{B}_{\mathit{m}}^{\mathrm{\ast }}\mathrm{)}\\ & & L_{\mathit{mm}}\mathrm{=}\mathrm{(}{D}_{\mathit{m}}\mathrm{\otimes }{D}_{\mathit{m}}^{\mathrm{\ast }}\mathrm{)}\mathrm{(}{P}_{\mathit{m}}\mathrm{\otimes }{P}_{\mathit{m}}^{\mathrm{\ast }}\mathrm{)}\end{array}$(8)where  ⊗  denotes an outer matrix product, and G_m,B_m,D_m and P_m are 2 × 2 Jones matrices denoting gain normalization, complex bandpass response, instrumental polarization, and parallactic angle feed rotation respectively. This follows the measurement equation formalism of Hamaker et al. (1996) but includes a self-noise term, as required in the autocorrelation data model considered here. In this circularly-polarized basis, the Jones matrices take the form, $\begin{array}{ccc}{G}_{\mathit{m}}\mathrm{=}\left(\begin{array}{c}\\ & \\ & \end{array}\right)& & \\ B_{\mathit{m}}\mathrm{=}\left(\begin{array}{c}\\ & \\ & \end{array}\right)& & \\ D_{\mathit{m}}\mathrm{=}\left(\begin{array}{c}\\ & \\ & \end{array}\right)& & \\ P_{\mathit{m}}\mathrm{=}\left(\begin{array}{c}\\ & \\ & \end{array}\right)& & \end{array}$K_mm is a therefore a diagonal 4 × 4 matrix containing the gain normalization ${\mathit{G}}_{\mathit{m}}^{\mathit{p}}$ and complex bandpass response terms ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, $K_{\mathit{mm}}\mathrm{=}\begin{array}{}\mathrm{diag}(\mathrm{\parallel }{\mathit{G}}_{\mathit{m}}^{\mathit{R}}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{R}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathit{,}\\ {\mathit{G}}_{\mathit{m}}^{\mathit{R}}{\mathit{G}}_{\mathit{m}}^{\mathit{L}\mathrm{\ast }}{\mathit{B}}_{\mathit{m}}^{\mathit{R}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{B}}_{\mathit{m}}^{\mathit{L}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{,}\\ {\mathit{G}}_{\mathit{m}}^{\mathit{L}}{\mathit{G}}_{\mathit{m}}^{\mathit{R}\mathrm{\ast }}{\mathit{B}}_{\mathit{m}}^{\mathit{L}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{B}}_{\mathit{m}}^{\mathit{R}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{,}\\ \mathrm{\parallel }{\mathit{G}}_{\mathit{m}}^{\mathit{L}}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{L}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}})\mathit{.}\end{array}$(13)The net instrumental polarization response L_mm including D-terms and parallactic angle rotation is the 4 × 4 matrix, $L_{\mathit{mm}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& {\mathit{D}}_{\mathit{m}}^{\mathit{R}\mathrm{\ast }}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{j\alpha }}& {\mathit{D}}_{\mathit{m}}^{\mathit{R}}{\mathrm{e}}^{\mathrm{2}\mathit{j\alpha }}& \mathrm{\parallel }{\mathit{D}}_{\mathit{m}}^{\mathit{R}}{\mathrm{\parallel }}^{\mathrm{2}}\\ {\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathrm{\ast }}& {\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{j\alpha }}& {\mathit{D}}_{\mathit{m}}^{\mathit{R}}{\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathrm{\ast }}{\mathrm{e}}^{\mathrm{2}\mathit{j\alpha }}& {\mathit{D}}_{\mathit{m}}^{\mathit{R}}\\ {\mathit{D}}_{\mathit{m}}^{\mathit{L}}& {\mathit{D}}_{\mathit{m}}^{\mathit{L}}{\mathit{D}}_{\mathit{m}}^{\mathit{R}\mathrm{\ast }}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{j\alpha }}& {\mathrm{e}}^{\mathrm{2}\mathit{j\alpha }}& {\mathit{D}}_{\mathit{m}}^{\mathit{R}\mathrm{\ast }}\\ \mathrm{\parallel }{\mathit{D}}_{\mathit{m}}^{\mathit{L}}{\mathrm{\parallel }}^{\mathrm{2}}& {\mathit{D}}_{\mathit{m}}^{\mathit{L}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{j\alpha }}& {\mathit{D}}_{\mathit{m}}^{\mathit{L}\mathrm{\ast }}{\mathrm{e}}^{\mathrm{2}\mathit{j\alpha }}& \mathrm{1}\end{array}\right)\mathit{.}$(14)

2.6.2. Cross-correlation data model

The cross-correlation data model does not include a self-noise term; in this case Eq. (8) takes the form, $\begin{array}{ccc}& & \\ & & K_{\mathit{mn}}\mathrm{=}\mathrm{(}{G}_{\mathit{m}}\mathrm{\otimes }{G}_{\mathit{n}}^{\mathrm{\ast }}\mathrm{)}\mathrm{(}{B}_{\mathit{m}}\mathrm{\otimes }{B}_{\mathit{n}}^{\mathrm{\ast }}\mathrm{)}\\ & & L_{\mathit{mn}}\mathrm{=}\mathrm{(}{D}_{\mathit{m}}\mathrm{\otimes }{D}_{\mathit{n}}^{\mathrm{\ast }}\mathrm{)}\mathrm{(}{P}_{\mathit{m}}\mathrm{\otimes }{P}_{\mathit{n}}^{\mathrm{\ast }}\mathrm{)}\end{array}$(15)where X is a diagonal 4 × 4 matrix term for cross-correlation coherence losses; this term does not apply to autocorrelation data. X does not contribute to differential polarization amplitude offsets, and forms part of the calibration of the absolute flux density scale.

Here, J_mn is the true source visibility cross-power spectrum, $J_{\mathit{mn}}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{J}}_{\mathit{mn}}^{\mathit{RR}}\\ {\mathit{J}}_{\mathit{mn}}^{\mathit{RL}}\\ {\mathit{J}}_{\mathit{mn}}^{\mathit{LR}}\\ {\mathit{J}}_{\mathit{mn}}^{\mathit{LL}}\end{array}\right)\mathrm{\left(}\mathit{u,v,\omega }\mathrm{\right)}$(16)as a function of (u,v)-spacing and frequency ω.

2.7. Sampling and quantization

The voltages in each nominal receptor polarization are digitally sampled and quantized by the data acquisition sub-systems at each antenna. Digital realizations ${\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}$ of the total-power and cross-power spectra are obtained by subsequent cross-correlation. We denote the relationship between the measured spectra ${\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}$ and their true analog counterparts ${\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}$ used in earlier sections as a general transformation function f_d: $\begin{array}{ccc}& & {\mathit{f}}_{\mathit{d}}\mathrm{:}{\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}\mathrm{\to }{\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}\\ & & {\mathit{f}}_{\mathit{d}}^{-1}\mathrm{:}{\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}\mathrm{\to }{\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}\mathit{.}\end{array}$The form of this relationship depends on quantization level, sampler voltage threshold stability and accuracy, and correlator architecture, and is reviewed by Thompson et al. (2001) and references therein. This transformation function does not always take closed analytic form.

The specific case of sampling and quantization effects for the VLBA correlator is considered by Kogan (1998). This correlator has an FX architecture, originally hardware-based but more recently upgraded to a VLBA-DiFX software implementation (Deller et al. 2007). Both one- and two-bit sampling are supported. An offset in sampler threshold voltage can produce an associated amplitude error in the output VLBA correlator spectra (Kogan 1998), particularly for two-bit sampling. If the mapping of individual sampler modules to sets of baseband converters mirrors receptor polarization baseband assignment – as occurs naturally for several common baseband converter configurations – these amplitude offsets will translate to differential polarization amplitude errors. The correction of quantization and sampler threshold errors is discussed in further detail in the following section.

3. Data reduction methods

In this section we assess the data reduction methods used to measure and correct the instrumental and propagation terms in the data model enumerated in Sect. 2. Consistent with earlier practice, we describe here only differences with the data reduction methodology presented by Kemball et al. (1995) and use identical notation wherever possible. This section follows the proposed data reduction sequence, which is shown in flowchart form in Fig. 2.

We have implemented these algorithms within a development framework that brings together several large community codes used for radio-astronomical imaging (Kemball et al. 2008). For consistency with earlier algorithmic work in Kemball et al. (1995) we have primarily used an adapted implementation of the Astronomical Image Processing System² within this software framework for the current work. Our intention is to make these revisions publicly available in the future.

Fig. 2 Flowchart depicting the overall data reduction sequence described in Sect. 3.

3.1. Sampling and quantization

In this section, we review the correction of digital sampling and quantization effects in VLBI circular polarimetry at millimeter wavelengths. We do not find it necessary to revise existing algorithmic practice, as described in the analysis below.

The total-power and cross-power spectra ${\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}$ from the VLBA FX correlator are first corrected for quantization effects and then scaled in amplitude by polarization-independent factors known a-priori for the array (the latter known historically as b-factors). The quantization correction element of the transformation function f_d defined in Sect. 2.7 above is defined as a relation for each quantization mode (one- or two-bit, η_Q = 11 or η_Q = 22 respectively) between a measured correlation ρ_m and the true underlying correlation ρ, both defined in the conjugate Fourier delay lag domain. In lag space, the correlation is of order $\mathit{O}{}^{\mathrm{\right(}}\frac{{\mathit{T}}_{\mathit{B}}}{{\mathit{T}}_{\mathit{B}}\mathrm{+}{\mathit{T}}_{\mathit{N}}}^{\mathrm{\left)}}$ for cross-power spectra; for the VLBA system temperatures listed in Table 1 and for typical SiO maser spectral flux densities, lag-domain correlation values are small |ρ_m| ≲ 0.1. In this case, the relation between ρ_m and ρ is in the highly linear regime (Kogan 1998) and the quantization correction can be applied as a direct linear scaling of cross-power spectral values ${\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}\mathrm{=}{\mathit{\alpha }}_{\mathit{Q}}{\mathit{r}}_{\mathit{d,mn}}^{\mathit{pq}}$ with numerical coefficient α_Q (Kogan 1993a). As this takes the form of a linear, polarization-independent correction, it cannot produce an artificial instrumental signature mimicking Stokes V. Cross-polarized autocorrelation spectra ${\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}\mathrm{(}\mathit{p}\ne \mathit{q}\mathrm{)}$ have lag domain correlation values of order $\mathit{O}\left(\frac{{\mathit{T}}_{\mathit{B}}\mathrm{+}\sqrt{\mathrm{2}}\mathrm{\parallel }\mathit{D}\mathrm{\parallel }{\mathit{T}}_{\mathit{N}}}{{\mathit{T}}_{\mathit{B}}\mathrm{+}{\mathit{T}}_{\mathit{N}}}\right)$, as indicated by the first-order form of Eq. (8), and similarly fall in the low-correlation regime. Accordingly they have the same quantization scaling correction as cross-power spectra.

Parallel-hand autocorrelation spectra ${\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}\mathrm{(}\mathit{p}\mathrm{=}\mathit{q}\mathrm{)}$ have lag-domain correlation values of unity at zero delay, and require Fourier transformation to and from the lag domain to allow application of the full non-linear correction between ρ_m and ρ (Kogan 1993a). As a result, this does not reduce to a linear scaling relationship for the total-power spectra ${\mathit{r}}_{\mathit{mm}}^{\mathit{pq}}\mathrm{(}\mathit{p}\mathrm{=}\mathit{q}\mathrm{)}$ in the frequency domain. Two effects limit the accuracy of the inverse transform to the lag domain. The VLBA FX correlator records only positive spectral channels (Kogan 1993a), omitting the opposite sideband which contains digitization noise (Thompson et al. 2001). In addition, an FX correlator requires equal-length zero padding to reconstruct the associated correlation function over the full range of sampled delay lags (Granlund 1986; Thompson et al. 2001).

Sampler threshold voltages define transition points between discrete quantization states. The non-linear mapping between ρ_m and ρ depends functionally on the value of the sampler threshold voltages relative to their optimal values (Kogan 1995). For two-bit sampling in the low-correlation limit, the relative error in FX spectral output depends linearly on the deviation in threshold voltage level. In one-bit sampling, the relative error in spectral output depends only in second order on the error in the threshold clipping voltage (Kogan 1995). In both VLBA quantization modes, the correction in output spectral amplitude can be derived from the measured mean deviation of the total-power spectra from the unit mean power level predicted for an ideal digitizer, independent of spectral shape (Kogan 1995), and we adopt that approach here.

The autocorrelation data are used in integrated template-fitting as part of calibration, as described in further detail below.

3.2. Bandpass calibration

Following Kemball et al. (1995), and the Jones matrix formalism in Eq. (10) in the current work, we denote the cross-power complex bandpass response as ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$; this is expanded as ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\parallel }{\mathrm{e}}^{\mathit{j}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}$. The autocorrelation bandpass amplitude response, denoted $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$, differs from the cross-power bandpass amplitude response $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ however, due to an unavoidable level of irreducible aliasing in the net system bandpass response. As the continuum calibrators have flat spectra across the baseband bandwidth they are used to solve for ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ and $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ using cross- and total-power continuum calibrator data respectively.

3.2.1. Bandpass frequency frame

The VLBA correlator produces output cross- and total-power spectra in a geocentric J2000.0 coordinate reference frame. In contrast, the bandpass response function ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ is defined in the data acquisition coordinate reference frame at each antenna, which is an apparent topocentric frame. These two frequency frames are offset by the time-variable natural geometric fringe rate at antenna m, denoted  △ ω_m(t). The angular frequency of channel number l in the recorded cross- and total-power visibility spectra (which are in a geocentric frame) is denoted as ω_l. The autocorrelation bandpass amplitude response ${\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ at antenna m in receptor polarization p ∈  {R,L} can be solved for by minimizing: $\begin{array}{ccc}{\mathrm{[}{{\mathit{\chi }}_{\mathit{m}}^{\mathit{pp}}}^{\mathrm{]}}}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{t}}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{c}}}{\mathrm{[}\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathrm{l}}\mathrm{-}\mathrm{△}{\mathit{\omega }}_{\mathit{m}}\mathrm{(}{\mathit{t}}_{\mathit{k}}\mathrm{\left)}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}{\begin{array}{}\mathrm{˜}\\ {\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\end{array}}^{\mathrm{]}}}^{\mathrm{2}}\\ \begin{array}{}\mathrm{˜}\\ {\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\end{array}& \mathrm{=}& \mathrm{\parallel }{\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }\mathit{/}\overline{\mathrm{\parallel }{\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }}\end{array}$where ${\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}$ is the (real-valued) continuum calibrator total-power spectrum in parallel-hand polarization correlation pp at antenna m over a pre-average integration interval centered at time t_k, the number of frequency channels is N_c, and the number of integration intervals is N_t. The pre-average interval is fixed and short relative to $\frac{\mathrm{d}\mathrm{△}\mathit{\omega }}{\mathrm{d}\mathit{t}}$. Pre-averaging improves the statistical robustness of the bandpass response solver, as this is a least-squares minimization problem. The autoscaling of each total-power spectrum by the instantaneous mean spectral amplitude, $\overline{\mathrm{\parallel }{\mathit{V}}_{\mathit{mm}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }}$ , ensures that the spectra do not need to be calibrated in amplitude before solving for ${\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$. The resulting solution for $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ is normalized to unit mean power $\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{c}}}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{k}\mathrm{=}{\mathit{N}}_{\mathrm{c}}}\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{=}\mathrm{1}$.

Bandpass calibration forms part of visibility amplitude calibration and accurate estimation of Stokes V requires a correspondingly accurate solution for $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }$. For example, a spurious bandpass amplitude spike, even if in an off-source spectral region, will bias the mean amplitude normalization of the bandpass over the region of source emission.

The measured visibility spectra and the unknown bandpass response are not in a common frequency frame; the channel shift corresponding to  △ ω_m is a non-negligible effect at the wavelengths of λ = 3 mm and λ = 7 mm considered in this paper. If $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ is parametrized discretely as a set of unknown values at each frequency channel (in a stationary topocentric frame) then this can be estimated by direct integration of the total-power, parallel-hand spectra if they are first shifted to the stationary reference frame by using a discrete Fourier transform pair including a channel shift. Alternatively, the unshifted visibility spectra can be fit directly if $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ is parametrized as a polynomial series expansion over frequency; this parametrization allows continuous variation of  △ ω_m(t) in the formulation of the chi-square minimization problem. This method was developed by Kemball & Diamond (1997), using a Chebyshev polynomial expansion for the bandpass response. The polynomial bandpass method avoids digital signal processing artifacts associated with the Fourier transform shift. In addition, as the polynomial expansion order can be significantly lower than the number of frequency channels N_c, this approach reduces the number of free parameters over the discrete parametrization case, thus improving the signal-to-noise ratio (SNR) of the bandpass solution. For these reasons (amongst others discussed below), we adopt the bandpass polynomial solution method. In a Chebyshev expansion, the autocorrelation bandpass response takes the form: $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }\mathrm{=}\frac{{\mathit{c}}_{\mathit{a,}\mathrm{0}}}{\mathrm{2}}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{p}}}{\mathit{c}}_{\mathit{a,j}}{\mathit{T}}_{\mathit{j}}\mathrm{\left(}\mathit{x}\mathrm{\right(}\mathit{\omega }\mathrm{\left)}\mathrm{\right)}$(21)where T_j(x) is the Chebyshev polynomial of degree j (Press et al. 2007), c_a,j are real coefficients in the series expansion for the autocorrelation bandpass amplitude, and the transformed coordinate is x(ω) = (2ω − ω_a − ω_b)/(ω_b − ω_a). The topocentric channel range for the bandpass solution covered by the measured visibility spectra is  [ω_a = 1 − max(△ ω_m),  ω_b = N_c − min(△ ω_m)] . In the current work, we found it important to fit over the full topocentric range  [ω_a,ω_b]  as opposed to only  [1,N_c]  as used in past practice, so that there is the greatest possible chi-square constraint on the bandpass solution at the edge of the frequency channel range. The derived solution coefficients c_a,j resulting from the chi-square minimization over  [ω_a,ω_b]  are then transformed to ${\mathit{c}}_{\mathit{a,j}}^{\mathrm{\prime }}$ over  [1,N_c]  as (Press et al. 2007), $\begin{array}{ccc}& & {\mathit{\omega }}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{\omega }}_{\mathit{b}}\mathrm{-}{\mathit{\omega }}_{\mathit{a}}}{\mathrm{2}}\cos \left(\frac{\mathit{\pi }\mathrm{(}\mathit{k}\mathrm{-}\mathrm{0.5}\mathrm{)}}{{\mathit{N}}_{\mathrm{c}}}\right)\mathrm{+}\frac{{\mathit{\omega }}_{\mathit{a}}\mathrm{+}{\mathit{\omega }}_{\mathit{b}}}{\mathrm{2}}\\ & & {\mathit{f}}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\frac{{\mathit{c}}_{\mathit{a,}\mathrm{0}}}{\mathrm{2}}\mathrm{+}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{p}}}{\mathit{c}}_{\mathit{a,j}}{\mathit{T}}_{\mathit{j}}\mathrm{\left(}\mathit{x}\mathrm{\right(}{\mathit{\omega }}_{\mathit{k}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right)}\\ & & {\mathit{c}}_{\mathit{a,j}}^{\mathrm{\prime }}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{p}}}{\mathit{f}}_{\mathit{k}}^{\mathrm{\prime }}\cos \left(\frac{\mathit{\pi }\mathrm{(}\mathit{j}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{0.5}\mathrm{)}}{{\mathit{N}}_{\mathrm{c}}}\right)\mathrm{·}\end{array}$During this coefficient transformation, if ω_b < N_c or ω_a > 1 then the outer ${\mathit{f}}_{\mathit{k}}^{\mathrm{\prime }}$ are extrapolated horizontally to cover the complete range  [1,N_c] .

We note that Eq. (19) is formulated for the case of a time-invariant bandpass at each antenna. This is an appropriate instrumental assumption given the expected stability of the net VLBA frequency response, and improves the signal-to-noise ratio of the resulting bandpass solution. In addition a single bandpass solution provides greater control over the calibration of overall differential polarization R/L amplitude and phase offsets.

In the most direct formulation, the cross-power bandpass response ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ can be solved for from the parallel-hand continuum calibrator cross-power spectra ${\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}$ by an analogous formulation of Eq. (19), $\begin{array}{ccc}{\left[{\mathit{\chi }}^{\mathit{pp}}\right]}^{\mathrm{2}}& \mathrm{=}& \\ \begin{array}{}\mathrm{˜}\\ {\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\end{array}& \mathrm{=}& {\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\mathit{/}\overline{{\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}}\\ {\mathit{\omega }}_{\mathit{lm}}& \mathrm{=}& {\mathit{\omega }}_{\mathrm{l}}\mathrm{-}\mathrm{△}{\mathit{\omega }}_{\mathit{m}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\\ {\mathit{\omega }}_{\mathit{ln}}& \mathrm{=}& {\mathit{\omega }}_{\mathrm{l}}\mathrm{-}\mathrm{△}{\mathit{\omega }}_{\mathit{n}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}\end{array}$(25)where $\overline{{\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}}$ is the complex normalization factor to auto-scale the cross-power spectrum in integration t_k over frequency to unit mean amplitude and zero mean phase. The auto-scaling applied during the chi-square minimization removes the requirement for amplitude calibration of ${\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}$, as in the case for the autocorrelation bandpass ${\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, but for the cross-power bandpass solution, the visibility data need to be calibrated in the phase domain for residual group delay, fringe-rate, phase, and differential polarization phase offset, as described by Kemball et al. (1995). The resulting solutions for ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ are normalized over frequency to unit mean power $\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{c}}}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{c}}}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ and zero mean phase. The polynomial bandpass method is further favored in the cross-power case as there is no unique shift defined for V_mn(t_k) if instead of a polynomial bandpass expansion, the visibility data are shifted using a Fourier transform.

3.2.2. Bandpass aliasing correction

The solutions obtained for the cross-power bandpass response ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ from the direct formulation in Eq. (25) are not optimal however for precise amplitude calibration at the millimeter wavelengths (λ = 3 mm and λ = 7 mm) considered here. At these wavelengths, the amplitude of the continuum calibrator cross-power spectra on the longer baselines invariably have insufficient SNR to allow estimation of $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ with low mean-square error (MSE) at the outlying antennas. The autocorrelation spectra have significantly higher SNR however, and offer a preferred path for the solution of the cross-power bandpass amplitude response $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ if corrected for aliasing, as found necessary in the current work. We adopt an aliasing model: $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{=}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{+}{\mathit{A}}_{\mathit{b}}^{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{(}\mathrm{2}{\mathit{N}}_{\mathrm{c}}\mathrm{-}\mathit{k}\mathrm{)}}\mathrm{\right)}$(26)where A_b(ω) is a model for the cross-power bandpass amplitude response at and above the upper end of the band (in a frequency range ω_k0 > βω_c, where ω_c is the bandpass cutoff frequency and β lies in the approximate range β ~  {0.7−0.9}). This function can be extrapolated and folded to approximate the aliased autocorrelation response. We find a 12-node Butterworth function (Bianchi & Sorrentino 2007) to give a satisfactory empirical fit to the upper-bandpass region ω_k0 > βω_c, with little sensitivity to the value of β adopted to define the channel range of the fit: ${\mathit{A}}_{\mathit{b}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{=}\mathit{a}\sqrt{\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{(}{\frac{{\mathit{\omega }}_{\mathit{k}}}{{\mathit{\omega }}_{\mathrm{c}}}}^{\mathrm{)}}}^{\mathrm{2}\mathit{n}}}}\mathrm{·}$(27)The fit minimizes: $\begin{array}{ccc}\mathrm{[}{\mathit{\chi }}^{\mathit{pp}}{\mathrm{]}}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{k}\mathrm{=}\mathit{k}\mathrm{0}}^{{\mathit{N}}_{\mathrm{c}}}{\left[\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}{\mathrm{(}{\mathit{A}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}{\mathrm{)}}^{\mathrm{)}}}^{\mathrm{2}}\right]}^{\mathrm{2}}\\ & & \mathrm{+}\sum _{\mathit{k}\mathrm{=}\mathit{k}\mathrm{0}}^{{\mathit{N}}_{\mathrm{c}}\mathrm{-}\mathrm{1}}\mathrm{[}\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathit{k}}\mathrm{)}{\mathrm{\parallel }}^{\mathrm{2}}\\ & & \mathrm{-}{\mathrm{(}{\mathit{A}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathrm{(}\mathrm{2}{\mathit{N}}_{\mathrm{c}}\mathrm{-}\mathit{k}\mathrm{)}}{\mathrm{)}}^{\mathrm{)}}}^{\mathrm{2}}{\mathrm{]}}^{\mathrm{2}}\end{array}$(28)where the number of nodes is fixed at n = 12, and factor a and the cutoff frequency ω_c are the only fitted parameters. The direct cross-power solutions for $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ obtained from Eq. (25) do however have sufficient SNR to allow a solution for A_b(ω) in Butterworth form (27). The solution for A_b(ω) allows the high-SNR autocorrelation bandpass amplitude solutions $\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }$ to be then transformed to cross-power form: $\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}{\mathit{A}}_{\mathit{b}}^{\mathrm{2}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{(}\mathrm{2}{\mathit{N}}_{\mathrm{c}}\mathrm{-}\mathit{k}\mathrm{)}}\mathrm{\right)}}\mathit{.}$(29)In this hybrid approach, the phase ${\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ of the cross-power bandpass response is solved for by minimizing the phase-only analog of Eq. (25): $\begin{array}{ccc}{\left[{\mathit{\chi }}^{\mathit{pp}}\right]}^{\mathrm{2}}& \mathrm{=}& \sum _{\mathit{m,n}\mathit{>}\mathit{m}}^{{\mathit{N}}_{\mathit{a}}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{t}}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{c}}}[{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{(}{\mathit{\omega }}_{\mathrm{l}}\mathrm{-}\mathrm{△}{\mathit{\omega }}_{\mathit{m}}\mathrm{(}{\mathit{t}}_{\mathit{k}}\mathrm{\left)}\mathrm{\right)}\\ & & \end{array}$(30)where $\begin{array}{}{}_{\mathrm{˜}}\\ {\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\end{array}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{l}}\mathit{,}{\mathit{t}}_{\mathit{k}}\mathrm{\right)}$ are the pre-averaged continuum calibrator cross-power visibility spectra, normalized to unit mean amplitude and zero mean phase, and calibrated in delay, fringe-rate, and phase as described above. In practice, the peak-to-peak residual phase error over frequency for the net VLBA baseband response after correction for phase and group delay is typically of order 5 degrees, and we fit ${\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ at each antenna in each receptor polarization with a low-order polynomial in order to maximize SNR. The complex cross-power bandpass response ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ is therefore represented by a separate Chebyshev polynomial in each of amplitude and phase, with a lower-order polynomial in phase compared to that in amplitude. The final cross-power bandpass solution is therefore constructed as, ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{\parallel }{\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{\prime }}{\mathrm{e}}^{\mathit{j}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}$(31)applying the unit-mean power and zero-mean phase normalization to the bandpass solutions as described above.

3.2.3. Reference antenna differential polarization bandpass phase response

We note that the solution for ${\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ is derived by solving a self-calibration problem that is linear in baseline-based phase, and so is known only relative to the phase at receptor polarization p for an adopted reference antenna (subscript zero), i.e. the determined bandpass phase solution takes the form ${\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{-}{\mathit{\varsigma }}_{\mathrm{0}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, as opposed to ${\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$, the argument of ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{\parallel }{\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\parallel }{\mathrm{e}}^{\mathit{j}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}$ in the Jones matrix of Eq. (10). A full solution, allowing correction of cross-polarized visibility spectra therefore requires an independent estimate of the differential polarization bandpass phase response at the reference antenna ${\mathit{\varsigma }}_{\mathrm{0}}^{\mathit{p}\mathrm{-}\mathit{q}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$. This is directly analogous to the correction of parallel-hand phase solutions for differential polarization phase and delay offsets at the reference antenna, as described by Kemball et al. (1995). The correction for ${\mathit{\varsigma }}_{\mathrm{0}}^{\mathit{p}\mathrm{-}\mathit{q}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ has not traditionally been applied to bandpass corrections in spectral-line VLBI but is relevant to the autocorrelation polarimetry described in subsequent sections of this paper, so was implemented in the current work. Without this correction, for example, there would be no bandpass phase correction applied to the cross-polarized autocorrelation spectra at the reference antenna.

3.2.4. Bandpass correction

The basic algebra for applying the autocorrelation bandpass response ${\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ and ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ to correct the visibility data is described by Kemball et al. (1995) and Diamond (1989). These references describe a traditional $\frac{\mathrm{ON}\mathrm{-}\mathrm{OFF}}{\mathrm{OFF}}\mathrm{=}{}^{\mathrm{(}}\frac{\mathrm{ON}}{\mathrm{OFF}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}$ heuristic for autocorrelation bandpass correction, to minimize the residual total-power offset. However, the complete autocorrelation data model used here (Eq. (8)) does not include this term – nor is it necessary – so we do not subtract 1 when correcting the autocorrelation bandpass response in the current work. Additionally, during bandpass application, the fringe-rate shift is applied when computing the polynomial expansion of ${\mathit{B}}_{\mathit{a,m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ or ${\mathit{B}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$.

3.3. Amplitude calibration

Amplitude calibration of the correlation spectra ${\mathit{r}}_{\mathit{mn}}^{\mathit{pq}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ in units of spectral flux density ${\mathit{J}}_{\mathit{mn}}^{\mathit{pq}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ (Jy) requires independent measurement of the total system temperature ${\mathit{T}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}$ (K) and point-source sensitivity ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}$ (Jy K^-1) both at each antenna m and in each receptor polarization p throughout the course of the observation. We denote the system equivalent flux density (SEFD) as ${\mathit{J}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}\mathrm{=}{\mathit{P}}_{\mathit{m}}^{\mathit{p}}{\mathit{T}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}$ (Jy). In the formalism used in this paper (see Eq. (13) and following), ${\mathit{G}}_{\mathit{m}}^{\mathit{p}}\mathrm{=}\mathrm{1}\mathit{/}\sqrt{{\mathit{J}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}}$. The ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}$ are time-variable due to changes in both atmospheric attenuation and antenna gravitational deformation, both as a function of time and antenna pointing position in local horizon coordinates.

Sections 3.3.1 and 3.3.2 contain a review and analysis of current practice in VLBI continuum and spectral-line amplitude calibration, with a specific focus on circular polarimetry at millimeter wavelengths. We describe the innovations in amplitude calibration introduced in the current work in Sect. 3.3.3.

3.3.1. Continuum amplitude calibration

The VLBA records integrated system temperatures every two minutes at each antenna in each receptor polarization p ∈  {R,L} obtained using an underlying switched noise calibration system (Thompson 1995). The VLBA project also publishes opacity-corrected gain curves for each antenna, which provide ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}\mathit{\tau }\mathrm{=}\mathrm{0}\mathit{,z}\mathrm{)}$ as a polynomial function of zenith angle z; these curves are obtained from analyses of regular single-dish service observations, separately scheduled (Walker 1999).

In standard a priori amplitude calibration of VLBI continuum observations, the ${\mathit{G}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ are computed as $\mathrm{1}\mathit{/}\sqrt{{\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{T}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}}$ from the measured system temperature values and published antenna point-source sensitivity curves. At the observing wavelengths λ = 7 mm and λ = 3 mm considered in this paper the opacity-corrected ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}\mathit{\tau }\mathrm{=}\mathrm{0}\mathit{,z}\mathrm{)}$ provided by the VLBA require correction for atmospheric attenuation e^−τ (Walker 1999; and see Eq. (5)). As described by Leppänen (1993), these corrections can be obtained by fitting the measured system temperatures over the course of the observations as a function of zenith angle z against a simplified form of Eq. (5): ${\mathit{T}}_{\mathrm{sys}}^{\mathit{p}}\mathrm{=}{\mathit{T}}_{\mathit{a}}^{\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{+}{\mathit{T}}_{\mathrm{spill}}\mathrm{+}{\mathit{T}}_{\mathrm{rx}}^{\mathit{p}}\mathrm{+}{\mathit{T}}_{\mathrm{atm}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}$(32)where ${\mathit{T}}_{\mathit{a}}^{\mathit{p}}$ is the antenna temperature contribution from the source. The number of free parameters can be controlled by adopting a plane-parallel atmosphere, with τ = τ₀secz for zenith opacity τ₀, extrapolating T_atm from ground-level metrology, and adopting an empirical model for the spill-over noise contribution for VLBA antennas (Leppänen 1993). In this case, only receiver temperature T_rx and zenith opacity τ₀ remain as free parameters, and once solved for using chi-square minimization, allow opacity correction in the form: $\begin{array}{ccc}& & {\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{t}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{T}}_{\mathrm{sys}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{-}{\mathit{T}}_{\mathrm{atm}}\mathrm{-}{\mathit{T}}_{\mathrm{rx}}^{\mathit{p}}\mathrm{-}{\mathit{T}}_{\mathrm{spill}}}{{\mathit{T}}_{\mathit{a}}^{\mathit{p}}\mathrm{-}{\mathit{T}}_{\mathrm{atm}}}\\ & & {\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{\tau ,t}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathit{\tau }\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}\mathit{\tau }\mathrm{=}\mathrm{0}\mathit{,t}\mathrm{)}\mathit{.}\end{array}$The a priori amplitude calibration methods outlined above are limited in their intrinsic accuracy by several sources of error. These include systematic errors in the measured or adopted power level of the noise calibration sources, intrinsic statistical error in the sampled switched T_sys measurements, and the stability over time of hardware elements affecting amplitude calibration, including receiver gain and the stability of the noise calibration sources. In addition, the opacity solution is limited by the assumption of a stable, plane-parallel atmosphere; this does not hold at low elevations or if local weather conditions vary significantly over the course of the observations (i.e. τ₀ = τ₀(t)). Other sources of systematic error include the assumption of an isothermal T_atm extrapolated from the measured ground air temperature, and the relatively coarse empirical model adopted for antenna spill-over noise contributions.

We estimate the overall accuracy of a priori VLBA calibration to be  ~10% at λ = 7 mm and  ~15% at λ = 3 mm, excluding low-elevation data.

A priori amplitude calibration, with or without opacity correction, can be refined using amplitude self-calibration. Historically, two assumptions are often made when performing amplitude self-calibration using parallel-hand visibility data: i) that Stokes V is identically zero; and ii) that amplitude calibration is completely separable from instrumental polarization calibration. The inapplicability of these assumptions when measuring non-zero fractional circular polarization m_c ≲ 0.5% in compact extragalactic continuum sources is considered in detail by Homan & Wardle (1999), Homan et al. (2001), and Homan & Lister (2006). Reduction methods for spectral-line sources with non-zero m_c are described by Kemball et al. (1995). As described in the latter reference, in the presence of circular polarization, independent calibration of the parallel-hand correlations r^RR(u,v,ω) and r^LL(u,v,ω) against a total intensity source model that assumes V = 0 will redistribute circularly-polarized emission in the image, including introducing positional offsets in the centroids of individual circularly-polarized components. As a result, Kemball et al. (1995) used only a single reference receptor polarization p ∈  {R,L} in phase-related self-calibration, and used measured R − L phase offsets to transfer solutions to the orthogonal receptor polarization. Similarly, amplitude self-calibration is avoided, for the reasons noted above.

Homan & Wardle (1999) introduced a method for calibrating continuum sources with low m_c while still allowing amplitude self-calibration; their method calibrates (r^RR + r^LL)/2 against a common total intensity source model I_mod, then solves for long-term residual differential gain errors at each antenna under the Stokes V = 0 assumption, i.e. r^RR = I_mod and r^LL = I_mod. The latter measurement of the differential polarization gain factors relies on the mean circular polarization of the observed ensemble of sources being close to zero. This approach has been demonstrated successfully at wavelengths longer than λ ≳ 1 cm (Homan et al. 2009). Simulation studies for this technique indicate a 1 − σ uncertainty in measured circular polarization of m_c ~ 0.1% for source brightness values exceeding 1 Jy/beam (Homan & Lister 2006). Uncertainties in amplitude gain calibration predominate over D-term errors or thermal noise contributions in this method (Homan & Lister 2006).

The coupling of amplitude and instrumental polarization is evident in the form of Eq. (15), which contains both G and D terms. Neglecting this coupling is equivalent to linearization of Eq. (15), which truncates terms that include instrumental polarization D-terms from the parallel-hand equations (Roberts et al. 1994; Homan & Wardle 1999). This can be addressed by an iterative approach to amplitude and polarization calibration (Homan & Wardle 1999) but with a modest impact on the derived location and magnitude of circularly-polarized components (Homan & Lister 2006).

3.3.2. Spectral-line amplitude calibration

The template-fitting method derives amplitude calibration from parallel-hand autocorrelation spectra, and was first introduced by Reid et al. (1980). The source total power spectrum usually has far greater frequency structure than the continuum noise terms in Eq. (5), and the two can therefore be treated as sufficiently orthogonal in functional form to allow a basis decomposition in terms of the scaled true source spectrum and a residual continuum term varying more slowly with frequency. The basic method derives time-variable gain normalization factors ${\mathit{G}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ by fitting baseline- and bandpass-corrected parallel-hand autocorrelation spectra ${\mathit{r}}_{\mathit{mm}}^{\mathit{pp}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ to a well-characterized source total power spectrum J^pp(ω) derived from a template scan at a single antenna; this template scan is chosen based on SNR and the quality of available a priori SEFD calibration information for the scan. The method is readily applied to polarization observations by fitting ${\mathit{G}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ separately in each receptor polarization p against the associated parallel-hand template spectra J^pp(ω) and J^qq(ω) (Kemball et al. 1995). These authors also introduced a robust method to fit a composite baseline during the template-fit itself, avoiding the need for prior baseline removal in the source and template spectra individually. Given the uncertainties in a priori SEFD calibration discussed above, a secondary correction is needed to determine and remove the differential polarization amplitude gain g^R/L tying the template spectra in each sense of parallel-hand receptor polarization. In the current paper, g^R/L is used to refer to a ratio of terms  ∥ G^R ∥  and  ∥ G^L ∥  in the Jones matrix of Eq. (9). At centimeter wavelengths, Kemball et al. (1995) solved for this differential polarization gain from relative ratios of the cross-power spectral amplitudes on continuum extragalactic calibrator sources, which were assumed to have Stokes V ~ 0. We note that in the case of maser sources it is necessary to account for the average antenna noise contribution from the source $\mathit{T̅}\begin{array}{}\mathit{p}\\ \mathit{a}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathit{\delta \omega }}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{\delta \omega }}{\mathit{T}}_{\mathit{a}}^{\mathit{p}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{d}\mathit{\omega }$ when transferring the continuum polarization gain ratio to the spectral-line source (Kemball et al. 1995).

The template fitting method has several intrinsic advantages at the wavelengths λ = 7 mm and λ = 3 mm considered in this paper (Kemball & Diamond 1997). The method implicitly includes opacity corrections in the determined amplitude scaling factors, and is therefore insensitive to systematic errors in either the opacity parametrization of Eq. (5) or in the published a priori antenna gain curves. In addition, the method intrinsically tracks short-term variations in the amplitude gain in each receptor polarization at each antenna that result from factors such as pointing errors in the presence of beam squint, or rapid changes in atmospheric attenuation due to highly time-variable local weather conditions. Such short-term amplitude gain fluctuations are significantly more pronounced at wavelengths shorter than  cm. The gain factor information is also encoded in the autocorrelation spectra sampled at the correlator integration rate, as opposed to to a separate sampling interval for switched-noise T_sys measurement. The template-fitting method is more statistically robust at measuring short-term amplitude gain fluctuations, such as those caused by pointing errors, than a method based on detecting associated changes in the measured T_sys values, as the latter requires separating out changes in two unknown continuum quantities, $\mathit{T̅}\begin{array}{}\mathit{p}\\ \mathit{a}\end{array}$ and the baseline noise terms. The template-fitting method exploits spectral structure in this decomposition.

Kemball & Diamond (1997) used two methods at λ = 7 mm to measure the differential polarization amplitude gain g^R/L for the template spectra: i) comparison of amplitude self-calibration gain corrections derived from continuum calibrators assuming Stokes V ~ 0, obtained for the reference antenna near the time of the template spectrum observation; and, ii) cross-fitting the template spectra in each polarization J^pp(ω) and J^qq(ω) to each other using the template-fitting algorithm to derive the relative differential polarization amplitude gain g^R/L directly. Method (ii) assumes that the small non-zero integrated Stokes V component of SiO emission does not significantly bias the estimated g^R/L, i.e. that .

In this paper we introduce several refinements to the template-fitting amplitude calibration method to improve its statistical performance in the millimeter wavelength range λ = 3 mm and λ = 7 mm considered here. These modifications are described in the following section.

3.3.3. Autocorrelation polarization self-calibration

VLBA system performance degrades sharply toward λ = 3 mm, consistent with array design specifications; this requires that we consider enhancements in the statistical performance of the template-fitting method in the low-SNR regime. In addition, the high representative fractional linear polarization of SiO masers m_l ~ 10−30% described earlier suggests that higher-order terms, such as O(D.(Q + jU)), that arise in the full non-linear coupling of amplitude and instrumental polarization calibration contained in Eq. (8), may need to be assessed when using the template-fitting method acting on the parallel-hand autocorrelation spectra.

To explore both concerns, we have implemented an iterative self-calibration method to derive the amplitude gains G(t), instrumental polarization D, and the true source correlation spectrum J(ω) from the measured un-calibrated autocorrelation spectra R(ω,t) at all antennas. The initial estimate of the true source spectrum J(ω) is obtained from an weighted average of baseline-corrected autocorrelation spectra over a restricted range of low zenith angle and over a subset of antennas with high site elevation or known low mean precipitable atmospheric water vapor. The spectra are calibrated a priori by the measured ${\mathit{T}}_{\mathrm{sys}}^{\mathit{p}}$ and opacity-corrected gain curves ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{(}\mathit{\tau }\mathrm{=}\mathrm{0}\mathit{,z}\mathrm{)}$ provided by the VLBA, and then converted to, and averaged separately in Stokes  {I,Q,U,V} form. The spectra are pre-averaged at each antenna over a short interval before baseline subtraction, and weighted by the inverse mean-square residual error of the baseline polynomial fit. Spectra with a completeness fraction below a specified fractional threshold of the pre-average interval are rejected in the global average used as an estimator for J(ω). The averaging of N high-quality spectra from a subset of antennas reduces statistical noise in the resulting estimate of J(ω) by approximately $\sqrt{\mathit{N}}$, an important improvement in the low-SNR regime, offset only by residual systematic error contributed primarily by the baseline removal process. The baseline is modeled as a low-order polynomial and removed in a fit to designated off-source spectral regions. Optionally, auxiliary piece-wise polynomial baselines are fit and removed above and below the outermost off-source regions; this approach minimizes the polynomial order required for the primary spectral baseline. The absolute flux density scale of J(ω) is subject to the uncertainty in a priori VLBA calibration, estimated above to be  ~10−15% at millimeter frequencies, and to the approximation of zero atmospheric opacity τ ~ 0 used in deriving the initial source correlation spectrum. However, our goal here is accurate measurement of the degree of circular polarization m_c; this does not require comparable accuracy in the absolute flux density scale. For this reason we enforce the constraint ${}^{\mathrm{\int }}\mathit{J}^{\mathit{pp}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{d}\mathit{\omega }\mathrm{=}{}^{\mathrm{\int }}\mathit{J}^{\mathit{qq}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{d}\mathit{\omega }$ over the two senses of parallel-hand correlation in the averaged estimate for J(ω), with corresponding scaling in cross-polarized spectra J^pq(ω). To establish the global differential polarization amplitude gain g^R/L connecting the mean template spectra J^pp(ω) and J^qq(ω), we apply the amplitude gains obtained from template-fitting to these spectra to a compact continuum calibrator and solve for g^R/L by minimizing (Kemball et al. 1995): $\sum _{{\mathit{t}}_{\mathit{k}}}\sum _{\mathit{m,n}\mathit{>}\mathit{m}}\sum _{{\mathit{\omega }}_{\mathit{k}}}{\left[\frac{{\mathit{g}}^{\mathit{R}\mathit{/}\mathit{L}}\mathrm{\parallel }{\mathit{V}}_{\mathit{mn}}^{\mathit{LL}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathit{,}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }\mathrm{-}\mathrm{\parallel }{\mathit{V}}_{\mathit{mn}}^{\mathit{RR}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathit{,}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }}{{\mathit{g}}^{\mathit{R}\mathit{/}\mathit{L}}\mathrm{\parallel }{\mathit{V}}_{\mathit{mn}}^{\mathit{LL}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathit{,}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }\mathrm{+}\mathrm{\parallel }{\mathit{V}}_{\mathit{mn}}^{\mathit{RR}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathit{,}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}\mathrm{\parallel }}\right]}^{\mathrm{2}}$(35)The parallel-hand cross-correlation data ${\mathit{V}}_{\mathit{mn}}^{\mathit{pp}}\mathrm{\left(}{\mathit{t}}_{\mathit{k}}\mathit{,}{\mathit{\omega }}_{\mathit{k}}\mathrm{\right)}$ are calibrated in group delay, fringe rate, and by the template-fitted amplitude gains, then pre-averaged over time and frequency sub-intervals in order to suppress outliers. The statistic in Eq. (35) is by design not sensitive to un-modeled calibrator spatial structure. Models of intrinsic circular polarization mechanisms in continuum extra-galactic radio sources (Homan et al. 2009) suggest intrinsic circular polarization m_c ≪ 0.1% at the millimeter wavelengths considered here.

The parallel-hand auto-correlation spectra are real-valued, unlike the complex cross-polarized autocorrelation spectra, which have non-zero phase and form a conjugate pair J^pq(ω) = J^qp ∗ (ω). Using bandpass solutions corrected for the differential polarization bandpass phase response at the reference antenna ${\mathit{\varsigma }}_{\mathrm{0}}^{\mathit{p}\mathrm{-}\mathit{q}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ described earlier, the instrumental phase correction for the autocorrelation cross-polarized spectra (Eq. (13)) at antenna m can be parametrized in the form, $\begin{array}{ccc}\arg \mathrm{\left(}K\mathrm{\right)}& \mathrm{=}& \mathrm{diag}\mathrm{(}\mathrm{0}\mathit{,}\mathrm{\Phi }\mathit{,}\mathrm{-}\mathrm{\Phi }\mathit{,}\mathrm{0}\mathrm{)}\\ \mathrm{\Phi }& \mathrm{=}& {\mathit{\phi }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}\mathrm{+}\mathit{\omega }{\mathit{\tau }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}\mathrm{+}\mathrm{(}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{R}}\mathrm{-}{{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{L}}}^{\mathrm{)}}\end{array}$where ${\mathit{\phi }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ is the differential polarization phase offset, ${\mathit{\tau }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ is the differential polarization group delay, and $\mathrm{(}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{R}}\mathrm{-}{\mathit{\varsigma }}_{\mathit{m}}^{\mathit{L}}\mathrm{)}$ is the bandpass phase correction, all at antenna m. The differential polarization phase offset and group delay at each antenna are assumed constant over the course of the observations. The variation of ${\mathit{\phi }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ by antenna across the array introduces amplitude decorrelation when averaging the autocorrelation spectra to estimate J(ω), as described above. Accordingly, during the averaging process, we normalize the mean phase (over frequency) of each pre-averaged spectrum to zero before it is added to the global average estimate of J(ω). This is equivalent to accepting an unknown additive term in the absolute EVPA of the true source correlation spectrum J(ω); this does not affect our calibration method, and would traditionally be measured using an astronomical calibrator of known absolute EVPA.

As evident from Eqs. (8) and (14) the amplitude (and for the cross-polarized spectra, phase) of the measured autocorrelation spectra R change over both time and frequency due to the combined effect of the polarization leakage terms acting on both the true source correlation spectrum J(ω) and the real-valued self-noise N, with the former terms modulated by parallactic angle variation. Our initial estimate of the source correlation spectrum J(ω) is accordingly imperfect, as it has been corrected only for an empirical estimate of ${\mathit{\phi }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ at each antenna, as described above.

Following traditional self-calibration practice, we use our initial estimate of the source spectrum J(ω) to solve for the instrumental calibration terms. We first hold the initial amplitude calibration in place, and solve for the instrumental polarization terms D at each antenna. This is solved as a non-linear chi-square problem, minimizing the norm $\mathrm{\parallel }R\mathrm{-}{}^{\mathrm{\right(}}\mathrm{KL}^{\mathrm{\left)}}{}^{\mathrm{(}}J\mathrm{+}{N}^{\mathrm{)}}\mathrm{\parallel }$ using full-polarization autocorrelation spectra on the target spectral line source, summed over pre-averaging intervals  △ t over which the parallactic angle does not vary significantly, and over a spectral channel range encompassing the source emission. This problem is completely separable for each antenna, as no cross-correlations between antennas are involved, and is therefore performed separately and sequentially for each antenna. The free parameters over p ∈  {R,L} are the complex ${\mathit{D}}_{\mathit{m}}^{\mathit{p}}$, assumed time- and frequency-invariant, the differential polarization quantities ${\mathit{\phi }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ and ${\mathit{\tau }}_{\mathit{m}}^{\mathit{R}\mathrm{-}\mathit{L}}$ used to specify Φ (Eq. (37)), and the self-noise contributions . The latter are modeled as low-order polynomials over frequency, in order to accommodate residual errors in bandpass correction, with independent coefficients per per-average interval due to time-variable atmospheric terms in Eq. (5). The total number of free real-valued parameters in the chi-square minimization per antenna is of order O(10²) with the greatest contribution from the independent baseline terms in each pre-average interval. Accordingly we use two unconstrained minimization methods suited for the large-parameter case. Given the Gaussian noise in the measured data, we use the limited-memory BFGS quasi-Newton method described by Nocedal (1980) and Liu & Nocedal (1989) to find an initial estimate of the minimum, and a truncated Newton method to refine the solution (Nash 1984). Both methods work well in the large-parameter case.

Given an estimate for ${\mathit{D}}_{\mathit{m}}^{\mathit{p}}$, Φ_m, and obtained from the chi-square solution, we correct the autocorrelation spectra by inverting Eq. (8), as, $J^{\mathrm{\prime }}\mathrm{=}{\mathrm{(}{\mathrm{KL}}^{\mathrm{)}}}^{-1}R\mathrm{-}N\mathit{.}$(38)The corrected J′ = J_m(ω,t) at each antenna is sampled at the original correlation integration interval and original frequency spacing. For perfect instrumental calibration we would expect an identical estimate of the true source spectrum at each antenna and at all times. However, in common with all self-calibration methodologies we only approach this incrementally over succeeding iterations. Accordingly, we update our estimate of the true source correlation spectrum J(ω) by averaging J′(ω) over a high-quality subset of data, exactly as described above in connection with the derivation of the initial estimate. Two sub-iterations of polarization calibration are performed before an outer iteration of amplitude calibration using the template-fitting method. In the amplitude calibration step, the corrected J′(ω) for the current iteration are fit against the current estimate of the true source correlation spectrum J(ω) obtained by averaging followed by an estimate of the differential polarization amplitude gain g^R/L as described using Eq. (35). This method was found to converge rapidly, generally within one or two outer iterations, to provide a joint estimate for the true source correlation spectrum J(ω) and the instrumental parameters $\mathrm{\parallel }{\mathit{G}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\parallel }$, ${\mathit{D}}_{\mathit{m}}^{\mathit{p}}$, Φ_m, and .

3.4. Interferometric phase and polarization calibration

The cross-correlation phase calibration method used follows that described in Kemball et al. (1995); namely phase calibration in a reference receptor polarization, combined with the determination of differential polarization phase offsets between the two orthogonal polarization receptors. Interferometric polarization self-calibration is performed as described by Kemball & Diamond (1997). Following Homan & Wardle (1999), we add a final phase self-calibration step after interferometric polarization calibration.

Fig. 3 Scan-preaveraged parallel-hand (RR) autocorrelation spectra for the two continuum calibrator scans at antenna LA with extremal values of geometric fringe rate over the course of the observing run. Fringe rates are given in units of channel width, and frequency channels 1–128 span the baseband bandwidth of 4 MHz.

4. Results

In this section we present results obtained by applying the data reduction techniques outlined above to representative observational data, chosen here to be a single epoch from a full-polarization VLBA monitoring campaign of the circumstellar v = 1,  J = 1−0 SiO maser emission in the λ = 7 mm band toward the Mira variable TX Cam, scheduled as VLBA project code BD46. The results of this synoptic imaging campaign have been published in total intensity by Diamond & Kemball (2003) and Gonidakis et al. (2010), and in linear polarization by Kemball et al. (2009). We choose epoch BD46AQ as a representative dataset for this work, as it contains polarization EVPA reversals at the circumstellar boundary and is therefore of independent scientific interest, but for no other special technical reason. The data were reduced for this epoch following the heuristics outlined in Sect. 3 and we highlight results relevant to the current data reduction method in this section.

Full details of the observing configuration for the project BD46AQ are described in the original references above; we provide only a concise synopsis here. The effective observing time was 6.5 h, allocated between the target source TX Cam and extra-galactic continuum calibrators 3C454.3, J0359+509, and J0609-157. The v = 1,  J = 1−0 SiO transition was centered in a 4 MHz baseband and cross-correlated in full polarization over 128 frequency channels with a correlator accumulation interval of 4.99 s. The array comprised all ten VLBA antennas and a single antenna from the VLA. The VLBA antenna at Mauna Kea (MK) did not contribute data due to operational difficulties.

The data were sampled in one-bit quantization; digital corrections for sampling and quantization where derived using methods described above. The derived quantization corrections have low magnitude  ≤ 0.5% and are stable in time, as expected for one-bit quantization. After these corrections, the antenna-based autocorrelation bandpass amplitude responses were solved for in Chebyshev polynomial form using the methods of Sect. 3.2. The bandpass frequency-frame considerations discussed in that section are illustrated in Fig. 3; this figure shows two parallel-hand continuum calibrator autocorrelation spectra produced by the VLBA correlator for a single antenna, here Los Alamos (LA), pre-averaged over time in the bandpass solver. These spectra are in a geocentric reference frame, and this figure shows the two scans that have extremal natural fringe-rate offsets between the topocentric and geocentric frequency reference frames over the course of the observation. As noted in Sect. 2, the polynomial bandpass solver accommodates the frequency shift algebraically in the Chebyshev expansion for the bandpass response function, without needing to interpolate the sampled data in frequency, and can use fewer free parameters than the total number of frequency channels N_c as a result. The matching full topocentric channel range needs to be used in the fit however, as opposed to  [1..N_c] , in order to provide a maximal constraint on χ² across the band. An expanded view of the frequency shifts in the scan pre-averaged data near the upper spectral roll-off is shown in Fig. 4; this depicts a subset of scans spanning the full range of natural fringe rate in the data for the Los Alamos antenna.

Fig. 4 Scan-preaveraged parallel-hand (RR) autocorrelation spectra in the upper spectral roll-off region of the baseband for a subset of continuum calibrator scans at antenna LA spanning the range of geometric fringe-rate over the course of the observing run. Fringe rates are given in units of channel width, and frequency channels 1–128 span the baseband bandwidth of 4 MHz.

As described in Sect. 3.2, a correction for aliasing is required if the high-SNR autocorrelation bandpass amplitude response functions are to be used to correct the cross-correlation data. The results of the aliasing correction are presented for a selection of representative antennas in Fig. 5; here shown in the RCP bandpass response functions for this antenna subset. For each antenna, a separate amplitude response is plotted for: i) the original autocorrelation bandpass; ii) a cross-correlation bandpass estimate; and; iii) the alias-corrected response. Note the low SNR of the cross-correlation bandpass response obtained from Eq. (25), for reasons described in that section.

The alias-corrected amplitude response preserves the SNR of the autocorrelation bandpass solution while correcting both the mid-band mean amplitude and upper band-edge frequency response, so that it can be used to correct cross-correlation data more accurately. This improvement in performance when correcting cross-correlation amplitudes is shown in Fig. 8. The final component of bandpass determination is the correction for the differential polarization phase response at the reference antenna, derived from the cross-polarized autocorrelation spectra as described in Sect. 3. The R-L bandpass phase response as solved for from the RL and LR autocorrelation spectra at the reference antenna LA is shown in Fig. 6. The corrected cross-polarized RL autocorrelation spectrum at LA, after correction for the R-L bandpass phase response is plotted in Fig. 7.

The data were then calibrated in amplitude using the autocorrelation polarization self-calibration algorithm described in Sect. 3.3.3. This method refines the standard template-fitting method to allow the use of a mean multi-antenna template spectrum with higher SNR, depicted here in Fig. 9 for both parallel hands. The method also provides an independent calibration of the differential polarization amplitude gain g^R/L (Eq. (35)), and uses the full non-linear autocorrelation polarization data model including self-noise, as described in Eq. (8). An example of fitting this equation to the autocorrelation spectra at the reference antenna LA is shown in Fig. 10. As described earlier, the chi-square minimization in this fit is performed over all Stokes autocorrelation polarization pairs; however for clarity of presentation we show the real-valued Stokes Q spectra pre-averaged over each scan in the form of both the original data and the fitted model. Each scan is annotated by the corresponding parallactic angle. The time-variability of the Stokes Q spectra across the observation arises from both the parallactic angle terms in Eq. (8) as well as the D-term mediated corruptions from the autocorrelation spectra in other Stokes polarization pairs. As described in Sect. 3.3.3, an iterative polarization self-calibration method is used to refine both the true source model correlation spectrum and the instrumental terms, including D-terms and the differential polarization phase and group-delay offsets. In Fig. 11 we show the rapid iterative convergence of the RL source model correlation spectrum over three iterations: the initial and final model spectra are plotted here.

Fig. 5 Bandpass amplitude response function solutions in RCP for a subset of representative antennas obtained from: i) autocorrelation calibrator data (dotted line) by solving Eq. (19); ii) cross-correlation calibrator data (dashed line) by solving Eq. (25); and iii) the autocorrelation bandpass solutions corrected for a fitted aliasing model (solid line), as defined by Eq. (26).

As a test of net error in Stokes V the amplitude gain factors derived from TX Cam using the preceding spectral-line calibration sequence were scaled by an additional constant user-specified R/L amplitude gain factor and applied to the continuum calibrator J0359+509. The J0359+509 data were then self-calibrated in phase only, an independent D-term solution derived using the similarity polarization approximation (Kemball 1999), and imaged in Stokes  {I,Q,U,V}. These images were produced over a range of such user-specified R/L amplitude gain factors. This particular continuum calibrator was observed most frequently during the schedule, and amongst the calibrators has the smallest angular separation from TX Cam on the sky. The integrated degree of circular polarization m_c was then measured by integrating Stokes V and Stokes I over a rectangular region tightly enclosing the central emission in the calibrator image of J0359+509. The resultant plot of m_c versus the additional user-specified R/L gain factor is shown in Fig. 12. The measured fractional circular polarization of the continuum calibrator for unit additional R/L amplitude gain factor provides an estimate of the error in the data reduction algorithm presented here. A linear fit produces a value m_c = −0.3% at the nominal value where the user-specified R/L factor is unity.

Fig. 6 R-L differential polarization bandpass phase response at the reference antenna LA, as measured from the cross-polarized RL and LR autocorrelation spectra. The amplitude response shown here is derived from the parallel-hand autocorrelation data.

Fig. 7 RL autocorrelation spectrum at the reference antenna LA, corrected by the bandpass depicted in Fig. 6.

Fig. 8 Continuum calibrator cross-power spectral amplitudes on representative baselines to the reference antenna LA, calibrated in amplitude, group delay, and fringe rate, and averaged over the duration of the observing run; calibrated additionally either with the original autocorrelation bandpass response (solid line) or the alias-corrected bandpass response (dotted line).

Fig. 9 Mean autocorrelation template spectra in RR and LL, plotted on a common y-axis of spectral flux density (Jy). Frequency channels 1–128 span the 4 MHz baseband bandwidth.

5. Discussion

General-purpose radio-interferometric arrays are optimized in design to provide radiometric accuracy matched primarily to the requirements of their key science goals. The a priori amplitude calibration information provided by the VLBA, in the form of measured system temperatures ${\mathit{T}}_{\mathrm{sys}\mathit{,m}}^{\mathit{p}}$ and point source sensitivities ${\mathit{P}}_{\mathit{m}}^{\mathit{p}}\mathrm{\left(}\mathit{z}\mathrm{\right)}$ is not the limiting final accuracy of amplitude calibration; in many cases it can be refined to greater accuracy (<5%) using careful subsequent amplitude self-calibration; and the absolute flux density scale can be set to comparable accuracy using compact amplitude calibrators of known brightness.

Fig. 10 Scan pre-averaged autocorrelation Stokes Q spectra before polarization self-calibration (dotted line) and as fitted against Eq. (8) including polarization terms (solid line). Each scan pre-average is labeled at left by parallactic angle α (in deg).

Fig. 11 Soure RL correlation total-power spectra for iteration one (S1) and three (S3) of the autocorrelation polarization self-calibration algorithm.

Fig. 12 Measured circular polarization mc of the continuum calibrator J0359+509 after application of the TX Cam spectral–line amplitude calibration and an additional user-specified R/L amplitude gain factor. The continuum calibrator circular polarization measurement for unit additional R/L gain factor provides a single measurement of the error in mc for the data reduction algorithm presented here. This intercept value is mc = −0.3%.

Amplitude self-calibration refinement in this broader context, traditionally applied to continuum parallel-hand data, requires imposition of an external constraint on Stokes V however. Invariably, Stokes V ~ 0 is assumed for individual continuum calibrator or target sources. However, as noted above, the more selective constraint that an ensemble of continuum target sources should have mean (Homan & Wardle 1999) allows non-zero m_c to be measured for individual sources to higher accuracy. For spectral-line sources, an analogous assumption of zero mean circular polarization across frequency could be made for an individual source, however in the latter case this must always be in essence an ad hoc assumption. The mean net for a given SiO maser source depends on a number of unknown physical properties of the source, e.g. magnetic field orientation and the morphological distribution of maser components across the source.

In the current paper we have shown that it is possible to measure circular polarization in spectral-line VLBI observations at millimeter wavelengths to a far higher level of accuracy (<1%) than suggested by the absolute accuracy of the a priori amplitude calibration information (10−15%), but without the need to perform amplitude self-calibration. This is possible by relying instead on the amplitude gain information encoded in the high-SNR autocorrelation spectra, which are additionally constrained to measure a common (antenna-independent) source correlation spectrum J(ω). This principle was recognized early in the development of the template-fitting method of amplitude calibration for spectral line VLBI (Reid et al. 1980); in the current work we have shown that it is possible to enhance this approach sufficiently to extend it to accurate circular polarization measurements at millimeter observing wavelengths. This approach has the advantage of avoiding amplitude self-calibration, which is statistically less robust when solving for antenna amplitude gain factors across sparse arrays when observing sources with complex spatial structure.

In Fig. 13, we show a plot of the polarization ratio of template-fitted amplitude gains for each antenna in each of the two orthogonal receptor polarizations against the value of the RR template-fitted gain derived as part of the same calibration. The amplitude gain factors are expected to increase with increasing atmospheric airmass along the line of sight or with increasing pointing error. These two effects are also partially coupled due to the expected increase in residual pointing rms at low elevation for mechanical reasons. The structure in this figure, especially when examined for individual antennas, is consistent with differential polarization amplitude gain effects due to pointing errors in the presence of beam squint (see Eq. (3)), and suggests the template-fitting method is correcting the associated amplitude gain effects. This effect is more pronounced at 3 mm, as expected.

Fig. 13 Parallel-hand template-fitted amplitude gain factors for each antenna derived from the mean template spectra (Fig. 9), plotted as differential polarization amplitude gain ratios against the RR amplitude gain derived by the same method.

We next consider the applicability of the algorithm we have presented in this paper. It does both require sufficient SNR in the autocorrelation spectra, and adequate spectral structure in the source emission so that the source spectrum is sufficiently mathematically non-degenerate relative to continuum baseline terms of polynomial form. However, these conditions are broadly met in observations of strong maser emission in these bands.

Our implementation of autocorrelation polarization self-calibration here has allowed an assessment the impact of second-order instrumental polarization terms on the parallel-hand autocorrelation spectra, and on the resultant amplitude gain factors estimated using the revised template-fitting method. Our conclusion from the current study is that this correction is helpful if the autocorrelation data have sufficient SNR and an adequate range of parallactic angle coverage. The SNR condition typically holds at λ = 7 mm but may not be routinely applicable at λ = 3 mm, due to the significantly higher SEFD in the latter band. Given the low D-term magnitude at the VLBA, in the low-SNR case the autocorrelation polarization self-calibration step may not be of net benefit and can be omitted.

Where applicable in terms of SNR however, we find the autocorrelation polarization method to be strongly convergent, especially when deriving an initial estimate of the cross-polarized autocorrelation source model using zero mean phase averaging, as noted above. This convergence is apparent in Fig. 11.

We note that the computational cost of the algorithms described in this paper is modest, broadly comparable to regular self-calibration.

Further work would be beneficial in several key areas. An exploration is needed of optimal polynomial bases for the expansion of complex bandpass responses, as a possible alternative to the Chebyshev basis used here, and especially for the low-SNR regime at these observing wavelengths. In some low-SNR cases the Chebyshev expansion will over-fit the noise and introduce low-level quasi-sinusoidal ripple. Opimized spectral baseline removal methods for the low-SNR case would likely improve both the template-fitting and autocorrelation polarization self-calibration elements. In addition, and most importantly, the data reduction technique described here needs to be applied to a greater sample of sources in order to set both the required SNR limit more precisely and to obtain a more comprehensive measurement of the statistical performance of the algorithm as an estimator for m_c and Stokes V over a greater source sample.

From the current work, we estimate the accuracy of the full algorithm presented in this paper as m_c ≤ 0.5% or better at λ = 7 mm, and m_c ≤ 0.5−1% or better at λ = 3 mm. We base this estimate on measurements of Stokes V in associated continuum calibrators, as shown in Fig. 12 for the test data BD46AQ considered here. These figures are also broadly consistent with sub-sampling bootstrap resampling tests on the variance of g^R/L estimated using Eq. (35). We note that both these tests tend to over-estimate the error in m_c as both are affected by interpolation errors from the position of the target source to the position of the continuum calibrator. These interpolation errors include differences in pointing errors and airmass along the lines of sight to the two separate sources. As noted above, further tests with larger source samples would be helpful in this regard.

The technique presented here however has significant scientific application in the study of circular polarization of astrophysical masers in these observing bands. This is especially true for SiO masers, for which accurate Stokes V measurements provide strong constraints on models of polarized maser propagation, and correspondingly of associated estimates of astrophysical magnetic fields.

6. Conclusions

The conclusions of our work are as follows:

• 1.

  We have examined the sources of calibration error in spectral-line VLBI imaging observations at λ = 7 mm and λ = 3 mm sensitive to circular polarization. This analysis was performed with a specific emphasis on SiO maser observations using the VLBA.

• 2.

  A algorithm is presented to provide accurate calibration of circular polarization without using amplitude self-calibration. This method is an enhancement of existing spectral-line VLBI calibration methods based on autocorrelation data, but with several innovations in bandpass estimation, autocorrelation polarization self-calibration, and the adaption of techniques for the low-SNR regime applicable at millimeter wavelengths.

• 3.

  We demonstrate an example reduction at λ = 7 mm and provide an estimate of circular polarization accuracy of m_c ≤ 0.5% or better at λ = 7 mm and m_c ≤ 0.5−1% or better at 3 mm. These estimates are based on Stokes V imaging of associated continuum calibrators and a statistical analysis of uv-data differential polarization amplitude ratios.

Acknowledgments

We thank the anonymous referee for their insightful comments on this paper. This material is based upon work partially supported by the National Science Foundation under grant AST-0507473. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

References

All Tables

All Figures

Fig. 1 Plot of normalized primary beam patterns in 1-d cross-section for receptor polarizations  {R,L} with a beam squint equal to 5% of the full-width at half-maximum θb; the x-ordinate for the beam pattern plot is in units of θb. The ratio of the two beam patterns is also plotted in the upper part of the diagram as a function of pointing error ϵp, with x-ordinate $\frac{{\mathit{ϵ}}_{\mathrm{p}}}{{\mathit{\theta }}_{\mathit{b}}}$).

In the text

	Fig. 2 Flowchart depicting the overall data reduction sequence described in Sect. 3.
In the text

	Fig. 3 Scan-preaveraged parallel-hand (RR) autocorrelation spectra for the two continuum calibrator scans at antenna LA with extremal values of geometric fringe rate over the course of the observing run. Fringe rates are given in units of channel width, and frequency channels 1–128 span the baseband bandwidth of 4 MHz.
In the text

	Fig. 4 Scan-preaveraged parallel-hand (RR) autocorrelation spectra in the upper spectral roll-off region of the baseband for a subset of continuum calibrator scans at antenna LA spanning the range of geometric fringe-rate over the course of the observing run. Fringe rates are given in units of channel width, and frequency channels 1–128 span the baseband bandwidth of 4 MHz.
In the text

	Fig. 5 Bandpass amplitude response function solutions in RCP for a subset of representative antennas obtained from: i) autocorrelation calibrator data (dotted line) by solving Eq. (19); ii) cross-correlation calibrator data (dashed line) by solving Eq. (25); and iii) the autocorrelation bandpass solutions corrected for a fitted aliasing model (solid line), as defined by Eq. (26).
In the text

	Fig. 6 R-L differential polarization bandpass phase response at the reference antenna LA, as measured from the cross-polarized RL and LR autocorrelation spectra. The amplitude response shown here is derived from the parallel-hand autocorrelation data.
In the text

	Fig. 7 RL autocorrelation spectrum at the reference antenna LA, corrected by the bandpass depicted in Fig. 6.
In the text

	Fig. 8 Continuum calibrator cross-power spectral amplitudes on representative baselines to the reference antenna LA, calibrated in amplitude, group delay, and fringe rate, and averaged over the duration of the observing run; calibrated additionally either with the original autocorrelation bandpass response (solid line) or the alias-corrected bandpass response (dotted line).
In the text

	Fig. 9 Mean autocorrelation template spectra in RR and LL, plotted on a common y-axis of spectral flux density (Jy). Frequency channels 1–128 span the 4 MHz baseband bandwidth.
In the text

	Fig. 10 Scan pre-averaged autocorrelation Stokes Q spectra before polarization self-calibration (dotted line) and as fitted against Eq. (8) including polarization terms (solid line). Each scan pre-average is labeled at left by parallactic angle α (in deg).
In the text

	Fig. 11 Soure RL correlation total-power spectra for iteration one (S1) and three (S3) of the autocorrelation polarization self-calibration algorithm.
In the text

Fig. 12 Measured circular polarization mc of the continuum calibrator J0359+509 after application of the TX Cam spectral–line amplitude calibration and an additional user-specified R/L amplitude gain factor. The continuum calibrator circular polarization measurement for unit additional R/L gain factor provides a single measurement of the error in mc for the data reduction algorithm presented here. This intercept value is mc = −0.3%.

In the text

	Fig. 13 Parallel-hand template-fitted amplitude gain factors for each antenna derived from the mean template spectra (Fig. 9), plotted as differential polarization amplitude gain ratios against the RR amplitude gain derived by the same method.
In the text