Identifying Companions in Pulsar Binary Systems via Gaia Data

1. Introduction

Pulsars are generally recognized as rapidly rotating neutron stars. Millisecond pulsars (MSPs), characterized by even faster rotation periods, are widely considered to be a class of recycled pulsars. These objects originate in binary star systems and have undergone a process of angular momentum transfer through mass accretion from their companion stars. When mass from the companion star escapes the Roche lobe under the star’s gravitational influence and flows towards the neutron star, it forms an accretion disk that releases immense energy. During this process, the angular momentum of the accreting material is transferred to the neutron star, accelerating its rotation to millisecond periods. Meanwhile, the strong magnetic field is temporarily weakened by the accreted material. At this stage, the neutron star typically manifests as a low-mass X-ray binary (LMXB) [1,2,3], emitting accretion energy in the form of X-ray radiation. As the accretion rate declines, the magnetic field re-emerges and activates radio pulsations, forming canonical MSPs. Compared to non-recycled pulsars, MSPs emit more frequent pulses with narrower profiles, enabling more precise measurements of pulse times of arrival (ToAs) [4]. For example, The timing precision for PSR J0437−4715 reaches 30 nanoseconds [5].

MSPs act as cosmic clocks and have extremely high spin rates. The tiny residuals in ToAs may contain gravitational wave signals. Using pulsar timing arrays (PTAs), we can detect gravitational waves in the nanohertz range, which are emitted by galaxy clusters containing supermassive black hole binaries. This ongoing exploration has significantly advanced the observational discovery of MSP populations [5]. The classification of spider pulsars is based on two categories, redback and black widow, which are determined by the mass and evolutionary state of their companion stars. A black widow system consists of an MSP and a very low-mass companion star with a mass of less than 0.1 $M_{\odot }$[6], from which pulsar winds strip away surface material, causing radio eclipses. These systems typically have orbital periods of less than 24 hours [7], with radio pulse profiles that vary significantly with frequency. Redback systems also comprise an MSP and a companion star with a higher mass (approximately 0.2–0.4 $M_{\odot }$) [6], likely representing a transitional phase between LMXBs and ordinary MSPs. They exhibit strong X-ray emission, and optical observations suggest that the companion star may be close to filling its Roche lobe.

Due to the extremely weak optical radiation emitted by pulsars, it is difficult to observe them directly in the optical wavelength range. However, pulsar binary systems present a different situation: their companion stars are typically much easier to detect in the optical range. The Gaia Data Release, published by the European Space Agency (ESA) in the form of incremental updates to the Gaia mission, provides a high-precision astronomical dataset containing positional, distance and spectroscopic information for billions of stars [8]. By providing astrometric parameters such as positions, proper motions and parallaxes with unparalleled accuracy, Gaia data have greatly enhanced the precision of astrometric measurements of pulsars by observing their companions in binary systems. When combined with radio observations of pulsars, this allows for more accurate constraints on their space velocities and orbital dynamics — critical factors for PTA in the detection of gravitational waves. More precise information on the positions and distances of pulsars can reduce noise in timing residuals more effectively, thereby enhancing the ability to identify ultra-low-frequency gravitational wave signals.

