Binary Neutron Star Mergers: Multi-Messenger Systematics and Prospects with Next-Generation Facilities

TL;DR

This paper develops an end-to-end, forward-modeling framework that connects binary neutron-star progenitor evolution to multi-messenger observables—gravitational waves, kilonovae, and gamma-ray burst afterglows—using COMPAS, gwfast, MOSFiT, and KN/afterglow models. By convolving metallicity-dependent star-formation history with population-synthesis outputs, the authors generate synthetic BNS populations and forecast joint GW+EM detectability for future third-generation networks, revealing that typical events yield $ ho \sim 10$–$20$ and localization areas of order $\Omega_{90} \sim 10\ \mathrm{deg}^2$, with KN peak magnitudes in the $i$-band around $m_i \sim 23$–$33$; yet, the resulting joint yields are highly sensitive to population-physics parameters such as the common-envelope efficiency $\alpha_{\rm CE}$. They show how including ET and larger CE configurations improves localization and reduces parameter degeneracies, but also demonstrate that modeling systematics from progenitor evolution dominate future inferences as detector sensitivities grow. The key contribution is highlighting multi-messenger modeling systematics and mapping how progenitor uncertainties propagate into GW/EM observables and joint detection rates, emphasizing the need for integrated, cross-disciplinary analyses to maximize science in the era of thousands of BNS detections per year. The work thus provides a blueprint for designing robust multi-messenger strategies and for interpreting large future catalogs in terms of astrophysical and nuclear-physics insights.

Abstract

Multi-messenger astronomy was galvanized by the detection of gravitational waves (GWs) from the binary neutron star (BNS) merger GW170817 and electromagnetic (EM) emission from the subsequent kilonova and short gamma ray burst. Maximizing multi-messenger constraints on these systems requires combining models of the progenitors and products of BNS mergers within a single framework. Motivated by GW170817, we create a combined model that relate the progenitor astrophysics of a BNS population with their GW observability and localizability, kilonova light curves, gamma-ray burst afterglow flux, and kilonova remnant evolution. We compute the BNS merger rate by convolving metallicity-dependent star-formation history with population-synthesis predictions, and we sample realistic populations to evaluate their GW and EM observables and joint detection rates. We find that next-generation detectors will typically observe BNS mergers with GW network signal-to-noise ratios of $\sim$ 10 to 20, 90th-percentile sky areas of order $\sim$ 10 deg$^2$, and kilonova $i$-band magnitudes spanning $\sim$ 23 to 33. The variation of the merger rate with respect to the common-envelope efficiency is shown in the GW and EM observables and the resulting multi-messenger detection yield, demonstrating how uncertainties propagate into all stages of joint GW+EM forecasting. Across the models examined, no more than $\sim$ 4% of BNS mergers are detectable simultaneously by a two-Cosmic-Explorer plus one-Einstein-Telescope network and by both Roman (in a $K$-like band) and Rubin ($i$ and $g$ bands). These results show that assumptions underlying the combination of progenitor evolution and source observables will constitute key multi-messenger modeling systematics for inference of astrophysical, nuclear, and fundamental physics from future datasets.

Figures (7)

