MeerKAT Pulsar Timing Array parallaxes and proper motions

We have determined positions, proper motions, and parallaxes of $77$ millisecond pulsars (MSPs) from $\sim3$ years of MeerKAT radio telescope observations. Our timing and noise analyses enable us to measure $35$ significant parallaxes ($12$ of them for the first time) and $69$ significant proper motions. Eight pulsars near the ecliptic have an accurate proper motion in ecliptic longitude only. PSR~J0955$-$6150 has a good upper limit on its very small proper motion ($<$0.4 mas yr$^{-1}$). We used pulsars with accurate parallaxes to study the MSP velocities. This yields $39$ MSP transverse velocities, and combined with MSPs in the literature (excluding those in Globular Clusters) we analyse $66$ MSPs in total. We find that MSPs have, on average, much lower velocities than normal pulsars, with a mean transverse velocity of only $78(8)$ km s$^{-1}$ (MSPs) compared with $246(21)$ km s$^{-1}$ (normal pulsars). We found no statistical differences between the velocity distributions of isolated and binary millisecond pulsars. From Galactocentric cylindrical velocities of the MSPs, we derive 3-D velocity dispersions of $\sigma_{\rho}$, $\sigma_{\phi}$, $\sigma_{z}$ = $63(11)$, $48(8)$, $19(3)$ km s$^{-1}$. We measure a mean asymmetric drift with amplitude $38(11)$ km s$^{-1}$, consistent with expectation for MSPs, given their velocity dispersions and ages. The MSP velocity distribution is consistent with binary evolution models that predict very few MSPs with velocities $>300$ km s$^{-1}$ and a mild anticorrelation of transverse velocity with orbital period.


INTRODUCTION
Studying the transverse velocities of pulsars enables us to better understand their origins in the Milky Way, the evolutionary scenarios leading up to their birth, and the kick velocities that they receive when they are born (Lyne et al. 1982;Cordes & Chernoff 1997).Distances and velocities of pulsars also enable us to constrain the dynamical models of supernovae (e.g.Gaensler & Johnston 1995) and understand the kinematics of binary evolution (e.g.Gonzalez et al. 2011;Freire et al. 2011).
Transverse velocities for radio pulsars are primarily determined using two different methods: firstly, measurements of pulsar proper motions and parallaxes through pulsar timing or Very Long Baseline Interferometry (VLBI) (e.g.Reardon et al. 2021;Ding et al. 2023), and secondly, using interstellar scintillation patterns in a pulsar's dynamic spectrum (e.g.Bogdanov et al. 2002;Reardon et al. 2020).Slow (or young) pulsars possess timing noise that makes simple fits for proper motion unreliable and most of the first proper motions ★ E-mail:msh.ph.ir@gmail.comdetermined were reliant on interferometric techniques (e.g., Lyne et al. 1982;Bailes et al. 1990;Harrison et al. 1993).
A breakthrough study of pulsar transverse velocities was made by Hobbs et al. (2005).They collated all the pulsar proper motions from the literature and significantly added to them with a new timing technique from the Jodrell Bank timing data to explore a total sample of 233 pulsar proper motions.They found the two-dimensional velocities of young pulsars range from a few tens to ≈ 1600 km s −1 .They used a novel deconvolution technique to derive a mean 3D pulsar birth velocity of approximately 400 ± 40 km s −1 for young (<3 Myr) pulsars, suggesting that pulsars receive large kicks at birth.Such high kick velocities imply that the local convective instability in the collapsed stellar core (Lai et al. 2001) is unlikely to be a pulsar kick mechanism (the mechanism that causes the neutron star to get kicked with a different velocity compared to its progenitor star after the supernova explosion), and they favoured more energetic mechanisms such as global asymmetric perturbations or neutrino emission in the supernovae.They found that a single Maxwellian distribution adequately fit the distribution of velocities, and derived a one-dimensional root-mean-square (rms) velocity of 265 km s −1 , much higher than their progenitors, the OB stars which typically only possess velocities of a few tens of km s −1 .
Slowly rotating pulsars cannot be timed sufficiently accurately to yield reliable timing parallaxes and so for many years pulsar parallaxes were rare.In the 2000s, dedicated very long baseline arrays such as the Very Long Baseline Array and efficient pipelines made pulsar proper motion and parallax studies more routine, and several large-scale surveys were conducted on both slow-spinning pulsars and MSPs that greatly expanded our knowledge of their kinematics and distances (Chatterjee et al. 2009;Deller et al. 2019;Ding et al. 2023).Millisecond pulsar timing arrays have also started significantly contributing to the number of MSPs with measured timing parallaxes.PTAs now routinely obtain sub-microsecond timing precision, and this is permitting the measurement of many MSP parallaxes via the pulse timing method (e.g., Desvignes et al. 2016;Arzoumanian et al. 2018;Reardon et al. 2021).In some special cases, the non-zero apparent orbital period derivative of pulsars enables a very model independent estimate of their distance (e.g., Bell & Bailes 1996).
In this work, we concentrate on measurements of transverse velocities for MSPs, and substantially increase the sample with which to study their distribution.For the purposes of this paper, MSPs are defined as having a spin period of  ≤ 30 ms and a spin-down rate of  ≤ 10 −17 .
MSPs are well known to have extremely precise spin periods and are accurate clocks.Their short spin periods and timing stability enhances the reliability of timing models.We are often able to predict and measure the arrival times of pulses (ToAs) of MSPs with submicrosecond precision.By having such precise ToAs, we can use the pulsar timing technique to measure the positions of MSPs often to sub-milliarcsecond accuracies.Having several years of data allows us to also measure the proper motions with the precision of better than milliarcseconds per year and even reach sub-milliarcsecond parallaxes.Using these proper motions and parallaxes, we can then derive the distances and transverse velocities with a precision of sometimes the order of a few percent and few km s −1 , respectively (e.g., Matthews et al. 2016;Reardon et al. 2021).To do so, we need to observe pulsars regularly over the course of at least one year in order to correct for the Rømer delay, which causes a nearly sinusoidal variation of the pulse's travel time due to the Earth's orbit around the Sun and allows us to measure the position of the pulsar.A transverse motion of the pulsar on the sky causes a linearly increasing sinusoid with a period of one year to appear in the timing residuals, and usually after about 3 years the proper motions can be determined to high accuracy.The curvature of the pulse train's wavefront due to its finite distance of origin allows us to measure a pulsar's parallax, provided it is not too far from the ecliptic plane -unfortunately an MSP at the ecliptic pole has almost no discernable parallax timing signature.One pulsar in our sample, (PSR J0711-6830) is almost 83 • from the ecliptic plane and has never had its timing parallax determined.
The applications enabled by possessing pulsar parallaxes are many.For example, combining a parallax-derived distance and the pulsar's dispersion measure (as will be explained in Section 3.1) can be used for developing Galactic electron density distribution models such as TC93 (Taylor & Cordes 1993), NE2001 (Cordes & Lazio 2002, 2003), and YMW16 (Yao et al. 2017).These models allow us to make estimates of pulsar distances when parallaxes are unavailable, using the DM as a proxy for distance.When combined with a proper motion, these yield an estimate of the pulsar's transverse velocity.In addition, the Galactic magnetic field parallel to the line of sight can be mapped out through the measurements of the distances, DMs and the Faraday rotation of pulsars (e.g., Han et al. 2006;Spiewak et al. 2022).Moreover, pulsar distances can be used to correct for the Shklovskii effect, which otherwise is the main contaminant to the observed period and orbital period derivatives from their intrinsic values (Shklovskii 1970).This leads to improved estimates of the pulsar's characteristic age, magnetic field strength and in tests of General Relativity that require the orbital period derivative (Camilo et al. 1994;Bell & Bailes 1996).
The focus of this paper is on MSP parallaxes, proper motions and hence distances and transverse velocities.We use these to assemble the largest-ever sample of MSPs with transverse velocities to study how they might differ from the slowly rotating pulsars, and also as a population.Determining pulsar transverse velocities requires measurement of pulsar distances and parallaxes, and consequently, we have selected a sample of MSPs observed by the MeerKAT timing program (MeerTime, Bailes et al. 2018) aimed at improving as many proper motions and parallaxes as possible.
In Section 2, we explain how the MeerKAT observations were performed and how the ToAs were determined.In Section 3, we describe our methods for timing each pulsar and the modelling of timing noise, and the methods employed for the measurement of parallaxes and proper motions.In Section 4, we present the positions, proper motions, and parallaxes.In Section 5, we derive distances and velocities for our sample.We also present the velocity distributions and the dispersions of velocity components in Galactocentric coordinates.In Section 6, we compare our results to the previous work and provide discussions about the velocity distributions and the velocity dispersions.In Section 7, we list our conclusions and discuss the implications of our findings.

