The MPIfR-MeerKAT Galactic Plane Survey II. The eccentric double neutron star system PSR J1208 − 5936 and a neutron star merger rate update

The MPIfR-MeerKAT Galactic Plane survey at L -band (MMGPS-L) is the most sensitive pulsar survey in the Southern Hemi-sphere, providing 78 discoveries in an area of 900sq.deg. Here, we present a follow-up study of one of these new discoveries, PSRJ1208 − 5936, a 28.71-ms recycled pulsar in a double neutron star system with an orbital period of P b = 0 . 632days and an eccentricity of e = 0 . 348, merging within the Hubble time. Through timing of almost one year of observations, we detected the relativistic advance of periastron ( ˙ ω = 0 . 918(1)degyr − 1 ), resulting in a total system mass of M t = 2 . 586(5) M (cid:12) . We also achieved low-signiﬁcance constraints on the amplitude of the Einstein delay and Shapiro delay, in turn yielding constraints on the pulsar mass ( M p = 1 . 26 + 0 . 13 − 0 . 25 M (cid:12) ), the companion mass ( M c = 1 . 32 + 0 . 25 − 0 . 13 M (cid:12) ), and the inclination angle ( i = 57 ± 12 ◦ ). This system is highly eccentric compared to other Galactic ﬁeld double neutron stars with similar periods, possibly hinting at a larger-than-usual supernova kick during the formation of the second-born neutron star. The binary will merge within 7 . 2(2)Gyr due to the emission of gravitational waves, making it a progenitor of the neutron star merger events seen by ground-based gravitational wave observatories. With the improved sensitivity of the MMGPS-L, we updated the Milky Way neutron star merger rate to be R newMW = 25 + 19 − 9 Myr − 1 within 90% credible intervals, which is lower than previous studies based on known Galactic binaries owing to the lack of further detections despite the highly sensitive nature of the survey. This implies a local cosmic neutron star merger rate of R newlocal = 293 +


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 dis-ruption 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 −1 (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. 2005a;Tauris et al. 2017).
Table 1: All known pulsars in Galactic DNS systems and candidates with their spin parameters, orbital parameters, and component masses ordered by merger times τ m , including J1208−5936 (highlighted).The first solid line splits systems merging within the Hubble time from those that do not.The last two entries have no total mass measurements, and therefore their DNS nature is uncertain.Nonetheless, they are promising candidates due to their eccentricity.2021), (15) Martinez et al. (2015), ( 16) Corongiu et al. (2007), (17) last published work (Janssen et al. 2008) and section 8.3 of Tauris et al. (2017), (18) Swiggum et al. (2023), ( 19) Swiggum et al. (2015), (20) Ng et al. (2018), (21) Keith et al. (2009) Nonetheless, if all of these steps are completed without disruption, then a new, eccentric DNS system is born.

