2. Observations and correlation

The pulsar astrometric analysis presented here uses two VLBA data sets, project codes: (1) BD174, which is being analysed and presented for the first time, and (2) BD152, presented earlier in Deller et al. (Reference Deller2019). The BD152 pulsar astrometric campaign observed 57 pulsars at 1.6 GHz in 8/9 epochs (similar to a typical VLBI astrometric campaign) spanning more than two years (from January 2011 to December 2013). BD174 extends this with another more than 20 observations of two bright pulsars, B0329+54 and B1133+16, lasting less than a year. The BD174 observing parameters are presented in Table 1. These observations were designed to better understand the contributions of systematic uncertainties to astrometric parameters, and were conducted in ‘filler’ time, with short observing duration and (sometimes) fewer antennas available.

Table 1. Details of the BD174 observations.

The following description of observing strategy and parameters pertains only to the BD174 data. Each observation is broken into four blocks, where a block consists of two scans on the target field (each of five minutes) interleaved and bracketed by scans on the Phase Reference Calibrator (PRC; each one minute), which are followed by a Fringe Finder Calibrator (FFC; one minute) scan. The first block is at L band, followed by an S band block, a second L band block, and a second S band block (Figure 1). In total, $4 \times 5$ min is thus obtained on the target field at each band. We have excluded the S-band data from the analysis since it has root-mean-square (RMS) noise almost ten times higher than the theoretical expectation due to the huge radio frequency interference (RFI) in this band (for details, see VLBA Scientific Memo 41Footnote ^b ). Consequently, our analysis focuses solely on the L-band data. The target field, with the pulsar positioned at the pointing centre, included three In-Beam Calibrators (IBCs). Details of observed pulsars with their period, DM, gated equivalent flux density at 1.4 GHz, and calibrators with their separation from the pulsar, along with their flux densities at 1.4 and 1.6 GHz for calibrators, are listed in Table 2.

Figure 1. Illustration of the BD174 observing strategy. Each observation is divided into four blocks, two at L-band and two at S-band, as shown by rectangular boxes. Within each block, two target scans (labelled as ‘T’) are observed in an interleaved manner, with each target scan bracketed by scans of the Phase Reference Calibrator (labelled as ‘C’).

Table 2. Parameters of target pulsars, as well as their associated bandpass and phase reference calibrators.

¹ https://www.atnf.csiro.au/research/pulsar/psrcat/; (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005).

² Deller et al. (Reference Deller2019).

³ Angular separation between the target and the PRC.

⁴ Angular separation between the target and the IBC.

⁵ Flux densities estimated using JMFIT on global model images of the PRC and IBC, obtained by combining all epochs.

The data were correlated using DiFX software correlator (Deller et al. Reference Deller, Tingay, Bailes and West2007, Reference Deller2011). Each observing epoch is processed through five correlation passes: three for IBCs and two for pulsar, one in ungated and the other in gated pulsar mode.Footnote ^c For pulsar gating, the DiFX software correlator utilises a ‘polyco’ file, which contains pulse longitude information as a function of frequency and observation time. For this astrometric campaign, the ‘polyco’ files were generated utilising timing data obtained using the Lovell Telescope. The first IBC pass also includes processing of PRC and FFC scans. The correlation positions of the PRCs, listed in Table 3, are taken from the ‘.sum’ filesFootnote ^d . The table also includes the updated positions of these sources, taken from the Radio Fundamental Catalogue (RFC; Petrov & Kovalev Reference Petrov and Kovalev2024). Our PRCs and IBCs do not have obvious calibrator structural variability on the scale of a significant fraction of the synthesised beam or larger.

Table 3. Correlation and RFC positions of the PRCs.

¹ PRCs coordinates are taken from ‘.sum’ files.

² PRCs coordinates are taken from RFC (Petrov & Kovalev Reference Petrov and Kovalev2024).

3. Data reduction and imaging

