© The Authors 2023

1. Introduction

Double neutron star (DNS) binaries are the evolutionary endpoint of massive binary systems (M_stars > 8 M_⊙) that survive two supernovae while remaining bound (e.g. Tauris et al. 2017; Vigna-Gómez et al. 2018; Chattopadhyay et al. 2020). In this picture, the primary star undergoes a Type Ib/c core collapse supernova, becoming the first-born neutron star (NS). Subsequently, the system undergoes a common-envelope phase as the secondary evolves out of the main sequence (e.g. van den Heuvel 2019), after which the outer shells of the evolved secondary are expelled, turning the system into a circular, compact NS–He star binary (e.g. Chattopadhyay et al. 2020). As the secondary continues to evolve, a case BB Roche-lobe overflow (RLO) ensues in which the primary NS sustains partial recycling (Tauris et al. 2015), burying its magnetic field and spinning up to a rotation period of a few tens of milliseconds (Bhattacharya & van den Heuvel 1991). Finally, the system must avoid the disruption during the second, ultra-stripped supernova in which the secondary NS is born with the expulsion of 0.1 to 1 M_⊙ from the system, and with a maximum kick velocity of ≈100 m s⁻¹ (Tauris et al. 2015). The nature of this supernova is still a matter of discussion as it can be triggered by either electron capture in the degenerate ONeMg core or iron core collapse, depending on how much stripping has been suffered by the He star, with more massive remnants undergoing the latter channel and suffering greater mass losses and supernova kicks (Tauris et al. 2015). The distinction between these two formation mechanisms is blurry due to the randomness of the resulting orbital parameters arising from the unpredictability of the supernova kick direction and magnitude, but the remaining mass of the second-born NS may be a tell-tale sign, as it has been seen that higher eccentricity systems correlate with more massive second-born NSs and less rapidly spinning recycled pulsars (e.g. Sengar et al. 2022; Andrews & Mandel 2019; Faulkner et al. 2005; Tauris et al. 2017). Nonetheless, if all of these steps are completed without disruption, then a new, eccentric DNS system is born.

Observable radio pulsars discovered in DNS systems are amongst the most useful astrophysical tools. With dedicated follow-up and observing campaigns using radio telescopes, their orbital parameters and component masses can be measured to high precision through the technique of pulsar timing (e.g. Lorimer & Kramer 2005), allowing for tests of formation channels and even of fundamental physics. This is demonstrated in the literature on formation mechanisms (e.g. Tauris et al. 2017), breakthrough tests of gravity (e.g. Taylor et al. 1979; Taylor & Weisberg 1982; Kramer et al. 2021), and the validation or exclusion of dense matter models (e.g. Özel & Freire 2016; Hu et al. 2020). Beyond radio observations, the discovery of pulsars in DNS systems has also had an indirect impact on other fields of astrophysics. The observation of the orbital decay of B1913+16 (Taylor et al. 1979) experimentally demonstrated the existence of gravitational waves as predicted by general relativity (GR) for the first time and supported the prediction of its merger in 301 Myr, providing scientific justification for the construction of ground-based gravitational-wave observatories. Indeed, with knowledge of the pulsar DNS population in the Milky Way, their rates of orbital decay, and of the sensitivity of blind surveys on the sky, predictions of the observed cosmic rate of NS mergers can be made (e.g. Kim et al. 2003, 2010, 2015; Pol et al. 2019, 2020; Grunthal et al. 2021). The most recent estimate (Grunthal et al. 2021) provides an upper limit of ℛ_local ≤ 597 yr⁻¹, which is consistent with the observed rates at ground-based gravitational wave detectors (Abbott et al. 2021), highlighting the synergy between radio and gravitational wave observations. Finally, the first detection of a gravitational wave signal emitted by the coalescence of two NSs, GW170817 (Abbott et al. 2017a), opened a new era for multi-wavelength and multi-messenger astronomy, highlighting DNS systems among other astrophysical systems even further (Abbott et al. 2017b).

However, observable pulsars in DNS systems are are a rarity. At the time of writing only 17 Galactic DNS systems are confirmed, of which only nine are competitive for tests of gravity and NS merger rate predictions (Table 1). Indeed, simulations predict a Galactic DNS formation rate of just 5−31 Myr⁻¹, with only ∼0.13% of massive binaries surviving the DNS formation channels (Vigna-Gómez et al. 2018). The number of discoveries is also limited by observational reasons. With the exceptions of J0737−3039B, J1906+0746 and perhaps J1755−2550 (Table 1), typically only the first-born, recycled NS is observable due to their longer-lasting emission and larger beaming fractions, while the second-born NS is set to spin down and cross its pulsar death line a few Myr after birth (Lorimer & Kramer 2005). Limitations are also computational, as the extreme orbital motion of pulsars in compact, relativistic DNS systems hampers the effectiveness of traditional periodicity searches from pulsar surveys due to the Doppler smearing of the signal across an observation (Johnston & Kulkarni 1991). It is, therefore, necessary to implement computationally expensive and sophisticated pulsar searches to increase the sample of known DNS systems (and hence maximise their science output). These generally include acceleration searches (e.g. Keith et al. 2010), jerk searches (e.g. Andersen & Ransom 2018; Eatough et al. 2021; Suresh et al. 2022) or even template bank searches (e.g. Allen et al. 2013; Balakrishnan et al. 2022).

The MPIfR-MeerKAT Galactic Plane survey at L-band (MMGPS-L; Kramer et al. 2016; Stappers & Kramer 2016; Padmanabh et al. 2023) has been highly successful in this objective. Designed to discover massive, compact binaries in the southern Galactic plane with a time resampling-based acceleration search, it has yielded the discovery of two previously unknown DNS systems, PSR J1155−6529 (Berezina et al., in prep.) and PSR J1208−5936 (J1208−5936 from now on). Additionally, the MMGPS-L is the most sensitive pulsar survey in the southern sky, making its large coverage area highly relevant to estimate the rate of NS mergers in the Milky Way. In this work, we present the follow-up study of J1208−5936, a 28.71-ms recycled pulsar in a highly eccentric orbit with another NS, with whom it is bound to merge in 7.2 ± 0.2 Gyr due to the emission of gravitational waves. In addition, we also use this discovery and the improvement in depth of sky coverage provided by the MMGPS-L to recompute the predicted NS merger rate.

This paper is structured as follows: in Sect. 2, we present the discovery and the derivation of an early orbital solution. In Sect. 3, we show a phase-coherent timing solution and perform preliminary mass measurements with almost one year of timing data. In Sect. 4, we study its general properties, including its pulse profile, a comparison with other pulsars in DNS systems and a discussion of its formation channel. In Sect. 5, we search for pulsations from its companion. In Sect. 6, we investigate the implications of the discovery and the performance of the MMGPS-L for estimations of the local NS merger rate, and the detection of these events by ground-based gravitational-wave observatories. And finally, in Sect. 7 we discuss possible biases and prospects for future improvements.

2. Discovery and orbital solution

Led by the Max Planck Institute for Radioastronomy¹ (MPIfR) in collaboration with the South African Radio Astronomy Observatory² (SARAO), the MMGPS-L was a fast Fourier transform (FFT)-based pulsar search survey with the MeerKAT radio interferometer³ (Jonas & MeerKAT Team 2016). Covering 900 sq. deg. on the southern Galactic plane with the primary aim to discover faint relativistic binaries, the MMGPS-L implemented a time-domain acceleration search (Padmanabh et al. 2023) with the PEASOUP⁴ pipeline (Barr 2020).

J1208−5936 was discovered on 30 May 2021 in an observation from May 6th as a topocentric 28.706-ms signal in the FFT with signal-to-noise ratio S/N = 15.1, a line-of-sight acceleration of −8.86 m s⁻², DM = 344.2 pc cm⁻³ and a PulsarX⁵ (Men et al. 2023) fold of S/N ≈ 19. Upon discovery, we performed a first 20 min-long follow-up observation at the same position on 4 June 2021. This yielded a re-detection with S/N ≈ 15, and a PRESTO/prepfold⁶ fold of S/N ≈ 12 with a changed period of P = 28.697 ms, thus confirming both its existence and its binary nature.

