2.1 MWA

The MWA (Tingay et al. Reference Tingay2013; Wayth et al. Reference Wayth2018) is a low-frequency (80–300 MHz) precursor to the SKAO located at the Inyarrimanha Ilgari Bundara, the CSIRO’s Murchison Radio-astronomy Observatory, in Western Australia. One project carried out with the MWA is the Galactic Plane Monitor (GPM). These observations were taken between June–September 2022 and June–September 2024 under project code G0080 (see Methods of Hurley-Walker et al. Reference Hurley-Walker2023). The prime objective of the campaign was to search for long-period transients by scanning the Galactic plane over $285^\circ \lt l \lt 65^\circ$ , $|b|\lt\pm15^\circ$ on a bi-weekly cadence at 185–215 MHz; a full description will be released by Hurley-Walker et al. (in preparation). Observations were undertaken in five-minute ‘snapshots’ and mosaics were formed for each night of observing. For any given region, the effective integration time was 30 min for the data taken in 2022, and 45 min in 2024. The four sources discussed in this paper were discovered in these data as clearly varying between ‘on’ and ‘off’ states; the survey and selection process details will be presented in a separate publication (Hurley-Walker et al. in preparation).

2.2 uGMRT

The uGMRT (Swarup Reference Swarup, Cornwell and Perley1991; Gupta et al. Reference Gupta2017) comprises of 30 antennas with a 45-m diameter, located in Pune, India. Our observations were taken between October 2023 and December 2024, using band-4 (550–750 MHz). We aimed to schedule the observations just before or after the candidate pulsars reach superior conjunction relative to their companions (i.e. as the candidate pulsar enters or exits eclipse), so that we can study the eclipsing material. Each observation commenced with a $\sim$ 5 min scan of a bright, compact radio source (e.g. 3C286) for flux and bandpass calibration. The targets were observed for $\sim$ 30 min each, bookended with a $\sim$ 5 min nearby phase calibrator scan. Calibration and flagging were done using Common Astronomy Software Applications (CASA; CASA Team et al. Reference Team2022). Imaging was done using WSClean (Offringa et al. Reference Offringa2014), performing multi-frequency synthesis imaging, with a Briggs $-$ 0.5 weighting scheme (Briggs Reference Briggs1995) and joint-channel deconvolution. The bandwidth was split between 2–10 channels, depending on the image’s signal-to-noise ratio (S/N). Aegean Footnote ^a (Hancock et al. 2012; Hancock et al. Reference Hancock, Trott and Hurley-Walker2018) was used to extract the flux densities and associated errors. When the source is eclipsed, the flux density is reported as an upper limit by measuring the pixel value at the location of the source.

Pulsar search mode data were also recorded simultaneously with the interferometric visibility data, using the phased-array observing mode of uGMRT. These data utilised 4 096 channels sampling 200 MHz bandwidth centred at 650 MHz, and a sampling time of 81.96 $\unicode{x03BC}$ s. RFIClean (Maan, van Leeuwen, & Vohl Reference Maan, van Leeuwen and Vohl2021) was used for the first level of RFI excision, as well as conversion of the data to SIGPROC filterbank format. A periodicity search for pulsed signals was then conducted using PRESTO (Ransom Reference Ransom2011) on the data recorded in October and November 2023; see Section 3.2 for more details.

2.3 ASKAP

The Australian Square Kilometre Array Pathfinder (ASKAP) (Hotan et al. Reference Hotan2021) Variables and Slow Transients Survey (VAST) (Murphy et al. Reference Murphy2021), is a radio survey designed to detect transient sources in the image domain. The Galactic components of VAST observes 42 fields over $|b|\lt6^\circ$ and $|l|\lt10^\circ$ totalling 1 260 deg $^2$ . It operates at a central frequency 888 MHz, taking observations roughly every two weeks since November 2022. The integration time for each observation is $\sim$ 12 min, with a typical sensitivity of 0.24 mJybeam $^{-1}$ . Science-ready data products are made available through CSIRO ASKAP Science Data Archive (CASDA) (Huynh et al. Reference Huynh, Dempsey, Whiting, Ophel, Ballester, Ibsen, Solar and Shortridge2020) and are produced using the ASKAPsoft pipeline (Guzman et al. 2019).

2.4 Parkes Murriyang

Murriyang, the Parkes 64-m radio telescope, is a 64-m single-dish radio telescope situated in New South Wales. Our observations utilised the MEDUSA back-end and were taken during available Director’s Time (under project code PX095) in the Ultra Wide Bandwidth (UWL) receiver (Hobbs et al. Reference Hobbs2020), which offers a continuous frequency range of 704–4 032 MHz. The observations were recorded at 64- $\unicode{x03BC}$ s time resolution, with 128 channels per subband and 26 subbands across the UWL band. Each source was observed for approximately 15 min. The data was processed using standard processing software, PRESTO performing a periodicity and acceleration search (that is described in detail in Section 3.2).

2.5 ATCA

The Australia Telescope Compact Array (ATCA; Wilson et al. Reference Wilson2011), consists of six 22-m antennas located in Narrabri, New South Wales. We took observations using the 16-cm receiver, utilising the Compact Array Broadband Backend (CABB) in the L-band (1 100–3 100 MHz). The observations were recorded with 1 MHz channels and a correlator integration time of 10 s, under the project code CX544. Each observation observed PKS B1934-638 for $\sim$ 10 min for flux density and bandpass calibration. A $\sim$ 2-min scan of a nearby phase calibrator was taken between each $\sim$ 15-min source scan. The data was flagged and calibrated using standard procedures with CASA and imaged using the task tclean, selecting a Briggs weighting and robust parameter $-$ 0.5.

2.6 MeerKAT

