© The Authors 2024

1 Introduction

Millisecond pulsars (MSPs) refer to recycled fast-spinning neutron stars, which exhibit unparalleled spin stability compared to other pulsars (Hobbs et al. 2010). Using the pulsar timing technique (e.g., Detweiler 1979), which records and models the pulse arrival times, astronomers have used MSPs to deliver the most stringent tests of gravitational theories (e.g., Freire & Wex 2024). Collectively, an array of MSPs scattered across the sky, known as a pulsar timing array (PTA), can be used to directly probe the stochastic gravitational-wave background (GWB) in the nanohertz (nHz) regime (Sazhin 1978).

Recent years have seen the major PTA consortia get closer than ever before to achieving high-significance detections of a homogeneous GWB (Agazie et al. 2023; Reardon et al. 2023; Antoniadis et al. 2023; Xu et al. 2023). Despite the breakthrough, deepening our understanding of the sources of the GWB still requires continuous improvement of the PTA sensitivities. The optimal strategy to sustain PTA sensitivity enhancement is to regularly add new MSPs to the PTAs (Siemens et al. 2013), as carried out by the MeerKAT Pulsar Timing Array (MPTA; Miles et al. 2023). However, it is generally difficult to quantify the red timing “noises” (in which the GWB signal resides) for a relatively shortly timed (≲3000 days) MSP; one way to overcome this difficulty is to incorporate independent astrometric measurements (i.e., sky position, proper motion, and parallax) into the inference of timing parameters (Madison et al. 2013).

The very long baseline interferometry (VLBI) technique can provide precise, robust, and model-independent astrometric measurements for MSPs, and can achieve a certain astrometric precision in much shorter time span compared to using pulsar timing (e.g., Brisken et al. 2002; Chatterjee et al. 2009; Deller et al. 2019). Therefore, the incorporation of precise VLBI astrometric measurements into the timing analysis of MSPs is essential for testing gravitational theories (e.g., Deller et al. 2008; Kramer et al. 2021), and may substantially enhance the PTA sensitivities (Madison et al. 2013).

However, this incorporation is technically challenging. First, incorporating precise VLBI proper motion and parallax into a timing analysis can be limited by the potential temporal structure evolution of the reference sources used in VLBI astrometry. Second, incorporating a VLBI pulsar position into a timing analysis necessitates a good understanding of the transformation between the two distinct kinds of reference systems used by VLBI astrometry and pulsar timing. Though both VLBI astrometry and pulsar timing are usually presented with respect to the barycenter of the Solar System, VLBI astrometry is conducted in the kinematic reference frame established with remote AGNs that are quasi-static on the sky (making them robust to inaccurate planetary ephemerides), while pulsar timing studies are carried out in the dynamic reference frame, which requires reliable planetary ephemerides to convert Earth-based pulse arrival times to the barycenter of the Solar System. The dynamic coordinate is floated to a kinematic coordinate system through observations of common objects, such as differential observations of asteroids with respect to stars. VLBI observations of MSPs with respect to AGNs will allow us to determine a rotation of the dynamic coordinate system defined by planetary ephemerides with respect to the inertial coordinate system (based on VLBI observations of AGNs) (Madison et al. 2013; Wang et al. 2017; Liu et al. 2023a).

The residuals in pulsar positions from VLBI and timing observations after a subtraction of the rotation will allow us to provide an independent assessment of pulsar timing errors and validate the PTA error model. We consider this question as a matter of great importance because any claim that a PTA has detected GWB is based upon a model of pulsar timing errors.

So far, planetary-ephemeris-dependent frame rotation remains poorly constrained. This is mainly because of the relatively large (≳1 mas) VLBI position uncertainties of MSPs (e.g., Liu et al. 2023a) as compared to the ≲0.2 mas timing position uncertainties for the best-timed MSPs (e.g., Perera et al. 2019). According to Liu et al. (2023b), 50 VLBI-based MSP positions with the current precision level would constrain the frame rotation to the 0.3 mas level; this constraint would be improved to the 0.1 mas level if VLBI could determine MSP positions to 0.3 mas precision. Therefore, reducing the VLBI position uncertainties of MSPs is key to building the planetary-ephemeris-dependent frame tie, examining the quality of planetary ephemerides, and hence facilitating the quantification of the timing noises resulting from inaccurate planetary ephemerides.

In this paper, we introduce and test a novel method that allows us to significantly improve the precision of pulsar VLBI positions. Throughout the paper, uncertainties are provided at 68% confidence, unless otherwise stated; the mathematical expressions (including the subscripts and superscripts) defined throughout this paper are universally valid.

2 A novel observing strategy and test observations

2.1 The PINPT strategy

As radio pulsars are generally faint (∼1 mJy at 1.4 GHz, and with steep spectra), the standard approach for directly determining positions in the quasi-inertial VLBI frame involving the measurement of group delays is not feasible. Instead, the absolute positions of pulsars have been determined using differential astrometry with respect to relatively bright reference sources, which are normally AGNs that are not point-like. As the position of a suitable nearby AGN can be determined using standard absolute astrometry techniques, such relative position measurements allow a connection between the pulsar position and the quasi-inertial frame. However, standard absolute astrometry techniques do not account for any structure in the AGN, and so the position extracted for any one of these sources is effectively that of the peak brightness in the image. The brightest spot in the 2D brightness distribution (or simply image) of an AGN also serves as the reference point for differential astrometry, which is usually the optically thick jet core (as long as the AGN is not flaring). As the AGN images normally vary with observing frequency ν, the reference point (or the jet core) evolves with ν as well; additionally, frequency-dependent image models of the reference sources are required for pulsar astrometry.

Three main error sources contribute to the error budget of the absolute pulsar position (see Sect. 3.2 of Ding et al. 2020a), namely (i) the uncertainty in the absolute position of the primary phase calibrator (derived and registered in the Radio Fundamental Catalogue¹, or the ICRF3 catalog; Charlot et al. 2020), (ii) the unknown frequency-dependent core shifts (hereafter simply referred to as core shifts; e.g., Bartel et al. 1986; Lobanov 1998) of the reference sources (or more generally, the frequency evolution of reference source structures), and (iii) differences in the line-of-sight propagation delay between the direction of the calibrator source (where it has been solved) and the direction of the target. At L band, core shifts amounting to ∼1.2 mas (Sokolovsky et al. 2011) usually dominate the error budget; second to that, the propagation-related systematic error of the absolute pulsar position is also prominent, given the relatively large (≳1 deg) separation between the pulsar and its primary phase calibrator.