• Figure 1: Properties of our fiducial BNS population from the COMPAS rapid binary population synthesis model. Top left panel: the masses of the neutron stars in each BNS with green x's (orange o's) for binaries that do (do not) undergo an episode of Common Envelope Evolution during their progenitor evolution history. Top right panel: the chirp mass $\mathcal{M}_{\rm BNS}$, mass ratio $q_{\rm BNS}$, and symmetric mass ratio $\eta_{\rm BNS}$ of each BNS. Bottom left panel: the correlation between the progenitor initial mass $m_{\rm ZAMS}$ with the mass of the neutron star in the BNS that it formed $m_{\rm BNS}$ for the primary (blue o's) and secondary (red x's) stars. Bottom right panel: the distributions of BNS total mass $M_{\rm BNS}$, orbital semi-major axis $a_{\rm BNS}$, and merge timescale due to GW emission $t_{\rm merge}$.
• Figure 2: For a single BNS of the fiducial population with parameters $m_{1, \rm BNS} = 1.3 \rm \,M_{\odot}$, $m_{2, \rm BNS} = 1.7 \rm \,M_{\odot}$, and $t_{\rm merge} = 0.4$ Myr and assuming redshift $z = 0.1$ and inclination $\cos(\iota) = \pi/4$, the left and right panels show the resulting GW polarization states $h_+(t)$ and $h_x(t)$ and the light curves from the kilonova explosion in the $i$, $g$, and $K$ bands, respectively. Horizontal lines in the right panel correspond to the $i$ and $g$ 180 second exposure with Rubin Andreoni2024a and the $F213$ filter (similar to the $K$ band) 55 second exposure with Roman Andreoni2024b.
• Figure 3: The output of the multi-messenger pipeline for the fiducial BNS population: the BNS masses $m_{1, \rm BNS}$, $m_{2, \rm BNS}$ from COMPAS, the peak magnitudes for the KN optical $i$ and infrared $K$ bands (with the number days $t_{\rm i}$ and $t_{\rm K}$ it is within 1 magnitude of the peak) from MOSFiT, the KNR timescale $t_{\rm equiv}$ in days, and the SNR $\rho$ and sky area $90^{\rm th}$ percentile $\Omega_{90}$ from gwfast for GW detector network CE 40km + CE 20km + ET 20 km. In this case, we assume a single redshift ($z = 0.1$) and inclination angle ($\iota=\pi/4$) for each BNS, implying that the variance in the GW and EM observables is due to only the BNS masses. As a single distance is assumed here, the 1 keV X-ray afterglow at 0.1 days after the BNS merger is $\Phi_{1 \rm keV} \approx 5.02\,\mu{\rm Jy}$ for all the BNSs. The sGRB afterglow and KNR assume the same ISM density $\rho_{\rm ISM, fid} = 1.5e-25 {\rm g}/{\rm cm}^3$.
• Figure 4: Histograms of the timescale $t_{\rm equiv}$ over which an amount of ISM mass equivalent to the total $m_{\rm ejecta}$ of the KN is swept-up by the KNR for the fiducial population of BNS mergers, i.e. see Fig.'s \ref{['F:evolutions']}-\ref{['F:cornerFiducial']}. The colors blue, red, and green correspond to the three densities $\rho_{\rm ISM} = \rho_{\rm ISM, fid}$, $10\rho_{\rm ISM, fid}$, $100\rho_{\rm ISM, fid}$, and the $m_{\rm ejecta}$ of each BNS is computed with the MOSFiT BNS kilonova model.
• Figure 5: The output of the multi-messenger pipeline for the fiducial BNS population: the BNS chirp mass $\mathcal{M}_{\rm BNS}$ and symmetric mass ratio $\eta_{\rm BNS}$, the peak magnitudes $m_{\rm i}$ and $m_{\rm K}$ for the KN optical i and infrared K bands from MOSFiT, the flux density in $\mu{\rm Jy}$ of the 1 keV GRB afterglow $\Phi_{1 \rm keV}$ at 0.1 days after the BNS mergers, the KNR timescale $t_{\rm equiv}$ in days, and the GW SNR $\rho$, sky area $90^{\rm th}$ percentile $\Omega_{90}$, and fractional uncertainties on the chirp mass $\Delta\mathcal{M}/\mathcal{M}$, symmetric mass-ratio $\Delta\eta/\eta$, luminosity distance $\Delta d_{\rm L}/ d_{\rm L}$, and inclination $\Delta\iota/\iota$ from gwfast for detector networks corresponding to colors: LVK (green), CE 40km (blue), CE 40km + CE 20km (orange), CE 40km + ET 20km (red), CE 40km + CE 20km + ET 20 km (purple). We uniformly assign redshift $0.01 \leq z \leq 2$ and inclination $0.1 \leq \cos(\iota) \leq 0.9$ to each BNS, implying that the variance in the GW and EM observables is due to a complicated interplay of parameters. The sGRB afterglow and KNR assume the same ISM density $\rho_{\rm ISM, fid} = 1.5e-25 {\rm g}/{\rm cm}^3$. Also, the statistical correlations between data points displayed here are weighed by $z^2$ to mimic a more realistic population of sources.