VLBA Astrometry of the Galactic Double Neutron Stars PSR J0509+3801 and PSR J1930$-$1852:
A Preliminary Transverse Velocity Distribution of Double Neutron Stars and Its Implications

Double neutron stars (DNSs) are valuable “laboratories” for probing theories of gravity, close binary star evolution, supernovae (SNe), and unveiling the composition of neutron stars (NSs). The gravitational-wave (GW) detection of the DNS merger event GW170817 (Abbott et al., 2017a), in conjunction with extensive observations across the electromagnetic spectrum (e.g. Abbott et al., 2017b; Mooley et al., 2018), has offered rich insights into the NS interior (Annala et al., 2018). A short $\gamma$-ray burst (SGRB) (Goldstein et al., 2017) coincident with the GW170817 merger event supports the widely believed association between DNS mergers and SGRBs (e.g. Fox et al., 2005; Coward et al., 2012; O’Connor et al., 2022; Fong et al., 2022), although most identified SGRBs are beyond the horizon of current GW detectors. Additionally, light curve monitoring and analysis after the GW170817 event suggest heavy elements were produced during the DNS merger (Drout et al., 2017; Watson et al., 2019), which reinforces the belief that DNS mergers are a prime source of r-process elements (Eichler et al., 1989; Korobkin et al., 2012). Nonetheless, it remains uncertain whether DNS mergers dominate the production of r-process elements.

Before the violent merger of a DNS is the phase of steady inspiraling. In this phase, Galactic DNSs that host an observable pulsar (hereafter referred to as DNS pulsars) have been extensively studied (e.g. Hulse & Taylor, 1975; Stairs et al., 2002; Faulkner et al., 2005; Kramer et al., 2006; Jacoby et al., 2006; Cameron et al., 2018; Stovall et al., 2018) with pulsar timing, a technique that measures and models the pulse time-of-arrivals (ToAs) from a pulsar. Despite being in shallower gravitational potentials compared to DNS mergers, DNS pulsars have offered some of the most stringent tests on gravitational theories in the strong-field regime with long-term timing observations (e.g. Wex, 2014; Fonseca et al., 2014; Weisberg & Huang, 2016; Ferdman et al., 2020; Kramer et al., 2021a). These tests are made by comparing observed post-Keplerian (PK) parameters, which quantify effects beyond a simple Keplerian model of motion, to the predictions of a specific gravitational theory, e.g. the general theory of relativity (GR).

There are two distinct channels of DNS formation: the predominant “isolated” channel and the “dynamical” channel (see Tauris & van den Heuvel, 2023 and references therein). The former channel gives rise to field DNSs, namely DNS systems born in isolated environments, while both channels can occur in dense regions such as globular clusters (Zevin et al., 2019). At the time of writing, only 19 known DNS pulsars and 2 suspected ones have been discovered from pulsar surveys (Sengar et al., 2022), including two found in globular clusters.

Precise proper motion and distance measurements for Galactic DNSs are essential for the studies related to DNSs. The orbital period derivative $\dot{P}_{\mathrm{b}}$ (or orbital decay) of a DNS system measured by pulsar timing is a PK parameter biased by the Shklovskii effect (Shklovskii, 1970), an apparent acceleration due to the pulsar’s transverse motion. $\dot{P}_{\mathrm{b}}$ is also affected by the radial acceleration caused by the overall gravity of the Galaxy (see Zhu et al., 2019 and references therein); this acceleration depends on our distance to the pulsar. Correcting the measured values for these effects requires an accurate measurement of the distance and proper motion. As the precision of the observed orbital decay $\dot{P}_{\mathrm{b}}^{\mathrm{obs}}$ improves as $t^{-2.5}$ (where $t$ refers to the observing time), the uncertain Shklovskii contribution to $\dot{P}_{\mathrm{b}}$ can quickly become the limiting factor of the $\dot{P}_{\mathrm{b}}$ test on gravitational theories (e.g. Deller et al., 2018; Ding et al., 2021).