MeerKAT S-band observations were undertaken under proposal code SCI-20241101-NH-01, using the S4 window at 2 625–3 500 MHz to minimise the effects of scattering. At the time of observing, this sub-band had by far the cleanest RFI environment,Footnote ^b ensuring that less data is lost to flagging. The observations included typical bandpass, polarisation, and phase calibrators, as well as phase-up and test pulsar observations. Each target source was observed for two 30-min integrations; we attempted to schedule the observations such that the sources would be observed as close to the candidate pulsars’ inferior conjunctions as possible (i.e. least likely to be in eclipse). This was done to increase the likelihood of detecting pulsations. As well as correlator observations undertaken at 8-s/854.492-kHz resolution, we also employed the Pulsar Timing User Supplied Equipment (PTUSE; Bailes et al. Reference Bailes2020) in search mode, using 37.45- $\unicode{x03BC}$ s sampling.

The correlator data were calibrated using the standard SARAO SDP calibration pipeline and imaged using WSClean. We imaged each 30-min scan separately to yield good signal-to-noise for source position measurements, while also yielding some information about the time-variability of the sources at S-band. A periodicity search for pulsation signal was conducted on the PTUSE data using the PRESTO software (see Section 3.2 for details).

Figure 1. Panels (i),(ii),(iii): DECaPS z-band images of the fields centred on the radio positions of GPM J1723 $-$ 33, GPM J1734-28, and GPM J1752 $-$ 30. The seeing in each image is about 1 arcsec. The radio positions are marked with a red ellipse corresponding to their positional uncertainties, except in the case of GPM J1815 $-$ 14 where the positional uncertainty is smaller than one pixel in the image. Optical sources in the DECaPS2 band-merged catalogue within 1.5 arcsec of the radio positions are marked by blue circles to guide the eye. The letter appearing next to each optical position corresponds to the catalogue entry in Table 3. Panel (iv): Pan-STARRS1 z band deep stacked image of the field centred on GPM J1815 $-$ 14. The radio position is marked by a red cross, as its positional uncertainty is the size of about one pixel in this image. Each cutout is $10 \times 10$ arcsec. One faint potential optical counterpart is detected for GPM J1723 $-$ 33, multiple are detected for GPM J1734 $-$ 28 and GPM J1752 $-$ 30, while none are detected for GPM J1815 $-$ 14.

2.7 Optical and near-infrared

We searched for optical and near-IR counterparts to the four radio positions among existing survey data. The deepest optical imaging of the fields around the first three sources is available through the Dark Energy Camera Plane Survey (DECaPS2; Saydjari et al. Reference Saydjari2023), which used 30 s exposures with the 4.0-m Blanco telescope in Chile in 2016 and 2019. As shown in the coadded fields in Figure 1, there is one potential counterpart consistent with the radio position of GPM J1723 $-$ 33. We identified multiple potential optical counterparts from the DECaPS imaging that are compatible with the radio positions of GPM J1734-28 and GPM J1752 $-$ 30. We report the DECaPS2 band-merged catalogue measurements for these optical sources in Section 3.4. However, due to significant interstellar extinction in the direction of the sources ( $A_z \sim 1.7$ and $2.3$ mag respectively for a distance of 8 kpc, Zucker et al. Reference Zucker2025), it is possible that the true optical counterparts to the radio sources are obscured.

We checked for potential near-IR counterparts among the VISTA Variables in the Via Lactea (VVV) DR5 catalogue and imaging (Minniti et al. Reference Minniti2010), which was performed between 2010–2015 and targeted central parts of the Galaxy primarily in the $K_s$ band with the 4.1-m VISTA telescope in Chile. No catalogue sources are detected within 1.5 arcsec of these positions. Although extinction is less of an obstacle at longer wavelengths, the exposure time used by the VVV survey imaging is quite short ( $\sim4$ sec). Visual inspection of the fields revealed only marginal evidence for flux at the locations of the optical sources in the DeCAPS imaging.

The field of GPM J1815 $-$ 14 is not covered by the DECaPS2 footprint. To search for nearby optical sources, we inspected deep stack images from Pan-STARRS1, which used a 1.8-m telescope in Hawai’i (Chambers et al. Reference Chambers2016; Flewelling et al. Reference Flewelling2020). In the fourth panel of Figure 1, we show a z band stack of 30 and 60 s exposures with a total exposure time of 780 s. The radio position of GPM J1815 $-$ 14 is close to a bright $z \sim 14.3$ mag foreground star. We note that there is one detection in the Pan-STARRS1 catalogue near this radio position in a single epoch, in the g band with $g_{PSF} = 21.7 \pm 0.2$ (were $g_{PSF}$ denotes the Point Spread Function magnitude) offset from the radio position by 0.8 arcsec observed in June 2012. The quality factor for the detection is 0.84, indicating that the detection may be suspect. No other optical sources compatible with the position are detected. The extinction in this direction is extremely high, and so the optical counterpart is likely entirely obscured by dust ( $A_z \sim 4.8$ mag at 8 kpc, Green Reference Green2018; Green et al. Reference Green, Schlafly, Zucker, Speagle and Finkbeiner2019). In the near-IR, the field of GPM J1815 $-$ 14 was covered by the VVV extended survey in 2016–2019 (VVVX; Saito et al. Reference Saito2024), and no sources are detected in the catalogue or upon inspection of the imaging down to a mean sensitivity of $\sim18.5$ mag in the $K_s$ band.

3.1 Fermi gamma-ray associations

As pulsars are known to emit $\gamma$ -ray emission, we cross-match our sources with the Fermi Fourth Full Catalogue of LAT Sources (4FGL-DR4; Abdollahi et al. Reference Abdollahi2022). The nearest Fermi-4FGL sources found, and corresponding angular separations for our sources, are provided in Table 4 and Figure 2.