The visibility data, resulting from DiFX processing, have been analysed primarily using AIPS (version 31DEC23; Greisen Reference Greisen and Heck2003), utilising ParselTongue (Kettenis et al. Reference Kettenis, van Langevelde, Reynolds, Cotton, Gabriel, Arviset, Ponz and Enrique2006), a Python interface, that facilitates AIPS functionality. The calibration pipeline, VLBI Data Calibration and Imaging pipeline in Python (VDCIpy),Footnote ^e is used in the present work (similar to Deller et al. Reference Deller2019). A brief summary of the calibration pipeline, including AIPS tasks, is as follows. It first performs a priori calibration, including ionospheric dispersive delay and Faraday rotation corrections using TEC maps (TECOR) and adjustments for the Earth orientation parameters using the latest Earth Rotation Parameter (ERP) file (CLCOR). The visibility amplitude corrections are applied through ACCOR and APCAL, which do the quantisation correction and convert the cross-correlation coefficients to Jy-unit, respectively. The primary beam corrections are performed using a ParselTongue script, which compensates for primary beam attenuation and beam squint (Middelberg et al. Reference Middelberg2011). The single-band delay, phase, and bandpass calibration solutions are derived using scans of the calibrators, listed in Table 2. Bandpass and time-independent single-band delay solutions are determined on the FFC scan using the BPASS and FRING tasks, respectively, and applied to all sources. The fringe-fitting solutions are computed on each band and on each scan of the PRC. Phase and amplitude self-calibration solutions are also obtained from PRC scans using the CALIB task. The calibration solutions obtained on the PRC are applied to both the PRC and target field sources. A calibrator model is used while performing the BPASS/FRING/CALIB task, which is generated by concatenating all epochs at the respective frequency bands.

To minimise the systematic errors introduced by the difference in the atmosphere sampled by the target and the PRC, we employ the inverse in-beam phase referencing technique (Imai et al. Reference Imai, Sakai, Nakanishi, Sakanoue, Honma and Miyaji2012; Deller et al. Reference Deller2019) since the target pulsars are adequately bright and can serve as the secondary phase calibrator. This technique determines the residual delays on the pulsar by incorporating the phase self-calibration (CALIB task, with a pulsar model constructed after primary calibration) and applies them to the IBC. To keep the residual delays at a minimum and avoid phase wraps, we have shifted the pulsar model position at each epoch based on the preliminary model of the pulsar motion. The preliminary model accounts for the pulsar position shift due to its proper motion, which can be taken from pulsar timing measurements or previous astrometric studies. The model position of the pulsar at $i{\rm th}$ epoch is $\text{P}_i^j\ =\ \text{P}_0^j\ +\ \mu_j t$ , where $\text{P}_0$ is the position at the model epoch, t is the time since the model epoch, and j represents either R.A. or Dec.

For each pulsar, the first epoch of the BD174 campaign is excluded from the analysis since it does not have the FFC scan. To study the frequency-dependent effects, we have analysed the 1.4 and 1.6 GHz data sets independently. For the 1.6 GHz data, the BD152 calibrator models are used in calibration, as both are observed at nearly the same central frequency (1.651 GHz versus 1.646 GHz). The average beam sizes during the BD174 observations of B0329+54 and B1133+16 are mentioned in Table 4. At some epochs, the beam is even broader, especially when a few antennas are missing during the observation or flagged in the calibration process. This can be due to the RFI or poor bandpass solutions at that antenna.

Table 4. Average synthesised beam size ( $\theta$ ) during the BD174 observations at 1.4 and 1.6 GHz.

The frequency-averaged calibrated visibility data for all sources are stored. Particularly, the in-beam calibrator u-v data is divided by its model to convert it to a pseudo-point source to improve the S/N. After that, the Stokes I visibility data for all sources are Fourier transformed using DIFMAP (version 2.5p; Shepherd Reference Shepherd, Hunt and Payne1997). The source images have been created using natural weighting, and are of a size of $1024 \times 1024$ pixels with each pixel size of 0.75 and 0.85 mas for 1.6 and 1.4 GHz data sets, respectively. The clean image of the IBC at each epoch is loaded into AIPS to fit an elliptical Gaussian centred at the peak using the JMFIT task. The position of the source can be determined by fitting a model directly to the visibility data. However, the reliability of the position uncertainties in this approach depends heavily on the scaling of the visibility weights, which may not fully account for calibration errors. On the other hand, image-plane fitting is insensitive to any constant scale factor errors in the weighting. For this reason (and to provide consistency with previous analyses of this VLBI data on these pulsars), we opted for image-domain fitting over visibility-domain fitting.