Collectively, the measurements of DNS proper motions and distances allow a sample study of space velocities (also known as peculiar velocities or systemic velocities) of DNSs, which, combined with a DNS delay time distribution and galactic gravitational potentials, can be used to predict the spatial distribution of DNS merger sites within, or nearby, a host galaxy (e.g. Bloom et al., 1999; Voss & Tauris, 2003; Fong & Berger, 2013; Beniamini & Piran, 2016; Tauris et al., 2017; Vigna-Gómez et al., 2018; Andrews & Zezas, 2019; Beniamini & Piran, 2019; Zevin et al., 2022). Such a prediction, once made, can be compared to DNS mergers (Abbott et al., 2017b) and SGRBs localized with respect to their host galaxies, in order to probe the poorly defined selection biases (on both sides of the comparison), and identify exotic transient events. Most pinpointed SGRBs are found far from the centers of their host galaxies and star-forming regions (e.g. Fox et al., 2005; Fong & Berger, 2013; Tunnicliffe et al., 2014; O’Connor et al., 2022; Fong et al., 2022). Assuming SGRBs are predominantly generated by DNS mergers, the transverse space velocities $v_{\perp}^{\mathrm{SGRB}}$ of $\approx 20$–140 $\mathrm{km~{}s^{-1}}$ with a median of $\approx 60$ $\mathrm{km~{}s^{-1}}$ were inferred from the SGRB displacements from their expected birth sites (Fong & Berger, 2013), in conjunction with the indicative DNS delay time estimated by Leibler & Berger (2010). On the other hand, as the direct link between SGRB localizations and DNS formation theories, the space velocities of Galactic DNSs have been previously investigated with $\leq 9$ proper motion and distance measurements (Wong et al., 2010; Tauris et al., 2017; Haniewicz et al., 2021): the preliminary observational constraints suggest overall consistency with $v_{\perp}^{\mathrm{SGRB}}$.

Furthermore, a space velocity distribution of DNSs is crucial for understanding the source of excess r-process elements in ultra-faint dwarf galaxies (UFDGs; Ji et al., 2016; Hansen et al., 2017, 2020) where the escape velocities ($\sim 15$ $\mathrm{km~{}s^{-1}}$) can be easily surpassed due to the NS natal kicks received at supernovae (Beniamini et al., 2016; Safarzadeh et al., 2019). On the other hand, if a considerable fraction (e.g., $\gtrsim 10$%) of DNSs have a space velocity of $\leq 15$ $\mathrm{km~{}s^{-1}}$, then it is plausible that the r-process elements detected in UFDGs came from DNS mergers bound by the gravity of the host UFDGs. Otherwise, the r-process elements in UFDGs must be produced through fast-merging DNSs (Safarzadeh et al., 2019) or other mechanisms.

Finally, the DNS formation is not yet fully understood: for example, it remains debated under which conditions an electron-capture supernova (ECSN), instead of a conventional iron core-collapse supernova (CCSN), can give rise to the second-born NS (e.g. Podsiadlowski et al., 2004, 2005; Knigge et al., 2011; Jones et al., 2016; Giacobbo & Mapelli, 2019; Tauris & Janka, 2019). As ECSNe are expected to give smaller natal kicks than CCSNe (e.g. Kitaura et al., 2006; Dessart et al., 2006; Gessner & Janka, 2018), DNS space velocities can be used to probe the formation channels of field DNSs. Specifically, if both the CCSN and ECSN channels contribute to the field DNS productions, one would likely find a bimodal feature in the DNS space velocity distribution. However, the picture is more blurred because i) the narrow progenitor mass range leading to ECSNe depends on metallicity (Podsiadlowski et al., 2004; Siess & Lebreuilly, 2018), and UFDGs are very metal-poor (Fu et al., 2023), and ii) ultra-stripped SNe (Tauris et al., 2013, 2015), occurring when the last-formed NS is produced in DNS progenitor systems, may also produce NSs with small kicks in both CCSNe and ECSNe (e.g. Tauris et al., 2017; Müller et al., 2019). Additionally, a positive correlation has been proposed between the second-born NS masses and SN kick magnitudes, and thus space velocities of DNSs (Tauris et al., 2017; Burrows et al., 2023), and possibly also between second-born NS masses and orbital eccentricities, which can also be examined with more proper motion and distance measurements of field DNSs.

Pulsar timing and VLBI astrometry can both provide high-precision proper motions and potentially distances for DNSs. However, the lower timing precision attainable for the moderately recycled pulsars typically found in DNS systems (compared to fully recycled “millisecond” pulsars with periods of $\lesssim 10$ ms) means that long campaigns – often 10+ years – are required to detect proper motion and annual geometric parallax with adequate precision. VLBI astrometry, on the other hand, generally only requires two well-separated observations to detect a high-precision proper motion, and O(10) observations over a 1 year (or longer) span to measure annual geometric parallax with a precision of tens of microarcseconds, irrespective of the pulse period.