The interferometric nature of MeerKAT, and therefore the multibeam pattern of the follow-up observation allowed for precise localisation of J1208−5936 with the SeeKAT⁷ software (Bezuidenhout et al., in prep.), which compares the S/N of detections in different beams and does a maximum-likelihood estimate of the best position of the source, taking into account beam positions and their point-spread functions as derived by the Mosaic⁸ software (Chen et al. 2021a). This improved the sky position to an uncertainty of just ∼5 arcsec, much better than the ∼20 arcsec precision provided by the beam size alone, and which was only 3.3 arcsec away from the timing-derived position (Table 2).

From 4 June 2021 to 31 August 2021, we typically scheduled dedicated sessions twice a week, consisting of two 20-min observations stored as filterbank data in the accelerated pulsar search user supplied equipment cluster in South Africa (APSUSE-mode; Padmanabh et al. 2023). Each session had the two observations conducted just a few hours apart. We folded each observation into archives with the dspsr⁹, where cycles of the pulsar are stacked every 8 or 10 seconds. We then cleaned these archives with the clfd¹⁰ radio-frequency interference (RFI) excision software (Morello et al. 2019). These resulting data allow for accurate tracking of the evolution or drift of the pulse over time, leading to measurements of the barycentric period (P_bary) at each observing epoch (t_obs) with the PSRCHIVE/pdmp¹¹ software (Hotan et al. 2004), but also with the timing software TEMPO2¹² (Hobbs et al. 2006; Edwards et al. 2006) for observations with large line-of-sight acceleration or jerk.

We implemented a modified version of the roughness estimate algorithm (REA) from Bhattacharyya & Nityananda (2008) to solve for the orbital period, available in estimateOrbit.py¹³. Given adequate orbital coverage, the REA is more sensitive to orbital periods in systems with non-zero eccentricity than the Lomb–Scargle periodogram. It performs a search in the orbital period P_b space by folding the P_bary(t_obs) series and evaluating the smoothness of the resulting curve with

(1)

where t_i = t_obs,i/P_b, trial − mod(t_obs, i, P_b, trial) are the folded observation epochs and ϕ corresponds to the orbital phase (mean anomaly). We then select P_b, trial values that minimise R, as they produce a smooth fold and are likely to correspond to the true P_b value. The normalisation by the difference in orbital phase Δϕ is a modification of the original REA presented in Bhattacharyya & Nityananda (2008) that gives higher significance to data points that are closer in orbital phase, dealing with possible gaps in the orbital coverage.

The REA found a significant signal at P_b, trial = 0.632 days. We then used pyfitorbit¹⁴ software to fit for the remaining parameters using the best REA P_b, trial as a first guess, which resulted in a confirmation of the orbital period and a fit for the remaining Keplerian parameters (Fig. 1), with a eccentricity e = 0.34 and a mass function f_M = 0.208 M_⊙.

Fig. 1. Orbital solution of J1208−5936 based on the Pbary(tobs) time series between 5 May 2021 and 31 August 2021. Observations with precise acceleration and/or jerk measurements were used to extract several data points (5 May: −10.5 m s−2; 11 June: +22.4 m s−2, and 16 August: +7.2 m s−2, −3.6 × 10−3 m s−3; at ϕ ≈ 0.72, 0.95 and 0.15). Assuming a canonical pulsar mass of Mp = 1.35 M⊙, the minimum companion mass computed from the mass function is Mc ≳ 1.1 M⊙ (see Table 2).

3. Timing and mass measurements

After August 2021, follow-up observations were scheduled to focus on the timing analysis. The observing time was dropped to monthly 5-min-long observations, but in exchange for the reduced observational time, observations were moved to being recorded with the dedicated MeerKAT pulsar timing user supplied equipment (PTUSE; Bailes et al. 2020) and APSUSE simultaneously. The advantages of PTUSE over APSUSE are GPU-based coherent de-dispersion, recording of full-Stokes information, real-time data folding and a finer sampling resolution of 9 μs at the same radio band, thus enabling high precision pulsar timing.

To achieve a phase-coherent timing solution, we used dracula2.py¹⁵, an implementation of the dracula algorithm (Freire & Ridolfi 2018) with python and TEMPO2. dracula searches for timing solutions assuming different combinations of phase wraps in between the observations, finding the unique solution that accounts for all the rotations of the pulsar from start to end. For this, we used PTUSE pulsar archives from the newer observations and APSUSE archives produced from older search data with dspsr, both of them with 128 phase bins across the profile. Those archives were de-dispersed and scrunched in frequency with the PSRCHIVE/pam command, and three Times of Arrival (ToAs) were produced per observation with PSRCHIVE/paas, using a single, narrow von Mises function fitted with PSRCHIVE/pat as an initial timing template. dracula.py quickly converged into a single solution without the need of fitting for the sky position in the process, indicative of the quality of the SeeKAT multibeam localisation.

Given the massive, eccentric and compact nature of the system, we performed a semi-coherent orbital campaign with 14 observations at selected orbital phases from 2 March 2022 to 7 March 2022, accumulating a total of 11 h. Using our phase-connected timing model as a predictor, the orbital phases of the observations were chosen to cover features of the unabsorbed signal of a potential Shapiro delay signal (Freire & Wex 2010) and the passage of periastron.

Together with an additional observation from 8 April 2022, the added PTUSE folded pulse profile with 1024 bins presents S/N ≈ 200, while the added 512-bin APSUSE profile presents S/N ≈ 170, the main difference in S/N coming from the improvement provided by coherent de-dispersion of PTUSE data. These two high S/N profiles were scrunched into four frequency sub-bands for frequency-resolved timing of the entire data set in order to account for DM evolution. Frequency-resolved analytical timing templates were produced with paas both for the PTUSE and APSUSE data sets. Subsequently and with dspsr, all PTUSE data were folded into 1024-bin, four-channel archives, and APSUSE data (when PTUSE data were not available) were folded into 512-bin, four-channel archives, with sub-integration times of at most 10 min in length in both cases. These files had been excised of RFI a priori at full time and frequency resolution using first the clfd¹⁶ software (Morello et al. 2019) and then manually with PSRCHIVE/paz upon inspection of the remaining RFI. Finally, pat was used with the Fourier Phase Gradient algorithm (Taylor 1992) to generate 384 frequency-resolved ToAs.

Assuming GR to be the correct theory of gravity, we analysed the produced ToAs to constrain the pulsar mass (M_p), the companion mass (M_c) and the inclination angle of the orbit (i) through two different methods. The first one is by measuring the theory-independent post-Keplerian (PK) parameters that arise from relativistic corrections of the orbital motion, and then deriving probability densities for M_p, M_c and i on which we quote median values and 1σ uncertainties from the two adjacent 34.1% percentiles (68.2% credible intervals). In this work, we measure the advance of periastron passage (), the amplitude of the Einstein delay (γ_E), and the orthometric parameters of Shapiro delay: the third-order orthometric amplitude (h₃) and the orthometric ratio (ς) as defined in Freire & Wex (2010). In the first post-Newtonian approximation of GR, arises from the rotation of the Keplerian ellipse of the orbit in the direction of the orbital motion, depending on the total mass M_t = M_p + M_c as

(2)

while h₃ and ς parameterise the unabsorbed part of the Shapiro delay (Freire & Wex 2010). These describe the delay of the pulses propagating through the companion’s gravitational field, depending on M_c and i as

(3)

and

(4)

for cos i > 0, which is applicable to our study as we do not have any constrain on the angle of the ascending node (Ω_A, see Sect. 7.2 for a discussion on this). It is worth noting that we use the h₃ and ς parametrisation introduced in Freire & Wex (2010) instead of the classic range (r) and shape (s) parameters used in other works (e.g. Kramer et al. 2021) because they better describe systems with low i and are less correlated with each other. And finally, the Einstein delay is caused by the periodic modulation of the relativistic time dilation due to the pulsar’s changing orbital velocity and its motion across the companion’s gravitational field, its amplitude γ_E being modelled by

(5)

These parameters are implemented in DDH, which is an extended version of the Damour-Deruelle pulsar timing model (DD; Damour & Deruelle 1986) that implements the orthometric parameters h₃ and ς instead of r and s. The fit was done with TEMPO2/TempoNEST¹⁷, a multi-nested Bayesian sampling plug-in to TEMPO2 (Lentati et al. 2014), in order to find stable values for both h₃ and ς with realistic uncertainties.

The other method is to assume GR from the start with the DDGR model (Taylor & Weisberg 1989), which measures M_t and M_c directly to model the PK effects. For this measurement we implemented a common likelihood approach (first introduced by Splaver et al. 2002, see therein for more details) with chi2Map.py¹⁸, computing the likely constraints from the quality of the TEMPO2 fits of the DDGR model in a uniformly spaced grid on the M_t − cos i plane, which produces agnostic prior that follows the random distribution of i values of binary systems in the sky.