Figure 2. MeerKAT S-band images for each source, overlaid with the Fermi 95 $\%$ error ellipse (shown in yellow) of the closest unassociated $\gamma$ -ray source. The angular separation between the radio and $\gamma$ -ray sources is indicated. See Section 3.1 for more details.

Possible counterparts for GPM J1723 $-33$ and GPM J1815 $-14$ lie within the Fermi 95 $\%$ error ellipse. The nearest 4FGL sources to GPM J1734 $-28$ and GPM J1752 $-30$ lie outside their respective ellipses, at 5 $\sigma$ and 3 $\sigma$ , respectively.

While this may indicate a low probability of association for the last two sources, Abdollahi et al. (Reference Abdollahi2022) notes that $\gamma$ -ray source localisation in the Galactic Plane is challenging due to high background emission. To assess the probability that GPM J1734 $-28$ and GPM J1752 $-30$ are associated with their nearest 4FGL sources even though they lie outside the 95 $\%$ confidence ellipse, we calculate the fractional sky area covered by Fermi error ellipses at 5 $\sigma$ and 3 $\sigma$ , respectively. We consider a 100 deg $^2$ region around our source. The percentage of sky occupied by the ellipses for GPM J1734 $-28$ and GPM J1752 $-30$ is found to be 13 $\%$ and 5 $\%$ , respectively. The results indicate that the likelihood of our source falling that close to a Fermi ellipse by chance is relatively low. Therefore, we still consider these as potential counterparts.

The implications of these results are discussed in detail in Section 4.7.

3.2 Pulsation searches

Pulsation searches were carried out with three radio telescopes: uGMRT (550–750 MHz), Parkes (704–4 032 MHz) and MeerKAT (2 625–3 500 MHz), covering a broad frequency range. A periodicity search was performed on the PTUSE data using the PRESTO software suite on the Setonix cluster hosted by Pawsey. Following radio-frequency interference (RFI) mitigation, the data were dedispersed over trial dispersion measures (DMs) from 0–5 000 pc cm $^{-3}$ . Each DM trial was searched for periodic signals using Fourier-based acceleration searches with harmonic summing up to 8 harmonics and a maximum acceleration parameter of $z_{\text{max}} = 300$ . Harmonic summing up to 8 was chosen for the acceleration search, as it is generally sufficient for detecting MSPs, particularly in systems such as spiders, which are known to be MSPs and typically do not require higher harmonic summing due to their narrow pulse profiles.

The acceleration parameter z denotes the maximum number of Fourier frequency bins a pulsar signal can drift due to a constant line-of-sight acceleration during the observation. It is given by:

(1) \begin{align}z = \left( \frac{Gm_c^3}{m_{\text{tot}}^2} \right)^{1/3} \left( \frac{2\pi}{P_{{\text{orb}}}} \right)^{4/3}\end{align}

where $m_c$ is the mass of the companion, $m_{\text{tot}}$ is the total mass of the system and $P_{\text{orb}}$ is the orbital period of the system. The maximum detectable constant line-of-sight acceleration is given by

(2) \begin{equation} a=\frac{z_{\text{max}}\ cP}{T_{\text{obs}}}\end{equation}

where c is the speed of light, P is the pulsar spin period, and $T_{\text{obs}}$ is the observation duration.

We arbitrarily choose a sufficiently high $z_{\text{max}}$ value of 300 for this search, as the expected z for spider pulsars is much smaller. For a example, a RB system with $m_c = 0.2\ \text{M}_\odot$ and $m_p = 1.4\ \text{M}_\odot$ , and a 12 hr orbital period, has $z\approx5.7$ , corresponding to an acceleration of 16 ms $^{-2}$ .

The PTUSE data were then folded to the period and DM of the candidates found in the search, and the diagnostic plots produced were manually inspected.

We also divided the Parkes UWL band into three segments (704–1 344, 1 344–2 368, 2 368–4 032 MHz) and searched each independently. High-frequency searches help reduce the potential for pulse smearing due to scattering, as the scattering timescale scales with frequency as $\tau_s \propto \nu^{-4}$ , while low-frequency searches are particularly effective for detecting pulsars due to their steep negative spectral indices. All candidates were manually inspected in the uGMRT and Parkes data, and no significant pulsar signals were identified.

An 8 $\sigma$ detection of GPM J1723 $-33$ was made using MeerKAT (see Figure 3). The spin period is 2.11030426(16) ms and the DM is 469(2) pc cm $^{-3}$ . The acceleration search measured $z=11$ , corresponding to a line-of-sight acceleration of 2.1493840(2) ms $^{-2}$ (for $T_{\text{obs}}=1\,800\,\mathrm{s}$ ).

Figure 3. Pulsation detection of the MSP GPM J1723 $-$ 33 with MeerKAT radio telescope. left: The main panel shows the evolution of pulsations over time, folded on its spin period of 2.11030426(16) ms. The top subpanel displays the integrated pulse profile, frequency-scrunched to enhance gaussian significance, while the side panel shows the reduced $\chi^2$ , indicating the significance of the detection. right: The upper plot shows pulse phase as a function of frequency, confirming the broadband nature of the signal. The bottom plot displays the dispersion measure (DM) against $\chi^2$ , highlighting the optimal DM solution. This figure is derived from the analysis described in Section 3.2.

All candidates from the MeerKAT observations were manually inspected; however, no convincing pulsar signals were detected. Potential reasons for the non-detection of radio pulses are discussed in Section 5.2.

3.3 Radio imaging

Despite the non-detection of pulsations, we can still obtain valuable insight into these systems through image-domain observations. All four sources were detected as variable sources in the image-domain by multiple radio surveys (MWA, VAST, uGMRT, ATCA, and MeerKAT). The variability observed in the time-integrated snapshot surveys such as VAST and MWA is typically characterised by a bimodal switch between ‘on’ and ‘off’ (see Figure 4 for an example).