To date, only 5 DNSs have been well measured astrometrically with VLBI (Deller et al., 2018; Kirsten et al., 2014; Kramer et al., 2021a; Ding et al., 2021, 2023b), including 4 field DNSs (i.e. PSR B1913$+$16, PSR B1534$+$12, PSR J0737-3039A/B and PSR J1518$+$4904) and a DNS in the globular cluster M15. In this work, we carried out VLBI astrometry on two field DNSs — PSR J0509+3801 and PSR J1930$-$1852 (see Section 6 for their introductions). In Sections 2 through 4, we describe the VLBI observation strategy and data reduction procedures. In particular, a novel astrometric tactic utilising multiple reference sources is introduced in Section 3.1. We present the calculation of distances and transverse space velocities for 11 field DNSs in Section 5. The obtained sample of transverse space velocities is discussed thoroughly in Section 7, while individual DNS systems are discussed in Section 6. Throughout this paper, uncertainties are provided at 68% confidence, unless otherwise stated.

The in-beam calibrator data have been used to refine the phase calibration for the target source, eliminating temporal interpolation along with reducing spatial interpolation of calibrator solutions. In a normal procedure, self-calibration is performed on one (or sometimes multiple) in-beam calibrator sources (e.g. Deller et al., 2019); the acquired solutions are then applied to all sources in the target field. The multi-source self-calibration technique (Middelberg et al., 2013; Radcliffe et al., 2016) is especially useful when all in-beam calibrators are slightly too faint for self-calibration. On the other hand, in the opposite scenario where at least two in-beam calibrators are bright enough for self-calibration (or, alternatively, at least two relatively faint in-beam calibrators are identified around a sufficiently bright target; see, e.g., Section 4.1.2 of Ding et al., 2023b), it is possible to further enhance the astrometric precision by doubling the number of relative position measurements. Specifically, self-calibration solutions can be derived independently on each sufficiently bright in-beam calibrator. As the self-calibration solutions are applied to the target, the position of the target is measured with respect to the in-beam calibrator. Namely, one target position can be measured with respect to each independently calibrated in-beam calibrator.

However, systematic errors in the position offsets measured against different in-beam calibrators are generally at least partially correlated, due to the predominance of direction-dependent calibration errors in the systematic error budget. As an extreme and straightforward example, two in-beam calibrators situated at the same sky position would render almost perfectly correlated self-calibration solutions. Assuming direction-dependent terms dominate the self-calibration solutions, and that these terms change linearly with the sky position of the in-beam calibrator, the solutions obtained with two in-beam calibrators are largely independent of each other only when the two in-beam calibrators are in perpendicular directions as viewed from the target (and even in this case, other sources of systematic error, such as variations taking place between adjacent calibration solutions, remain correlated). Such a pair of in-beam calibrators are hereafter referred to as orthogonal in-beam calibrators, or simply orthogonal calibrators. The chance of having a pair of physical orthogonal in-beam calibrators is rather small. However, with the help of the 1D interpolation technique (e.g. Fomalont & Kopeikin, 2003; Doi et al., 2006; Ding et al., 2020b), the target can be measured with respect to a virtual calibrator (see Ding et al., 2020b for more explanations), whose effective position can be manipulated (either along the line connecting two self-calibratable in-beam calibrators or along the line connecting the out-of-beam main phase calibrator and a self-calibratable in-beam calibrator) to form an orthogonal in-beam calibrator pair with a physical in-beam calibrator. For brevity, we hereafter refer to this novel strategy of VLBI astrometry as the “2X” (read as “two cross”) strategy.

Interestingly, both DNSs possess at least two self-calibratable in-beam calibrators (see Table 1). Therefore, the 2X strategy is potentially applicable to the astrometry of both DNSs. Based on the brightness and compactness of the in-beam calibrators as well as their angular distances from the target, we adopted NVSS J051132$+$380913 and NVSS J193138$-$185258 as the primary in-beam calibrators for PSR J0509+3801 and PSR J1930$-$1852, respectively. In-beam calibrators other than the primary in-beam calibrators are referred to as secondary in-beam calibrators. Unfortunately, both NVSS J050958$+$381358 and NVSS J193102$-$185830 display a resolved jet feature in the right ascension direction. Any time dependence to the jet brightness profile along the jet direction may corrupt the parallax measurement (e.g. Deller et al., 2013), given that the parallax signature of either DNS is predominantly revealed in the right ascension direction. In such a case, the systematic error budget captured by the self-calibration solution might no longer be dominated by direction-dependent effects, breaking the underlying assumption of the 2X calibration strategy. Therefore, NVSS J050958$+$381358 and NVSS J193102$-$185830 are not used as secondary in-beam calibrators. Accordingly, we only implemented the 2X tactic on PSR J0509+3801 with its primary in-beam calibrator NVSS J051132$+$380913 and the secondary in-beam calibrator NVSS J050940$+$381301.

