# The MSPSRπ catalogue: VLBA astrometry of 18 millisecond pulsars

H. Ding, A.T. Deller, B.W. Stappers, T.J.W. Lazio, D. Kaplan, S. Chatterjee, W. Brisken, J. Cordes, P.C.C. Freire, E. Fonseca, et al.

► To cite this version:

H. Ding, A.T. Deller, B.W. Stappers, T.J.W. Lazio, D. Kaplan, et al. The MSPSRπ catalogue: VLBA astrometry of 18 millisecond pulsars. Monthly Notices of the Royal Astronomical Society, 2023, 519 (4), pp.4982-5007. 10.1093/mnras/stac3725. hal-03925284

HAL Id: hal-03925284

https://hal.science/hal-03925284

Submitted on 22 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# The MSPSRπ catalogue: VLBA astrometry of 18 millisecond pulsars

H. Ding<sup>1,2,*</sup> A. T. Deller<sup>1,2*</sup> B. W. Stappers<sup>3</sup> T. J. W. Lazio,<sup>4</sup> D. Kaplan<sup>5</sup> S. Chatterjee<sup>1,6</sup>
W. Brisken,<sup>7</sup> J. Cordes,<sup>6</sup> P. C. C. Freire<sup>8</sup> E. Fonseca,<sup>9,10</sup> I. Stairs,<sup>11</sup> L. Guillemot,<sup>12,13</sup> A. Lyne,<sup>3</sup>
I. Cognard,<sup>12,13</sup> D. J. Reardon<sup>1,2</sup> and G. Theureau<sup>12,13,14</sup>

Accepted 2022 December 12. Received 2022 November 10; in original form 2022 June 24

## ABSTRACT

With unparalleled rotational stability, millisecond pulsars (MSPs) serve as ideal laboratories for numerous astrophysical studies, many of which require precise knowledge of the distance and/or velocity of the MSP. Here, we present the astrometric results for 18 MSPs of the 'MSPSRπ' project focusing exclusively on astrometry of MSPs, which includes the re-analysis of three previously published sources. On top of a standardized data reduction protocol, more complex strategies (i.e. normal and inverse-referenced 1D interpolation) were employed where possible to further improve astrometric precision. We derived astrometric parameters using *sterne*, a new Bayesian astrometry inference package that allows the incorporation of prior information based on pulsar timing where applicable. We measured significant (>$3\sigma$) parallax-based distances for 15 MSPs, including $0.81 \pm 0.02$ kpc for PSR J1518+4904 – the most significant model-independent distance ever measured for a double neutron star system. For each MSP with a well-constrained distance, we estimated its transverse space velocity and radial acceleration. Among the estimated radial accelerations, the updated ones of PSR J1012+5307 and PSR J1738+0333 impose new constraints on dipole gravitational radiation and the time derivative of Newton's gravitational constant. Additionally, significant angular broadening was detected for PSR J1643−1224, which offers an independent check of the postulated association between the HII region Sh 2-27 and the main scattering screen of PSR J1643−1224. Finally, the upper limit of the death line of γ-ray-emitting pulsars is refined with the new radial acceleration of the hitherto least energetic γ-ray pulsar PSR J1730−2304.

**Key words:** gravitation – stars: kinematics and dynamics – pulsars: individual: PSR J0030+0451, PSR J0610−2100, PSR J0621+1002, PSR J1024−0719, PSR J1537+1155, PSR J1853+1303, PSR J1910+1256, PSR J1918−0642, PSR J1939+2134 – gamma-rays: stars – radio continuum: stars.

## 1 INTRODUCTION

### 1.1 Millisecond pulsars: a key for probing theories of gravity and detecting the gravitational-wave background

Pulsars are an observational manifestation of neutron stars (NSs) that emit non-thermal electromagnetic radiation while spinning (Gold 1968; Pacini 1968; Hewish et al. 1969). Over 3000 radio pulsars

have been discovered to date throughout the Galaxy and the nearest members of the Local Group (Manchester et al. 2005). Due to the large moment of inertia of pulsars, the pulses we receive on the Earth from a pulsar exhibit highly stable periodicity. By measuring a train of pulse time-of-arrivals (ToAs) of a pulsar and comparing it against the model prediction, a long list of model parameters can be inferred (e.g. Detweiler 1979; Helfand et al. 1980). This procedure to determine ToA-changing parameters is known as pulsar timing, hereafter referred to as timing.

In the pulsar family, recycled pulsars (commonly refereed to as millisecond pulsars, or MSPs) have the shortest rotational periods.

<footnote n="1"><sup>1</sup>Centre for Astrophysics and Supercomputing, Swinburne University of Technology, John St, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="2"><sup>2</sup>ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Australia</footnote>

<footnote n="3"><sup>3</sup>Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="4"><sup>4</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Blvd, Pasadena, CA 91109, USA</footnote>

<footnote n="5"><sup>5</sup>Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA</footnote>

<footnote n="6"><sup>6</sup>Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="7"><sup>7</sup>National Radio Astronomy Observatory, PO Box O, Socorro NM 87801, USA</footnote>

<footnote n="8"><sup>8</sup>Max-Planck-Institut für Radioastronomie, auf dem Hügel 69, Bonn D-53121, Germany</footnote>

<footnote n="9"><sup>9</sup>Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506-6315, USA</footnote>

<footnote n="10"><sup>10</sup>Center for Gravitational Waves and Cosmology, Chestnut Ridge Research Building, Morgantown, WV 26505, USA</footnote>

<footnote n="11"><sup>11</sup>Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="12"><sup>12</sup>Observatoire Radioastronomique de Nançay, Observatoire de Paris, Université PSL, Université d'Orléans, CNRS, F-18330 Nançay, France</footnote>

<footnote n="13"><sup>13</sup>Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, Université d'Orléans/CNRS, F-45071 Orléans Cedex 02, France</footnote>

<footnote n="14"><sup>14</sup>LUTH, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, F-92195 Meudon, France</footnote>

<footnote>* E-mail: hdingastro@hotmail.com (HD); adeller@astro.swin.edu.au (ATD)</footnote>
They are believed to have been spun-up through the accretion from their donor stars during a previous evolutionary phase as a low-mass X-ray binary (LMXB) (Alpar et al. 1982). As the duration of the recycling phase (and hence the degree to which the pulsar is spun-up) can vary depending on the nature of the binary, there is no clear spin period threshold that separates MSPs from canonical pulsars. In this paper, we define MSPs as pulsars with spin periods of $\le 40$ ms and magnetic fields $\le 10^{10}$ G. This range encompasses most partially recycled pulsars with NS companions, such as PSR J1537+1155 (also known as PSR B1534+12) and PSR J1518+4904. Compared to non-recycled pulsars, ToAs from MSPs can be measured to higher precision due to both the narrower pulse profiles and larger number of pulses. Additionally, MSPs exhibit more stable rotation (e.g. Hobbs, Lyne & Kramer 2010); both factors promise a lower level of random timing noise. Consequently, MSPs outperform non-recycled pulsars in the achievable precision for probing theories underlying ToA-changing astrophysical effects. In particular, MSPs provide the hitherto most precise tests for gravitational theories (e.g. Freire et al. 2012; Zhu et al. 2019; Kramer et al. 2021). Einstein's theory of general relativity (GR) is the simplest form among a group of possible candidate post-Newtonian gravitational theories. The discovery of highly relativistic double neutron star (DNS) systems (e.g. Hulse & Taylor 1975; Wolszczan 1991; Burgay et al. 2003; Lazarus et al. 2016; Cameron et al. 2018; Stovall et al. 2018), and their continued timing have resulted in many high-precision tests of GR and other gravity theories (Fonseca, Stairs & Thorsett 2014; Weisberg & Huang 2016; Ferdman et al. 2020, and especially Kramer et al. 2021). The precise timing, optical spectroscopy and VLBI observations of pulsar-white-dwarf (WD) systems have, in addition, achieved tight constraints on several classes of alternative theories of gravity (Deller et al. 2008; Lazaridis et al. 2009; Freire et al. 2012; Antoniadis et al. 2013; Ding et al. 2020b; Guo et al. 2021; Zhao et al. 2022).

Gravitational Waves (GWs) are changes in the curvature of spacetime (generated by accelerating masses), which propagate at the speed of light. Individual GW events in the Hz-kHz range have been detected directly with GW observatories (e.g. Abbott et al. 2016; see the third Gravitational-Wave Transient Catalog<fnref n="1" />), and indirectly using the orbital decay of pulsar binaries (e.g. Taylor & Weisberg 1982; Weisberg & Huang 2016; Kramer et al. 2021; Ding et al. 2021a). Collectively, a gravitational wave background (GWB) formed with primordial GWs and GWs generated by later astrophysical events (Carr 1980) is widely predicted, but has not yet been confirmed by any observational means. In the range of $10^{-9}$–0.1 Hz, supermassive black hole binaries are postulated to be the primary sources of the GWB (Sesana, Vecchio & Colacino 2008). In this nano-hertz regime, the most stringent constraints on the GWB are provided by pulsar timing (Detweiler 1979).

To enhance the sensitivity for the GWB hunt with pulsar timing, and to distinguish GWB-induced ToA signature from other sources of common timing ‘noise’ (e.g. Solar system planetary ephemeris error, clock error and interstellar medium, Tiburzi et al. 2016), a pulsar timing array (PTA), composed of MSPs scattered across the sky (see Roebber 2019 for spatial distribution requirement), is necessary (Foster & Backer 1990). After two decades of efforts, no GWB has yet been detected by a PTA, though common steep-spectrum timing noise (in which GWB signature should reside) has already been confirmed by several radio PTA consortia (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al. 2021; Antoniadis et al. 2022). At $\gamma$-rays, a competitive GWB amplitude upper limit was recently

achieved using the Fermi Large Area Telescope with 12.5 years of data (Fermi-LAT Collaboration 2022).

## 1.2 Very long baseline astrometry of millisecond pulsars

In timing analysis, astrometric information for an MSP (reference position, proper motion, and annual geometric parallax) can form part of the global ensemble of parameters determined from ToAs. However, the astrometric signatures can be small compared to the ToA precision and/or covariant with other parameters in the model, especially for new MSPs that are timed for less than a couple of years (Madison, Chatterjee & Cordes 2013). Continuing to add newly discovered MSPs into PTAs is considered the best pathway to rapidly improve the PTA sensitivity (Siemens et al. 2013), and is particularly important for PTAs based around newly commissioned high-sensitivity radio telescopes (e.g. Bailes et al. 2020). Therefore, applying priors to the astrometric parameters can be highly beneficial for the timing of individual MSPs (especially the new ones) and for enhancing PTA sensitivities (Madison et al. 2013).

Typically, the best approach to independently determine precise astrometric parameters for MSPs is the use of phase-referencing (e.g. Lestrange et al. 1990; Beasley & Conway 1995) very long baseline interferometry (VLBI) observations, which can achieve submas positional precision (relative to a reference source position) for MSPs in a single observation. By measuring the sky position of a Galactic MSP a number of times and modelling the position evolution, VLBI astrometry can obtain astrometric parameters for the MSP. Compared to pulsar timing, VLBI astrometry normally takes much shorter time to reach a given astrometric precision (e.g. Brisken et al. 2002; Chatterjee et al. 2009; Deller et al. 2019).

One of the limiting factors on searching for the GWB with PTAs is the uncertainties on the Solar system planetary ephemerides (SSEs) (Vallisneri et al. 2020), which are utilized to convert geocentric ToAs to ones measured in the (Solar system) barycentric frame (i.e. the reference frame with respect to the barycentre of the Solar system). Various space-mission-driven SSEs have been released mainly by two SSE providers – the NASA Jet Propulsion Laboratory (e.g. Park et al. 2021) and the IMCCE (e.g. Fienga, Avdellidou & Hanuš 2020). In pulsar timing analysis adopting different SSEs may lead to discrepant timing parameters (e.g. Wang et al. 2017). On the other hand, VLBI astrometry measures offsets with respect to a source whose position is measured in a quasi-inertial (reference) frame defined using remote quasars (e.g. Charlot et al. 2020). Although VLBI astrometry also relies on SSEs to derive annual parallax, it is robust against SSE uncertainties. In other words, for VLBI astrometry, using different SSEs in parameter inference would not lead to a noticeable difference in the inferred parameters. Therefore, VLBI astrometry of MSPs can serve as an objective standard to be used to discriminate between various SSEs. Specifically, if an SSE is inaccurate, the barycentric frame based on the SSE would display rotation with respect to the quasar-based frame. This frame rotation can be potentially detectable by comparing VLBI positions of multiple MSPs against their timing positions (Chatterjee et al. 2009; Wang et al. 2017). By eliminating inaccurate SSEs, VLBI astrometry of MSPs can suppress the SSE uncertainties, and hence enhance the PTA sensitivities.

Besides the GWB-related motivations, interferometer-based astrometric parameters (especially distances to MSPs) have been adopted to sharpen the tests of gravitational theories for individual MSPs (e.g. Deller, Bailes & Tingay 2009; Deller et al. 2018; Guo et al. 2021; Ding et al. 2021a). Such tests are normally made by comparing the model-predicted and observed post-Keplerian (PK) parameters

<footnote n="1">¹https://www.ligo.org/science/Publication-O3aFinalCatalog/</footnote>
that quantify excessive gravitational effects beyond a Newtonian description of the orbital motion. Among the PK parameters is the orbital decay $\dot{P}_b$ (or the time derivative of orbital period). The intrinsic cause of $\dot{P}_b$ in double neutron star systems is dominated by the emission of gravitational waves, which can be predicted using the binary constituent masses and orbital parameters (e.g. Lazaridis et al. 2009; Weisberg & Huang 2016). To test this model prediction, however, requires any extrinsic orbital decay $\dot{P}_b^{\text{ext}}$ due to relative acceleration between the pulsar and the observer to be removed from the observed $\dot{P}_b$. Such extrinsic terms depend crucially on the proper motion and the distance of the pulsar, however these (especially the distance) can be difficult to estimate from pulsar timing. Precise VLBI determination of proper motions and distances can yield precise estimates of these extrinsic terms, and therefore play an important role in orbital-decay tests of gravitational theories. Likewise, Gaia astrometry on nearby pulsar-WD systems can potentially serve the same scientific goal though the method is only applicable to a small number of pulsar-WD systems, where the WDs are sufficiently bright for the Gaia space observatory (see Section 5.2).

Last but not least, pulsar astrometry is crucial for understanding the Galactic free-electron distribution, or the Galactic free-electron number density $n_e(\vec{x})$ as a function of position. An $n_e(\vec{x})$ model is normally established by using pulsars with well determined distances as benchmarks. As the pulsations from a pulsar allow precise measurement of its dispersion measure (DM), the average $n_e$ between the pulsar and the Earth can be estimated given the pulsar distance. Accordingly, a large group of such benchmark pulsars across the sky would enable the establishment of an $n_e(\vec{x})$ model. In a relevant research field, extragalactic fast radio bursts (FRBs) have been used to probe intergalactic medium distribution on a cosmological scale (e.g. Macquart et al. 2020; Mannings et al. 2021), which, however, demands the removal of the DMs of both the Galaxy and the FRB host galaxy. The Galactic DM cannot be determined without a reliable $n_e(\vec{x})$ model, which, again, calls for precise astrometry of pulsars across the Galaxy.

## 1.3 The MSPSRπ project

Using the Very Long Baseline Array (VLBA), the PSRπ project tripled the sample of pulsars with precisely measured astrometric parameters (Deller et al. 2019), but included just three MSPs. The successor project, MSPSRπ, is a similarly designed VLBA astrometric program targeting exclusively MSPs. Compared to canonical pulsars, MSPs are generally fainter. To identify MSPs feasible for VLBA astrometry, a pilot program was conducted, which found 31 suitable MSPs. Given observational time constraints, we selected 18 MSPs as the targets of the MSPSRπ project, focusing primarily on sources observed with pulsar timing arrays. The 18 MSPs are listed in Table 1 along with their spin periods $P_s$ and orbital periods $P_b$ (if available) that have been obtained from the ATNF Pulsar Catalogue<fnref n="2" /> (Manchester et al. 2005). The astrometric results for three sources (PSR J1012+5307, PSR J1537+1155, PSR J1640+2224) involved in the project have already been published (Vigeland et al. 2018; Ding et al. 2020b, 2021a). In this paper, we present the astrometric results of the remaining 15 MSPs studied in the MSPSRπ project. We also re-derived the results for the three published MSPs, in order to ensure consistent and systematic astrometric studies.

Along with the release of the catalogue results, this paper covers several scientific and technical perspectives. First, this paper explores novel data reduction strategies such as inverse-referenced 1D phase interpolation (see Section 3.2). Second, a new Bayesian astrometry inference package is presented (see Section 4). Third, with new parallax-based distances and proper motions, we discriminate between the two prevailing $n_e(\vec{x})$ models (see Section 6.1.1), and investigate the kinematics of MSPs in Section 6.2. Fourth, with new parallax-based distances of two MSPs, we re-visit the constraints on alternative theories of gravity (see Section 7). Finally, discussions on individual pulsars are given in Section 8, which includes a refined ‘death line’ upper limit of γ-ray pulsars (see Section 8.7). The study of SSE-dependent frame rotation, which depends on an accurate estimation of the reference points of our calibrator sources in the quasi-inertial VLBI frame, requires additional multifrequency observations and will be presented in a follow-up paper.

Throughout this paper, we abide by the following norms unless otherwise stated. (1) The uncertainties are provided at 68 per cent confidence level. (2) Any mention of flux density refers to unresolved flux density $S_{\text{unres}}$ in our observing configuration (e.g. a 10-mJy source means $S_{\text{unres}} = 10$ mJy). (3) All bootstrap and Bayesian results adopt the 50th, 16th, and 84th percentile of the marginalized (and sorted) value chain as, respectively, the estimate and its $1-\sigma$ error lower and upper bound. (4) Where an error of an estimate is required for a specific calculation but an asymmetric error is reported for the estimate, the mean of upper and lower errors is adopted for the calculation. (5) VLBI positional uncertainties will be broken down into the uncertainty of the offset from a chosen calibrator reference point, and the uncertainty in the location of that chosen reference point. This paper focuses on the relative offsets, which are relevant for the measurement of proper motion and parallax, and the uncertainty in the location of the reference source is presented separately.

## 2 OBSERVATIONS AND CORRELATION

As is mentioned in Section 1.2 to achieve high-precision pulsar astrometry requires the implementation of a VLBI phase referencing technique. There are, however, a variety of such techniques, including the normal phase referencing, relayed phase referencing, inverse phase referencing, and interpolation. These techniques are described and discussed in Chapter 2 of Ding 2022. Generally, a given phase referencing approach and hence observational set-up maps directly to a corresponding data reduction procedure though occasionally other data reduction opportunities could arise by chance (see Section 3).

The MSPSRπ project systematically employs the relayed phase referencing technique, in which a secondary phase reference source (explained in Chapter 2 of Ding 2022) very close to the target on the sky is observed to refine direction-dependent calibration effects. The observing and correlation tactics are identical to those of the PSRπ project (Deller et al. 2019). All MSPs in the MSPSRπ catalogue (see Table1) were observed at L-band with the VLBA at 2-Gbps data rate (256 MHz total bandwidth, dual polarization) from mid-2015 to no later than early 2018. To minimize radio-frequency interference (RFI) at L-band, we used eight 32 MHz subbands with central frequencies of 1.41, 1.44, 1.47, 1.50, 1.60, 1.66, 1.70, and 1.73 GHz, corresponding to an effective central frequency of 1.55 GHz. The primary phase calibrators were selected from the Radio Fundamental Catalogue.<fnref n="3" /> The secondary phase calibrators were identified from the FIRST (Faint Images of the Radio Sky at

<footnote n="2">²https://www.atnf.csiro.au/research/pulsar/psrcat/</footnote>

<footnote n="3">³astrogeo.org/rfc/</footnote>
<table><caption>**Table 1.** List of pulsars and phase calibrators.</caption><thead><tr><th>PSR</th><th>P<sub>s</sub><br>(ms)</th><th>P<sub>b</sub><br>(d)</th><th>Gating<br>gain</th><th>Project codes</th><th>Primary phase calibrator</th><th>ΔPC-IBC<br>(deg)</th><th>f</th><th>Secondary phase calibrator</th><th>IBC<sup>a</sup><br>code</th><th>Δp<sub>pr-IBC</sub><sup>b</sup><br>(arcmin)</th><th>S<sup>b</sup>IBC<br>units<br>(mJy)</th></tr></thead><tbody><tr><td>J0030+0451</td><td>4.87</td><td>—</td><td>1.75</td><td>BD179B, BD192B</td><td>ICRF J002945.8+055440</td><td>0.97</td><td>—</td><td>FIRST J003054.6+045908</td><td>00027</td><td>10.1</td><td>20.6</td></tr><tr><td>J0610-2100</td><td>3.86</td><td>0.29</td><td>1.85</td><td>BD179C, BD192C</td><td>—</td><td>—</td><td>—</td><td>NVSS J061002-211538<sup>c</sup></td><td>00238</td><td>15.4</td><td>112.8</td></tr><tr><td>J0621+1002</td><td>28.85</td><td>8.3</td><td>2.15</td><td>BD179D, BD192D</td><td>ICRF J061909.9+073641</td><td>2.84</td><td>—</td><td>NVSS J062153+102206</td><td>00303</td><td>21.0⇒2.9<sup>d</sup></td><td>18.1</td></tr><tr><td>J1012+5307**</td><td>5.26</td><td>0.60</td><td>1.43</td><td>BD179E, BD192E</td><td>ICRF J095837.8+503957</td><td>3.38</td><td>—</td><td>NVSS J101307+531233</td><td>00462</td><td>7.51</td><td>20.8</td></tr><tr><td>J1024-0719</td><td>5.16</td><td>—</td><td>1.69</td><td>BD179E, BD192E</td><td>ICRF J102838.7-084438</td><td>1.59</td><td>—</td><td>FIRST J102526.3-072216</td><td>00529</td><td>12.2</td><td>11.7</td></tr><tr><td>J1518+4904</td><td>40.93</td><td>8.6</td><td>2.27</td><td>BD179F, BD192F</td><td>ICRF J150644.1+493355</td><td>1.78</td><td>—</td><td>NVSS J151733+491626</td><td>00691</td><td>13.8⇒6.5<sup>g</sup></td><td>35.6</td></tr><tr><td>J1537+1155**</td><td>37.90</td><td>0.42</td><td>1.62</td><td>BD179F, BD192F, BD229</td><td>ICRF J154049.4+144745</td><td>3.18</td><td>—</td><td>FIRST J153746.2+114215</td><td>00840</td><td>16.3</td><td>19.2</td></tr><tr><td>J1640+2224**</td><td>3.16</td><td>175</td><td>1.90</td><td>BD179G, BD192F</td><td>ICRF J164125.2+225704</td><td>0.79</td><td>—</td><td>NVSS J164018+221203</td><td>00920</td><td>12.1</td><td>98.0</td></tr><tr><td>J1643-1224</td><td>4.62</td><td>147</td><td>1.56</td><td>BD179G, BD192G</td><td>ICRF J163845.2-141550</td><td>2.49</td><td>—</td><td>NVSS J164515-122013</td><td>01120</td><td>24.3</td><td>6.0</td></tr><tr><td>J1721-2457</td><td>3.50</td><td>—</td><td>1.78</td><td>BD179H, BD192H</td><td>ICRF J172658.9-225801</td><td>2.47</td><td>—</td><td>NVSS J172129-250538</td><td>01188</td><td>10.1</td><td>7.4</td></tr><tr><td>J1730-2304</td><td>8.12</td><td>—</td><td>1.86</td><td>BD179H, BD192H</td><td>ICRF J172658.9-225801</td><td>0.76</td><td>—</td><td>NVSS J172932-232722</td><td>01220</td><td>25.5</td><td>37.9</td></tr><tr><td>J1738+0333</td><td>5.85</td><td>0.35</td><td>1.98</td><td>BD179I, BD192I</td><td>ICRF J174037.1+031147</td><td>0.66</td><td>—</td><td>NVSS J173823+033305</td><td>01313</td><td>7.5</td><td>6.5</td></tr><tr><td>J1824-2452A</td><td>3.05</td><td>—</td><td>1.53</td><td>BD179Q, BD192Q</td><td>ICRF J182057.8-252812</td><td>0.61</td><td>—</td><td>NVSS J182301-250438</td><td>01433</td><td>23.8⇒3.7<sup>d</sup></td><td>3.3</td></tr><tr><td>J1853+1303</td><td>4.09</td><td>116</td><td>1.64</td><td>BD179K, BD192K</td><td>ICRF J185250.5+142639</td><td>1.51</td><td>—</td><td>NVSS J185456+130110</td><td>01535</td><td>14.6</td><td>5.9</td></tr><tr><td>J1910+1256</td><td>4.98</td><td>58.5</td><td>2.65</td><td>BD179L, BD192L</td><td>ICRF J191158.2+161146</td><td>3.16</td><td>—</td><td>NVSS J190957+130434</td><td>01769</td><td>8.7</td><td>4.9</td></tr><tr><td>J1911-1114</td><td>3.63</td><td>2.7</td><td>1.67</td><td>BD179M, BD192M</td><td>ICRF J190528.5-115332</td><td>1.76</td><td>—</td><td>NVSS J191233-113327</td><td>01816</td><td>21.9</td><td>23.4</td></tr><tr><td>J1918-0642</td><td>7.65</td><td>10.9</td><td>2.12</td><td>BD179M, BD192M</td><td>ICRF J191207.1-080421</td><td>2.14</td><td>—</td><td>NVSS J191731-062435</td><td>01846</td><td>26.1</td><td>50.7</td></tr><tr><td>J1939+2134</td><td>1.56</td><td>—</td><td>1.35</td><td>BD179K, BD192K</td><td>ICRF J193510.4+203154</td><td>1.88</td><td>—</td><td>NVSS J194104+214913<sup>e</sup></td><td>01647</td><td>24.5⇒4.1</td><td>10.6</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>NVSS J194106+215304<sup>e</sup></td><td>01648</td><td>27.4⇒2.7</td><td>1.8</td></tr></tbody></table>

Notes. •The image models for the primary and secondary calibrators listed here are publicly available<sup>6</sup>.

*P<sub>s</sub> and P<sub>b</sub> stand for spin period and orbital period, respectively.

**Unresolved flux intensity of the secondary phase calibrator at 1.55 GHz.

***Published in Ding et al. (2020b), Ding et al. (2021a), Vigeland et al. (2018).

aSecondary phase calibrators are named IBCXXXX in the BD179 and BD192 observing files, where 'XXXX' represents a unique 5-digit IBC code.

bAngular separation between target and secondary calibrator.

cNVSS J061002-211538, close to the pulsar on the sky, is bright enough to serve as primary phase calibrator.

dAs ID interpolation is applied, the pulsar-to-virtual-calibrator separation is also provided after '⇒' (see Section 3.1).

eHere, inverse phase referencing is adopted, where the 'secondary phase calibrators' are ultimately the targets (see Section 3.2).

fAngular separation between primary and secondary calibrator.

gThe NVSS J151815+491105, a 4.5-mJy-bright source 6.5 away from the pulsar, is used as the final reference source (see Section 3).
Twenty-cm) catalogue (Becker, White & Helfand 1995) or the NVSS (NRAO VLA sky survey) catalogue (Condon et al. 1998) (for sky regions not covered by the FIRST survey) using a short multifield observation. Normally, more than one secondary phase calibrators were observed together with the target. Among them, a main one that is preferably the brightest and the closest to the target is selected to carry out self-calibration; the other secondary phase calibrators are hereafter referred to as redundant secondary phase calibrators. The primary and the main secondary phase calibrators for the astrometry of the 18 MSPs are summarized in Table 1, alongside the project codes. At correlation time, pulsar gating was applied (Deller et al. 2011) to improve the S/N on the target pulsars. The median values of the gating gain, defined as $(S/N)_{gated}/(S/N)_{ungated}$, are provided in Table 1.