The resulting probability distributions, derived from multiplying the measured PK parameters from DDH and the χ² values from DDGR, were then marginalised into one-dimensional probability densities for M_p, M_c and cos i, on which we quote median values and 1σ uncertainties from the two adjacent 34.1% percentiles (68.2% credible intervals). The resulting TempoNEST fit of the DDH model and the chi2Map.py exploration of the DDGR χ² space are consistent with each other. The most significant PK parameter is , detected in DDH with ≈900σ significance, yielding a highly precise measurement of the total system mass at

(6)

This is in excellent agreement with the

(7)

given by the direct DDGR TEMPO2 fit. We note that ς and γ_E are detected with low significance, and there is an important upper limit of Shapiro delay amplitude h₃. The derived 1σ constraints from each DDH PK parameter and the χ² DDGR fits are overlayed in Fig. 2, with the DDGR contour tracing very strictly the limits imposed by , and being consistent with the loose limits imposed by γ_E, h₃ and ς. The most likely M_p, M_c and i values derived from the DDH and DDGR models are in very good consistency with each other, with

(8)

Fig. 2. Corner plot showing the constraints on the inclination angle and component masses from PK parameters and the χ2 mapping of DDGR solutions. Central plots: mass–mass diagrams portraying the DDH PK parameters constraints, with each color corresponding to a different parameter as indicated in the legend (solid lines: nominal value, dashed lines: 1σ limits), and the 1σ limits from DDGR χ2 mapping in black. The shaded grey area on the right plot represents the region excluded by the mass function (i > 90), while the shaded areas on the left plot represent the areas outside the prior for Mp (outside of 0.6 < Mp < 1.5 M⊙). The explored regions of the Mp and Mc were decided based on the limits given by the mass function and Mt (from ). Outer plots: marginalised one-dimensional probability densities for Mp, Mc and cos i from DDGR χ2 mapping, showcasing the median value (black solid line) and the 31.4%, 47.4% and 49.9% percentiles on both sides (shaded areas). The resulting mass constraints are consistent with a pair of NSs.

from multiplying probability distributions given by the PK parameter limits, and

(9)

from the χ² mapping of DDGR. Additionally, as a final consistency check, we also re-derived the likely values of γ_E, h₃ and ς by calculating their marginal one-dimensional probability densities from the DDGR χ² mapping experiment, on which we quote median values and 1σ uncertainties from the two adjacent 34.1% percentiles, and attest that they are in good consistency with their DDH limits (Table 2).

We also attempted to fit any orbital period derivative Ṗ_b or apparent change of the projected semi-major axis ẋ, but none of them yield significant limits with DDH. From the DDGR χ² mapping limits on M_p, M_c and i and assuming GR, we predict them to have values of s s⁻¹ and ls s⁻¹, but the current timing sensitivity is not high enough to yield a proper measurement. Nonetheless, Ṗ_b is not likely to be detected without the addition of another decade of timing and it may be contaminated by effects introduced by the galactic acceleration field and proper motion in the sky, while ẋ is likely to be dominated by proper motion effects (Sect. 7.1, Table 5).

4. Properties of J1208–5936

4.1. General characteristics

Table 2 presents the timing parameters resulting from fits to the ToAs with the DDH and the DDGR models, and Fig. 3 shows the timing residuals. Both models are a good description of the data, with reduced χ² ≈ 0.95, and with spin, astrometric, Keplerian and PK parameters in very good consistency.

Fig. 3. Timing residuals against observing time, mean anomaly, and radio frequency of the ToAs extracted from the TempoNEST DDH fit, presented in Table 2. No significant trends are seen in the data, which is indicative of the quality of the fit. The DDGR fit presents virtually the same residuals.

The results of the timing analysis are consistent with J1208−5936 being the first-born NS in a DNS system. It is a mildly recycled pulsar (P₀ = 28.71 ms, Ṗ ≲ 4 × 10⁻²⁰ s s⁻¹; Fig. 4) in an eccentric, compact orbit (e = 0.3480, P_b = 0.6316 days) in a binary with both masses being consistent with those of NSs (M_t = 2.586(6) M_⊙, , ). Similar to most massive Galactic pulsar binaries, the system lies close to the Galactic plane, at a galactic latitude of b = 2.813°. With DM ≈ 344.4 cm⁻³ pc, and assuming the NE2001¹⁹ (Cordes 2004) or the YMW16²⁰ (Yao et al. 2017) models of Galactic electron density, the corresponding DM-inferred distances from Earth are d ≈ 8.2(1.6) kpc or d ≈ 8.5(1.7) kpc respectively and with 20% uncertainties, placing it as possibly the furthest known Galactic DNS system.

Fig. 4. P0 − Ṗ diagram showcasing all isolated and binary pulsars, and highlighting all known DNS systems and candidates (Table 1). The dashed lines indicate the characteristic age τc, the surface magnetic field Bsurf and the spin-down luminosity. J1208−5936 falls at the bottom of the mildly recycled population. All data points except J1208−5936 have been retrieved from the Australia Telescope National Facility (ATNF) website.

Due to the low significance of the spin frequency derivative , both the characteristic age τ_c and surface magnetic field B_surf of J1208−5936 are poorly constrained. Furthermore, the observed magnitude of is likely dominated by contributions from the Shklovskii effect and the Galactic acceleration field (Sect. 7.1, Table 5). However, the current constraint points towards J1208−5936 possibly having the lowest spin-down rate among all DNS systems, indicating a weak magnetic field.

With the constrained masses and orbital parameters, we numerically integrate Eqs. (5.6) and (5.7) in Peters (1964) to compute a merger time of τ_m = 7.2 ± 0.2 Gyr, where the uncertainty arises from the ranges of the individual M_p and M_c values. Therefore, J1208−5936 joins the family of DNS systems merging within the Hubble time due to the orbital decay because of gravitational wave radiation (Table 1).

4.2. Profile properties

The emission of J1208−5936 exhibits no detectable linear polarisation. Considering possible Faraday smearing, we searched the Rotation Measure (RM) space in the range of −20 000 < RM < 20 000 rad m⁻² using a 1024-channel frequency-resolved profile from the PTUSE data of the semi-coherent orbital campaign. We used both psrchive/rmfit and RMcalc²¹ but found no detection, which may be caused by the scattering of the pulse at low frequencies.

As seen in Fig. 5, a secondary component trails the brighter, primary component of the pulse at high frequencies with a relative separation between peaks of approximately 12.7 deg. This secondary component disappears in the lower half of the MeerKAT band, being absorbed into the scattered tail of the main component, but becomes more prominent at high frequencies (> 1700 MHz). In addition to the MeerKAT observations, we also have eight fold-mode Parkes/Murriyang²² observations totaling more than 15 h. They were recorded with the Ultra-Wideband (UWL) receiver, which covers the band between 704 MHz and 4032 MHz (Hobbs et al. 2020). We calibrated the bandpass with the psrchive/pac command and the standard candle observations provided by CSIRO²³, cleaned them from RFI with clfd, and then took the pulse profile between 1715 and 2305 MHz by adding all the observations. Figure 6 shows the resulting pulse profile, where the secondary component is seen to gain prominence.

Fig. 5. De-dispersed (DM = 344.258 pc cm−3) PTUSE profiles at the upper and lower half of the bands of the emission of J1208−5936, with full Stokes resolution at RM = 0 rad m−2, derived with psrchive/psrplot. No linearly polarised emission has been detected.

Fig. 6. De-dispersed (DM = 344.308 pc cm−3) UWL profile from scrunching the emission of J1208−5936 between 1715 and 2305 MHz. Only total intensity data is available.

We find good evidence for scattering in the low-frequency band of MeerKAT observations (< 1000 MHz). We divide the MeerKAT band into 32 sub-bands and fit the profile in the 10 lowest frequency sub-bands with a Gaussian function convolved with an exponential

(10)

where b stands for bins, b₀ for the Gaussian centre, Δb for the standard deviation thereof, and τ_s for the scattering time-scale. Table 3 shows the τ_s value for each frequency, and Fig. 7 shows the fit at f = 877 MHz as an example, where it is clear that the pulse shape is well described by a scattered Gaussian function. This results in a significant trend of τ_s decreasing in frequency f in the shape of the power law