For each DNS, self-calibration was first performed on the primary in-beam calibrator; the acquired solutions were applied to all sources in the target field. In this way, the target position is measured with respect to the primary in-beam calibrator. After the application of in-beam phase calibration solutions to all target-field sources, we imaged both the target sources and the in-beam calibrator sources after dividing by the calibrator source model using natural weighting. From each image per source per epoch, we extracted a position and uncertainty using the image-plane fitting task JMFIT in AIPS. The positions so obtained are essentially anchored to the primary in-beam calibrator.

For PSR J0509+3801 only, we went one step further in order to fulfill the 2X strategy: this time PSR J0509+3801 needs to be phase-referenced to a virtual calibrator that forms orthogonal calibrators with the primary in-beam calibrator NVSS J051132$+$380913. In practice, we implemented 1D interpolation along the line connecting the NVSS J051132$+$380913 and NVSS J050940$+$381301 (see Ding et al., 2020b for the detailed procedure), then applied the solutions to only PSR J0509+3801. By placing the virtual calibrator 1.27 times further away from NVSS J051132$+$380913 compared to NVSS J050940$+$381301, the virtual calibrator and NVSS J051132$+$380913 form an orthogonal pair around PSR J0509+3801 (see Figure 1). In this instance, the virtual calibrator is further away from PSR J0509+3801 than NVSS J050940$+$381301 (which may not be the case given a calibrator plan of another target); this means that the systematic errors (of the target positions) resulting from direction-dependent propagation effects are accordingly amplified. Despite this small setback, the overall astrometric precision is expected to be improved (see Section 4).

We estimated the absolute reference positions of PSR J0509+3801 and PSR J1930$-$1852 in the same way as described in Section 3.2 of Ding et al. (2020a), except that, instead of the bootstrap uncertainties of the relative reference positions with respect to the primary in-beam calibrator, the Bayesian ones were adopted here as part of the error budget of the absolute reference positions. At the reference epoch MJD 57381, the absolute position of PSR J0509+3801 is $\alpha_{\mathrm{J0509+3801}}=05^{\rm h}09^{\rm m}31\fs 78871\pm 0.4\,\mathrm{mas}\pm 0.8\,\mathrm{mas}$, $\delta_{\mathrm{J0509+3801}}=+38\arcdeg 01^{\prime}18\farcs 0730\pm 0.6\,\mathrm{mas}\pm 0.8\,\mathrm{mas}$, while the absolute position of PSR J1930$-$1852 is $\alpha_{\mathrm{J1930-1852}}=19^{\rm h}30^{\rm m}29\fs 7156\pm 2.2\,\mathrm{mas}\pm 0.8\,\mathrm{mas}$, $\delta_{\mathrm{J1930-1852}}=-18\arcdeg 51^{\prime}46\farcs 164\pm 8.1\,\mathrm{mas}\pm 0.8\,\mathrm{mas}$ at MJD 57433. For both positions, the second error terms are indicative of the frequency-dependent core shift of the off-beam phase calibrator (Sokolovsky et al., 2011, also see the explanation in Section 3.2 of Ding et al., 2020a); the first error terms consist of 1) the Bayesian uncertainty of the relative reference position measured with respect to the primary in-beam calibrator, 2) the position scatter of the primary in-beam calibrator (across all epochs) with respect to the off-beam phase calibrator, and 3) the uncertainty of the off-beam phase calibrator position routinely updated in the Radio Fundamental Catalog⁴⁴4astrogeo.org/rfc/ (RFC).

Where parallax measurements are available, the parallax-based distances are adopted as the DNS distances $D$ (see Table 3). Otherwise, $D$ is approximated by the weighted mean of the two DM-based distances $d_{\mathrm{DM}}$, where the uncertainty of $D$ includes the indicative 20% fractional uncertainty as well as the standard deviation of the two $d_{\mathrm{DM}}$. We note that the DM-based distances might be subject to extra systematic uncertainties (see Section 6); nevertheless, they will get increasingly reliable with time, as $n_{\mathrm{e}}$ models become better constrained with growing number of independent pulsar distance measurements. The parallax-based distances of PSR J0509+3801 and PSR J1930$-$1852 were estimated in the same way detailed in Section 6 of Ding et al. (2023b), while the parallax-based distances of the four other field DNSs with well constrained parallaxes are directly quoted from respective publications (see Table 3 for the references).