## 3 DATA REDUCTION AND FIDUCIAL SYSTEMATIC ERRORS

We reduced all data with the psrvlbireduce pipeline<fnref n="4" /> written in parseltongue (Kettenis et al. 2006), a python-based interface for running functions provided by AIPS (Greisen 2003) and DIFMAP (Shepherd, Pearson & Taylor 1994). The procedure of data reduction is identical to that outlined in Ding et al. (2020b), except for four MSPs – PSR J1518+4904, PSR J0621+1002, PSR J1824–2452A, and PSR J1939+2134. For PSR J1518+4904, the self-calibration solutions acquired with NVSS J151733+491626, a 36-mJy secondary calibrator 13:8 away from the pulsar, are extrapolated to both the pulsar and NVSS J151815+491105 – a 4.5-mJy source about a factor of two closer to PSR J1518+4904 than NVSS J151733+491626. The positions relative to NVSS J151815+491105 are used to derive the astrometric parameters of PSR J1518+4904. For the other exceptions, the data reduction procedures as well as fiducial systematics estimation are described in Sections 3.1 and 3.2.

At the end of the data reduction, a series of positions as well as their random errors $\sigma_i^R$ (where $i = 1, 2, 3, ...$ refers to right ascension or declination at different epochs) are acquired for each pulsar. For each observation, on top of the random errors due to image noise, ionospheric fluctuations would introduce systematic errors that distort and translate the source, the magnitude of which generally increases with the angular separation between a target and its (secondary) phase calibrator (e.g. Chatterjee et al. 2004; Pradel, Charlot & Lestrade 2006; Kirsten et al. 2015; Deller et al. 2019). We estimate fiducial values for these systematic errors $\sigma_i^S$ of pulsar positions using the empirical relation (i.e. equation 1 of Deller et al. 2019) derived from the whole PSR$\pi$ sample. While this empirical relation has proven a reasonable approximation to the actual systematic errors for a large sample of sources, for an individual observational set-up $\sigma_i^S$ may overstate or underestimate the true systematic error (see Section 4). We can account for our uncertainty in this empirical estimator by re-formulating the total positional uncertainty as

$$ \sigma_i(\eta_{\text{EFAC}}) = \sqrt{(\sigma_i^R)^2 + (\eta_{\text{EFAC}} \cdot \sigma_i^S)^2}, \quad (1) $$

where $\eta_{\text{EFAC}}$ is a positive correction factor on the fiducial systematic errors. In this work, we assume $\eta_{\text{EFAC}}$ stays the same for each pulsar throughout its astrometric campaign. The inference of $\eta_{\text{EFAC}}$ is described in Section 4. We reiterate that the target image frames have been determined by the positions assumed for our reference

<table><caption>**Table 2.** An example set of astrometric results for J1738+0333, where the presented uncertainty excludes the calibrator reference point uncertainty as described in the text.</caption><thead><tr><th>obs. date<br>(yr)</th><th>α<sub>J2000</sub> (RA.)</th><th>δ<sub>J2000</sub> (Decl.)</th></tr></thead><tbody><tr><td>2015.6166</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>969242(3|5)</td><td>03°33'10''90430(9|17)</td></tr><tr><td>2015.8106</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>969329(3|6)</td><td>03°33'10''90491(9|18)</td></tr><tr><td>2016.6939</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>969726(5|6)</td><td>03°33'10''90981(16|21)</td></tr><tr><td>2017.1304</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970000(6|7)</td><td>03°33'10''91262(21|25)</td></tr><tr><td>2017.2068</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970040(2|4)</td><td>03°33'10''91217(7|15)</td></tr><tr><td>2017.2860</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970078(3|5)</td><td>03°33'10''91307(11|17)</td></tr><tr><td>2017.2997</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970062(17|17)</td><td>03°33'10''91272(59|61)</td></tr><tr><td>2017.7232</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970208(15|16)</td><td>03°33'10''91484(64|74)</td></tr><tr><td>2017.7669</td><td>17<sup>h</sup>38<sup>m</sup>53<sup>s</sup>970248(7|8)</td><td>03°33'10''91466(27|33)</td></tr></tbody></table>

Notes. • This table is compiled for PSR J1738+0333.

• The values on the left and the right side of ‘|’ are, respectively, statistical errors given in J1738+0333.pmpar.in.preliminary<fnref n="5" />, and systematics-included errors provided in J1738+0333.pmpar.in<fnref n="5" />.

sources (or virtual calibrators, see Section 3.1), and that any change in the assumed reference source position would transfer directly into a change in the recovered position for the target pulsar. Accordingly, the uncertainty in the reference source position must be accounted for in the pulsar’s reference position error budget, after fitting the pulsar’s astrometric parameters.

All pulsar positions and their error budgets are provided in the online<fnref n="5" /> ‘pmpar.in.preliminary’ and ‘pmpar.in’ files. The only difference between ‘pmpar.in.preliminary’ and ‘pmpar.in’ (for each pulsar) files are the position uncertainties: ‘pmpar.in.preliminary’ and ‘pmpar.in’ offer, respectively, position uncertainties $\sigma_i(0) = \sigma_i^R$ and $\sigma_i(1) = \sqrt{(\sigma_i^R)^2 + (\sigma_i^S)^2}$. As an example, the pulsar positions for PSR J1738+0333 are presented in Table 2, where the values on the left and right side of the ‘|’ sign stand for, respectively, $\sigma_i(0)$ and $\sigma_i(1)$. Additionally, to facilitate reproducibility, the image models for all primary and secondary phase calibrators listed in Table 1 are released<fnref n="6" /> along with this paper. Following Deller et al. (2019), Ding et al. (2020b), the calibrator models were made with the calibrator data concatenated from all epochs in an iterative manner.

### 3.1 1D interpolation on PSR J0621+1002 and PSR J1824–2452A

One can substantially reduce propagation-related systematic errors using 1D interpolation with two phase calibrators quasi-colinear with a target (e.g. Fomalont & Kopeikin 2003; Ding et al. 2020a). After 1D interpolation is applied, the target should in effect be referenced to a ‘virtual calibrator’ much closer (on the sky) than either of the two physical phase calibrators, assuming the phase screen can be approximated by a linear gradient with sky position (Ding et al. 2020a).

According to Table 1, seven secondary phase calibrators (or the final reference sources) are more than 20’ away from their targets, which would generally lead to relatively large systematic errors (e.g. Chatterjee et al. 2004; Kirsten et al. 2015; Deller et al. 2019). Fortunately, there are three MSPs – PSR J0621+1002, PSR J1824–2452A, and PSR J1939+2134, for which the pulsar and its primary and secondary phase calibrators are near-colinear (see online<fnref n="5" /> calibrator plans as well as Fig. 1). Hence, by applying 1D

<footnote n="4">⁴available at https://github.com/dingswin/psrvlbireduce</footnote>

<footnote n="5">⁵https://github.com/dingswin/publication_related_materials</footnote>

<footnote n="6">⁶https://github.com/dingswin/calibrator_models_for_astrometry</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Left: The calibrator plan for VLBI astrometry of PSR J1939+2134 (see Table 1 for full source names), where PSR J1939+2134 serves as the secondary phase calibrator and J1935 is the primary phase calibrator. Right: The zoomed-in field for reference sources as well as the virtual calibrator (VC) along the J1935-to-pulsar line. For the inverse 1D interpolation on PSR J1939+2134, we used the VC location that forms the largest included angle (65.7) with the two reference sources (see Section 3.2 for explanation), which corresponds to $\Delta_{\text{VC}} - \text{PC}/\Delta_{\text{PC}} - \text{psr} = 1.2836$ (i.e. the VC-to-J1935 separation is 1.2836 times the J1935-to-pulsar separation).</figcaption></figure>

interpolation, each of the 3 ‘1D-interpolation-capable’ MSPs can be referenced to a virtual calibrator much closer than the physical secondary phase calibrator (see Table 1).

We implemented 1D interpolation on PSR J0621+1002 and PSR J1824-2452A in the same way as the astrometry of the radio magnetar XTE J1810-197 carried out at 5.7 GHz (Ding et al. 2020a). None the less, due to our different observing frequency (i.e. 1.55 GHz), we estimated $\sigma_i^S$ differently. The post-1D-interpolation systematic errors should consist of (1) first-order residual systematic errors related to the target-to-virtual-calibrator offset $\Delta_{\text{psr-VC}}$ and (2) higher-order terms. Assuming negligible higher-order terms, we approached post-1D-interpolation $\sigma_i^S$ with equation (1) of Deller et al. (2019), using $\Delta_{\text{psr-VC}}$ as the calibrator-to-target separation. The assumption of negligible higher-order terms will be tested later and discussed in Section 4.1.3.

## 3.2 Inverse-referenced 1D interpolation on PSR J1939+2134

For PSR J1939+2134, normal 1D interpolation (Fomalont & Kopeikin 2003; Ding et al. 2020a) with respect to the primary phase calibrator ICRF J193510.4+203154 (J1935), and the brightest secondary reference source NVSS J194104+214913 (J194104) is still not the optimal calibration strategy. The ≈10-mJy (at 1.55 GHz) PSR J1939+2134 is the brightest MSP in the northern hemisphere, and only second to PSR J0437-4715 in the whole sky. After pulsar gating, PSR J1939+2134 is actually brighter than J194104. PSR J1939+2134 is unresolved on VLBI scales, and does not show long-term radio feature variations (frequently seen in quasars), which makes it an ideal secondary phase calibrator. Both factors encouraged us to implement the inverse-referenced 1D interpolation (or simply inverse 1D interpolation) on PSR J1939+2134, where PSR J1939+2134 is the de-facto secondary phase calibrator and the two ‘secondary phase calibrators’ serve as the targets. To avoid confusion, we refer to the two ‘secondary phase calibrators’ for PSR J1939+2134 (see Table 1) as secondary reference sources or simply reference sources.

Though inverse phase referencing (without interpolation) has been an observing/calibration strategy broadly used in VLBI astrometry (e.g. Imai et al. 2012; Yang et al. 2016; Li et al. 2018; Deller et al. 2019), inverse interpolation is new with the 2D approach of Hyland et al. (2022) at 8.3 GHz being a recent and independent development. We implemented inverse 1D interpolation at 1.55 GHz on PSR J1939+2134 in three steps (in addition to the standard procedure) detailed as follows.

### 3.2.1 Tying PSR J1939+2134 to the primary-calibrator reference frame

Inverse 1D interpolation relies on the residual phase solutions $\Delta\phi_n(\vec{x}, t)$ of self-calibration on PSR J1939+2134 (where $\vec{x}, t$, and $n$ refers to, respectively, sky position, time, and the $n$-th station in a VLBI array), which, however, change with $\Delta\vec{x}_{\text{psr}}$ – the displacement from the ‘true’ pulsar position to its model position. When $|\Delta\vec{x}_{\text{psr}}|$ is much smaller than the synthesized beam size $\theta_{\text{syn}}$, the changes in $\Delta\phi_n(\vec{x}, t)$ would be equal across all epochs, hence not biasing the resultant parallax and proper motion. However, if $|\Delta\vec{x}_{\text{psr}}| \gtrsim \theta_{\text{syn}}$, then the phase wraps of $\Delta\phi_n(\vec{x}, t)$ would likely become hard to uncover. The main contributor of considerable $|\Delta\vec{x}_{\text{psr}}|$ is an inaccurate pulsar position. The proper motion of the pulsar would also increase $|\Delta\vec{x}_{\text{psr}}|$ with time, if it is poorly constrained (or neglected). For PSR J1939+2134, the effect of proper motion across our observing duration is small ($\lesssim 1$ mas across the nominal observing span of 2.5 years; see the timing proper motion in Section 5) compared to $\theta_{\text{syn}} \sim 10$ mas.

In order to minimize $|\Delta\vec{x}_{\text{psr}}|$, we shifted the pulsar model position on an epoch-to-epoch basis by $\Delta\vec{x}_{\text{cor}}$ (which ideally should approximate $-\Delta\vec{x}_{\text{psr}}$) to the position measured in the J1935 reference frame (see Section 4.1 of Ding et al. 2020a for explanation of ‘reference frame’). This J1935-frame position was derived with the method for determining pulsar absolute position (Ding et al. 2020b) (where J194104 was used temporarily as the secondary phase calibrator) except that there is no need to quantify the position uncertainty.
We typically found $|\Delta\vec{x}_{psr}| \sim 50$ mas, which is well above $\theta_{syn} \sim 10$ mas. After the map centre shift, PSR J1939+2134 becomes tied to the J1935 frame.

### 3.2.2 1D interpolation on the tied PSR J1939+2134

The second step of inverse 1D interpolation is simply the normal 1D interpolation on PSR J1939+2134 that has been tied to the J1935 frame as described above (in Section 3.2.1). When there is only one secondary reference source, optimal 1D interpolation should see the virtual calibrator moved along the interpolation line (that passes through both J1935 and PSR J1939+2134) to the closest position to the secondary reference source (e.g. Ding et al. 2020a). However, there are two reference sources for PSR J1939+2134 (see Table 1), and the virtual calibrator point can be set at a point that will enable both of them to be used.

After calibration, a separate position series can be produced for each reference source. While we used each reference-source position series to infer astrometric parameters separately, we can also directly infer astrometric parameters with the combined knowledge of the two position series (which can be realized with *sterne*<fnref n="7" />). If the errors in the two position series are (largely) uncorrelated, this can provide superior astrometric precision. Since position residuals should be spatially correlated, we would ideally set the virtual calibrator at a location such that the included angle between the two reference sources is 90°. While achieving this ideal is not possible, we chose a virtual calibrator location that forms the largest possible included angle (65°7) with the two reference sources to minimise spatially correlated errors (see Fig. 1). This virtual calibrator is 1.2836 times further away from J1935 than PSR J1939+2134. Accordingly, the $\Delta\phi_{ij}(\vec{x}, t)$ solutions (obtained from the self-calibration on the tied PSR J1939+2134) were multiplied by 1.2836, and applied to the two reference sources.

### 3.2.3 De-shifting reference source positions

After data reduction involving the two steps outlined in Sections 3.2.1 and 3.2.2, one position series was acquired for each reference source. At this point, however, the two position series are not yet ready for astrometric inference, mainly because both proper motion and parallax signatures have been removed in the first step (see Section 3.2.1) when PSR J1939+2134 was shifted to its J1935-frame position. Therefore, the third step of inverse 1D interpolation is to cancel out the PSR J1939+2134 shift (made in the first step) by moving reference source positions by $-1.2836 \cdot \Delta\vec{x}_{cor}$, where the multiplication can be understood by considering fig. 1 of Ding et al. 2020a. This de-shifting operation was carried out separately outside the data reduction pipeline⁴. After the operation, we estimated $\sigma_{ij}^S$ of the reference sources (where $j = 1, 2$ refers to an individual reference source) following the method described in Section 3.1. The final position series of the reference sources are available online⁵. The astrometric parameter inference based on these position series is outlined in Section 4.

## 4 ASTROMETRIC INference METHODS AND QUASI-VLBI-ONLY ASTROMETRIC RESULTS

After gathering the position series⁵ with basic uncertainty estimation (see Section 3), we proceed to infer the astrometric parameters.

The inference is made by three different methods: (a) direct fitting of the position series with pmpar,<fnref n="8" /> (b) bootstrapping (see Ding et al. 2020b), and (c) Bayesian analysis using *sterne*<fnref n="7" /> (see Ding et al. 2021a). The two former methods directly adopt $\sigma_i(1) = \sqrt{(\sigma_i^R)^2 + (\sigma_i^S)^2}$ as the position errors. In Bayesian analysis, however, we inferred $\eta_{EFAC}$ along with other model parameters using the likelihood terms

$$P_1 \propto \left( \prod_i \sigma_i \right)^{-1} \exp \left[ -\frac{1}{2} \sum_i \left( \frac{\Delta \epsilon_i}{\sigma_i} \right)^2 \right], \quad (2)$$

where $\sigma_i = \sigma_i(\eta_{EFAC})$ obeys equation (1); $\Delta\epsilon_i$ refers to the model offsets from the measured positions. As is discussed in Section 4.2, Bayesian inference outperforms the other two methods, and is hence consistently used to present final results in this work. In all cases, the uncertainty in the reference source position should be added in quadrature to the uncertainty in the pulsar's reference position acquired with any method (of the three), in order to obtain a final estimate of the absolute positional uncertainty of the pulsar.

To serve different scientific purposes, we present two sets of astrometric results in two sections (i.e. Sections 4 and 5), which differ in whether timing proper motions and parallaxes are used as prior information in the inference.

### 4.0.4 Priors of canonical model parameters used in Bayesian analysis