THE DATA SETS
We used the first MeerKAT pulsar timing array (MPTA) data release provided by Miles et al. (2023) as the foundation of our study.This data set contains 78 MSPs observed with approximately a two week cadence with the 64-dish MeerKAT radio telescope over a period of ∼ 2.5 years, starting from early 2019.We augmented the timing baseline of the data set to ∼ 3 years by adding all data taken after the first data release and up to May 2022.
The observations were made with the L-band receiver, operating at the centre frequency of 1284 MHz with 856 MHz of frequency bandwidth (Bailes et al. 2020), and recorded with 1024 frequency channels.Following Spiewak et al. (2022), we removed the top and bottom 48 channels (∼ 10% of the total bandwidth) from the L-band receiver where response roll-off affects the signal-to-noise ratio, and the remaining 928 channels were averaged in time and polarization to form Stokes I profiles.The typical minimum length of the observations was 4 minutes, but depending on the brightness of the pulsars and the orbital phase of the binary pulsars, they could be up to 4 hours long.
Timing baselines of observations ranged from 2.33 yr (for J1652−4838) to 3.26 yr (for PSR J1909−3744).We used only 1.98 years of PSR J1713+0747 observations taken prior to a change in its pulse shape between April 16 and17, 2021 (e.g., Jennings et al. 2022) that spoil its timing properties.Channels affected by radio frequency interference (RFI) were zero-weighted using MeerGuard, a modified version of CoastGuard (Lazarus et al. 2016) for removing RFI from the MeerKAT observations.Almost no observations were deleted due to RFI, although a small percentage had to be removed because of observatory set-up issues such as poor phasing of the array.These outliers were very obvious and did not present any analysis issues.
We formed sub-banded pulse profiles by averaging frequency channels by a factor of 58 in order to detect any frequency-dependent variation of dispersion measure (DM) across the band (Keith et al. 2013;Jones et al. 2017;Donner et al. 2020;Nobleson et al. 2022) in addition to any frequency-dependent trends that could be attributed to the intrinsic profile shape changes across frequency (Kramer et al. 1999;Ahuja et al. 2007;Xu et al. 2021).This left us with 16 subbands that could be used to model dispersion measure variations.
We produced smoothed, high signal-to-noise, standard 'template' profiles for each pulsar from the sub-banded pulse profiles for use in measuring ToAs.Each template represents a set of denoised analytical pulse profiles that evolve across frequency (Pennucci et al. 2014;Pennucci 2019), using the PulsePortraiture software package1 .Employing a 2D portrait enabled us to calculate ToAs for each subband, while simultaneously examining profile stability and correct for the frequency-dependent trends (if they are significant) due to time-dependent dispersion measure.
The most recently available ephemerides for the pulsars (Spiewak et al. 2022, private communication) were used to initiate the timing analysis.Using the pat software tool in the psrchive package (van Straten et al. 2012), we measured the ToAs for all the observations, setting the threshold for timing at a signal-to-noise ratio SNR > 10.The number of ToAs ranged from 269 for PSR J1327−0757 to 3074 for PSR J1909−3744.We describe the timing analysis of this set of ToAs in the next section.

Timing and Noise analysis
The astrometric, period, spin-down, DM, and (where necessary) binary parameters of every pulsar are described in a timing model.Using the pulsar timing software tempo2 (Hobbs et al. 2006;Edwards et al. 2006), we refined the timing models by fitting for the par-file parameters.tempo2 uses a linear least-squares approach to fitting parameter values, such that some of the parameters that have non-linear forms need to be fit multiple times in order to yield the best timing residuals.We also include systematic time jumps in our timing models derived from the entire population.We used the JPL DE440 ephemeris as the model of the Earth's orbit around the Solar System barycentre (Park et al. 2021), and the ToAs were referred to the TT(BIPM2020)2 , and the Barycentric Coordinate Time (TCB) was used as the coordinate time standard for calculating the orbits in the Solar System (Seidelmann & Fukushima 1992).
Once accurate timing models had been produced for each pulsar, they were then employed for the noise analysis using the Bayesian inference software temponest (Lentati et al. 2014).The advantage of using temponest in our analysis is that it uses Multinest (Feroz & Hobson 2008;Feroz et al. 2009) to sample both timing parameter and noise parameter spaces efficiently.The noise parameters were modeled by three white noise terms, an achromatic red noise term, and a dispersion noise (i.e.due to DM variations) term.
• White noise: there are three contributions for the white noise as follows: (i) EFAC: while calculating the ToAs described in Section 2, the radiometer noise is a certain contributor to the ToA uncertainties  i .EFAC is modelled using a parameter, which is multiplied by the ToA uncertainties and scales them.An EFAC estimate is made for each pulsar/receiver combination.If EFAC deviates too far from unity the error estimates are not well understood.(ii) EQUAD: The EFAC parameter might not be able to fully account for the actual error, and we can broaden the errors by defining EQUAD parameter which is added in quadrature to EFAC.As for the EFAC parameter, an EQUAD estimate is made for each pulsar/receiver combination.Again if EQUAD is much greater than the rms residual, this is concerning as it means there are systematic errors in the arrival times that are of unknown origin.(iii) ECORR: for sub-banded data, there is likely to be a correlation between ToAs at different frequencies that are collected simultaneously because of issues such as pulse jitter (NANOGrav Collaboration et al. 2015), and this introduces another noise term that can be described using ECORR.In this analysis, we used the TECORR estimator, as described in Bailes et al. (2020, Eq. 1).ECORR is also estimated for each pulsar/receiver combination.
• Red noise: this is thought to be related to the irregularities of the rotation of neutron stars (Shannon & Cordes 2010).Red noise is a stochastic, achromatic term with a power-law spectrum of the form where  red (  ) is the power spectral density of red noise with fluctuation frequency  ,  red is the amplitude in s yr 1/2 , and  red is the spectral index.
• Dispersion noise: The ToAs are affected by frequency dispersion caused by the ionized material in the interstellar medium (ISM), and ToAs at lower observing frequencies arrive with a delay with respect to the higher ones.This delay, measured in seconds, is equal to 1/(2.41 × 10 −4 ) DM/ 2 , where DM is dispersion measure (in pc cm −3 ), and  is the observing frequency (in MHz) (Lorimer & Kramer 2005).Due to the space motions of the Earth and pulsars, as well as the bulk motion of the ISM, DM shows a temporal variation (Keith et al. 2013;Jones et al. 2017).This variation introduces a dispersion noise term that can be modelled using a power law with the same form as the red noise.The only difference is the frequencydependent amplitude.
For every pulsar, we performed four noise analyses with the noise models including (1) DM, red, and white noise parameters (2) DM and white noise parameters only (3) red and white noise parameters only and (4) white noise parameters only.We calculated a log Bayesian evidence for every noise model.By subtracting the log Bayesian evidence values of each of two noise models we can measure the log Bayes factor of ln B. If B > 10 (corresponding to the log Bayes factor of 2.3) (Kass & Raftery 1995) the noise model with higher Bayesian evidence is more significant than the other and is preferred as the final noise model.The prior ranges for the noise parameters are listed in Table 1.