Investigating the kinematics of Galactic DNSs ideally requires the knowledge of the 3D space velocity, $v^{\mathrm{birth}}$, of each DNS system immediately after its birth (which is inherently linked to the formation history and the SN properties of the system, e.g. Tauris et al., 2017). Here, “space velocity” refers to the peculiar velocity with respect to the local stellar neighbourhood of the DNS (see discussions in Section 6 of Ding et al., 2023b). In theory, $v^{\mathrm{birth}}$ of a DNS system can be derived from its current full kinematic information (including the 3D velocity and the 3D position) and the age, $t_{\mathrm{age}}$, of the second-born NS, by applying a Galactic potential. This pathway to $v^{\mathrm{birth}}$ is, however, impractical, due to a highly uncertain value of $t_{\mathrm{age}}$ and the absence of the radial velocity, $v_{\mathrm{r}}$, for all known Galactic DNSs. With measured proper motion and 3D position (from sky coordinates and distance) of a DNS, we can unfortunately only constrain its current transverse velocity, i.e. the 2D velocity along the plane of the sky, and thus only infer its 2D peculiar velocity, $v_{\perp}$, relative to the stellar neighbourhood of the DNS system.

With no better option, this work uses $v_{\perp}$ as a proxy of $v_{\perp}^{\mathrm{birth}}$, i.e., the transverse velocity at birth, to probe the DNS kinematics. This approach hinges on the underlying assumption that $v_{\perp}$ is largely correlated with $v_{\perp}^{\mathrm{birth}}$, which is the first caveat to our analysis of DNS kinematics. Conceivably, this assumption holds more truth for DNSs with low space velocities, as their space velocities change less over time in the Galactic potential. Therefore, we expect that a $v_{\perp}$ distribution compiled from a large number of DNSs should be close to the DNS $v_{\perp}^{\mathrm{birth}}$ distribution at the low-space velocity end, while being more dispersed than the $v_{\perp}^{\mathrm{birth}}$ distribution at the high-space velocity end. Assuming natal kicks are equally possible in all directions (isotropic kicks), it is easy to calculate that on average $v_{\perp}^{\mathrm{birth}}=\left(\pi/4\right)\cdot v^{\mathrm{birth}}$. Accordingly, we approximate the $v^{\mathrm{birth}}$ distribution by $\left(4/\pi\right)\cdot v_{\perp}$ in this work. To carefully justify the proposed $v_{\perp}$ to $v_{\perp}^{\mathrm{birth}}$ mapping (as well as the coefficient $4/\pi$) requires much more data and a dedicated population analysis, which is, however, not possible at this stage.

The $v_{\perp}$ of the DNSs were calculated in a consistent manner following Section 6 of Ding et al. (2023b). The mathematical recipe of the calculation has been detailed in Verbunt et al. (2017). We evaluated the $v_{\perp}$ uncertainties using Monte Carlo simulations, where proper motions and distances with asymmetric uncertainties were approached with split normal distributions (e.g. Possolo et al., 2019). In the $v_{\perp}$ estimation, we approximated the motion of the stellar neighbourhood of each DNS by a circular Galactocentric motion parallel to the Galactic plane. This treatment is conceivably less reliable for DNSs $\gtrsim 1$ kpc away from the Galactic plane. Any anticorrelation (or correlation) between the Galactic height $|z|$ and $v_{\perp}$ would suggest that the approximation is oversimplified and inappropriate. Though no correlation between $|z|$ and $v_{\perp}$ is identified (see the upper panel of Figure 3) with an insignificant Pearson correlation coefficient of $0.35^{+0.23}_{-0.24}$ (let alone a high false alarm probability of 0.29, see Table 4), we cautiously note that extra systematics of $v_{\perp}$ are likely introduced by the approximation for the relatively high-$|z|$ DNSs. Due to the overly uncertain distance of PSR J1757$-$1854 (see Table 3, or Figure 13 of Cameron et al., 2023), we did not infer $v_{\perp}$ for PSR J1757$-$1854. Neither did we use it in the sample studies discussed in Section 7. The other 11 DNSs used in the sample studies are hereafter referred to as the 11 DNSs.