In the case of inverse phase referencing, we get a position-time series of the IBC, which is converted to the pulsar position-time series by subtracting the IBC position offset (P $^i_{\text{IBC}}$ - P $^0_{\text{IBC}}$ ; IBC position at $i{\rm th}$ epoch minus the mean position of the IBC) from the pulsar model position, and stored in pmpar.in.preliminary file.Footnote ^f The mean IBC position is estimated from its position-time series, which is obtained from the calibrated IBC data set, without applying pulsar self-calibration solutions. The pmpar file includes only the thermal uncertainties to the IBC positions. The pulsar model position uncertainty should also be added in quadrature, but since our pulsars are bright, their contribution to the overall positional uncertainties is negligible. Just for completeness, one can also fit astrometric parameters directly to the IBC position-time series, and can then force the in-beam source to have zero proper motion and parallax, which corresponds to subtracting its proper motion and parallax from the a priori pulsar model.

4. Astrometric fitting

For each pulsar, three position-time series are obtained: two from the BD174 data set at 1.4 and 1.6 GHz and one from the BD152 data set (Figure 2). These position-time series are phase-referenced to the nearest IBC (Table 2) to minimise ionospheric effects. However, the reference sources –IBCs, typically Active Galactic Nuclei (AGNs) – exhibit a frequency-dependent position shift known as the core-shift (Bartel et al. Reference Bartel, Herring, Ratner, Shapiro and Corey1986; Lobanov Reference Lobanov1998). If a constant source model is assumed, this shift introduces a systematic offset between the 1.4 and 1.6 GHz position-time series and propagates to pulsar positions. More generally, any difference between the actual and the modelled angular brightness distribution of the calibrator source at any frequency introduces a position shift in the calibration terms at that frequency that propagates to the pulsar position.

Figure 2. Position evolution of B0329+54 (left two plots) and B1133+16 (right two plots), including error bars solely based on thermal noise. Data from BD152 are shown in blue, while BD174 data are represented by black (1.4 GHz) and green (1.6 GHz). The plotted positions include the contributions from both the proper motion and the parallax. The best-fit reference epoch (MJD: 56300) position is subtracted from the pulsar position-time series.

To correct for this misalignment, we adjust the 1.4 GHz position-time series by subtracting the reference epoch positions offset (P $_{\text{ref}}^{1.4}$ - P $_{\text{ref}}^{1.6}$ ) from it. This approach is justified because 1.6 GHz astrometry generally provides better position accuracy due to reduced ionospheric distortions and improved resolution. The reference epoch is chosen in the middle of the astrometric campaign to minimise the propagation of proper motion uncertainties to the reference epoch position. The reference epoch positions of B0329+54 and B1136+16 at 1.4 and 1.6 GHz differ by approximately 2.2 and 2.3 mas, respectively. This position difference arises due to the independent phase referencing at each frequency. Once aligned, the corrected BD174 position-time series are merged with BD152 data to determine the final astrometric parameters. Since the 1.6 GHz and BD152 data sets were observed at nearly the same frequency and calibrated using identical calibrator models, as expected, no significant systematic frequency-dependent offset was detected. Nevertheless, we applied the same merging procedure as used for 1.4 and 1.6 GHz. The time offset between the BD174 and BD152 observations implies that any systematic frequency-dependent offset between the calibrator model used in BD152 and that in the 1.6 GHz observations would mostly be absorbed into the proper motion fit. However, since the 1.6 GHz and BD152 data sets were observed at nearly the same frequency, using an identical model is unlikely to introduce significant unmodeled core-shift errors.

During the astrometric parameters fitting, we have excluded a few epochs, particularly one epoch (ae) for B0329+54 and three epochs (bk, bl, and bm) for B1133+16, because they exhibit significant deviations from the model (described by Equation 1). Although we could not identify the exact reason, but as a part of our diagnostic process, we have calculated the summed chi-square values ( $\displaystyle\sum_{i=1}^{3} \dfrac{(\text{P}_i - \text{P}_{\text{model}})^2}{\sigma_i^2}$ , $\text{P}_{\text{model}}$ is expressed by Equation 1) for each epoch, considering all three IBCs. The summed chi-square values are around 40, 60, 250, and 130 for the mentioned epochs in that order, while the average chi-square value for the rest of the epochs is around five. The very large number of epochs in the BD174 observing campaign makes it feasible to exclude a few epochs, while it would generally not be feasible for a campaign with a typical number of epochs, $\sim10$ . The final position-time series from BD174 and BD152 data sets (phase-referenced to the closest IBC), excluding the mentioned epochs, are shown in Figure 2.