Figure 4. left: VAST detection of GPM J1815 $-14$ on 2022-12-21T03:20:25.5 UTC. right: VAST non-detection of GPM J1815 $-14$ on 2022-11-14T06:57:53.7 UTC.

Figure 5. left: Light curve of GPM J1815 $-14$ using the MWA GPM data at 200 MHz. right: The corresponding folded light curve on the orbital period, $P_\mathrm{orbit}= $ 9.81969(2) h. The plots demonstrate the increasing eclipse duration over time; see Section 4.1 for details.

The ASKAP positions of the sources were used for initial follow-up observations, determined using weighted averages with an accuracy of $\sim1$ arcsec. We refined the positions using MeerKAT observations, which offer improved angular resolution due to their longer baselines and higher sensitivity. We measured the positions using Aegean, achieving sub-arcsecond precision. The final positions and associated uncertainties are listed in Table 1. We adopted the MeerKAT positions for our analysis due to their improved accuracy and resolution.

The MeerKAT 30-min integration images had RMS noise levels between 8–9 $\unicode{x03BC}$ Jy beam $^{-1}$ . GPM J1723 $-33$ was detected at a flux density of 0.11(2) mJy in one scan, with a non-detection in the other. GPM J1734 $-28$ was detected with a flux density of 60(9) $\unicode{x03BC}$ Jy in both scans. GPM J1752 $-30$ was not detected in the first scan but was detected with a flux density of 40(8) $\unicode{x03BC}$ Jy in the second. GPM J1815 $-14$ was measured with a flux density is 0.14(1) mJy in both scans.

The uGMRT observations were scheduled to coincide with the time of superior conjunction and successfully captured data during this phase for GPM J1734 $-28$ , GPM J1752 $-30$ , and GPM J1815 $-14$ . The wide-bandwidth of uGMRT enables the data to be segmented by frequency, allowing an investigation of frequency-dependent effects.

We aimed to observe the sources with ATCA at higher radio frequencies and at orbital phases where the candidate pulsars are in inferior conjunction to their companion. GPM J1723 $-33$ was detected at a flux density of 0.35(6) mJy and GPM J1752 $-30$ at 0.48(6) mJy. However, the S/N ratios ( $\sim$ 6 and $\sim$ 8, respectively) were too low to allow for time or frequency segmentation of the data for orbital period analysis. ATCA did not detect GPM J1734 $-28$ and GPM J1815 $-14$ due to noise levels, where the background root-mean-square (RMS) is 30 and 70 $\unicode{x03BC}$ Jy beam $^{-1}$ , respectively.

3.4 Search for optical and near-IR counterparts

Optical observations of RBs are used to characterise the companion and the binary orbital configuration through detailed light curve modelling and radial velocity studies. We searched for optical and near-IR counterparts to the four new RB candidates in survey imaging and identified potential counterparts to three shown in Figure 1.

Only one faint source in the band-merged DECaPS2 catalogue is detected within 1.5 arcsec of GPM J1723-33. For GPM J1734-28, there are four faint optical sources potentially consistent with the radio position, two of which are visually detected in Figure 1. For GPM J1752 $-$ 30, there are three faint optical sources consistent with the radio position, two of which are marginally visible upon inspection of the DECaPS images. We summarise the DECaPS2 catalogue measurements of these optical sources in Table 3. No near-IR sources within 1.5 arcsec of the radio positions are detected in the VVV DR5 catalogue. Visual inspection of the stacked and epoch imaging showed only marginal evidence for excess brightness at the positions of the DECaPS2 sources, and we found no evidence for variability.

Ground-based follow-up of these fields to determine the correct optical counterparts is challenging due to the high extinction and intrinsic faintness of the targets. Optical observations with larger telescopes or deeper observations in the near-IR bands are needed to detect the sources and variability signatures, which would solidify any association with the radio sources and allow detailed modelling of the binary configurations. With the initial orbital periods determined in Section 4.2, observations could be strategically planned to sample orbital modulation.

The search for optical counterparts to GPM J1815 $-$ 14 did not yield any optical sources consistent with its radio position. This field is strongly affected by dust extinction, and it is unlikely that the optical counterpart will ever be accessible. Future deep observations in the near-IR may reveal a counterpart.

4.1 Dynamic evolution

Figure 5 shows the MWA GPM data for GPM J1815 $-14$ , revealing that the system exhibits significant variability over time. Notably, the source is not detected for most of the 2024 data. RB systems are known to evolve dynamically, as the companion star evaporates due to intense radiation from the pulsar. As the density of the ablated material varies over month–timescales, the amount of flux density that is absorbed or scattered may change, resulting in variations in the observed eclipse duration. This likely explains the decreasing detectability of GPM J1815 $-14$ over time. Therefore, systems like this one may experience varying detectability over time, making consistent detection challenging.

We also considered refractive interstellar scintillation (RISS) as a possible explanation for the observed variability. Hancock et al. (Reference Hancock, Charlton, Macquart and Hurley-Walker2019) created an all-sky model for RISS, which can predict variability on a variety of timescales, survey locations, and observing frequencies. We use the available code RISS19 Footnote ^c to obtain an order-of-magnitude estimate for the RISS at the location of GPM J1815 $-14$ . We find the scintillation strength is very high near the galactic plane, but the expected flux density modulation at 200 MHz is only a few per cent, with a characteristic timescale of hundreds of years. This suggests that RISS is unlikely to produce significant short-term variability in our source. Intrinsic pulsar variability is another potential cause, but the observed flux density appears correlated with the orbital phase, suggesting an eclipse-related origin rather than random pulsar behaviour.

4.2 Orbital period analysis