PSR J0509+3801 was discovered in the Green Bank Northern Celestial Cap Pulsar Survey at 350 MHz (Stovall et al., 2014), and later identified as a DNS system in a highly eccentric ($e=0.586$) 9.1-hr orbit (Lynch et al., 2018). The pulsar mass and the companion mass of the system are well determined thanks to the significant measurements of the periastron advance and the gravitational redshift (Lynch et al., 2018). Among the 12 field DNSs having proper motion measurements, PSR J0509+3801 possesses the second highest orbital eccentricity. Therefore, the relatively precise determination of $v_{\perp}$ for PSR J0509+3801 (see Table 2) is crucial for probing the $v_{\perp}$-to-$e$ correlation (see Section 7.1). By comparing the parallax-based distance of PSR J0509+3801 to the two DM distances, we find that, at $\approx 2\,\sigma$ confidence, both DM distances (based on the NE2001 and the YMW16 $n_{\mathrm{e}}$ models) are underestimated.

PSR J1930$-$1852 is the sixth pulsar discovered by students attending the Pulsar Search Collaboratory (Rosen et al., 2013) from the drift-scan observations made with the Green Bank Telescope at 350 MHz (Boyles et al., 2013; Lynch et al., 2013). PSR J1930$-$1852 is an exceptional DNS system: it challenges the DNS formation theory with the longest $P_{\mathrm{s}}$ (185.5 ms) and the widest ($P_{\mathrm{b}}=45$ d) orbit among the Galactic field DNSs (Swiggum et al., 2015, also see Table 3). Accordingly, the astrometric determination of PSR J1930$-$1852 is important for investigating the $v_{\perp}$-to-$P_{\mathrm{s}}$ and $v_{\perp}$-to-$P_{\mathrm{b}}$ correlations (see Section 7.1). Unfortunately, we did not achieve significant parallax measurement for PSR J1930$-$1852 due to its radio weakness: random errors dominate the error budget of the position series. In this regard, the parallax precision of PSR J1930$-$1852 can be potentially enhanced by a factor of $\gtrsim 4$ in future with high-sensitivity VLBI arrays like the High Sensitivity Array.

The 3D systemic space velocity, $v^{\mathrm{birth}}$, at DNS birth, right after the second SN, is acquired from the combination of the anisotropic nature (i.e. natal momentum kick) and instantaneous mass loss of the SN producing the second-born NS. Therefore, it has been proposed that $v^{\mathrm{birth}}$ is likely correlated with a few parameters which include the orbital eccentricity, $e$, the orbital period, $P_{\mathrm{b}}$, and the companion mass $m_{\mathrm{c}}$ (e.g. Tauris et al., 2017; Haniewicz et al., 2021). Despite the caveats of $v_{\perp}$ estimation (see Section 5.2), we investigate if $v_{\perp}$ is correlated with the aforementioned parameters. In addition to $v_{\perp}$, another important indicator of DNS kinematics is the DNS distance $|z|$ from the Galactic plane, which has also been refined with VLBI astrometry of DNSs. Assuming DNS systems were born in the Galactic plane and received natal kicks in random directions, DNSs with large space velocities are more likely found at high $|z|$. To comprehensively examine the correlation between the two DNS kinematics indicators (i.e. $v_{\perp}$ and $|z|$) and other parameters, we carried out linear and power-law correlation analyses on 7 parameters, including $v_{\perp}$, $|z|$, $e$, $P_{\mathrm{b}}$, spin period (of the first-born and recycled NS) $P_{\mathrm{s}}$, $m_{\mathrm{c}}$ and pulsar mass, $m_{\mathrm{p}}$. The calculated Pearson correlation coefficients $\rho$ and their false alarm probabilities (or p-values) $P_{\mathrm{F}}$ are summarized in Table 4, where the two DNS systems (i.e., PSRs J1930$-$1852 and J1411$+$2551) that have not had significant measurements of the pulsar mass are excluded from the calculation of mass-related Pearson correlation coefficients.