To suppress the aforementioned positional uncertainties, we designed a special observing strategy of differential astrometry, which extends the core-shift-determining method pioneered by Voitsik et al. (2018), and combines the method with the multiview (referring specifically to 2D interpolation throughout this paper) strategy (e.g., Rioja et al. 2017; Hyland et al. 2023). The proposed observing strategy, referred to as the PINPT (Phasescreen Interpolation plus frequeNcy-dePendent core shifT correction; read as “pinpoint”) strategy, requires a group of less than or approximately equal to 6 observations for absolute position determination: ≲3 L-band Multi-View observations of the pulsar (pulsar sessions), and ≲3 core-shift-determining observations on nearby AGNs (which must include the three calibrators used in the L-band Multi-View sessions), each at different observing frequency ν. Where suitable L-band in-beam calibrators are identified, the pulsar sessions can also be used as the L- band core-shift-determining session, hence reducing the required observing time.

Our observing strategy of multifrequency observations with multiple phase calibrators has three advantages. First, with the multifrequency observations of AGNs, the core shifts of the AGNs will be well determined, which will significantly reduce the core-shift-related errors of the absolute pulsar position. Second, using the Multi-View strategy will remove propagation- related systematic errors to at least the first order (e.g., Ding et al. 2020b). Finally, with the three phase calibrators of the MultiView setup, the pulsar position uncertainties due to uncertain reference source positions will drop by a factor of ≤ √3 (see Sect. 4.4 of Ding et al. 2020b) compared to using only one phase calibrator.

2.2 Observations

To test the PINPT strategy, we ran four observing sessions in December 2021 and January 2022 using the Very Long Baseline Array (VLBA) on PSR J2222-0137, a millisecond pulsar that has been well determined astrometrically by previous campaigns (Deller et al. 2013; Guo et al. 2021, hereafter referred to as D13 and G21). The observations, carried out under the project code BD244, include two 2 h L-band Multi-View observations and two 2hr core-shift-determining observations at S/X band and Ku (∼15 GHz) band, respectively. The first L-band observation was scheduled within 2 days of the S/X- and Ku-band observations (in order to minimize structure evolution of the phase calibrators), and served as both pulsar session and core-shift-determining session; whereas the second L-band observation was solely a pulsar session.

More specifically, in the L-band observations, a “target pointing” covers PSR J2222-0137 and two in-beam (at L band) calibrators identified by D13, namely FIRST J222112-012806 (J2221) and FIRST J222201-013236 (J2222); scans on this field were interleaved with scans on three brighter but off-beam AGNs (hereafter referred to as off-beam calibrators), namely ICRF J221852.0-033536 (J2218), ICRF J221947.2-005132 (J2219), and ICRF J222646.5+005211 (J2226), with a cycle time of 5 minutes. In the S/X- and Ku-band observations, J2221 and the three off-beam calibrators (hereafter referred to as the core-shift-probing AGNs) were observed alternately (with 5min and 2min cycles, respectively). All of the four core-shiftprobing AGNs are selected to have displayed resolved jet-core radio features in VLBI data, which eases the core-shift determination (see Sect. 4.2.2). Being fainter at higher observing frequencies, PSR J2222-0137 was skipped in the S/X- and Ku- band observations. The calibrator plans of the observations are displayed in Fig. A.1. In addition to the phase calibrators relatively close to PSR J2222-0137 on the sky, the bright blazar ICRF J214805.4+065738, further away from the sky region of interest, was used as the fringe finder to correct instrumental delays and filter bandpass. The unresolved flux densities and the angular separations of the aforementioned sources can be found in Table A.1.

To precisely determine the core shifts, a wide variety of observing frequencies are required. The S/X-band observation simultaneously covers two frequency bands, at around ∼2.3 and ∼8.4 GHz. Additionally, the L-band observations were designed to cover two separate frequency ranges centered around 1.44 and 1.76 GHz, respectively. Altogether, we sampled five distinct frequency ranges for the purpose of core-shift determination. The observing configurations, the phase calibration strategies (see Sect. 3), and the purposes of the four BD244 sessions are summarized in Table 1. The VLBA data prior to data reduction can be accessed with the project code BD244 at the NRAO Archive Interface².

3 Data reduction and direct results

All data reduction was performed with the psrvlbireduce³ pipeline, which runs functions of AIPS (Greisen 2003) through ParselTongue (Kettenis et al. 2006), and images sources with DIFMAP (Shepherd et al. 1994). The data reduction follows the standard workflow described in Ding et al. (2020a), except for the phase calibration. Despite the total observing time of only 8 hours, the phase calibration of this work is sophisticated. There are two different procedures of phase calibration, and which is used depends on the purpose of the session, and more precisely on whether it is a target session or a core-shift-determining session. As noted in Sect. 2.2, the first L-band observation serves as both a target session and a core-shift-determining session, and is therefore reduced twice in two distinct procedures. The phase calibration strategies applied to the four observing sessions are summarized in Table 1 and described as follows.

3.1 Target sessions

For the data reduction of target sessions, we did not split the data into two by observing frequency (one at around 1.44 GHz and the other at around 1.76 GHz). By doing so, the average central frequency 1.6 GHz agrees with previous astrometric campaigns of PSR J2222-0137 (D13, G21), and the acquired image S/N (and hence the positional precision) is not lowered. As mentioned in Sect. 2.2, the two L-band observations involve five phase calibrators, including off-beam calibrators and two in-beam calibrators (also see Fig. A.1 and Table A.1). We implemented the phase calibration of PSR J2222-0137 in two different ways, depending on the astrometric goal.

3.1.1 In-beam astrometry

A previous astrometric campaign targeting PSR J2222-0137 was carried out between October 2010 and June 2012, spanning 1.7 yr (D13, G21). The two new target sessions extend the astrometric time baseline by a factor of 6.6 to 11.3 yr, promising higher astrometric precision (especially for proper motion). To capitalize on the long time baseline, we phase-referenced PSR J2222-0137 to the same reference source (i.e., J2222) used in the previous campaign following the data reduction procedure of D13 (including using the same image model of J2222). The updated results of in-beam astrometry are reported in Sect. 4.1.

3.1.2 Multi-View (2D interpolation)

As one form of the Multi-View strategy, 2D interpolation uses ≥3 reference sources to derive the phase solution at the sky position of the target (e.g., Rioja et al. 2017), which can at least remove propagation-related systematic errors to first order. In order to determine the precise absolute position of PSR J2222-0137, we applied Multi-View (as part of the PINPT strategy) with the off-beam calibrators, all of which have well-determined absolute positions¹. We reiterate that realizing the PINPT strategy requires a combination of Multi-View session(s) and core-shiftdetermining session(s); these two kinds of sessions need to be arranged close to each other to minimize the effects of structure evolution. The second L-band observation is ≈40 days apart from the core-shift-determining sessions. Therefore, MultiView was applied to the first L-band observation, but not the second one.