We searched for periodic variability in the MWA and VAST datasets independently for each source, as spider systems often exhibit frequency-dependent eclipses. Using the MJD, flux density, and corresponding flux density errors, we applied the Lafler-Kinman’s String Length (LKSL) method (Clarke Reference Clarke2002), using the P4J Footnote ^d package. LKSL is well-suited for identifying periodicities in sparse, unevenly sampled light curves. It works by folding the data on a range of trial periods and calculating the string length at each trial; the minima in string length correspond to the best trial period. The P4J implementation inverts this measure so that peaks in the periodogram (LKSL power versus period) correspond to best candidate periods.

We searched trial periods between 0.5 and 24 h, corresponding to a frequency range of approximately 1 to 0.04167 cycles/day. The frequency grid was sampled at a resolution of $10^{-4}$ , and the 12 most significant local peaks were selected. We independently searched for periodicity in the MWA and VAST data and then compared candidate periods to identify consistencies.

Figure 6. The normalised light curves, where the flux densities for the MWA (red points, 200 MHz) and VAST (blue points, 888 MHz) data are folded on the best-fit orbital periods of the systems. The flux densities for each dataset have been normalised to their maxima to remove the effect of the steep spectral indices of the source.

Using this subset of best candidate periods, we phase-folded the flux density measurements from both the VAST and MWA datasets, incorporating the MeerKAT detection and non-detection times, as well as the uGMRT observation window, to assess whether these observations aligned consistently in phase space. We determined the best candidate orbital period as the one most compatible with all the data.

Spurious peaks in the periodogram can arise from measurement errors, noise, or random fluctuations in the data. To assess the significance of our candidate period, we applied a bootstrap-based approach. We generated 200 surrogate (bootstrap) light curves by randomly resampling the original data, then computed a periodogram for each, and extracted the 20 highest maxima. An extreme value probability density function (PDF) was then fitted to the distribution of these maxima. From the fitted PDF, we derived 90 and 99 $\%$ confidence thresholds, which are shown in the periodograms in Figure A1. The candidate period peaks for GPM J1734 $-28$ , GPM J1752 $-30$ , and GPM J1815 $-14$ exceed the 90 $\%$ threshold, allowing us to confidently conclude that the periodicities are genuine.

Table 1. The table presents the positional coordinates of the sources discussed in this paper (see Section 3.3), along with the 200 MHz flux density at the pulsar’s inferior conjunction, obtained from model fitting of Equation (3). It also includes the orbital period $P_\mathrm{orb}$ , and the reference MJD $T_{0}$ (see Section 4.2).

Table 2. Summary of the follow-up observations outlined in Section 2.

Table 3. Optical sources in the DECaPS2 band-merged catalogue that are compatible with the radio positions of GPM J1723-33, GPM J1734-28, and GPM J1752 $-$ 30. The Label column indicates the marker used for the optical sources in Figure 1. Unique object identifiers, offsets from the radio positions, and mean grizY background-corrected magnitudes from DECaPS2 are given. All catalogue entries here are $\gt5\sigma$ detections.

Table 4. Details of the potential Fermi sources associated with each source in this paper. The table lists the separation between the radio and Fermi 4FGL sources, along with the semi-major and semi-minor axes of the 95 $\%$ error ellipses for the 4FGL sources supplied by Fermi.

Figure 6 presents the folded light curves on the best candidate orbital period found using the LKSL method. All systems are found to have short orbital periods ( $\lt$ 1 day) and exhibit eclipsing behaviour. The eclipsing variability is not clearly visible in the VAST GPM J1723 $-33$ data due to large errors caused by deconvolution issues. The final orbital periods and associated uncertainties, derived from the full width at half maximum of the periodogram peak, are given in Table 1.

Figure 7. The folded light curves for each source observed with the MWA (200 MHz), uGMRT (550–750 MHz), and VAST (888 MHz), fitted with either Equation (3) or (4) at each frequency. The flux densities are normalised to the maximum flux density ( $S_0$ ) derived from the fits. For clarity, a vertical offset of 0.5 has been added between successive frequencies relative to the MWA baseline in panels (c) and (d).

4.3 Eclipse properties

The properties of the eclipse can be quantified by fitting a double Fermi-Dirac function (see Zic et al. Reference Zic2024):

(3) \begin{equation}S(\phi) = S_0 \left[ \left( \frac{1}{1 + e^{\frac{\phi - \phi_i}{w_i}}} - \frac{1}{1 + e^{\frac{\phi - \phi_e}{w_e}}} \right) \right]\end{equation}

where $\phi$ is the phase, $\phi_i$ and $\phi_e$ are the phases at the eclipse ingress and egress, $w_i$ and $w_e$ are the widths of the ingress and egress, and $S_0$ is the flux density at the pulsar’s inferior conjunction (i.e. out of eclipse).

For the uGMRT data, and for the VAST data of GPM J1734 $-28$ and GPM J1752 $-30$ , which sparsely sampled the orbital phase, we fitted a modified version of Equation (3) to model only the ingress of the eclipse:

(4) \begin{equation}S(\phi) = S_0 \left[1 + \left( \frac{1}{1 + e^{\frac{\phi - \phi_i}{w_i}}} \right) \right]\end{equation}

Equations (3) and (4) were fitted to the folded light curves from the MWA, VAST, and uGMRT observations (see Figure 7) using non-linear least-squares fitting via scipy.curve_fit, with bounds and initial guesses chosen to ensure numerical stability. The flux densities are normalised by the maximum flux density obtained from the fit, denoted $S_0$ in Equation (3) and (4). The folded light curves are scaled to a reference MJD, $T_{0}$ (given in Table 1), corresponding to the time when the candidate pulsars are at superior conjunction relative to their companions. The uncertainty in $T_{0}$ arises from the propagation of errors from the ingress and egress phase measurements, derived from the covariance matrix of the fitted models.