(11)

Fig. 7. Fit of Eq. (10) (red line) at the bottom-most sub-band of 877 MHz (cyan line). The fit is effective at getting the general shape of the pulse profile, with little structure left in the residuals (dark blue line). Table 3 lists the best scattering timescales and Gaussian widths for all fitted frequencies.

The scattering index 2.8  ±  0.2 is lower than the 4.0 typically measured in single-component scattered pulsars (e.g. Oswald et al. 2021), which is likely indicative of a bias introduced by the double-component nature of J1208−5936. Nonetheless, our analysis provides robust evidence in favour of the presence of scattering in the pulse. Finally, from the PTUSE 1024-channel profile, the best DM measurement from the 2 March 202 to 8 March 2022 observational campaign is 344.298 ± 0.054 pc cm⁻³ from the pdmp fit at full frequency resolution.

4.3. Comparison with the known DNS population

J1208−5936 falls on the lower side of the expected P₀ and P_b relationship, having a spin period below the average (Fig. 8), but it is still consistent with the rest of the DNS population. This relationship is explained through the less efficient recycling of the primary NS resulting from a delayed RLO accretion onset in longer-orbit progenitor NS–He star systems (Tauris et al. 2017).

Fig. 8. P0 − Pb diagram of the known recycled pulsars in DNS systems population. The lines represent linear regression fits to the data points in the log space, added to aid in the visualisation of the trend.

A relationship between e and P_b has been postulated on a similar basis (Tauris et al. 2015, 2017). In this case, wide NS–He star progenitors undergo a reduced mass transfer due to the delayed RLO onset, which results in a less stripped He star and therefore in an increased mass-loss during the second supernova. Therefore, large orbital periods are associated with large eccentricities. Figure 9 shows that J1208−5936 is a high-eccentricity outlier amongst DNS systems with e < 0.5 with short orbital periods. Such outliers are to be expected due to the large spread in outcome eccentricities introduced by supernova kicks (Tauris et al. 2017), but it could also be a hint towards larger supernova kick (see Sect. 4.4).

Fig. 9. e − Pb diagram of the known pulsars in DNS systems population. The lines represent linear fits to the data points at e < 0.5 and e > 0.5, with the extended fit including unconfirmed DNS systems.

Finally, we also look at the postulated e − P₀ relationship in Galactic DNS systems. This one goes in parallel with the e − P_b relationship, and is explained from the P₀ − P_b relationship: pulsars in longer orbital periods undergo less efficient recycling, the He star undergoes greater mass loss during the supernova, and therefore high eccentricities should imply longer spin periods. This was first observed by Faulkner et al. (2005), and then theorised by Dewi et al. (2005) under the assumption of supernova kicks smaller than 50 km s⁻¹. Figure 10 compares the updated DNS population with the simulated results from Dewi et al. (2005), and it shows that J1208−5936 is once again at the high-eccentricity end but still lies within the distribution. As recently reported by Sengar et al. (2022), the simulations from Dewi et al. (2005) do not coincide with the observed high-eccentricity end of the Galactic DNS population, which would favour the different formation channel hypothesis brought forward in Andrews & Mandel (2019).

Fig. 10. e − P0 diagram of the known recycled pulsars in DNS systems population. The spread of values in observed systems is larger than predicted in some simulations (Dewi et al. 2005), with J1208−5936 staying slightly above the expectation.

4.4. Formation channel

We note in Fig. 9 that J1208−5936 has a particularly high orbital eccentricity compared to other DNS systems with e < 0.5. This could be indicative that J1208−5936 has been formed through a different channel than other systems in this group. We also pay attention to the division between e > 0.5 and e < 0.5 systems in the Galactic DNS populations in the e − P_b space (Fig. 9). According to Andrews & Mandel (2019), the high-eccentricity population can be explained by larger supernova kicks from heavier He stars progenitors of the second-born NS. Tauris et al. (2015) also explores the possible evolutionary origin of this division. In most known DNS systems, the He star loses enough mass during binary interaction to stop the possibility of nuclear fusion in the ONeMg core. In these cases, the He star with a core mass ≲1.43 M_⊙ would implode through the electron-capture process instead of following shell burning until reaching the iron core collapse (Tauris et al. 2015). On the other hand, He stars that keep a core with ≳1.43 M_⊙ would be able to reach iron core collapse, possibly leaving a heavier NS behind with a larger supernova kick and with greater mass loss (Tauris et al. 2015). In this picture, most known DNS systems that do not get disrupted form through the electron-capture channel, constituting the e < 0.5 population, while the DNS systems forming through the iron core collapse channel could create a heavier e > 0.5 population. This picture is in principle corroborated by Fig. 11, where the high eccentricity systems are also shown to contain more massive second-born NSs, consistent with the idea of larger supernova kicks being associated with larger masses.

Fig. 11. Mc − e diagram of second-born NSs with measured masses (quoting the mass derived from the χ2 mapping of DDGR solutions for J1208−5936). The line represents linear fits to the data points excluding the companion of J1208−5936, which does not yet have a sufficiently constrained mass. For PSR J1906−0746 we quote the pulsar mass instead of the companion mass because it is believed to be the second-born, un-recycled NS.

It is tempting to see J1208−5936 as a system bridging the gap between these two postulated populations. While its current eccentricity is consistent with the tail of the eccentricity distribution with a second supernova kick of 50 km s⁻¹ and a progenitor He stars mass of 3 M_⊙ like other low-eccentricity DNS systems (see Fig. 10 in Tauris et al. 2017), it is also plausible that its companion has formed through the iron core collapse channel, or that it has at least suffered a stronger supernova quick than in other e < 0.5 DNS systems formed through the electron capture process owing to a more massive He star progenitor. However, in Fig. 11 it is shown that the current uncertainties on M_c are too large to discriminate whether J1208−5936 lies within any of these two populations. Further timing in the following years will constrain the mass of the companion of J1208−5936, providing clarification on its formation channel.

5. Search for pulsations from the companion

We also searched for pulsations from the companion of J1208−5936. Detecting them would not only increase the sample of known young pulsars in DNS systems, but it would also provide much more precise mass measurements and even gravity tests in the future through the inclusion of the mass ratio as an extra constraint on top of the PK parameters, such as in the case of PSR J0737−3039A/B (Lyne et al. 2004).

We performed deep searches on de-modulated APSUSE data from 30-min and 60-min long observations from the Shapiro delay semi-coherent orbital campaign. We began by cleaning the data with presto/rfifind, de-dispersing it and integrating it into a barycentric time series at DM = 344.43 pc cm⁻³ with presto/prepdata. This DM value is slightly offset from our current best estimate presented in Sect. 4.2, owing to the fact that we are dealing with APSUSE data instead of PTUSE. Nonetheless, small discrepancies in DM do not affect the search for pulsations at P > 10 ms. We de-modulated the time series with python software pysolator²⁴, which undoes the orbital motion of the companion assuming a chosen mass ratio, allowing us to search for pulsations in its rest frame of the companion. For this, we assumed several system mass ratios chosen from 27 different inclination angles between i = 45° and i = 85° uniformly spaced in sin i, assuming the total mass of M_t = 2.586 M_⊙ from the measurement of (Sect. 3, Table 2). The spacing in sin i was chosen on the basis of introducing a maximum line-of sight accelerations discrepancy of at most 1 m s⁻² between the different de-modulations based on different sin i values.

Subsequently, two different searches were performed on this data set. Firstly, we computed the FFT of each de-modulated time series with presto/realfft and then de-reddened them with presto/rednoise. The de-reddened FFT spectra were searched with presto/accelsearch with z_max = 2. For a 30-min-long and a 60-min-long observations, this corresponds to acceleration ranges of ±1.85 m s² and ±0.46 m s² for 10-ms pulsar, with increasing ranges for increasing spin period. All candidates were then refolded with presto/prepfold.

The second search implemented a fast Folding algorithm (FFA) instead, a phase-coherent periodicity search technique on evenly sampled time series (Staelin 1969). The FFA has the advantage of being more sensitive to long-period, narrow signals than the FFT (e.g. Cameron et al. 2017). Additionally, since it does not compute a Fourier Transform, de-reddening is typically performed with a running-median filter on the time series, which is less harmful towards real low-frequency signals than removing power from Fourier bins (Singh et al. 2022). Therefore, the FFA is well-suited for searching for a non-recycled second-born pulsar in the properly de-modulated time series. We restricted our FFA search to periods above 100 ms, on the basis that shorter periods would have easily been detected by the PRESTO FFT search. We used the RIPTIDE²⁵ software (Morello et al. 2019) to perform a non-accelerated FFA search on the de-modulated time series and fold the resulting candidates. For that we implemented the following steps with our presto wrapper pipelines demodulate-search²⁶.