In practice, we realized the Multi-View in two different approaches, as described in Appendix B. Both approaches consistently render the position 22^h22^m05^s99997 ± 0.1 mas, –01°37′15.″7825±0.2 mas for PSR J2222-0137 at MJD 59554.0, and 22^h22^m01^s.373131 ±0.03 mas, –01°32′36.″797654 ±0.06mas for J2222. For further examination, we also applied the two approaches of Multi-View to the X-band data (where the chance of phase wrap errors is much smaller than at L band), and confirmed with J2221 that the two approaches give almost identical positions. The availability of the J2222 position, as well as the J2221 position 22^h21^m12^s680887 ± 0.01 mas, –01°28′06.″730985 ± 0.02 mas, offers a distinct pathway to the absolute position of PSR J2222-0137 (see Sect. 5). We note that the three positions presented above can only be considered relative positions with respect to the off-beam calibrators; and the positional uncertainties only include the statistical component related to random (thermal) noise in the VLBI image of PSR J2222-0137.

3.2 Core-shift-determining sessions

The determination of AGN core shifts requires a wide coverage of observing frequencies. As mentioned in Sect. 2.2, the core-shift-determining sessions cover five frequency ranges. We first split the S/X-band data into S-band data and X-band data, and likewise split the wide-L-band data into two datasets – one at ~l.44GHz and the other at ~1.76 GHz. Thus, we acquired five datasets altogether, each taken at a different observing wavelength: λ – 2, 4, 13, 17, and 21 cm.

From each of the five datasets, we measured the J2221 position with respect to each off-beam calibrator by (1) phasereferencing J2221 to each off-beam calibrator, and (2) dividing the phase-referenced J2221 data by its final image model obtained at λ, and (3) fitting the J2221 position. In phase referencing, the final image model of an off-beam calibrator at λ was applied to the phase calibration of the off-beam calibrator; subsequently, the acquired phase solution was applied to both the off-beam calibrator and J2221. To make the final image models of the off-beam calibrators and J2221, we first applied 13 cm models to the datasets of respective sources at all λ during phase calibration (or self-calibration for J2221). In this way, we aligned the reference points of each core-shift-probing AGN across λ. Thereafter, we singled out the aligned datasets, and remade the models of J2221 and off-beam calibrators at each λ. The final image models of in-beam and off-beam calibrators at all λ have been made publicly available⁴ to allow the reproduction of our results.

4 Astrometric parameters and core shifts

The data reduction described in Sect. 3 produced (a) two pulsar positions measured with respect to J2222, (b) Multi-View positions of the pulsar, J2221, and J2222, and (c) 5 × 3 = 15 J2221 positions from the core-shift-determining sessions. The products (a) and (c) can provide stringent constraints on astrometric parameters and core shifts, respectively.

4.1 Astrometric inference

We added the two new pulsar positions measured with respect to J2222 to previous pulsar positions (measured with respect to J2222), and remade the Bayesian inference as described in G21. The resultant astrometric parameters are provided in Table 2. For comparison, the previous VLBI results reported by G21 are also listed in Table 2. Unsurprisingly, the proper motion improves substantially, namely by a factor of approximately 7, while the parallax is also enhanced by ≈14%, corresponding to a refined trigonometric distance of ${268.6}_{-0.9}^{+1.0}$ pc for PSR J2222-0137. On the other hand, the 3σ discrepancy in µ_δ reveals an underestimated systematic uncertainty in Declination reported by previous works, which is likely due to the vertical J2222 structure evolution that has a much greater impact on astrometric results determined over short (≲2 yr) time span.

4.2 Core-shift determination

Following the pioneering work of Voitsik et al. (2018), we developed new packages to infer the core shifts of the four coreshift-probing AGNs from the 15 J2221 positions (with respect to three off-beam calibrators at five λ) in a Bayesian style. The core-shift inference requires three ingredients – (1) a prescription of the systematic errors of the 15 J2221 positions, (2) a mathematical description of the core shifts, and (3) an underlying mathematical relation between core shift and λ (or observing frequency ν). These three ingredients are addressed below.

4.2.1 Systematic errors on measured J2221 positions

Along with the 15 J2221 positions, their random errors due to noise in the J2221 images were also obtained with data reduction (described in Sect. 3.2). Additionally, atmospheric propagation effects will introduce systematic errors, which change with ν and the angular separation between J2221 and the reference source.

We quantified the systematic errors by $$\frac{{\sigma }_{ijk}^{\mathcal{S}}}{1\text{\hspace{0.17em}mas}}={\eta }_{\text{EFAC}}\cdot {\eta }_0\cdot \frac{s_i}{1\text{\hspace{0.17em}deg}}\cdot \widehat{l}(ϵ)\cdot \frac{{\sigma }_0\left(v_j\right)}{1\text{\hspace{0.17em}mas}}\cdot {\widehat{\Theta }}_k,$$(1)

where ${\sigma }_{ijk}^{\mathcal{S}}$ denotes systematic errors; the subscript i = J2218, J2219, J2226 specifies the reference source; the second subscript j = 1–5 corresponds to one of the five observing frequencies (see Sect. 3.2); the last subscript k = α, δ refers to Right Ascension (RA) or Declination (Dec); η_EFAC is a scaling factor (for the estimated systematic uncertainty) to be determined with Bayesian analysis; η₀ is the initial scaling factor that brings ${\sigma }_{ijk}^{\mathcal{S}}$ to a reasonable value (∼1 mas) before inferring η_EFAC, which eases the inference of η_EFAC at a later stage; s_i represents the angular separation between J2221 and the reference source; ${\widehat{\text{Θ}}}_k$ denotes the fractional synthesized beam size (projected to RA or Dec); and $\widehat{l}$ and σ₀ stand for the two terms changing with antenna elevations ϵ and ν, respectively. In this work, we set η₀ ≡ 1, and the adopted $\widehat{l}(ϵ)$ and σ₀(ν) are described in Appendix C.

4.2.2 The directions of core shifts

Regarding the ingredient (2), we describe the core shift of the i′- th (i′ = J2218, J2219, J2221, J2226) AGN with two parameters – the core-shift magnitude $r_{i^{\prime }}$, and the direction of core shift ${\theta }_{i^{\prime }}$. As noted in Voitsik et al. (2018), it is difficult to infer both $r_{i^{\prime }}$ and ${\theta }_{i^{\prime }}$ for each of the four core-shift-probing AGNs from the 15 J2221 positions without prior constraints, which has been taken into account in the design of the observations. As mentioned in Sect. 2.2, all of the core-shift-probing AGNs are selected to have displayed a clear jet-core radio feature on VLBI scales. Following Voitsik et al. (2018), we assumed that core-shift directions are aligned with the respective AGN jet directions, thereby gaining prior knowledge of the core-shift directions with analysis of VLBI images of the AGNs. The determination of AGN jet directions is detailed in Appendix D.

4.2.3 The relation between core shift and observing frequency

The relation between core shift and observing frequency is given by Konigl (1981) as $r\propto v^{-1/k_r},$, where the index k_r is equal to 1 when synchrotron self-absorption (that leads to the jet core radio emissions) is in equipartition (Blandford & Königl 1979). For the i′ -th AGN, we initially adopted an equivalent formalism, $$r_{i^{\prime }}=r_{0i^{\prime }}{\left(\frac{v}{v_0}\right)}^{-{\beta }_{i^{\prime }}},$$(2)