A correlation is normally considered tentative when $0.5\lesssim|\rho|\lesssim 0.75$ and $P_{\mathrm{F}}<0.05$, while being strong if $\rho\gtrsim 0.75$ and $P_{\mathrm{F}}<0.05$. Unsurprisingly, two of the examined parameters (i.e., $P_{\mathrm{s}}$ and $|z|$) are well correlated with $P_{\mathrm{b}}$. These correlations are thought to have an origin in binary stellar evolution. For example, Tauris et al. (2015, 2017) demonstrated that the wider the pre-SN binaries orbit is, the less efficient the final Case BB mass transfer to the first-born NS is (the reason being a shorter lifetime prior to the core collapse of more evolved helium star donors, i.e. the progenitors of the second-born NSs). Therefore, after the second SN, wider-orbit DNS systems have, on average, slower spinning recycled NSs compared to more tight-orbit DNS systems. This is confirmed from our analysis which yields a Pearson correlation coefficient of $\rho=0.92$. Due to the unique kick magnitude and direction of each DNS system, some scatter is expected in the final correlation between $P_{\mathrm{b}}$ and $P_{\mathrm{s}}$. Tauris et al. (2015, 2017) also argued for weak correlations between $e$ and both $P_{\mathrm{b}}$ and $P_{\mathrm{s}}$. We find that the inclusion of PSR J1930$-$1852, which has exceptionally long $P_{\mathrm{b}}$ and $P_{\mathrm{s}}$ (Swiggum et al., 2015) among Galactic field DNSs, plays an essential role in dismissing the $e$-to-$P_{\mathrm{b}}$ correlation. But even without PSR J1930$-$1852, a tentative $e$-to-$P_{\mathrm{b}}$ correlation is not revealed by our small sample. On the other hand, a tentative logarithmic correlation between $e$ and $P_{\mathrm{s}}$ is identified, even when PSR J1930$-$1852 with the long $P_{\mathrm{s}}$ of 185.5 ms is included.

Between the two kinematics indicators, $v_{\perp}$ and $|z|$, we did not find any significant correlation, which is the precondition of the $v_{\perp}$ estimation (see Section 5.2). However, $v_{\perp}$ and $|z|$ apparently show correlations with the other parameters: $v_{\perp}$ is tentatively correlated with $e$ and $m_{\mathrm{c}}$ (see Figure 3), while $|z|$ is well correlated with $P_{\mathrm{b}}$. In addition, both kinematics indicators are correlated with $P_{\mathrm{s}}$ (see Figure 3). Large SN kicks will, in general, result in relatively large eccentricities. Hence, the former correlation between $v_{\perp}$ and $e$ is perhaps not surprising. There are several components that affect the SN outcome: large kicks will disrupt the binaries. Because of the strong dependence on the stochastic kick direction in the second SN, as well as the different masses of the naked helium stars resulting from the previous common-envelope phase, it is non-trivial to directly link kinematic parameters ($v_{\perp}$ and $|z|$) of DNS systems to their orbital parameters and component masses. A proposed correlation between the SN kick velocity magnitude and the resulting NS mass (Tauris et al., 2017) would translate into a (weaker) correlation between e.g. $v_{\perp}$ or $|z|$ and $m_{\mathrm{c}}$. Our analysis confirms the postulated correlation between $v_{\perp}$ and $m_{\mathrm{c}}$, while finding no correlation between $|z|$ and $m_{\mathrm{c}}$. On the other hand, any correlation between $|z|$ (or $v_{\perp}$) and $P_{\mathrm{s}}$ is somewhat unexpected and non-trivial to explain.

We note that great caution must be taken with postulating correlations involving parameters that evolve with time, such as $P_{\mathrm{s}}$ and $e$. To properly analyse the DNS population, time evolution must be included, e.g. Tauris et al. (2017). Future analysis with more DNS systems will likely improve the statistical significance of the tentative correlations suggested in this work. Given no indication of insufficient $e$ or $P_{\mathrm{s}}$ coverage in the pulsar search programs (that led to the discoveries of the Galactic DNSs), no selection bias is uncovered for the $v_{\perp}$ distribution. Among the four above-mentioned correlations between the two kinematics indicators and other parameters, linear correlations are either preferred over or comparable to the power-law counterparts by $\rho$ and $P_{\mathrm{F}}$ (see Table 4). This conclusion can also be drawn for the whole Table 4, with the {$e,P_{\mathrm{s}}$} pair being the only exception.