4.4 Radio spectra

Figure 8 shows the spectra of the sources when they are at inferior conjunction with their companions. For MWA, VAST, and uGMRT, the flux density at inferior conjunction is taken as $S_{0}$ , which is derived from the fitted functions (see Section 4.3). The values for $S_{0}$ at 200 MHz are provided in Table 1, with the uncertainty on each value assumed to be 20 $\%$ , accounting for typical calibration and flux scale errors. The data points for ATCA and MeerKAT are taken from the measured flux density at the sources’ inferior conjunction.

Figure 8. Spectra of the four sources fitted with Equations (5)–(7), using flux densities corresponding to the inferior conjunction of the candidate pulsars. The resulting model parameters are listed in Table 5.

The spectra were fitted with the following models, taken from Swainston et al. (Reference Swainston, Lee, McSweeney and Bhat2022):

Simple power law

(5) \begin{equation}S_\nu=S_0\left( \frac{\nu}{\nu_0} \right)^\alpha\end{equation}

where $S_0$ is the reference flux density at frequency $\nu_0$ , and $\alpha$ is the spectral index.

Log parabola

(6) \begin{equation}S_\nu = S_0\space\left( \frac{\nu}{\nu_0} \right)^{a+b\log \left(\frac{\nu}{\nu_0}\right)}\end{equation}

where a is the spectral index and b is the curvature parameter.

Broken power Law

(7) \begin{equation}S_\nu =\begin{cases}S_0\space \left( \frac{\nu}{\nu_0} \right)^{a_1} & \text{for} \ \nu \leq \nu_b \\S_0 \space \left( \frac{\nu}{\nu_0} \right)^{a_2} \left( \frac{\nu_b}{\nu_0} \right)^{a_1-a_2} & \text{for} \ \nu \gt \nu_b\end{cases}\end{equation}

where $a_1$ and $a_2$ are the spectral indices below and above the break frequency $\nu_b$ .

The spectral models were fitted using curve-fit, with uncertainties on the parameters derived from the covariance matrix of the fit. The resulting fits are shown in Figure 8 and the corresponding parameter values are given in Tabel 5. The best model was determined using Akaike information criterion (AIC) (Jankowski et al. Reference Jankowski, van Straten, Keane, Bailes, Barr, Johnston and Kerr2018), which estimates the information lost for a given model. The best models are highlighted in the table, and all sources are found to have steep negative spectral indices.

4.5 Companion modelling

We compute the range of companion mass ( $m_c$ ) and semi-major axis (a) for a binary system using the mass function

(8) \begin{equation}f(m_p,m_c) = \frac{4\pi^2}{G}\frac{(a \sin i)^2}{P_\mathrm{orbit}^2} = \frac{(m_csini)^3}{(m_p+m_c)^2}\end{equation}

where $m_p$ is the fixed mass of the neutron star assumed to be $1.4\ M_{sun}$ (which is typically seen for pulsars Zhang et al. Reference Zhang2011), the inclination angle i is fixed to $\pi/2$ (edge-on orbit), and G is the gravitational constant. From Section 3.2, we calculated an acceleration for GPM J1723 $-33$ of 2.1493840(2) ms $^{-2}$ . Using Equation (3) from Ng et al. (Reference Ng2015), we calculated a mass function of 5.7 $\times 10^{-6}$ and a corresponding minimum companion mass of 0.023 M $_\odot$ . Assuming the same acceleration value for the other three sources, we derived the following minimum companion masses GPM J1815 $-14$ : 0.05 M $_\odot$ , GPM J1734 $-28$ : 0.15 M $_\odot$ , and GPM J1752 $-30$ : 0.063 M $_\odot$ . These values are used as lower limits for the minimum companion mass in Equation (8).

Figure 9 plots the minimum companion mass against semi-major axis for known binary systems in the ATNF Pulsar Catalogue Footnote ^e (Catalogue Version 2.6.0). The range of $m_c$ values for our four systems is plotted, showing a distribution along the line consistent with the observed population. Assuming an upper limit of 0.4 M $_{\odot}$ for the companion mass, as typically observed in RB systems (Roberts Reference Roberts and van Leeuwen2013), we can estimate an upper limit for the semi-major axis of our systems. The upper limit for semi-major axis values calculated are as follows, GPM J1723 $-33$ : $0.410\,{\rm R}_{\odot}$ , GPM J1734 $-28$ : $1.01\,{\rm R}_{\odot}$ , GPM J1752 $-30$ : $0.67\,{\rm R}_{\odot}$ , and GPM J1815 $-14$ : $0.63\,{\rm R}_{\odot}$ .

Using these upper limits on semi-major axis, the Roche lobe radius can be computed using the Eggleton approximation (Eggleton Reference Eggleton1983). The approximate values computed are as follows, GPM J1723 $-33$ : $0.11\,{\rm R}_{\odot}$ , GPM J1734 $-28$ : $0.28\,{\rm R}_{\odot}$ , GPM J1752 $-30$ : $0.19\,{\rm R}_{\odot}$ , and GPM J1815 $-14$ : $0.17\,{\rm R}_{\odot}$ .

Figure 9. Minimum companion masses against the projected semi-major axes for known binary pulsars from the ATNF catalogue, along with the range of companion masses allowed for each of our targets, calculated using Equation (8). The colour map indicates the orbital period in hours. The bottom panel provides a zoomed-in view of the top panel. The upper limits on the minimum companion mass for our sources are marked with an ’X’ along each line. Binary MSPs are categorised by companion type: Neutron Star (NS), Helium White Dwarf (He), Carbon-Oxygen White Dwarf (CO), and Ultra-light companion (UL).