as the relation between $r_{i^{\prime }}$ and ν, where ${\beta }_{i^{\prime }}$ is a power index to be determined in Bayesian analysis; and v₀ and $r_{0i^{\prime }}$ refer to the reference frequency and the core-shift magnitude at the reference frequency, respectively.

Using the Bayesian inference described in Appendix E, we derive ${\beta }_{i^{\prime }}$ along with other model parameters, which are provided in Table E.1. For the derivation, we choose v₀ = 2.27 GHz, the median of the five central frequencies (see Sect. 3.2 and Fig. 1). The $r_{0i^{\prime }}$ at a different v0 can be calculated using Eq. (2). When ${\beta }_{i^{\prime }}$ are included in the inference, the reduced chi-square ${\chi }_v^2$ of inference is 1.9. Although ${\beta }_{i^{\prime }}$ is significantly (≳3 σ) determined for 2 of the 4 AGNs, it is consistent with 1 in all cases. In comparison, when performing the inference with all ${\beta }_{i^{\prime }}$ fixed to 1, we acquired consistent results with generally higher precision. Moreover, the ${\chi }_v^2$ decreases to 1.2, suggesting ${\beta }_{i^{\prime }}\equiv 1$ as an appropriate assumption for this work. We adopt the results derived assuming ${\beta }_{i^{\prime }}\equiv 1$ for further analysis. The adopted core-shift model is illustrated in Fig. 1.

Fig. 1 Jet core relative positions of FIRST J222112-012806 measured with respect to the jet cores of three reference sources (i.e., ICRF J221852.0-033536, ICRF J221947.2-005132, and ICRF J222646.5+005211) at five observing wavelengths. Systematic uncertainties calculated with Eq. (1) have been included in the positional uncertainties. Models of the combined core shifts are derived assuming ${\beta }_{i^{\prime }}\equiv 1$ with Bayesian inference (see Appendix E), and displayed with straight lines.

5 An ultraprecise MSP position obtained with VLBI

With $r_{0i^{\prime }}\text{ and }{\theta }_{i^{\prime }}$ determined for the three off-beam calibrators and J2221, we calculated the core shifts of the four AGNs at 1.6 GHz, and proceeded to derive the absolute position of PSR J2222-0137. The derivation was realized with three different approaches. In the first approach, the pulsar position obtained with Multi-View (presented in Sect. 3.1.2) was corrected for the core-shift refinements in the assumed positions of the off- beam calibrators (as described by Eq. (F.3) in Appendix F). We obtained the pulsar position ${\alpha }_1^{\mathcal{G}}={22}^{\text{h}}{22}^{\text{m}}{05}^{\text{s}}.99993\pm 0.17\text{\hspace{0.17em}mas},{\delta }_1^{\mathcal{G}}=-{01}^{\circ }{37}^{\prime }{15}^{\prime \prime }7841\pm 0.32\text{\hspace{0.17em}mas}$ in the geocentric reference frame at the reference epoch t_ref of MJD 59554.0, which corresponds to ${\alpha }_1^{ℬ}={22}^{\text{h}}{22}^{\text{m}}06.00015\pm 0.17\text{ mas, }{\delta }_1^{ℬ}=-{01}^{\circ }{37}^{\prime }{15}^{\prime }\cdot 7827\pm 0.32\text{\hspace{0.17em}mas}$ with respect to the barycenter of the Solar System after removing the parallax effects. Here, the positional uncertainty is the addition in quadrature of the statistical uncertainty (provided in Sect. 3.1.2) and the ∆x_target uncertainty described in Appendix F.

In the second approach, we first acquired the J2221 position x_J2221 at 1.6 GHz by correcting the X-band J2221 position (obtained with Multi-View, see Sect. 3.1.2) using the same process employed in the previous approach for the target source PSR J2222-0137. Subsequently, we reprocessed all VLBA data (including the two new pulsar sessions and the historical ones) of PSR J2222-0137, phase-referencing PSR J2222-0137 to J2221. With the position series obtained from the data reduction, we inferred the reference pulsar position $x_{\text{psr}}^{\text{J}2221}$ (measured with respect to J2221) and its uncertainty at t_ref = MJD 59554.0 (along with other astrometric parameters) using the Bayesian analysis described in Ding et al. (2023). Combining x_J2221, the image model position $x_{\text{J}2221}^{\text{model }}$ of J2221, and $x_{\text{psr}}^{\text{J}2221}$, the absolute position of PSR J2222-0137 was determined with Eq. (F.6) to be ${\alpha }_2^{ℬ}={22}^{\text{h}}{22}^{\text{m}}{066}^{\text{s}}.00015\pm 0.19\text{\hspace{0.17em}mas},{\delta }_2^{ℬ}=-{01}^{°}{37}^{\prime }{15}^{\prime \prime }\cdot 7825\pm 0.42\text{\hspace{0.17em}mas}$, where the error budget of the position is detailed in Appendix F.1 and Table F.2.

The last approach is essentially the same as the second one, except that J2222 is used instead of J2221. Accordingly, the absolute pulsar position was calculated with the J2222 position x_J2222 at 1.6 GHz, the image model position $x_{\text{J}2222}^{\text{model }}$ of J2222, and the reference pulsar position $x_{\text{psr}}^{\text{J}2222}$ measured with respect to J2222. We obtained ${\alpha }_3^{ℬ}={22}^{\text{h}}{22}^{\text{m}}06.00015\pm 0.14\text{\hspace{0.17em}mas},{\delta }_3^{ℬ}=-{01}^{\circ }{37}^{\prime }{15}^{\prime \prime }.7828\pm 0.27\text{ mas}.$.

The pulsar positions derived with the three approaches are summarized in the upper part of Table 3. The intermediate results leading to the three pulsar positions are provided in Table F.2. Though all three approaches render highly consistent absolute pulsar positions, the first approach offers a snapshot (on the timescale of AGN structure evolution) localization (of PSR J2222-0137) that is independent of the previous astrometric observations of PSR J2222-0137, as opposed to the two other approaches (which use the multifrequency multi-source observations to perfect the position and structure of one of the in-beam calibrators – rather than the pulsar – at that snapshot in time, before standard differential astrometry is performed using the frozen model of the astrometric reference source). Therefore, we report $\left({\alpha }_1^{ℬ},{\delta }_1^{ℬ}\right)$ at tref = MJD 59554.0 as the primary absolute position of PSR J2222–0137.

6 Comparison to timing positions

A priori, we expect the position difference between the VLBI PINPT measurement and that made by pulsar timing to be consistent. With a sample size of 1, it is difficult to ascribe any discrepancy to one or the other of the measurement techniques, or to a systematic difference between dynamic and kinematic frames – but we can use the level of agreement to set a probabilistic upper limit on the error contribution from any of these three potential sources. The most precise published timing-based position of PSR J2222-0137 is reported in G21. In Table 3 of G21, three positions are provided for t_ref = MJD 55743, which include one position derived assuming a proper motion determined with VLBI, and two positions derived without using any VLBI prior (hereafter referred to as timing-only positions).