We have employed three statistical inferences to estimate five astrometric parameters (position, $\alpha_{\text{J2000}}$ , $\delta_{\text{J2000}}$ , proper motion in R.A. and Dec., $\mu_{\alpha}$ , $\mu_{\delta}$ , and parallax, $\varpi$ ) from the concatenated position-time series.

4.1. Least-square fitting

The PMPAR Footnote ^g (version 1.8) software is used to estimate astrometric parameters that are best fit to the position-time series of our pulsars based on a model of the apparent path of a star, which can be modelled as,

(1) \begin{equation}\left\{\begin{aligned}\alpha(t) &= \alpha_{\text{J2000}} + \mu_{\alpha t} + \varpi f_{\alpha(t;} \alpha_{\text{J2000}}, \delta_{\text{J2000}}) \\\delta(t) &= \delta_{\text{J2000}} + \mu_{\delta t }+ \varpi f_{\delta(t;} \alpha_{\text{J2000}}, \delta_{\text{J2000}}),\end{aligned}\right.\end{equation}

where t is the time since the reference epoch, ( $\alpha_{\text{J2000}}$ , $\delta_{\text{J2000}}$ ) are the position of the pulsar at that epoch, and $f_{\alpha}$ , $f_{\delta}$ represent the components of the apparent elliptical path associated with parallax in units of 1 arcsec (which is the parallax value for a source at one pc distance) in R.A. and Dec., respectively,

\[\left\{\begin{aligned}f_{\alpha} &= \dfrac{1}{15}\sec\delta (X\sin\alpha - Y\cos\alpha) \\f_{\delta} &= X \cos\alpha \sin\delta - Y \sin\alpha \sin\delta - Z \cos\delta,\end{aligned}\right.\]

where X, Y, and Z are the epoch-dependent coordinates of the Earth’s centre in the equatorial coordinate system, with its origin at the Solar system barycentre (Green Reference Green1985).

The parameter estimation using the least-squares method is always strongly influenced by the uncertainties in the input data. The reduced chi-square for astrometric parameters fitting generally exceeds unity (for example, the reduced chi-square is eleven for the data set on B0329+54), indicating a significant underestimation of position uncertainties. In the JMFIT task, the estimation of position uncertainties is based solely on thermal noise in the image. However, the total error budget should include systematic errors like atmospheric contributions, calibrator structural evolution, if any, etc. The atmosphere relevant to us comprises the troposphere and ionosphere, which introduce the nondispersive ( $\propto\nu$ ) and dispersive ( $\propto\nu^{-1}$ ) phase errors, respectively. The dominant systematic error contribution in L-band astrometry comes from the ionosphere (Chatterjee et al. Reference Chatterjee, Cordes, Vlemmings, Arzoumanian, Goss and Lazio2004; Brisken et al. Reference Brisken, Benson, Beasley, Fomalont, Goss and Thorsett2000), which mainly depends on the separation between the phase calibrator and the target (Chatterjee et al. Reference Chatterjee, Cordes, Vlemmings, Arzoumanian, Goss and Lazio2004; Deller et al. Reference Deller2019), emphasising the importance of choosing a suitable reference frame (IBC) close to the target.

Table 5. Astrometric parameters obtained using all three statistical approaches. For comparison, Deller et al. (Reference Deller2019) estimates for these pulsars are also included. To facilitate direct comparison, the position of each pulsar at the reference epoch (MJD: 56000) has been aligned with that of Deller et al. (Reference Deller2019) and is therefore not listed in this table.

Deller et al. (Reference Deller2019) have proposed the following empirical relation to compute the systematic error based on calibrator-target separation (s; in arcminutes), antenna elevation (el) during the observation, and S/N on the calibrator (S),

(2) \begin{equation} \sigma_{\text{sys}} = A \times \dfrac{s \times \sum_{a=1}^{N_{\textrm{ant}}} \sum_{o=1}^{M} \csc(el_{a, o})}{M \times N_{\textrm{ ant}}} + \dfrac{B}{S}, \end{equation}

here $\sigma_{\text{sys}}$ is the systematic error in units of the width of the synthesised beam; $el_{a,o}$ is the elevation of antenna a in scan o on the target pulsar; $N_{\text{ant}}$ and M are the number of antennas and number of target scans, respectively and coefficients A and B are 0.001 and 0.6 respectively. These coefficients are estimated using PSRPI observations (Deller et al. Reference Deller2019). Therefore, the empirical relation provides the best estimate of the systematic errors at a mean frequency of 1.6 GHz and the typical ionospheric condition similar to that time. We estimate the systematic error using Equation (2), add it to the random error in quadrature, and store it in pmpar.in file.Footnote ^f The position uncertainty estimates then become more realistic, and the reduced chi-square, as expected, approaches unity. We find that in the case of B0329+54, the reduced chi-square turns out to be unity, while for B1133+16, the reduced chi-square is improved from twelve to three. The least-squares fitting estimates for our pulsars are mentioned in Table 5.

4.2. Bootstrap fitting

As noted above, the least-squares fitting is not the optimal technique to estimate astrometric parameters, particularly their uncertainties, when uncertainties in input data are not properly accounted. In astrometry, systematic error contributions to the uncertainties in the measured positions are well appreciated but difficult to quantify. Hence, we appeal to the bootstrap fitting method (Efron & Tibshirani Reference Efron and Tibshirani1991), which is less sensitive to input uncertainties and thus provides a more reliable estimate of the fitting parameters and their uncertainties.

Here, a large number of bootstrap sample sets (100 000; each of size equal to observation epochs, N $_e$ ) are generated by randomly choosing any of the epoch entries with replacement from the given input data set, where some of the epoch entries from the original set might repeat and/or be absent. If N $_e$ is small, then a bootstrap sample may have fewer than three unique epochs, which are required to solve for the aforementioned five astrometric parameters. We have set the required minimum unique epochs to five to reduce the biasing in the fit due to the chance selection of unique epochs at one side of the parallax ellipse, and all with a small parallax contribution to estimate the underlying parallax swing. If a bootstrap sample has fewer than the minimum unique epochs, then the sample is not considered in the bootstrap fitting. For each valid bootstrap sample set, astrometric parameters are inferred using PMPAR software, and a probability density distribution is examined for each parameter from which the median value of the parameter and the $84^{\rm th}$ and $16^{\rm th}$ percentile values ( $1\sigma$ uncertainty) are estimated. The bootstrap fitting estimates for both pulsars are shown in Figures 3 and 4 and tabulated in Table 5. For comparison, we also include the astrometric estimates for these pulsars from Deller et al. (Reference Deller2019). However, they reported the most probable values of the parameters and used the most compact 68% confidence interval around the parameter value to estimate the $1\sigma$ error bars, which is less robust, particularly for asymmetric distributions. Our bootstrap parameter estimates from the combined data set (BD152+BD174) are consistent with those of Deller et al. (Reference Deller2019) within the reported uncertainties, but with a significant improvement in precision (Table 5). As expected, the enhancement is more pronounced in proper motion than in parallax since the proper motion precision grows with observation time as $t^{3/2}$ , whereas parallax precision grows with $t^{1/2}$ (assuming all observations are of equal sensitivity and equally spaced in time).

Figure 3. Illustration of bootstrap fitting estimates for B0329+54. Top: pulsar position offsets (R.A., left and Dec., right) from the reference epoch (MJD: 56300) position with error bars, plotted as a function of epochs after removing the best-fit proper motion. The first nine data points (from left) are from BD152 observations, while the remaining data points are from BD174 observations. Each BD174 epoch includes two data points, representing the 1.4 and 1.6 GHz datasets. Each light blue line represents the parallax signature fitted to a respective bootstrap sample, and the band of these lines indicates the spread in the fitted parallax signature. The increased spread in the parallax signature curves for the Dec. offset, compared to the R.A. offset, can be attributed to a combination of a) sampling epochs that sample R.A. offsets closer to the parallax signature extremes in earlier epochs (i.e., BD152), and b) relatively larger error bars in Dec. offset estimates. Middle and bottom: probability density functions for the parallax and the proper motion in R.A. and Dec. are shown for both the bootstrap in black colour and least-squares results (Gaussian function with the parameters estimated using PMPAR after adding the systematic error contribution to obtain the reduced chi-square to be one, plotted in blue colour). The dotted vertical lines indicate the most probable value (black) and the median value (red), with green dashed-dotted lines marking the 68% confidence interval ( $\pm 1 \sigma$ ) estimated from the bootstrap fit.

