Multi-messenger observations of double neutron stars in Galactic disk with gravitational and radio waves

Binary pulsars emit both radio and gravitational waves and provide laboratories for testing gravity, exemplified by the Hulse–Taylor pulsar and the Double Pulsar. Current radio surveys struggle to discover compact DNSs with orbital periods below ~1 hour due to Doppler smearing and computational limits, whereas upcoming space-borne GW detectors (LISA, TianQin, Taiji) can access systems with orbital periods of a few minutes to ~1 hour. Motivated by the gap between radio-detected DNSs (>1 hr periods) and ground-based GW detections of merging DNSs, this study investigates whether DNSs identified via mHz-band GWs can be followed up in radio, enabling multi-messenger studies that can improve constraints on DNS demographics and neutron star equation of state.

Past work estimated LISA-detectable DNS counts using various merger rates and population models. Seto (using higher merger rates) predicted a few DNSs detectable in the Local Group at SNR >10 for 10-year observations. Lau et al. (COMPAS synthesis; 4-year LISA; merger rate 33 Myr−1) estimated ~35 detections (mostly Galactic). Andrews et al. (merger rate 210 Myr−1) found LISA could detect ~240 (330) Galactic DNSs in 4 (8) years at SNR ≥7, plus a few in nearby galaxies. Kyutoku et al. proposed LISA-informed radio searches for systems with Pb <10 min, showing large reductions in telescope pointings and orbital-modulation trials. Thrane et al. explored LISA+SKA constraints on the neutron star EOS via Lense–Thirring precession, finding potential mass–radius precision of ~0.2% over 10 years, better than current X-ray and GW-only constraints. These studies highlight both the detection prospects and scientific payoff of GW–radio multi-messenger observations, but a systematic assessment incorporating TianQin, joint networks, parameter estimation forecasts, and realistic radio follow-up yields was lacking.

• DNS population synthesis: Adopt a Milky Way merger rate R_MW = 210 Myr−1 (from LIGO/Virgo DNS rate), simulating N_DNS = 2100 systems that will merge within 10 Myr. Sources are distributed in the Galactic disk with an exponential radial profile (scale length L=2.5 kpc) and sech^2 vertical profile (scale height β=0.2 kpc), and random azimuth. Component masses are fixed at 1.4 M⊙ each (M_c = 1.22 M⊙). Orbital evolution follows Peters’ equations for eccentric binaries; initial reference near merger is based on PSR B1913+16 (Pb0=0.2 s, e0=5×10−6). Draw merger times uniformly over the past 10 Myr and compute current Pb and e by integrating the evolution equations.
• GW signal and detector modeling: Use quadrupole-order eccentric inspiral waveforms with harmonics (n=1–4) including inclination, periastron phase, polarization, annual Doppler modulation, and, when measurable, frequency derivative (Ṗb) effects. Detector responses for two independent Michelson channels are included for LISA and TianQin in the long-wavelength limit. Non-sky-averaged sensitivity curves S_L(f), S_T(f) are adopted from mission references; a combined sensitivity S_T+L(f)=S_T S_L/(S_T+S_L) is used for joint analysis (uncorrelated noises). Galactic confusion noise is included for LISA; omitted for TQ in the relevant band.
• Detection and parameter estimation: Assume 4-year observations, random orientations (cos ι ~ U[−1,1], ψ ~ U[0,π], Φ_P ~ U[0,2π]) and SNR threshold 7. For sources with f above the chirp-measurement threshold, include Ṗb in the parameter set; otherwise exclude it. Compute SNRs and Fisher Information Matrices (time-domain approximation) with central finite differences and sampling at 4 f0. Perform 100 Monte Carlo realizations to obtain distributions of detectable counts and parameter errors. Sky localization error ΔΩ is derived from the covariance matrix.
• Radio follow-up modeling: Assume each DNS contains one recycled MSP and one normal pulsar (NP). Draw MSP spin periods uniformly from 10–30 ms; NP periods from a lognormal distribution. Beam geometry: adopt period-dependent beam radius ρ with scatter and random magnetic/line-of-sight angles; pulsar detectable if |ξ−α|<ρ. Pseudo-luminosities at 1400 MHz are drawn from a lognormal distribution (applied to both MSPs and NPs). Convert to flux via S_1400=L_1400/D^2. Compute integrated-profile SNR via the radiometer equation including effective width broadening by sampling, dispersion smearing, and scattering (YMW16 electron density model for DM). Adopt SNR threshold 9 for radio detection. Telescope configurations considered: Parkes, FAST, SKA1-Mid, SKA-Mid with respective gains, bandwidths, and sky-visibility constraints.
• Timing precision: For radio-detected pulsars, estimate TOA rms uncertainty as quadrature sum of phase jitter and radiometer noise, using effective width and observing setup. Aggregate results over 100 Monte Carlo realizations.