Position, Proper Motion, and Parallax Measurements
After determining a preferred noise model, we measured the pulsar's astrometric parameters using the ToAs, and ascertained which are measured with sufficient statistical significance for use in the later tangential velocity analysis.For each pulsar the astrometric parameters are the equatorial coordinates for the position (, ), the proper motions (  ≡  cos ,   ) and the parallax.Due to the high [−20,−10]  DM [0,6] sensitivity of MeerKAT observations, the position of pulsars can be obtained extremely well, and we did not need to sample for the position parameter; so, these remained fixed.We performed three analyses using temponest sampling for (1) proper motion, (2) parallax, and (3) proper motion and parallax.Bayes factors produced by temponest showed that, for all pulsars, the proper motion and parallax were detectable (with B > 10) and superior to models without them.In addition, we used ecliptic coordinates for the position (, ) and proper motion (  ≡  cos ,   ), and repeated the analyses and compared the Bayesian evidence for these solutions with the equatorial ones, and found them to be similar.Note that covariance between position and proper motion parameters was minimized by setting the epoch of the positions to be near the middle of the data set's baseline in time.Matthews et al. (2016) demonstrated that positional uncertainties reported in ecliptic coordinates are less than the equatorial coordinates.Therefore, we used timing analyses in ecliptic coordinates to derive the Galactic position (, ) and Galactic proper motion (  ≡  cos ,   ) of all pulsars using the Astropy package (Astropy Collaboration et al. 2013).

Distances from parallax measurement
The most direct method for measuring distance is parallax  =  −1 , where  is trigonometric parallax in mas and  is in kpc.This method does not depend on any intrinsic physical property of the object and is purely geometric, and thus model-independent.For pulsars, the two most common methods of measuring parallax are using pulsar timing solutions and VLBI observations.As will be seen, these two methods dominate the distances used in section 5 for our transverse velocity analysis.

Velocities
The full 3D velocity of a celestial object relative to the Sun can be broken into two components: a radial velocity component and a transverse (or tangential) velocity component.The radial velocity for stars in the Milky Way is typically estimated using the Doppler effect, although this is not possible for pulsars, as no spectral lines are available unless they are in a binary system with a star.The transverse velocity (in km s −1 ), is the velocity perpendicular to the line of sight and is given by where  is the total proper motion, and  is the distance (Lorimer et al. 1997).

Correcting for the Local Standard of Rest and the Galactic differential rotation
Proper motions are usually measured relative to the Solar System's barycentre, but if we want to understand their origins as Milky Way objects, we need to adjust for the reference frame.To obtain the velocity of pulsars with respect to the Galactic Center (GC), we need to correct the relative proper motions for two effects: (1) the peculiar velocity of the Sun with respect to the local standard of rest (LSR) and ( 2) differential rotation of the Galaxy (Hobbs et al. 2005).There are two types of LSR, both of which are defined as reference frames at the location of the Sun.The first is the 'kinematic LSR' which orbits the Galaxy with the mean velocity of the neighbouring stars ( v  ).The second is the 'dynamic LSR' which rotates around the Galaxy in a perfectly circular orbit with a velocity that is exerted by the gravitational potential of the Galaxy ( c ).The difference between the two ( c − v  ) is defined as the "asymmetric drift" (Golubov et al. 2013;Ferreras 2019).
The Galactic disk rotates differentially, and the Galactic orbital periods of stars are dependent on their distances from the GC.The circular rotation speed of stars at distance  from the GC defines the rotation curve of the Galaxy.In the last decade, there have been a number of studies to improve the measured rotation curve of the Galaxy using a variety of different methods (e.g., Bovy et al. 2012;Bhattacharjee et al. 2014;Reid et al. 2014Reid et al. , 2019)).In this work, we used the results and rotation curve model from Reid et al. (2019).
To correct for differential rotation of the Galaxy, we calculated the distances of the pulsars from the Galactic centre, and used the universal rotation curve of the Galaxy from Reid et al. (2019, Tab. 4) to calculate the  and  components of the Galactic rotation at the position of MSPs.Finally, we subtracted these from the velocities of MSPs and derived their transverse velocities with respect to their own LSR.All velocities reported in this work are thus corrected for the peculiar velocity of the Sun and the rotation curve of the Galaxy.

RESULTS: ASTROMETRY
In this section we report our astrometric results, including positions, proper motions, parallaxes for our sample pulsars, and compare with MSPs for which these have been previously measured in the literature.