Figure 4. Bootstrap fitting estimates for B1133+16. In the top panel plots, the first eight data points (from left) are from BD152 observations, and the rest of the data points are from BD174 observations. For details, see the caption of Figure 3.

The BD174 data set exhibits comparatively larger error bars (Figures 3 and 4) due to a couple of reasons: (a) higher thermal noise resulting from less than half the target scan time compared to BD152, (b) observations taken during the solar maxima; the typical mean Vertical TEC (VTEC) value during BD174 observations $\sim25$ TECU ( $1\,\textrm{ TECU} = 10^{16}\,\textrm{electrons}\,{\rm m}^{-2}$ ), which is considerably higher than the $\sim17$ TECU seen during BD152 observations (estimated by randomly sampling 5–10 epochs from each observation series), (c) bigger beam size due to missing antennas (sometimes) and low-frequency observation in the case of 1.4 GHz data set.

Figure 5. Error ellipses and marginalised histograms are plotted for $\eta_{\textrm{ EFAC}}$ and five astrometric parameters of B0329+54 using the Bayesian approach. The contours show the $1\sigma$ , $2\sigma$ , and $3\sigma$ confidence intervals, and the position offset ( $\Delta\alpha$ and $\Delta\delta$ ) is with respect to the reference epoch position (MJD: 56300), while dashed lines in histograms mark the median value of astrometric parameters and $1\sigma$ deviations.