No significant pulsations from the companion were found. This suggests that the companion has either crossed its pulsar death line, that its radio emission is beamed away from Earth, or that the pulses are too faint to be detected with the available data and telescope sensitivity.

6. Implications for NS merger rate

The discovery of J1208−5936 and, more significantly, the increase in explored depth of the southern Galactic plane provided by the sensitivity of the MMGPS-L survey are a great opportunity to update the prediction of the observed NS coalescence events with gravitational wave observatories such as LIGO and Virgo. This is possible assuming that the observed population of pulsars in DNS systems is representative of the broad Galactic DNS population, with each new survey and discovery providing a more accurate sampling of it.

Following the methodology established in Kim et al. (2003) and also implemented in Kim et al. (2010, 2015), Pol et al. (2019, 2020), Grunthal et al. (2021), we use PsrPopPy2²⁷ to populate the Galactic field with DNS systems and simulate blind surveys on the sky. In this case, we simulate the known DNS systems merging within the Hubble time. For each known DNS system j, we seek to find the proportionality constant α_j in

(12)

where λ_j is the number of discovered systems and N_tot, j is the number of simulated systems with similar characteristics each. This α_j value is then used to compute a probability distribution for the merger rate ℛ_j

(13)

where τ_life, j = τ_age + τ_obs is the lifetime of the system, composed by its age and future observable time, and f_b is beaming factor of the pulsar (inverse of the beam sweep fraction). Typically, τ_age is defined either by the characteristic age (τ_c) of the pulsar or the time after it has exited a recycling fiducial line, while τ_obs consists either on the remaining time until the pulsar crosses the death-line or the merger time (τ_m). Finally, the rate of NS merger events in the Milky Way (ℛ_MW) is then the sum of all the ℛ_j computed for all DNS systems merging within the Hubble time, which for the probability distribution translates to the convolution of all probability distributions.

The latest update on the NS merger rate prediction (Grunthal et al. 2021) includes the following blind pulsar surveys: the Pulsar Arecibo L-band Feed Array survey (PALFA; Cordes et al. 2006), the Low-latitude High Time-Resolution Universe pulsar survey (HTRU-Low; Keith et al. 2010), the Parkes High-latitude pulsar survey (Burgay et al. 2006), the Parkes Multibeam Survey (Manchester et al. 2001), the PSR B1534+12 discovery survey with Arecibo (Wolszczan 1991) and the Green Bank North Celestial Cap survey (GBNCC; Stovall et al. 2014). These surveys used to provide a realistic coverage of the explored pulsar sky, and they all have in common that they were performed with single-dish telescopes and therefore single-beam pointings on the sky, which are straightforward to model in PsrPopPy2 as circular beams with a Gaussian sensitivity function. Now we update this picture with the MMGPS-L, the most sensitive survey in the southern Galactic plane. Its interferometric nature requires an extra layer of complexity in the simulation of pointings in the sky, which we implement in PsrPopPy2 with the addition of an extra degradation factor for coherent beams (see Appendix B for details).

The biggest uncertainty in accounting for J1208−5936 in the simulations is the computation of τ_{life, J1208}. Thanks to the preliminary mass estimates presented in Sect. 3, we compute a merger time of τ_m = 7.2 ± 0.2 Gyr for J1208−5936. However, unlike with the rest of pulsars merging within the Hubble time, there is little information on the characteristic age τ_c of J1208−5936 due to the poorly constrained Ṗ. Despite this, a reasonable assumption for τ_age can be made without sacrificing accuracy. Due to its large merger time, J1208−5936 is already the pulsar with the smallest contribution to the estimated galactic merger rate, and therefore we do not expect the choice to have a significant impact when all contributions are added. Given its small Ṗ, it is likely a pulsar comparable to PSR J1913+1102 (Table 1). We therefore pick a realistic recycling age of 2.5 Gyr. For the beaming factor, we choose f_b, J1208 = 4.59, consistent with the average of pulsars in DNS systems with actually measured f_b, i values (Pol et al. 2019), and model the pulse duty fraction of δ_J1208 = 8% from the part of the profile that was visible above the noise during the discovery. Finally, we compute its Doppler degeneration factor 0 ≤ γ_2m ≤ 1 in all surveys, assuming acceleration search for each of them (see Appendix B for more details).

For consistency with Pol et al. (2019, 2020), Grunthal et al. (2021), we simulate pulsars with a mean luminosity of ⟨log₁₀(L/[mJy kpc²])⟩= − 1.1 with a standard deviation of σ_log10L = 0.9, and a mean spectral index of ⟨Γ⟩ = 1.4 with a deviation of σ_Γ = 1. The height scale of the simulated population is set to z₀ = 0.33 kpc, with pulsars following a density distribution above and below the Galactic plane of f(z) = exp( − z/z₀). For PSR J1906+0746, we use the beam shape modelling implemented in Sect. 4 of Grunthal et al. (2021). From each DNS system j, we start with simulating only N_tot, j = 100 pulsars on the sky and simulate the surveys to see the number of discoveries λ_j. Then, on each step, the number of pulsars is increased by ΔN_tot, j = 100 and the surveys are repeated until we reach N_tot, j = 4000. Within each step, the population is simulated 100 different times, and in each simulation, the surveys are performed 100 different times, producing a loop with 10 000 iterations. All of this is done to get an averaged measurement of α_j from Eq. (12).

The updated contributions to the Galactic merger rate, along with the parameters used in the simulations and computations, are presented in Fig. 12 and Table 4. PSR J1906+0746 remains the most impactful contribution due to its short lifetime, while the contribution of J1208−5936 is the smallest due to its large lifetime. The added Galactic merger rate results in

(14)

Fig. 12. Updated probabilities distribution of the merger rates after the inclusion of J1208−5936 and the MMGPS-L survey. Top: individual contributions from each DNS population, derived from the values in Table 4 and Eq. (13). The 90% credible intervals are quoted in Table 4. Bottom: galactic merger rate probability distribution from the convolution of individual DNS distributions, with a highlighted 90% credible interval.

with a 90% confidence interval. This result is shifted downwards with respect to the Myr⁻¹ presented by Grunthal et al. (2021), owing to the inclusion of the MMGPS-L survey, because despite it being the most sensitive blind pulsar survey on the southern Galactic plane, it has only added one DNS system merging within Hubble time, therefore implying a reduction of the expected number of unseen merging DNS systems in the Milky Way from the lack of new detections.

We pay attention to the increasingly constrained NS merger rates in the Milky Way based on the observation of Galactic binaries, as well as the decreasing trend of the estimated NS merger rates in the Milky Way based on the observation of Galactic binaries. As shown in Fig. 13, the discovery of PSR J0737−3039 with its relatively short merger time brought forward a drastic increase of the predicted Galactic merger rate (Burgay et al. 2003; Kalogera et al. 2004; Kim et al. 2010) with respect to previous estimates based only on PSR B1913+16 and PSR B1534+12 (e.g., Kim et al. 2003). However, further discoveries have not increased the estimated rate in recent years, as better modelling of beam shapes and sky coverage have allowed for an increase inn the accuracy of the estimates. Pol et al. (2019) presented Myr⁻¹ at the 90% confidence with the inclusion of at-the-time newly discovered systems such as of PSR J1757−1854 and the highly relativistic PSR J1946+2052, the most constrained estimate up to that date owing to the increased coverage of the sky. Then, in Pol et al. (2020) the sensitivity of the newly included the GBNCC survey (Stovall et al. 2014) outweighed the inclusion of the newly discovered relativistic J0509+3801 DNS system, decreasing the estimated number of unseen binaries in the sky and reducing the N merger rate estimate to Myr⁻¹, with more constrained uncertainties. Grunthal et al. (2021) implemented a modified beam shape in their simulations, leading to another reduction of the total value and uncertainties. Our estimate presents a continuation of this trend, being our estimate once again reduced in comparison to the most recent estimates owing to the sensitivity of the MMPGS-L, the most important newly contributing factor. Therefore, our work confirms that NS merger rate estimates based on know electromagnetic binaries are converging into more constrained, lower values as pulsar surveys scout the sky at larger depths.