PSR
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. 2003Kim et al. , 2010Kim et al. , 2015;;Pol et al. 2019Pol et al. , 2020;;Grunthal et al. 2021).The most recent estimate (Grunthal et al. 2021) provides an upper limit of R local ≤ 597 yr −1 , 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. 2017b), opened a new era for multi-wavelength and multi-messenger astronomy, highlighting DNS systems among other astrophysical systems even further (Abbott et al. 2017a).
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 −1 , 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 observa-tion (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 Section 2, we present the discovery and the derivation of an early orbital solution.In Section 3, we show a phase-coherent timing solution and perform preliminary mass measurements with almost one year of timing data.In Section 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 Section 5, we search for pulsations from its companion.In Section 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 gravitationalwave observatories.And finally, in Section 7 we discuss possible biases and prospects for future improvements.

Discovery and orbital solution
Led by the Max Planck Institute for Radioastronomy 1 (MPIfR) in collaboration with the South African Radio Astronomy Observatory 2 (SARAO), the MMGPS-L was a fast Fourier transform (FFT)-based pulsar search survey with the MeerKAT radio interferometer 3 (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 4 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 S /N = 15.1, a line-of-sight acceleration of −8.86 m s −2 , DM = 344.2pc cm −3 and a PulsarX 5 fold of S /N ≈ 19.Upon discovery, we performed a first 20 minute-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 6 fold of   2).
S /N ≈ 12 with a changed period of P = 28.697ms, 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 SeeKAT7 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 Mosaic8 software (Chen et al. 2021).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 20minute 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 dspsr9 , where cycles of the pulsar are stacked every 8 or 10 seconds.We then cleaned these archives with the clfd10 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/pdmp11 software (Hotan et al. 2004), but also with the timing software TEMPO212 (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.py13 .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 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 pyfitorbit14 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 ⊙ .

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.py15, 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 hours.Using our phaseconnected 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 subintegration times of at most 10 minutes in length in both cases.These files had been excised of RFI a priori at full time and frequency resolution using first the clfd16 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 3 ) 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 while h 3 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 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 Pulsar The explored regions of the M p and M c were decided based on the limits given by the mass function and M t (from ω).Outer plots: marginalised one-dimensional probability densities for M p , M c 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.
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 Section 7.2 for a discussion on this).It is worth noting that we use the h 3 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 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 3 and ς instead of r and s.The fit was done with TEMPO2/TempoNEST17 , a multi-nested Bayesian sampling plugin to TEMPO2 (Lentati et al. 2014), in order to find stable values for both h 3 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.py18, 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 χ 2 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.pyexploration of the DDGR χ 2 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 This is in excellent agreement with the 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 3 .The derived 1σ constraints from each DDH PK parameter and the χ 2 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 3 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 −0.12 M ⊙ , and i DDH = 59 ± 9 deg (8) from multiplying probability distributions given by the PK parameter limits, and from the χ 2 mapping of DDGR.Additionally, as a final consistency check, we also re-derived the likely values of γ E , h 3 and ς by calculating their marginal one-dimensional probability densities from the DDGR χ 2 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 χ 2 mapping limits on M p , M c and i and assuming GR, we predict them to have values of ṖGR b = −1.225+0.026 −0.009 × 10 −13 s s −1 and ẋGR = −6.37 +0.14  −0.05 × 10 −18 ls s −1 , 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 (Section 7.1, Table 5).