Figure 6. Error ellipses and marginalised histograms are plotted for $\eta_{\textrm{ EFAC}}$ and five astrometric parameters of B1133+16 using the Bayesian approach. For details, see the caption of Figure 5.

Figure 7. The left plot shows the difference between the traditional and Petrov23 mapping functions (equivalent to STEC/VTEC), while the right plot illustrates the position offset caused by residual dispersive delay at 1.6 GHz when the delay is computed using the respective mapping function. A nominal TEC value of 10 TECU is assumed, and the IBC is assumed at six arcminutes above the target.

4.3. Bayesian approach

We have also inferred astrometric parameters based on Bayesian statistics using STERNE Footnote ^h package (Ding et al. 2021, Reference Ding2022). This tool is specifically designed for VLBI astrometry and allows estimation of additional parameters, such as orbital parameters in the case of a binary system and $\eta_{\text{EFAC}}$ , a scale factor to account for the underestimation or overestimation of systematic uncertainties. The $\eta_{\text{EFAC}}$ factor is assumed to be one in the least-squares fitting and bootstrap approach. The total error at the $i{\rm th}$ epoch is

(3) \begin{equation} \sigma^i({\eta_{\text{EFAC}}}) = \sqrt{(\sigma_R^i)^2 + (\eta_{\text{EFAC}} \cdot \sigma_{\text{sys}}^i)^2}, \end{equation}

where $\sigma_R$ and $\sigma_{\text{sys}}$ are the random and systematic errors.

The priors for astrometric parameters are assumed uniformly distributed, $\mathcal{U}(X^{\text{DF}} - 20 \widetilde{\sigma}^{\text{DF}}_X, X^{\text{DF}} + 20 \widetilde{\sigma}^{\text{DF}}_X)$ , where $X^{\text{DF}}$ and $\widetilde{\sigma}^{\text{DF}_X}$ stand for astrometric parameters and associated scaled uncertainties ( $\widetilde{\sigma}^{\text{DF}}_X = \sqrt{\chi^2}\sigma^{\text{DF}}_X$ ), taken from the least-squares fitting, and $\mathcal{U}(0, 15)$ for $\eta_{\text{EFAC}}$ . The uniform distribution avoids biasing, as we have used PMPAR results obtained from the same data to set the priors. Figures 5 and 6 illustrate the inferred astrometric parameters for B0329+54 and B1133+16 using the Bayesian approach alongside $\eta_{\textrm{ EFAC}}$ , and the median value of astrometric parameters along with their associated uncertainties are tabulated in Table 5. The $\eta_{\text{EFAC}}$ value of unity and two for B0329+54 and B1133+16, respectively, indicates that the initially estimated systematic uncertainties are underestimated for the latter.