Fig. 13. Evolution of estimated NS merger rates in the Milky Way based on observable Galactic binaries during the last two decades, quoting the 90% or 95% confidence intervals (Kim et al. 2003, 2010, 2015; Kalogera et al. 2004; Pol et al. 2019, 2020; Grunthal et al. 2021).

We transform this rate into a local cosmic merger rate density and a prediction for the LIGO event detection rates by assuming that the number of DNS systems merging within the Hubble time in a galaxy is proportional to its total B-band luminosity, as it is a tracer of the star-formation rate (Kopparapu et al. 2008). This results in a local merger rate density of

(15)

which results on an upper limit of ℛ_local ≤ 515 Gpc⁻³ yr⁻¹. Assuming that LIGO has sensitivity to a range distance of D_r = 130 Mpc during the O3 run²⁸ (see Chen et al. 2021b for a definition of D_r), we predict a LIGO NS merger detection rate of

(16)

Based on the observed events between April 1st, 2019 and October 1st, 2019 catalogued in the Gravitational-Wave Transient Catalogue-2.1, a NS merger rate density of Gpc⁻³ yr⁻¹ is computed, which for detections translates into yr⁻¹ (Abbott et al. 2021), in good consistency with our estimate. This shows that the rates computed from electromagnetic binaries are in good consistency with gravitational wave observatories (see also Pol et al. 2020; Grunthal et al. 2021). With this in mind, we assume D_r = 175 Mpc to make a prediction of the rate NS merger events seen by the LIGO-Virgo-KAGRA O4 run²⁸, resulting in

(17)

which implies a prediction of the detection of events during the 18 months of the O4 run²⁸, or at most 18 events within 90% credible intervals.

7. Biases and future prospects

7.1. Proper motion and galactic field

There are three main contributions that could be biasing our mass constraints. Firstly, we consider the Shklovskii and the Galactic acceleration field effects (Shklovskii 1970; Damour & Taylor 1991). These two combine into an apparent evolution of any periodicity or length L = {P₀, P_b, x}, expressed with the equation

(18)

where D is the Doppler factor, K₀ is the unit vector from the Solar System barycentre (SSB) to the pulsar system, a_PSR and a_SSB are the Galactic acceleration vectors at the pulsar system and at the SSB, V_t is the magnitude of the transverse velocity of the system with respect to the Sun, and d is the distance to the source. Secondly, we also consider the effects on and ẋ introduced by the proper motion (PM) vector μ_t = (μ_RA, μ_Dec), along with the angle of the ascending node Ω_a and inclination angle i, expressed in Kopeikin (1996) as

(19)

and

(20)

In our case, since we do not have a measurement of the proper motion or Ω_a, these effects are not yet quantifiable, but we can nonetheless estimate the likely magnitude of the contributions. From the Galactic coordinates and DM distance of J1208−5936, we take the Galactic rotation velocity curve from Sofue (2020) and predict V_t ≃ 240 km s⁻¹ and |μ_t|≃6 mas s⁻¹. This assumption is taken based on the small magnitude of the supernova kicks introduced during the formation of the system, as otherwise the binary would have been disrupted (Tauris et al. 2017), making the Galactic flow the more dominant component of the acceleration. For the Galactic acceleration field, we use the Galactic mass distribution model presented in McMillan (2017) and extract a_PSR and a_SSB.

Our estimation of the magnitude of these contributions are listed in Table 5. We expect the Shklovskii and Galactic contributions to have similar orders of magnitude but of opposite sign, likely cancelling each other out. However, the exact values will only be computable when precise measurements of PM and d are available, and it may well be possible that either of them is dominant over the other. Comparing it with the values listed in Table 2 (Ṗ = 2.6 ± 1.0 s s⁻¹), we see that the measured Ṗ may also be heavily contaminated by either contribution. Nonetheless, it is also clear that our mass constraints should not be biased by the proper motion contribution to , as the maximum expected contribution is still three orders of magnitude smaller than our current measurement uncertainty, while the effect on Ṗ_b will eventually be able to constrain the distance through comparison with the prediction given by GR (Table 2, s s⁻¹).

7.2. Future prospects for timing

We estimate the prospects for improved mass measurements of J1208−5936 by simulating new ToAs with the TEMPO2/fake plug-in. We emulate the observing cadence outside of the Shapiro delay campaign (4 ToAs/month), our current timing precision (37 μs) and add an assumed orbital decay of Ṗ_b = −1.247533 × 10⁻¹³ s s⁻¹ from the crossing of the current values of and γ_E. We also keep h₃ and ς fixed to their current values. As seen in Fig. 2, is already providing a very tight constraint, so the evolution of the uncertainties on M_p, M_c and i is expected to be dependant on the evolution of the γ_E measurement. From Fig. 14, we expect the precision of γ_E to surpass a significance of 10σ at T = 3 yr, decreasing the individual uncertainties of M_p and M_c to ±0.1 M_⊙, and that of i to ±5 deg. Beyond this point, the uncertainties on PK parameters are likely to become dominated by correlated spin and DM noise, but assuming only white noise or a good modelling with Bayesian sampling of correlated noise as offered by TempoNEST (Lentati et al. 2014), the first reliable measurement of Ṗ_b at 3σ is reached at T = 11 yr, while 5σ is reached at T = 14 yr, and 10σ at T = 19 yr. At this time, the effects of Shklovskii and the Galactic acceleration field will already have an impact on the measurement, which would only be accountable from a good measurement of the PM. Depending on the magnitude of the PM, an independent estimate of the distance may be performed at that time. The uncertainties on the masses at those times are ±0.01, ±0.008 and ±0.005 M_⊙, and for i they are ±1 deg, ±0.5 deg and ±0.4 deg. The parameters and γ_E always remain the most constraining ones, with γ_E dominating the uncertainty and Ṗ_b providing an independent estimate of true distance by forcing it to be consistent with GR and the measured PM, which in turn will help constrain the true value of Ṗ.

Fig. 14. Fractional uncertainty evolution with timing baseline of the measured , and Ṗb PK parameters from the TEMPO2/fake simulation (solid lines: uncertainty/measurement, cyan dashed line: uncertainty/simulation input value, grey dashed lines: significance thresholds). For Ṗb, the dashed line is taken as a reference due to the low significance of the measurement. and scale with T−3/2, while Ṗb scales with T−5/2.

It is unlikely that the geometry of the system will be well determined in the future. Even at T = 20 yr, the proper motion contribution to (Eq. (19), Table 5) is expected to have a maximum possible value which is one order of magnitude below the uncertainty with our assumed PM magnitude. The value of ẋ is also unlikely to be constrained in the following decades, as it is highly degenerate with γ_E, resulting in a spurious measurement until the precession angle of ω is large enough to break the degeneracy (for a discussion of this phenomenon, the reader is encouraged to consult Ridolfi et al. 2019). Furthermore, even a detection of the PM contribution to ẋ (Eq. (20), Table 5), which is expected to be three orders of magnitude larger than the GR-predicted one (Table 2), will not break the sign degeneracy of i owing to the large distance to the system.

7.3. Future profile changes

The two components present in the profile (Sect. 4.2, Fig. 5) are likely to change both in relative amplitude and in phase separation due to the geodetic precession of the pulsar spin axis around the orbital angular momentum vector during the following years (Damour & Ruffini 1974). In GR and assuming M_p = M_c, the rate of the spin-orbit coupling-induced precession is proportional to (Barker & O’Connell 1975) as in

(21)

which in the case of J1208−5936 results in an expected value of Ω_g ≈ 0.268 deg yr⁻¹. Assuming a cone-shaped pulse from which we are seeing the maximal cross-section, a spin axis perpendicular to the angular momentum and the pulse being emitted from the equator, this would give us a minimum timescale of ≈24 yr before the pulses precess out of view. However, a different spin-orbit orientation would only extend this range upwards and the primary component will most likely be visible for a longer period of time. A detailed study of profile changes during the following years will enable a more accurate prediction from the detection of any change in the profile or lack thereof.

8. Conclusions