General characteristics
Table 2 presents the timing parameters resulting from fits to the ToAs with the DDH and the DDGR models, and Fig. 4 shows the timing residuals.Both models are a good description of the data, with reduced χ 2 ≈ 0.95, and with spin, astrometric, Keplerian and PK parameters in very good consistency.
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 0 = 28.71ms, Ṗ ≲ 4 × 10 −20 s/s, Fig. 3) 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 ⊙ , M p = 1.26 +0.13  −0.25 M ⊙ , M c = 1.32 +0.25 −0.13 M ⊙ ).Similar to most massive Galactic pulsar binaries, the system lies close to the Galactic plane, at a galactic latitude of b = 2.813 degrees.With DM ≈ 344.4 cm −3 pc, and assuming the NE2001 19 (Cordes 2004) or the YMW16 20 (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.
19 https://pypi.org/project/pyne2001/ 20http://119.78.162.254/dmodel/index.phpDue 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 (Section 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 equations 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).

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 −20000 < RM < 20000 rad m −2 using a 1024-channel frequency-resolved profile from the PTUSE data of the semi-coherent orbital campaign.We used both psrchive/rmfit and RMcalc21 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/Murriyang22 observations totaling more than 15 hours.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,23 cleaned them from RFI with clfd, and then took the pulse profile between 1715 and 2305 MHz by adding all the observations.Fig. 6 shows the resulting pulse profile, where the secondary component is seen to gain prominence.
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 where b stands for bins, b 0 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 a Taken from the direct TEMPO2 DDGR fit.Variations from the χ 2 mapping (Section 3) are only of the order of 1/369 ≈ 0.003 (Reduced χ 2 = χ 2 / degrees of freedom).b Derived through equation 2. In the DDGR case, this simple derivation has been chosen instead of the marginalisation of the χ 2 mapping due to it being much more constraining.c The first value is the direct DDGR fit, with a Gaussian uncertainty.The second one results from the χ 2 mapping of solutions (Section 3).Spin period derivative, P (s s 1 ) 1 0 4 y r 1 0 6 y r 1 0 8 y r 1 0 1 0 y r 1 0 3 0 e r g / s 1 0 3 3 e r g / s 1 0 3 6 e r g / s Fig. 3: P 0 − Ṗ 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 B surf 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.

Comparison with the known DNS population
J1208−5936 falls on the lower side of the expected P 0 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).
A relationship between e and P b has been postulated on a similar basis (Tauris et al. 2015(Tauris et al. , 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.Fig. 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 Section 4.4).
Finally, we also look at the postulated e − P 0 relationship in Galactic DNS systems.This one goes in parallel with the e − P b

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 popula-  tion.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.
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 distribu-   ), 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.

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-minute and 60-minute 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.43pc cm −3 with presto/prepdata.This DM value is slightly offset from our current best estimate presented in Section 4.2, owing to the fact that we are dealing with APSUSE  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.
data instead of PTUSE.Nonetheless, small discrepancies in DM do not affect the search for pulsations at P > 10 ms.We demodulated the time series with python software pysolator24 , 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 ω (Section 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 −2 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-minute-long and a 60-minute-long observations, this corresponds to acceleration ranges of ±1.85 m s 2 and ±0.46 m s 2 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 riptide25 software (Morello et al. 2019) to perform a non-accelerated FFA search on the demodulated time series and fold the resulting candidates.For that we implemented the following steps with our presto wrapper pipelines demodulate-search26 .
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.

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. (2010Kim et al. ( , 2015)); Pol et al. (2019Pol et al. ( , 2020)); Grunthal et al. (2021), we use PsrPopPy227 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 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 R j 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 (R MW ) is then the sum of all the R 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 Table 4: Pulsars used in the NS merger rate computation simulation, along with their used lifetime, beaming fraction, pulse duty (δ, used only during the PsrPopPy2 runs), and the derived α i and individual merger rate contributions.For PSR J1906+0746, the lifetime is shorter than the merger time due to it crossing the pulsar deadline in ∼60 Myr.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 Section 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. (2019Pol et al. ( , 2020)); Grunthal et al. (2021), we simulate pulsars with a mean luminosity of  4 and equation 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.
⟨log 10 L/ mJy kpc 2 ⟩ = −1.1 with a standard deviation of σ log 10 L = 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 = 0.33 kp, with pulsars following a density distribution above and below the Galactic plane of f (z) = exp (−z/z 0 ).
For PSR J1906+0746, we use the beam shape modelling implemented in Section 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 = 4 000.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 equation 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 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;Kalogera et al. 2004;Kim et al. 2010Kim et al. , 2015;;Pol et al. 2019Pol et al. , 2020;;Grunthal et al. 2021).
lifetime.The added Galactic merger rate results in with a 90% confidence interval.This result is shifted downwards with respect to the R 2021 MW = 32 +19 −9 Myr −1 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 (eg., 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 R 2019 MW = 42 +30 −14 Myr −1 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 R 2020 MW = 37 +24 −11 Myr −1 , 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.
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 −237 Gpc −3 yr −1 is computed, which for detections translates into R GRTC-2.1  LIGO,O3 = 2.36 +4.69 −2.18 yr −1 (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, 28 resulting in R new LIGO,O4 = 6.73 +5.12 −2.37 yr −1 , ( which implies a prediction of the detection of 10 +8 −4 events during the 18 months of the O4 run, 28 or at most 18 events within 90% credible intervals.

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   28 https://observing.docs.ligo.org/plan/1991).These two combine into an apparent evolution of any periodicity or length L = {P 0 , P b , x}, expressed with the equation where D is the Doppler factor, K 0 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 and 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 −1 and |µ t | ≃ 6 mas s −1 .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 −1 ), 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, ṖGR b = 1.225 +0.026 −0.009 s s −1 ).

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 −13 s/s from the crossing of the current values of ω and γ E .We also keep h 3 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 For Ṗb , the dashed line is taken as a reference due to the low significance of the measurement.ω and γE scale with T −3/2 , while Ṗb scales with T −5/2 .Fig. 14, we expect the precision of γ 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 Ṗ.
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 ω (Equation 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.ẋ 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, the PM contribution to ẋ (equation 20, Table 5) is expected to have a maximum possible value one hundred times smaller than the GR-predicted one (Table 2).Therefore, a measurement of the ascending node Ω a and a breaking of the sign degeneracy of i will remain unlikely unless a large, unexpected PM magnitude is detected.

Future profile changes
The two components present in the profile (Section 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 which in the case of J1208−5936 results in an expected value of Ω g ≈ 0.268 deg yr −1 .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 spinorbit 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.

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 M p = 1.26 +0.13 −0.25 M ⊙ , M c = 1.32 +0.25  −0.13 M ⊙ and i = 57 ± 12 deg from the mapping χ 2 mapping of DDGR solutions, with the tightest constraint coming from the 900σ measurement of the periastron advance ω.The measurement of Ṗ < 4 × 10 −4 s s −1 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 −1 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. 2017b), 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 R new MW = 25 +19 −9 Myr −1 and local cosmic merger rate of R new local = 293 +222 −103 Gpc −3 yr −1 within a 90% confidence interval, smaller than the limits provided by previous studies on Galactic DNS systems (Grunthal et al. 2021;Pol et al. 2019Pol et al. , 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 10 +8 −4 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.
where the 480 coherent beams are uniformly distributed.Otherwise, the area taken by the tiling is A = π × (0.95 × r c arcmin) 2 , (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 code30 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 For the MMGPS-L, we only consider as discoveries signals with S /N pulsar ≥ 9.

Fig. 2 :
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 M p (outside of 0.6 < M p < 1.5 M ⊙ ).The explored regions of the M p and M c were decided based on the limits given by the mass function and M t (from ω).Outer plots: marginalised one-dimensional probability densities for M p , M c 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.

Fig. 4 :
Fig. 4: Timing residuals against observing time, mean anomaly, and radio frequency of the ToAs extracted from the TempoNEST DDH fit, presented in Table2.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.

Fig. 5 :
Fig. 5: De-dispersed (DM = 344.258pc 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. 7 :
Fig. 7: Fit of equation 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.

Fig. 8 :
Fig. 8: P 0 -P b 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.

Fig. 9 :Fig. 10 :
Fig. 9: e -P b 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.

Fig. 11 :
Fig. 11: M ce 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, whichdoes 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.

Fig. 12 :
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 Table4and equation 13.The 90% credible intervals are quoted in Table4.Bottom: Galactic merger rate probability distribution from the convolution of individual DNS distributions, with a highlighted 90% credible interval.
which results on an upper limit of R local ≤ 515 Gpc −3 yr −1 .Assuming that LIGO has sensitivity to a range distance of D r = 130 Mpc during the O3 run 28 (seeChen et al. 2021  for a definition of D r ), we predict a LIGO NS merger detection rate of R new LIGO,O3 = 2.76 +2.10 −0.97 yr −1 .(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 R GRTC-2.1 local = 256 +510

Fig. 14 :
Fig.14: Fractional uncertainty evolution with timing baseline of the measured ω, γE 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 γE scale with T −3/2 , while Ṗb scales with T −5/2 .

Table 2 :
Timing parameters resulting from the DDH and DDGR fits.Values without brackets come from the TempoNEST fit of DDH or the direct TEMPO2 fit of DDGR.Values with square brackets [...] are derived from the TempoNEST DDH parameters or the direct TEMPO2 DDGR parameters (M t for DDH and ω for DDGR).Values with curly brackets {...} are derived from the multiplication of PK density probabilities in DDH or the χ 2 mapping probability in DDGR.All Keplerian, spin, position and DM parameters presented for DDGR are taken from the direct TEMPO2 fit.These fits also include a time jump between the APSUSE and PTUSE data-sets, as well as a DM jump (see Appendix A for details).

Table 3 :
Gaussian width ∆b and scattering timescale τ ms , with both units of bins (out of 1024) and ms.Aside from a couple of outliers, there is a clear increase of τ s with decreasing f .
Oswald et al. 2021)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 −3 from the pdmp fit at full frequency resolution.
Table 3 lists the best scattering timescales and Gaussian widths for all fitted frequencies.

Table 5 :
Possible biases onto spin a PK parameters measurements introduced by the Shklovskii effect and the acceleration in the Galactic field (eq.18) and from proper motion (eq.19 and 20), assuming d = 8.2 kpc and |µ t | = 6 mas s −1 (see text in Section 5).