Positions and proper motions
In Table 2, positions from timing analyses in equatorial and ecliptic coordinates for 77 MSPs are given.These were not sampled using temponest and only fitted using tempo2.In addition, the positions in Galactic coordinates are derived from ecliptic coordinates using Astropy.All uncertainties shown are computed via the median with 68.3% equally-tailed credible intervals.
In Table 3, we show proper motions from timing analyses in both equatorial and ecliptic coordinates.In this table, the proper motions in Galactic coordinates are derived from the ecliptic proper motions using Astropy.The reported values and uncertainties are the lower bound (16%), median (50%), and upper bound (84%) of marginalized posterior probability distributions.The timing baseline of observations for each pulsar, as well as the best previous measurements of proper motions in equatorial coordinates and their corresponding references, are listed.In our study, our limited timing baseline means that our proper motion results are not superior to all previous results, mainly because the relative error is proportional to the observing span  −3/2 , and some of the pulsars have been timed for over a decade with other instruments.
Out of 77 pulsars, 69 had proper motions with statistical significance greater than 3, and are shown in the upper portion of Table 3.The remaining 8 pulsars (i.e.those listed in the lower portion of the Table ) consist of an isolated MSP (J1721−2457) with relatively poor rms timing residuals of 3.7s and less than two deg from the ecliptic, and seven binary MSPs (J0955−6150, J1022+1001, J1327−0755, J1653−2054, J1705−1903, J1802−2124, J1811−2405) that have either poor timing residuals or small proper motions (as low as 0.21(15) mas yr −1 for J0955−6150).PSR J1022+1001 is near the ecliptic plane, and as expected had a very accurate proper motion in ecliptic longitude, but a poor constraint in ecliptic latitude.For this pulsar, we set the proper motion and orbital parameters from the high accuracy VLBI study of Deller et al. (2016) and derived a parallax of 1.2(4) mas.
We typically require at least two years of data in order to separate the annual parallax from the proper motion of a pulsar, due to the motion of the pulsar against the background of the sky and the orbit of the Earth around the Sun.We note that it is not always simple to break this degeneracy for the pulsars with high rms timing residuals, and worse for pulsars in (wide) binary systems.A solution to this is to implement VLBI techniques, which can often measure pulsar parallaxes significantly better, e.g.Deller et al. (2008) measured the parallax for PSR J0437−4715 to be 6.396(54) mas with the Australian Long Baseline Array.
All but two of the proper motions listed in Table 3 are consistent with previous results within 3.One component of the proper motion of PSRs J1446−4701 is just over 3 inconsistent with the values in Reardon et al. (2021), but the total proper motion is the same to within 8%.The origin of this discrepancy is unclear.We checked the use of the Solar System ephemeris that Reardon et al. (2021) used for their timing analysis, and our parallax remained the same to within errors.This inconsistency might be due to the different noise models that were adopted or be due to different noise properties at different epochs.
Figure 1 presents the trajectories of the 69 MSPs since 5 Myr ago until the present in Galactic longitude () and Galactic latitude (), assuming zero radial velocity for the MSPs.Both effects of the peculiar velocity of the Sun and the Galactic differential rotation are subtracted from the observed proper motions of MSPs.In addition, the Galactic gravitational potential is assumed to be zero, because the oscillation period in the galactic potential (> 100Myr) is much larger than 5 Myr, and so the MSPs are simply moving along their corrected proper motion vectors.The area in gray indicates the part of the sky that the MeerKAT radio telescope is not able to reach.From the diagram it is obvious that the MSPs are both moving towards and away from the galactic plane, consistent with a relaxed ancient population.

Parallaxes
The parallax measurements made in this study, as well as the best previous measurements of those pulsars in the literature and their corresponding measurement techniques (timing or VLBI), are listed in Table 4 and shown in Figure 2 (> 3 detections).The parallax values and their uncertainties are derived from the marginalized posterior probability distributions from temponest.There are 35 of our MSPs that had parallax measurements with > 3 significance.12 of these 35 are new (PSRs J0125−2327, J0614−3329, J0636−3044, J1421−4409, J1629−6902, J1652−4838, J1658−5324, J1757−5322, J1801−1417, J1811−2405, J1933−6211 and J1946−5403).We improved the precision of 4 other parallaxes (of the 35 MSPs) by a factor of 1.25 to 6 including PSRs J1446−4701, J1455−3330, J2124−3358 and J2322−2650.For the MSPs with less than 3 significance in parallax, we calculated the 68% confidence upper level and listed them in the lower portion of Table 4.
To gain confidence in our parallax measurements, we compared them with their values derived from other authors in Figure 2 as 22 of our significant (> 3) parallaxes have previous measurements.For PSR J1713+0747, we measured the parallax of 0.88(10) mas using only 1.98 yr of data as its subsequent timing was compromised by a sudden step-change in its pulse profile and gradual decay back to the original profile.In addition to the previously measured parallax value of 0.763(21) mas by Reardon et al. (2021)  Of the 35 MSPs with significant measurements, only our parallaxes for PSRs J1643−1224 and J2322−2650 (shown in Figure 2) are markedly inconsistent with the best previous measurements.We discuss each of these two pulsars in turn below.Overall the standard deviation of the differences in relative parallaxes for the 23 pulsars with previous measurements in literature is only 10 percent, and if we exclude the two outliers discussed above, only 7 percent.