To test the PINPT strategy with pulsar timing, we re-derived $\left({\alpha }_2^{ℬ},{\delta }_2^{ℬ}\right)$ and $\left({\alpha }_3^{ℬ},{\delta }_3^{ℬ}\right)$ at tref = MJD 55743 (same as that of G21) following the method described above in Sect. 5, under the assumption that the structures of J2221 and J2222 do not evolve with time. The results are displayed in Table 3. Intermediate results for the calculations are provided in Table F.2. Compared to the absolute pulsar position determined for MJD 59554 using either the J2221 or the J2222 approach, the pulsar position derived for MJD 55743 with the same approach is reported with smaller uncertainty. This is because the reference pulsar position $x_{\text{psr}}^{\text{PR}}$ is better constrained at MJD 55743 (see Appendix F.1 and Table F.2), around which the pulsar was observed more densely with VLBI.

To minimize the correlation between the timing and VLBI positions, we only compare the VLBI positions to the two timing-only positions (of G21), which are listed in Table 3 as the “DMX” and “non-DMX” positions. Here, “DMX”, named by G21, refers to the pulsar timing model that describes DM with a piecewise constant function (Demorest et al. 2013; G21), while “non-DMX” refers to the timing model that approximates DM variations with a cubic function (G21). The DMX and non-DMX timing positions are marginally consistent to within 2σ (here and hereafter, σ refers to the addition in quadrature of the uncertainties of the two compared sides), with the uncertainty of the DMX position being more conservative than the non-DMX one.

On the other side of the comparison, ${\delta }_2^{ℬ}$ and ${\delta }_3^{ℬ}$ are consistent with each other. However, ${\alpha }_2^{ℬ}$ is smaller than ${\alpha }_3^{ℬ}$ at 3 σ significance, which is associated with a ~3 σ discrepancy between the two µ_α measured with respect to J2221 ($44{.755}_{-0\text{.014 }}^{\text{+0}.017}\text{mas yr}-1$) and J2222 (44.707 ± 0.005 mas yr^–1, as reported in Table 2), respectively. The discrepancy between ${\alpha }_2^{ℬ}$ and ${\alpha }_3^{ℬ}$ likely indicates the violation of the assumption that the structures of J2221 and J2222 do not evolve with time. The jet direction of J2221 is almost aligned with the RA direction (see Fig. D.1). As noted in D13, the relatively severe structure evolution of J2221 is believed to result in (1) a less precise parallax derived with respect to J2221 (as the parallax magnitude in the RA direction is ≈2.4 times larger than in the Dec. direction), and (2) a biased µ_α determination. Specifically, the discrepancy between ${\alpha }_2^{ℬ}$ and ${\alpha }_3^{ℬ}$ (or the µ_α discrepancy) can be explained by a 0.56 mas larger (or fractionally 18% higher according to Table E.1) r_J2221 (1.6 GHz) at MJD 55743 compared to at MJD 59554. On the other hand, thanks to the horizontal jet of J2221, ${\delta }_2^{ℬ}$ is almost free from the impact of structure evolution (of J2221), and is therefore more favorable than ${\delta }_3^{ℬ}$. The favorability is further strengthened by the indication of vertical structure evolution of J2222 mentioned in Sect. 4.1. Therefore, we adopt $\left({\alpha }_3^{ℬ},{\delta }_2^{ℬ}\right)$ as the absolute pulsar position at MJD 55743. A more sophisticated treatment in the future could incorporate positions in both axes from both sources, weighting them by some measure of expected reliability and taking into account the covariance between the two measurements.

When comparing $\left({\alpha }_3^{ℬ},{\delta }_2^{ℬ}\right)$ to the two timing-only positions, we first find good consistency between $\left({\alpha }_3^{ℬ},{\delta }_2^{ℬ}\right)$ and the DMX position, with ${\alpha }_3^{ℬ}$ and the DMX RA being consistent to within 0.3 σ, and ${\alpha }_2^{ℬ}$ and the DMX Dec consistent to within 1.1σ. On the other hand, despite the 1.5σ consistency between ${\delta }_2^{ℬ}$ and the non-DMX Dec., ${\alpha }_3^{ℬ}$ is smaller than the non-DMX RA at >3σ significance. It is not impossible that the >3σ RA discrepancy is caused by the use of different kinds of reference systems. Therefore, we cannot yet conclude that the new PINPT results favor the DMX timing model over the non-DMX model with just one MSP.

Moreover, we reiterate that ${\alpha }_3^{ℬ}$ is subject to additional errors induced by the potential horizontal structure evolution of J2222. A more robust test of the PINPT strategy with pulsar timing can be achieved by comparing $\left({\alpha }_1^{ℬ},{\delta }_1^{ℬ}\right)$ to timing positions (based on different timing models) derived at MJD 59554. A preliminary timing analysis, incorporating a large number of new datasets, shows reasonably good consistency with $\left({\alpha }_1^{ℬ},{\delta }_1^{ℬ}\right)$. Detailed results will be reported as part of an upcoming timing paper. Additionally, applying the PINPT strategy to just a few well-timed MSPs will allow a much more stringent test of the new strategy.

7 Summary

Using the PINPT strategy, we determined an MSP position with a precision of ~0.2 mas using VLBI, which improves on the previous precision level by a factor of approximately 5 (e.g., Ding et al. 2020a; Liu et al. 2023a), and is comparable with the precision level of the timing positions of the best-timed MSPs (e.g., Perera et al. 2019). According to Liu et al. (2023a), applying the PINPT strategy to 50 MSPs promises <0.1 mas precision for the connection between the kinematic and dynamic reference frames, which is more than three times more precise than the previous expectation. Considering that systematic errors of calibrator source positions are in the range of 0.05–0.2 mas¹ and are unlikely to be improved within several decades, our strategy provides a position accuracy that approaches the practical limit.

In general, the PINPT strategy can be used in a broader context: it can sharpen the VLBI localization of any compact steep-spectrum radio source, which could facilitate studies of, for example, fast radio bursts (e.g., Bhandari et al. 2023) or mergers of two neutron stars (e.g., Mooley et al. 2022). Furthermore, it is important to stress that the PINPT strategy is not only meant to enhance the precision of absolute positions, but could have a fundamental impact on differential astrometry as well. Due to the temporal structure evolution of the reference source(s) (which includes core-shift variations; Plavin et al. 2019), proper motions and parallaxes measured with respect to AGNs could be potentially biased (e.g., D13, also see Sect. 6), which is believed to have driven the occasional inconsistencies between VLBI and timing proper motions for the longest timed MSPs (Ding et al. 2023). However, traditional differential astrometry with one calibrator is incapable of quantifying the positional variations due to structure evolution, and is therefore subject to extra systematic errors. In the scenario of precise in-beam astrometry, the variability in the core shifts of the reference sources is a leading source of systematic error. Multi-epoch PINPT observations, on the other hand, can constrain the coreshift variabilities, thus minimizing their detrimental influence on astrometric results. As pulsar astrometry is considered, joining the core-shift-determining sessions to in-beam or Multi-View astrometry can provide accurate pulsar proper motions and parallaxes that remain unbiased by core-shift variabilities, which can then be safely incorporated into pulsar timing analysis.