As the systems are no longer accreting, we can reasonably assume that the companions are not overflowing their Roche lobes. Thus, these radii represent lower limits on the sizes of the companion stars. Assuming a lower-limit effective temperature of 3 000 K, we use the Stefan–Boltzmann law to estimate the corresponding bolometric luminosities, and from there calculate lower-limit apparent magnitudes (assuming distances determined in Section 5.2). The calculated approximate lower limit on apparent magnitude are as follows, GPM J1723 $-33$ : 22 , GPM J1734 $-28$ : 17, GPM J1752 $-30$ : 19, and GPM J1815 $-14$ : 19. Comparing these estimates with the DECaPS2 photometry in Table 3, we find that our lower-limit apparent magnitude estimates for GPM J1723 $-33$ and GPM J1752 $-30$ are consistent with sources listed in the table, specifically, those with obj_id 4522440858633389827 and 4954927160361131656, respectively.

4.6 Absorption models

RB pulsars are seen to ablate their companion stars. This ablated material surrounds the binary system and can absorb or scatter the pulsar’s radio emission, often leading to long-duration eclipses in the observed radio emission. We explore the eclipse mechanisms suggested by Thompson et al. (Reference Thompson, Blandford, Evans and Phinney1994) that could be responsible for the observed eclipses. This is done using the flux density from the fitted Fermi-Dirac models (see Section 4.3) for MWA, VAST, and uGMRT. The optical depth is then computed using Equations (9)–(12), as outlined by Thompson et al. (Reference Thompson, Blandford, Evans and Phinney1994).

4.6.1 Cyclotron absorption

One possible mechanism is cyclotron absorption, where the companion star and pulsar wind are in a magnetic field. Following the same approach as Thompson et al. (Reference Thompson, Blandford, Evans and Phinney1994), we equate the magnetic pressure $P_{B}=\frac{B_{E}^2}{8\pi}$ of the plasma to the pulsar wind pressure $U_{E}=\frac{\dot{E}}{4\pi ca^2}$ . Taking the spin-down power of the pulsar $\dot{E}\sim10^{34}$ (Kansabanik et al. Reference Kansabanik, Bhattacharyya, Roy and Stappers2021), we can calculate the characteristic magnetic field $B_E$ , the values are given in (10). The optical depth for cyclotron absorption is given by (Thompson et al. Reference Thompson, Blandford, Evans and Phinney1994):

(9) \begin{equation}\tau_{cyc}(\nu)= \frac{\pi}{2} \frac{m^{m+1}}{m!} \left(\frac{mk_{B}T}{2m_{e}c^2}\right)^{(m-1)} \frac{n_{e}e^2L_{B}}{m_{e}c\nu}\end{equation}

where $m=\frac{\nu}{\nu_{B}}$ is the cyclotron harmonic, where $\nu_{B}=\frac{eB}{2\pi m_{e}c}$ is the cyclotron frequency. T is the temperature in Kelvin, and $L_B=|\frac{ds}{dlnB}$ | is the scale length of the magnetic field in cm. We assume that the scale length of electron density variations and magnetic field variations are similar for simplicity. Therefore, we take $n_{e}L_{B}$ to be equal to the average electron column density $N_{e}$ (Thompson et al. Reference Thompson, Blandford, Evans and Phinney1994). For Cyclotron absorption to be valid, the temperature must satisfy $T\lt \frac{m_ec^2}{2k_{B}m^3}$ .

To calculate the optical depth for our systems, we allow $N_e$ to vary between $10^{4}$ and $10^{30} {\rm cm}^{-2}$ , T between $10^{2}$ and $10^{9}$ K, and B between $10^{-2}$ and $10^{2}$ G, following the approach of Zic et al. (Reference Zic2024). The minimum temperature found to produce eclipses at 888 MHz is $\gt10^8$ K, for each system.

4.6.2 Compton scattering

For Compton scattering, the optical depth is given as (Thompson et al. Reference Thompson, Blandford, Evans and Phinney1994):

(10) \begin{equation}\tau_{ind}(\nu)= 4 \times 10^{6} \frac{N_{\text{e}} S_{\nu}}{\nu^{2}}\langle\,f(\phi) \rangle \, |\alpha + 1| \cdot \left(\frac{d_{kpc}}{a}\right)^2 M \end{equation}

where $S_{\nu}$ is the mean flux density of the pulsar in mJy at frequency $\nu$ MHz, $\alpha$ is the spectral index of the incident radiation, $d_{kpc}$ is the distance to the scattering centre from the observer in kiloparsecs, and a is the distance between the companion and the pulsar in cm. $\langle\,f(\phi) \rangle$ is the angular factor which is averaged out over the scattering region and M is the demagnefication factor (Kumari et al. Reference Kumari2024) which can be between 0 and 1. The distance of GPM J1723 $-$ 33 is $\sim$ 9.7 kpc obtained using the YMW16 model, and we have assumed a distance of 8 kpc for the other three sources (see Section 3.2 for details). The maximum $\tau_{ind}$ found for each system is much less than 1.

4.6.3 Free-Free absorption

For free-free absorption, the optical depth is given by (Thompson et al. Reference Thompson, Blandford, Evans and Phinney1994):

(11) \begin{equation}\tau_{ff}(\nu) \approx 3.1 \times 10^{-8} \cdot \frac{f_{\text{cl}} N_e^2}{T^{3/2} \nu^{2} L} \cdot \ln \left( 5 \times 10^{10} \cdot \frac{T^{3/2}}{\nu} \right)\end{equation}

where $f_{cl}$ is the clumping factor, T is the temperature, and L is the absorption length. We fix the absorption length as $L=2R_{E}$ , where the eclipse radius is defined as $R_E=\pi a\Delta\phi$ . Where, a is the semi-major axis calculated in Section 4.5, and $\Delta\phi=\phi_{e}-\phi_{i}$ is the eclipse duration orbital phase.