• GW-detectable DNS counts (4-yr, SNR≥7): mean 217 (TQ), 368 (LISA), 429 (TQ+LISA) over 100 Monte Carlo runs.
• Parameter estimation (TQ+LISA): Errors generally follow power laws with SNR. Typical accuracies can reach ΔPb/Pb ~ 10−6, ΔΩ ~ 100 deg2 (improving to ~1 deg2 at SNR ~100), Δe/e ~ 0.3 (∼0.1 at SNR ~100), and ΔṖb/Ṗb ~ 0.02 (∼0.004 at SNR ~100). The orbital frequency derivative Ṗb is measurable for ~17% of sources. Strong anti-correlations exist between amplitude and inclination and between polarization and periastron phase; fixing one of a correlated pair using EM info tightens constraints on the other.
• SNR and harmonic content: For the low to moderate eccentricities produced (max e≈0.18), the n=2 harmonic dominates the SNR; calculations include n=1–4 (neglecting cross-terms).
• Derived quantities: For chirping sources, propagated errors yield small fractional uncertainties in chirp mass and distance when ΔPb, ΔṖb, and Δe are well measured, enabling distinction between DNS and typical WD–NS systems if ΔM_c/M_c ≲1%.
• Radio follow-up yields for TQ+LISA DNSs (mean over realizations): Parkes 8, FAST 10, SKA1-Mid 43, SKA-Mid 87 pulsar binaries (sums over classes MSP–NP, MSP–NS, NP–NS consistent with Table II). LISA alone yields ~50–70% more than TQ; TQ+LISA improves over LISA by up to ~17%.
• Distances of jointly detected systems (mean; kpc): Increase with telescope sensitivity. For TQ+LISA, typical means are: Parkes ~1–3.6 kpc across classes; FAST ~2.1–4.5 kpc; SKA1-Mid ~3.0–5.9 kpc; SKA-Mid ~4.7–6.8 kpc (consistent with Table III).
• Timing precision: For SKA-Mid, MSP TOA rms can be as low as ~70 ns, with most probable around ~100 μs; NPs can reach ~5 μs minimum, with most around ~130 μs. MSPs generally have TOA rms smaller by a factor ~1.3–8.2 than NPs due to narrower pulses.

The results demonstrate that space-borne mHz GW detectors will bridge the current observational gap for short-period Galactic DNSs (Pb ~3–60 min), enabling a substantial sample: hundreds detectable with LISA/TianQin over 4 years. Fisher forecasts show that critical orbital parameters (Pb, e, Ṗb) and sky location can be constrained well enough—especially for higher SNRs—to guide efficient radio searches. Incorporating GW-informed ephemerides can drastically reduce radio search parameter space and computational cost, overcoming Doppler smearing limitations that hinder blind radio searches of compact binaries. Joint GW–radio detections at the levels predicted for FAST and SKA phases will allow detailed population studies (e.g., orbital period and eccentricity distributions), improving our understanding of formation channels. Accurate measurements of Ṗb and harmonic content may also enable estimates of chirp mass, distance, and total mass (via periastron advance signatures), thereby constraining neutron star equations of state when combined with radio timing and, potentially, X-ray data. The timing precision achievable, particularly for MSPs with SKA-Mid, underpins sensitive tests of relativistic gravity and dense matter properties. Overall, the multi-messenger strategy effectively addresses the challenges of discovering and characterizing ultra-compact DNSs and maximizes their scientific return.

We synthesized a Galactic DNS population consistent with LIGO/Virgo merger rates and forecasted detections with TianQin, LISA, and their joint network over 4 years. Hundreds of DNSs are detectable in GWs, with parameter accuracies following power-law scaling with SNR; Ṗb is measurable for a non-negligible fraction. Using realistic beam geometry and empirical spin and luminosity distributions, we predict that radio follow-ups with Parkes, FAST, SKA1-Mid, and SKA-Mid can detect from a handful to dozens of these GW-identified DNSs, at typical distances of 1–7 kpc, and achieve TOA precisions down to tens of nanoseconds for the best MSPs. These results support the feasibility and promise of GW–radio multi-messenger studies of ultra-compact DNSs, enabling population mapping, precision gravity tests, and improved constraints on neutron star structure. Future work should incorporate additional low-frequency facilities (SKA-Low and pathfinders), refine pulsar spin and luminosity distributions in DNSs as new systems are discovered, and exploit higher harmonics and periastron advance to extract total masses and component properties.

Key assumptions may impact quantitative predictions: (1) The orbital reference at merger (based on PSR B1913+16) sets the eccentricity evolution; adopting PSR J0737−3039 parameters yields generally smaller eccentricities (e≲0.04) and increases Δe/e by a factor ~5 while preserving power-law trends. (2) MSP spin periods in DNSs are assumed uniform in 10–30 ms due to limited observed samples; this will be refined as more DNS MSPs are discovered. (3) Pulsar pseudo-luminosity distribution is assumed identical to normal pulsars; long evolutionary timescales could dim DNS pulsars relative to this model. (4) Only SKA1-Mid and SKA-Mid are modeled; low-frequency arrays (SKA-Low and precursors) may significantly aid searches for steep-spectrum pulsars and improve timing via DM and solar-wind mitigation. (5) Simplifications include neglect of cross-harmonic terms in SNR/FIM and omission of LISA confusion noise for TQ in the relevant band; real-data systematics and duty-cycle effects are not included.