In this study, we use Gaia Data Release 3 (Gaia DR3) to search for pulsar companions. Sample data were obtained from the Australia Telescope National Facility (ATNF) Pulsar Catalogue (https://www.atnf.csiro.au/research/pulsar/psrcat/). Section 2 outlines our sample selection and methodology. Section 3 presents our findings. Section 4 provides a discussion and summary of the results.

2. Data and Methods

We obtained pulsar data from the ATNF Pulsar Catalogue [9], which contains a total of 3,781 pulsars at present. For this research, our focus was exclusively on pulsars located in the galactic field, with those found in globular clusters being excluded. Because Gaia faces significant challenges in distinguishing the true companions of pulsars due to the extremely high stellar density in globular clusters [10]. Moreover, pulsars within globular clusters listed in the ATNF Pulsar Catalogue are usually detected at high frequencies, and the astrometric precision is far less precise than that achieved by Gaia. Hence, the cross-identification of pulsars in globular clusters is less reliable. As our primary interest lies in identifying astrometrically close pairs, we applied an initial screening condition of positional accuracy better than 0.5 arcseconds, reducing the sample size to 1,527.

Our study uses Gaia DR3 data (reference epoch J2016.0) and requires positional alignment with the ATNF Pulsar Catalogue. We propagate the positions of pulsars to the Gaia epoch using their proper motions in right ascension and declination, as well as their parallax measurements from the ATNF Catalogue. We use Equation (1) to align the positions of the pulsars:

$$\begin{array}{c}\hfill \Delta \alpha ={\mu }_{\alpha }\Delta t+\pi f_{\alpha }(\alpha ,\delta ,\Delta t)\\ \hfill \Delta \delta ={\mu }_{\delta }\Delta t+\pi f_{\delta }(\alpha ,\delta ,\Delta t)\end{array}$$

(1)

In the equation, $\alpha$ and $\delta$ represent the right ascension and declination of the pulsar, respectively, while ${\mu }_{\alpha }$ and ${\delta }_{\delta }$ denote the proper motions in the two directions. $\pi$ is the parallax of the pulsar. $f_{\alpha }$ and $f_{\delta }$ are the parallax displacements of the pulsar in right ascension and declination, which can be expressed as:

$$\begin{array}{c}\hfill f_{\alpha }=sin\lambda cos\alpha cos\epsilon -cos\lambda sin\alpha \\ \hfill f_{\delta }=sin\lambda sin\epsilon cos\delta -sin\lambda cos\epsilon sin\alpha sin\delta -cos\lambda cos\alpha sin\delta \end{array}$$

(2)

Here, $\lambda$ represents the solar ecliptic longitude, and $\epsilon$ denotes the obliquity of the ecliptic. We then used the TOPCAT tool [11] to cross-match our sample with Gaia DR3. TOPCAT is primarily designed for processing and analyzing tabular data, particularly astronomical catalogues. It facilitates cross-matching between catalogues from various sources, based on the coordinates of celestial objects. Studies conducted by Antoniadis [12] and Liu et al. [13] have shown that the angular separations between pulsars and their optical companions in the Gaia data are usually small. Based on this, we conduct a cross-match within a $3^{\prime \prime }$ radius, yielding 615 matched objects. We then apply the following criteria: (1) the targets must have parallax measurements in Gaia DR3 and (2) the Gaia-measured proper motions in right ascension and declination must agree with the ATNF values within 3 $\sigma$ error margins.

An alternative approach is to use the photometric bands from the Gaia data for selection purposes:

$$\begin{array}{c}\hfill M_G\le 2.7(G_{BP}-G_{RP})+10.5\\ \hfill M_G>3.7(G_{BP}-G_{RP})+2.4\end{array}$$

(3)

Here, $M_G$ represents the absolute magnitudes of the pulsars. Meanwhile, $G_{BP}$ and $G_{RP}$ denote the mean magnitudes in Gaia’s BP and RP photometric bands. We use the formula $M_G=m_g-5{log}_{10}d+5$ to roughly estimate the absolute magnitude of pulsars. Here, $m_g$ and d represent the apparent magnitudes of the pulsars and the distance to the pulsars (in parsecs), respectively. However, the absolute magnitudes calculated using this method are not accurate. For instance, the known binary system PSR J0437−4715, paired with a helium white dwarf with a mass of $\sim 0.2$$M_{\odot }$ [14], yields an absolute magnitude estimate of between 6.7 and 13.6 using Equation (3). However, Antoniadis et al. [12] previously calculated its absolute magnitude to be approximately 14.4. If this method were used as a selection criterion, PSR J0437−4715 (and similarly PSR J1024−0719) would erroneously be excluded from our sample, despite being confirmed binaries. These obtained results are clearly incorrect. In Section 4, we will analyse the limitations of this method and present the correct formula for calculating absolute magnitude. This further emphasizes the reliability of using Gaia proper motions and parallaxes for sample selection.

3. Results

Using our methodology, we identified a total of 58 qualifying objects and 3 candidate objects. Figure 1 shows the spatial distribution of 58 pulsars. Red triangles represent redback systems, green squares represent black widow systems, black circles represent newly identified pulsars, and blue pentagrams represent non-spider pulsars. Based on their companion types, we classify identified Gaia companions and candidate objects as main-sequence stars, neutron stars, white dwarfs, ultra-light companion or planet. Among them, there are 23 spider pulsars, accounting for approximately 40%; main-sequence stars and white dwarfs account for about 14% and 12%, respectively; neutron stars and ultra-light companion stars together make up roughly 3% (Note: For example, PSR J1048+2339 is a spider pulsar with a main-sequence star as its companion, and we classify it as spider pulsars). As can be seen clearly in Figure 1, the newly identified pulsars are predominantly concentrated in the southern hemisphere, which further expands the spatial distribution of pulsars in this region. This is primarily due to the results obtained from Meer Karoo Array Telescope (MeerKAT) observations. Among the 9 newly identified pulsars, 8 are redback pulsars, which is significant for studying the evolution of MSPs and the composition of their companion stars. Table 1 provides key information on all 58 astrometric pairs, including Gaia DR3 designations, the angular separation between ATNF and Gaia positions, Gaia parallax values, proper motion from the ATNF, proper motion from Gaia, spin period, orbital period, and apparent magnitude. Given the continuous updates to the Gaia data employed, particularly the variations in the precision of parallax and proper motion measurements compared to those at the time of previous identifications, we have compiled the parameters for all 58 pulsars. Figure 2 displays the parallax values obtained from the ATNF and Gaia for the newly identified pulsars. The proper motion discrepancies of these pulsars are shown in Figure 3. We can observe that these pulsars exhibit small differences in proper motion in both directions. This further supports their identification as pulsars with companion stars. The next part of this section will provide detailed information on each of the 9 newly confirmed pulsars and the 3 targets currently designated as temporary candidates.

3.1. Main-sequence Star

PSR J2108+4516: We identified this pulsar in Gaia DR3 as source 2162555482829978496, with an angular separation $\theta$=0.103 arcseconds. The radio pulses of J2108+4516 were first detected by CHIME/FRB in 2018 [15], and subsequent observations by CHIME/Pulsar confirmed an orbital period of 269 days. Eclipsing events occur between orbital phases 0.06-0.76, lasting 33%-69% of the orbital period, potentially caused by a protoplanetary disk around the companion star or an extremely dense stellar wind [15]. To achieve precise localization, Andersen et al. [15] conducted observations using the Very Large Array (VLA), clearly identifying the companion as EM*UHA 138 - a bright star located in the North America Nebula region of the Cygnus constellation.

PSR J1036−4353, J1628−3205, and J1803−6707 have main-sequence star companions and belong to redback systems, which will be described in Section 3.5.

3.2. Neutron Star

PSR J2150+3427: In their study, Wu et al. [16] proposed two possibilities for the companion type of PSR J2150+3427: a double neutron star (DNS) system or a neutron star–white dwarf (NS–WD) system. Its orbital period of 10.592 days, eccentricity of 0.601 and total system mass of 2.59 ± 0.13 $M_{\odot }$ align with known DNS parameters. Meanwhile, the lower limit of the companion mass (>0.98 $M_{\odot }$) overlaps with the white dwarf mass range in certain NS-WD systems (e.g. B2303+46) [16]. However, since no second neutron star signal has been detected and key observational data such as the Shapiro delay are lacking, neither possibility can be fully ruled out. Future high-precision measurements of relativistic effects and improved distance modelling are required to definitively determine the companion type. Here, we list Gaia source 1948939718066589952 as a candidate.

3.3. White Dwarf

PSR J1910+1256: PSR J1910+1256 is a binary system consisting of a MSP and a low-mass white dwarf companion with a minimum mass of ∼0.18 $M_{\odot }$ [17]. Mingarelli et al. [18] derived the distance by combining PTA parallax measurements with Gaia companion parallax data. In Gaia DR2, source 4314046781982561920 was identified as the companion; however, its parallax measurement yielded a negative value of −0.42 ± 0.80 mas with a low signal-to-noise ratio of approximately 0.5. This anomaly persists in Gaia DR3 and requires verification with future data releases. We concur with Mingarelli et al. [18] and classify this system as a candidate.

3.4. Ultra-light Companion or Planet

This type of pulsar is defined by a companion mass criterion of less than 0.08 $M_{\odot }$, with our sample containing only PSR J1908+2105.

PSR J1908+2105: This system represents a transitional MSP system exhibiting characteristics of both black widow and redback spider pulsars. With a short orbital period of 3.51 hours and a companion mass of ∼0.055 $M_{\odot }$ [19], it aligns with black widow properties, yet displays anomalous high-frequency eclipsing behavior typical of redbacks - eclipses persist for >30% of the orbital period and extend up to 4 GHz, the highest frequency band of the Parkes Ultra-Wideband receiver. This makes it one of only three MSPs (alongside PSR J1723−2837 and J1731−1847) observed to eclipse at frequencies >3 GHz [20]. We tentatively list the Gaia source 4519819661567533696 as a candidate.

3.5. Spiders

PSR J0838−2827: This pulsar belongs to the redback system, which has an orbital period of 0.21 days. This is consistent with the short orbital periods characteristic of redback systems. Its optical companion is a low-mass M-type dwarf star that orbits the primary star every 5.15 hours. The mass of the primary star is comparable to that of a neutron star [21]. Initially, the pulsar’s orbital period was unknown due to time constraints in the observations [22], but its redback nature was later confirmed through improved photometric methods. We identify Gaia source 5645504747023158400 as the pulsar’s companion, with an angular separation of 0.002 arcseconds. Thongmeearkom et al. [23] employed Monte Carlo techniques to estimate the pulsar’s mass and found a companion mass of 0.13 $M_{\odot }$ for J0838−2827, which is slightly below the typical redback range of 0.2–0.4 $M_{\odot }$. However, the mass ratio q = 7.06 ($M_{psr}/M_c$) aligns with the high mass ratio characteristic of redback systems.

PSR J0955−3947: We identified its companion in Gaia 5419965878188457984, with an angular separation of 0.003 arcseconds. In a study by Drake et al. [24], an optical variable with a period of 9.3 hours was considered a potential redback binary system, while Li et al. [25] provided further empirical support through the analysis of X-ray variability and spectroscopic features. Thongmeearkom et al. [23] determined the mass of the companion of J0955−3947 to be 0.30 $M_{\odot }$, which is within the typical range for a redback system, and found that the mass ratio q = 4.43 was consistent with that of a neutron star-led binary system. Through two MeerKAT observations, Thongmeearkom et al. [23] reported that J0955−3947 exhibited several brief eclipsing events (lasting ∼40 seconds) near orbital phase 0.1, caused by companion material occulting the pulsar signal — a hallmark feature of redback systems.

PSR J1036−4353: Clark et al. [26] identified PSR J1036−4353 as a redback system with a 6.23 hour orbital period and a minimum companion mass of 0.23 $M_{\odot }$. During their observations, the researchers detected eclipsing events in which the pulsar became undetectable for several orbital phases when positioned behind the companion. This indicates that the eclipses are periodic and caused by companion material, which is a hallmark of redback systems. PSR J1036−4353 exhibits similarities in eclipsing behavior, orbital parameters and optical characteristics to PSR J1803−6707 (described below). We identified its companion as Gaia 5367876720979404288, with an angular separation of 0.001 arcseconds, and it is classified as a main-sequence star.

PSR J1628−3205: We identified its companion as Gaia 6025344817107454464, with an angular separation of 0.125 arcseconds. Li et al. [27] obtained a minimum companion mass of 0.16 $M_{\odot }$ and an orbital period of 5 hours for PSR J1628−3205 through three observations carried out in 2013 and 2014. Additionally, light-curve modelling revealed that the companion’s effective temperature ranges from 3560 K to 4670 K, exhibiting a lower density and a higher temperature than main-sequence stars of a similar mass (possibly due to pulsar irradiation or an evolutionary influence), which are characteristic features of redback companions.

PSR J1803−6707: Clark et al. [26] established an orbital timing solution for PSR J1803−6707. This revealed a ∼9 hour orbital period and a minimum companion mass of 0.26 $M_{\odot }$. However, significant orbital period variations were also observed, requiring multiple higher-order orbital frequency derivative terms to fully describe the orbital phase evolution across the 1-year timing dataset. During two dedicated observations in 2020, pulsed emission was detected for approximately 30 minutes only once, which is characteristic of the broad eclipses typical of redback systems. Clark et al. [26] confirmed the optical counterpart using ultra-high-frequency MeerKAT observations, demonstrating a 1.5-magnitude light curve with a single-peaked structure caused by pulsar irradiation heating (similar to that observed in PSR J2215+5135 by Schroeder [28]). Notably, the optical counterpart remains detectable across all orbital phases — a defining feature of redback companions. We identified the companion as Gaia 6436867623955512064, with an angular separation of $\theta$ = 0.002 arcseconds.

PSR J1910−5320: PSR J1910−5320 was discovered during observations of unassociated Fermi-LAT gamma-ray sources with the MeerKAT telescope [29]. Dodge et al. [30] obtained optical light curves using the ULTRACAM high-speed multi-band photometer combined with MeerKAT timing and SOAR/Goodman spectroscopy. This confirmed its binary nature through sinusoidal radial velocity curves with an orbital period of 0.35 days. The detection of radio pulsations as part of the Fermi-LAT source survey [26] validated the prediction of a redback system made by Au et al [29].

PSR J2055+1545: PSR J2055+1545 was discovered by Lewis et al. [31] through radio observations. It exhibits a spin period of 2.16 milliseconds and an orbital period of 4.8 hours. Radio signals from the pulsar are eclipsed by companion material for 36% of the orbital period, and the companion is estimated to have a mass of 0.29 $M_{\odot }$. Using the Gran Telescopio Canarias (GTC), Turchetta et al. [32] obtained optical spectra, revealing that the companion is consistent with a G-dwarf star. Temperature variations across orbital phases — from 5749 ± 34 K at inferior conjunction to 6106 ± 35 K at superior conjunction — support the hypothesis that the companion is irradiated by the pulsar wind, a characteristic feature of redback systems. We identified the companion as Gaia 1763537692275731328, with an angular separation of $\theta$ = 0.02 arcseconds.

PSR J2333−5526: Similar to PSR J1910−5320, PSR J2333−5526 was identified as a redback candidate through searches for optical and X-ray counterparts to unassociated Fermi-LAT gamma-ray sources [33]. We identified its companion in Gaia 6496325574947304448. Thongmeearkom et al. [23] confirmed an orbital period of 6.9 hours and a companion mass of 0.38 $M_{\odot }$ through MeerKAT radio searches. They noted that the companion resides in a nearly edge-on orbit [34]. Clark et al. [35] validated this scenario by detecting a gamma-ray eclipse. For PSR J2333−5526, the eclipse duration spans 0.045–0.07 orbital phases [23].

4. Discussions and Conclusions

PSR J0337+1715 is a highly hierarchical triple star system with nearly coplanar, nearly circular orbits. It consists of a millisecond pulsar and white dwarfs in the inner and outer orbits [36,37]. It originated from a close triple main-sequence star system with initial masses of approximately 10.0, 1.10, and 1.30 $M_{\odot }$. Following common envelope evolution, the inner and outer orbits contracted significantly (with the inner orbit shortening from several hundred days to approximately 2.5 days and the outer orbit from several thousand days to approximately 17 days) [38]. During the shell helium burning phase, the core of the primary star underwent mass transfer and subsequently exploded as a supernova, forming a neutron star. Following this, the system underwent two LMXB phases. During the outer LMXB phase, the tertiary star (initially 1.30 $M_{\odot }$) transferred mass to the neutron star, forming a 0.410 $M_{\odot }$ white dwarf and substantially widening the outer orbit. During the inner LMXB phase, the secondary star (initially 1.10 $M_{\odot }$) transferred mass to the neutron star, resulting in the formation of a 0.197 $M_{\odot }$ white dwarf and the recycling of the neutron star into a millisecond pulsar. For this pulsar, the parallax and proper motion measurements obtained by Gaia are more precise than those obtained from ATNF. We may potentially employ Very Long Baseline Interferometry (VLBI) observations to further refine the astrometric parameters of this star.

PSR J1024−0719 was initially classified as an isolated millisecond pulsar with a 5.2 ms spin period and period derivative of 1.8×${10}^{-20}$ [39]. Kaplan et al. [40] later identified it as a long-period binary (2–20 kyr orbital period) with a low-mass (0.4 $M_{\odot }$), low-metallicity (−0.9 dex) main-sequence companion. The extreme orbital length results in a significant second period derivative due to acceleration variations from orbital motion [40]. The leading hypothesis is dynamical exchange in a globular cluster: an original compact binary (possibly containing a white dwarf) interacted with a third star, ejecting the white dwarf and capturing the companion into a wide orbit [41]. Alternatively, the system may have formed from a hierarchical triple in the galactic disk that underwent supernova-driven acceleration [42]. However, the companion’s low metallicity—consistent with globular cluster stars—supports the former scenario. High stellar densities in globular clusters promote such dynamical interactions, enabling the formation of rare binaries like this one. For such long-period pulsars, Gaia measurements are less precise than those obtained through pulsar timing. When using the astrometric parameters of these pulsars, it may be worth considering using the timing results instead of Gaia’s parameters.

In Section 2, we used the apparent magnitude and distance of the pulsar to make a rough estimate of its absolute magnitude. However, this approach excludes some well-known pulsar binary systems due to errors in distance calculations between the pulsar and Earth. Pulsar distance measurements currently commonly employ the Dispersion Measure (DM) method, which is based on the frequency-dependent dispersion effect that causes delays to electromagnetic waves from the pulsar as they propagate through the interstellar medium. This can be expressed mathematically as:

$$\begin{array}{c}\hfill \Delta t={\displaystyle \frac{DM}{2.41\times {10}^{-4}}}\left({\displaystyle \frac{1}{v_1^2}}-{\displaystyle \frac{1}{v_2^2}}\right)\end{array}$$

(4)

$$\begin{array}{c}\hfill DM={\int }_0^d{n}_e\left(l\right)dl\end{array}$$

(5)

Where $\Delta t$ represents the time delay between frequencies ${\nu }_1$ and ${\nu }_2$ (where ${\nu }_1$ < ${\nu }_2$), and $n_e\left(l\right)$ denotes the electron density along the path. Using Equation (5), we can estimate the distance as $d\approx {\displaystyle \frac{DM}{\left\{n_e\right\}}}$ if we determine the average electron density $n_e$ along the line-of-sight path from Earth to the pulsar. However, the distribution of free electrons in the interstellar medium is inhomogeneous, clumpy and complex, so calculations must rely on electron density models. The YMW16 model [43] and the NE2001 model [44] are currently the most widely used. Nevertheless, estimation accuracy depends entirely on the accuracy and suitability of these models, as they can never perfectly describe local complex structures. For example, if the line of sight passes through an unmodelled, small, high-density cloud, significant errors will be incurred in the distance estimation, depending on the specific sightline direction, with errors typically being stochastic. Furthermore, the real interstellar medium contains turbulence and small-scale structures that these models inevitably smooth out. Therefore, it is not accurate to use the distance d to calculate the absolute magnitude.

Based on the above analysis, we precisely calculate the pulsar’s absolute magnitude using Equation (6):

$$\begin{array}{c}\hfill H_g=m_g+5{log}_{10}{\displaystyle \frac{\mu }{\mathrm{mas}{\mathrm{yr}}^{-1}}}-10+A_g\\ \hfill H_g\equiv M_g+5{log}_{10}{\displaystyle \frac{v_{tr}}{4.74057\mathrm{km}{\mathrm{s}}^{-1}}}\end{array}$$

(6)

Here, $m_g$ and $A_g$ represent the apparent magnitude and extinction of the pulsar in the Gaia data, respectively. The proper motion magnitude is given by $\mu =\sqrt{{\mu }_{\alpha }^2+{\mu }_{\delta }^2}$, and ${\nu }_{tr}$ denotes the transverse velocity of the source. Figure 4 presents a Hertzsprung-Russell (HR) diagram with $G_{BP}-G_{RP}$ on the x-axis and absolute magnitude $M_G$ on the y-axis. $G_{BP}$ and $G_{RP}$ correspond to the magnitudes in the blue and red photometric bands, respectively, and their difference serves as a color index. This index increases for redder (cooler) objects and decreases for bluer (hotter) ones. The dashed line represents the empirical cutoff defined by Equation (3). As shown in Figure 4, the majority of the stars in our identified sample lie along the main sequence band that extends from the top left to the bottom right. This distribution is consistent with the theoretical HR diagram.

For the majority of pulsars among all the samples we included, their proper motion differences are generally around 1 mas/yr. However, for some pulsars, such as PSR J2150+3427, this difference can be several tens of mas/yr. This discrepancy arises because Gaia’s astrometric solutions are computed under the assumption that all observed stars are single systems. Consequently, the proper motions derived for certain pulsars through Gaia astrometry differ significantly from those obtained via timing methods [45]. When selecting our pulsar samples, we initially excluded 4 pulsars: PSR J0210+5845, J1529−5609, J1903+0327 and J1937+1658. PSR J1903+0327 and PSR J1937+1658 were excluded due to their slightly larger angular distances ($\theta$ = 2.25 and 2.84 arcsec, respectively), while the remaining two were excluded due to the limited information currently available on them. Figure 5 briefly presents the magnitude and proper motion information for the 3 candidate pulsars and the 4 initially excluded pulsars. We have identified certain samples among the pulsars listed in Section 3 that can serve as targets for subsequent research endeavors. For example, we could use the Australian Long Baseline Array (LBA) network to conduct further observations of PSR J1803−6707 and PSR J1910−5320. PSR J1910+1256 is a pulsar included in the International Pulsar Timing Array (IPTA) DR2. Observations using the Very Long Baseline Array (VLBA) would make a contribution to the work of Guo et al. [46] on linking pulsar reference frames. We find that the flux densities of the pulsars listed in Section 3 are all notably low. PSR J2108+4516, for example, has an average flux density of merely 3 mJy at 600 MHz, which suggests that high-precision timing observations using the Five-hundred-meter Aperture Spherical radio Telescope (FAST) would be necessary to detect it successfully. We will conduct follow-up observations and continuous monitoring of these pulsars, as well as the samples presented in Figure 5, to further confirm their relevant information.

We conducted companion searches for pulsars in the ATNF pulsar catalogue using high-precision astrometric data from Gaia DR3. Our method identified 58 binary systems, including 9 newly confirmed pulsars. The ESA Gaia mission’s fourth data release (DR4) is scheduled for completion in 2026. It will contain 66 months of observational data covering approximately two billion stars. This release will include additional galaxies and quasars and will significantly improve data precision and depth. In the future, the Chinese VLBI Network (CVN), comprising the FAST, the Shanghai 65-Metre Radio Telescope, the Nanshan 26-Metre Radio Telescope, the Miyun 50-Metre Radio Telescope, and so on, will provide robust support for pulsar position and motion measurements thanks to its high-precision positioning capabilities [47]. Integrating the CVN’s radio observation capabilities with the optical observation data from Gaia DR4 will enable multi-wavelength collaborative observations to be conducted in order to study the radiation characteristics and variability patterns of pulsars and their binary systems across different wavebands.