To constrain the physical parameters for free-free absorption in Equation (11), we performed a Markov Chain Monte Carlo (MCMC) analysis. The details of the MCMC setup are provided in Section 4.6.5.

4.6.4 Synchrotron absorption

Synchrotron absorption arises from a population of non-thermal electrons, often assumed to possess a power-law energy distribution $n(E)=n_{0}E^{-p}$ . The optical depth is given by (Thompson et al. Reference Thompson, Blandford, Evans and Phinney1994):

(12) \begin{equation}\tau_{\text{sync}(\nu)} = \left(\frac{ 3^\frac{(p+1)}{2} \Gamma \left( \frac{3p + 2}{12} \right) \, \Gamma \left( \frac{3p + 22}{12} \right)}{4}\right) \left( \frac{\sin \theta}{m} \right)^{\frac{p+2}{2}} \frac{n_0 e^2}{m_e c \, \nu } L\end{equation}

Where L is the absorption length, p is the power-law index, $\theta$ is the angle of the magnetic field lines with our line of sight, $n_0$ is the non-thermal electron density, $m_e$ is the mass of the electron, m is the cyclotron harmonic at frequency $\nu$ , e is the charge on the electron, and c is the speed of light. We fix the absorption length as $L=2R_{E}$ (as in Equation 11), and a fixed viewing angle of $\pi/3$ . An MCMC analysis is performed to constrain the physical parameters in Equation (12), with further details presented in Section 4.6.5.

4.6.5 MCMC

We perform an MCMC analysis to constrain the physical parameters for free-free and synchrotron absorption, using the emcee sampler. We modelled the observed eclipse durations as a function of frequency using the models given in Equations (11) and (12). Our MCMC setup and parameter constraints follow those used in Zic et al. (Reference Zic2024).

Sampling was performed in logarithmic space, using uniform priors over the following ranges: for free-free absorption, $\log(f_{cl}) \in (0, 10)$ , $\log(T) \in (0, 9)$ , and $\log(N_{e}) \in (10, 20)$ ; for synchrotron absorption, $\log(B) \in (\!-\!2, 2)$ , $\log(p) \in (0, 8)$ , and $\log(n_{0}) \in (0, 8)$ .

The log-likelihood function compares how well the modelled optical depths, $\tau_{\mathrm{ff}(\nu)}$ and $\tau_{\mathrm{sync}(\nu)}$ , fits the observed optical depth measurements $\tau(\nu)$ , assuming gaussian errors $\tau_{\text{err}(\nu)}$ .

We calculated $\tau(\nu)$ at frequencies corresponding to MWA (200 MHz), uGMRT (500–600 MHz), and VAST (888 MHz) at the phase just after the eclipse. This was done using the radiative transfer equation $S=S_{0}(\nu)e^{-\tau(\nu)}$ , where $S_{0}(\nu)$ is the flux density at the pulsars’ inferior conjunctions (from Equation 3).

The log-likelihood is then given by

(13) \begin{equation}\log L = -\frac{1}{2} \sum \left( \frac{ \tau(\nu) - (\tau_{\mathrm{ff}}(\nu)/\tau_{\mathrm{sync}}(\nu)) }{ \tau_{\mathrm{err}}(\nu) } \right)^2\end{equation}

and the total log-probability is the sum of the log-prior and the log-likelihood. Figures 12 and 13 present the marginal posterior distributions for the free-free and synchrotron absorption parameters, the results are discussed in Section 5.3.

4.7 Fermi gamma-ray properties

Figure 10. left: Radio flux density at 1 400 MHz ( $S_{1\,400}$ ) plotted against Fermi $\gamma$ -ray energy flux ( $E_{100}$ ) at 100 MeV. There is no evident correlation between radio and $\gamma$ -ray flux. All four of our systems lie within the observed range of both $E_{100}$ and $S_{1\,400}$ values for the current population of binary pulsars. right: $\gamma$ -ray variability index against spectral curvature significance. The plots display the four sources from this paper alongside ATNF binary pulsars with Fermi associations (see Section 4.7).

Section 3.1 revealed possible Fermi $\gamma$ -ray associations for each source. The two primary requirements for $\gamma$ -ray emission from pulsar candidates are: they have a spectral curvature significance higher than $3\sigma$ and the variability index is lower than 27.7 (Abdollahi et al. Reference Abdollahi2022). Pulsars are expected to be steady $\gamma$ -ray emitters, hence, their variability index should be low. The spectral curvature significance quantifies how much the $\gamma$ -ray spectrum deviates from a simple power law, with pulsars expected to have curved spectra. Figure 10 plots the $\gamma$ -ray spectral curvature significance against the $\gamma$ -ray variability index (taken from the 4FGL catalogue), for $\gamma$ -ray binary pulsars taken from the ATNF pulsar catalogue (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005) along with the four sources in this paper. All four of our systems align with the characteristics of the currently known binary pulsar population.

We compiled $\gamma$ -ray fluxes at 100 MeV (E $_{100}$ ) and radio flux densities at 1 400 MHz (S $_{1\,400}$ ) for known $\gamma$ -ray binary pulsars in the ATNF catalogue. The E $_{100}$ values are taken from the 4FGL catalogue (Abdollahi et al. Reference Abdollahi2022) and S $_{1\,400}$ from the ATNF pulsar catalogue (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005). In Figure 10 we show S $_{1\,400}$ as a function of E $_{100}$ for all $\gamma$ -ray binary pulsars, together with our four systems. There is no evident correlation between radio and $\gamma$ -ray flux. All four of our systems lie within the observed range of both E $_{100}$ and S $_{1\,400}$ values for the current population of binary pulsars.