Acknowledgements

H.D. appreciates the EACOA Fellowship awarded by the East Asia Core Observatories Association. H.D. thanks Wei Zhao for the discussion about jet direction determination. This work is mainly based on observations with the Very Long Baseline Array (VLBA), which is operated by the National Radio Astronomy Observatory (NRAO). The NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Pulse ephemerides of PSR J2222–0137 were made for the purpose of pulsar gating, using data from the Effelsberg 100-meter telescope of the Max-Planck-Institut für Radioastronomie. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under license (Deller et al. 2011).

Appendix A Supporting materials for Sect. 2.2

Fig. A.1 The calibrator plans for the PINPT observations of this work, which involve two in-beam phase calibrators and three out-of-beam calibrators of PSR J2222–0137 (shown in red). More source information can be found in Table A.1. Among the two in-beam calibrators, FIRST J222201–013236, marked by golden circle, was used in previous astrometric campaigns (Deller et al. 2013; Guo et al. 2021); FIRST J222112–012806, marked by blue diamond, is used for core-shift determinations (see Sect. 4.2.2). The adopted phase calibration strategies are described in Sect. 2.2 and summarized in Table 1.

Appendix B Two approaches of Multi-View (2D interpolation)

As mentioned in Sect. 3.1.2, Multi-View was realized in two approaches. In both approaches, J2219 — the closest to PSR J2222–0137 among the off-beam calibrators, was used as the main phase calibrator. In the first approach, we made a phase calibration with J2219, then interpolated the phase solution to the approximate sky position of PSR J2222–0137 by performing two 1D interpolation operations (e.g., Ding et al. 2020b) one after another. Specifically, the phase solution was first extrapolated along the J2219-to-J2226 line to the intersection with the straight line connecting J2218 and PSR J2222–0137. The derived phase solution was then extrapolated again along the J2218-to-pulsar line, to the position of PSR J2222–0137.

In the second approach, the phase solution acquired with J2219 was passed to J2218 and J2226. Subsequently, selfcalibration was performed with both J2218 and J2226. The acquired incremental phase solutions ϕ_J2218 and ϕ_J2226 were then linearly added together as c_J2218⋅ϕ_J2218 + c_J2226⋅ϕ_J2226 (here, c_J2218 and c_J2226 are constants described in Appendix F), which essentially moves the virtual calibrator (explained in Ding et al. 2020b) to the approximate sky position of the target (i.e., PSR J2222–0137, J2221, or J2222).

Appendix C î(є) and σ₀(v)

In Eq. 1, there are two functions — î(є) and σ₀(v). The former describes the fractional atmospheric path length as a function of the antenna elevation є, while the latter characterizes the evolution of the systematic error with respect to the observing frequency v. We define î(є) as l/R_E, where l and R_E are, respectively, the atmospheric path length and the radius of the Earth (not including the atmosphere). It is easy to calculate that $$\widehat{l}(ϵ)=\overline{\sqrt{{\sin }^2ϵ+2{\eta }_{\text{H}}+{\eta }_{\text{H}}^2}}-\sin ϵ,$$(C.1)

where the overline commands averaging over the observation; η_H = h_atmo/R_E is the atmosphere thickness divided by R_E. We adopted η_Η = 0.15 calculated with the upper height of the ionosphere (≈ 965 km) and the average Earth radius of 6371 km.

Using simulations, the relation between propagation-related positional error and v has been studied by Martí-Vidal et al. (2010), and is provided in Fig. 6 of Martí-Vidal et al. (2010). As the simulations of Martí-Vidal et al. (2010) assume a 5° angular separation between the target and the reference source, we divided the results of Martí-Vidal et al. (2010) by a factor of 5 (to reach the systematics per degree separation), and adopted the divided results as σ₀(v). In Eq. 1, we assume ${\sigma }_{ijk}^{\mathcal{S}}\propto \widehat{l}(ϵ)\cdot {\sigma }_0\left(v_j\right)$, while having little knowledge about the coefficient. Therefore, the nuisance parameter η_EFAC (to be determined in Bayesian inference) is essential for completing Eq. 1, and recovering the true magnitude of ${\sigma }_{ijk}^S$.

Appendix D The determination of AGN jet directions

AGN jet directions ϕ_jet and their uncertainties have been directly estimated from their VLBI images (e.g., Kovalev et al. 2017). Likewise, we determined ϕ_jet of the 4 core-shift-probing AGNs from the final image models⁴ of the 4 sources using the newly developed package arcfits⁵. For the ϕ_jet determination, we prioritized the use of image models obtained at 15 GHz because, compared to at lower frequencies, 1) the jet core is closer to the central engine, and 2) the spatial resolution is higher. In case extended jet features were not identified by arcfits, we moved on to the image model of the next highest frequency. The obtained ϕ_jet results are illustrated alongside the AGN images in Fig. D.1, and provided in Table D.1. arcfits derives ϕ_jet and their uncertainties in an automatic way described as follows.

D.1 The noise level of the residual map

After reading a VLBI image, a 2D array of flux densities can be obtained. From this 2D map, the noise level rms of the residual map (i.e., the image after removing all detected source components) can be derived, which is a preparation for measuring ϕ_jet. We established the residual map and estimated the rms in the following way. We first marginalized the 2D array of flux densities into a 1D array, and calculated the standard deviation of the flux densities. All flux densities higher than 7 times the standard deviation were considered “detected”, and were removed from the 1D array. We repeated this standard deviation calculation and removal of detected points until we reached the residual map, in which no more detection can be identified. The standard deviation of this residual map was adopted as the residual map noise level rms.

D.2 Elliptical cuts on the inner regions of VLBI images

The hitherto most advanced method of measuring the position angles of radio features outside the compact radio core is by 1) converting the positions in a VLBI image from Cartesian coordinate to a polar coordinate centered around the brightest spot of the image, and 2) applying circular cuts to the VLBI image (e.g., Cui et al. 2023). As far as core-shift studies are concerned, a ϕ_jet determined close to the compact radio core is expected to better approximate the core-shift direction, as compared to a ϕ_jet determined afar (e.g., Konigl 1981). Due to intrinsic structure and/or scatter broadening, the compact radio core is usually larger than the synthesized beam, while normally resembling the beam (e.g., Fig. D.1). In other words, there is no prior information about the size of the compact radio core. As a result, when the synthesized beam is non-circular, circular cuts applied in close proximity to the compact radio core might lead to the misidentification of the compact radio core (stretching along the major axis) as radio features outside the compact core. Although this anomaly can be corrected by human intervention, other systematic errors might be introduced during the human intervention.