PSR J1643−1224
For PSR J1643−1224, we obtained the parallax value of  = 2.6 +0.5 −0.6 mas which is quite inconsistent with the VLBI value of  = 1.1(1) mas obtained by Ding et al. (2023) and also has the lower formal error of the measurements.
Matthews et al. ( 2016) measured a value of 0.7(6) mas using the NANOGrav nine-year data set and a detailed model for DM variation that is known to be significant.They compared their value with the measurement of Verbiest et al. (2009), who obtained a value of 2.2(4) mas which did not include a complex model for DM variation.Matthews et al. (2016) hence argued that the limited DM variation modelling might be the possible source of the inconsistency between their measurements.
Two more measurements come from Reardon et al. (2021), who obtain the parallax to be 0.82(17) mas using the 24 yr observations of the Parkes Pulsar Timing Array data set and Desvignes et al. (2016), who obtain a parallax of 1.17( 26   dence for chromatic variations (Miles et al. in prep), probably due to time-dependent scattering in the HII region which are not modelled by our techniques.Ocker et al. (2020) reported that there is a significant excess DM for PSR J1643−1224 consistent with its distance beyond the HII region.
We examined whether our parallax might be radio-frequency dependent, by measuring the parallax using just the top half of our band (1284-1712 MHz), and we obtained a smaller value of 1.7(4) mas, which is in better agreement (1.5) with the measurement of Ding et al. (2023).In our noise analysis for this pulsar, using the entire band, the noise model with the white, red, and DM noise parameters had the highest Bayesian evidence; however, using only the top part of the band the noise model with the white and red noise parameters had the highest Bayesian evidence.We know that this pulsar is significantly affected by interstellar scattering, but this is not currently included in our noise model.Weaknesses in our noise model such as this may lead to small biases in the parallax measurement.This would be good topic for a new study.For the rest of this paper, we adopt the VLBI parallax.

PSR J2322−2650
For PSR J2322−2650, we have measured a parallax of 1.3(2) mas, which is very different from the previously measured value of 4.4(12) mas by Spiewak et al. (2018) but since it has 6 times lower uncertainty we favour our value.Our parallax distance for PSR J2322−2650 is thus 0.80 +0.18 −0.12 kpc which is much nearer the YMW16 distance of 0.76 kpc than the 0.23 kpc suggested by Spiewak et al. (2018).We speculate that the inconsistency between our values is due to the absence of any DM noise model in the earlier analysis of Spiewak et al. (2018) and their poorer timing residuals.

KINEMATICS OF THE MSP SAMPLE
In this section we study the kinematics of our MSP sample, using only those with parallax-based distances and velocities.While models of the ISM can be used to estimate distances to pulsars, it is well known that this can lead to large uncertainties.We choose here to base our kinematic analysis on MSPs with high quality distances (> 3 significance) only.
For a detailed analysis of the uncertainties involved in the use of DM-based distances via ISM models, we refer the reader to Appendix A.

Distances
We used the posterior samples of the parallax measurements to derive the posterior distribution of the distances.The distances of the 35 MSPs with > 3 significance are listed in the upper portion of Table 4.For the MSPs with < 3 significance in parallax, we calculated 68% confidence lower limit on distances and listed them in the lower portion of Table 4. Lutz & Kelker (1973) pointed out a bias in the parallax of stars and showed that the bias depends on the significance of the parallax with low-significance parallaxes more likely to be over-estimated due to volume effects.This has the effect of pushing low-significance pulsar parallax objects further from the Sun.Verbiest et al. (2010Verbiest et al. ( , 2012) ) confirmed the existence of the Lutz-Kelker (LK) bias in observations,  and they provided a method in order to correct the values of parallaxes and distances.For our discussions concerning pulsar distances and velocities, we will not adopt the LK corrections for individual objects to ease comparison with other recent authors (e.g., Matthews et al. 2016;Arzoumanian et al. 2018), but will discuss what effects it would have on the median velocities of the population.One problem with the LK correction is that the magnitude of the correction depends upon the assumed model for the spatial (i.e.Galactic) distribution of the pulsars, and also the luminosity function of MSPs.For highprecision parallaxes, the effects are negligible.

Tangential Velocities
From the proper motions and distances of the MSPs, we calculated the barycentric transverse velocities,  ⊥ , and then using the method explained in Section 3.4.1 we corrected them for the peculiar velocity of the Sun and the rotation curve of the Galaxy, assuming zero radial velocity with respect to our line of sight.In Table 5, a summary of 39 velocities for two subsets of our MSPs is given.The first subset is those MSPs with > 3 significance for their parallax and proper motion measurements.The second subset is those MSPs with > 3 significance for their proper motions measurements but < 3 significance for their parallax measurements.For MSPs in the second subset namely: PSRs J0610−2100, J1045−4509, J1545−4550, J1643−1224, J1843−1113, and J2234+0944, we used parallaxes measured by other authors listed in the fifth column of Table 4 to calculate their transverse velocities (these velocities are denoted by a † after the pulsar name).We used the positions, proper motions, and parallaxes of other MSPs in the pulsar catalogue (Manchester et al. 2005) and selected those with significant parallaxes (> 3) and calculated their transverse velocities using our method described in Section 3 in order to increase our MSP velocity sample.

Isolated versus Binary Millisecond Pulsars
We divided the MSPs into two subsets depending upon whether they are isolated or binary MSPs.The velocity distributions for the isolated and binary MSPs are presented in Figure 3.We obtained a mean velocity of 67(12) km s −1 for 16 isolated MSPs, and a mean veloc-ity of 72(8) km s −1 for the 49 binary MSPs.Between these MSPs, PSR J1024−0719 is a special case with an extremely wide binary with the orbital period of 2-20 kyr (Kaplan et al. 2016) (presented with a black bin in Figure 3).Also, PSR J1300+1240 is a planetary system (Wolszczan 1990) (presented with an olive bin in Figure 3).We did not include them for calculating the mean velocity of the binary MSPs as their evolutionary history is probably quite different from the other systems.The two-sample Kolmogorov-Smirnov test was performed for comparing the distributions of velocity samples.The maximum absolute difference between the cumulative distribution functions of the two samples was 0.17 and the -value was 0.87.This means with our limited sample that there is no statistical evidence that the two samples are drawn from different distributions.Whatever causes isolated MSPs to be created does not have an appreciably different effect on their velocities from that of the binaries.

Dispersion of the Velocity Components
We obtained cylindrical Galactocentric velocity components of   (radial),   (rotational), and  z (perpendicular) and then calculated the components that are at least 70 • from the line of sight (i.e. for the components that are nearly perpendicular to the line of sight direction).This condition satisfies the 68% confidence level condition for the measures of the velocity components.The angles with respect Table 5. Pulsar transverse velocities.The second column shows the velocities that are corrected for the peculiar velocity of the Sun and the Galactic differential rotation ( ⊥,LSR ).The third column shows the observed transverse velocities without any corrections ( ⊥,obs ).The fourth column shows the best previous reported velocities, some of which did not correct for the local standard of rest.The size of the error bars is showing 68% confidence interval from posterior probability distributions of velocities.
Figure 3. Histograms of transverse velocities for the isolated MSPs (red) and the binary MSPs (blue).PSR J1024−0719 is shown in black due to its extremely long orbital period (2-20 kyr (Kaplan et al. 2016)).PSR J1300+1240 is shown in olive as it is a planetary system (Wolszczan 1990).These two pulsars are not included in obtaining the velocity dispersion of binary MSPs as the evolutionary history is quite unlike that of the others in the sample.

Asymmetric drift in the Galaxy
An asymmetric drift is predicted to be observed in normal pulsars with characteristic ages greater than 10 7 yr (Hansen & Phinney 1997).
MSPs have characteristic ages greater than 10 8 yr and are expected to be old enough to have reached a virialized state; so they are also expected to demonstrate such a drift.Cordes & Chernoff (1997) predicted an asymmetric drift of 13 km s −1 for the case of a uniform surface density model, and 25 km s −1 for the case of an exponential surface density model, both with the MSP velocity dispersion of 60 km s −1 .In Figure 5, we can clearly see most of the MSPs moving opposite to the direction of Galactic rotation, which is the signature of the asymmetric drift.Utilising the rotational velocities,   , we obtained a median rotational velocity of 38( 16) km s −1 and the mean rotational velocity of 45(13) km s −1 .By removing the outliers (as defined in 5.3), we obtained the median value of 33( 14) km s −1 and the mean value of 38(11) km s −1 .Our values are higher than the 25(12) km s −1 obtained by Toscano et al. (1999).Our median value is closer to 25 km s −1 , the predicted value corresponding to the exponential surface density model of Cordes & Chernoff (1997) although the MSPs in our sample have higher mean velocities than theirs.

Millisecond Pulsars vs. Normal Pulsars
Soon after their discovery, a population study by Gunn & Ostriker (1970), based on data from just 41 young pulsars, argued that pulsars are probably born in a disk distribution with an initial scale height of ∼ 80 pc (similar to the normal OB stars), and move with average velocities of ∼ 100 km s −1 .In the first large scale study of pulsar proper motions, Lyne et al. (1982) studied a sample of 26 and showed that the pulsars with Galactic scale height greater than a few tens of parsecs tend to be moving away from the Galactic plane, probably due to a velocity kick at birth.They found pulsars had an rms velocity of about 210 km s −1 .Later, Lyne & Lorimer (1994)  youngest pulsars with proper motions, and found they had a mean space velocity at birth of 450(90) km s −1 .
By modelling the kinematics of the spatial distribution of MSPs, Cordes & Chernoff (1997) estimated that such pulsars on average receive a -velocity kick at birth of just 52 +17 −11 km s −1 , much less that than of the normal pulsars.They suggested that the rms speed of young pulsars is ∼ 5 time larger than that of MSPs, and a significant contribution to the observed -velocity (the velocity component along Galactocentric ) of MSPs originates from the diffusive processes that affect all old stars in the disk.
We plotted the histograms of the transverse velocities of the normal pulsars from the pulsar catalogue and the MSPs in Figure 7, to highlight how different the two distributions are.In this Figure, the top, middle, and bottom panels show the distributions of MSPs from this work, MSPs from the combination of this work and the pulsar catalogue, and normal pulsars from the pulsar catalogue, respectively.For the 87 normal pulsars in the pulsar catalogue with significant parallaxes and proper motions (> 3), we found a mean velocity of 246(21) km s −1 .Comparing the mean velocities showed that the normal pulsars seem to be faster than MSPs by a factor of ∼ 3.2.The characteristic age of MSPs is generally higher than that of normal pulsars, and accordingly, the old pulsars seem to be slower than young pulsars.If we restrict our attention to the pulsars with characteristic ages less than 3 Myr, we find a mean transverse velocity of 283(44) km s −1 , ∼ 3.7 times faster than the MSPs.
The velocities of the MSPs and normal pulsars as a function of Galactocentric  are shown in Figure 8.This figure demonstrates that pulsars with velocities < 100 km s −1 are more concentrated in the Galactic plane and are less scattered compared to the pulsars with velocities > 100 km s −1 .In addition, on average, normal pulsars seem to have a wider -distribution around the Galactic plane (i.e.larger  scale height) compared to the MSPs consistent with their higher velocities.This figure provides additional evidence for the results from the numerical studies of Bhattacharya & van den Heuvel (1991, and references therein); Tauris & Bailes (1996); Cordes & Chernoff (1997) that suggested MSPs have lower velocities compared to the young, long-period pulsars.The higher velocities of normal pulsars reinforce the idea that the high recoil velocity that pulsars receive during their birth in supernova explosions (Burrows 2013;Verbunt et al. 2017;Deller et al. 2019), and are in agreement with some of the core-collapse supernovae simulations performed by Wongwathanarat et al. (2012); Müller (2020).On the other hand, some of the low velocities of normal pulsars might be due to the weak supernova kicks for a sub-population of the pulsars (Willcox et al. 2021).Some caution in interpreting Figure 8 is required, as many selection effects are at play.

Velocity dispersions for millisecond pulsars
The dispersions of the 3 velocity components that Matthews et al. (2016) obtained from the analysis of NANOGrav nine-year MSP timing data set, after excluding outliers, are   = 46 km s −1 ,   = 40 km s −1 , and  z = 24 km s −1 .We obtained   = 63(11) km s −1 ,   = 48(8) km s −1 , and  z = 19(3) km s −1 .Our   and  z are thus comparable with their values.Also, they calculated the dispersion of velocity components, using the fit equations to the local stellar data provided by Aumer & Binney (2009), to be   = 34 km s −1 ,   = 22 km s −1 , and  z = 18 km s −1 for the characteristic age of  ∼ 5 Gyr, and   = 42 km s −1 ,   = 28 km s −1 , and  z = 24 km s −1 for  ∼ 10 Gyr.This is showing that our  z is consistent with the model fit corresponding to the characteristic age of  ∼ 5 Gyr, but our   and   are closer to the model fit corresponding to the characteristic age of  ∼ 10 Gyr.Both comparisons show that the MSPs are consistent with having been drawn from the old disk stellar population, and that they are subject to kick velocities at birth, as discussed in the next section.
The fact that the -components are the lowest is perhaps easiest to understand.MSPs with low -velocities spend more time near that galactic disk (and hence the Sun) and are preferentially detected by MSP surveys, plus as mentioned before, the -velocity is on average lower than the birth -velocity as it exhibits simple harmonic motion in the Galactic potential.Hobbs et al. (2005) studied the kinematics of 233 pulsars and found the mean two-dimensional speeds of 211(18) km s −1 for all pulsars in their sample, 307(47) km s −1 for pulsars with characteristic ages less than 3 Myr, 87(13) km s −1 for recycled pulsars, and 209(19) km s −1 for the normal pulsars with characteristic ages greater than 3 Myr.Our mean velocity of 78(8) km s −1 for MSPs is consistent with their mean two-dimensional speed of 87(13) km s −1 for recycled pulsars.Also, Hobbs et al. (2005) obtained the mean 2D speeds of 77( 16) km s −1 for seven isolated MSPs and 89(15) km s −1 for 28 binary MSPs.These are comparable with the mean velocities of 67(12) km s −1 and 77(9) km s −1 for the 16 solitary and 49 binary MSPs in our sample, respectively, excluding PSR J1024−0719.Johnston et al. (1998), by studying scintillation parameters for a sample of 49 pulsars, suggested that binary MSPs have higher velocities compared with the isolated MSPs.Toscano et al. (1999) calculated velocities for a sample of 23 MSPs and presented that, on average, the binary MSPs are one-third faster than isolated MSPs.However, Hobbs et al. (2005) reported that there is no significant difference between the mean velocities of both.Our results are also showing that the mean velocities of both are not significantly different.Using simulations, Tauris & Bailes (1996) predicted a mild inverse correlation between the recoil velocity (the post-explosion velocity that a pulsar receives at birth) and the orbital period of binary MSPs, depending upon the birth component masses of the binary, the separation between the two objects, and the evolution of the system during the common-envelope and mass transfer stages.Toscano et al. (1999) did not find any significant correlation with orbital period using a sample of 23 binary MSPs.In addition, Hobbs et al. (2005) looked for this correlation using the binary MSPs in their pulsar sample, but they did not find any either.We also investigated the relationship between transverse velocity and orbital period of binary MSPs in Figure 9, and measured the correlation coefficient to be −0.24(excluding the very long orbital period pulsar PSRs J1024−0719, the planet pulsar J1300+1240, and the eccentric MSP J2234+0611which probably has a complex evolutionary history) showing a weak anti-correlation.We need to keep in mind that many selection effects are influencing our observed MSP distribution.For instance, binary MSPs with short orbital periods and heavy companions might be less likely to appear in pulsar surveys unless acceleration searches have been undertaken.Nevertheless, Tauris & Bailes (1996) predicted an approximate velocity range of 25-270 km s −1 for binary MSPs and that those above 300 km s −1 should be very rare.This is borne out by our observations of our sample's transverse velocities.In the stan- dard model MSPs obtain their velocities as the vector sum of three distinct mechanisms.First, MSP progenitors, the binaries containing OB stars, are formed in large stellar nurseries where they obtain velocities of 10s of km s −1 due to stellar interactions.Second, when the neutron star is produced there is a Blaauw momentum kick to the binary as a result of the mass loss of the supernova (Blaauw 1961), and possibly an asymmetric kick imparted to the neutron star.To get very high velocities requires very compact binaries, large mass loss and asymmetric retrograde kicks that do not disrupt the orbit.These must be rare in MSP formation.As we only measure the 2D velocities the presence of some very low velocities is not unexpected and could be a projection effect, but MSP velocities can help constrain models for neutron star formation in population synthesis codes such as COMPAS (Riley et al. 2022) and STARTRACK (Belczynski et al. 2002).

Velocity distributions for millisecond pulsars
Any correlation between the minimum masses of companions of binary MSPs and their transverse velocities of binary MSPs was explored.We used the pulsar catalogue to find the minimum mass of the companions and plotted them as a function of transverse velocity in Figure 10.The minimum masses were calculated assuming the orbital inclination angle to be 90 • and the mass of the pulsar to be 1.35 M ⊙ .We measured the correlation coefficient between the minimum mass and the transverse velocity (excluding PSRs J1024−0719, J1300+1240, and J2234+0611) to be only −0.10, showing at best a very weak anti-correlation.Transverse velocity versus orbital period for binary MSPs.Velocities calculated using MSPs in this work are in blue and from the pulsar catalogue are in red.The correlation coefficient between velocity and binary period (excluding PSRs J1024−0719, J1300+1240, and J2234+0611 for reasons discussed in the text) was found to be −0.27,indicating a weak anti-correlation, similar to that predicted.

CONCLUSIONS
We have undertaken a study on the parallaxes and proper motions of 77 MSPs, observed with the MeerKAT radio telescope, over a time span of ∼ 3 years.Out of the 77 MSPs, 35 and 69 had significant parallaxes and proper motions (above 3), respectively.We calculated the transverse velocities for MSPs that had both significant parallaxes and proper motions.For the MSPs for which we did not have significant parallaxes, we used parallaxes taken from the literature.Pulsars with significant parallaxes and proper motions from the pulsar catalogue psrcat were also added to our MSP sample to produce the largest MSP velocity sample to date.We found that the transverse velocity of MSPs has a mean of 78(8) km s −1 and an rms scatter of 101 km s −1 .
The Lutz-Kelker effect on the derived transverse velocities was investigated for a sample of 35 pulsars with parallaxes.The median transverse velocities increased by ∼ 2% and ∼ 13% for parallaxes with > 10 and < 10 significance, respectively.The median of the whole population increased by ∼ 7%.
Out of our 69 MSPs with significant (> 3) proper motions, the long orbital period pulsar J1024−0719 had the highest transverse velocity of 278 +44 −32 km s −1 , and the 12.4 d orbital period binary pulsar PSR J1652−4838 had the lowest transverse velocity of just 5 +7 −2 km s −1 .
Although the mean velocity of 72(8) km s −1 for 49 binary MSPs is slightly faster than that of 67(12) km s −1 for 16 isolated ones, there was no evidence that their distributions are statistically significantly different.In comparison to the normal pulsars in the pulsar catalogue with the mean transverse velocity of 246(21) km s −1 , MSPs had a much lower mean transverse velocity of 78(8) km s −1 .
The velocity dispersions of the pulsars in the (cylindrical) Galac- tocentric velocity components (radial, rotational, and perpendicular) were measured to be   = 63(11) km s −1 ,   = 48(8) km s −1 , and  z = 19(3) km s −1 after removal of a small number of outliers.The lower -velocity component is likely a consequence of the fact that high velocity MSPs will spend less time near the Sun suppressing their representation in pulsar surveys.
The expected asymmetric drift was clearly seen in the rotational component of the velocities, and the mean value was found to be 38(11) km s −1 after removing outliers.This is consistent with the number predicted by Cordes & Chernoff (1997), who model the formation of MSPs in our galaxy.They predicted an asymmetric drift of approximately 25 km s −1 for their more realistic (exponential surface density) Milky Way model.
The substantially increased sample of MSPs with good parallaxes and proper motions presented in this study bolsters the case that MSPs arise out of the old disk of the Milky Way and, are subject to kick velocities at birth significantly smaller than those seen for young pulsars (Hobbs et al. 2005).We explored the distribution of MSP velocities versus orbital period and found a weak anti-correlation, in agreement with the predictions of Tauris & Bailes (1996), based on simulations of the recoil velocities of MSPs in various binary stellar evolutionary scenarios.
MSPs from the pulsar catalogue (Manchester et al. 2005) with > 3 significance in parallax into our MSP sample and calculated their DM-based distances.We excluded all MSPs in globular clusters as they have different origins and their velocities are dominated by their host clusters.
The top and bottom panels of Figure A1 show the comparison of parallax distances with DM distances obtained from the NE2001 and YMW16 models, respectively.In this Figure, points in red show a DM distance overestimation (i.e.MSPs are closer than what the models predict) and points in blue show a DM distance underestimation (i.e.MSPs are further than what the models predict).The ratios of  DM /  and   / DM are represented by the area of the circles in red and blue for the cases of overestimation and underestimation, respectively.The black circles in the legend of Figure A1 provide the standard circle sizes for the ratios of 10, 5, 2, and 1.We found the maximum underestimation based on the YMW16 model for PSR J1737−0811 with the distance ratio of   / DM ≈ 5.9 (considering the lower limit of its parallax distance), and the maximum overestimation was for PSR J1652−4838 with the distance ratio of  DM /  ≈ 10.5.
The 16%, 50%, and 84% of the underestimation and overestimation of the YMW16 distances ratios were obtained to be 1.4 +1.1 −0.3 and 1.2 +1.1 −0.1 , respectively.In addition, the underestimation and overestimation of the NE2001 distances were obtained to be 1.6 +1.1 −0.5 and 1.3 +1.8 −0.2 , respectively.Also, the distance ratios of   / DM were obtained for YMW16 model to be 1.1 +0.7 −0.3 and for NE2001 to be 1.2 +1.3  −0.4 .The comparison of distance ratios highlighted that there are obviously patches in Figure A1 that the DM-based distances are systematically overestimated or underestimated by the two models.

A2 Velocities
The top panel of Figure A2 presents the histogram of transverse velocities for the total sample that possesses parallax measurements.The mean transverse velocity is 78(8) km s −1 .The fastest moving MSP for these parallax-based velocities is the planet pulsar PSR J1300+1240 (Wolszczan 1990), which has a velocity of 330 +20 −18 km s −1 , and the slowest-moving, PSR J1652−4838, has a velocity of just 5 +7 −2 km s −1 .Using the code 4 provided by Antoniadis (2021), we explored whether the LK corrections might affect our velocities by breaking our sample into two groups, those with > 10 parallaxes and those with lower significance.The median velocities of the first group before and after the LK correction were 67 and 68 km s −1 respectively and the second group 43 and 49 km s −1 .Clearly, the LK correction raises the velocities of the second group and brings their median closer to that of the (more accurate) group, as might be predicted by the LK selection effect.Overall, the median velocities of the entire population before and after the LK correction were 46 and 50 km s −1 , an increase of ∼ 7 %.
To illustrate the effect of using DM-based distances rather than our preferred parallax-based distances, we have derived transverse velocities using DM-based distances obtained from NE2001 and YMW16 models.The middle and the bottom panels of Figure A2 present the histograms of transverse velocities for the MSPs that their velocities are calculated from YMW16 and NE2001 models, respectively.
Using the YMW16 model, the mean transverse velocity is 4 https://github.com/jantoniadis/GAIA_pulsars_eDR3/blob/main/GAIA_pulsars_eDR3.ipynbIt is curious to note that distance errors in YMW16 are responsible for a marked increase in the number of very low velocity MSPs (see Figure A2, lower two panels).In reality there are much fewer low velocity MSPs in the sample based on parallax distances only (Figure A2, upper panel), and this is merely an artifact of moving the MSPs too close to the Sun, reducing their apparent speeds.
We have performed a two-sample Kolmogorov-Smirnov test for comparing the continuous distributions of the parallax-based velocities and either of the two DM-based velocities to see if their macroscopic properties would differ.For the YMW16 model, the maximum absolute differences between the cumulative distribution functions were 0.11 and the p-value was 0.86.For the NE2001 model, the maximum absolute differences between the cumulative distribution functions were 0.17 and the p-values was 0.32.There is thus no statistical evidence that the parallax-based distances and DM-based distances differ that markedly in  ⊥ .
For the 35 pulsars with significant parallaxes (> 3), ∼ 2/3 of their YMW16 DM-based distances are more than 1 different from the parallax distances.This affects any velocities calculated from DM- based distances and makes the distribution of the velocities broader.
It is beyond the scope of this current work to examine each DM model and search for reasons why they predict the wrong distances but using either the NE2001 or YMW16 models for calculating distances for finding velocities of individual objects is subject to large uncertainties and for the rest of the paper we only consider the transverse velocities of pulsars that have reliable parallax measurements.It is very difficult to calibrate distance models at high galactic latitudes in the absence of pulsars with parallaxes.
This paper has been typeset from a T E X/L A T E X file prepared by the author.
using timing analysis of the Parkes pulsar timing array second data release, a parallax of 0.82(3) mas was obtained by Arzoumanian et al. (2018) using timing analysis of the NANOGrav 11-year data set.Using the VLBI technique, Chatterjee et al. (2009) measured a value of 0.95 +0.06 −0.05 mas.Desvignes et al. (2016) obtained a value of 0.90(3) mas, very close to our value.All measurements are in agreement with our value to within about 1.
) mas.Clearly, the timing of this pulsar presents challenges.PSR J1643−1224 has been shown to be located behind the HII region Sh 2−27 by Harvey-Smith et al. (2011) and undergoes annual DM variations.Its timing with MeerKAT also exhibits strong evi-