In this work we report the discovery and follow-up study of the MMGPS-L discovery J1208−5936. Spinning at 28.71 ms and in close orbit with another NS, this is in the tenth known Galactic DNS merging within the Hubble time. We have constrained the masses and inclination angle to , and i = 57 ± 12 deg from the mapping χ² mapping of DDGR solutions, with the tightest constraint coming from the 900σ measurement of the periastron advance . The measurement of Ṗ < 4 × 10⁻⁴ s s⁻¹ is consistent with a mildly recycled pulsar and makes J1208−5936 the pulsar in a DNS with the smallest period derivative. However, the value is likely to be biased by the Shklovskii and Galactic acceleration field contamination. Its high eccentricity is still consistent with the tail of eccentricity distribution arising from a 50 km s⁻¹ supernova kick during the formation of the companion NS (Tauris et al. 2017), but at the same time it could be indicative of a larger supernova kick caused by a massive He star. A more precise measurement of M_c in the future may clarify which is the case and help confirm the idea of two main formation channels for Galactic DNS systems depending on the supernova type. We have been unable to detect polarised emission from J1208−5936, but we have observed a faint, secondary leading component to the main pulse that becomes more prominent at high frequencies, and significantly scattered at low frequencies. We have also found robust evidence for scattering, with a scattering index of 2.8 ± 0.2, even though a better modelling may be able to provide a better measurement of the scattering index.

The merger time of τ_m = 7.2 ± 0.2 Gyr adds J1208−5936 to the family of DNS systems merging withing the Hubble time, therefore making it a progenitor of NS merger events seen by gravitational-wave observatories such as the landmark GW170817 event (Abbott et al. 2017a), making it relevant for predictions of the cosmic NS merger rate based on Galactic binaries. The performance of the MMGPS-L, the most sensitive survey in the southern sky, encouraged us to revisit these predictions. The end result provides an updated merger rate of Myr⁻¹ and local cosmic merger rate of Gpc⁻³ yr⁻¹ within a 90% confidence interval, smaller than the limits provided by previous studies on Galactic DNS systems (Grunthal et al. 2021; Pol et al. 2019, 2020) owing to the fact that the despite the high sensitivity of the MMGPS-L only one new system merging within the Hubble time has been discovered, reducing the expected number of unseen DNS systems. This continues the trend of more constrained, decreasing estimates over time as the depth of pulsar surveys increase and the modelling of pulsar beam shapes improves. The resulting prediction for the LIGO-Virgo-KAGRA O4 run is the observation of NS merger events within 90% credible intervals.

We expect the mass constraints in this system to improve significantly in the following years. Through simulations, we predict the masses and inclination angle uncertainty to be reduced to ±0.1 M_⊙ and ±5 deg with only two extra years of timing. After two decades, mass and inclination angle uncertainties can be reduced down to ±0.005 M_⊙ and ±0.4 deg, with the uncertainty always being dominated by the precision in the measurement of the Einstein delay amplitude. An eventual detection of Ṗ_b will help constrain the true distance to the system by forcing consistency with GR and the PM.

Deep surveys with new sensitive facilities such as MeerKAT, FAST or the SKA in the future will continue to provide new systems similar to J1208−5936, and increase the discovery rate of DNS systems with improved sensitivity and search algorithms. The MPIfR-MeerKAT Galactic Plane survey at S-band (Kramer et al. 2016; Padmanabh et al. 2023) will probe deep into southern Galactic plane at high radio frequencies, allowing the discovery of even more distant and faint compact pulsar binaries without being hampered by propagation effects introduced by the interstellar medium. Further in time, space-based gravitational-wave observatories like LISA will to probe tens of electromagnetically dark DNS systems with orbital periods of one hour or less (Lau et al. 2020), constraining estimates of the merger rate in the Milky Way even further.

Acknowledgments