As a novel method that is i) largely free of human intervention, and ii) dedicated to ϕ_jet determination in close proximity to the compact radio core, we cut an image model with a number of increasingly large ellipses that a) have the same axis ratio and position angle as the synthetic beam, and b) are centered at the image pixel of the highest flux density. Given that the position angle and axis ratio are both constant, the n-th (n = 1,2,3,…) ellipse can be characterized by its semi-major axis a_n ; a position on this ellipse can be defined with one additional parameter — the position angle (east of north) ϕ. In this work, we universally adopted $$a_n=a_{\text{beam}}+(n-1)\cdot \left(a_{\text{beam}}/2\right)$$(D.1)

for the ellipses, where a_beam refers to the semi-major axis of the synthesized beam.

Among the ellipses ascertained by Eq. D.1, the innermost ellipse outside the compact radio core was used to derive ϕ_jet; the position angle corresponding to the maximum flux density on the ellipse was adopted as the ϕ_jet. An ellipse is considered outside the compact radio core, when it meets the following criteria:

• (i)

  the maximum flux density on the ellipse is > 7 rms (see Appendix D.1 for the meaning and calculation of rms), and

• (ii)

  the median flux density on the ellipse is < 3 rms, and

• (iii)

  ϕ_jet uncertainty can be calculated with the method detailed in Appendix D.3.

If no ellipse meets the criteria, the AGN is considered a compact source without extended radio features.

D.3 The uncertainty on the AGN jet direction

In Appendix D.2, we adopted the position angle of the maximum flux density 𝒮_max on the innermost elliptical cut outside the compact radio core as the ϕ_jet. Random noises in the image might distort the flux density distribution, and deviate the position of the maximum flux density on the ellipse. The degree of flux density change due to random noises is limited: the chance of large flux density changes due to random noises is lower than smaller flux density changes. Given a flux density drop δ𝒮 induced by random noises, one can estimate the chance of the drop by counting flux densities $<{\overline{\mathcal{S}}}_{\text{R}}-\text{Δ}\mathcal{S}$ in the residual map (see Appendix D.1 for explanation of the residual map), where ${\overline{\mathcal{S}}}_{\text{R}}$ stands for the mean flux density of the residual map. Reversely, provided a confidence level where the flux densities in the residual map is $\ge {\overline{\mathcal{S}}}_{\text{R}}-\text{Δ}\mathcal{S}$, we can derive the δ𝒮 corresponding to the confidence level. In this way, we estimated Δ𝒮 corresponding to 68% confidence level from the residual map. At 68% confidence, we expect ${S^{\prime }}_{\max }>S_{\max }-\text{Δ}\mathcal{S}$, where ${S^{\prime }}_{\max }$ isthe maximum flux density changed by random noises.

On the ellipse (where ϕ_jet is determined), we identify all position angles ϕ where the flux densities equal to S_max − δ𝒮. When exactly two ϕ are acquired, the two ϕ are adopted as the 1σ uncertainty interval of ϕ_jet. Otherwise (which is rare), we consider the ellipse likely still intersects the compact radio core, and move onto the next ellipse further afield.

Fig. D.1 AGN contour plots obtained from the core-shift-determining sessions, and the AGN jet directions determined with arcfits. For the jet direction determinations, image models obtained at higher observing frequencies are preferred, unless extended jet features are not identified by arcfits. The flux density on the m-th (m = 1,2,3,…) contour equals to $3\cdot \text{rms}\cdot {(\sqrt{2})}^{m-1}$ , where rms is 0.18 mJy beam−1, 7.3 mJy beam−1, 0.50 mJy beam−1 and 0.65 mJy beam−1 for J2221, J2218, J2219 and J2226, respectively. Δα and Δδ are, respectively, relative right ascension (left to the east) and relative declination in the unit of image pixel. The reference position of each image is provided over the image. In each panel, the gray ellipse shows the size of the synthesized beam; the AGN jet direction is marked with the orange arrow, while the uncertainty of the direction is illustrated with the two dashed lines.

Appendix E The Bayesian inference of core shifts

We derived the core-shift model with quartet⁶, a newly developed package dedicated to inferring core-shift-related parameters in a Bayesian manner, which is explained as follows.

E.1 Mathematical formalism

Apart from the aforementioned parameters r_i, θ_i, β_i and η_EFAC, the reference positions $x_{ik}^{*}$ are also required in the model of core shifts. For the Bayesian inference, the likelihood function is $$P_{CS}\propto {\left({\displaystyle \prod _i}{\displaystyle \prod _j}{\displaystyle \prod _k}{\sigma }_{ijk}\right)}^{-1}\exp \left[-\frac{1}{2}{\displaystyle \sum _i}{\displaystyle \sum _j}{\displaystyle \sum _k}{\left(\frac{x_{ijk}-{\tilde{x}}_{ijk}}{{\sigma }_{ijk}}\right)}^2\right],$$(E.1)

where x_ijk and ${\tilde{x}}_{ijk}$ refer to, respectively, the observed and the modeled J2221 positions with respect to the i-th reference source at observing frequency ν_j; the total positional uncertainty σ_ijk is the addition-in-quadrature of random and systematic errors. In Eq. E.1, the modeled J2221 positions ${\tilde{x}}_{ijk}$ follow the relation $$\begin{array}{cc}& {\tilde{x}}_{ij\alpha }=\left[r_{\text{J}2221}\left(v_j\right)\cos {\theta }_{\text{J}2221}-r_i\left(v_j\right)\cos {\theta }_i\right]+x_{i\alpha }^{\ast }\\ & {\tilde{x}}_{ij\delta }=\left[r_{\text{J}2221}\left(v_j\right)\sin {\theta }_{\text{J}2221}-r_i\left(v_j\right)\sin {\theta }_i\right]+x_{i\delta }^{\ast },\end{array}$$(E.2)

where the formalism of r_i′ (v_j) is given by Eq. 2.

E.2 Priors

For the Bayesian inference, we adopted the following prior information:

• (i)

  When the inference of β_i′ is requested, the prior constraints of β_i′ follow a uniform distribution between 0.3 and 3, which is denoted as β_i′ ~ 𝓤 (0.3,3).

• (ii)

  r_0i′ ~ 𝓤 (0,5), where the unit is mas.

• (iii)

  $x_{ik}^{*}~\mathcal{U}\left({\overline{x_{ij}}|}_k-{3{\sigma }_{x_{yj}}|}_k,{\overline{x_{ij}}|}_k+{3{\sigma }_{x_{ij}}|}_k\right)$, where ${\overline{x_{ij}}|}_k$ and σ_xij|_k stand for the average and standard deviation of x_ijk over all observing frequencies j; one lower limit and one upper limit are universally used for all observing frequencies.

• (iv)

  η_EFAC ~𝓤 (0, 20).