To facilitate reproduction of our Bayesian results, the priors (of Bayesian inference) we use for canonical model parameters and $\eta_{EFAC}$ are detailed as follows. Priors for the two orbital parameters can be found in Section 4.3. We universally adopt the prior uniform distribution $U(0, 15)$ (i.e. uniformly distributed between 0 and 15) for $\eta_{EFAC}$. This prior distribution can be refined for future work with an ensemble of results across many pulsars. With regard to the canonical astrometric parameters (seven parameters for PSR J1939+2134 and five for the other pulsars), we adopt $U(X_0^{(\text{DF})} - 20, \bar{\sigma}_X^{(\text{DF})}, X_0^{(\text{DF})} + 20, \bar{\sigma}_X^{(\text{DF})})$ for each X, where X refers to one of $\alpha_{\text{ref}}, \delta_{\text{ref}}, \mu_{\alpha}, \mu_{\delta}$, and $\bar{\sigma}$. Here, $X_0^{(\text{DF})}$ stands for the direct-fitting estimate of X; $\bar{\sigma}_X^{(\text{DF})}$ represents the direct-fitting error corrected by the reduced chi-square $\chi_v^2$ (see Table 3) with $\bar{\sigma}_X^{(\text{DF})} = \sigma_X^{(\text{DF})} \cdot \sqrt{\chi_v^2}$. The calculation of prior range of X is made with the function *sterne.priors(generate_initsfile*<fnref n="7" />. We note that the adopted priors are relaxed enough to ensure robust outcomes: shrinking or enlarging the prior ranges by a factor of two would not change the inferred values. Meanwhile, the specified prior ranges are also reasonably small so that the global minimum of equation (2) can be reached.

## 4.1 Astrometric inference disregarding orbital motion

### 4.1.1 Single-reference-source astrometric inferences

All MSPs (in this work) excepting PSR J1939+2134 have only one reference source. For each of these single-reference-source MSPs, we fit for the five canonical astrometric parameters, i.e. reference position ($\alpha_{\text{ref}}$ and $\delta_{\text{ref}}$), proper motion ($\mu_{\alpha} \equiv \dot{\alpha} \cos \delta$ and $\mu_{\delta}$), and parallax ($\bar{\sigma}$). In the Bayesian analysis alone, $\eta_{EFAC}$ is also inferred alongside the astrometric parameters. At this stage, we neglect any orbital reflex motion for binary pulsars – the effects of orbital reflex motion are addressed in Section 4.3. The proper motions

<footnote n="7">⁷https://github.com/dingswin/sterne</footnote>

<footnote n="8">⁸https://github.com/walterfb/pmpar</footnote>
<table><caption>**Table 3.** Proper motion ($\mu_{\alpha}$, $\mu_{\delta}$) and parallax $\sigma$ from astrometry inferences disregarding orbital motion.</caption><thead><tr><th>PSR</th><th>μ<sub>α</sub><sup>(DF)</sup><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><sup>(DF)</sup><br>(mas yr<sup>-1</sup>)</th><th>σ<sup>(DF)</sup><br>(mas)</th><th>χ<sub>ν</sub><sup>2</sup></th><th>μ<sub>δ</sub><sup>(Bo)</sup><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><sup>(Bo)</sup><br>(mas yr<sup>-1</sup>)</th><th>σ<sup>(Bo)</sup><br>(mas)</th><th>μ<sub>α</sub><sup>(Ba)</sup><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>α</sub><sup>(Ba)</sup><br>(mas yr<sup>-1</sup>)</th><th>σ<sup>(Ba)</sup><br>(mas)</th><th>η<sub>EF&amp;C</sub></th><th>η<sub>orb</sub></th><th>P<sub>b</sub><br>(d)</th></tr></thead><tbody><tr><td colspan="15">Non-1D-interpolated results</td></tr><tr><td>J0030+0451</td><td>-6.13(4)</td><td>0.33(9)</td><td>3.02(4)</td><td>1.4</td><td>-6.1(1)</td><td>0.34<sup>+0.09</sup><sub>-0.08</sub></td><td>3.02<sup>+0.09</sup><sub>-0.08</sub></td><td></td><td>-6.13(7)</td><td>0.34<sup>+0.15</sup><sub>-0.16</sub></td><td>3.02(7)</td><td>1.35<sup>+0.45</sup><sub>-0.52</sub></td><td>-</td><td>-</td></tr><tr><td>J0610-2100</td><td>9.11(8)</td><td>15.9(2)</td><td>0.74(7)</td><td>1.0</td><td>9.10<sup>+0.06</sup><sub>-0.05</sub></td><td>15.9(2)</td><td>0.73<sup>+0.05</sup><sub>-0.04</sub></td><td></td><td>9.1(1)</td><td>15.96<sup>+0.25</sup><sub>-0.24</sub></td><td>0.73(10)</td><td>1.1<sup>+0.4</sup><sub>-0.2</sub></td><td>3 × 10<sup>-3</sup></td><td>0.29</td></tr><tr><td>J0621+1002</td><td>3.51(9)</td><td>-1.32(16)</td><td>0.88(7)</td><td>3.7</td><td>3.5<sup>+0.3</sup><sub>-0.4</sub></td><td>-1.3(2)</td><td>0.9<sup>+0.2</sup><sub>-0.3</sub></td><td></td><td>3.4(2)</td><td>-1.3(5)</td><td>0.85(21)</td><td>2.33<sup>+0.56</sup><sub>-0.44</sub></td><td>0.4</td><td>8.3</td></tr><tr><td>J1012+5307*</td><td>2.68(3)</td><td>-25.38(6)</td><td>1.17(2)</td><td>1.9</td><td>2.67<sup>+0.13</sup><sub>-0.06</sub></td><td>-25.39(12)</td><td>1.18<sup>+0.05</sup><sub>-0.06</sub></td><td></td><td>2.67(5)</td><td>-25.39<sup>+0.14</sup><sub>-0.15</sub></td><td>1.17<sup>+0.04</sup><sub>-0.05</sub></td><td>1.7<sup>+0.6</sup><sub>-0.4</sub></td><td>0.1</td><td>0.60</td></tr><tr><td>J1024-0719</td><td>-35.32(4)</td><td>-48.2(1)</td><td>0.94(3)</td><td>1.1</td><td>-48.2(1)</td><td>-48.2(1)</td><td>0.94<sup>+0.07</sup><sub>-0.06</sub></td><td></td><td>-35.32(7)</td><td>-48.1(2)</td><td>0.94(6)</td><td>1.2<sup>+0.4</sup><sub>-0.3</sub></td><td>-</td><td>-</td></tr><tr><td>J1518+4904</td><td>-0.69(2)</td><td>-8.54(4)</td><td>1.24(2)</td><td>1.2</td><td>-35.32<sup>+0.05</sup><sub>-0.04</sub></td><td></td><td>1.25(3)</td><td></td><td>-0.69(3)</td><td>-8.52(7)</td><td>1.24(3)</td><td>1.2<sup>+0.4</sup><sub>-0.3</sub></td><td>4.9</td><td>8.6</td></tr><tr><td>J1537+1155*</td><td>1.51(2)</td><td>-25.31(5)</td><td>1.06(7)</td><td>0.81</td><td>-0.69(3)</td><td>-8.53<sup>+0.10</sup><sub>-0.07</sub></td><td></td><td>1.07<sup>+0.09</sup><sub>-0.08</sub></td><td></td><td>1.51(3)</td><td>-25.30<sup>+0.05</sup><sub>-0.06</sub></td><td>1.06<sup>+0.11</sup><sub>-0.10</sub></td><td>0.54<sup>+0.57</sup><sub>-0.38</sub></td><td>0.24</td><td>0.42</td></tr><tr><td>J1640+2224*</td><td>2.199(56)</td><td>-11.29(9)</td><td>0.676(46)</td><td>1.1</td><td>-0.69(3)</td><td>-25.31<sup>+0.04</sup><sub>-0.05</sub></td><td></td><td>0.68<sup>+0.07</sup><sub>-0.06</sub></td><td></td><td>2.20(9)</td><td>-11.29<sup>+0.16</sup><sub>-0.14</sub></td><td>0.68(7)</td><td>1.29<sup>+0.66</sup><sub>-0.54</sub></td><td>3.1</td><td>175</td></tr><tr><td>J1643-1224</td><td>6.2(2)</td><td>3.3(5)</td><td>1.3(1)</td><td>0.8</td><td>6.1<sup>+0.7</sup><sub>-0.1</sub></td><td>3.3(4)</td><td>1.3<sup>+0.2</sup><sub>-0.5</sub></td><td></td><td>6.2(2)</td><td>3.3(6)</td><td>1.33(18)</td><td>1.0<sup>+0.3</sup><sub>-0.2</sub></td><td>1.1</td><td>147</td></tr><tr><td>J1721-2457</td><td>2.5(1)</td><td>-1.9(3)</td><td>0.02(7)</td><td>4.4</td><td>2.5(6)</td><td>-1.7<sup>+0.4</sup><sub>-0.7</sub></td><td></td><td>0.2<sup>+0.2</sup><sub>-0.6</sub></td><td></td><td>2.5(3)</td><td>-1.9(9)</td><td>0.0(2)</td><td>3.1<sup>+0.8</sup><sub>-0.3</sub></td><td>-</td><td>-</td></tr><tr><td>J1730-2304</td><td>20.3(1)</td><td>-4.79(26)</td><td>1.56(9)</td><td>1.7</td><td></td><td>-4.80<sup>+0.33</sup><sub>-0.35</sub></td><td></td><td>1.56<sup>+0.15</sup><sub>-0.17</sub></td><td></td><td>20.3(2)</td><td>-4.8(5)</td><td>1.57(18)</td><td>1.4<sup>+0.3</sup><sub>-0.2</sub></td><td>-</td><td>-</td></tr><tr><td>J1738+0333</td><td>6.97(4)</td><td>5.18(7)</td><td>0.51(3)</td><td>1.8</td><td>7.00<sup>+0.06</sup><sub>-0.11</sub></td><td>5.2(1)</td><td>0.50<sup>+0.07</sup><sub>-0.06</sub></td><td></td><td>6.98(8)</td><td>5.18(16)</td><td>0.50(6)</td><td>1.9<sup>+0.7</sup><sub>-0.6</sub></td><td>0.02</td><td>0.35</td></tr><tr><td>J1824-2452A</td><td>-0.03(49)</td><td>-6.6(1.2)</td><td>-0.25(37)</td><td>0.8</td><td>-0.03<sup>+0.28</sup><sub>-0.03</sub></td><td>-6.8<sup>+1.3</sup><sub>-1.5</sub></td><td>-0.20<sup>+0.54</sup><sub>-0.09</sub></td><td></td><td>-0.1(6)</td><td>-6.6(1.5)</td><td>-0.22(48)</td><td>0.9<sup>+0.4</sup><sub>-0.3</sub></td><td>-</td><td>-</td></tr><tr><td>J1853+1303</td><td>-1.37(9)</td><td>-2.8(2)</td><td>0.49(6)</td><td>1.0</td><td>-1.39<sup>+0.28</sup><sub>-0.28</sub></td><td>-2.8(2)</td><td>0.48<sup>+0.13</sup><sub>-0.13</sub></td><td></td><td>-1.36(10)</td><td>-2.8(2)</td><td>0.50(6)</td><td>0.45<sup>+0.36</sup><sub>-0.30</sub></td><td>1.4</td><td>116</td></tr><tr><td>J1910+1256</td><td>0.49(8)</td><td>-6.85(15)</td><td>0.26(7)</td><td>0.15</td><td>0.49<sup>+0.07</sup><sub>-0.08</sub></td><td>-6.85(6)</td><td>0.26<sup>+0.06</sup><sub>-0.02</sub></td><td></td><td>0.50(4)</td><td>-6.85(9)</td><td>0.25(3)</td><td>0.19<sup>+0.15</sup><sub>-0.12</sub></td><td>0.87</td><td>58.5</td></tr><tr><td>J1911-1114</td><td>-13.8(1)</td><td>-10.3(2)</td><td>0.38(9)</td><td>0.9</td><td></td><td>-10.3<sup>+1.1</sup><sub>-0.3</sub></td><td>0.36(8)</td><td></td><td>-13.8(2)</td><td>-10.3(4)</td><td>0.38<sup>+0.13</sup><sub>-0.14</sub></td><td>1.1<sup>+0.4</sup><sub>-0.3</sub></td><td>3 × 10<sup>-4</sup></td><td>2.7</td></tr><tr><td>J1918-0642</td><td>-7.12(8)</td><td>-5.7(2)</td><td>0.60(7)</td><td>1.1</td><td>-13.76<sup>+0.06</sup><sub>-0.08</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>J1939+2134<sup>(i)</sup></td><td>0.07(14)</td><td>-0.24(24)</td><td>0.35(10)</td><td>0.2</td><td>-7.1(1)</td><td>-5.8<sup>+0.3</sup><sub>-0.2</sub></td><td>0.60<sup>+0.12</sup><sub>-0.08</sub></td><td></td><td>-7.1(1)</td><td>-5.7(3)</td><td>0.60(12)</td><td>1.2<sup>+0.3</sup><sub>-0.2</sub></td><td>0.27</td><td>10.9</td></tr><tr><td colspan="15">Single-reference-source 1D-interpolated results</td></tr><tr><td>J0621+1002</td><td>3.68(6)</td><td>-1.33(9)</td><td>0.94(4)</td><td>6.3</td><td>3.53<sup>+0.36</sup><sub>-0.35</sub></td><td>-1.3(2)</td><td>0.9<sup>+0.2</sup><sub>-0.3</sub></td><td></td><td>3.5(2)</td><td>-1.37(35)</td><td>0.86(15)</td><td>3.7<sup>+1.0</sup><sub>-0.7</sub></td><td>0.4</td><td>8.3</td></tr><tr><td>J1824-2452A</td><td>0.3(3)</td><td>-3.7(6)</td><td>0.1(3)</td><td>1.5</td><td>0.4<sup>+0.5</sup><sub>-0.5</sub></td><td>-4<sup>+1</sup><sub>-2</sub></td><td>0.1<sup>+0.7</sup><sub>-0.4</sub></td><td></td><td>0.3(6)</td><td>-3.9<sup>+1.2</sup><sub>-1.3</sub></td><td>0.1(5)</td><td>1.6<sup>+0.8</sup><sub>-0.7</sub></td><td>-</td><td>-</td></tr><tr><td>J1939+2134<sup>(ii)</sup></td><td>0.08(4)</td><td>-0.45(6)</td><td>0.38(3)</td><td>1.3</td><td>0.08<sup>+0.07</sup><sub>-0.07</sub></td><td>-0.44<sup>+0.09</sup><sub>-0.02</sub></td><td>0.36<sup>+0.05</sup><sub>-0.05</sub></td><td></td><td>0.08(7)</td><td>-0.44(11)</td><td>0.380<sup>+0.049</sup><sub>-0.049</sub></td><td>1.4<sup>+0.7</sup><sub>-0.2</sub></td><td>-</td><td>-</td></tr><tr><td>J1939+2134<sup>(iii)</sup></td><td>0.3(3)</td><td>-0.3(2)</td><td>0.36(19)</td><td>1.3</td><td>0.2<sup>+1.0</sup><sub>-0.2</sub></td><td>-0.3(4)</td><td>0.36<sup>+0.35</sup><sub>-0.58</sub></td><td></td><td>0.3(4)</td><td>-0.3(4)</td><td>0.38<sup>+0.29</sup><sub>-0.28</sub></td><td>3.7<sup>+3.4</sup><sub>-2.6</sub></td><td>-</td><td>-</td></tr><tr><td colspan="15">Multireference-source 1D-interpolated results</td></tr><tr><td>J1939+2134<sup>(iv)</sup></td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td></td><td>0.08(7)</td><td>-0.43(11)</td><td>0.384<sup>+0.048</sup><sub>-0.046</sub></td><td>1.5<sup>+0.7</sup><sub>-0.6</sub></td><td>-</td><td>-</td></tr></tbody></table>

Notes. *DF', 'Bo', and 'Ba' stands for, respectively, direct fitting, bootstrap and Bayesian inference. $\chi^2$ is the reduced chi-square of direct fitting using pmpars⁸.

*The middle and top block presents, respectively, 1D-interpolated (see Sections 3.1 and 3.2) and non-1D-interpolated results.

• The bottom entry for PSR J1939+2134 shows the result of multireference-source astrometry inference (see Section 4.1.2).

• $P_b$ represents orbital period (see Table 4 for their references). $\eta_{orb}$ is defined in equation (3).

* Already published in Vigeland et al. (2018), Ding et al. (2020b), Ding et al. (2021a).

(i) Based on (non-1D-interpolated) J194104 positions inverse-referenced to PSR J1939+2134.

(ii) Using 1D-interpolated J194104 positions inverse-referenced to PSR J1939+2134.

(iii) Using 1D-interpolated J194106 positions inverse-referenced to PSR J1939+2134.

(iv) Based on 1D-interpolated J194104 and J194106 positions inverse-referenced to PSR J1939+2134 (see Section 3.2).
<figure><img src="image_2.png" /><figcaption>**Figure 2.** The common parallax signature of PSR J1939+2134 revealed by the position measurements of both reference sources (see Table 1). In both panels, the best-fit proper motion has been subtracted. The magenta curve in each panel represents the best-inferred astrometric model. The fuzzy region around the curve consists of various Bayesian simulations, the scatter of which can visualize the uncertainty level of the underlying model parameters (see Section 4.1.2). As a result of the inverse referencing, the common parallax revealed here is actually the negative of the PSR J1939+2134 parallax presented in Table 3.</figcaption></figure>

and parallaxes derived with single-reference-source astrometry and disregarding orbital motion are summarized in Table 3. The reference positions are presented in Section 4.4.

### 4.1.2 Multisource astrometry inferences

When multiple sources share proper motion and/or parallax (while each source having its own reference position), a joint multisource astrometry inference can increase the degrees of freedom of inference (i.e. the number of measurements reduced by the number of parameters to infer), and tighten constraints on the astrometric parameters. Multisource astrometry inference has been widely used in maser astrometry (where maser spots with different proper motions scatter around a region of high-mass star formation, Reid et al. 2009), but has not yet been used for any pulsar, despite the availability of several bright pulsars with multiple in-beam calibrators (e.g. PSR J0332+5434, PSR J1136+1551) in the PSRπ project (Deller et al. 2019).

PSR J1939+2134 is the only source (in this work) that has multiple (i.e. two) reference sources, which provides a rare opportunity to test multireference-source astrometry. We assumed that the position series of J194104 is uncorrelated with that of NVSS J194106+215304 (hereafter J194106), and utilized sterne⁷ to infer the common parallax and proper motion, alongside two reference positions (one for each reference source). The acquired proper motion and parallax are listed in Table 3. As inverse phase referencing is applied for PSR J1939+2134, the parallax and proper motion of PSR J1939+2134 are the inverse of the direct astrometric measurements. For comparison, the proper motion and parallax inferred solely with one reference source are also reported in Table 3. Due to the relative faintness of J194106 (see Table 1), the inclusion of J194106 only marginally improves the astrometric results (e.g. $\sigma$) over those inferred with J194104 alone.

The constraints on the parallax (as well as the proper motion) are visualized in Fig. 2. The best-inferred model (derived from the

J194104 and J194106 positions) is illustrated with a bright magenta curve, amidst two sets of Bayesian simulations – each set for a reference source. Each simulated curve is a time series of simulated positions, with the best-inferred reference position ($\alpha_{\text{ref},j}$ and $\delta_{\text{ref},j}$, where j refers to either J194104 or J194106) and proper-motion-related displacements (i.e. $\mu_{\alpha}\Delta t$ and $\mu_{\delta}\Delta t$, where $\Delta t$ is the time delay from the reference epoch) subtracted. As the simulated curve depends on the underlying model parameters, the degree of scatter of simulated curves would increase with larger uncertainties of model parameters. Though sharing simulated parallaxes and proper motions with J194104, the simulated curves for J194106 exhibits broader scatter (than the J194104 ones) owing to more uncertain reference position (see Section 4.4 for $\alpha_{\text{ref},J194106}$ and $\delta_{\text{ref},J194106}$). The large scatter implies that the J194106 position measurements impose relatively limited constraints on the common model parameter (i.e. parallax and proper motion), which is consistent with the findings from Table 3.

### 4.1.3 Implications for 1D/2D interpolation

On the three 1D-interpolation-capable MSPs, we compared astrometric inference with both the 1D-interpolated and non-1D-interpolated position series (one at a time). For PSR J1939+2134, the $\eta_{\text{EFAC}}$ of the three 1D-interpolated realizations are consistent with each other, but larger than the non-1D-interpolated counterpart. This post-1D-interpolation inflation of $\eta_{\text{EFAC}}$ also occurs to the other two 1D-interpolation-capable pulsars (see Table 3), which suggests the post-1D-interpolation fiducial systematic errors $\sigma_i^S$ might be systematically under-estimated. One obvious explanation for this under-estimation is that the higher-order terms of systematic errors are non-negligible (as opposed to the assumption we started with in Section 3.1): they might be actually comparable to the first-order residual systematic errors (that are related to $\Delta_{\text{psr-VC}}$) at the ~1.55 GHz observing frequencies.

On the other hand, the astrometric results based on the non-1D-interpolated J194104 positions inverse-referenced to PSR J1939+2134 are less precise than the 1D-interpolated counterpart by ≈40 per cent, as is also the case for PSR J0621+1002 (see Table 3). Moreover, the post-1D-interpolation parallax of PSR J1824–2452A becomes relatively more accurate than the negative parallax obtained without applying 1D interpolation. All of these demonstrate the utility of 1D/2D interpolation, even in the scenario of in-beam astrometry that is already precise. In the remainder of this paper, we only focus on the 1D-interpolated astrometric results for the three 1D-interpolation-capable MSPs.

## 4.2 Bayesian inference as the major method for MSPSRπ

We now compare the three sets of astrometric parameters (in Table 3) obtained with different inference methods, and seek to proceed with only one set in order to simplify the structure of this paper. Among the three inference methods we use in this work, direct least-square fitting is the most time-efficient, but is also the least robust against improperly estimated positional uncertainties. Conversely, the other two methods (i.e. bootstrap and Bayesian methods) do not rely solely on the input positional uncertainties, and can still estimate the model parameters, and their uncertainties $\sigma_X^{(Y)}$ ($X = \mu_{\alpha}, \mu_{\delta}$ or $\sigma$; $Y = \text{'Bo'}$ or $\text{'Ba'}$) more robustly in the presence of incorrectly estimated positional errors.

Generally speaking, $\sigma_X^{(Y)}$ inferred from a pulsar position series are expected to change with the corresponding $\chi^2$-corrected direct-
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Bootstrap (denoted as ‘Bo’) and Bayesian (‘Ba’) errors (of three inferred parameters) divided by the corresponding $\chi_v^2$-corrected direct-fitting errors. Here, $\tilde{\sigma}_X^{(\text{DF})} = \sigma_X^{(\text{DF})} \cdot \sqrt{\chi_v^2}$ represents the $\chi_v^2$-corrected errors of direct fitting, where X stands for one of the $\mu_\alpha, \mu_\delta$, and $\omega$ groups. The dimensionless $\tilde{\sigma}_X^{(\text{DF})}$ is defined as an individual $\tilde{\sigma}_X^{(\text{DF})}$ divided by the standard deviation $s_X^{(\text{DF})}$ for all $\tilde{\sigma}_X^{(\text{DF})}$ of the group X. The grey and orange shaded regions show, respectively, the standard deviation of $\sigma_X^{(\text{Bo})}/\tilde{\sigma}_X^{(\text{DF})}$ and $\sigma_X^{(\text{Ba})}/\tilde{\sigma}_X^{(\text{DF})}$ across all of the three groups (i.e. $\mu_\alpha, \mu_\delta$, and $\omega$) around the respective mean value outlined with the grey and orange dashed lines. Both bootstrap and Bayesian errors are generally slightly higher than the level of direct-fitting errors illustrated with the cyan dashed line, and are well consistent with each other as anticipated. Despite the consistency, bootstrap errors show larger scatter than Bayesian ones.</figcaption></figure>

fitting error $\tilde{\sigma}_X^{(\text{DF})} \equiv \sigma_X^{(\text{DF})} \cdot \sqrt{\chi_v^2}$. In order to investigate the relation between $\sigma_X^{(\text{Y})}$ and $\tilde{\sigma}_X^{(\text{DF})}$, we divided $\sigma_X^{(\text{Y})}$ by $\tilde{\sigma}_X^{(\text{DF})}$ for each pulsar entry in the top block of Table 3. The results are displayed in Fig. 3. For the convenience of illustration, we calculated the dimensionless $\tilde{\sigma}_X^{(\text{DF})}$ defined as $\tilde{\sigma}_X^{(\text{DF})}/s_X^{(\text{DF})}$ (where $s_X^{(\text{DF})}$ represents the standard deviation of $\tilde{\sigma}_X^{(\text{DF})}$ over the group X), which allows all the three sets (i.e. $\mu_\alpha, \mu_\delta$, and $\omega$) of dimensionless $\tilde{\sigma}_X^{(\text{DF})}$ to be horizontally more evenly plotted in Fig. 3.

Across the entire MSPSRπ sample, we see that $\sigma_X^{(\text{Y})}$ scales with $\tilde{\sigma}_X^{(\text{DF})}$ in a near-linear fashion. The mean scaling factors across all of the three parameter groups (i.e. $\mu_\alpha, \mu_\delta$, and $\omega$) are $\langle \sigma_X^{(\text{Bo})}/\tilde{\sigma}_X^{(\text{DF})} \rangle = 1.67 \pm 0.85$ and $\langle \sigma_X^{(\text{Ba})}/\tilde{\sigma}_X^{(\text{DF})} \rangle = 1.49 \pm 0.24$ (see Fig. 3). The two mean scaling factors show that parameter uncertainties inferred using either a bootstrap or Bayesian approach will be slightly higher (and on average, consistent between the two approaches) than would be obtained utilizing direct-fitting (illustrated with the cyan dashed line in Fig. 3).

The more optimistic uncertainty predictions of $\tilde{\sigma}_X^{(\text{DF})}$ can be understood as resulting from two causes: first, it neglects both the finite width and the skewness of the $\chi^2$ distribution, and second, to achieve the expected $\chi^2$ it scales the total uncertainty contribution at each epoch, rather than the systematic uncertainty contribution alone. When (as is typical for pulsar observations) the S/N and hence statistical positional precision can vary substantially between observing epochs, this simplified approach preserves the relative weighting between epochs, whereas increasing the estimated systematic uncertainty contribution acts to equalise the weighting between epochs (by reducing the position precision more for epochs, where the pulsar was bright and the statistical precision high than for epochs, where the pulsar was faint and the statistical precision is already low).

While the consistency between $\langle \sigma_X^{(\text{Bo})}/\tilde{\sigma}_X^{(\text{DF})} \rangle$ and $\langle \sigma_X^{(\text{Ba})}/\tilde{\sigma}_X^{(\text{DF})} \rangle$ suggests that both approaches can overcome this shortcoming in the direct fitting method, $\sigma_X^{(\text{Bo})}/\tilde{\sigma}_X^{(\text{DF})}$ shows a much larger scatter (3.5 times) compared to $\sigma_X^{(\text{Ba})}/\tilde{\sigma}_X^{(\text{DF})}$ (see Fig. 3). To determine which approach best represents the true (and unknown) parameter uncertainties, it is instructive to consider the outliers in the bootstrap distribution results.

First, consider cases where the bootstrap results in a lower uncertainty than $\tilde{\sigma}_X^{(\text{DF})}$. For the reasons noted above, we expect $\tilde{\sigma}_X^{(\text{DF})}$ to yield the most optimistic final parameter uncertainty estimates, and yet the bootstrap returns a lower uncertainty than $\tilde{\sigma}_X^{(\text{DF})}$ in a number of cases. Second, the cases with the highest values of $\sigma_X^{(\text{Bo})}/\tilde{\sigma}_X^{(\text{DF})}$ reach $\gtrsim 3$ on a number of occasions, which imply an extremely large (or very non-Gaussian) systematic uncertainty contribution, which would lead (in those cases) to a surprisingly low-reduced $\chi^2$ for the best-fitting model. Given the frequency with which these outliers arise, we regard it likely that bootstrap approach mis-estimates parameter uncertainties at least occasionally, likely due to the small number of observations available.

Therefore, we consider the Bayesian method described in this paper as the preferred inference method for the MSPSRπ sample, and consistently use the Bayesian results in the following discussions. We note that as continued VLBI observing campaigns add more results, the systematic uncertainty estimation scheme applied to Bayesian inference can be further refined in the future.

## 4.3 Astrometric inference accounting for orbital motion

For some binary pulsars, VLBI astrometry can also refine parameters related to the binary orbit on top of the canonical astrometric parameters. The orbital inclination $i$ and the orbital ascending
node longitude $\Omega_{asc}$ have been previously constrained for a few nearby pulsars, such as PSR J1022+1001, PSR J2145-0750, and PSR J2222-0137 (Deller et al. 2013; Deller et al. 2016; Guo et al. 2021). To assess the feasibility of detecting orbital reflex motion with VLBI, we computed

$$ \eta_{\text{orb}} \equiv \frac{2a_1}{1 \text{ au}} \cdot \frac{\sigma}{\sigma_{\bar{\sigma}}} = 2a_1 \cdot \left( \frac{1 \text{ au}}{\sigma} \right)^{-1} \cdot \frac{1}{\sigma_{\bar{\sigma}}} = \frac{2a_1}{D} \cdot \frac{1}{\sigma_{\bar{\sigma}}}, \quad (3) $$

where $D$ and $a_1 \equiv a \sin i$ stands for, respectively, the distance (to the pulsar) and the orbital semi-major axis projected onto the sightline. On the other hand, $\tilde{\theta}_{\text{orb}} \equiv 2a/D$ reflects the apparent angular size of orbit. Provided the parallax uncertainty $\sigma_{\bar{\sigma}}$, $\tilde{\theta}_{\text{orb}}/\sigma_{\bar{\sigma}}$ quantifies the detectability of orbital parameters using VLBI astrometry. Hence,

$$ \frac{\tilde{\theta}_{\text{orb}}}{\sigma_{\bar{\sigma}}} \equiv \frac{2a}{D} \cdot \frac{1}{\sigma_{\bar{\sigma}}} \geq \eta_{\text{orb}}. \quad (4) $$

Since $i$ is usually unknown, the $\eta_{\text{orb}}$ defined in equation (3) serves as a lower limit for $\tilde{\theta}_{\text{orb}}/\sigma_{\bar{\sigma}}$, and is used in this work to find out pulsar systems with $i$ and $\Omega_{asc}$ potentially measurable with VLBI observations. In general, the orbital reflex motion should be negligible when $\eta_{\text{orb}} \ll 1$, easily measurable when $\eta_{\text{orb}} \gg 1$, and difficult to constrain (but non-negligible) when $\eta_{\text{orb}} \sim 1$. By way of comparison, Guo et al. (2021) were able to firmly constrain $\Omega_{asc}$ and $i$ for PSR J2222-0137 ($\eta_{\text{orb}} = 10.2$), while Deller et al. (2016) could place weak constraints for PSR J1022+1001 and PSR J2145-0750 ($\eta_{\text{orb}} = 3.2$ and 1.6, respectively).

Accordingly, in this work, we fit for orbital reflex motion if all the following conditions are met:

(i) $a_1$ is well determined with pulsar timing;

(ii) $\eta_{\text{orb}} > 1$;

(iii) the orbital period $P_b < 2$ yr, where 2 yr is the nominal time span of an MSPSR$\pi$ astrometric campaign.

For the calculation of $\eta_{\text{orb}}$, we simply use the direct-fitting parallax $\sigma^{(\text{DF})}$ for $\sigma$, and its $\chi_v^2$-corrected uncertainty $\sigma_{\bar{\sigma}}^{(\text{DF})} \cdot \sqrt{\chi_v^2}$ for $\sigma_{\bar{\sigma}}$ (see Table 3). We note that this choice of parallax and its uncertainty would generally lead to slightly larger $\eta_{\text{orb}}$ compared to using $\sigma^{(\text{Ba})}$ and $\sigma_{\bar{\sigma}}^{(\text{Ba})}$, according to Fig. 3 and the discussion in Section 4.2. Nevertheless, the choice (1) enables the comparison with $\eta_{\text{orb}}$ of the historically published pulsars (that do not have $\sigma^{(\text{Ba})}$ and $\sigma_{\bar{\sigma}}^{(\text{Ba})}$), (2) simplifies the procedure of analysis, (3) facilitates the reproduction of $\eta_{\text{orb}}$ by other researchers, and (4) is more conservative in the sense that more candidates with $\eta_{\text{orb}} > 1$ would be found. The calculated $\eta_{\text{orb}}$ as well as $P_b$ are summarized in Table 3. Among the 18 MSPSR$\pi$ pulsars, PSR J1518+4904, PSR J1640+2224, PSR J1643-1224, and PSR J1853+1303 meet our criteria (see Table 3), where PSR J1518+4904 is a DNS system, and the others are pulsar-WD binaries. Hereafter, the four pulsars are referred to as the ‘8P’ pulsars for the sake of brevity, as we would perform 8-parameter (i.e. the five canonical astrometric parameters and $\eta_{\text{EFAC}}$ plus $i$ and $\Omega_{asc}$) inference on them.

For the 8-parameter inference, prior probability distributions of the canonical parameters and $\eta_{\text{EFAC}}$ are described in Section 4.0.1. Both $i$ and $\Omega_{asc}$ are defined in the TEMPO2 (Edwards, Hobbs & Manchester 2006) convention. The prior probability distribution of $\Omega_{asc}$ follows $U(0, 360^{\circ})$. Sine distribution $S(0, 180^{\circ})$ is used for $i$ of the four 8P pulsars (i.e. the probability density $p(i) \propto \sin i$, $i \in [0, 180^{\circ}$]). Where available, tighter constraints are applied to $i$ in accordance with Table 4 (also see the descriptions in Section 8).

Moreover, extra prior constraints can be applied to $i$ and $\Omega_{asc}$ based on $\dot{a}_1$, the time derivative of $a_1$ (e.g. Nice, Splaver & Stairs 2001;

<table><caption>Table 4. Prior constraints on $i$ and $\Omega_{asc}$.</caption><thead><tr><th>PSR</th><th>dot;<sub>1</sub><br>(10<sup>-15</sup> It-s s<sup>-1</sup>)</th><th>dot;<sub>1</sub>/a<sub>1</sub><br>(10<sup>-15</sup> s<sup>-1</sup>)</th><th>i</th><th>&Omega;<sub>asc</sub><br>(deg)</th></tr></thead><tbody><tr><td>J1518+4904</td><td>-11(3)<sup>a1</sup></td><td>-0.55(15)</td><td>sin i &le; 0.73<sup>a2</sup></td><td>-</td></tr><tr><td>J1640+2224</td><td>12(1)<sup>b</sup></td><td>0.22(2)</td><td>sin i = 0.973(9)<sup>b</sup></td><td>-</td></tr><tr><td>J1643-1224</td><td>-49.7(7)<sup>c</sup></td><td>-1.98(3)</td><td>-</td><td>-</td></tr><tr><td>J1853+1303</td><td>14(2)<sup>b</sup></td><td>0.34(5)</td><td>85(14)<sup>c d</sup></td><td>-</td></tr></tbody></table>

Notes. $a_1$ Janssen et al. (2008); $a_2$ inferred from the non-detection of Shapiro delay effects.

$^b$Perera et al. (2019); $^c$Reardon et al. (2021).

$^d$based on Shapiro delay measurements (Faisal Alam et al. 2021).

Deller et al. 2016; Reardon et al. 2021). As $a_1 \equiv a \sin i$,

$$ \frac{\dot{a}_1}{a_1} = \frac{\dot{a}}{a} + \frac{\partial i}{\partial t} \cot i \approx \frac{\partial i}{\partial t} \cot i. \quad (5) $$

Here, the $\dot{a}/a$ term reflects the intrinsic variation of the semimajor axis $a$ due to GR effects (Peters 1964), which is however ~8 and ~5 orders smaller than $\dot{a}_1/a_1$ for the 8P WD-pulsar systems and the DNS system PSR J1518+4904, respectively (see Nice et al. 2001 for an analogy). Accordingly, the apparent $\dot{a}_1/a_1$ is predominantly caused by apparent $i$ change as a result of the sightline shift (Kopeikin 1996). When proper motion contributes predominantly to the sky position shift (as is the case for the 8P pulsars),

$$ \frac{\partial i}{\partial t} = \mu \sin(\theta_{\mu} - \Omega_{\text{asc}}), \quad (6) $$

where $\theta_{\mu}$ refers to the position angle (east of north) of the proper motion $\mu$ (Kopeikin 1996; Nice et al. 2001). We incorporated the $\dot{a}_1/a_1$ measurements (with equations 5 and 6) on top of other prior constraints, and inferred $i$, $\Omega_{asc}$, $\eta_{\text{EFAC}}$, and the canonical five astrometric parameters for the 8P pulsars with sterne$^7$, following similar approaches taken by Deller et al. (2016), Guo et al. (2021).

While we ultimately did not significantly constrain $i$ or $\Omega_{asc}$ for any pulsar, including their non-negligible reflex motion in the inference is still necessary for correctly inferring the uncertainties of the non-orbital model parameters. The non-orbital inferred parameters are provided in Section 4.4 below, along with all the non-8P pulsars. As we found minimal differences between the constraints obtained on orbital parameters with or without the adoption of priors based on pulsar timing, we defer the presentation of the posterior constraints on orbital inclinations and ascending node longitudes (of the 8P pulsars) to Section 5 in order to avoid repetition.

## 4.4 The quasi-VLBI-only astrometric results

To wrap up this section, we summarize in Table 5 the full (including $\alpha_{\text{ref}}$ and $\delta_{\text{ref}}$) final astrometric results obtained with no exterior prior proper motion or parallax constraints, which we simply refer to as quasi-VLBI-only astrometric results (we add ‘quasi’ because timing constraints on two orbital parameters, i.e. $i$ and $\dot{a}_1$, have already been used for the 8P pulsars). These quasi-VLBI-only results are mainly meant for independent checks of timing results (which would enable the frame connection mentioned in Section 1.2), or as priors for future timing analyses. For the most precise possible pulsar parallaxes and hence distances, we recommend the use of the ‘VLBI+timing’ results presented in Section 5.

The reference positions $\alpha_{\text{ref}}$ and $\delta_{\text{ref}}$ we provide in Table 5 are precisely measured, but only with respect to the assumed location of the in-beam calibrator source for each pulsar. In all cases, the uncertainties on the in-beam source locations (also shown in Table 5)
<table><caption>**Table 5.** Final astrometric models inferred without using timing proper motions as priors.</caption><thead><tr><th>PSR</th><th>t<sub>ref</sub><br>(MJD)</th><th>α<sub>ref</sub> (J2000)*</th><th>σ<sub>αref</sub><br>(mas)</th><th>δ<sub>ref</sub> (J2000)*</th><th>σ<sub>δref</sub><br>(mas)</th><th>μ<sub>α</sub><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>σ<sub>μ</sub><br>(mas)</th><th>η<sub>EFAC</sub></th><th>ρ<sub>μα,σ</sub></th><th>ρ<sub>μδ,σ</sub></th></tr></thead><tbody><tr><td>J030+0451</td><td>57849</td><td>00<sup>b</sup>30<sup>m</sup>27<sup>s</sup>42502</td><td>0.06[±0.3±0.2±0.8]</td><td>04°51′39.7159</td><td>0.2[±0.8±0.3±0.8]</td><td>-6.13(7)</td><td>0.34<sup>+0.15</sup><sub>-0.16</sub></td><td>3.02(7)</td><td>1.35<sup>+0.45</sup><sub>-0.32</sub></td><td>0.39</td><td>0.08</td></tr><tr><td>J0610-2100</td><td>57757</td><td>06<sup>b</sup>10<sup>m</sup>13<sup>s</sup>60053</td><td>0.08[±0±3.5±0.8]</td><td>-21°00′27.7923</td><td>0.2[±0±15.2±0.8]</td><td>9.1(1)</td><td>15.96<sup>+0.25</sup><sub>-0.24</sub></td><td>0.73(10)</td><td>1.1<sup>+0.4</sup><sub>-0.2</sub></td><td>0.51</td><td>0.05</td></tr><tr><td>J0621+1002</td><td>57685</td><td>06<sup>b</sup>21<sup>m</sup>22<sup>s</sup>11617</td><td>0.12[±2.1±0.5±0.8]</td><td>10°02′38.7261</td><td>0.3[±2.8±0.7±0.8]</td><td>3.5(2)</td><td>-1.37(35)</td><td>0.86(15)</td><td>3.7<sup>+0.7</sup><sub>-0.7</sub></td><td>0.60</td><td>-0.01</td></tr><tr><td>J1012+5307</td><td>57700</td><td>10<sup>b</sup>12<sup>m</sup>33<sup>s</sup>43991</td><td>0.04[±1.2±0.3±0.8]</td><td>53°07′02<sup>s</sup>1110</td><td>0.1[±2.8±0.3±0.8]</td><td>2.67(5)</td><td>-25.39<sup>+0.14</sup><sub>-0.15</sub></td><td>1.17<sup>+0.04</sup><sub>-0.05</sub></td><td>1.7<sup>+0.6</sup><sub>-0.4</sub></td><td>0.37</td><td>-0.06</td></tr><tr><td>J1024-0719</td><td>57797</td><td>10<sup>b</sup>24<sup>m</sup>38<sup>s</sup>65725</td><td>0.06[±0.6±0.8±0.8]</td><td>-07°19′19′8014</td><td>0.2[±1.4±1.5±0.8]</td><td>-35.32(7)</td><td>-48.1(2)</td><td>0.94(6)</td><td>1.2<sup>+0.3</sup><sub>-0.4</sub></td><td>0.23</td><td>2 × 10<sup>-3</sup></td></tr><tr><td>J1518+4904<sup>b</sup></td><td>57795</td><td>15<sup>b</sup>18<sup>m</sup>1679817</td><td>0.04[±0.5±0.4±0.8]</td><td>49°04′34.1132</td><td>0.1[±1.4±0.4±0.8]</td><td>-0.69<sup>+0.04</sup><sub>-0.03</sub></td><td>-8.53<sup>+0.07</sup><sub>-0.09</sub></td><td>1.238(36)</td><td>1.4<sup>+0.6</sup><sub>-0.4</sub></td><td>-0.19</td><td>-0.68</td></tr><tr><td>J1537+1155</td><td>57964</td><td>15<sup>b</sup>37<sup>m</sup>09<sup>s</sup>96347</td><td>0.06[±1.2±0.1±0.8]</td><td>11°55′55<sup>s</sup>0274</td><td>0.1[±2.7±0.1±0.8]</td><td>1.51(3)</td><td>-25.30<sup>+0.05</sup><sub>-0.06</sub></td><td>1.06<sup>+0.11</sup><sub>-0.10</sub></td><td>0.5<sup>+0.6</sup><sub>-0.4</sub></td><td>-0.17</td><td>0.04</td></tr><tr><td>J1640+2224<sup>b</sup></td><td>57500</td><td>16<sup>b</sup>40<sup>m</sup>1674587</td><td>0.07[±0.4±0.3±0.8]</td><td>22°24′08<sup>s</sup>7642</td><td>0.1[±0.3±0.6±0.8]</td><td>2.19(9)</td><td>-11.30<sup>+0.13</sup><sub>-0.16</sub></td><td>0.68(8)</td><td>1.3<sup>+0.7</sup><sub>-0.4</sub></td><td>-0.57</td><td>-0.02</td></tr><tr><td>J1643-1224<sup>b</sup></td><td>57700</td><td>16<sup>b</sup>43<sup>m</sup>38<sup>s</sup>16407</td><td>0.1[±2.0±0.1±0.8]</td><td>-12°24′58<sup>s</sup>6531</td><td>0.4[±6.2±0.1±0.8]</td><td>6.2(2)</td><td>3.3(6)</td><td>1.31<sup>+0.17</sup><sub>-0.18</sub></td><td>1.0<sup>+0.5</sup><sub>-0.4</sub></td><td>-0.43</td><td>0.04</td></tr><tr><td>J1721-2457</td><td>57820</td><td>17<sup>b</sup>21<sup>m</sup>05<sup>s</sup>49936</td><td>0.2[±3.0±0.3±0.8]</td><td>-24°57′06<sup>s</sup>2210</td><td>0.6[±7.3±0.6±0.8]</td><td>2.5(3)</td><td>-1.9(9)</td><td>0.0(2)</td><td>3.1<sup>+0.8</sup><sub>-0.2</sub></td><td>-0.13</td><td>-0.01</td></tr><tr><td>J1730-2304</td><td>57821</td><td>17<sup>b</sup>30<sup>m</sup>21<sup>s</sup>67969</td><td>0.2[±0.5±0.3±0.8]</td><td>-23°04′31<sup>s</sup>1749</td><td>0.5[±1.1±0.6±0.8]</td><td>20.3(2)</td><td>-4.8(5)</td><td>1.57(18)</td><td>1.4<sup>+0.3</sup><sub>-0.7</sub></td><td>-0.38</td><td>0.01</td></tr><tr><td>J1738+0333</td><td>57829</td><td>17<sup>b</sup>38<sup>m</sup>5397001</td><td>0.06[±0.2±0.3±0.8]</td><td>03°33′10<sup>s</sup>9124</td><td>0.1[±0.6±0.6±0.8]</td><td>6.98(8)</td><td>5.18(16)</td><td>0.50(6)</td><td>1.9<sup>+0.6</sup><sub>-0.7</sub></td><td>-0.53</td><td>4 × 10<sup>-4</sup></td></tr><tr><td>J1824-2452A</td><td>57836</td><td>18<sup>b</sup>24<sup>m</sup>32<sup>s</sup>00791</td><td>0.4[±2.0±0.1±0.8]</td><td>-24°52′10<sup>s</sup>9121</td><td>0.1[±5.2±0.2±0.8]</td><td>0.3(6)</td><td>-3.9<sup>+1.2</sup><sub>-1.3</sub></td><td>0.1(5)</td><td>1.6<sup>+0.7</sup><sub>-0.7</sub></td><td>-0.65</td><td>-3 × 10<sup>-3</sup></td></tr><tr><td>J1853+1303<sup>b</sup></td><td>57846</td><td>18<sup>b</sup>53<sup>m</sup>57<sup>s</sup>31785</td><td>0.06[±0.8±0.2±0.8]</td><td>13°03′44<sup>s</sup>0471</td><td>0.1[±1.6±0.4±0.8]</td><td>-1.4(1)</td><td>-2.8(2)</td><td>0.49(7)</td><td>0.5<sup>+0.6</sup><sub>-0.3</sub></td><td>-0.37</td><td>0.26</td></tr><tr><td>J1910+1256</td><td>57847</td><td>19<sup>b</sup>10<sup>m</sup>09<sup>s</sup>70165</td><td>0.03[±1.0±0.1±0.8]</td><td>12°56′25<sup>s</sup>4316</td><td>0.06[±2.5±0.1±0.8]</td><td>0.50(4)</td><td>-6.85(9)</td><td>0.254(35)</td><td>0.19<sup>+0.15</sup><sub>-0.12</sub></td><td>-0.48</td><td>0.03</td></tr><tr><td>J1911-1114</td><td>57768</td><td>19<sup>b</sup>11<sup>m</sup>49<sup>s</sup>27544</td><td>0.1[±0.9±0.3±0.8]</td><td>-11°14′22<sup>s</sup>5547</td><td>0.3[±2.1±0.5±0.8]</td><td>-13.8(2)</td><td>-10.3(4)</td><td>0.38<sup>+0.13</sup><sub>-0.14</sub></td><td>1.1<sup>+0.4</sup><sub>-0.3</sub></td><td>-0.39</td><td>-0.02</td></tr><tr><td>J1918-0642</td><td>57768</td><td>19<sup>b</sup>18<sup>m</sup>48<sup>s</sup>02959</td><td>0.1[±0.9±0.2±0.8]</td><td>-06°42′34<sup>s</sup>9335</td><td>0.2[±2.2±0.4±0.8]</td><td>-7.1(1)</td><td>-5.7(3)</td><td>0.60(12)</td><td>1.2<sup>+0.3</sup><sub>-0.2</sub></td><td>-0.39</td><td>-0.01</td></tr><tr><td>J1939+2134<sup>a</sup></td><td>57850</td><td>19<sup>b</sup>39<sup>m</sup>38<sup>s</sup>56134</td><td>[0.07[±0.8±0.1±0.8]</td><td>21°34′59<sup>s</sup>1233</td><td>0.2[±1.9±0.1±0.8]</td><td>0.08(7)</td><td>-0.43(11)</td><td>0.384<sup>+0.048</sup><sub>-0.046</sub></td><td>1.5<sup>+0.7</sup><sub>-0.6</sub></td><td>-0.62</td><td>-0.06</td></tr></tbody></table>

Notes: * α<sub>ref</sub> and δ<sub>ref</sub> refer to the reference position at reference epoch *t*<sub>ref</sub>. The error budgets of the reference positions are provided in the adjacent columns, which include from left to right, the error of relative reference position with respect to the reference point, the uncertainty of the reference point with regard to the main phase calibrator (estimated with Equation 1 of Deller et al. 2019), the position uncertainty of the main phase calibrator, and the typical (0.8 mas in each direction, Sokolovsky et al. 2011) frequency-dependent core shift (e.g. Bartel et al. 1986; Lobanov 1998) between 1.55 GHz and ~8 GHz. We note that the errors outside ‘[’ are obtained with Bayesian inference, while the errors inside ‘]’ are only indicative. To properly determine the absolute pulsar position and its uncertainty requires the procedure described in Section 3.2 of Ding et al. (2020b). This analysis will be made and presented in an upcoming paper.

* ρ<sub>μα,σ</sub> and ρ<sub>μδ,σ</sub> stand for correlation coefficients between σ<sub>μ</sub> and the two proper motion components.

*a* Since inverse referencing is applied for PSR J1939+2134, the two reference sources are the de-facto targets. Accordingly, the proper motion and parallax are the negative values of the direct measurements out of inverse referencing. For the original astrometric model, the reference positions for the two reference sources J194104 and J194106 are, respectively, 19<sup>b</sup>41<sup>m</sup>04<sup>s</sup>319769(2) + 21°49′13′'19731(7) and 19<sup>b</sup>41<sup>m</sup>06<sup>s</sup>86774(1) + 21°53<sup>s</sup>04<sup>s</sup>9594(2), where the uncertainties do not contain those of the reference point (i.e. the inside-the-bracket terms of σ<sub>αref</sub> and σ<sub>δref</sub>). The reference position of PSR J1939+2134 presented in the table is estimated, using normal phase referencing with respect to J194104.

*b* Results of the 8-parameter Bayesian inference are reported here; the constraints on *i* and Ω<sub>asc</sub> are described in Section 5 (see Section 4.3 for explanations).
dominate the total uncertainty in the pulsar's reference position. A future work, incorporating additional multifrequency observations of the in-beam calibrations, will enable significantly more precise pulsar reference positions to be obtained, as is discussed in Section 1.3.

# 5 VLBI + TIMING ASTROMETRIC RESULTS

In Bayesian inference, the output of a model parameter $X_j$ (where $j$ refers to various model parameters) hinges on its prior probability distribution: generally speaking, tighter prior constraints (on $X_j$) that are consistent with data (in the sense of Bayesian analysis) would sharpen the output $X_j$. In cases where a strong correlation between $X_j$ and another model parameter $X_k$ is present, tighter prior $X_j$ constraints that are consistent with the data would potentially sharpen both the output $X_j$ and the output $X_k$.

As noted in Section 1.2, VLBI astrometry serves as the prime method to measure parallaxes of Galactic pulsars. A VLBI astrometric campaign (on a Galactic pulsar) normally spans ~2 yr, as a substantial parallax can likely be achieved in this time span. On the other hand, most MSPSRπ pulsars have been timed routinely for >10 yr, which allows their proper motions to be precisely determined, as the precision on proper motion grows with $t^{3/2}$ (see, e.g. Section 4.4 of Ding, Deller & Miller-Jones 2021b) for a regularly observed pulsar. In Table 6, we collect one timing proper motion (denoted as $\mu_{\alpha}^{(\text{Ti})}$ and $\mu_{\delta}^{(\text{Ti})}$) and one timing parallax ($\bar{\omega}^{(\text{Ti})}$) for each MSPSRπ pulsar. Among the published timing results, we select the timing proper motions measured over the longest time span, and the $\bar{\omega}^{(\text{Ti})}$ having the smallest uncertainties. According to Tables 5 and 6, most timing proper motions are more precise than the quasi-VLBI-only counterparts. On the other hand, timing parallaxes are mostly less precise than the quasi-VLBI-only counterparts. Nevertheless, adopting appropriate timing parallaxes as priors can still effectively lower parallax uncertainties.

The precisely measured $\mu_{\alpha}^{(\text{Ti})}$ and $\mu_{\delta}^{(\text{Ti})}$ provide the opportunity to significantly refine the quasi-VLBI-only proper motions. Furthermore, as shown with the Pearson correlation coefficients (Pearson 1895), $\rho_{\mu_{\alpha},\bar{\omega}}$ and $\rho_{\mu_{\delta},\bar{\omega}}$, that we summarized in Table 5, large correlation between parallax and proper motion is not rare for VLBI astrometry. Therefore, using the $\mu_{\alpha}^{(\text{Ti})}$ and $\mu_{\delta}^{(\text{Ti})}$ measurements as the prior proper motion constraints in Bayesian inference can potentially refine both proper motion and parallax determination.

The astrometric results inferred with timing priors, hereafter referred to as VLBI + timing results, are reported in Table 6. To differentiate from the notation of quasi-VLBI-only astrometric parameter $Y$, we denote a VLBI + timing model parameter in the form of $Y'$. Comparing Tables 5 and 6, we find almost all VLBI + timing proper motions and parallaxes more precise than the quasi-VLBI-only counterparts; the most significant parallax precision enhancement occurs to PSR J1918−0642 (by 42 per cent), followed by PSR J1939+2134 (by 36 per cent) and PSR J1537+1155 (by 33 per cent). Hence, we use the VLBI + timing results in the remainder of this paper.

In 7 cases (i.e. PSR J0610−2100, PSR J1643−1224, PSR J1730−2304, PSR J1738+0333, PSR J1853+1303, PSR J1824−2452A, PSR J1910+1256), one of $\mu_{\alpha}^{(\text{Ti})}$, $\mu_{\delta}^{(\text{Ti})}$, or $\bar{\omega}^{(\text{Ti})}$ is more than $2\sigma$ discrepant from the quasi-VLBI-only counterpart. Using such timing priors may widen the uncertainties of resultant model parameters, as $\eta_{\text{EFAC}}$ would be lifted to counter-balance the increased $\chi^2$. Without any indication that the discrepant timing values are less reliable, we use them as priors regardless. However, we caution the use of these seven sets of VLBI + timing results, and

would recommend the quasi-VLBI-only results to be considered if our adopted timing priors are proven inaccurate in future.

We also now consider any possible effects that could, despite our best efforts to characterise all sources of position noise, bias the fitted VLBI positions. For any given VLBI calibrator source, evolution in the source structure can lead to a detectable position offset (e.g. Perger et al. 2018; Zhang, An & Frey 2020) that is then transferred to the target pulsar. Due to the long time-scales of AGN structure evolution, over the ~2-yr time-scale of the MSPSRπ observations, this error may be quasi-linear in time and be absorbed into the pulsar proper motion (e.g. Deller et al. 2013). Redundant secondary calibrators can be used to probe the astrometric effect of structure evolution. However, with small numbers of redundant calibrator sources, such probes are hardly conclusive, as the structure evolution of the redundant calibrators would also be involved. Among the seven pulsars showing >2 $\sigma$ discrepancy between quasi-VLBI-only and timing results (see Table 6), PSR J0030+0451, PSR J1643−1224, PSR J1730−2304, PSR J1738+0333, and PSR J1824−2452A either display no relative motion between the redundant secondary calibrators, and the main secondary calibrators or do not have any redundant calibrator (i.e. PSR J1643−1224), although the suboptimal main secondary calibrators of PSR J1643−1224 and PSR J1824−2452A (see Sections 8.5 and 8.9) may likely affect the astrometric performance. For PSR J1853+1303, the main secondary calibrator has a clear jet aligned roughly with the right ascension (RA) direction, and thus source structure evolution is potentially significant. The two redundant calibrators for PSR J1853+1303 do display a relative proper motion of up to 0.2 mas yr$^{-1}$ with respect to the main secondary calibrator, so while the mean relative motion seen between the two redundant secondary calibrators is small, calibrator structure evolution remains a possible explanation for the VLBI-timing discrepancy. Finally, the main secondary calibrator of PSR J1910+1256 also exhibits a jet structure at a position angle of ~45°. When using the only redundant calibrator of PSR J1910+1256 as the reference source, we obtained the VLBI-only result $\mu_{\alpha} = 0.25 \pm 0.06$ mas yr$^{-1}$, $\mu_{\delta} = -7.3 \pm 0.1$ mas yr$^{-1}$, and $\bar{\omega} = 0.61 \pm 0.05$ mas with Bayesian inference, where $\mu_{\alpha}$ becomes consistent with $\mu_{\alpha}^{(\text{Ti})}$, but $\mu_{\delta}$ and $\bar{\omega}$ are further away from the timing counterparts. The $\mu_{\alpha}$ consistency between VLBI and timing indicates that structure evolution in our chosen calibrator is likely contributing to the VLBI-timing discrepancy. However, as the redundant calibrator is both fainter and further away from PSR J1910+1256 (compared to the main secondary calibrator), we do not use this source as the final reference source.

## 5.1 The posterior orbital inclinations and ascending node longitudes

For the four 8P pulsars, orbital inclinations $i'$, and ascending node longitudes $\Omega'_{\text{asc}}$ are also inferred alongside the five canonical parameters and $\eta_{\text{EFAC}}$ (see Section 4.3). The full 8D corner plots out of the 8-parameter inferences are available online⁵. Prior constraints on $i'$ and $\Omega'_{\text{asc}}$ have been provided in Section 4.0.1. Owing to bi-modal features of all 1D histograms of $i'$, no likelihood component is substantially favoured over the other. Hence, no tight posterior constraint on $i'$ is achieved for any 8P pulsar. Likewise, all 1D histograms of $\Omega'_{\text{asc}}$ show multimodal features, which precludes stringent constraints on $\Omega'_{\text{asc}}$.

## 5.2 Comparison with Gaia results

From the Gaia Data Release 2 (Gaia Collaboration et al. 2018), Gaia counterparts for pulsars with optically bright companions
<table><caption>**Table 6.** Left of the dashed line: proper motions {$\mu'_{\alpha}$, $\mu'_{\delta}$} and parallaxes $\sigma'$ inferred with timing proper motion and parallax priors. Right: distances *D* and transverse space velocities $v_{\perp}$ based on $\mu'_{\alpha}$, $\mu'_{\delta}$, and $\sigma'$.</caption><thead><tr><th>PSR</th><th>μ<sup>(T)</sup><sub>α</sub><br>(mas yr<sup>-1</sup>)</th><th>μ<sup>(T)</sup><sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>σ<sup>(T)</sup><sub>α</sub><br>(mas)</th><th>μ'<sub>α</sub><br>(mas yr<sup>-1</sup>)</th><th>μ'<sub>δ</sub><br>(mas yr<sup>-1</sup>)</th><th>σ'<sub>α</sub><br>(mas)</th><th>η'<sub>EFAC</sub></th><th>DM<br>(pc cm<sup>-3</sup>)</th><th>d<sup>(NE)</sup><sub>DM</sub><br>(kpc)</th><th>d<sup>(YMW)</sup><sub>DM</sub><br>(kpc)</th><th>b<br>(deg)</th><th>D<br>(kpc)</th><th>v<sub>⊥</sub><sup>**</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0030+0451</td><td>-6.2(1)</td><td>0.5(3)<sup>d</sup></td><td>3.08(8)<sup>d</sup></td><td>-6.15(5)</td><td>0.37(14)</td><td>3.04(5)</td><td>1.3<sup>+0.4</sup><sub>-0.3</sub></td><td>4.3</td><td>0.32(6)</td><td>0.35(7)</td><td>-57.6</td><td>0.329<sup>+0.006</sup><sub>-0.005</sub></td><td>15.4(2)</td></tr><tr><td>J0610-2100</td><td>9.04(8)</td><td>"16.7(1)<sup>c</sup></td><td>—</td><td>9.06(7)</td><td>16.6(1)</td><td>0.72(11)</td><td>1.4<sup>+0.4</sup><sub>-0.3</sub></td><td>60.7</td><td>3.5(7)</td><td>3.3(7)</td><td>-18.2</td><td>1.5<sup>+0.3</sup><sub>-0.2</sub></td><td>120<sup>+25</sup><sub>-17</sub></td></tr><tr><td>J0621+1002</td><td>"3.2(1)</td><td>"0.6(5)<sup>c</sup></td><td>—</td><td>3.27(9)</td><td>-1.1(3)</td><td>0.74(14)</td><td>3.9<sup>+1.1</sup><sub>-0.8</sub></td><td>36.5</td><td>1.4(3)</td><td>0.42(8)</td><td>-2.0</td><td>1.6<sup>+0.5</sup><sub>-0.3</sub></td><td>20<sup>+8</sup><sub>-5</sub></td></tr><tr><td>J1012+5307</td><td>"2.61(1)</td><td>-25.49(1)<sup>c</sup></td><td>"0.9(2)<sup>c</sup></td><td>2.61(1)</td><td>-25.49(1)</td><td>1.14(4)</td><td>1.7<sup>+0.5</sup><sub>-0.4</sub></td><td>9.02</td><td>0.41(8)</td><td>0.8(2)</td><td>50.9</td><td>0.877(35)</td><td>95(4)</td></tr><tr><td>J1024-0719</td><td>-35.270(17)</td><td>-48.22(3)<sup>b</sup></td><td>0.83(13)<sup>b</sup></td><td>-35.27(2)</td><td>-48.22(3)</td><td>0.93(5)</td><td>1.2<sup>+0.4</sup><sub>-0.3</sub></td><td>6.5</td><td>0.39(8)</td><td>0.38(8)</td><td>40.5</td><td>1.08(6)</td><td>300(20)</td></tr><tr><td>J1518+4904</td><td>-0.67(4)</td><td>-8.53(4)<sup>b</sup></td><td>—</td><td>-0.683(26)</td><td>-8.528(36)</td><td>1.237<sup>+0.035</sup><sub>-0.031</sub></td><td>1.2<sup>+0.6</sup><sub>-0.4</sub></td><td>11.61</td><td>0.6(1)</td><td>1.0(2)</td><td>54.3</td><td>0.81(2)</td><td>16.0(6)</td></tr><tr><td>J1537+1155</td><td>1.482(7)</td><td>-25.285(12)<sup>g</sup></td><td>0.86(18)<sup>g</sup></td><td>1.484(7)</td><td>-25.286(11)</td><td>1.07(7)</td><td>0.5<sup>+0.5</sup><sub>-0.3</sub></td><td>11.62</td><td>1.0(2)</td><td>0.9(2)</td><td>48.3</td><td>0.94<sup>+0.07</sup><sub>-0.06</sub></td><td>102<sup>+8</sup><sub>-7</sub></td></tr><tr><td>J1640+2224</td><td>2.08(1)</td><td>-11.34(2)<sup>c</sup></td><td>0.6(4)<sup>c</sup></td><td>2.08(1)</td><td>-11.34(2)</td><td>0.73(6)</td><td>1.3<sup>+0.7</sup><sub>-0.6</sub></td><td>18.43</td><td>1.2(2)</td><td>1.5(3)</td><td>38.3</td><td>1.39<sup>+0.13</sup><sub>-0.11</sub></td><td>53<sup>+5</sup><sub>-4</sub></td></tr><tr><td>J1643-1224</td><td>"5.970(18)</td><td>3.77(8)<sup>b</sup></td><td>"0.82(17)<sup>b</sup></td><td>5.97(2)</td><td>3.76(8)</td><td>1.1(1)</td><td>1.0<sup>+0.3</sup><sub>-0.2</sub></td><td>62.3</td><td>2.4(5)</td><td>0.8(2)</td><td>21.2</td><td>0.95<sup>+0.15</sup><sub>-0.11</sub></td><td>41<sup>+4</sup><sub>-4</sub></td></tr><tr><td>J1721-2457</td><td>1.9(1.2)</td><td>"25(16)<sup>f</sup></td><td>—</td><td>2.5(3)</td><td>-1.9(9)</td><td>0.0(2)</td><td>3.1<sup>+0.8</sup><sub>-0.6</sub></td><td>48.3</td><td>1.3(3)</td><td>1.4(3)</td><td>6.8</td><td>&gt; 1.5 ▲</td><td>—</td></tr><tr><td>J1730-2304</td><td>"20.06(12)</td><td>-4(2)<sup>b</sup></td><td>"2.11(11)<sup>b</sup></td><td>20.1(1)</td><td>-4.7(6)</td><td>2.0(1)</td><td>1.6<sup>+0.4</sup><sub>-0.3</sub></td><td>9.62</td><td>0.5(1)</td><td>0.5(1)</td><td>6.0</td><td>0.51(3)</td><td>54<sup>+3</sup><sub>-2</sub></td></tr><tr><td>J1738+0333</td><td>7.037(5)</td><td>5.073(12)<sup>c</sup></td><td>"0.68(5)<sup>c</sup></td><td>7.036(5)</td><td>5.07(1)</td><td>0.589(46)</td><td>2.3<sup>+0.8</sup><sub>-0.6</sub></td><td>33.8</td><td>1.4(3)</td><td>1.5(3)</td><td>17.7</td><td>1.74<sup>+0.15</sup><sub>-0.13</sub></td><td>90<sup>+7</sup><sub>-6</sub></td></tr><tr><td>J1824-2452A</td><td>-0.25(4)</td><td>"" - 8.6(8)<sup>b</sup></td><td>—</td><td>-0.25(4)</td><td>-7.8(8)</td><td>0.4(5)</td><td>2.5<sup>+0.8</sup><sub>-0.6</sub></td><td>119.9</td><td>3.1(6)</td><td>3.7(7)</td><td>-5.6</td><td>&gt;0.5 ▲</td><td>—</td></tr><tr><td>J1853+1303</td><td>" - 1.63(2)</td><td>-2.96(4)<sup>c</sup></td><td>0.48(14)<sup>c</sup></td><td>-1.62(2)</td><td>-2.96(4)</td><td>0.53(7)</td><td>0.9(5)</td><td>30.6</td><td>2.1(4)</td><td>1.3(3)</td><td>5.4</td><td>2.0(3)</td><td>14<sup>+3</sup><sub>-2</sub></td></tr><tr><td>J1910+1256</td><td>"""0.21(3)</td><td>" - 7.04(5)<sup>c</sup></td><td>0.1(3)<sup>c</sup></td><td>0.24(3)</td><td>-7.03(5)</td><td>0.36(6)</td><td>0.8<sup>+0.3</sup><sub>-0.2</sub></td><td>38.1</td><td>2.3(5)</td><td>1.5(3)</td><td>1.8</td><td>3.4<sup>+1.0</sup><sub>-0.6</sub></td><td>55<sup>+16</sup><sub>-9</sub></td></tr><tr><td>J1911-1114</td><td>-13.7(2)</td><td>" - 9.3(9)<sup>c</sup></td><td>—</td><td>-13.7(1)</td><td>-10.2<sup>+0.4</sup><sub>-0.3</sub></td><td>0.37<sup>+0.13</sup><sub>-0.14</sub></td><td>1.1<sup>+0.4</sup><sub>-0.3</sub></td><td>31.0</td><td>1.2(2)</td><td>1.1(2)</td><td>-9.6</td><td>&gt; 1.3 ▲</td><td>—</td></tr><tr><td>J1918-0642</td><td>-7.15(2)</td><td>-5.94(5)<sup>c</sup></td><td>"0.8(1)<sup>c</sup></td><td>-7.15(2)</td><td>-5.93(5)</td><td>0.71(7)</td><td>1.1<sup>+0.3</sup><sub>-0.2</sub></td><td>6.1</td><td>1.2(2)</td><td>1.0(2)</td><td>-9.1</td><td>1.48<sup>+0.19</sup><sub>-0.14</sub></td><td>44<sup>+6</sup><sub>-5</sub></td></tr><tr><td>J1939+2134</td><td>0.074(2)</td><td>-0.410(3)<sup>c</sup></td><td>"0.28(5)<sup>c</sup></td><td>0.074(2)</td><td>-0.410(3)</td><td>0.35(3)</td><td>1.2(6)</td><td>71.1</td><td>3.6(7)</td><td>2.9(6)</td><td>-0.3</td><td>2.9<sup>+0.3</sup><sub>-0.2</sub></td><td>80<sup>+9</sup><sub>-8</sub></td></tr></tbody></table>

Notes. •T" denotes historical pulsar timing results. ρ stand for correlation coefficients of Bayesian inference without using timing proper motion priors.

•The 't's before $\mu'_{\alpha}$, $\mu'_{\delta}$ and $\sigma'$ convey the significance of the offset from the quasi-VLBI-only counterparts. In specific, the 't' repetition number $N_t$ means the offset significance is between $N_t$ and $(N_t+1)\sigma$. VLBI + timing results obtained with timing proper motion or parallax priors that are more than 2 $\sigma$ away from the quasi-VLBI-only counterparts (i.e. $\mu'_{\alpha}$, $\mu'_{\delta}$ and $\sigma'$) marked with $\ge 2$ 't's) should be used with caution.

•$\mu'_{\alpha}$, $\mu'_{\delta}$, and $\sigma'$ designate results of Bayesian inference using timing proper motion priors.

•For each pulsar, we present the most precise timing estimates published in $^{a}$Faisal Alam et al. (2021), $^{b}$Reardon et al. (2021), $^{c}$Perera et al. (2019), $^{d}$Arzoumanian et al. (2018), $^{e}$Freire et al. (2012), $^{f}$Desvignes et al. (2016), $^{g}$Fonseca et al. (2014), and $^{h}$Janssen et al. (2008).

•For comparison, we list $d_{\text{DM}}^{(NE)}$ and $d_{\text{DM}}^{(YMW)}$ derived with pygedm<fnref n="9" /> given the sky position and the DM of a pulsar, based on the two latest $n_e(\bar{x})$ models (Cordes & Lazio 2002; Yao, Manchester & Wang 2017).

*20 per cent relative uncertainties are assumed for all DM-based distances (i.e. $d_{\text{DM}}^{(NE)}$ and $d_{\text{DM}}^{(YMW)}$).

**Tangential space velocities corrected for the differential rotation of the Galaxy (see 6.2).

▲ The reciprocal of the 3 $\sigma$ upper limit of the parallax is adopted as the lower limit of the distance.

<footnote n="9">⁹https://github.com/FRBs/pygedm</footnote>
<table><caption>**Table 7.** Gaia astrometric results.</caption><thead><tr><th>PSR</th><th>Gaia DR3<br>source ID</th><th>μ<sub>α</sub><sup>(G)</sup><br>(mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub><sup>(G)</sup><br>(mas yr<sup>-1</sup>)</th><th>ω<sub>1</sub><sup>(G)</sup><br>(mas)</th><th>GoF.*</th></tr></thead><tbody><tr><td>J1012+5307</td><td>851610861391010944</td><td>2.7(3)</td><td>¹.25.9(3)</td><td>¹.7(3)</td><td>-1.5</td></tr><tr><td>J1024-0719</td><td>3775277872387310208</td><td>-35.5(3)</td><td>-48.35(36)</td><td>0.86(28)</td><td>0.4</td></tr><tr><td>J1910+1256</td><td>4314046781982561920<sup>?</sup></td><td>!!!-2.3(6)</td><td>¹.6.1(6)</td><td>-0.1(8)</td><td>1.9</td></tr></tbody></table>

Notes. •Sources marked with ‘?’ are tentative Gaia counterpart candidates.

●Values marked with N?!’s are Nσ − (N + 1)σ offset from the VLBI + timing counterparts

*Goodness of fitting, a parameter (of Gaia data releases) approximately following N(0, 1) distribution. A GoF closer to zero indicates better fitting performance.

have been identified and studied by Jennings et al. (2018), Mingarelli et al. (2018), Antoniadis (2021). In the MSPSRπ sample, PSR J1012+5307 and PSR J1024−0719 have secure Gaia counterparts, while PSR J1910+1256 has a proposed Gaia counterpart candidate (Mingarelli et al. 2018). In Table 7, we updated the Gaia results for these three Gaia sources to the Gaia Data Release 3 (DR3, Gaia Collaboration et al. 2022).

For PSR J1024−0719, the Gaia proper motion {$\mu_{\alpha}^{(G)}$, $\mu_{\delta}^{(G)}$} and parallax $\omega_{1}^{(G)}$ are highly consistent with the VLBI + timing ones, which further strengthens the proposal that PSR J1024−0719 is in an ultra-wide orbit with a companion star (Bassa et al. 2016; Kaplan et al. 2016, also see Sections 6.2 and 7.2). The Gaia proper motion and parallax of PSR J1012+5307 is largely consistent with the VLBI + timing counterparts. The > 1 σ discrepancy between $\mu_{\delta}^{(G)}$ and $\omega_{1}^{(G)}$ and the respective VLBI + timing counterparts can be explained by non-optimal goodness of (Gaia astrometric) fitting (GoF) (see Table 7). On the other hand, the Gaia counterpart candidate for PSR J1910+1256 (proposed by Mingarelli et al. 2018) possesses a $\mu_{\alpha}^{(G)}$ 4 σ discrepant from the VLBI + timing one. Though this discrepancy is discounted by the relatively bad GoF by roughly a factor of 1.9 (see Table 7), the connection between the Gaia source and PSR J1910+1256 remains inconclusive. We note that the parallax zero-points $\omega_{0}^{(G)}$ (Lindegren et al. 2021) of the three Gaia sources are negligible and hence not considered, as $\omega_{0}^{(G)}$ is small ($|\omega_{0}^{(G)}| \lesssim 0.02$ mas, Ding et al. 2021b) compared to the uncertainty of $\omega_{1}^{(G)}$ (see Table 7).

# 6 DISTANCES AND SPACE VELOCITIES

In this section, we derive pulsar distances *D* from parallaxes $\omega'$ (see Section 5), and compare them to the dispersion-measure-based distances. Incorporating the proper motions {$\mu'_{\alpha}$, $\mu'_{\delta}$} (see Section 5), we infer the transverse space velocity $v_{\perp}$ (i.e. the velocity with respect to the stellar neighbourhood) for each pulsar in an effort to enrich the sample of ~40 MSPs with precise $v_{\perp}$ (Hobbs et al. 2005; Gonzalez et al. 2011), and refine the $v_{\perp}$ distributions of MSP subgroups such as binary MSPs and solitary MSPs.

## 6.1 Parallax-based distances

Inferring a source distance from a measured parallax requires assumptions about the source properties, for which a simple inversion implicitly makes unphysical assumptions (e.g. Bailer-Jones et al. 2021). Various works (e.g. Lutz & Kelker 1973; Verbiest et al. 2012; Bailer-Jones 2015; Igoshev, Verbunt & Cator 2016) have contributed to developing and consolidating the mathematical formalism of parallax-based distance inference, which we briefly recapitulate as follows, in order to facilitate comprehension and ready the mathematical formalism for further discussion.

A parallax-based distance *D* can be approached from the conditional probability density function (PDF)

$$p(D|\omega', l, b) \propto p(\omega'|D)p(D, l, b), \quad (7)$$

where *l* and *b* stands for Galactic longitude and latitude, respectively; $\omega' = \omega'_{0} \pm \sigma_{\omega'}$. The first term on the right takes the form of

$$p(\omega'|D) \propto \exp \left[ -\frac{1}{2} \left( \frac{1/D - \omega'_{0}}{\sigma_{\omega'}} \right)^{2} \right], \quad (8)$$

assuming $\omega'_{0}$ is Gaussian-distributed, or more specifically, $\omega'_{0} \sim N(1/D, \sigma_{\omega'}^{2})$. The second term on the right side of equation (7) can be approximated as $p(D, l, b) \propto D^{2}$, when the parent population $\Psi$ of the target celestial body is uniformly distributed spatially (Lutz & Kelker 1973). Given a postulated (Galactic) spatial distribution $\rho(D, l, b)$ of $\Psi$, $p(D, l, b) \propto D^{2}\rho(D, l, b)$. Hence,

$$p(D|\omega', l, b) \propto D^2 \rho(D, l, b) \exp \left[ -\frac{1}{2} \left( \frac{1/D - \omega'_{0}}{\sigma_{\omega'}} \right)^2 \right]. \quad (9)$$

We join Verbiest et al. (2012) and Jennings et al. (2018) to adopt the $\rho(D, l, b)$ (of the ‘Model C’) determined by Lorimer et al. (2006) for Galactic pulsars. While calculating the $\rho(D, l, b)$ with Equations 10 and 11 of Lorimer et al. (2006), we follow Verbiest et al. (2012) and Jennings et al. (2018) to increase the scale height (i.e. the parameter ‘E’ of Lorimer et al. 2006) to 0.5 kpc to accommodate the MSP population. In addition, the distance to the Galactic Centre (GC) in Equation 10 of Lorimer et al. 2006 is updated to $d_{\odot} = 8.12 \pm 0.03$ kpc (Gravity Collaboration et al. 2018). We do not follow Verbiest et al. (2012), Igoshev et al. (2016) to use pulsar radio fluxes to constrain pulsar distances, as pulsar luminosity is relatively poorly constrained.

Using the aforementioned mathematical formalism, we calculated $p(D|\omega', l, b)$ for each MSPSRπ pulsar, and integrated it into the cumulative distribution function (CDF) $\Phi(D|\omega', l, b) = \int_{0}^{D} p(D'|\omega', l, b) dD'$. The $p(D|\omega', l, b)$ and $\Phi(D|\omega', l, b)$ is plotted for each pulsar and made available online⁵. An example of these plots are presented in Fig. 4. The median distances $D_{median}$ corresponding to $\Phi(D|\omega', l, b) = 0.5$ are taken as the pulsar distances, and summarized in Table 6. The distances matching $\Phi(D|\omega', l, b) = 0.16$ and $\Phi(D|\omega', l, b) = 0.84$ are respectively used as the lower and upper bound of the 1 σ uncertainty interval.

### 6.1.1 Comparison with DM distances

As mentioned in Section 1.2, the precise DM measured from a pulsar can be used to assess the pulsar distance, provided an $n_{e}(\bar{x})$ model. Using pygeadm⁹, we compile into Table 6 the DM distances (i.e. $d_{DM}^{(NE)}$ and $d_{DM}^{(YMW)}$) of each pulsar based on the two latest realizations of $n_{e}(\bar{x})$ model – the NE2001 model (Cordes & Lazio 2002) and the YMW16 model (Yao et al. 2017). For all the DM distances, we
<figure><img src="image_4.png" /><figcaption>**Figure 4.** An example posterior probability density function $p(D|\omega', l, b)$ (of distance) and its cumulative distribution function $\Phi(D|\omega', l, b) = \int_0^D p(D'|\omega', l, b) dD'$. The vertical dashed lines correspond to $\Phi(D|\omega', l, b) = 0.16$ and $\Phi(D|\omega', l, b) = 0.84$, which are, respectively, used as the lower and upper bound of the $1\sigma$ uncertainty interval. The mode distance $D_{\text{mode}}$ and median distance $D_{\text{median}}$ are marked with dot-dashed blue line and dotted cyan line, respectively. Plots of this kind are also made for other MSPSR$\pi$ pulsars, and made available online$^5$. Staying in line with the norm (see Section 1.2) of this paper, we universally adopt $D_{\text{median}}$ as the distances (i.e. $D$ in Table 6) for all MSPSR$\pi$ pulsars in this paper.</figcaption></figure>

adopt typical 20 per cent fractional uncertainties. We have obtained significant ($\ge 3\sigma$) parallax-based distances $D$ for 15 out of 18 MSPSR$\pi$ pulsars. These distances enable an independent quality check of both $n_e(\bar{x})$ models.

Among the 15 pulsars with parallax-based distance measurements, YMW16 is more accurate than NE2001 in three cases (i.e. PSR J1012+5307, PSR J1643-1224, and PSR J1939+2134), but turns out to be the other way around in four cases (i.e. PSR J0621+1002, PSR J1853+1303, PSR J1910+1256, and PSR J1918-0642). In other eight cases, the $D$ cannot discriminate between the two models. The small sample of 15 $D$ measurements shows that NE2001 and YMW16 remain comparable in terms of outliers. In two (out of the 15) cases (i.e. PSR J0610-2100, PSR J1024-0719), $D$ is about $2.6\sigma$ and $6.8\sigma$ away from either DM distance, which reveals the need to further refine the $n_e(\bar{x})$ models. Such a refinement can be achieved, with improved pulsar distances including the ones determined in this work.

## 6.2 Transverse space velocities

Having determined the parallax-based distances $D$ and the proper motions $\{\mu'_\alpha, \mu'_\delta\}$, we proceed to calculate transverse space velocities $v_\perp$ for each pulsar, namely the transverse velocity with respect to the neighbouring star field of the pulsar. In estimating the transverse velocity of a pulsar neighbourhood, we assume the neighbourhood observes circular motion about the axis connecting the North and South Galactic Poles, which is roughly valid given that all MSPSR$\pi$ pulsars with significant ($> 3\sigma$) $D$ share a median $|z| = D \sin|b|$ of 0.3 kpc. Using the Galactic rotation curve from Reid et al. (2019) and the full circular velocity of the Sun $247 \pm 1$ km s$^{-1}$, we derived the apparent transverse velocity of the neighbourhood $v_\perp$, $N$, thus obtaining $v_\perp$ by subtracting the apparent transverse velocity of the pulsar by $v_\perp$, $N$. Here, the full circular velocity (denoted as $\Theta_0 + V_\odot$)

in Reid et al. 2019) is calculated with $d_\odot = 8.12 \pm 0.03$ kpc (Gravity Collaboration et al. 2018) and the proper motion of Sgr A* from Reid et al. (2019).

To estimate the uncertainty of $v_\perp$, we simulated a chain of 50 000 distances for each pulsar based on the $p(D|\omega', l, b)$ that we have obtained in Section 6.1. Besides, we also acquired chains of 50 000 $\mu'_\alpha$ and $\mu'_\delta$ given the VLBI + timing proper motions of Table 6, assuming $\mu'_\alpha$ and $\mu'_\delta$ follow Gaussian distributions. With these chains of $D$, $\mu'_\alpha$, and $\mu'_\delta$, we calculated 50 000 $v_\perp$ values, which form a PDF of $v_\perp$ for each pulsar. The $v_\perp$ inferred from the PDFs are summarized in Table 6.

In Fig. 5, we illustrate the $v_\perp$ in relation to $|z|$ for 16 pulsars with precise distance estimates. Among the 16 pulsars, only PSR J1824-2452A does not have a significant parallax-based distance. Nevertheless, its $v_\perp$ can be inferred by incorporating its proper motion with the astrometric information (i.e. distance and proper motion) of its host globular cluster (see Section 8.9). No clear correlation is revealed between $v_\perp$ and $|z|$, which reinforces our decision to treat all MSPSR$\pi$ pulsars across the $|z| \lesssim 1$ kpc regime equally. By concatenating the simulated $v_\perp$ chains, we acquired the PDF for the 16 MSPs (see Fig. 5), which gives $v_\perp^{(\text{MSP})} = 53^{+48}_{-37}$ km s$^{-1}$.

Amongst the MSPSR$\pi$ sources, PSR J1024-0719 is an obvious outlier, with a velocity of $\sim 300$ km s$^{-1}$ that is $3\sigma$ above the mean. As proposed by Bassa et al. (2016) and Kaplan et al. (2016), PSR J1024-0719 is theorized to have been ejected from a dense stellar region, thus possibly following a different $v_\perp$ distribution from typical field MSPs (isolated along with their respective companions throughout their lives). In this regard, we turn our attention to the binary sample of pulsars with well determined orbital periods $P_b$ (see $P_b$ of Table 3), and obtain $v_\perp^{(\text{BI})} = 50^{+49}_{-34}$ km s$^{-1}$ for field binary MSPs. Based on this small sample, we do not find the $v_\perp$ of the three solitary MSPs (i.e. PSR J0030+0451, PSR J1730-2304, and PSR J1939+2134) to be inconsistent with $v_\perp^{(\text{BI})}$. Neither are the two DNSs (i.e. PSR J1518+4904 and PSR J1537+1155). If we exclude the two DNSs from the binary sample, we would come to $v_\perp^{(\text{WD})} = 50^{+46}_{-31}$ km s$^{-1}$ for the MSPSR$\pi$ pulsars with WD companions, which is highly consistent with $v_\perp^{(\text{BI})}$ and $v_\perp^{(\text{MSP})}$.

Compared to $113 \pm 17$ km s$^{-1}$ previously estimated for a sample of $\sim 40$ MSPs (Gonzalez et al. 2011), our $v_\perp^{(\text{MSP})}$ is largely consistent but on the smaller side. Boodram & Heinke (2022) recently shows that MSP space velocities have to be near zero to explain the Galactic Centre $\gamma$-ray excess (e.g. Abazajian & Kaplinghat 2012). Interestingly, the $v_\perp$ PDF based on our small sample of 16 shows a multimodal feature, with the lowest mode consistent with zero. Specifically, the seven MSPSR$\pi$ pulsars with the smallest $v_\perp$ share an equally weighted mean $v_\perp$ of only 25 km s$^{-1}$, which suggests MSPs with extremely low-space velocities are not uncommon. Accordingly, we suspect the MSP origin of the GC $\gamma$-ray excess can still not be ruled out based on our sample of $v_\perp$.

## 7 RADIAL ACCELERATIONS OF PULSARS AND ORBITAL-DECAY TESTS OF GRAVITATIONAL THEORIES

As described in Section 1.2, VLBI astrometry of pulsars, in conjunction with pulsar timing, can enhance the orbital-decay tests of gravitational theories. For binary systems involved in this work, the observed orbital decay has three significant components:

$$ \dot{P}_b^{\text{obs}} = \dot{P}_b^{\text{GW}} + \dot{P}_b^{\text{Shk}} + \dot{P}_b^{\text{Gal}}, \quad (10) $$

where $\dot{P}_b^{\text{GW}}$ reflects the effect of gravitational-wave damping intrinsic to a binary system, while $\dot{P}_b^{\text{Shk}}$ and $\dot{P}_b^{\text{Gal}}$ are both extrinsic contribu-
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Upper: The transverse space velocities $v_⊥$ versus the Galactic vertical heights $|z| = D \sin|b|$ of the 16 MSPSR$\pi$ pulsars with significant (> 3$\sigma$) distance measurements (including 15 parallax-based distances and a globular cluster distance). Lower: The probability density function (PDF) of $v_⊥$ for the 16 MSPs. The median of the $v_⊥$ PDF is marked with the dashed line, while the 1$\sigma$ error interval is shown with the shaded region.</figcaption></figure>

tions caused, respectively, by relative line-of-sight accelerations (of pulsars) $A_{\text{Shk}}$ and $A_{\text{Gal}}$. Specifically, $\dot{P}_{\text{b}}^{\text{Shk}} = A_{\text{Shk}}/c \cdot P_{\text{b}} = \mu^2 D/c \cdot P_{\text{b}}$ (where $\mu^2 = \mu_{\alpha}^{\prime 2} + \mu_{\delta}^{\prime 2}$) is the radial acceleration caused by the tangential motion of pulsars (Shklovskii 1970), which becomes increasingly crucial for pulsars with larger $\mu$ (e.g. PSR J1537+1155, Ding et al. 2021a), as $A_{\text{Shk}} \propto \mu^2$. On the other hand,

$$ \dot{P}_{\text{b}}^{\text{Gal}} = \frac{A_{\text{Gal}}}{c} P_{\text{b}} = \frac{[-\nabla\varphi(\vec{x})]_{\mid_{\vec{x}_{\odot}}} \cdot \vec{e}_{\text{r}}}{c} P_{\text{b}} \quad (11) $$

is a consequence of the gravitational pull (or push) exerted by the Galaxy. Here, $\varphi(\vec{x})$ and $\vec{e}_{\text{r}}$ are, respectively, the Galactic gravitational potential (as a function of Galactic position $\vec{x}$) and the unit vector in the Earth-to-pulsar direction.

In order to test any theoretical prediction of $\dot{P}_{\text{b}}^{\text{GW}}$, it is necessary to estimate $A_{\text{Shk}}$ and $A_{\text{Gal}}$ and remove their effect on $\dot{P}_{\text{b}}^{\text{obs}}$. Besides this impact, the radial accelerations $A_{\text{Shk}}$ and $A_{\text{Gal}}$ would, more generally, affect the time derivative of all periodicities intrinsic to a pulsar system, which include the pulsar spin period derivative $\dot{P}_{\text{s}}$. Similar to $\dot{P}_{\text{b}}^{\text{Shk}}$ and $\dot{P}_{\text{b}}^{\text{Gal}}$, $\dot{P}_{\text{s}}^{\text{Shk}} = A_{\text{Shk}}/c \cdot P_{\text{s}}$ and $\dot{P}_{\text{s}}^{\text{Gal}} = A_{\text{Gal}}/c \cdot P_{\text{s}}$ (where $P_{\text{s}}$ stands for the spin period of a pulsar). As MSPs consist of nearly half of the $\gamma$-ray pulsar population, determining the extrinsic terms of $\dot{P}_{\text{s}}$ and the intrinsic spin period derivative $\dot{P}_{\text{s}}^{\text{int}} = \dot{P}_{\text{s}}^{\text{obs}} - \dot{P}_{\text{s}}^{\text{Shk}} - \dot{P}_{\text{s}}^{\text{Gal}}$ is essential for exploring the ‘death line’ (i.e. the lower limit) of high-energy emissions from pulsars (e.g. Guillemot et al. 2016).

Sections 7.1 and 7.2, we evaluate $A_{\text{Shk}}$ and $A_{\text{Gal}}$ one after another. The evaluation only covers pulsars with significant $D$, as both $A_{\text{Shk}}$ and $A_{\text{Gal}}$ are distance-dependent.

## 7.1 Shklovkii effects

We estimate the model-independent $A_{\text{Shk}}$ in a way similar to the estimation of $v_⊥$ (see Section 6.2). Three chains of 50 000 $\mu_{\alpha}^{\prime}$, $\mu_{\delta}^{\prime}$, and $D$ were simulated from their respective PDFs. Using the relation $A_{\text{Shk}} = (\mu_{\alpha}^{\prime 2} + \mu_{\delta}^{\prime 2}) D$, 50 000 $A_{\text{Shk}}$ were calculated to assemble the PDF of $A_{\text{Shk}}$ for each pulsar with significant $D$. The $A_{\text{Shk}}$ inferred from the PDFs are compiled in Table 8 along with their resultant $\dot{P}_{\text{s}}^{\text{Shk}}$ and $\dot{P}_{\text{b}}^{\text{Shk}}$.

## 7.2 Relative radial accelerations due to Galactic gravitational pull

We estimate $A_{\text{Gal}}$ in the same way as Ding et al. (2021a), following the pioneering work of Zhu et al. (2019). To briefly demonstrate this method, we present, in Table 9, the $A_{\text{Gal}}$ based on five different $\varphi(\vec{x})$ models for the 15 pulsars with significant $D$ measurements. The five $\varphi(\vec{x})$ models are denoted as NT95 (Nice & Taylor 1995), DB98 (Dehnen & Binney 1998), BT08 (Binney & Tremaine 2011), P14 (Piffl et al. 2014), and M17 (McMillan 2017), in this paper. The
<table><caption>**Table 8.** Extrinsic terms of $\hat{P}_s$ and $\hat{P}_b$ for 15 pulsars with significant D.</caption><thead><tr><th>PSR</th><th>A<sub>Shk</sub><br>(pm s<sup>-2</sup>)*</th><th>A<sub>Gal</sub><br>(pm s<sup>-2</sup>)</th><th>P<sub>s</sub><br>(ms)</th><th>P&#770;<sub>Shk</sub><br>(zs s<sup>-1</sup>)*</th><th>P&#770;<sub>Gal</sub><br>(zs s<sup>-1</sup>)</th><th>P&#770;<sub>obs</sub><br>(zs s<sup>-1</sup>)</th><th>P&#770;<sub>int</sub><br>(zs s<sup>-1</sup>)</th><th>P&#770;<sub>Shk</sub><br>(fs s<sup>-1</sup>)</th><th>P&#770;<sub>Gal</sub><br>(fs s<sup>-1</sup>)*</th><th>P&#770;<sub>obs</sub><br>(fs s<sup>-1</sup>)</th><th>P&#770;<sub>int</sub><br>(fs s<sup>-1</sup>)</th><th>P&#770;<sub>b</sub><sup>GW</sup><br>(fs s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0030+0451</td><td>9.1(2)</td><td>-33.0(3.7)</td><td>4.87</td><td>0.148(3)</td><td>-0.54(6)</td><td>10.2</td><td>10.59(6)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J0610-2100</td><td>3.9<sup>+0.8</sup><sub>-0.6</sub> × 10<sup>2</sup></td><td>-9(2)</td><td>3.86</td><td>5.0<sup>+1.0</sup><sub>-0.7</sub></td><td>-0.12(3)</td><td>12.3<sup>a</sup></td><td>7.4(9)</td><td>32<sup>+6</sup><sub>-5</sub></td><td>-0.72(17)</td><td>-70(30)<sup>a</sup></td><td>-101(31)</td><td>~4.6<sup>a</sup></td></tr><tr><td>J0621+1002</td><td>14<sup>+4</sup><sub>-3</sub></td><td>23.8(4.5)</td><td>28.85</td><td>1.3<sup>+0.4</sup><sub>-0.3</sub></td><td>2.3(4)</td><td>47.3</td><td>43.7(6)</td><td>33<sup>+11</sup><sub>-7</sub></td><td>57(11)</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1012+5307</td><td>419<sup>+17</sup><sub>-15</sub></td><td>-23.5(2.4)</td><td>5.26</td><td>7.3(3)</td><td>-0.41(4)</td><td>17.1</td><td>10.2(3)</td><td>73(3)</td><td>-4.1(4)</td><td>61(4)</td><td>-7.9(5.0)</td><td>-13(1)<sup>b</sup></td></tr><tr><td>J1024-0719</td><td>2.8(2) × 10<sup>3</sup></td><td>-40(3)</td><td>5.16</td><td>48(3)</td><td>-0.69(5)</td><td>18.6</td><td>-29(3)**</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1518+4904</td><td>43(1)</td><td>-48.5(3.2)</td><td>40.93</td><td>5.9(2)</td><td>-6.6(4)</td><td>27.2</td><td>27.9(5)</td><td>107(3)</td><td>-120(8)</td><td>2.4(2.2) × 10<sup>2</sup></td><td>2.6(2.2) × 10<sup>2</sup></td><td>~1.2<sup>e</sup></td></tr><tr><td>J1537+1155</td><td>4.4(3) × 10<sup>2</sup></td><td>-42(3)</td><td>37.90</td><td>55.6<sup>+4.0</sup><sub>-3.5</sub></td><td>-5.3(4)</td><td>2422.5</td><td>2372(4)</td><td>53.3<sup>+3.8</sup><sub>-3.3</sub></td><td>-5.1(4)</td><td>-136.6(3)</td><td>-185(4)</td><td>-192.45(6)<sup>c</sup></td></tr><tr><td>J1640++2224</td><td>1.3(1) × 10<sup>2</sup></td><td>-48.5(4.3)</td><td>3.16</td><td>1.4(1)</td><td>-0.51(5)</td><td>2.8</td><td>1.9(1)</td><td>6.8<sup>+0.6</sup><sub>-0.5</sub> × 10<sup>3</sup></td><td>-2.45(22) × 10<sup>3</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1643-1224</td><td>35<sup>+5</sup><sub>-4</sub></td><td>1.0(1.7)</td><td>4.62</td><td>0.53<sup>+0.08</sup><sub>-0.06</sub></td><td>2(3) × 10<sup>-2</sup></td><td>18.5</td><td>17.95(7)</td><td>1.5(2) × 10<sup>3</sup></td><td>41(72)</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1730-2304</td><td>158<sup>+9</sup><sub>-8</sub></td><td>11.7(9)</td><td>8.12</td><td>4.3(2)</td><td>0.32(2)</td><td>20.2</td><td>15.6(2)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1738++0333</td><td>95<sup>+8</sup><sub>-7</sub></td><td>-6.2(1.6)</td><td>5.85</td><td>1.86<sup>+0.16</sup><sub>-0.14</sub></td><td>-0.12(3)</td><td>24.1</td><td>22.4(2)</td><td>9.7<sup>+0.9</sup><sub>-0.7</sub></td><td>-0.64(16)</td><td>-17(3)</td><td>-26.1(3.1)</td><td>-27.7<sup>+1.5</sup><sub>-1.9</sub></td></tr><tr><td>J1853+1303</td><td>17<sup>+3</sup><sub>-2</sub></td><td>-16(5)</td><td>4.09</td><td>0.23<sup>+0.04</sup><sub>-0.03</sub></td><td>-0.22(7)</td><td>8.7</td><td>8.69(8)</td><td>5.6<sup>+0.9</sup><sub>-0.7</sub> × 10<sup>2</sup></td><td>-5(2) × 10<sup>2</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1910+1256</td><td>1.2<sup>+0.4</sup><sub>-0.2</sub> × 10<sup>2</sup></td><td>-35(17)</td><td>4.98</td><td>2.0<sup>+0.6</sup><sub>-0.4</sub></td><td>-0.6(3)</td><td>9.7</td><td>8.3(6)</td><td>2.0<sup>+0.6</sup><sub>-0.4</sub> × 10<sup>3</sup></td><td>-6(3) × 10<sup>2</sup></td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1918-0642</td><td>93<sup>+12</sup><sub>-9</sub></td><td>6.8(1.5)</td><td>7.65</td><td>2.4<sup>+0.3</sup><sub>-0.2</sub></td><td>0.17(4)</td><td>25.7</td><td>23.1(3)</td><td>2.9<sup>+0.4</sup><sub>-0.3</sub> × 10<sup>2</sup></td><td>21(5)</td><td>-</td><td>-</td><td>-</td></tr><tr><td>J1939+2134</td><td>0.37<sup>+0.04</sup><sub>-0.03</sub></td><td>-60(10)</td><td>1.56</td><td>1.9(2) × 10<sup>-3</sup></td><td>-0.31(5)</td><td>105.1</td><td>105.41(5)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Notes. •Unless otherwise specified, we adopted $P_s$, $\hat{P}_{\text{obs}}$, $P_b$, and $\hat{P}_{\text{obs}}$ from the ATNF Pulsar Catalogue².

•Other references: a van der Wateren et al. (2022); b Ding et al. (2020b) and references therein; c Ding et al. (2021a) and references therein; d Freire et al. (2012); e Janssen et al. (2008).

*'zs' meanszeptosecond, which is 10⁻²¹s. Besides, 'pm' and 'fs' represent 10⁻¹²m and 10⁻¹⁵s, respectively.

**The negative $\hat{P}_s$ is the result of the radial acceleration caused by a far-away companion (Bassa et al. 2016; Kaplan et al. 2016).
<table><caption>**Table 9.** Radial accelerations due to Galactic gravitational pull based on different models of Galactic gravitational potential.</caption><thead><tr><th>PSR</th><th>A<sub>Gal</sub><sup>(NT95)</sup><br>(pm s<sup>-2</sup>)</th><th>A<sub>Gal</sub><sup>(DB98)</sup><br>(pm s<sup>-2</sup>)</th><th>A<sub>Gal</sub><sup>(BT08)</sup><br>(pm s<sup>-2</sup>)</th><th>A<sub>Gal</sub><sup>(P14)</sup><br>(pm s<sup>-2</sup>)</th><th>A<sub>Gal</sub><sup>(M17)</sup><br>(pm s<sup>-2</sup>)</th><th>A<sub>Gal</sub><sup>(GR)</sup><sup>*</sup><br>(pm s<sup>-2</sup>)</th></tr></thead><tbody><tr><td>J0030+0451</td><td>-29(3)</td><td>¹-37.0(2)</td><td>¹-27.3(2)</td><td>-35.3(2)</td><td>-32.5(3)</td><td>-</td></tr><tr><td>J0610-2100</td><td>¹-12(1)</td><td>-8.6<sup>+1.2</sup><sub>-0.8</sub></td><td>¹-6.0<sup>+1.0</sup><sub>-0.5</sub></td><td>-10.9<sup>+0.7</sup><sub>-0.3</sub></td><td>-8.8<sup>+1.0</sup><sub>-0.4</sub></td><td>-</td></tr><tr><td>J0621+1002</td><td>24(4)</td><td>23(5)</td><td>24(5)</td><td>24(5)</td><td>25(5)</td><td>-</td></tr><tr><td>J1012+5307</td><td>!!!-32.0(6)</td><td>-24.0(2)</td><td>¹-19.44(9)</td><td>-24.80(9)</td><td>-25.80(6)</td><td>05(29)</td></tr><tr><td>J1024-0719</td><td>¹-45.1(9)</td><td>-38.5(6)</td><td>¹-35.6(9)</td><td>-42(1)</td><td>-43(1)</td><td>-</td></tr><tr><td>J1518+4904</td><td>-47.5(5)</td><td>-48.1(5)</td><td>¹-44.8(7)</td><td>-50.4(7)</td><td>¹-51.9(7)</td><td>-</td></tr><tr><td>J1537+1155</td><td>!!!-29(1)</td><td>-42(1)</td><td>-39(2)</td><td>-43(2)</td><td>-45(2)</td><td>21<sup>+28</sup><sub>-31</sub></td></tr><tr><td>J1640+2224</td><td>!!!-33(1)</td><td>-46(3)</td><td>-45(3)</td><td>-50(3)</td><td>-52(4)</td><td>-</td></tr><tr><td>J1643-1224</td><td>¹¹0(3)</td><td>¹-1.2<sup>+0.8</sup><sub>-0.6</sub></td><td>¹ 3.2<sup>+0.7</sup><sub>-0.5</sub></td><td>0.6<sup>+0.7</sup><sub>-0.6</sub></td><td>1.3<sup>+0.6</sup><sub>-0.4</sub></td><td>-</td></tr><tr><td>J1730-2304</td><td>13.2(8)</td><td>10.8(6)</td><td>12.1(6)</td><td>11.5(6)</td><td>12.5(7)</td><td>-</td></tr><tr><td>J1738+0333</td><td>!!!!!!10.1(1.8)</td><td>-6.4<sup>+0.4</sup><sub>-0.6</sub></td><td>-4.2<sup>+0.8</sup><sub>-1.1</sub></td><td>-7.5<sup>+0.7</sup><sub>-1.0</sub></td><td>-6.9<sup>+0.8</sup><sub>-1.2</sub></td><td>9(35)</td></tr><tr><td>J1853+1303</td><td>-13(3)</td><td>-13<sup>+3</sup><sub>-4</sub></td><td>-13<sup>+3</sup><sub>-5</sub></td><td>-19<sup>+4</sup><sub>-5</sub></td><td>-16<sup>+4</sup><sub>-5</sub></td><td>-</td></tr><tr><td>J1910+1256</td><td>-35(13)</td><td>-29<sup>+8</sup><sub>-16</sub></td><td>-31<sup>+10</sup><sub>-18</sub></td><td>-42<sup>+12</sup><sub>-21</sub></td><td>-36<sup>+10</sup><sub>-20</sub></td><td>-</td></tr><tr><td>J1918-0642</td><td>¹¹4(2)</td><td>5.9(5)</td><td>¹'8.7(5)</td><td>¹'5.0(2)</td><td>7.4(3)</td><td>-</td></tr><tr><td>J1939+2134</td><td>-64(8)</td><td>-53(8)</td><td>-56(8)</td><td>-67(9)</td><td>-63(9)</td><td>-</td></tr></tbody></table>

Notes. ●NT95, DB98, BT08, P14, and M17 refer to five different $\varphi(\vec{x})$ models (see Section 7.2 for the references).

●The '!'s indicate, in the same way as Table 6, the significance of the offset between the $A_{\text{Gal}}$ in Table 8 and that of each $\varphi(\vec{x})$ model.

*$A_{\text{Gal}}^{(\text{GR})} = (\dot{P}_{\text{b}}^{\text{obs}} - \dot{P}_{\text{b}}^{\text{GW}} - \dot{P}_{\text{b}}^{\text{Shk}})c/P_{\text{b}}$ is based on the assumption that GR is correct.

results obtained with NT95, which uses a simple analytical approach, are frequently discrepant compared to the other four $\varphi(\vec{x})$ models. Accordingly, and following Ding et al. (2021a), we exclude it and use the remaining four models to derive the estimate for $A_{\text{Gal}}$ and its uncertainty, which we present in Table 8 (along with $\dot{P}_{\text{b}}^{\text{Gal}}$ and $\dot{P}_{\text{s}}^{\text{Gal}}$).

Incorporating the $\dot{P}_{\text{s}}^{\text{Shk}}$ derived in Section 7.1, we calculated the intrinsic spin period derivative $\dot{P}_{\text{s}}^{\text{int}} = \dot{P}_{\text{s}}^{\text{obs}} - \dot{P}_{\text{s}}^{\text{Shk}} - \dot{P}_{\text{s}}^{\text{Gal}}$. We note that the negative $\dot{P}_{\text{s}}^{\text{int}}$ of PSR J1024-0719 is probably the consequence of radial acceleration induced by a putative companion in an extremely wide orbit with PSR J1024-0719 (Bassa et al. 2016; Kaplan et al. 2016, also see Section 5.2). In addition to $\dot{P}_{\text{s}}^{\text{int}}$, $\dot{P}_{\text{b}}^{\text{int}} = \dot{P}_{\text{b}}^{\text{obs}} - \dot{P}_{\text{b}}^{\text{Shk}} - \dot{P}_{\text{b}}^{\text{Gal}}$ are estimated for four pulsar systems with reported $\dot{P}_{\text{b}}^{\text{obs}}$. The improved PSR J1738+0333 parallax as well as the re-assessed PSR J1012+5307 parallax calls for an update to the constraint on alternative theories of gravity (e.g. Freire et al. 2012; Zhu et al. 2019; Ding et al. 2020b), which is discussed in Section 7.3.

While performing the $A_{\text{Gal}}$ analysis, we found an error in the code that had been used to implement the calculation of equation (11) for the Ding et al. (2021a) work (which, to be clear, is not an error in the GalPot<fnref n="10" /> package that provides the $\varphi(\vec{x})$ models). Therefore, we note that the $\dot{P}_{\text{b}}^{\text{Gal}}$ of PSR J1537+1155 in Table 8 is a correction to the Ding et al. (2021a) counterpart. Further discussions on PSR J1537+1155 can be found in Section 8.3.

Last but not the least, assuming GR is correct, the approach taken above can be inverted to infer $A_{\text{Gal}}^{(\text{GR})} = (\dot{P}_{\text{b}}^{\text{obs}} - \dot{P}_{\text{b}}^{\text{GW}} - \dot{P}_{\text{b}}^{\text{Shk}})c/P_{\text{b}}$, which can be used to constrain Galactic parameters for the local environment (of the Solar system) (Bovy 2020), or probe the Galactic dark matter distribution in the long run (Phillips et al. 2021). The $A_{\text{Gal}}^{(\text{GR})}$ for the three viable pulsars are listed in Table 9.

## 7.3 New constraints on alternative theories of gravity

In the GR framework, the excess orbital decay $\dot{P}_{\text{b}}^{\text{ex}} = \dot{P}_{\text{b}}^{\text{int}} - \dot{P}_{\text{b}}^{\text{GW}}$ is expected to agree with zero. However, some alternative theories of gravity expect otherwise due to their predictions of non-zero dipole gravitational radiation and time-varying Newton's gravitational constant $G$. Both phenomena are prohibited by GR. Namely, in GR, the dipole gravitational radiation coupling constant $\kappa_D = 0$, and $\dot{G}/G = 0$. The large asymmetry of gravitational binding energy of pulsar-WD systems makes them ideal testbeds for dipole gravitational emissions (e.g. Eardley 1975). In an effort to test (and possibly eliminate) alternative theories of gravity, increasingly tight constraints on $\kappa_D$ and $\dot{G}/G$ have been placed using multiple pulsar-WD systems (Deller et al. 2008; Freire et al. 2012; Zhu et al. 2019; Ding et al. 2020b).

With the reassessed astrometric results of PSR J1012+5307, the $\dot{P}_{\text{b}}^{\text{ex}}$ of PSR J1012+5307 changes from $10.6 \pm 6.1 \text{ fs s}^{-1}$ in Ding et al. (2020b) to $5.1 \pm 5.1 \text{ fs s}^{-1}$. This change is mainly caused by three reasons: (1) priors are placed on the proper motion during inference in this work (but not in Ding et al. 2020b); (2) a Bayesian framework is applied in this work (while Ding et al. 2020b reported bootstrap results); (3) this work adopts PDF medians as the estimates (while Ding et al. 2020b used PDF modes). Though barely affecting this work (see Fig. 4), the choice between PDF mode and median makes a difference to Ding et al. (2020b) given that their parallax PDF is more skewed (see fig. 4 of Ding et al. 2020b). After employing the new VLBI + timing distance, the $\dot{P}_{\text{b}}^{\text{ex}}$ of PSR J1738+0333 has shifted from $2.0 \pm 3.7$ (Freire et al. 2012) to $1.6 \pm 3.5 \text{ fs s}^{-1}$. More discussions on PSR J1738+0333 can be found in Section 8.8.

With the new $\dot{P}_{\text{b}}^{\text{ex}}$ of PSR J1012+5307 and PSR J1738+0333, we updated the constraints on $\kappa_D$ and $\dot{G}/G$ in exactly the same way as Ding et al. (2020b). The prerequisites of this inference are reproduced in Table 10, where the two underlined $\dot{P}_{\text{b}}^{\text{ex}}$ are the only difference

<footnote n="10">¹⁰https://github.com/PaulMcMillan-Astro/GalPot</footnote>
<table><caption>**Table 10.** Excess orbital decay $\dot{P}_b^{\text{ex}} = \dot{P}_b^{\text{obs}} - \dot{P}_b^{\text{Shk}} - \dot{P}_b^{\text{Gal}}$ and other pre-requisites for constraining $\dot{G}/G$ and $\kappa_D$.</caption><thead><tr><th>PSR</th><th><i>P</i><sub>b</sub><br>(d)</th><th><i>P</i>'<sub>b</sub><sup>ex</sup><br>(fs s<sup>-1</sup>)</th><th><i>m</i><sub>p</sub><br>(M<sub>☉</sub>)</th><th><i>m</i><sub>c</sub><br>(M<sub>☉</sub>)</th><th><i>q</i></th></tr></thead><tbody><tr><td>J0437-4715</td><td>5.74</td><td>12(32)</td><td>1.44(7)</td><td>0.224(7)</td><td>-</td></tr><tr><td>J1012+5307</td><td>0.60</td><td><u>5.1(5.1)</u></td><td>-</td><td>0.174(11)</td><td>10.44(11)</td></tr><tr><td>J1713+0747</td><td>67.83</td><td>30(150)</td><td>1.33(10)</td><td>0.290(11)</td><td>-</td></tr><tr><td>J1738+0333</td><td>0.35</td><td><u>1.6(3.5)</u></td><td>1.46(6)</td><td>-</td><td>8.1(2)</td></tr></tbody></table>

Note. •*m*<sub>p</sub>, *m*<sub>c</sub> and *q* stand for, respectively, pulsar mass, companion mass and mass ratio (i.e. *m*<sub>p</sub>/*m*<sub>c</sub>). See Ding et al. (2020b) for their references.

from the Table 6 of Ding et al. (2020b). We obtained

$$ 
\begin{aligned}
\dot{G}/G &= -1.6^{+5.3}_{-4.8} \times 10^{-13} \text{ yr}^{-1}, \\
\kappa_D &= -1.1^{+2.4}_{-0.9} \times 10^{-4}.
\end{aligned} 
\qquad (12) $$

Compared to Ding et al. (2020b), $\kappa_D$ becomes more consistent with zero, while the new uncertainties of $\kappa_D$ and $\dot{G}/G$ remain at the same level.

## 8 INDIVIDUAL PULSARS

In this section, we discuss the impacts of the new astrometric measurements (particularly the new distances) on the scientific studies around individual pulsars. Accordingly, special attention is paid to the cases where there is no published timing parallax $\sigma^{(\text{Ti})}$. In addition, we also look into the two pulsars (i.e. PSR J1721-2457 and PSR J1824-2452A) that have $\sigma'$ consistent with zero in an effort to understand the causes of parallax non-detection.

### 8.1 PSR J0610-2100

PSR J0610-2100 is the third black widow pulsar discovered (Burgay et al. 2006), which is in a 7-hr orbit with an extremely low-mass ($\approx 0.02$ M<sub>☉</sub>, Pallanca et al. 2012) star. Adopting a distance of around 2.2 kpc, van der Wateren et al. (2022) obtained a $\gamma$-ray emission efficiency $\eta_\gamma \equiv 4\pi F_\gamma D^2/\dot{E}^{\text{int}}$ in the range of 0.5–3.7, where $\dot{E}^{\text{int}}$ and $F_\gamma$ are, respectively, the intrinsic NS spin-down power and the $\gamma$-ray flux above 100 MeV.

In addition, van der Wateren et al. (2022) estimated a mass function

$$ f(m_{\text{p}}, q) = m_{\text{p}} \frac{\sin^3 i}{q(q+1)^2} = \frac{4\pi^2 a_1^3}{G P_b^2} \qquad (13) $$

of $5.2 \times 10^{-6}$ M<sub>☉</sub> for the PSR J0610-2100 system (where $q \equiv m_{\text{p}}/m_{\text{c}}$). Besides, they determined the irradiation temperature (of the companion) $T_{\text{irr}} = 2820 \pm 190$ K as well as the projected orbital semi-major axis $a_1 = 7.3 \times 10^{-3}$ lt-s. Combining these three estimates, we calculated the heating luminosity

$$ 
\begin{aligned}
L_{\text{irr}} &\equiv 4\pi \left[ \frac{a_1(1+q)}{\sin i} \right]^2 \sigma_{\text{SB}} T_{\text{irr}}^4 \\
&\approx 4\pi a_1^2 \left[ \frac{m_{\text{p}}}{f(m_{\text{p}}, q)} \right]^{2/3} \sigma_{\text{SB}} T_{\text{irr}}^4 \\
&\sim 9.1 \times 10^{32} \left( \frac{m_{\text{p}}}{1.4 \text{ M}_{\odot}} \right)^{2/3} \text{erg s}^{-1},
\end{aligned} 
\qquad (14) $$

where $\sigma_{\text{SB}}$ represents the Stefan-Boltzmann constant.

Our new distance $D = 1.5^{+0.3}_{-0.2}$ kpc to PSR J0610-2100 is less than half the DM-based distances (see Table 6), and significantly below

that assumed by van der Wateren et al. (2022). Assuming a NS moment of inertia $I_{\text{NS}} = 10^{45}$ g cm<sup>2</sup>, the $\dot{P}_s^{\text{int}}$ of PSR J0610-2100 (see Table 8) corresponds to an intrinsic spin-down power

$$ \dot{E}^{\text{int}} \equiv 4\pi^2 I_{\text{NS}} \frac{\dot{P}_s^{\text{int}}}{P_s^3} \qquad (15) $$

of $(5.1 \pm 0.5) \times 10^{33}$ erg s<sup>-1</sup>, which is roughly twice as large as the $\dot{E}^{\text{int}}$ range calculated by van der Wateren et al. (2022). In conjunction with a smaller $\gamma$-ray luminosity $L_\gamma = 4\pi F_\gamma D^2$ (due to closer distance), the $\dot{E}^{\text{int}}$ reduced $\eta_\gamma$ to around 0.37 (from 0.5 < $\eta_\gamma$ < 3.7 estimated by van der Wateren et al. 2022), disfavoring unusually high-$\gamma$-ray beaming towards us. Moreover, the heating efficiency $\epsilon_T$ drops to ~0.17 (from 0.15 < $\epsilon_T$ < 0.77 evaluated by van der Wateren et al. 2022), disfavoring the scenario where the NS radiation is strongly beamed towards the companion.

### 8.1.1 On the DM discrepancy

In Section 6.1.1, we noted that our VLBI parallax-derived distance and the DM model-inferred distance to this pulsar differed substantially. Specifically, PSR J0610-2100 has a measured DM = $60.7$ pc cm<sup>-3</sup> while the NE2001 model predicts $27.5$ pc cm<sup>-3</sup> for a line of sight of length 1.5 kpc. We attribute this discrepancy to thermal plasma or ‘free electrons’ along the line of sight that is not captured fully by a ‘clump’ in the NE2001 model. The NE2001 model includes this ‘clump’ to describe the effects due to the Mon R2 H II region, centred at a Galactic longitude and latitude of (214°, -12.6°), located at an approximate distance of ~0.9 kpc (Herbst 1975). However, the WHAM survey shows considerable $H\alpha$ in this direction, extending over tens of degrees. Lines of sight close to the pulsar show changes in the $H\alpha$ intensity by factors of two, but an approximate value toward the pulsar is roughly 13 Rayleighs, equivalent to an emission measure EM = $29$ pc cm<sup>-6</sup> (for a temperature $T = 8000$ K). Using standard expressions, as provided in the NE2001 model, to convert EM to DM, there is sufficient $H\alpha$ intensity along the line of sight to account for the excess DM that we infer from the difference between our parallax-derived distance and the NE2001 model distance.

## 8.2 PSR J1518+4904

The 41-ms PSR J1518+4904, discovered by Sayer, Nice & Taylor (1997), is one of the only two DNSs in the current sample. According to Janssen et al. (2008), the non-detection of Shapiro delay effects suggests $\sin i \le 0.73$ at 99 per cent confidence level. Accordingly, we adopted 0.73 as the upper limit of $\sin i$, and carried out 8-parameter Bayesian inference, which led to a bi-modal posterior PDF on $i'$ and a multimodal PDF on $\Omega'_{\text{asc}}$ (see the online corner plot⁵). The predominant constraints on both $i'$ and $\Omega'_{\text{asc}}$ come from the $\dot{a}_1$ measurement (Janssen et al. 2008 or see Table 4). Though there are three major likelihood peaks for the $\Omega'_{\text{asc}}$, two of them gather around 171°, making the PDF relatively more concentrated. When a much more precise $\dot{a}_1$ measurement is reached with new timing observations, the existing VLBI data will likely place useful constraints on $i'$ and $\Omega'_{\text{asc}}$. So will additional VLBI observations.

In addition to $i'$ and $\Omega'_{\text{asc}}$, the 8-parameter Bayesian inference also renders a $40$ σ parallax $\sigma'$, which becomes the most significant parallax achieved for a DNS. In contrast, to detect a timing parallax $\sigma^{(\text{Ti})}$ for PSR J1518+4904 would take $\gtrsim 600$ yr (Janssen et al. 2008) due to its relatively high-ecliptic latitude of 63°.
## 8.3 PSR J1537+1155

PSR J1537+1155, also known as PSR B1534+12, is the second discovered DNS (Wolszczan 1991). The DNS displays an exceptionally high-proper motion amongst all Galactic DNSs (see Table 3 of Tauris et al. 2017), leading to an unusually large Shklovskii contribution to observed timing parameters. Therefore, precise astrometry of the DNS plays an essential role in its orbital decay test of GR. The most precise astrometric parameters of PSR J1537+1155 are provided by Ding et al. (2021a) based on the same data set used in this work, which result in $\dot{P}_b^{\text{Shk}} = 53 \pm 4$ fs s$^{-1}$. Subsequently, Ding et al. (2021a) estimated $\dot{P}_b^{\text{Gal}} = -1.9 \pm 0.2$ fs s$^{-1}$, and concluded $\dot{P}_b^{\text{int}}/\dot{P}_b^{\text{GW}} = 0.977 \pm 0.020$.

In this work, we inferred $\eta_{\text{EFAC}}$ on top of the canonical astrometric parameters, which is the only difference from the Bayesian method of Ding et al. (2021a). Despite this difference, our astrometric results of PSR J1537+1155 remain almost the same as Ding et al. (2021a). So is our re-derived $\dot{P}_b^{\text{Shk}} = 53.3_{-3.3}^{+3.8}$ fs s$^{-1}$. However, as is mentioned in Section 7.2, the $\dot{P}_b^{\text{Gal}}$ estimated by Ding et al. (2020a) is incorrect due to a coding error. After correction, $\dot{P}_b^{\text{Gal}}$ drops to $-5.1 \pm 0.4$ fs s$^{-1}$ (see Table 8). Consequently, we obtained $\dot{P}_b^{\text{int}}/\dot{P}_b^{\text{GW}} = 0.96 \pm 0.02$.

As Ding et al. (2021a) have pointed out, the limiting factor of the GR orbital decay test using PSR J1537+1155 remains the distance precision, which generally improves as $t^{-1/2}$ with additional observations, but can be accelerated if more sensitive instrumentation can be deployed. On the other hand, the extremely high-braking index of 157 (two orders higher than the normal level) calculated from the rotational frequency $v_s \equiv 1/P_s$, its first derivative $\dot{v}_s$, and its second derivative $\ddot{v}_s$ (Fonseca et al. 2014) indicate likely timing noise contributions that may affect the observed orbital period derivative to some degree. This will be clarified with continued timing observations and refined timing analysis.

## 8.4 PSR J1640+2224

PSR J1640+2224 is a 3.2-ms MSP (Foster et al. 1995) in a wide ($P_b = 175$ d) orbit with a WD companion (Lundgren et al. 1996). The MSPSR$\pi$ results for PSR J1640+2224 have been determined using bootstrap and published in Vigeland et al. (2018), which are highly consistent with our re-assessed quasi-VLBI-only results (see Table 2 of Vigeland et al. 2018 and Table 5), and also agree with the VLBI + timing ones (see Table 6). Our 8-parameter Bayesian inference renders a 1D histogram of $\Omega'_{\text{asc}}$ with four likelihood components at 0°, 140°, 200°, and 320°, which is predominantly shaped by the prior on $\dot{a}_1$ from pulsar timing (see Table 4).

## 8.5 PSR J1643−1224

PSR J1643−1224 is a 4.6-ms pulsar in a 147-d orbit with a WD companion (Lorimer et al. 1995). As a result of multipath propagation due to inhomogeneities in the ionised interstellar medium (IISM), the pulse profiles of PSR J1643−1224 are temporally broadened (e.g. Lentati et al. 2017). As the Earth-to-pulsar sightline moves through inhomogeneous scattering ‘screen(s)’ (in the IISM), the temporal broadening $\tau_{sc}$ varies with time; at 1 GHz, $\tau_{sc}$ fluctuates up and down by $\lesssim 5$ µs on a yearly time-scale (Lentati et al. 2017). Meanwhile, the moving scattering ‘screen(s)’ would also change the radio brightness of the pulsar. This effect, as known as pulsar scintillation, is used to constrain the properties of both PSR J1643−1224, and

the scattering screen(s) between the Earth and the pulsar (Mall et al. 2022). The scintillation of PSR J1643−1224 has previously been modelled with both isotropic and anisotropic screens (Mall et al. 2022). The isotropic model renders a pulsar distance $D = 1.0 \pm 0.3$ kpc and locates the main scattering screen at $D_{sc} = 0.21 \pm 0.02$ kpc; in comparison, the anisotropic model yields a pulsar distance $D = 1.2 \pm 0.3$ kpc, and necessitates a secondary scattering screen $0.34 \pm 0.09$ kpc away (from the Earth) in addition to a main scattering screen at $0.13 \pm 0.02$ kpc distance (Mall et al. 2022). On the other hand, the HII region Sh 2-27 in front of PSR J1643−1224 is suspected to be associated with the main scattering screen of the pulsar. This postulated association is strengthened by the agreement between the distance to the main scattering screen (based on the two-screen anisotropic model, Mall et al. 2022) and the distance to the HII region (i.e. $112 \pm 3$ pc, Ocker, Cordes & Chatterjee 2020).

### 8.5.1 Independent check on the postulated association between the HII region Sh 2-27 and the main scattering screen

Besides the pulse broadening effect, the scattering by the IISM would lead to apparent angular broadening of the pulsar, which has been detected with VLBI at $\gtrsim 8$ GHz (e.g. Bower et al. 2014). By the method outlined in Appendix A of Ding et al. (2020a), we measured a semi-angular-broadened size $\theta_{sc} = 3.65 \pm 0.43$ mas for PSR J1643−1224, which is the only significant $\theta_{sc}$ determination in the MSPSR$\pi$ catalogue. Likewise, the secondary in-beam calibrator of PSR J1643−1224 is also scatter-broadened, which may likely introduce additional astrometric uncertainties (see more explanation in Section 8.6).

As both pulse broadening and angular broadening are caused by the IISM deflection, $\theta_{sc}$, $\tau_{sc}$, the pulsar distance $D$, and the distance(s) $D_{sc}$ to the scattering screen(s) are geometrically related. Assuming there is one dominant thin scattering screen, we make use of following approximate relation

$$ \frac{\theta_{sc}^2}{2\pi\tau_{sc}} = \frac{1}{D_{sc}} - \frac{1}{D} \quad (\text{when } \theta_{sc} \lesssim 1^{\circ}), \qquad (16) $$

where c stands for the speed of light in vacuum.

To estimate the unknown $\tau_{sc}$ at our observing frequency of $\sim 1.55$ GHz, we used the data spanning MJD 54900–57500 from the PPTA second data release (Kerr et al. 2020). We analysed the dynamic spectra of observations centred around 3.1 GHz and recorded with the PDFB4 processor, using the scintools<fnref n="11" /> package (Reardon et al. 2020). A model was fit to the auto-correlation function of each dynamic spectrum, which has an exponential decay with frequency (Reardon et al. 2019). The characteristic scintillation scale (in frequency) $\Delta v_d$ is related to the scattering time-scale with $\tau_{sc} = 1/(2\pi\Delta v_d)$. We found the mean temporal broadening $\tau_{sc} = 103 \pm 25$ ns at 3.1 GHz, with fluctuations of $\lesssim 60$ ns (see Fig. 6). To convert this $\tau_{sc}$ to our observing frequency 1.55 GHz, we compare the maximum degree (i.e. 60 ns) of fluctuations at 3.1 GHz to that (i.e. 5 µs, Lentati et al. 2017) at 1 GHz, and acquired an indicative scaling relation

$$ \tau_{sc} \propto v^{-3.9}, \qquad (17) $$

where v is the observing frequency. This relation reasonably agrees with the scaling relation $\tau_{sc} \propto v^{-11/3}$ associated with the Kolmogorov

<footnote n="11">¹¹https://github.com/danielreardon/scintools</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Temporal broadening $\tau_{sc}$ of PSR J1643−1224 at 3.1 GHz. The solid red line and the dashed red line show the mean temporal broadening and a 68 per cent confidence interval, respectively.</figcaption></figure>

turbulence (e.g. Armstrong, Rickett & Spangler 1995). With the indicative scaling relation, we calculated $\tau_{sc} = 1.54 \pm 0.37$ µs. It is timely to note that $\theta_{sc}^2/\tau_{sc}$ (on the left side of equation 16) is frequency-independent. By combining equations (16) and (17), we come to another equivalent indicative scaling relation

$$ \theta_{sc} \propto \nu^{-1.95}. \qquad (18) $$

Substituting $\tau_{sc} = 1.54 \pm 0.37$ µs, $\theta_{sc} = 3.65 \pm 0.43$ mas and $D = 0.95_{-0.11}^{+0.15}$ kpc into equation (16), we obtained $D_{sc} = 86_{-24}^{+30}$ pc, where the uncertainty is derived with a Monte-Carlo simulation. This $D_{sc}$ is consistent with the distance to the HII region Sh 2-27 (Ocker et al. 2020), hence independently supporting the association between the HII region and the main scattering screen of PSR J1643−1224.

### 8.5.2 Probing scintillation models

Apart from the above check on the connection between the HII region Sh 2-27 and the main scattering screen, the angular broadening of PSR J1643−1224 also promises a test of the aforementioned isotropic scintillation model proposed by Mall et al. (2022). To do so, we changed the pulsar distance to the one inferred with the model (i.e. $1.0 \pm 0.3$ kpc). With this change, we derive $D_{sc} = 86_{-24}^{+30}$ pc, which disagrees with $0.21 \pm 0.02$ kpc based on the isotropic model. To investigate the impact of a different scaling relation $\tau_{sc} \propto \nu^{-\alpha_{sc}}$, we inferred $\tau_{sc} = 4.3$ µs with both $D$ and $D_{sc}$ based on the isotropic model (Mall et al. 2022), which corresponds to an unreasonably large $\alpha_{sc} = 5.4$. Hence, we conclude that our $\theta_{sc}$ and $\tau_{sc}$ cannot easily reconcile with the one-screen isotropic model proposed by Mall et al. (2022).

Fundamentally, the irreconcilability implies a one-screen model might be incapable of describing both scintillation and angular broadening effects of PSR J1643−1224. In principle, it is possible to analyse a multiscreen model with a $\theta_{sc}(t)$ series (at various time $t$), and its associated $\tau_{sc}(t)$ instead of only using their mean values. However, this analysis is not feasible for this work, as $\tau_{sc}$ and $\theta_{sc}$ were not measured on the same days. None the less, we can still investigate whether our observations of PSR J1643−1224 can reconcile with the scintillation observations (Mall et al. 2022) in the context of a two-screen model.

In the scenario of two thin scattering screens, we derived the more complicated relation

$$ \begin{cases} 2c\tau_{sc} = k_1\beta_{sc}^2 + k_2\beta_{sc}\theta_{sc} + k_3\theta_{sc}^2 & (\theta_{sc} \lesssim 1^{\circ} \text{ and } \beta_{sc} \lesssim 1^{\circ}) \\ k_1 = \frac{(D - D_{sc2})(D - D_{sc1})}{D_{sc2} - D_{sc1}} \\ k_2 = -\frac{2D_{sc1}(D - D_{sc2})}{D_{sc2} - D_{sc1}} \\ k_3 = \frac{D_{sc1}D_{sc2}}{D_{sc2} - D_{sc1}}, \end{cases} \qquad (19) $$

where $D_{sc1}$ and $D_{sc2}$ are the distance to the first and the second scattering screen, respectively; $\beta_{sc}$ represents a half of the opening angle of the second scattering screen (closer to the pulsar) as seen from the pulsar. As a side note, equation (16) can be considered a special case (i.e. $D = D_{sc2}$) of equation (19). In equation (19), all parameters except $\beta_{sc}$ are known, either determined with the anisotropic two-screen model or obtained in this work. Hence, it is feasible to constrain the geometric parameter $\beta_{sc}$ with the known parameters.

However, equation (19) can yield unphysical solutions (i.e. $\beta_{sc} > 0$). We applied the simple condition

$$ \frac{\theta_{sc}^2}{2c\tau_{sc}} \le \frac{1}{D_{sc1}} - \frac{1}{D} \qquad (20) $$

to ensure equation (19) gives physical solutions of $\beta_{sc}$. This equation is equivalent to $D_{sc1} \le D_{sc}$, where $D_{sc}$ corresponds to the one-screen solution of equation (16). This is because $D_{sc1} > D_{sc}$ would always lead to longer routes, thus exceeding the $\tau_{sc}$ budget. It is important to note that equation (20) is valid for a model with any number of scattering screens. Hence, we recommend to use equation (20) in scintillation model inferences as a prior condition, to cater for the constraints imposed by $\theta_{sc}$ and $\tau_{sc}$ (and thereby truncate the parameter space of a scintillation model).

To test the anisotropic two-screen model (Mall et al. 2022) with our $\theta_{sc}$ and $\tau_{sc}$, we calculated $D_{sc} = 89_{-26}^{+33}$ pc with the pulsar distance (i.e. $D = 1.2 \pm 0.3$ kpc) based on the anisotropic two-screen model. This $D_{sc}$ is consistent with $D_{sc1} = 129 \pm 15$ pc (Mall et al. 2022) (therefore not ruling out $D_{sc1} < D_{sc}$). That is to say, our $\theta_{sc}$ and $\tau_{sc}$ measurements can loosely reconcile with the anisotropic two-screen model (Mall et al. 2022). In comparison, we reiterate our finding that a one-screen model is difficult to describe both scintillation, and angular broadening effects of PSR J1643−1224.

## 8.6 PSR J1721−2457

PSR J1721−2457 is a 3.5-ms solitary MSP discovered at intermediate Galactic latitudes (Edwards & Bailes 2001). The main secondary phase calibrator of PSR J1721−2457 (and indeed, all the sources near it on the plane of the sky) is heavily resolved due to IISM scattering, which leads to non-detections on the longest baselines and a lack of calibration solutions for some antennas, reducing the spatial resolution of the VLBI observations. The non-uniform IISM distribution also leads to refractive image wander as the line-of-sight to the pulsar changes (e.g. Kramer et al. 2021), which is most pronounced for heavily scatter-broadened sources such as the calibrator for PSR J1721−2457. In conjunction with the lower spatial resolutions, which reduces positional precision to begin with, this additional noise term likely results in the parallax non-detection (see Table 5).
## 8.7 PSR J1730−2304

PSR J1730−2304 is a solitary pulsar spinning at $P_s = 8.1$ ms (Lorimer et al. 1995). Being so far the least-energetic (in terms of $\dot{E}^{\text{int}}$) γ-ray pulsar (Guillemot et al. 2016), the pulsar plays a key role in exploring the death line of NS high-energy radiation. Substituting $\dot{P}_s^{\text{int}}$ and $P_s$ of equation (15) with values listed in Table 8, we substantially refined the $\dot{E}^{\text{int}}$ death line (of all γ-ray-emitting pulsars) to

$$ \dot{E}_{\text{death}} \le \dot{E}_{\text{J1730}}^{\text{int}} = (1.15 \pm 0.01) \times 10^{33} \left( \frac{I_{\text{NS}}}{10^{45} \text{ g cm}^2} \right) \times \text{erg s}^{-1}, \quad (21) $$

which is consistent with (but on the higher side of) the previous estimate ($8.4 \pm 2.2 \times 10^{32}$ erg s$^{-1}$) by Guillemot et al. (2016) (assuming the same $I_{\text{NS}}$). On the other hand, we updated the γ-ray luminosity (above 100 MeV) to $L_{\gamma} = 4\pi D^2 F_{\gamma} = (3.1 \pm 1.6) \times 10^{32}$ erg s$^{-1}$, where the precision is limited by the less precise $F_{\gamma}$ (Guillemot et al. 2016). Accordingly, we obtained $\eta_{\gamma} = 0.27 \pm 0.14$.

## 8.8 PSR J1738+0333

PSR J1738+0333, discovered from a 1.4-GHz high-Galactic-latitude survey with the 64-m Parkes radio telescope (Jacoby et al. 2009) is a 5.85-ms pulsar in a 8.5-hr orbit with a WD companion. Thanks to the relatively short $P_b$, the WD-pulsar system plays a leading role in constraining alternative gravitational theories that predict dipole gravitational radiation (Freire et al. 2012; Zhu et al. 2015).

Our VLBI-only $\sigma$ is 2.3 $\sigma$ away from the most precise $\sigma^{(\text{Ti})}$ (see Table 5 and 6). After adopting timing priors, $\sigma' = 0.589 \pm 0.046$ mas becomes closer to $\sigma^{(\text{Ti})} = 0.68 \pm 0.05$ mas (Freire et al. 2012), meaning that $\dot{P}_b^{\text{Shk}}$ is only 1.2 times larger than the previous estimate. On the other hand, our re-assessed $\dot{P}_b^{\text{Gal}}$, based on more realistic $\phi(\bar{x})$ models, is smaller than that estimated with the NT95 $\phi(\bar{x})$ model (Freire et al. 2012). Combining the unchanged $\dot{P}_b^{\text{obs}} = -17 \pm 3$ fs s$^{-1}$ the re-derived $\dot{P}_b^{\text{Int}} = -26.1 \pm 3.1$ fs s$^{-1}$ is almost the same as the previous estimate, as the change of $\dot{P}_b^{\text{Gal}}$ happens to nearly cancels out that of $\dot{P}_b^{\text{Shk}}$.

Future pulsar timing or VLBI investigation into the discrepancy between $\sigma^{(\text{Ti})}$ (Freire et al. 2012) and $\sigma$ is merited by the importance of the pulsar system. Specifically, if the true parallax turns out to be around 0.5 mas, it would not only mean that $\dot{P}_b^{\text{Shk}}$ is 1.4 times higher than the estimate by Freire et al. (2012), but also suggest the WD radius $R_{\text{WD}}$ to be 1.4 larger (as $R_{\text{WD}} \propto D$ according to Equation 1 of Antoniadis et al. 2012). A higher $R_{\text{WD}}$ would lead to lighter WD and NS (as the mass ratio is well determined), thus smaller $\dot{P}_{\text{GW}}$.

## 8.9 PSR J1824−2452A

PSR J1824−2452A is a 3-ms solitary pulsar discovered in the Globular cluster M28 (NGC 6626) (Lyne et al. 1987). The calibration configuration for this pulsar was suboptimal as the best in-beam phase calibrator for PSR J1824−2452A was both resolved and faint (3.3 mJy, see Table 1), leading to noisy solutions, especially on the longest baselines. This is likely responsible for the parallax non-detection (see Table 5), and indicates that higher sensitivity to enable improved calibration solutions would be advantageous in any future VLBI campaign.

The proper motion of M28 is estimated to be $\mu_{\alpha}^{\text{M28}} = -0.28 \pm 0.03$ mas yr$^{-1}$ and $\mu_{\delta}^{\text{M28}} = -8.92 \pm 0.03$ mas yr$^{-1}$ (Vasiliev & Baumgardt 2021) with Gaia Early Data Release 3 (EDR3). Hence, the relative proper motion of PSR J1824−2452A with respect to

M28 is $\Delta\mu_{\alpha} = 0.03 \pm 0.05$ mas yr$^{-1}$ and $\Delta\mu_{\delta} = 1.1 \pm 0.8$ mas yr$^{-1}$. Combining the M28 distance $D = 5.4 \pm 0.1$ kpc estimated by Baumgardt & Vasiliev (2021), we obtained the transverse space velocity $v_{\perp} = 28 \pm 20$ km s$^{-1}$ for PSR J1824−2452A, which is smaller than the typical escape velocity ($\approx 50$ km s$^{-1}$) of a globular cluster. Therefore, the pulsar is probably (as expected) bound to M28.

## 9 SUMMARY AND FUTURE PROSPECTS

In this MSPSRπ release paper, we have presented VLBI astrometric results for 18 MSPs, including a re-analysis of three previously published sources. From the 18 sources, we detect significant parallaxes for all but three. For each MSPSRπ pulsar, at least one self-calibratable in-beam calibrator was identified to serve as the reference source of relative astrometry. In three cases, 1D interpolation, a more complex observing and data reduction strategy, is adopted to further suppress propagation-related systematic errors. Among the three pulsars, PSR J1939+2134 is the brightest MSP in the northern hemisphere. Hence, we took one step further to perform inverse-referenced 1D interpolation using PSR J1939+2134 as the in-beam calibrator. Compared to the pioneering Multi-View study of Rioja et al. (2017) at 1.6 GHz, the larger number of observations and targets here provides more opportunities to characterize the interpolation performance, which is crucial for ultra-precise astrometric calibration schemes proposed for future VLBI arrays incorporating the Square Kilometre Array (SKA). Based on a small sample of three, we found that $\eta_{\text{EFAC}}$ has consistently inflated after applying 1D interpolation (see Section 4.1.3). This inflation implies that the higher-order terms of in the phase screen approximation may not be negligible, and could become the limiting factor of the ultra-precise SKA-based astrometry using the Multi-View strategy. Further investigations of the same nature, especially using temporally simultaneous (in-beam) calibrators, at low-observing frequencies are merited and encouraged.

In this paper, we present two sets of astrometric results – the quasi-VLBI-only results (see Section 4) and the VLBI + timing results (see Section 5). Both sets of astrometric results are inferred with the astrometry inference package `sterne`$^7$. The former set of results is largely independent of any input based on pulsar timing, making use only of orbital parameters as priors in the inference of orbital reflex motion, which affects only four pulsars from our sample and is near-negligible in any case. The latter, however, additionally adopts the latest available timing parallaxes and proper motions as priors of inference wherever possible, affecting all pulsars in our sample. While the latter approach typically gives more precise results, we note that this is dependent on the accuracy of the timing priors, and identify seven pulsars (PSR J0610−2100, PSR J1643−1224, PSR J1730−2304, PSR J1738+0333, PSR J1824−2452A, PSR J1853+1303, and PSR J1910+1256) for which disagreement between the quasi-VLBI-only and the timing priors mean that the VLBI + timing results should be treated with caution. From the VLBI perspective, we looked into possible causes of additional astrometric uncertainties, including non-optimal calibrator quality (see Sections 8.6 and 8.9) and calibrator structure evolution (see Section 5). In future, proper motion uncertainties (including any unaccounted systematic error due to calibrator structure evolution) can be greatly reduced with only $\lesssim 2$ extra observations per pulsar. For example, a 10-yr time baseline can improve the current VLBI-only proper motion precision by roughly a factor of 8.

From the VLBI + timing parallaxes $\sigma'$, we derived distances $D$ using equation (9). Incorporating the PDFs of $D$ and proper motions
$\{\mu'_\alpha, \mu'_s\}$, we estimated the transverse space velocities $v_\perp$ for 16 pulsars with significant distance determination, and found their $v_\perp$ to be generally on the smaller side of the previous estimates (Hobbs et al. 2005; Gonzalez et al. 2011). Boodram & Heinke (2022) propose that MSPs must have near-zero space velocities in order to explain the Fermi Galactic Centre excess. Our relatively small space velocities inferred for 16 MSPs suggest that MSPs may not be ruled out as the source of the Galactic $\gamma$-ray Centre excess. If the multimodal feature of the $v_\perp$ is confirmed with a sample of ~50 MSPs, it may serve as a kinematic evidence for alternative formation channels of MSPs (Bailyn & Grindlay 1990; Gautam et al. 2022, also see Ding et al. 2022 as an analogy).

In addition, we estimated the radial accelerations of pulsars with their distances and proper motions (see Section 7), which allows us to constrain the intrinsic spin period derivative $\dot{P}_s^{\text{int}}$ and the intrinsic orbital decay $\dot{P}_b^{\text{int}}$ (see Table 8). We used the improved $\dot{P}_s^{\text{int}}$ of PSR J1730−2304 to place a refined upper limit to the death line of $\gamma$-ray pulsars (see Section 8.7), and the $\dot{P}_b^{\text{int}}$ (of PSR J1012+5307 and PSR J1738+0333) to constrain alternative theories of gravity (see Section 7.3). As already noted by Ding et al. (2020b), the orbital decay tests (of gravitational theories) with the three viable MSPSRπ systems (i.e. PSR J1012+5307, PSR J1537+1155, and PSR J1738+0333) will be limited by the distance uncertainties, as parallax precision improves much slower than the $\dot{P}_b^{\text{obs}}$ precision (Bell & Bailes 1996).

Moreover, we detected significant angular broadening of PSR J1643−1224, which we used to (1) provide an independent check of the postulated connection between the HII region Sh 2-27 and the main scattering screen of PSR J1643−1224, and (2) test the scintillation models proposed by Mall et al. (2022). In future scintillation model inferences, angular broadening and temporal broadening measurements, where available, are suggested to be used as priors using equation (20), in order to achieve more reliable (and potentially more precise) scintillation model parameters. Such an inference would also complete the geometric information of the deflection routes (using equation (19), for example, in the two-screen case).

## DATA AND CODE AVAILABILITY

(i) The VLBA data can be downloaded from the NRAO Archive Interface at https://data.nrao.edu with the project codes in Table 1.

(ii) Image models of phase calibrators are provided at https://github.com/dingswin/calibrator_models_for_astrometry.

(iii) Supplementary materials supporting this paper can be found at https://github.com/dingswin/publication_related_materials.

(iv) The data reduction pipeline psrvlbireduce is available at https://github.com/dingswin/psrvlbireduce (version ID: b8ddafd).

(v) The astrometry inference package sterne can be accessed at https://github.com/dingswin/sterne.

## ACKNOWLEDGEMENTS

The authors thank Y. Kovalev, L. Petrov and N. Wex for useful discussions. HD is supported by the OzGrav scholarship through the Australian Research Council project number CE170100004. PF thanks the continued support by the Max-Planck-Gesellschaft. SC, JMC, EF, DK, and TJWL are members of the NANOGrav Physics Frontiers Center, which is supported by NSF award PHY-1430284. This work is 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. Pulsar research at Jodrell Bank Observatory is supported by a consolidated grant from STFC. Pulsar research at UBC is supported by an NSERC Discovery Grant and by the Canadian Institute for Advanced Research. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). We acknowledge financial support from the ‘Programme National de Cosmologie et Galaxies’ (PNCG) and ‘Programme National Hautes Energies’ (PNHE) of INSU,CNRS, France. The Parkes radio telescope (Murriyang) is part of the Australia Telescope, which is funded by the Commonwealth Government for operation as a National Facility managed by CSIRO. The Wisconsin Hα Mapper and its Hα Sky Survey have been funded primarily by the National Science Foundation. The facility was designed and built with the help of the University of Wisconsin Graduate School, Physical Sciences Lab, and Space Astronomy Lab. NOAO staff at Kitt Peak and Cerro Tololo provided on-site support for its remote operation. Data reduction and analysis was performed on OzSTAR, the Swinburne-based supercomputer. 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. sterne as well as other data analysis involved in this work made use of numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), astropy (The Astropy Collaboration et al. 2018) and the bilby package (Ashton et al. 2019).

## REFERENCES

Abazajian K. N., Kaplinghat M., 2012, Phys. Rev. D, 86, 083511

Abbott B. P. et al., 2016, Phys. Rev. Lett., 116, 061102

Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982, Nature, 300, 728

Antoniadis J., 2021, MNRAS, 501, 1116

Antoniadis J., Van Kerkwijk M., Koester D., Freire P., Wex N., Tauris T., Kramer M., Bassa C., 2012, MNRAS, 423, 3316

Antoniadis J. et al., 2013, Science, 340, 448

Antoniadis J. et al., 2022, MNRAS, 510, 4873

Armstrong J., Rickett B., Spangler S., 1995, ApJ, 443, 209

Arzoumanian Z. et al., 2018, ApJS, 235, 37

Arzoumanian Z. et al., 2020, ApJ, 905, L34

Ashton G. et al., 2019, ApJS, 241, 27

Bailer-Jones C. A., 2015, PASP, 127, 994

Bailer-Jones C., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147

Bailes M. et al., 2020, PASA, 37, e028

Bailyn C. D., Grindlay J. E., 1990, ApJ, 353, 159

Bartel N., Herring T. A., Ratner M. I., Shapiro I. I., Corey B. E., 1986, Nature, 319, 733

Bassa C. et al., 2016, MNRAS, 460, 2207

Baumgardt H., Vasiliev E., 2021, MNRAS, 505, 5957

Beasley A., Conway J., 1995, in ASP Conf. Ser. Vol. 82, Very Long Baseline Interferometry and the VLBANRAO. Astron. Soc. Pac., San Francisco, p. 327

Becker R. H., White R. L., Helfand D. J., 1995, ApJ, 450, 559

Bell J. F., Bailes M., 1996, ApJ, 456, L33

Binney J., Tremaine S., 2011, Galactic Dynamics. Princeton Univ. Press, Princeton

Boodram O., Heinke C. O., 2022, MNRAS, 512, 4239

Bovy J., 2020, preprint (arXiv:2012.02169)

Bower G. C. et al., 2014, ApJ, 780, L2

Brisken W. F., Benson J. M., Goss W., Thorsett S., 2002, ApJ, 571, 906

Burgay M. et al., 2003, Nature, 426, 531

Burgay M. et al., 2006, MNRAS, 368, 283

Cameron A. et al., 2018, MNRAS, 475, L57
Carr B. J., 1980, A&A, 89, 6

Charlot P. et al., 2020, A&A, 644, A159

Chatterjee S., Cordes J., Vlemmings W., Arzoumanian Z., Goss W., Lazio T., 2004, ApJ, 604, 339

Chatterjee S. et al., 2009, ApJ, 698, 250

Chen S. et al., 2021, MNRAS, 508, 4970

Condon J. J., Cotton W., Greisen E., Yin Q., Perley R. A., Taylor G., Broderick J., 1998, AJ, 115, 1693

Cordes J. M., Lazio T. J. W., 2002, preprint (arXiv:astro-ph/0207156)

Dehnen W., Binney J., 1998, MNRAS, 294, 429

Deller A., Verbiest J., Tingay S., Bailes M., 2008, ApJ, 685, L67

Deller A. T., Bailes M., Tingay S. J., 2009, Science, 323, 1327

Deller A. T. et al., 2011, PASP, 123, 275

Deller A., Boyles J., Lorimer D., Kaspi V., McLaughlin M., Ransom S., Stairs I., Stovall K., 2013, ApJ, 770, 145

Deller A. T. et al., 2016, ApJ, 828, 8

Deller A. T., Weisberg J. M., Nice D. J., Chatterjee S., 2018, ApJ, 862, 139

Deller A. T. et al., 2019, ApJ, 875, 100

Desvignes G. et al., 2016, MNRAS, 458, 3341

Detweiler S., 1979, ApJ, 234, 1100

Ding H., 2022, PhD thesis, Swinburne University of Technology

Ding H. et al., 2020a, MNRAS, 498, 3736

Ding H., Deller A. T., Freire P., Kaplan D. L., Lazio T. J. W., Shannon R., Stappers B., 2020b, ApJ, 896, 85

Ding H., Deller A. T., Fonseca E., Stairs I. H., Stappers B., Lyne A., 2021a, ApJ, 921, L19

Ding H., Deller A. T., Miller-Jones J. C. A., 2021b, PASA, 38, e048

Ding H., Deller A., Lower M., Shannon R., 2022, preprint (arXiv:2201.07376)

Eardley D. M., 1975, ApJ, 196, L59

Edwards R. T., Bailes M., 2001, ApJ, 553, 801

Edwards R. T., Hobbs G., Wolszczan A., Anderson S., 2006, MNRAS, 372, 1549

Faisal Alam M. et al., 2021, ApJS, 252, 4

Ferdman R. et al., 2020, Nature, 583, 211

Fermi-LAT Collaboration, 2022, Science, 376, 521

Fienga A., Avdellidou C., Hanuš J., 2020, MNRAS, 492, 589

Fomalont E. B., Kopeikin S. M., 2003, ApJ, 598, 704

Fonseca E., Stairs I. H., Thorsett S. E., 2014, ApJ, 787, 82

Foster R. S., Backer D. C., 1990, ApJ, 361, 300

Foster R., Cadwell B., Wolszczan A., Anderson S., 1995, ApJ, 454, 826

Freire P. C. C. et al., 2012, MNRAS, 423, 3328

Gaia Collaboration et al., 2018, A&A, 616, A1

Gaia Collaboration et al., 2022, A&A, preprint (arXiv:2208.00211)

Gautam A., Crocker R. M., Ferrario L., Ruiter A. J., Ploeg H., Gordon C., Macias O., 2022, Nat. Astron., 6, 703

Gold T., 1968, Nature, 218, 731

Goncharov B. et al., 2021, ApJ, 917, L19

Gonzalez M. et al., 2011, ApJ, 743, 102

Gravity Collaboration et al., 2018, A&A, 615, L15

Greisen E. W., 2003, in Heck A., ed., Astrophys. Space Sci. Libr. Vol. 285, Information Handling in Astronomy - Historical Vistas. Kluwer Academic Publishers, Dordrecht, p. 109

Guillemot L. et al., 2016, A&A, 587, A109

Guo Y. et al., 2021, A&A, 654, A16

Harris C. R. et al., 2020, Nature, 585, 357

Helfand D., Taylor J., Backus P., Cordes J., 1980, ApJ, 237, 206

Herbst W., 1975, AJ, 80, 503

Hewish A., Bell S., Pilkington J., Scott P., Collins R., 1969, Nature, 224, 472

Hobbs G., Lorimer D., Lyne A., Kramer M., 2005, MNRAS, 360, 974

Hobbs G., Lyne A., Kramer M., 2010, MNRAS, 402, 1027

Hulse R. A., Taylor J. H., 1975, ApJ, 195, L51

Hyland L. J., Reid M. J., Ellingsen S. P., Rioja M. J., Dodson R., Orosz G., Masson C. R., MacCallum J., 2022, ApJ, 932, 52

Igoshev A., Verbunt F., Cator E., 2016, A&A, 591, A123

Imai H., Sakai N., Nakanishi H., Sakanoue H., Honma M., Miyaji T., 2012, PASJ, 64, 142

Jacoby B., Bailes M., Ord S., Edwards R., Kulkarni S., 2009, ApJ, 699, 2009

Janssen G., Stappers B., Kramer M., Nice D., Jessner A., Cognard I., Purver M., 2008, A&A, 490, 753

Jennings R. J., Kaplan D. L., Chatterjee S., Cordes J. M., Deller A. T., 2018, ApJ, 864, 26

Kaplan D. L. et al., 2016, ApJ, 826, 86

Kerr M. et al., 2020, PASA, 37, e020

Kettenis M., van Langevelde H. J., Reynolds C., Cotton B., 2006, in Gabriel C., Arviset C., Ponz D., Enrique S., eds, ASP Conf. Ser. Vol. 351, Astronomical Data Analysis Software and Systems XV. Astron. Soc. Pac., San Francisco, p. 497

Kirsten F., Vlemmings W., Campbell R. M., Kramer M., Chatterjee S., 2015, A&A, 577, A111

Kopeikin S., 1996, ApJ, 467, L93

Kramer M. et al., 2021, Phys. Rev. X, 11, 041050

Lazaridis K. et al., 2009, MNRAS, 400, 805

Lazarus P. et al., 2016, ApJ, 831, 150

Lentati L., Kerr M., Dai S., Shannon R., Hobbs G., Oslowski S., 2017, MNRAS, 468, 1474

Lestrade J.-F., Rogers A., Whitney A., Niell A., Phillips R., Preston R., 1990, AJ, 99, 1663

Li Z. et al., 2018, MNRAS, 476, 399

Lindegren L. et al., 2021, A&A, 649, A4

Lobanov A. P., 1998, A&A, 330, 79

Lorimer D. et al., 1995, ApJ, 439, 933

Lorimer D. et al., 2006, MNRAS, 372, 777

Lundgren S., Cordes J., Foster R., Wolszczan A., Camilo F., 1996, ApJ, 458, L33

Lutz T. E., Kelker D. H., 1973, PASP, 85, 573

Lyne A., Brinklow A., Middleditch J., Kulkarni S., Backer D. A., Clifton T., 1987, Nature, 328, 399

Macquart J.-P. et al., 2020, Nature, 581, 391

Madison D. R., Chatterjee S., Cordes J. M., 2013, ApJ, 777, 104

Mall G. et al., 2022, MNRAS, 511, 1104

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ, 129, 1993

Mannings A. G. et al., 2021, ApJ, 917, 75

McMillan P. J., 2017, MNRAS, 465, 76

Mingarelli C. M. F., Anderson L., Bedell M., Spergel D. N., 2018, preprint (arXiv:e-print)

Nice D. J., Taylor J. H., 1995, ApJ, 441, 429

Nice D. J., Splaver E. M., Stairs I. H., 2001, ApJ, 549, 516

Ocker S. K., Cordes J. M., Chatterjee S., 2020, ApJ, 897, 124

Pacini F., 1968, Nature, 219, 145

Pallanca C., Mignani R., Dalessandro E., Ferraro F., Lanzoni B., Possenti A., Burgay M., Sabbi E., 2012, ApJ, 755, 180

Park R. S., Folkner W. M., Williams J. G., Boggs D. H., 2021, AJ, 161, 105

Pearson K., 1895, Notes on Regression and Inheritance in the Case of Two Parents Proceedings of the Royal Society of London, Vol. 58., p. 214

Perera B. et al., 2019, MNRAS, 490, 4666

Perger K. et al., 2018, MNRAS, 477, 1065

Peters P. C., 1964, Phys. Rev., 136, 1224

Phillips D. F., Ravi A., Ebadi R., Walsworth R. L., 2021, Phys. Rev. Lett., 126, 141103

Piffl T. et al., 2014, MNRAS, 445, 3133

Pradel N., Charlot P., Lestrade J.-F., 2006, A&A, 452, 1099

Reardon D., Coles W., Hobbs G., Ord S., Kerr M., Bailes M., Bhat N., Venkatraman Krishnan V., 2019, MNRAS, 485, 4389

Reardon D. J. et al., 2020, ApJ, 904, 104

Reardon D. et al., 2021, MNRAS, 507, 2137

Reid M., Menten K., Brunthaler A., Zheng X., Moscadelli L., Xu Y., 2009, ApJ, 693, 397

Reid M. et al., 2019, ApJ, 885, 131

Rioja M. J., Dodson R., Orosz G., Imai H., Frey S., 2017, AJ, 153, 105

Roebber E., 2019, ApJ, 876, 55

Sayer R., Nice D., Taylor J., 1997, ApJ, 474, 426

Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192

Shepherd M. C., Pearson T. J., Taylor G. B., 1994, Bull. Astron. Soc., 26, 987

Shklovskii I. S., 1970, Soviet Ast., 13, 562

Siemens X., Ellis J., Jenet F., Romano J. D., 2013, Class. Quantum Gravity, 30, 224015
Sokolovsky K. V., Kovalev Y. Y., Pushkarev A. B., Lobanov A. P., 2011, A&A, 532, A38

Stovall K. et al., 2018, ApJ, 854, L22

Tauris T. et al., 2017, ApJ, 846, 170

Taylor J. H., Weisberg J. M., 1982, ApJ, 253, 908

The Astropy Collaboration et al., 2018, AJ, 156, 123

Tiburzi C. et al., 2016, MNRAS, 455, 4339

Vallisneri M. et al., 2020, ApJ, 893, 112

van der Wateren E. et al., 2022, A&A, 661, A57

Vasiliev E., Baumgardt H., 2021, MNRAS, 505, 5978

Verbiest J., Weisberg J., Chael A., Lee K., Lorimer D., 2012, ApJ, 755, 39

Vigeland S. J., Deller A. T., Kaplan D. L., Istrate A. G., Stappers B. W., Tauris T. M., 2018, ApJ, 855, 122

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wang J. B. et al., 2017, MNRAS, 469, 425

Weisberg J. M., Huang Y., 2016, ApJ, 829, 55

Wolszczan A., 1991, Nature, 350, 688

Yang J., Paragi Z., van der Horst A., Gurvits L., Campbell R., Giannios D., An T., Komossa S., 2016, MNRAS, 462, L66

Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29

Zhang Y., An T., Frey S., 2020, Science Bulletin, 65, 525

Zhao J., Freire P. C., Kramer M., Shao L., Wex N., 2022, Class. Quantum Gravity, 39, 11LT01

Zhu W. W. et al., 2015, ApJ, 809, 41

Zhu W. et al., 2019, MNRAS, 482, 3249

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a TeX/LATEX file prepared by the author.