The MeerKAT telescope is operated by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation. The Parkes radio telescope is part of the ATNF, which is funded by the Australian Government for operation as a National Facility managed by the Commonwealth Scientific and Industrial Research Organisation. We acknowledge the Wiradjuri people as the Traditional Owners of the Observatory site. SARAO acknowledges the ongoing advice and calibration of GPS systems by the National Metrology Institute of South Africa (NMISA) and the time space reference systems department of the Paris Observatory. Observations used the FBFUSE and APSUSE computing clusters for data acquisition, storage and analysis. These clusters were funded and installed by the MPIfR and the Max-Planck-Gesellschaft (MPG). All authors affiliated with the MPG acknowledge its constant support. Marina Berezina acknowledges support from the Bundesministerium für Bildung und Forschung D-MeerKAT award 05A17VH3 (Verbundprojekt D-MeerKAT). Marta Burgay acknowledges support through the research grant ‘iPeska’ (PI: A. Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. Vivek Venkatraman Krishnan acknowledges financial support from the European Research Council (ERC) starting grant ‘COMPACT’ (grant agreement number: 101078094). We also thank Alessandro Ridolfi for providing a working version of pysolator.py and for his input in Sect. 5, Norbert Wex for his comments on the interpretations of the NS merger rate results in Sect. 6, and Livia Silva Rocha and Robert Main for their general feedback on this manuscript. The data underlying this work will be shared on reasonable request to the MMGPS Collaboration.

References

Appendix A: Implementation of the DM jump fit

One of the most relevant differences between APSUSE and PTUSE is the implementation of coherent de-dispersion during recording for PTUSE. This leads to a discrepancy between the best DM in the two data sets, as APSUSE data suffers from intra-channel smearing. While this hampers the S/N of APSUSE-derived ToAs, it also introduces a best-DM discrepancy between the APSUSE-derived ToAs and the PTUSE-derived ToAs, which leads to spurious DM trends in the time series when combined, and therefore to biased fit parameters and reduced quality fit. We implement a DM jump between the two data sets to enforce the continuity of the DM variations in the data, and to ensure a good quality fit of the DM and DM1 parameters in the DDH and DDGR models, taking advantage that in some epochs the data sets overlap.

Such DM jump is not implemented in TEMPO2, and therefore we implement a manual fitting instead. The process is as follows: firstly, a global fit is implemented on the entire data set, including DM and DM1. Such fit has an unreliable DM1 value because APSUSE data dominates the beginning of the time series, while PTUSE dominates the middle and later half. Then, we take the resulting model and fit only DM and DM1 for the APSUSE and PTUSE data sets individually. This creates two individual DM models, one for each data set, represented as

(A.1)

where and are the DM value at a reference time and its constant derivative, and t the instantaneous time, which covers a diferent range for each data set. Then, the average of ΔDM = DM^PTUSE − DM^APSUSE is computed in the overlapping range, and the DM value of the APSUSE template is shifted by −ΔDM. Subsequently, the APSUSE ToAs are re-processed with this fix. This cycle is repeated several times until the ΔDM value is seen to approach zero. In our implementation, the profiles were de-dispersed originally with DM = 344.312 pc cm⁻³, and the value that sets ΔDM ∼ 0 was found to be DM^APSUSE = 344.336 pc cm⁻³.

The setup for the DDH and DDGR fits and their measurements are listed in Table 2. The timing residuals are plotted in Fig. 3.

Appendix B: Implementation of the MerKAT coherent beam pattern in PsrPopPy2

Unlike the pulsar surveys included in the previous NS merger rate estimates, the MMGPS-L makes use of the MeerKAT multi-dish telescope. It is, therefore, necessary to implement the interferometric nature of its observations, which manifest as a tiling of coherently phases beams withing the main survey beam defined by the resolution of individual antennas (Chen et al. 2021a; Padmanabh et al. 2023).

Each MMGPS-L pointing consists of a survey beam with a Full Width at Half Maximum (FWHM) of 29.85 arcmin, within which a tiling of 480 coherently phased beams is embedded (Chen et al. 2021a; Padmanabh et al. 2023). The size and orientations of these coherent beams depends on the local sky position and the antenna array configuration, and therefore they can be considered independent of the RA and DEC coordinates (Chen et al. 2021a). Thus, we consider them random in nature, with a probability distribution in semi-major and semi-minor axis that mimics the true distribution of MMGPS-L size values. The probability distributions are derived by fitting skewed Gaussian functions to the full set of MMGPS-L pointing size values, and upon simulation of random values, we see that this method is successful in mimicking typical MMGPS-L coherent beam sizes (Fig. B.1).

Fig. B.1. Randomly computed coherent beam sizes for the PsrPopPy2 simulation of the MMGPS-L survey for the NS-NS merger rate computation (orange) against the true distribution of coherent beam sizes of the bulk of the MMGPS-L pointings (blue).

A second aspect of randomness in the simulation of MMGPS-L pointings is the size of the coherent beam tiling itself. Before October 5th, 2021, the coherent tiling was circular and its radius r_c depended on the local sky position (Padmanabh et al. 2023). In most of the pointings, it was held that r_c < FWHM/2, but occasionally the opposite was true. We find that a sum of two normal Gaussian functions with standard deviations of σ₁ = 0.898 and σ₂ = 2.264 arcmin, central values of ⟨x₁⟩ = 12.546 and ⟨x₂⟩ = 12.761 arcmin and a height ratio of A₁/A₂ = 2.942 describes the distribution of the distance of the furthest coherent beam from the centre of the tiling, and take it as a distribution for r_c. After October 5th, 2021, the coherent beam tiling was instead forced to adopt an hexagonal shape with an inner radius equal to the survey beam to provide a complete coverage of the sky, which allows us to model the tiling consistently across pointings (Padmanabh et al. 2023).

This parametrisation is implemented to PsrPopPy2 by modifying the doSurvey.py script with the computation of an extra degradation factor 0 ≤ D_coh ≤ 1 for the coherent beams on top of the Gaussian degradation factor of the survey beam 0 ≤ D_sur ≤ 1 used for single-dish surveys²⁹. The process goes as following: for each computed D_sur, we randomly select whether the pointing should be considered under the older or newer rules for tiling filling (before or after October 5th, 2021) with a 0.39 or 0.61 chance, corresponding to the fraction of pointings under the old and new tiling rules, respectively. In the case that the newer tiling rules are assigned, then we assume an efficient hexagonal tiling which takes an area of

(B.1)

where the 480 coherent beams are uniformly distributed. Otherwise, the area taken by the tiling is

(B.2)

where the 480 coherent beams are also uniformly distributed, and where the 0.95 factor comes from the fact that the outermost coherent beam is typically an outlier slightly outside of the main radius. This area is then divided between the 480 beams and a random position for a potential pulsar discovery is drawn within the square of the corresponding size. Then, D_coh is a second degradation factor computed from a 2D Gaussian function with the random dimensions drawn from the coherent beam tiling size distributions (Fig. B.1) centered at the square and aligned with its sides. However, for cases from the older tiling rules in which r_c < FWHM/2, if the offset from the centre of the survey beam is found to be larger than 1.05 × r_c, then we set D_coh = 0 instead.

Finally, like in the previous works (Grunthal et al. 2021; Pol et al. 2020), each pulsar in a DNS merging within the Hubble time has a computed Doppler degeneration factor that depends on the orbital parameters, the spin period, the duration of the observations and the applied search method based on the definitions given by Bagchi et al. (2013). Like the other surveys, the MMGPS-L implements an acceleration search, and therefore we computed the parameter 0 ≤ γ_2m ≤ 1 from equation (11) in Bagchi et al. (2013) using the code³⁰ prepared in Pol et al. (2019), under the assumption of 500-second long observations. Therefore, the S/N of a simulated pulsar in PsrPopPy2 is computed as

(B.3)

For the MMGPS-L, we only consider as discoveries signals with S/N_pulsar ≥ 9.

All Tables

All Figures

Fig. 1. Orbital solution of J1208−5936 based on the Pbary(tobs) time series between 5 May 2021 and 31 August 2021. Observations with precise acceleration and/or jerk measurements were used to extract several data points (5 May: −10.5 m s−2; 11 June: +22.4 m s−2, and 16 August: +7.2 m s−2, −3.6 × 10−3 m s−3; at ϕ ≈ 0.72, 0.95 and 0.15). Assuming a canonical pulsar mass of Mp = 1.35 M⊙, the minimum companion mass computed from the mass function is Mc ≳ 1.1 M⊙ (see Table 2).

In the text

Fig. 2. Corner plot showing the constraints on the inclination angle and component masses from PK parameters and the χ2 mapping of DDGR solutions. Central plots: mass–mass diagrams portraying the DDH PK parameters constraints, with each color corresponding to a different parameter as indicated in the legend (solid lines: nominal value, dashed lines: 1σ limits), and the 1σ limits from DDGR χ2 mapping in black. The shaded grey area on the right plot represents the region excluded by the mass function (i > 90), while the shaded areas on the left plot represent the areas outside the prior for Mp (outside of 0.6 < Mp < 1.5 M⊙). The explored regions of the Mp and Mc were decided based on the limits given by the mass function and Mt (from ). Outer plots: marginalised one-dimensional probability densities for Mp, Mc and cos i from DDGR χ2 mapping, showcasing the median value (black solid line) and the 31.4%, 47.4% and 49.9% percentiles on both sides (shaded areas). The resulting mass constraints are consistent with a pair of NSs.

In the text

	Fig. 3. Timing residuals against observing time, mean anomaly, and radio frequency of the ToAs extracted from the TempoNEST DDH fit, presented in Table 2. No significant trends are seen in the data, which is indicative of the quality of the fit. The DDGR fit presents virtually the same residuals.
In the text

Fig. 4. P0 − Ṗ diagram showcasing all isolated and binary pulsars, and highlighting all known DNS systems and candidates (Table 1). The dashed lines indicate the characteristic age τc, the surface magnetic field Bsurf and the spin-down luminosity. J1208−5936 falls at the bottom of the mildly recycled population. All data points except J1208−5936 have been retrieved from the Australia Telescope National Facility (ATNF) website.

In the text

	Fig. 5. De-dispersed (DM = 344.258 pc cm−3) PTUSE profiles at the upper and lower half of the bands of the emission of J1208−5936, with full Stokes resolution at RM = 0 rad m−2, derived with psrchive/psrplot. No linearly polarised emission has been detected.
In the text

	Fig. 6. De-dispersed (DM = 344.308 pc cm−3) UWL profile from scrunching the emission of J1208−5936 between 1715 and 2305 MHz. Only total intensity data is available.
In the text

	Fig. 7. Fit of Eq. (10) (red line) at the bottom-most sub-band of 877 MHz (cyan line). The fit is effective at getting the general shape of the pulse profile, with little structure left in the residuals (dark blue line). Table 3 lists the best scattering timescales and Gaussian widths for all fitted frequencies.
In the text

	Fig. 8. P0 − Pb diagram of the known recycled pulsars in DNS systems population. The lines represent linear regression fits to the data points in the log space, added to aid in the visualisation of the trend.
In the text

	Fig. 9. e − Pb diagram of the known pulsars in DNS systems population. The lines represent linear fits to the data points at e < 0.5 and e > 0.5, with the extended fit including unconfirmed DNS systems.
In the text

	Fig. 10. e − P0 diagram of the known recycled pulsars in DNS systems population. The spread of values in observed systems is larger than predicted in some simulations (Dewi et al. 2005), with J1208−5936 staying slightly above the expectation.
In the text

Fig. 11. Mc − e diagram of second-born NSs with measured masses (quoting the mass derived from the χ2 mapping of DDGR solutions for J1208−5936). The line represents linear fits to the data points excluding the companion of J1208−5936, which does not yet have a sufficiently constrained mass. For PSR J1906−0746 we quote the pulsar mass instead of the companion mass because it is believed to be the second-born, un-recycled NS.

In the text

Fig. 12. Updated probabilities distribution of the merger rates after the inclusion of J1208−5936 and the MMGPS-L survey. Top: individual contributions from each DNS population, derived from the values in Table 4 and Eq. (13). The 90% credible intervals are quoted in Table 4. Bottom: galactic merger rate probability distribution from the convolution of individual DNS distributions, with a highlighted 90% credible interval.

In the text

	Fig. 13. Evolution of estimated NS merger rates in the Milky Way based on observable Galactic binaries during the last two decades, quoting the 90% or 95% confidence intervals (Kim et al. 2003, 2010, 2015; Kalogera et al. 2004; Pol et al. 2019, 2020; Grunthal et al. 2021).
In the text

Fig. 14. Fractional uncertainty evolution with timing baseline of the measured , and Ṗb PK parameters from the TEMPO2/fake simulation (solid lines: uncertainty/measurement, cyan dashed line: uncertainty/simulation input value, grey dashed lines: significance thresholds). For Ṗb, the dashed line is taken as a reference due to the low significance of the measurement. and scale with T−3/2, while Ṗb scales with T−5/2.

In the text

	Fig. B.1. Randomly computed coherent beam sizes for the PsrPopPy2 simulation of the MMGPS-L survey for the NS-NS merger rate computation (orange) against the true distribution of coherent beam sizes of the bulk of the MMGPS-L pointings (blue).
In the text