• (v)

  The prior constraints on θ_i′ follow Gaussian distributions characterized by the AGN jet directions ϕ_jet in Table D.1. Namely, θ_i′ ~ 𝒢 (ϕ_jet|_i′ ,σϕ_jet|_i′). In the case of asymmetric ϕ_jet uncertainty, the larger side of the uncertainty is used as the σϕ_jet|_i′.

Appendix F The calculation of an absolute pulsar position and its uncertainty

The determination of the absolute position of PSR J2222–0137 through the PINPT strategy relies on well determined absolute positions of the off-beam calibrators. This work is based on the off-beam calibrator positions $x_i^{\text{RFC}}$ reported in the 2024A release of the Radio Fundamental Catalogue¹ (RFC). Generally speaking, for Multi-View carried out with 3 off-beam calibrators, the relation between the target position x_target and the off-beam calibrator positions x_q (q = 1,2, 3) is $$x_{\text{target}}={\displaystyle \underset{3}{\overset{q}{{\sum }^{\text{}}}}}{c}_q{x}_q\text{,}$$(F.1)

where x_target is known to a precision good enough to guide the Multi-View; c_q are coefficients that can be solved with the additional condition $${\displaystyle \sum _q^3}{c}_q=1.$$(F.2)

Therefore, in Multi-View, offsets in the off-beam calibrator positions would lead to offset in the target position (refined with Multi-View), as $$\text{Δ}{x}_{\text{target}}={\displaystyle \sum _q^3}{c}_q\cdot \text{Δ}{x}_q,$$(F.3)

where Δx_q result from 1) core shifts of the off-beam calibrators, and 2) inaccurate positions of the off-beam calibrator image models, and are calculated by $$\text{Δ}{x}_q=\left(x_q^{\text{model}}-x_q^{\text{RFC}}\right)+\left[r_q\left(v_{psr}\right)-r_q^{\text{RFC}}\right].$$(F.4)

Here, $x_q^{\text{model }}$ refers to the image model position of the q-th off- beam calibrator; r_q (v_psr) denotes the core shift of the q-th off- beam calibrator at the observing frequency of interest (which is 1.6 GHz, the observing frequency of the pulsar sessions, for this work); $r_q^{\text{RFC}}$ is the residual core shift of the RFC position.

Specific to this work, q = i = J2212, J2219, J2226. c_i for various targets are provided in Table F.1. According to Porcas (2009), $r_i^{\text{RFC}}=0$, when 1) $x_i^{\text{RFC}}$ is derived with group-delay astrometry (that removes the ionosphere-induced group delay with dual-band or multi-band geodetic observations), and 2) β_i = 1. Both conditions are met for this work: $x_i^{\text{RFC}}$ are estimated with group-delay astrometry; and we already assume that β_i ≡ 1 in the Bayesian analysis (see Sect. 4.2.3). Therefore, we adopt $r_i^{\text{RFC}}=0$ in this work.

The uncertainty on Δx_target is calculated by the addition-in- quadrature of the uncertainties on the components c_iΔx_i, where the uncertainty of Δx_i is further derived by the addition-in- quadrature of the uncertainties on $x_i^{\text{RFC}}$ and r_i.

F.1 Deriving an absolute pulsar position via a nearby AGN

When the target itself is an AGN (e.g., J2221, J2222), Eq. F.3 needs to be generalized to $$\text{Δ}{x}_{\text{target}}={\displaystyle \sum _q^3}{c}_q\cdot \text{Δ}{x}_q+\text{Δ}{r}_{\text{target}},$$(F.5)

where Δr_target is the core-shift difference between v_psr and the observing frequency v_IBC of the AGN target. Accordingly, the uncertainty on Δx_target is further added in quadrature by the uncertainty of Δr_target.

The availability of a nearby AGN (hereafter referred to as IBC) around a pulsar (or other targets of interest) provides an alternative pathway to the absolute position of the pulsar. After applying the position correction Δx_target calculated with Eq. F.5, the absolute position x_IBC of the IBC at v_psr can be derived. Provided i) the image model position $x_{\text{IBC}}^{\text{model}}$ of the IBC, and ii) the pulsar position $x_{\text{psr}}^{\text{PR}}$ determined with phase referencing with respect to the IBC, the absolute pulsar position can be calculated as $$x_{\text{psr}}=x_{\text{psr}}^{\text{PR}}+\left(x_{\text{IBC}}-x_{\text{IBC}}^{\text{model}}\right).$$(F.6)

This pathway of deriving the absolute position may become the only option in cases where Multi-View of the target of interest cannot be arranged (e.g., when the target is an unpredictable radio transient). The uncertainty of x_psr derived via the IBC is the addition-in-quadrature of the Δx_target uncertainty, the $x_{\text{psr}}^{\text{PR}}$ uncertainty (estimated with Bayesian inference in this work), and the statistical uncertainty of the IBC position obtained with Multi-View at v_IBC.

References

All Tables

All Figures

Fig. 1 Jet core relative positions of FIRST J222112-012806 measured with respect to the jet cores of three reference sources (i.e., ICRF J221852.0-033536, ICRF J221947.2-005132, and ICRF J222646.5+005211) at five observing wavelengths. Systematic uncertainties calculated with Eq. (1) have been included in the positional uncertainties. Models of the combined core shifts are derived assuming ${\beta }_{i^{\prime }}\equiv 1$ with Bayesian inference (see Appendix E), and displayed with straight lines.

In the text

Fig. A.1 The calibrator plans for the PINPT observations of this work, which involve two in-beam phase calibrators and three out-of-beam calibrators of PSR J2222–0137 (shown in red). More source information can be found in Table A.1. Among the two in-beam calibrators, FIRST J222201–013236, marked by golden circle, was used in previous astrometric campaigns (Deller et al. 2013; Guo et al. 2021); FIRST J222112–012806, marked by blue diamond, is used for core-shift determinations (see Sect. 4.2.2). The adopted phase calibration strategies are described in Sect. 2.2 and summarized in Table 1.

In the text

Fig. D.1 AGN contour plots obtained from the core-shift-determining sessions, and the AGN jet directions determined with arcfits. For the jet direction determinations, image models obtained at higher observing frequencies are preferred, unless extended jet features are not identified by arcfits. The flux density on the m-th (m = 1,2,3,…) contour equals to $3\cdot \text{rms}\cdot {(\sqrt{2})}^{m-1}$ , where rms is 0.18 mJy beam−1, 7.3 mJy beam−1, 0.50 mJy beam−1 and 0.65 mJy beam−1 for J2221, J2218, J2219 and J2226, respectively. Δα and Δδ are, respectively, relative right ascension (left to the east) and relative declination in the unit of image pixel. The reference position of each image is provided over the image. In each panel, the gray ellipse shows the size of the synthesized beam; the AGN jet direction is marked with the orange arrow, while the uncertainty of the direction is illustrated with the two dashed lines.

In the text