All astrometric parameters obtained using the three mentioned statistical techniques are consistent within the error bars. However, the bootstrap and Bayesian methods perform better, as they rely less on prior knowledge of uncertainties in the input data. In particular, the Bayesian approach proves more robust when reliable prior information is available for the fitting parameters. Therefore, we will only consider Bayesian estimates for the remainder of the discussion.

For each pulsar, the proper motion and parallax estimates obtained from 1.4 and 1.6 GHz data sets are consistent within their error bars, except for the parallax of B1133+16, which deviates at the $1.5 \sigma$ level (Table 7). As expected, the 1.6 GHz data set provides better astrometric precision, and the combination of 1.4 and 1.6 GHz further enhances the precision of astrometric parameters by approximately 1.2 times. Similarly, incorporating the BD174 data set alongside BD152 improves astrometric precision by increasing the number of observations and, more importantly, extending the observation time span. Comparing the Bayesian estimates of the inferred parameters from the BD152-only data set (Table 5, reference epoch MJD 56000) with those from the combined BD152+BD174 data set (Table 7, reference epoch MJD 56300), the proper motion and parallax estimates for both pulsars remain consistent within the respective error bars. Both analyses used the traditional mapping function. However, the combined data set yields almost threefold improvement in proper motion precision, and parallax precision improves by nearly a factor of two.

Table 6. Median of position-time series differences ( $\Delta\alpha$ , $\Delta\delta$ ) obtained using each mapping function (Figures 8 and 9). The second and third columns correspond to results from the BD174 data set.

Table 7. Bayesian estimates of astrometric parameters for B0329+54 and B1133+16 utilising each TEC mapping function to compute the ionospheric dispersive delays and different combinations of data sets included in the fitting.

¹ The position offsets ( $\Delta\alpha_{\text{J2000}}$ , $\Delta\delta_{\text{J2000}}$ ) are with respect to the mean reference position of B0329+54 and B1133+16, which are, R.A. = 03:32:59.4109883, Dec. = 54:34:43.31937 and R.A. = 11:36:03.1154509, Dec. = 15:51:14.48331, respectively. The reference positions are with respect to the IBC at a reference epoch, MJD: 56300, and associated error bars do not include the position uncertainties in the PRC and IBC and other systematic errors.

5. Effect of TEC mapping function on astrometric parameters

The TECOR task performs a priori ionospheric delay corrections by computing the excess path delay using global TEC maps. For both data sets, the ionospheric dispersive delays are estimated using TEC maps provided by the International Global Navigation Satellite System (GNSS) Service (centre code: IGSG). These maps represent the VTEC distribution on a global grid with a spatial and temporal resolution of $5^{\circ}\, \times\, 2.5^\circ$ ( $\textrm{ longitude}\, \times\, \textrm{ latitude}$ ) and 2 h, respectively. Co-located GNSS receivers at VLBI stations have also been characterised (Skeens et al. Reference Skeens, York, Petrov, Munton, Herrity, Ji-Cathriner, Bettadpur and Gaussiran2023), and may be useful in future for deriving site-specific TEC measurements and enabling real-time ionospheric corrections.

The TECOR task utilises a mapping function to convert the VTEC to Slant TEC (STEC). The traditional mapping function is the so-called single-layer model, $M(e)\, =\, 1/\cos(\textrm{ ZA})$ , where ZA is the Zenith angle at the ionosphere puncture point. Recently, Petrov (Reference Petrov2023) has studied a so-called thin-shell ionospheric mapping function, suggested by Schaer (Reference Schaer1999),

(4) \begin{equation}M(e) = k \dfrac{1}{\sqrt{1 - \left(\dfrac{R_{\oplus}}{R_{\oplus} + H_i + \Delta H}\right)^2 \cos^2\alpha e_{gc}}},\end{equation}