Figure 1 .
Figure 1.Celestial distribution of the 69 MSPs in this study, shown in Galactic longitude,  and Galactic latitude, .The curves are the paths of these MSPs from 5 Myr up until the present, and the current positions of these MSPs are shown with black stars.The effects of the peculiar motion of the Sun and the Galactic rotation curve have been subtracted from the observed proper motions.The black unit vectors are showing the proper motions with at least 3 significance, and the length of the brown vectors are showing the magnitudes of velocities with at least 3 significance in parallax detections.See the text for more details.

Figure 2 .
Figure 2. Comparison of our 23 parallax measurements with the measurements from literature.All of our measurements are consistent with previous work with only two exceptions, namely PSRs J1643−1224 and J2322−2650 that we discuss in Section 4.2.

Figure 4 .
Figure 4. Radial component of the cylindrical Galactocentric velocity (  ) versus distance from the Galactic centre () for MSPs whose radial velocity vectors are > 70 • from the line of sight direction.The dotted line indicates the assumed position of the Sun.Velocities calculated using MSPs in this work are in blue and from the pulsar catalogue are in red.The velocity dispersion of these MSPs is 63(11) km s −1 , if we exclude PSR J0125−2327 and J1300+1240 as outliers.

Figure 5 .
Figure 5. Rotational component of the cylindrical Galactocentric velocity (Δ  ) versus distance from the Galactic centre () for MSPs whose rotational velocity vector are > 70 • from the line of sight direction.The dotted line indicates the assumed position of the Sun.Velocities calculated using MSPs in this work are in blue and from the pulsar catalogue are in red.The velocity dispersion of these MSPs is this component is 46(8) km s −1 , excluding PSR J1909−3744 as an outlier.The signature of asymmetric drift is visible with the mean drift velocity of 39(11) km s −1 .