As noted in Section 5.2, a $v_{\perp}$ distribution based on a large sample of astrometrically measured DNSs is expected to better approximate the $v_{\perp}^{\mathrm{birth}}$ distribution at the low velocity end, while being more dispersed than the $v_{\perp}^{\mathrm{birth}}$ counterpart at the high velocity end. With this in mind, we compiled the $v_{\perp}$ results in Table 3 to a DNS $v_{\perp}$ distribution using Monte Carlo re-sampling. Specifically, we assume that all of the 11 field DNSs equally represent the whole population of field DNSs, and randomly drew 10,000 values of $v_{\perp}$ for each DNS from a split normal distribution. Subsequently, the $11\times 10,000$ randomly drawn values of the DNSs were concatenated together. From this re-sampled $v_{\perp}$ chain, we estimated $v_{\perp}$ of field DNSs to be $60^{+85}_{-44}$ $\mathrm{km~{}s^{-1}}$(see Figure 3), hence the 3D space velocities of field DNSs should be $\sim 4/\pi\times 60$ $\mathrm{km~{}s^{-1}}$$=76$ $\mathrm{km~{}s^{-1}}$ (see Section 5.2). Here, the median of the $v_{\perp}$ chain is adopted as the $v_{\perp}$ estimate, while the central 68% of the $v_{\perp}$ chain is used as the $v_{\perp}$ uncertainty interval. To mitigate binning effects, we smoothed the histogram of the $v_{\perp}$ simulations with kernel density estimation provided by scipy.stats.gaussian_kde. The resultant dashed curve (in the bottom panel of Figure 3) represents the estimated probability density distribution (PDF) of $v_{\perp}$.

In this work, we measured proper motions and parallaxes for two DNS systems PSRs J0509+3801 and J1930$-$1852 with VLBA observations. In a consistent manner, we estimated the two kinematic indicators, the transverse space velocity $v_{\perp}$ and the vertical height $|z|$, for all Galactic DNSs astrometrically measured (see Section 5). With this small sample of 11, we investigated the correlation between the kinematic indicators and 5 other parameters, including the orbital period $P_{\mathrm{b}}$, the spin period $P_{\mathrm{s}}$, the orbital eccentricity $e$, the companion mass $m_{\mathrm{c}}$ and the pulsar mass $m_{\mathrm{p}}$. Tentative correlations were identified between $v_{\perp}$ and $e$, $v_{\perp}$ and $P_{\mathrm{s}}$, $v_{\perp}$ and $m_{\mathrm{c}}$, $|z|$ and $P_{\mathrm{s}}$, and $|z|$ and $P_{\mathrm{b}}$ (see Section 7.1). By combining the $v_{\perp}$ estimates, we obtained the preliminary $v_{\perp}$ distribution of DNSs (see Section 7.2), which is used to pursue three scientific studies. Firstly, we found that the $v_{\perp}$ distribution is consistent with the indicative transverse velocities $v_{\perp}^{\mathrm{SGRB}}$ of short gamma-ray bursts (SGRBs) evaluated by Fong & Berger (2013), which suggests that the underlying DNS delay time provided by Leibler & Berger (2010) is reasonable. Secondly, according to the $v_{\perp}$ distribution and the typical escape velocity of ultra-faint dwarf galaxies (UFDGs), $\gtrsim 11$% of the DNSs born in UFDGs are gravitationally bound to UFDGs. We argue that a significant fraction of the DNSs gravitationally bound to UFDGs would eventually contribute to the r-process enrichment of the UFDGs. Therefore, DNS mergers cannot be ruled out as the source of r-process enrichment in UFDGs. Finally, we looked into the prospects of using the $v_{\perp}$ distribution to probe the formation channels of DNSs. Our simulation suggests that astrometric measurements of $\sim 120$ additional isolated DNSs were needed to rule out unimodality of the $v_{\perp}$ distribution at 90% confidence, if the $v_{\perp}$ distribution shown in Figure 3 is true.

1. (i)

   The availability of the three components — the $v_{\perp}$ distribution (directly approached with Galactic DNSs), the SGRB displacements from the expected birth sites in the host galaxies (e.g. O’Connor et al., 2022; Fong et al., 2022), the DNS delay time distribution (constrained by Galactic DNSs, SGRB host galaxies and GW events, e.g. Andrews & Zezas, 2019; Beniamini & Piran, 2019; Zevin et al., 2022), promises a new research opportunity: by combining the knowledge of the three components, one can uncover the selection biases (in the three components) to unprecedented details, and achieve the most reliable posterior distributions for the three components.

2. (ii)

   As mentioned in Section 7.2.3, ongoing and future relativistic binary pulsar search programs using high-sensitivity radio telescopes (or tied arrays) will drastically increase the number of Galactic DNSs, which will eventually enable the multimodality test of the $v_{\perp}$ distribution, and refine the estimation of the fraction of DNSs gravitationally bound to UFDGs (see Section 7.2.2). Meanwhile, high-sensitivity VLBI observations of DNSs will not only increase the number of DNSs well measured astrometrically, but also achieve higher astrometric precision for the DNSs.