here, $R_{\oplus}$ denotes the Earth’s radius, $e_{gc}$ is the geocentric elevation angle with respect to the radius vector between the geocentre and the station, and $H_i$ is the height of the ionosphere. Petrov (Reference Petrov2023) has recommended values for scaling factor (k), thickness of the ionospheric layer ( $\Delta H$ ), and fudge factor ( $\alpha$ ) to 0.85, 56.7 km, and 0.9782, respectively. The modified (Petrov23) mapping function aims to enhance the accuracy of ionospheric corrections based on global TEC maps, and has been implemented in the TECOR task of AIPS since 31DEC23. We emphasise that these global ionospheric maps only capture the bulk/smooth component of the ionosphere, due to their limited time and spatial resolution, and small-scale or rapidly varying ionospheric disturbances will remain uncorrected regardless of the TEC mapping function used.These two mapping functions exhibit an offset in STEC prediction that becomes increasingly significant at lower elevations. Even at an elevation of $90^{\circ}$ , the Petrov23 mapping function predicts a dispersive delay that is approximately 20% lower than that of the traditional mapping function for a given VTEC value (Figure 7). These differences affect the inferred source position and can introduce noticeable position shifts even across the small angular separations typical for differential astrometry. For example, consider a field in which the target is separated from a calibrator source by six arcminutes in declination, observed at an elevation of $35^{\circ}$ , at a frequency of 1.6 GHz, and with a maximum baseline length of 8000 km. When ionospheric delays are corrected using TEC maps (assuming a nominal TEC value of 10 TECU) utilising each of the mapping functions, the resulting position shift can be as large as $\sim0.1$ mas (Figure 7).

To examine the impact of TEC mapping functions on astrometry, we have processed the data twice, using each mapping function and compared the resulting position-time series. For each pulsar, the median of the difference of position-time series obtained using traditional and Petrov23 mapping functions is shown in Table 6. The quantity determining the shift in the target position, assuming the shift is only due to the switching between mapping functions, is the difference in dispersive delays (computed by TECOR) between the IBC and target directions when using the traditional and Petrov23 mapping functions. The obtained position-time series difference is within the expectation from the position difference shown in Figure 7 at higher elevations.

The astrometric parameters obtained using each mapping function are shown in Table 7. The difference in the fitted astrometric parameters is within the respective error bars. However, it is not apparent that the Petrov23 mapping function is definitely an improvement, by examining the precision of the astrometric parameters and the reduced chi-square values (in least-squares fitting) for different combinations of data sets. Although the reduced chi-square values are not appreciably better or worse for the full BD152+BD174 data in either case, nor are the BD174-only results consistently closer to the truth (where we assume that BD152+BD174 gives a reasonable approximation of truth, at least compared to BD174 alone) with one mapping function compared to the other. We conclude that significant model uncertainty remains in the application of ionospheric corrections based on global ionosphere models. In any case, ionospheric variations on small spatial and temporal scales that are not captured by global TEC maps (regardless of the mapping function used) may be more significant for differential astrometric observations than the difference between the mapping function used in global TEC maps, at least at small angular separations typical for in-beam calibrator observations.

6. Impact of frequency variations on astrometric parameters

To probe the ionospheric propagation effects on astrometry, we have analysed 1.4 and 1.6 GHz data sets independently. For each source, the fitted reference epoch positions have an offset at the 1.4 and 1.6 GHz (Table 7), which remain almost the same when switching between different TEC mapping functions. It seems this position offset might be due to the mis-modelled structure of the IBC and the residual ionospheric error, which are associated with the ionospheric fluctuations within the spatial and temporal resolution of the global TEC maps.

For each pulsar, the ratio of the uncertainties to the parallax at both bands is close to the ratio of the wavelength squared (almost 1.5; Table 7), which one can expect if the ionosphere is the dominant contributor to the position uncertainties. This underscores the notion that the dominant systematic error in position arises from ionospheric effects at low frequencies. A comparison of astrometric parameters at 1.4 and 1.6 GHz is presented in Table 7, confirming the anticipated superiority of astrometric results at 1.6 GHz, irrespective of the TEC mapping function.