Figure 6 .
Figure 6. component of the cylindrical Galactocentric velocity (  ) versus distance from the Galactic centre () for MSPs whose -velocity vector are nearly perpendicular to the line of sight direction (> 70 • ).The dotted line indicates the assumed position of the Sun.Velocities calculated using MSPs in this work are shown in blue and those taken from the pulsar catalogue (psrcat) are shown in red.The velocity dispersion along the  direction is 19(3) km s −1 , excluding 6 outliers namely PSRs J1801−1417, J0610−2100, J1911−1114, J1227−4853, J1417−4402 and J1431−4715.It is clear that near the Sun (at  ≈ 8.12 kpc)   is small, and it becomes larger for MSPs further away from the Sun.As mentioned in the text multiple selection effects are at play.

Figure 7 .
Figure 7. Histograms of transverse velocities for our MSP sample (top panel and in blue), our MSPs and MSPs in the pulsar catalogue (middle panel and in gray), and the normal pulsars in pulsar catalogue that have > 3 significance in parallaxes (bottom panel and in red).This plot indicates how different the distributions of MSP velocities and the normal pulsar velocities are.

Figure 8 .
Figure 8. Transverse velocities versus Galactocentric  coordinate of pulsars.This plots shows the difference between the velocities of MSPs and normal pulsars.Velocities of MSPs from this work, MSPs from the pulsar catalogue, and normal pulsars from pulsar catalogue(Manchester et al. 2005) are shown in blue, red, and olive, respectively.Pulsars from the pulsar catalogue are shown if they have at least 3 significance in parallax.The mean velocity of the normal pulsars is 246(21) km s −1 , ∼ 3.2 times the mean MSP transverse velocity of 77(8) km s −1 .Also, the velocity dispersion of the normal pulsars is 194(15) km s −1 , ∼ 3 times the MSP velocity dispersion of 65(6) km s −1 .Pulsars with velocities < 100 km s −1 are often closer to the Galactic plane.
Figure9.Transverse velocity versus orbital period for binary MSPs.Velocities calculated using MSPs in this work are in blue and from the pulsar catalogue are in red.The correlation coefficient between velocity and binary period (excluding PSRs J1024−0719, J1300+1240, and J2234+0611 for reasons discussed in the text) was found to be −0.27,indicating a weak anti-correlation, similar to that predicted.

Figure 10 .
Figure 10.Transverse velocity versus minimum companion mass assuming orbital inclination of 90 • and the pulsar mass of 1.35 M ⊙ for binary MSPs.Velocities calculated using MSPs in this work are in blue and from the pulsar catalogue are in red.The correlation coefficient between velocity and minimum companion mass (excluding PSRs J1024−0719, J1300+1240, and J2234+0611 for reasons discussed in the text) was found to be −0.10,indicating almost no correlation.

Figure A1 .
Figure A1.The ratios of the parallax derived distances and DM derived distances.Top panel: distance ratios for the NE2001 model.Bottom panel: Distance ratios for the YMW16 model.The blue points indicate DM model distance underestimation (  / DM ), and the red points indicate distance overestimation ( DM /  ).The sizes of the black circles in the legend indicate the distance ratios.The black crosses in bottom panel show the parallaxes used for establishing the YMW16 model that we might expect to fit well.

Figure A2 .
Figure A2.Histograms of parallax-based velocities compared to the DMbased velocities.The transverse velocities are calculated using distances derived from the parallax measurements in the top panel, YMW16 DM model in the middle panel, and NE2001 DM model in the bottom panel.As discussed in Section A2, there is an excess of pulsars with very low velocities (lower two panels) for which distances have been derived from an ISM model, compared to the sample compiled from pulsars with parallaxes only (upper panel).This is an artifact of these models placing pulsars too close to the Sun.The velocity standard deviation of  ⊥ for the YMW16 model compared to the one for NE2001 model shows the number of low velocities in NE2001 model are higher.

Table 1 .
Uniform prior distributions of noise parameters are adopted in the ranges given in the second column

Table 2 .
Pulsar positions in Equatorial (, ), Ecliptic (,), and Galactic (,) coordinates to 3 significant figures.The digits in parentheses are showing 1 uncertainty in the last significant digits.The last column shows the reference epochs (in MJD) at which the positions are referred to.

Table 3 .
Pulsar Proper Motions in Equatorial, Ecliptic, and approximate Galactic Coordinates.The Galactic proper motions are derived from the Ecliptic ones.The time span (column 2) of observations for each pulsar is provided in yr.The previous best measurements of proper motions in Equatorial coordinates and the corresponding references are shown in the final 3 columns.Proper motions with > 3 and < 3 detections are listed in the upper and lower portions of the table, respectively.The list of the references are provided at the end of the table.