EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 704 (December 2025) A&A, 704 (2025) A260 Full HTML
Open Access
Issue
A&A
Volume 704, December 2025
Article Number A260
Number of page(s) 18
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/202554326
Published online 15 December 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

On 17 August 2017, the first gravitational wave (GW) signal consistent with the coalescence of a binary neutron star (NSNS) system, GW170817, was detected (Abbott et al. 2017a, 2019) by the LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2015), and KAGRA (Aso et al. 2013) (LVK) GW detectors. Remarkably, less than two seconds later, a short gamma-ray burst (GRB), GRB170817A, was detected by the Fermi and INTEGRAL satellites (Abbott et al. 2017b), marking the beginning of the multi-messenger (MM) era with GWs (Abbott et al. 2017c). The presence of Virgo in the GW detector network enabled GW170817 to be localized within an area of 28 square degrees in the sky (Veitch et al. 2015; Abbott et al. 2019), prompting an extensive observing campaign that spanned the entire electromagnetic (EM) spectrum (Abbott et al. 2017c). Approximately 11 hours after the merger, a faint and rapidly evolving optical and near-infrared transient was discovered in the nearby galaxy NGC 4993 (Coulter et al. 2017). This transient was spectroscopically identified as a kilonova (KN; Pian et al. 2017), characterized by quasi-thermal emission from the merger expanding ejecta (Li & Paczyński 1998; Metzger 2019). Subsequently, a non-thermal broadband source (radio to X-rays) was detected at the same location and was identified months later as the afterglow of an off-axis relativistic jet, confirmed through very long baseline interferometry observations (Lazzati et al. 2018; Mooley et al. 2018; Ghirlanda et al. 2019).

In 2019, a second NSNS merger, GW190425, was discovered (Abbott et al. 2020a). However, no EM counterpart was identified, possibly due to the poorly constrained sky localization and the distance. No significant NSNS merger candidates have been reported in the first part of the fourth observing run (O4a) of the LVK detector network1, despite the improved GW detection range with respect to the previous runs. Therefore, although GW170817 showcased the potential of multi-messenger astronomy, it remains a single case so far. More recently, two long GRBs, GRB 211211A (Rastinejad et al. 2022; Troja et al. 2022; Mei et al. 2022) and GRB 230307A (Levan et al. 2024; Yang et al. 2024), exhibited KN signatures, suggesting an NSNS origin. Although these events were within the nominal detection horizon of GW detectors, the instruments were unfortunately offline for upgrades at the time.

Even black hole-neutron star (BHNS) mergers could be associated with EM counterparts. This is expected to occur when the distance, dtidal, at which the tidal disruption of the NS takes place is larger than the innermost stable circular orbit (ISCO) of the black hole (BH): dtidal > RISCO (Kawaguchi et al. 2015; Foucart et al. 2018, 2019; Barbieri et al. 2020). This requirement depends on several factors that are more favourable in the case of larger NS tidal deformability, higher BH spins and lower BH masses. The presence of matter outside the BH can potentially power a KN (Li & Paczyński 1998; Kawaguchi et al. 2015) and launch a relativistic jet, which may generate prompt and afterglow emission, contributing to a subclass of GRBs (Lattimer & Schramm 1974; Eichler et al. 1989; Mochkovitch et al. 1993).

To date, no EM counterpart has been conclusively associated with BHNS mergers, even though these events have been observed through GWs (Abbott et al. 2021; Abac et al. 2024). Notably, during O4a in April 2024, an event named GW230529 was reported (Abac et al. 2024). The mass of the BH in this event was 3 . 6 − 1.2 + 0.8 M ⊙ $ 3.6^{+0.8}_{-1.2}\,M_\odot $, placing it within the lower mass gap between NSs and BHs. Although the probability of EM emission from this event was limited (around 10%), the potential existence of systems with low-mass BHs increases the likelihood of EM counterparts (Xing et al. 2025), even in scenarios with non-spinning BHs or soft equations of state (EoSs), which appear to be the most favoured cases based on LVK constraints (Abbott et al. 2023).

The development of next-generation GW observatories, such as the Einstein Telescope (ET) (Punturo et al. 2010) and Cosmic Explorer (CE) (Abbott et al. 2017d; Reitze et al. 2019), is anticipated to lead to a monumental advancement in the ability to observe GW sources. These future facilities are set to revolutionize the field in the next decades by offering unprecedented detection capabilities far beyond what is achievable with current instruments. ET is planned to be constructed as an underground observatory to achieve the best insulation from seismic and environmental noise. Until recently, the leading design concept was that of a triangular ‘xylophone’ consisting of three pairs of low- and high-frequency interferometers with 10 km arms. A recent cost-benefit comparison of different designs (Branchesi et al. 2023) considered a different configuration consisting of two L-shaped interferometers with 15 km arms misaligned by 45 degrees among them and again featuring both a low- and a high-frequency instrument. The current CE design consists of an L-shaped surface-based interferometer with an arm length of 40 km, located in the United States. A possible second 20 km CE interferometer is being considered (Evans et al. 2023; Gupta et al. 2024).

The current second-generation GW detectors, including LVK and the upcoming LIGO-India, are projected to detect a few to several hundred NSNS mergers annually at design sensitivity, with a maximum observable redshift of approximately z ∼ 0.2 (Abbott et al. 2020b). By contrast, next-generation observatories such as ET and CE will significantly extend the observable volume of the Universe (Maggiore et al. 2020; Kalogera et al. 2021). ET is expected to detect NSNS mergers up to redshift z ∼ 4 (for optimally located and oriented systems), which is well beyond the peak of star formation, while CE will push the boundary even further, with a detection capability virtually extending to z ∼ 10 (Maggiore et al. 2020; Branchesi et al. 2023; Evans et al. 2021, 2023; Iacovelli et al. 2022a; Gupta et al. 2024). This remarkable increase in sensitivity and detection range will enable the observation of compact binary mergers across cosmic history up to the epoch of the first forming structures.

Beyond sheer detection rates, ET and CE will also deliver precise measurements of source parameters. For the most favourable systems, and in the presence of more than one detector, these will include improved estimates of sky localization, luminosity distance, and binary inclination angle, thereby enhancing our ability to study the astrophysical and cosmological properties of GW sources. Such advancements will not only enable more accurate source characterization but also facilitate more effective coordination with EM observatories, bolstering the field of multi-messenger astronomy (Maggiore et al. 2020; Branchesi et al. 2023).

The era of ET will coincide with significant advancements in EM facilities, enabling a deeper exploration of multi-messenger astrophysics. Among the most anticipated developments in the optical and near-infrared is the Vera C. Rubin Observatory (Ivezić et al. 2019; Andreoni et al. 2022), whose exceptional sensitivity and wide field of view make it uniquely suited to detect KNæ associated with GW events. Moreover, the introduction of highly sensitive spectroscopic telescopes, such as the Extremely Large Telescope (ELT; Marconi et al. 2022) and the Wide-field Spectroscopic Telescope project (WST; Mainieri et al. 2024), will be pivotal in studying the EM counterparts of GW events. In the radio domain, the Square Kilometre Array (SKA; Braun et al. 2019) and the Next Generation Very Large Array (ngVLA; Corsi et al. 2019) will provide remarkable sensitivity for detecting GRB afterglows. The detection of these counterparts will also benefit from new X-ray facilities such as NewAthena (Nandra et al. 2013). For the MeV and GeV bands, there are currently only proposed but not yet confirmed projects, such as the Transient High-Energy Sky and Early Universe Surveyor (THESEUS; Amati et al. 2021), aimed at enhancing high-energy transient monitoring. However, in this energy range, the new generation of AstroSats, such as the High-Energy Rapid Modular Ensemble of Satellites (HERMES; Fiore et al. 2020; Ghirlanda et al. 2024) could play a crucial role due to their lower costs and rapid development timelines, making them a valuable asset for future observations. Finally, in the very high energy (VHE; > 100 GeV) gamma-ray band, the advent of the Cherenkov Telescope Array (CTA; Cherenkov Telescope Array Consortium et al. 2019) will offer unprecedented opportunities to detect transient emission associated with the most extreme astrophysical events.

To optimize the scientific impact of next-generation GW detectors in a multi-messenger framework, it is crucial to assess the expected joint detections and the properties of the expected EM emission, in order to define the necessary instrument requirements and determine the most effective observation strategies. In this work, we present our projections for the rates and properties of NSNS and BHNS events that will be detected as multi-messenger sources in the ET era. This study is based on a state-of-the-art population synthesis model considering the emission from KN, GRB prompt, and GRB afterglow, including the VHE band, and thus it extends and improves on the methods presented in Colombo et al. (2022, 2024). The results are reported considering ET both as a standalone detector with one of the two most likely configurations and as part of a network with two CE detectors, representing the most optimistic scenario. Additionally, we also investigate the impact of the mass distribution and the NS EoS on the NSNS results, as well as the effects of the NS EoS and the BH spin on the BHNS outcomes.

The structure of the paper is as follows. In Section 2 we outline the distribution of binary parameters used in our NSNS and BHNS population models, describe the emission models for both GW and EM signals, and detail the representative multi-messenger detection limits assumed to compute the detection rates. In Section 3, we present our findings on multi-messenger detection rates as well as variations in the starting NSNS and BHNS populations. Additionally, we discuss the impact of the different network configurations of detectors. In Section 4, we analyse the expected properties of the various EM counterparts. Section 5 provides a summary discussion of our results, including a comparison with similar studies in the literature, and presents our conclusions. Throughout this study, we adopt a flat cosmology with parameters from Planck Collaboration VI (2020).

2. GW and EM population models

2.1. NSNS and BHNS populations

For the analysis conducted in this work, we considered two populations of merging compact objects, one of NSNS and one of BHNS. The NSNS population was the same as assumed in Colombo et al. (2022), in which the mass distribution was calibrated with current observational data from GW detections and Galactic NSNS systems (see Figure A.1). To analyse how the NSNS mass distribution impacts our results, we also considered two variations: a Gaussian and a uniform mass distribution (see also Schwab et al. 2010; Valentim et al. 2011; Kiziltan et al. 2013; You et al. 2025). The Gaussian has a mean value of 1.33 M⊙ with standard deviation of 0.09 M⊙ and it is based on the masses of galactic NSNS binaries (Özel et al. 2012; Özel & Freire 2016). The uniform mass distribution ranges between 1 M⊙ and MTOV, the latter representing the maximum mass for a non-rotating NS, whose value depends on the choice of the NS EoS. We assumed two different EoS models: the soft SFHo (Hempel et al. 2012) and the stiff DD2 EoS (Steiner et al. 2013). The SFHo EoS predicts a maximum non-rotating NS mass of MTOV = 2.06 M⊙ and a radius of R1.4 = 11.30 km for a 1.4 M⊙ NS, while the DD2 has MTOV = 2.46 M⊙ and R1.4 = 13.25 km. Hereafter, we call ‘fiducial’ the NSNS population constructed adopting the mass distribution from Colombo et al. (2022) and the SFHo EoS, treating the other combinations of mass distributions and EoS as variations over this fiducial model.

For all variations, we assumed the same cosmic merger rate density as in Colombo et al. (2022), derived by convolving a power law delay time distribution P(td)∝td−1 (with a minimal delay of td, min = 50 Myr, Mapelli & Giacobbo 2018; Safarzadeh & Berger 2019; Zevin et al. 2020) with the cosmic star formation rate from Madau & Dickinson (2014). This was normalized to align with a local rate density R 0 = 188 − 138 + 290 Gpc − 3 yr − 1 $ R_0 = 188_{-138}^{+290}\,\mathrm{{Gpc}}^{-3}\,\mathrm{yr}^{-1} $, obtained ensuring a consistent match with the observed frequency of significant NSNS mergers in LVK observing runs up to O4a. Given the absence of NSNS detections in the latter run, this is equivalent to the NSNS local rate density obtained in Colombo et al. (2022) (see Appendix A.1) with a correction factor that stems from the increased time-volume surveyed after O4a2.

For BHNS we considered the population described in Colombo et al. (2024). Specifically, we relied on the BH and NS mass distributions from the standard parameter set (model A) detailed in Broekgaarden et al. (2021). We also incorporated the fiducial metallicity-specific star formation rate density from the same study, grounded in the phenomenological model of Neijssel et al. (2019). As we did for the NSNS population, the rate density was normalized to be consistent with the merger rate density R 0 = 81 − 47 + 82 Gpc − 3 yr − 1 $ R_0 = 81^{+82}_{-47}\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1} $ at redshift z = 0, correcting the value presented in Colombo et al. (2024) to reflect all the BHNS events observed until the end of O4a (Abac et al. 2024)3.

The binary stellar evolution model behind our BHNS population adopts the prescriptions from the ‘delayed’ supernova explosion mechanism of Fryer et al. (2012). This results in a BH mass distribution that extends into the lower mass gap, particularly with BHs having masses below 5 M⊙. Prior to GW observations, the existence of such low-mass BHs was widely debated, largely due to X-ray binary observations in the Milky Way, which suggested a sharp cut-off at around 5 M⊙ for BH masses (Özel et al. 2010; Farr et al. 2011). However, recent GW detections indicate that the lower mass gap may not be as empty as once believed, with observations of a number of systems that likely contain a component within this mass range (Abbott et al. 2020c; Zevin et al. 2020; Abac et al. 2024; Xing et al. 2025).

As in Colombo et al. (2024), we considered two different configurations for the BH spin parameter χBH prior to the merger: a conservative approach with χBH = 0 for all binaries, and a more optimistic scenario with a uniform distribution in the range χBH ∈ [0, 0.5]. These values align with the typical spin values obtained in various simulations, assuming that the helium star progenitors of BHs that form second in the binary can spin-up through tidal interactions (Fuller & Ma 2019; Belczynski et al. 2020; Román-Garza et al. 2021a; Bavera et al. 2020, 2021, 2023). If the BH forms first, it is believed to have zero spin as a consequence of efficient angular momentum transport (Fragos & McClintock 2015; Qin et al. 2018; Fuller & Ma 2019; Belczynski et al. 2020).

In order to compute the NS compactness, we assumed the SFHo and the DD2 EoS, as we did for the NSNS population. In what follows, we adopt the most conservative setup with non-spinning BHs and the SFHo EoS as our fiducial population.

2.2. GW signal models and parameters estimation

For each merging event within our study, we computed the optimal matched-filter signal-to-noise ratio (S/N; see e.g. Chap. 7 of Maggiore 2007):

( S / N ) 2 = ∑ i ( S / N ) i 2 , ( S / N ) i 2 = 4 ∫ f min f cut | h ∼ ( i ) ( f ) | 2 S n , i ( f ) d f , $$ \begin{aligned} (\mathrm{S/N})^2 = \sum \nolimits _{i} (\mathrm{S/N})^2_i\,, \quad (\mathrm{S/N})^2_i = 4 \int _{f_{\rm min}}^{f_{\rm cut}} \dfrac{|\tilde{h}_{(i)}(f)|^2}{S_{n,i}(f)} \mathrm{d}f \,, \end{aligned} $$(1)

with the index, i, running over the considered detectors in the case of a network, fmin = 2 Hz for ET and fmin = 5 Hz for CE; fcut being a cut-off frequency depending on the events’ parameters; h ∼ ( i ) ( f ) $ \tilde{h}_{(i)}(f) $ denoting the Fourier-domain GW strain projected onto the detector, i; and Sn, i(f) as the noise power spectral density (PSD) of the ith interferometer. Additionally, we calculated the 90% credible area for sky localization, ΔΩ90%, for each detected signal. We assumed two different configurations for ET: a triangular configuration with 10 km arms (ETΔ), locating the detector in Sardinia, and a 2L-shaped interferometer design with 15 km arms misaligned by 45 deg among them (ET2L), locating one detector in Sardinia and the other one in the Meuse-Rhine region. For both the configurations, we also considered the possibility of ET operating in a global network with either one CE with 40 km arms or two CEs one with 40 km and one with 20 km arms (ETΔ+2CE, ET2L+2CE), located in the US. For ET, we adopted the same sensitivity curves as in Branchesi et al. (2023)4, while for CE we used the latest publicly available official PSDs5. Moreover, for each detector (and each interferometer in the case of a triangle) we incorporated a 85% uncorrelated duty cycle, aligning with the standard set in Branchesi et al. (2023). Given that the considered signals can stay in band as long as O ( 1 day ) $ {\cal O}(1\,\mathrm{day}) $ at 3G instruments, resulting in improved localization capabilities especially for a single detector, we include the effect of Earth’s rotation in the reconstruction as described in Iacovelli et al. (2022b,a). The main assumptions about the interferometers are summarized in Table 1.

Table 1.

Summary of the interferometer assumptions.

The computations of S/N and sky localization were performed via the public GWFAST package (Iacovelli et al. 2022b,a), adopting the IMRPhenomD_NRTidalv2 (Dietrich et al. 2019) waveform approximant for NSNS, and the IMRPhenomNSBH approximant (Pannarale et al. 2015; Dietrich et al. 2019) for BHNS. These waveform models depend on several parameters, including the detector-frame chirp mass, mass ratio, dimensionless spin parameters of the binary components, luminosity distance, sky position, binary inclination angle, polarization angle, time and phase of coalescence, and the NS tidal deformability (Iacovelli et al. 2022b). For parameters not explicitly discussed previously, values were drawn from non-informative priors within their physically relevant ranges, as detailed in Iacovelli et al. (2022b). The sky localization areas are computed within the Fisher-information-matrix formalism, valid in the high-S/N limit, in which the GW likelihood is approximated as a multivariate Gaussian near the peak (see e.g. Vallisneri (2008) for a comprehensive discussion of the formalism and its limitations). To avoid numerical instabilities in the sky position subspace due to degeneracies between some parameters (in particular distance and inclination for nearly face-on/-off systems, more likely associated with a GRB detection), we resort to a singular value decomposition as in Dupletsa et al. (2023), Ronchini et al. (2022) and eliminate from the inversion singular values below a threshold of 10−10.

2.3. Ejecta properties and EM emission models

Our framework for modelling ejecta properties and their EM emission is based on the methods presented in Colombo et al. (2022, 2024). A detailed description of how ejecta masses and velocities, and accretion disc masses are estimated for NSNS and BHNS systems using analytical formulae informed by numerical relativity is provided in these previous studies. The EM counterparts, such as KN and GRB prompt and afterglow, are modelled through semi-analytical approaches using the ejecta properties as inputs. Specifically, for each binary, we computed the KN light curves between 0.1 and 30 days in the g (484 nm), z (900 nm), and J (1250 nm) bands. The GRB afterglow light curves over 0.1–1000 days were computed in the radio, optical, and X-ray bands.

To explore the potential synergy between CTA and ET, we also calculated the afterglow light curves at 0.1 and 1 TeV over 0.01–1000 days. We assumed the VHE emission to be produced by the jet-driven external forward shock as a result of the synchrotron-self Compton (SSC) radiative process. The method we followed in computing such emission is described in Salafia et al. (2022). In brief, this model adopts a simplified description of the emissivity obtained from a delta function approximation of the average single-particle spectrum. Klein-Nishina effects were roughly accounted for by assuming a vanishing cross section for the electron-photon interaction beyond the Klein-Nishina limit. The observed flux was then computed accounting for Doppler beaming and for light-travel-time effects. The latter effects smooth out the spectra, reducing the impact of the employed approximations, producing results that deviate from a more exact treatment (e.g. Miceli & Nava 2022) only by factors of the order of one.

In computing the GRB prompt emission, we modified slightly some model parameters with respect to Colombo et al. (2022, 2024), for the following reasons. Our methodology, which is similar to that in Salafia et al. (2019), assumes that a constant fraction ηγ of the jet energy density, restricted to regions with a bulk Lorentz factor Γ ≥ 10, is emitted in the form of photons. The spectrum of these photons in the jet comoving frame is assumed to be the same at all angles. The observed time-integrated spectrum is then obtained by transforming the spectrum to the observer frame and integrating the emission over the jet’s solid angle. In our previous works, the photon flux spectrum was then calculated by dividing the time-integrated spectrum by a fixed duration for all bursts, regardless of the viewing angle. This means that the photon flux computed from such a model represents an average over the GRB duration, rather than a peak photon flux. In previous works we compared the distribution of such photon fluxes with the peak photon flux distribution of Fermi/GBM short GRBs based on light curves with a 64 ms binning. To avoid this inconsistency, in this work we compare the distribution of bolometric fluence Eiso(1 + z)/4πdL2 from our model population with that obtained from the Fermi/GBM sample of short GRBs with spectral information available in the online catalog (see Figure 1). To improve the agreement between the observed and model distributions, we decreased the ηγ parameter to 0.1 (it was 0.15 in Colombo et al. 2022, 2024). We also decreased the parameter that sets the efficiency of conversion of disc mass into jet energy (indicated with the symbol ϵ in equation B1 of Colombo et al. 2022) by a factor 2.4. The value remains well within the uncertainty limits of this parameter (Salafia & Giacomazzo 2021). We note that, while these adjustments allow for better consistency between our fiducial population and the observed GRB properties, the impact on the predictions is less than the uncertainties induced by the poorly constrained local rate density and by the large number of model parameters.

thumbnail Fig. 1.

Cumulative distribution of events above a threshold bolometric fluence. The solid line represents the distribution for Fermi/GBM observed SGRBs (with a duration below the customary 2 s threshold), constructed using the spectral information available in the online catalog. The coloured band shows the associated Poisson uncertainty. The dotted line shows our model distribution for comparison.

For NSNS systems observed within a viewing angle θv ≤ 60°, we also included a cocoon shock breakout component, modelled following the characteristics of GRB 170817A (Abbott et al. 2017b), namely a luminosity of LSB = 1047 erg/s and a cut-off power-law spectrum with νFν peak photon energy Ep, SB = 185 keV and a low-energy photon index of α = −0.62. In BHNS we did not consider this additional emission because of the lack of observing constraints from these kinds of sources.

In both the populations, we assumed a jet angular structure motivated by the one of GRB 170817A (Ghirlanda et al. 2019), characterized by a uniform core with a half-opening angle of θj = 3.4° (see again Colombo et al. 2022, for more details). However, BHNS and NSNS jets may differ due to the different environments in which they form. BHNS jets may experience less self-collimation because of the lower polar region density compared to NSNS systems (Bromberg et al. 2011; Duffell et al. 2015; Lazzati & Perna 2019; Urrutia et al. 2021; Salafia et al. 2020; Hamidani & Ioka 2021; Gottlieb & Nakar 2022; Salafia & Ghirlanda 2022). Additionally, NSNS mergers produce more isotropic ejecta and stronger post-merger winds (Foucart 2020; Kawaguchi et al. 2016; Fernández & Metzger 2013; Just et al. 2015). For this reason, in the BHNS population we also considered a broader jet opening angle of θj = 15°6, while keeping other parameters the same, adjusting only the jet core isotropic-equivalent energy Ec = E(0) to maintain the total jet energy constant.

2.4. Multi-messenger detection criteria

In order to represent the multi-messenger detection of our model sources in a simple and general way, avoiding facility-specific simulations, we opted for representing the detection condition as a threshold on integrated photon flux (for gamma-rays) or flux density (for radio and X-rays) and, equivalently, apparent magnitude for ultraviolet-to-infrared (UVOIR) bands. These thresholds are given below and are summarized in Table 2. For the GW signal detection, we imposed a network S/N threshold of 12. This stringent S/N threshold also enhances the reliability of the parameter estimation forecasts based on the GWFAST Fisher-information-matrix (Iacovelli et al. 2022a).

Table 2.

Detection limits and predicted detection rates for NSNS and BHNS assuming different detector network configurations.

Regarding EM follow-up, we anticipate substantial improvements in radio and optical search depths, owing to new instrumentation. In the radio spectrum, we expect a tenfold sensitivity increase to 0.01 mJy (Dobie et al. 2021) with respect to the representative limits we assumed for O4 (Colombo et al. 2024). Such a depth is achievable by next-generation instruments such as the SKA2 (Braun et al. 2019), Next-Generation VLA (Corsi et al. 2019), or DSA-2000 (Hallinan et al. 2019). In optical observations, advancements are anticipated with the coming on line of large FoV instruments such as the Vera Rubin Observatory (Ivezic et al. 2008). We considered magnitude thresholds of 26 in the g band and 24.4 in the z band, corresponding to expectations for Rubin Observatory’s target-of-opportunity program (Andreoni et al. 2022). For the X-ray band we assumed a flux density limit of 10−13 erg/cm2/s/keV at 1 keV, achievable by Swift/XRT.

For the GRB prompt emission, given the uncertainty on the future observational landscape, we conservatively assumed a Fermi/GBM such as instrument. We represent its sensitivity with a threshold on bolometric fluence of 3.09 × 10−7 erg cm−2, based on a visual comparison of the fluence distribution predicted by our model with the one observed by Fermi/GBM (see Figure 1). To determine the final GRB prompt emission detection rates, we accounted for the restricted field of view and duty cycle of Fermi/GBM by applying a correction factor of 0.60 (Burns et al. 2016).

Regarding the VHE afterglow band, we used the sensitivity curves of CTA North and South7 relative to photon energies of 0.1 and 1 TeV, assuming an integration time ranging from the minimum time of the light curve (0.01 days) up to 50 hours. We considered an event as detected if at least one point in its light curve was above this sensitivity threshold. We additionally assumed a 15% duty cycle to account for weather and moon constraints and a 50% reduction in sky visibility, considering that the sub-arrays are unable to observe the sky beyond a zenith angle of 60 degrees.

It is crucial to note that our analysis is predicated on the assumption that the GW sky localization areas for NSNS and BHNS mergers will be thoroughly surveyed to the outlined detection thresholds that can be obtained with extensive follow-up. A more comprehensive analysis of practical detection rates would require simulations mimicking the search strategies of individual observatories. Nevertheless, as we show, the possibility to cover the GW sky localization regions with a realistic investment of resources critically depends on the presence of multiple next-generation GW detectors in a global network.

3. Multi-messenger detection prospects

3.1. Bird’s eye view of the accessible population

Figure 2 represents a ‘bird’s eye view’ of the accessible multi-messenger populations of NSNS and BHNS mergers according to our fiducial models. The blue dot in the centre of the figure represents the Earth (not to scale). The distance scale is linear in the redshift z and two circles of constant distance (corresponding to z = 0.5 and z = 1.0) are represented with white dashed lines. Yellow (orange) dots in the left (right) half of the plot represent NSNS (BHNS) mergers that produce GW signals that exceed the assumed detection threshold in ETΔ. A red butterfly shaped symbol is drawn around GW-detectable mergers that produce a KN whose emission exceeds the assumed detection thresholds in at least one of the considered UVOIR bands; a light blue elongated jet-like symbol is drawn around GW-detectable mergers that produce a jet whose emission (prompt or afterglow) exceeds the assumed thresholds in at least one of the considered bands. The inclination of the jet and KN symbol symmetry axes with respect to the direction towards the Earth reflect those in the actual synthetic population. The number of sources represented corresponds to five hypothetical years of ETΔ operation and EM follow-up.

thumbnail Fig. 2.

Representation of two of our fiducial multi-messenger synthetic populations in a geocentric universe. Compact binaries that are detected by ETΔ are represented by yellow (NSNS – left semicircle) and orange (BHNS – right semicircle) dots. The distance of each dot from the Earth (blue circle) is proportional to the redshift of the corresponding compact binary. If the simulated merger produces a relativistic jet whose prompt or afterglow emission is detectable, according to the limits reported in Table 2, a cyan jet is plotted centered on the dot, with its axis inclined by the actual viewing angle with respect to the line of sight to the Earth. If a KN is also produced and if it is detectable, then a red butterfly shape is also plotted. The total number of binaries is representative of 5 years of ET operation.

3.2. Detection rates

In the left panel of Figure 3, we present our forecasts for the NSNS EM counterpart detection scenario in conjunction with ETΔ, based on the limits established in the previous section. The light grey line, labelled ‘All NSNS’, depicts the intrinsic cumulative merger rate within redshift z, with the corresponding light grey band illustrating its uncertainty. The latter originates from the assumed uncertainty in the intrinsic local NSNS merger rate density, and it propagates as a constant relative error to all subsequent rate estimates, displayed in the Figure as coloured bands. In this section, no other sources of error related to the initial population properties or the EM model are considered. These are discussed in detail in the following section.

thumbnail Fig. 3.

Cumulative multi-messenger detection rates as a function of redshift (luminosity distance) for our fiducial NSNS and BHNS population (SFHo EoS, non-spinning BHs), assuming the ET triangle 10 km configuration. Left panel: NSNS population. The light grey line (“All NSNS”) represents the intrinsic merger rate for the NSNS population, with the grey band showing its uncertainty due to that on the local merger rate. This uncertainty propagates as a constant relative error to all the other rates. The black (“GW ET”) line is the cumulative GW detection rate (events per year with network S/N ≥ 12). The blue (“Kilonova+GW”), red (“Afterglow+GW”), purple (“Afterglow VHE+GW”) and orange (“Prompt+GW”) lines are the cumulative detection rates for the joint detection of ET GW plus a KN (g band), a GRB afterglow (radio and VHE bands), or a GRB prompt (the orange and purple lines account, respectively, for the Fermi/GBM and CTA duty cycle and field of view). The dashed lines are the cumulative detection rates assuming only the binaries with ΔΩ90% < 100 deg2. The assumed thresholds or instruments sensitivity are shown in the legend. Right panel: Same as the left panel but for the BHNS population.

Each of the other solid lines represents a cumulative detection rate. The black line, labelled “GW ETΔ”, represents NSNS mergers with detectable GW assuming an ET triangle 10 km configuration. The other colours refer to NSNS mergers with detectable GW and one particular counterpart that exceeds the assumed detection threshold: KN in g band (“Kilonova+GW”, blue); GRB radio afterglow (“Afterglow+GW”, red); GRB prompt emission, assuming Fermi/GBM as a representative instrument (“Prompt+GW”, orange); GRB VHE afterglow, with limits representative of CTA (“Afterglows VHE+GW”). Results for other spectral bands are summarized in Table 2. The dashed lines represent the cumulative detection rates considering only binaries with a sky localization ΔΩ90% < 100 deg2.

We find that ET will detect 1 . 16 − 0.85 + 1.79 × 10 4 $ 1.16^{+1.79}_{-0.85}\times 10^4 $ NSNS mergers per year, with the 90th percentile of the redshift distribution of detectable events at z90% ∼ 1.0. Among these GW detectable binaries, 2 . 54 − 1.88 + 3.93 × 10 3 $ 2.54^{+3.93}_{-1.88}\times 10^3 $ also have a KN that exceeds, up to redshift z90% ∼ 0.4, the brightness thresholds at its peak. The large number of detectable events and the rapid evolution of these sources poses a significant challenge to the EM follow-up infrastructure. The ability of ET to localize the GW source conditions the effectiveness of the EM search for a counterpart, by determining the number of telescope pointings needed to adequately cover the GW localization region (Branchesi et al. 2023). For this reason, in Table 2 and Figure 3 we also report the detections rates for the best localized binaries, having ΔΩ90% < 100 deg2 or ΔΩ90% < 10 deg2. The rate of GW-detectable NSNS mergers with a potentially detectable KN that are also localized to within such an accuracy decreases to 104 − 76 + 160 $ 104^{+160}_{-76} $ yr−1, which represents a more realistically manageable number.

The high fraction (98%) of off-axis jets within the population makes the predicted joint GW and GRB rates comparatively lower, with a detection rate of 15 − 11 + 23 $ 15^{+23}_{-11} $ yr−1 for GRB radio afterglows and 21 − 13 + 53 $ 21^{+53}_{-13} $ yr−1 for GRB prompt. For events localized within 100 deg2, these rates decrease to 2 . 8 − 2.1 + 4.4 $ 2.8^{+4.4}_{-2.1} $ yr−1 and 1 . 3 − 1.0 + 2.1 $ 1.3^{+2.1}_{-1.0} $ yr−1, respectively. According to this analysis, the majority (∼54%) of short GRBs (SGRBs) detectable by Fermi/GBM will have an associated detectable GW signal, in agreement with the estimate from Ronchini et al. (2022). This makes the GRB prompt emission the most promising probe for multi-messenger astronomy in the distant Universe. In the cumulative distribution of the GRB prompt emission, an increase in rates can be observed at low redshifts, corresponding to events dominated by the cocoon shock breakout component. This component becomes relevant for events occurring at distances closer than about 100 Mpc.

Regarding the VHE afterglow band, the detection rates are not promising, with values around 0 . 03 − 0.02 + 0.04 $ 0.03^{+0.04}_{-0.02} $ yr−1. This low rate is primarily due to the abundance of off-axis events in the population, which produce flux levels that are too low compared to the projected sensitivity of CTA. Of course, this rate is dependent on the microphysical parameters of the afterglow model and the assumptions about the average density of the circum-burst medium. In particular, increasing the circum-burst medium density to n = 0.1 cm−3 produces detection rates higher by a factor of 10, as discussed in the following section.

In the right panel of Figure 3, we report the same cumulative detection rates for the BHNS population. We expect a GW detection rate of 1 . 39 − 0.81 + 1.41 × 10 4 $ 1.39^{+1.41}_{-0.81}\times 10^4 $ yr−1, again with z90% ∼ 1. In our fiducial population, only about 2% of the mergers produce some mass remnant (mout > 0), potentially powering an EM counterpart. This results in significantly lower rates of multi-messenger detectable events compared to the NSNS population. In particular, we predict a detectable KN rate of 20 − 11 + 20 $ 20^{+20}_{-11} $ yr−1, that decreases to 2 . 6 − 1.5 + 2.6 $ 2.6^{+2.6}_{-1.5} $ yr−1 considering only the events with ΔΩ90% < 100 deg2. The KN horizon is the same as the NSNS case, z90% ∼ 0.4, determined by the assumed sensitivity of EM facilities.

The GW-detectable systems with an observable radio afterglow and GRB prompt are forecasted to achieve a total rate of 0 . 21 − 0.12 + 0.21 $ 0.21^{+0.21}_{-0.12} $ yr−1 and 0 . 43 − 0.24 + 0.43 $ 0.43^{+0.43}_{-0.24} $ yr−1 ( 0 . 08 − 0.05 + 0.08 $ 0.08^{+0.08}_{-0.05} $ yr−1 and 0 . 06 − 0.04 + 0.06 $ 0.06^{+0.06}_{-0.04} $ yr−1 for the events with ΔΩ90% < 100 deg2). The GRB prompt horizon is smaller than the GW one. This occurs because only events with low BH masses, and hence an intrinsically fainter GW signal, lead to a NS disruption outside of the BH ISCO. Consequently, the combined GW+GRB detection horizon is defined by the GW detection of events involving BHs below a specific mass threshold. For radio afterglows, on the other hand, the value at which the curve saturates is still set by the assumed EM detection limit.

Although the detection rates for the BHNS population are significantly lower compared to NSNS cases, it is important to stress here that we are considering the most conservative population scenario, with non-spinning BHs and a soft EoS for NSs. In the next section we discuss some variations in the progenitor BHNS population to explore also more optimistic scenarios. Nevertheless, the detection rates predicted in our fiducial set up increase by more than a factor of 10 with respect to those of the O5 run (see Colombo et al. 2024). This suggests that we may have to wait for next-generation GW detectors for the first identification of EM counterparts from these types of sources: thus, ET could represent a pivotal advancement for multi-messenger astronomy in the context of BHNS observations.

3.3. Variations in the NSNS and BHNS populations

The results shown in Figure 3 refer to the NSNS and BHNS populations that we have defined as fiducial. In Figure 4 we examine the impact on the detection rates of changing some initial population assumptions. For clarity, we only show the 90th percentile of the redshift distributions of multi-messenger detectable events, assuming the same EM bands and detection limits as in the previous section, without constraints on the sky localization. The error bars always include the uncertainty due to the merger rate density. In some cases, described below, they also include some additional sources of uncertainty.

In the left panel of Figure 4, we show the variations in the NSNS population, assuming the ET triangular configuration with a 10 km arm length. The different markers represent different mass distribution choices: triangles denote the fiducial distribution used in Figure 3, while squares and circles represent a Gaussian distribution centered at 1.33 with a standard deviation of 0.9, and a uniform distribution between 1 and MTOV, respectively. A filled marker represents a result obtained assuming the SFHo EoS, while an empty marker indicates the DD2 EoS has been assumed.

thumbnail Fig. 4.

Predicted 90th percentile of the cumulative multi-messenger detection rates for our NSNS and BHNS population model variations. Left panel: NSNS fiducial population, assuming ET triangle 10 km configuration and EM bands and detection limits reported in Figure 3. Different colours refer to different counterparts as in Figure 3. Different marker shapes indicate different adopted mass distributions (triangle: fiducial; square: uniform; circle: Gaussian). Filled markers indicate the SFHo EoS, while empty markers are for the DD2 EoS. The error bars indicate the uncertainty on the local merger rate and for the “Afterglow VHE+GW” channel also a variation on the median circum-burst density. Right panel: BHNS fiducial population, assuming ET triangle 10 km configuration and EM limits reported in Table 2. Different marker shapes indicate different adopted BH spin distributions (triangle: χBH = 0; square: uniform between 0 and 0.5). Filled markers indicate the SFHo EoS, while empty markers are for the DD2 EoS. The error bars indicate the uncertainty on the local merger rate. For GRB afterglow and prompt they also take into account a variation on the jet core half-opening angle (θj = 15°) and for the Afterglow VHE also a variation on the median circum-burst density.

Regarding GW detections, represented by black markers, we observe an increase in detection rates as we move from the Gaussian distribution to the fiducial and uniform distributions, which correspond to more massive NSs on average. The variation in EoS results in only marginal changes, as expected for the inspiral signal.

For KN+GW detections (blue symbols), the only noticeable effect is an increase in the horizon, from z90% ∼ 0.4 to z90% ∼ 0.6, when assuming the uniform mass distribution, due to an increased fraction of events with a massive accretion disc, whose winds produce a strong KN emission.

In the Afterglow+GW and Prompt+GW channels (red and orange symbols, respectively) we observe differences of up to a few orders of magnitude in the detection rates depending on the population assumptions. For the SFHo EoS (filled symbols), the rates do not differ appreciably between the Gaussian and fiducial mass distributions, but they increase when assuming the uniform distribution, due to the increased fraction of events with massive accretion discs (see Figure A.1 for a comparison between mass distribution and ejecta and accretion disc masses). On the other hand, for the DD2 EoS (empty symbols), the rates vary significantly with different mass distributions. This is because our jet-launching conditions require both a non-zero accretion disc and the formation of a BH within a timescale much shorter than the accretion timescale. This latter condition is enucleated in the requirement that Mrem > 1.2 MTOV. In the case of DD2 we have MTOV = 2.46 M⊙, while for SFHo it is only MTOV = 2.06 M⊙. Therefore, when assuming DD2, the fiducial and Gaussian distributions mostly produce supramassive or stable NS remnants, leading to a reduced jet production rate (see Figure A.1, where the HMNS formation condition is indicated by a pink line).

A similar trend is visible in the VHE afterglow band (purple markers), with significant spread for the same reasons. In this case, the error bars also include the effects of a variation in the assumed median circum-burst density, from the fiducial value n = 5 × 10−3 to a higher n = 0.1 cm−3. The latter more optimistic assumption leads to an increase in detection rates by an order of magnitude.

It is important to note that our assumed jet structure and the model for the GRB prompt emission were calibrated to reproduce the bolometric fluence distribution observed by Fermi/GBM. The variations in the mass distribution and EoS naturally alter the fluence distribution, preventing a match with SGRB observations unless the jet and prompt emission models are adjusted. Therefore, the changes in the detection rates for jet-related emission should be understood as an exploration aimed at clarifying the model’s dependences and how these factors influence the results.

In the right panel of Figure 4, we show the variations in the BHNS mergers population. As in the previous case, filled and empty markers represent the SFHo and DD2 EoS, respectively, while triangles and circles indicate the assumed BH spin distribution, either χBH = 0 or a uniform distribution between 0 and 0.5. In this case, we did not consider variations in the BH mass distribution. While it is clear that different assumptions on the mass distribution can affect the multi-messenger detection rates, we focused on the uncertainties associated with the BH spin, the NS EoS, and the overall merger rate. These factors, as well the systematics in the modelling, tend to dominate the impact on the fraction of events producing EM emission, as shown in Broekgaarden et al. (2021), Román-Garza et al. (2021b), Biscoveanu et al. (2023), Chen et al. (2024), Xing et al. (2024). We note, however, that the role of BH spin must be interpreted carefully. While our models include a wide range of spin assumptions to bracket potential outcomes, realistic spin distributions from isolated binary evolution are expected to favour low BH spins, with only a small fraction of systems potentially reaching high spins due to formation order reversal (e.g. Fuller & Ma 2019; Xing et al. 2024). Therefore, the impact of BH spin on EM detectability may be somewhat overestimated under more optimistic spin assumptions. Additionally, we acknowledge that recent observational results challenge some long-standing assumptions about the BH mass distribution. In particular, the detection of GW230529 (Abac et al. 2024) not only questions the existence of a lower-mass gap–reproduced by the ‘rapid’ SN model of Fryer et al. (2012)–but also suggests that the typical BH mass in BHNS systems may be lower than previously assumed. While this potential shift in the BH mass distribution is not explicitly modelled in our current work, it is an important direction for future studies, since it could further enhance EM detectability in the BHNS population.

The GW detection rate is essentially insensitive to the variations considered in the EoS and χBH distribution. On the other hand, the impact of varying population assumptions on EM+GW detection rates is larger compared to the NSNS case. For KN+GW detections, we observe an increase in the rate by more than an order of magnitude, which also extends the redshift range at which these events can be detected. This is due both to the change in the EoS – where DD2 leads to more deformable NSs and, consequently, larger ejecta masses – and to the inclusion of the more optimistic spin assumption, which increases the fraction of events capable of producing an EM counterpart (see Figure A.2 for a comparison of the mass distribution with the ejecta and accretion disc masses). The error bars for these channels, in addition to the uncertainty related to the local merger rate, include the uncertainty due to variations in the jet opening angle, ranging from 3.4 to 15 degrees. In the case of the VHE afterglow band, we also include the effect of the variation in the median circum-burst density. Similar to the NSNS case, this allows for an increase in the detection rate by roughly a factor of ten.

A similar trend can be observed in the jet-related emission, where, unlike the NSNS population, there is no longer a dependence on the conditions for launching a jet, since all events result in the formation of a BH. The strong dependence of multi-messenger detection rates on the initial population properties demonstrates how future multi-messenger observations could provide important constraints on the population and, consequently, on its formation channels.

3.4. GW sky localization

Figure 5 displays the ΔΩ90% sky localization distribution as a function of redshift for the NSNS population (top panel, blue) and the BHNS population (bottom panel, green), considering only events capable of producing an EM counterpart. This condition requires the presence of non-zero ejecta and/or an accretion disk. The shaded regions show the extent of the 50%, 90%, and 99% confidence intervals at each redshift, with the median indicated by a solid line. Panels in different columns refer to different configurations: ET triangle configuration (ETΔ, first column), ET 2L configuration (ET2L, second column), ETΔ combined with two CE detectors (ETΔ+2CE, third column), and ET2L combined with two CE detectors (ETL+2CE, fourth column).

thumbnail Fig. 5.

Distribution of the extent of the 90% sky localization region of GW-detected NSNS events (upper row) and BHNS events (lower row) that are capable of producing an EM counterpart, as a function of redshift. Different columns refer to ETΔ (first column), ET2L (second column), ETΔ+2CE (third column), and ET2L+2CE (four column). The solid line represents the median, while the coloured bands encompass the 50%, 90%, and 99% credible interval at each fixed redshift.

The figure shows how sky localization deteriorates as redshift increases. For the triangle configuration, at redshifts below approximately z ≲ 0.04, more than 50% of events are localized within 100 deg2. However, beyond z ≳ 0.2, more than 50% of events exhibit sky localization areas exceeding 1000 deg2. In contrast, the 2L configuration generally shows improved localization capabilities, with median values reduced by nearly an order of magnitude. Specifically, at z ≲ 0.04, the majority of events are localized to within 1 to 10 deg2, while at redshifts greater than z ≳ 0.5, 50% of events have localization areas larger than 1000 deg2. The broader distribution of sky localization areas in the 2L configuration, compared to the triangle configuration, is due to instances where only a single detector is active, as dictated by the assumed duty cycle.

Considering ET as part of a network with two CE detectors, the sky localization improves even more significantly, as expected. Specifically, for both the NSNS and BHNS populations, 50% of the events within a redshift of z = 1 achieve a sky localization better than 100 square degrees. In this case, the differences between the two ET configurations become negligible.

3.5. ET configurations and detector network

In order to show how our predicted detection rates depend on the chosen EM detection thresholds, as well as different ET configurations and detector networks, we provide in Fig. 6, the detection rates as functions of EM detection limits for the joint KN+GW channels analysed in this study, for ETΔ, ET2L, ETΔ+2CE, and ET2L+2CE, corresponding respectively to the coloured, light grey, dark grey, and black lines. Panels in the top row refer to the NSNS population, while the bottom row refers to the BHNS population. The blue, green, and red colours correspond to the g, z, and J bands, respectively. The solid line considers all GW detected events, the dashed line represents events with a sky localization ΔΩ90% < 100 deg2, and the dotted line corresponds to events with ΔΩ90% < 10 deg2.

thumbnail Fig. 6.

Kilonova detection rate as a function of the EM detection limit threshold for our fiducial NSNS (upper panels) and BHNS (lower panels) populations. The blue, green, and red colours indicate the g, z, and J bands, respectively, assuming the ETΔ configuration. In each panel we also report in grey and black the ET2L and ETΔ+2CE configurations. The solid line indicates all the detectable binaries; the dashed and dotted lines show the detectable binaries with ΔΩ90% < 100 deg2 and the ones with ΔΩ90% < 10 deg2, respectively.

For a fixed detector network and band, one can observe an increase in the detection rate as the considered magnitude limit increases, until reaching saturation, which corresponds to the detection of all KNæ associated with those GW events. The plot also allows us to compare the impact of different ET configurations and detector networks on rates and sky localizations. Assuming a magnitude limit of 26 in the g band for the NSNS population, the rate of events with ΔΩ90% < 10(100) deg2 increases by approximately a factor of 3.5 (3) in the ET2L configuration. Assuming instead ETΔ in a network with two CEs, the rate increases by a factor of 400 (30). In particular, the differences between the two ET configurations become negligible in a network with two CEs, as the source distance is limited by the considered EM detection limit. For higher magnitude limits, the differences become more significant, as they correspond to larger EM horizons. In Appendix B we report the same figures for the joint GRB afterglow+GW and GRB prompt+GW channels, for which similar considerations apply.

4. EM properties

4.1. Kilonova

In Figure 7, we display the distribution of KN apparent AB magnitude as a function of time in the g band for ET-detectable (assuming ETΔ) binary systems in our NSNS (left column) and BHNS (right column) populations. Filled regions show the ranges encompassing 50%, 90%, and 99% of the brightness at each fixed time, considering all ET-detectable binaries (top row) or only the binaries with ΔΩ90% < 100 deg2 (bottom row). For comparison, we also report the observed data of AT2017gfo (Villar et al. 2017) at the median distance of the GW-detectable events: 3.5 Gpc for NSNS (1.0 Gpc for the well-localized events) and 2.4 Gpc for BHNS (1.3 Gpc for the well-localized events).

thumbnail Fig. 7.

Distribution of optical KN brightness as a function of time for ET-detectable (assuming ETΔ) events in our NSNS (left column) and BHNS (right column) populations. The shaded regions contain 50%, 90%, and 99% of the KN light curves in the g (484 nm) band. Coloured circles show extinction-corrected AT2017gfo data rescaled to the median distance of the considered populations (upper limits are marked with a downward arrow, data from Villar et al. 2017). The bottom row shows the result restricted to binaries with ΔΩ90% < 100 deg2.

For the NSNS population we find that, considering all the ET-detectable binaries, the apparent AB magnitude of the KN at peak spans from 23 down to 29, with 50% being concentrated in the 24.5–27.5 interval. When applying the constraint on the sky localization, the majority of the peaks are found in the 23.5–25.5 interval, making almost all the KNæ accessible to the Rubin Observatory.

Considering all simulated BHNS binary mergers detectable by ET which do not produce a direct plunge of the NS, the KN peak apparent AB magnitudes span the range 24.5 to 32, with 50% clustered between 27 and 29. When focusing on those with ΔΩ90% < 100 deg2, most of the peaks fall within the 25–27.5 range. Hence, we predict fainter KNæ on average with respect to NSNS. This is also made apparent by the comparison to GW170817 observational data in Figure 7.

The plot also shows the rapid KN brightness decline, particularly for the BHNS population. A possible strategy to counteract the challenge of rapidly dimming KNæ involves directly seeking out non-thermal counterparts, such as radio afterglows, using tiling instruments, as discussed in Colombo et al. (2024). This approach is well-suited for upcoming radio surveys, enabling another possible search for these transient events.

4.2. GRB afterglow

In Figure 8, we illustrate the properties of GRB afterglows associated with ET-detectable NSNS (left panel) and BHNS (right panel) binaries. The figure presents contours encompassing 50% (solid lines) and 90% (dashed lines) of the peaks of GRB afterglow light curves in the radio (1.4 GHz, red), optical (g band, green) and X-rays (1 keV, blue). Most of the peaks occur beyond 102 days. We note that light curve calculations were restricted between 10−1 and 103 rest-frame days. To enhance visualization of the underlying light curve behaviour, we included a sample of randomly chosen optical light curves (thin grey lines) in the background. For context, we also include GRB170817A data (Makhathini et al. 2021, small circles) rescaled to the median distance of the simulated populations.

thumbnail Fig. 8.

Distribution of brightness versus time for GRB afterglows associated with ET-detectable NSNS (left panel) and BHNS (right panel) mergers. Solid and dashed contours contain 50% and 90% of the peaks, respectively. Red, green, and blue colours refer to the radio (Fν at 1.4 × 109 Hz), optical (AB magnitude at 4.8 × 1014 Hz), and X-rays (νFν at 1 keV), respectively. The coloured circles are the observed data of GRB170817A (Makhathini et al. 2021) at the median distance of our population. The grey lines in the background are a random sample of optical light curves from the underlying population. We assume for both the populations the same jet half-opening angle, θj = 3.4 deg. In the BHNS panel, we also show with orange solid lines the 50% contour for radio peaks when assuming θj = 15 deg.

These observations are largely influenced by the strong dependency of GRB afterglow light curves on the viewing angle, combined with the viewing angle distribution skewed by GW detection (which favours somewhat smaller angles compared to an isotropic distribution, with a peak at about 30° – Schutz 2011). As a result, the majority of peaks occur several months to years post-GW, with a minority peaking earlier (around hours) in the optical and X-rays, exhibiting bright emission due to closer viewing angle.

In the right panel, referring to the BHNS population, we also include for comparison the region containing 50% of the afterglow peaks in the radio band assuming a jet opening angle of θj = 15 deg, represented by the solid orange line. In this case, more events are expected to be observed within the jet core or close to its border. As a result, the 50% region now comprises light curve peaks at < 1 day. The impact on the optical and X-ray bands (not shown here solely to avoid clutter) is qualitatively similar.

4.3. GRB prompt

In Figure 9, we present the distribution of rest-frame spectral energy distribution (SED) peak energy Epeak versus isotropic-equivalent energy Eiso for NSNS (left panel) and BHNS (right panel) events that meet our detection criteria for both the GW signal and the GRB prompt emission.

thumbnail Fig. 9.

Rest-frame SED peak photon energy Epeak versus the isotropic-equivalent energy, Eiso, for our NSNS (left panel) and BHNS (right panel) populations. The filled green coloured regions contain 50%, 90%, and 99% of the binaries both GRB prompt- and ET-detectable. The magenta lines contain 50%, 90%, and 99% (solid, dashed, and dotted, respectively) of the GRB prompt-detectable binaries. The black dashed line contains 90% of the ET-detectable binaries. The black dots with error bars represent an SGRB sample for comparison (Salafia et al. 2023). The orange dot is GRB170817A. We assumed for both populations the same jet half-opening angle θj = 3.4 deg. In the BHNS panel, we also show with orange solid lines the 50% contour for radio peaks when assuming θj = 15 deg.

The green shaded areas encompass 50%, 90%, and 99% of the binary systems detectable by both Fermi/GBM and the ET. The magenta lines, varying in style from solid to dashed to dotted, represent the 50%, 90%, and 99% containment levels, respectively, of binary systems that are detectable through GRB prompt emission without the requirement of ET detection.

To contextualize Fermi/GBM-detected GRB prompt events within the broader cosmological population, we incorporate a sample of SGRBs with known redshift, symbolized by grey diamonds (Salafia et al. 2023). The specific case of GRB 170817A is marked as an orange diamond on the plot, serving as a reference point.

A dashed black line is included to denote the 90% confidence region for binaries detectable by ET, independent of the GRB prompt emission detectability constraint. This curve displays the most significant differences between the two populations. For NSNS, the curve shows a bimodal distribution, as we also consider a cocoon shock breakout component for events with a viewing angle of less than 60 degrees, replicating the properties of GRB170817A. This component is significant for particularly close events (within 100 Mpc), thus it is not as relevant in comparison to the cosmological population of SGRBs. Since our BHNS population lacks this cocoon shock breakout component, the black curve does not intersect with the cocoon shock breakout event cluster and extends into lower Eiso values.

For the BHNS population, we also report a dashed orange line corresponding to 90% of binary systems that are detectable through GRB prompt emission, assuming a jet half-opening angle θj = 15 deg. Varying the jet core half-opening angle while keeping the total jet energy fixed primarily results in a shift of the Eiso distribution: wider jets correspond to lower on-axis Eiso values, as this quantity scales approximately as θj−2. As a result, populations with systematically larger opening angles would trace a distinct locus in the Epeak–Eiso plane compared to BNS-associated GRBs, assuming similar intrinsic emission properties.

5. Summary and conclusions

In this study, we have explored the multi-messenger detection prospects for NSNS and BHNS mergers in the era of ET. By leveraging state-of-the-art population synthesis models and incorporating both GW and EM emission, we have provided detailed forecasts of detection rates and outlined the observable characteristics of the EM counterparts, including KN, GRB prompt, and GRB afterglow. Our analysis considered ET as a standalone detector with different configurations as well as its potential integration within a global network that includes two CE observatories.

Our findings indicate that ET will significantly enhance the ability to observe NSNS mergers. We project that ET could detect more than 104 NSNS mergers per year, with thousands of these events also producing detectable KN emission.

The rate of events with detectable jet-related emission is on the order of a few tens per year, primarily due to the large number of off-axis events. Specifically, we find that most SGRBs currently observed by Fermi/GBM are likely to have a detectable GW counterpart. On the other hand, the prospects for VHE afterglow detection are not promising, with a rate of fewer than 10−1 events per year with our fiducial assumptions and with a CTA-like sensitivity.

For BHNS mergers, our results are less optimistic. In fact, only a small fraction of these events are expected to generate EM emission, primarily because tidal disruption of the NS occurs outside of the BH ISCO only under favourable conditions, such as low BH mass, high BH spin, and large NS deformability. Nevertheless, we anticipate a tenfold increase in the detection rate of BHNS mergers compared to current expectations for the O5 run, with a few KN detections per year and one GRB detection every few years, providing a unique opportunity to study these peculiar systems.

With an ETΔ detector alone, the above rates are reduced by a factor of 102 when considering only events with a sky localization uncertainty ΔΩ90% < 100 deg2. The fraction of reasonably well localized events increases by a factor of around four with a network of two L-shaped ET detectors, while a much more drastic improvement appears when considering ET in tandem with two CE detectors, bringing the fraction of well-localized events (for which a multi-wavelength follow-up can be realized with conceivable resources) up to more than 70% of the detected sources.

We examined the impact of variations in the NS EoS, BH spin distribution, and mass distributions on our detection rates. For NSNS systems, the primary source of errors in the detection rates is the uncertainty in the local merger rate, while variations in the mass distribution and EoS generally result in smaller deviations, except in certain cases involving GRB emission, due to the specific conditions assumed for the launch of relativistic jets. For BHNS systems, variations in the EoS and the BH spin distribution can lead to an increase in rates by more than an order of magnitude, surpassing the uncertainty associated with the local merger rate. These dependences demonstrate the potential of multi-messenger observations to constrain the astrophysical properties of compact objects and refine our understanding of their formation channels.

Prospects for the observation of EM counterparts in the era of ET have already been studied in the literature using different approaches, considering optical (Chen et al. 2021; Branchesi et al. 2023; Loffredo et al. 2025), high-energy (Ronchini et al. 2022; Hendriks et al. 2023), and very-high-energy (Banerjee et al. 2023) observations. In general, these studies focus on a specific band and the corresponding observational facilities, considering only counterparts from NSNS mergers. For BHNS, studies have considered observations of KN and afterglow (Zhu et al. 2021; Boersma & van Leeuwen 2022).

In this work, for the first time, we consider a population of both NSNS and BHNS with a consistent approach, simultaneously analysing the joint detection of KN, GRB afterglow, and GRB prompt emission across multiple EM bands. Additionally, we provide prospects for the VHE band of the afterglow assuming CTA, an estimate not yet present in the literature.

In conclusion, ET, particularly when operating as part of a network with other third-generation GW observatories, will open new frontiers in multi-messenger astronomy. It will allow for a comprehensive investigation of the origin and properties of compact binary mergers, advancing our knowledge of dense matter physics, jet dynamics, and heavy element nucleosynthesis. The next decade promises the arrival of transformative discoveries that will bridge GW and EM observations, providing a holistic view of some of the Universe’s most powerful events.

Data availability

The data produced in this work are publicly available on Zenodo through the link https://doi.org/10.5281/zenodo.15411012. The scripts and files to reproduce the main figures in the text are publicly available at https://github.com/acolombo140/ET_MM.

1

The plan for the observing run can be seen here: https://observing.docs.ligo.org/plan/

2

A detailed derivation of the correction factor for the local merger rate density is provided here: https://dcc.ligo.org/LIGO-P2400022/public

3

We assumed four BHNS events with a false alarm rate of less than one in 4 years, as reported in https://emfollow.docs.ligo.org/userguide/capabilities.html

4

The ET PSDs are publicly available at https://apps.et-gw.eu/tds/?content=3r=18213

5

The CE PSDs we used are publicly available at https://dcc.cosmicexplorer.org/CE-T2000017/public

6

The effect of this variation on the rates is included in the error bars in the right panel of Figure 4. The effect on the GRB properties is shown as an orange curve in the right panels of Figures 8 and 9.

7

The sensitivity curves are reported here: https://www.ctao.org/for-scientists/performance/

Acknowledgments

We acknowledge useful discussion with E. Loffredo, S. Ronchini, T. Fragos, P. Schmidt, F. Gabrielli. A.C. is supported by ERC grant GravitySirens 101163912. Funded by the European Union. Views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. M.M. is supported by the French government under the France 2030 investment plan, as part of the Initiative d’Excellence d’Aix-Marseille Université – A*MIDEX AMX-22-CEI-02. The work of F. I. is supported by the Swiss National Science Foundation, grant 200020_191957, by the SwissMap National Center for Competence in Research, and by a Miller Postdoctoral Fellowship. Part of the numerical computations made use of the Baobab cluster at the University of Geneva.

References

  1. Aasi, J., Abbott, B. P., Abbott, R., & the LIGO Scientific Collaboration 2015, Class. Quant. Grav., 32, 074001 [Google Scholar]
  2. Abac, A. G., Abbott, R., Abouelfettouh, I., et al. 2024, ApJ, 970, L34 [CrossRef] [Google Scholar]
  3. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 119, 161101 [Google Scholar]
  4. Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., et al. 2017b, ApJ, 848, L13 [CrossRef] [Google Scholar]
  5. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJ, 848, L12 [Google Scholar]
  6. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, Class. Quant. Grav., 34, 044001 [NASA ADS] [CrossRef] [Google Scholar]
  7. Abbott, B. P., LIGO Scientific Collaboration,& Virgo Collaboration 2019, Phys. Rev. X, 9, 011001 [NASA ADS] [Google Scholar]
  8. Abbott, B. P., LIGO Scientific Collaboration,& Virgo Collaboration 2020a, ApJ, 892, L3 [NASA ADS] [CrossRef] [Google Scholar]
  9. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020b, Liv. Rev. Relat., 23, 3 [NASA ADS] [Google Scholar]
  10. Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, ApJ, 896, L44 [Google Scholar]
  11. Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]
  12. Abbott, R., LIGO Scientific Collaboration,& Virgo Collaboration 2023, Phys. Rev. X, 13, 041039 [Google Scholar]
  13. Acernese, F., Agathos, M., Agatsuma, K., Aisa, D., et al. 2015, Class. Quant. Grav., 32, 024001 [NASA ADS] [CrossRef] [Google Scholar]
  14. Amati, L., O’Brien, P. T., Götz, D., et al. 2021, Exp. Astron., 52, 183 [NASA ADS] [CrossRef] [Google Scholar]
  15. Andreoni, I., Margutti, R., Salafia, O. S., et al. 2022, ApJS, 260, 18 [NASA ADS] [CrossRef] [Google Scholar]
  16. Aso, Y., Michimura, Y., Somiya, K., et al. 2013, Phys. Rev. D, 88, 043007 [NASA ADS] [CrossRef] [Google Scholar]
  17. Banerjee, B., Oganesyan, G., Branchesi, M., et al. 2023, A&A, 678, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Barbieri, C., Salafia, O. S., Perego, A., Colpi, M., & Ghirlanda, G. 2020, Eur. Phys. J. A, 56, 8 [NASA ADS] [CrossRef] [Google Scholar]
  19. Barbieri, C., Salafia, O. S., Colpi, M., Ghirlanda, G., & Perego, A. 2021, A&A, 654, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Bavera, S. S., Zevin, M., & Fragos, T. 2021, Res. Notes Am. Astron. Soc., 5, 127 [Google Scholar]
  22. Bavera, S. S., Fragos, T., Zapartas, E., et al. 2023, Nat. Astron., 7, 1090 [NASA ADS] [CrossRef] [Google Scholar]
  23. Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Biscoveanu, S., Landry, P., & Vitale, S. 2023, MNRAS, 518, 5298 [Google Scholar]
  25. Boersma, O. M., & van Leeuwen, J. 2022, A&A, 664, A160 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Branchesi, M., Maggiore, M., Alonso, D., et al. 2023, JCAP, 2023, 068 [CrossRef] [Google Scholar]
  27. Braun, R., Bonaldi, A., Bourke, T., Keane, E., & Wagg, J. 2019, arXiv e-prints [arXiv:1912.12699] [Google Scholar]
  28. Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, MNRAS, 508, 5028 [NASA ADS] [CrossRef] [Google Scholar]
  29. Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2011, ApJ, 740, 100 [NASA ADS] [CrossRef] [Google Scholar]
  30. Burns, E., Connaughton, V., Zhang, B.-B., et al. 2016, ApJ, 818, 110 [NASA ADS] [CrossRef] [Google Scholar]
  31. Chen, H.-Y., Cowperthwaite, P. S., Metzger, B. D., & Berger, E. 2021, ApJ, 908, L4 [NASA ADS] [CrossRef] [Google Scholar]
  32. Chen, S., Wang, L., Hayashi, K., et al. 2024, Phys. Rev. D, 110, 063016 [Google Scholar]
  33. Cherenkov Telescope Array Consortium, Acharya, B. S., Agudo, I., et al. 2019, Science with the Cherenkov Telescope Array (WORLD SCIENTIFIC) [Google Scholar]
  34. Colombo, A., Salafia, O. S., Gabrielli, F., et al. 2022, ApJ, 937, 79 [NASA ADS] [CrossRef] [Google Scholar]
  35. Colombo, A., Duqué, R., Salafia, O. S., et al. 2024, A&A, 686, A265 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Corsi, A., Lloyd-Ronning, N. M., Carbone, D., et al. 2019, arXiv e-prints [arXiv:1903.10589] [Google Scholar]
  37. Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556 [NASA ADS] [CrossRef] [Google Scholar]
  38. Dietrich, T., Samajdar, A., Khan, S., et al. 2019, Phys. Rev. D, 100, 044003 [NASA ADS] [CrossRef] [Google Scholar]
  39. Dobie, D., Murphy, T., Kaplan, D. L., et al. 2021, MNRAS, 505, 2647 [NASA ADS] [CrossRef] [Google Scholar]
  40. Duffell, P. C., Quataert, E., & MacFadyen, A. I. 2015, ApJ, 813, 64 [NASA ADS] [CrossRef] [Google Scholar]
  41. Dupletsa, U., Harms, J., Banerjee, B., et al. 2023, Astron. Comput., 42, 100671 [NASA ADS] [CrossRef] [Google Scholar]
  42. Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126 [NASA ADS] [CrossRef] [Google Scholar]
  43. Evans, M., Adhikari, R. X., Afle, C., et al. 2021, arXiv e-prints [arXiv:2109.09882] [Google Scholar]
  44. Evans, M., Corsi, A., Afle, C., et al. 2023, arXiv e-prints [arXiv:2306.13745] [Google Scholar]
  45. Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103 [NASA ADS] [CrossRef] [Google Scholar]
  46. Farrow, N., Zhu, X.-J., & Thrane, E. 2019, ApJ, 876, 18 [NASA ADS] [CrossRef] [Google Scholar]
  47. Fernández, R., & Metzger, B. D. 2013, MNRAS, 435, 502 [CrossRef] [Google Scholar]
  48. Fiore, F., Burderi, L., Lavagna, M., et al. 2020, SPIE Conf. Ser., 11444, 114441R [NASA ADS] [Google Scholar]
  49. Foucart, F. 2020, Front. Astron. Space Sci., 7, 46 [NASA ADS] [CrossRef] [Google Scholar]
  50. Foucart, F., Hinderer, T., & Nissanke, S. 2018, Phys. Rev. D, 98, 081501 [NASA ADS] [CrossRef] [Google Scholar]
  51. Foucart, F., Duez, M. D., Kidder, L. E., et al. 2019, Phys. Rev. D, 99, 103025 [NASA ADS] [CrossRef] [Google Scholar]
  52. Fragos, T., & McClintock, J. E. 2015, ApJ, 800, 17 [Google Scholar]
  53. Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91 [Google Scholar]
  54. Fuller, J., & Ma, L. 2019, ApJ, 881, L1 [Google Scholar]
  55. Ghirlanda, G., Salafia, O. S., Paragi, Z., et al. 2019, Science, 363, 968 [NASA ADS] [CrossRef] [Google Scholar]
  56. Ghirlanda, G., Nava, L., Salafia, O., et al. 2024, A&A, 689, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Gottlieb, O., & Nakar, E. 2022, MNRAS, 517, 1640 [NASA ADS] [CrossRef] [Google Scholar]
  58. Gupta, I., Afle, C., Arun, K. G., et al. 2024, Class. Quant. Grav., 41, 245001 [Google Scholar]
  59. Hallinan, G., Ravi, V., Weinreb, S., et al. 2019, Bull. Am. Astron. Soc., 51, 255 [Google Scholar]
  60. Hamidani, H., & Ioka, K. 2021, MNRAS, 500, 627 [Google Scholar]
  61. Hempel, M., Fischer, T., Schaffner-Bielich, J., & Liebendörfer, M. 2012, ApJ, 748, 70 [CrossRef] [Google Scholar]
  62. Hendriks, K., Yi, S.-X., & Nelemans, G. 2023, A&A, 672, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Iacovelli, F., Mancarella, M., Foffa, S., & Maggiore, M. 2022a, ApJ, 941, 208 [NASA ADS] [CrossRef] [Google Scholar]
  64. Iacovelli, F., Mancarella, M., Foffa, S., & Maggiore, M. 2022b, ApJS, 263, 2 [NASA ADS] [CrossRef] [Google Scholar]
  65. Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1 [Google Scholar]
  66. Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [Google Scholar]
  67. Just, O., Bauswein, A., Ardevol Pulpillo, R., Goriely, S., & Janka, H. T. 2015, MNRAS, 448, 541 [NASA ADS] [CrossRef] [Google Scholar]
  68. Kalogera, V., Sathyaprakash, B. S., Bailes, M., et al. 2021, arXiv e-prints [arXiv:2111.06990] [Google Scholar]
  69. Kawaguchi, K., Kyutoku, K., Nakano, H., et al. 2015, Phys. Rev. D, 92, 024014 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kawaguchi, K., Kyutoku, K., Shibata, M., & Tanaka, M. 2016, ApJ, 825, 52 [NASA ADS] [CrossRef] [Google Scholar]
  71. Kiziltan, B., Kottas, A., De Yoreo, M., & Thorsett, S. E. 2013, ApJ, 778, 66 [NASA ADS] [CrossRef] [Google Scholar]
  72. Krüger, C. J., & Foucart, F. 2020, Phys. Rev. D, 101, 103002 [CrossRef] [Google Scholar]
  73. Lattimer, J. M., & Schramm, D. N. 1974, ApJ, 192, L145 [NASA ADS] [CrossRef] [Google Scholar]
  74. Lazzati, D., & Perna, R. 2019, ApJ, 881, 89 [NASA ADS] [CrossRef] [Google Scholar]
  75. Lazzati, D., Perna, R., Morsony, B. J., et al. 2018, Phys. Rev. Lett., 120, 241103 [NASA ADS] [CrossRef] [Google Scholar]
  76. Levan, A. J., Gompertz, B. P., Salafia, O. S., et al. 2024, Nature, 626, 737 [NASA ADS] [CrossRef] [Google Scholar]
  77. Li, L.-X., & Paczyński, B. 1998, ApJ, 507, L59 [NASA ADS] [CrossRef] [Google Scholar]
  78. Loffredo, E., Hazra, N., Dupletsa, U., et al. 2025, A&A, 697, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415 [Google Scholar]
  80. Maggiore, M. 2007, Gravitational Waves. Vol. 1: Theory and Experiments (Oxford University Press) [Google Scholar]
  81. Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, JCAP, 2020, 050 [CrossRef] [Google Scholar]
  82. Mainieri, V., Anderson, R. I., Brinchmann, J., et al. 2024, arXiv e-prints [arXiv:2403.05398] [Google Scholar]
  83. Makhathini, S., Mooley, K. P., Brightman, M., et al. 2021, ApJ, 922, 154 [NASA ADS] [CrossRef] [Google Scholar]
  84. Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391 [NASA ADS] [CrossRef] [Google Scholar]
  85. Marconi, A., Abreu, M., Adibekyan, V., et al. 2022, SPIE Conf. Ser., 12184, 1218424 [NASA ADS] [Google Scholar]
  86. Mei, A., Banerjee, B., Oganesyan, G., et al. 2022, Nature, 612, 236 [NASA ADS] [CrossRef] [Google Scholar]
  87. Metzger, B. D. 2019, Liv. Rev. Relat., 23, 1 [Google Scholar]
  88. Miceli, D., & Nava, L. 2022, Galaxies, 10, 66 [NASA ADS] [CrossRef] [Google Scholar]
  89. Mochkovitch, R., Hernanz, M., Isern, J., & Martin, X. 1993, Nature, 361, 236 [NASA ADS] [CrossRef] [Google Scholar]
  90. Mooley, K. P., Frail, D. A., Dobie, D., et al. 2018, ApJ, 868, L11 [Google Scholar]
  91. Nandra, K., Barret, D., Barcons, X., et al. 2013, arXiv e-prints [arXiv:1306.2307] [Google Scholar]
  92. Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [Google Scholar]
  93. Özel, F., & Freire, P. 2016, ARA&A, 54, 401 [Google Scholar]
  94. Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918 [CrossRef] [Google Scholar]
  95. Özel, F., Psaltis, D., Narayan, R., & Santos Villarreal, A. 2012, ApJ, 757, 55 [Google Scholar]
  96. Pannarale, F., Berti, E., Kyutoku, K., Lackey, B. D., & Shibata, M. 2015, Phys. Rev. D, 92, 084050 [NASA ADS] [CrossRef] [Google Scholar]
  97. Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67 [Google Scholar]
  98. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Class. Quant. Grav., 27, 194002 [Google Scholar]
  100. Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Rastinejad, J. C., Gompertz, B. P., Levan, A. J., et al. 2022, Nature, 612, 223 [NASA ADS] [CrossRef] [Google Scholar]
  102. Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, Bull. Am. Astron. Soc., 51, 35 [Google Scholar]
  103. Román-Garza, J., Bavera, S. S., Fragos, T., et al. 2021a, ApJ, 912, L23 [CrossRef] [Google Scholar]
  104. Román-Garza, J., Bavera, S. S., Fragos, T., et al. 2021b, ApJ, 912, L23 [CrossRef] [Google Scholar]
  105. Ronchini, S., Branchesi, M., Oganesyan, G., et al. 2022, A&A, 665, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Safarzadeh, M., & Berger, E. 2019, ApJ, 878, L12 [Google Scholar]
  107. Salafia, O. S., & Ghirlanda, G. 2022, Galaxies, 10, 93 [NASA ADS] [CrossRef] [Google Scholar]
  108. Salafia, O. S., & Giacomazzo, B. 2021, A&A, 645, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Salafia, O. S., Ghirlanda, G., Ascenzi, S., & Ghisellini, G. 2019, A&A, 628, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Salafia, O. S., Barbieri, C., Ascenzi, S., & Toffano, M. 2020, A&A, 636, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Salafia, O. S., Ravasio, M. E., Yang, J., et al. 2022, ApJ, 931, L19 [NASA ADS] [CrossRef] [Google Scholar]
  112. Salafia, O. S., Ravasio, M. E., Ghirlanda, G., & Mandel, I. 2023, A&A, 680, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Schutz, B. F. 2011, Class. Quant. Grav., 28, 125023 [Google Scholar]
  114. Schwab, J., Podsiadlowski, P., & Rappaport, S. 2010, ApJ, 719, 722 [Google Scholar]
  115. Steiner, A. W., Lattimer, J. M., & Brown, E. F. 2013, ApJ, 765, L5 [Google Scholar]
  116. Troja, E., Fryer, C. L., O’Connor, B., et al. 2022, Nature, 612, 228 [NASA ADS] [CrossRef] [Google Scholar]
  117. Urrutia, G., De Colle, F., Murguia-Berthier, A., & Ramirez-Ruiz, E. 2021, MNRAS, 503, 4363 [CrossRef] [Google Scholar]
  118. Valentim, R., Rangel, E., & Horvath, J. E. 2011, MNRAS, 414, 1427 [NASA ADS] [CrossRef] [Google Scholar]
  119. Vallisneri, M. 2008, Phys. Rev. D, 77, 042001 [NASA ADS] [CrossRef] [Google Scholar]
  120. Veitch, J., Raymond, V., Farr, B., et al. 2015, Phys. Rev. D, 91, 042003 [NASA ADS] [CrossRef] [Google Scholar]
  121. Villar, V. A., Guillochon, J., Berger, E., et al. 2017, ApJ, 851, L21 [Google Scholar]
  122. Xing, Z., Bavera, S. S., Fragos, T., et al. 2024, A&A, 683, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  123. Xing, Z., Kalogera, V., Fragos, T., et al. 2025, ApJ, 989, 188 [Google Scholar]
  124. Yang, Y.-H., Troja, E., O’Connor, B., et al. 2024, Nature, 626, 742 [NASA ADS] [CrossRef] [Google Scholar]
  125. You, Z.-Q., Zhu, X., Liu, X., et al. 2025, Nat. Astron., 9, 552 [Google Scholar]
  126. Zevin, M., Spera, M., Berry, C. P. L., & Kalogera, V. 2020, ApJ, 899, L1 [Google Scholar]
  127. Zhu, J.-P., Wu, S., Yang, Y.-P., et al. 2021, ApJ, 917, 24 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Population model

Figure A.1 illustrates the mass distribution from Colombo et al. (2022) compared against observational data on the (M1, M2) plane. The plot also includes iso-contours of ejecta and accretion disc masses derived using the adopted fitting formulae (Krüger & Foucart 2020; Barbieri et al. 2021) and the equations of state (EoS) SFHo (left panel) and DD2 (right panel). These contours provide a visual representation of the lack of EM counterparts for events located in the upper right region of the plane and highlight the general trends of ejecta and disc mass distributions within the population. Additionally, the pink line indicates the condition for forming a HMNS remnant (Mrem > 1.2MTOV), which is also the threshold for launching a relativistic jet, together with the requirement of a non-negligible accretion disc (mdisk > 0).

thumbnail Fig. A.1.

Mass distribution of our NSNS population in the M1, M2 plane. The filled blue coloured regions contain 50%, 90%, and 99% of the binaries. The black dashed lines and the grey lines represent respectively the contours for the predicted dynamical ejecta and disc mass, assuming the SFHo EoS (left panel) and the DD2 EoS (right panel). The pink line indicates the condition for a HMNS remnant (Mrem > 1.2MTOV). Red stars and contours show the best fit and 90% credible regions for the known Galactic NSNS (Özel & Freire 2016; Farrow et al. 2019) systems that merge within a Hubble time. Yellow and aquamarine lines represent the 50% confidence regions for the component masses in GW190425 (Abbott et al. 2020a) and GW170817 (Abbott et al. 2019), both constructed using the publicly available low-spin-prior posterior samples.

thumbnail Fig. A.2.

Mass distribution of our BHNS population at redshift z = 0 (fiducial model in Broekgaarden et al. (2021)) in the MNS, MBH plane. The filled blue coloured regions contain 50%, 90%, and 99% of the binaries. The black dashed lines and the dark red lines represent respectively the contours for the predicted dynamical ejecta and disc mass, assuming the SFHo EoS (upper panel) and the DD2 EoS (lower panel). Violet and green lines represent the 50% and 90% confidence regions for the component masses in GW200115, GW200105 (Abbott et al. 2021) and GW230529, both constructed using the publicly available low-spin-prior posterior samples.

In Figure A.2, we display analogous information to that shown in Figure A.1, but for the BHNS population. The figure illustrates our fiducial mass distribution for BHNS binaries at redshift z = 0 (Broekgaarden et al. 2021) on the (MNS, MBH) plane. The shaded blue regions represent the areas containing 50%, 90%, and 99% of the binary systems. Iso-contours of ejecta and accretion disc masses are also included, computed using the fitting formulae from Krüger & Foucart (2020) and Kawaguchi et al. (2016) for two equations of state (EoS), SFHo (upper panel) and DD2 (right panel), under three different BH spin configurations: χBH = 0, χBH = 0.5, and χBH = 0.9. The figure highlights the general trends in the distribution of ejecta and disc masses within the population and emphasizes the strong dependence of ejecta production on the BH spin.

Appendix B: Detection rates as a function of the assumed EM detection threshold

Figures B.1 and B.2 present the same information as in Figure 6, but for GRB afterglow+GW events and for GRB prompt+GW events, respectively. In Figure B.1, the upper panel refers to the NSNS population, while the lower panel refers to the BHNS population. The red, green, and blue colours indicate the radio, optical, and X bands. The solid, dashed, and dotted lines correspond to all GW events, and to those with ΔΩ90% < 100 deg2 and ΔΩ90% < 10 deg2, respectively. The configurations ETΔ, ET2L, ETΔ+2CE, and ET2L+2CE are represented by the coloured, light grey, dark grey, and black lines, respectively.

In Figure B.2, all lines refer to the bolometric fluence. The left panel refers to the NSNS population and the right panel to the BHNS population.

thumbnail Fig. B.1.

Gamma-ray burst afterglow detection rate as a function of the EM detection limit threshold for our fiducial NSNS (upper panels) and BHNS (lower panels) populations. The red, green, and blue colours indicate the radio, optical, and X bands, respectively, assuming the ETΔ configuration. In each panel we also report in grey and black the ET2L and ETΔ+2CE configurations. The solid line indicates all the detectable binaries, the dashed and dotted lines the detectable binaries with ΔΩ90% < 100deg2 and the ones with ΔΩ90% < 10deg2, respectively.

thumbnail Fig. B.2.

Gamma-ray burst prompt detection rate as a function of the EM detection limit threshold for our fiducial NSNS (left panel) and BHNS (right panel) populations. We account for the Fermi/GBM duty cycle. In each panel we report in orange, grey and black the ETΔ, ET2L and ETΔ+2CE configurations, respectively. The solid line indicates all the detectable binaries, the dashed and dotted lines the detectable binaries with ΔΩ90% < 100deg2 and the ones with ΔΩ90% < 10deg2, respectively.

All Tables

Table 1.

Summary of the interferometer assumptions.

In the text
Table 2.

Detection limits and predicted detection rates for NSNS and BHNS assuming different detector network configurations.

In the text

All Figures

thumbnail Fig. 1.

Cumulative distribution of events above a threshold bolometric fluence. The solid line represents the distribution for Fermi/GBM observed SGRBs (with a duration below the customary 2 s threshold), constructed using the spectral information available in the online catalog. The coloured band shows the associated Poisson uncertainty. The dotted line shows our model distribution for comparison.
In the text
thumbnail Fig. 2.

Representation of two of our fiducial multi-messenger synthetic populations in a geocentric universe. Compact binaries that are detected by ETΔ are represented by yellow (NSNS – left semicircle) and orange (BHNS – right semicircle) dots. The distance of each dot from the Earth (blue circle) is proportional to the redshift of the corresponding compact binary. If the simulated merger produces a relativistic jet whose prompt or afterglow emission is detectable, according to the limits reported in Table 2, a cyan jet is plotted centered on the dot, with its axis inclined by the actual viewing angle with respect to the line of sight to the Earth. If a KN is also produced and if it is detectable, then a red butterfly shape is also plotted. The total number of binaries is representative of 5 years of ET operation.
In the text
thumbnail Fig. 3.

Cumulative multi-messenger detection rates as a function of redshift (luminosity distance) for our fiducial NSNS and BHNS population (SFHo EoS, non-spinning BHs), assuming the ET triangle 10 km configuration. Left panel: NSNS population. The light grey line (“All NSNS”) represents the intrinsic merger rate for the NSNS population, with the grey band showing its uncertainty due to that on the local merger rate. This uncertainty propagates as a constant relative error to all the other rates. The black (“GW ET”) line is the cumulative GW detection rate (events per year with network S/N ≥ 12). The blue (“Kilonova+GW”), red (“Afterglow+GW”), purple (“Afterglow VHE+GW”) and orange (“Prompt+GW”) lines are the cumulative detection rates for the joint detection of ET GW plus a KN (g band), a GRB afterglow (radio and VHE bands), or a GRB prompt (the orange and purple lines account, respectively, for the Fermi/GBM and CTA duty cycle and field of view). The dashed lines are the cumulative detection rates assuming only the binaries with ΔΩ90% < 100 deg2. The assumed thresholds or instruments sensitivity are shown in the legend. Right panel: Same as the left panel but for the BHNS population.
In the text
thumbnail Fig. 4.

Predicted 90th percentile of the cumulative multi-messenger detection rates for our NSNS and BHNS population model variations. Left panel: NSNS fiducial population, assuming ET triangle 10 km configuration and EM bands and detection limits reported in Figure 3. Different colours refer to different counterparts as in Figure 3. Different marker shapes indicate different adopted mass distributions (triangle: fiducial; square: uniform; circle: Gaussian). Filled markers indicate the SFHo EoS, while empty markers are for the DD2 EoS. The error bars indicate the uncertainty on the local merger rate and for the “Afterglow VHE+GW” channel also a variation on the median circum-burst density. Right panel: BHNS fiducial population, assuming ET triangle 10 km configuration and EM limits reported in Table 2. Different marker shapes indicate different adopted BH spin distributions (triangle: χBH = 0; square: uniform between 0 and 0.5). Filled markers indicate the SFHo EoS, while empty markers are for the DD2 EoS. The error bars indicate the uncertainty on the local merger rate. For GRB afterglow and prompt they also take into account a variation on the jet core half-opening angle (θj = 15°) and for the Afterglow VHE also a variation on the median circum-burst density.
In the text
thumbnail Fig. 5.

Distribution of the extent of the 90% sky localization region of GW-detected NSNS events (upper row) and BHNS events (lower row) that are capable of producing an EM counterpart, as a function of redshift. Different columns refer to ETΔ (first column), ET2L (second column), ETΔ+2CE (third column), and ET2L+2CE (four column). The solid line represents the median, while the coloured bands encompass the 50%, 90%, and 99% credible interval at each fixed redshift.
In the text
thumbnail Fig. 6.

Kilonova detection rate as a function of the EM detection limit threshold for our fiducial NSNS (upper panels) and BHNS (lower panels) populations. The blue, green, and red colours indicate the g, z, and J bands, respectively, assuming the ETΔ configuration. In each panel we also report in grey and black the ET2L and ETΔ+2CE configurations. The solid line indicates all the detectable binaries; the dashed and dotted lines show the detectable binaries with ΔΩ90% < 100 deg2 and the ones with ΔΩ90% < 10 deg2, respectively.
In the text
thumbnail Fig. 7.

Distribution of optical KN brightness as a function of time for ET-detectable (assuming ETΔ) events in our NSNS (left column) and BHNS (right column) populations. The shaded regions contain 50%, 90%, and 99% of the KN light curves in the g (484 nm) band. Coloured circles show extinction-corrected AT2017gfo data rescaled to the median distance of the considered populations (upper limits are marked with a downward arrow, data from Villar et al. 2017). The bottom row shows the result restricted to binaries with ΔΩ90% < 100 deg2.
In the text
thumbnail Fig. 8.

Distribution of brightness versus time for GRB afterglows associated with ET-detectable NSNS (left panel) and BHNS (right panel) mergers. Solid and dashed contours contain 50% and 90% of the peaks, respectively. Red, green, and blue colours refer to the radio (Fν at 1.4 × 109 Hz), optical (AB magnitude at 4.8 × 1014 Hz), and X-rays (νFν at 1 keV), respectively. The coloured circles are the observed data of GRB170817A (Makhathini et al. 2021) at the median distance of our population. The grey lines in the background are a random sample of optical light curves from the underlying population. We assume for both the populations the same jet half-opening angle, θj = 3.4 deg. In the BHNS panel, we also show with orange solid lines the 50% contour for radio peaks when assuming θj = 15 deg.
In the text
thumbnail Fig. 9.

Rest-frame SED peak photon energy Epeak versus the isotropic-equivalent energy, Eiso, for our NSNS (left panel) and BHNS (right panel) populations. The filled green coloured regions contain 50%, 90%, and 99% of the binaries both GRB prompt- and ET-detectable. The magenta lines contain 50%, 90%, and 99% (solid, dashed, and dotted, respectively) of the GRB prompt-detectable binaries. The black dashed line contains 90% of the ET-detectable binaries. The black dots with error bars represent an SGRB sample for comparison (Salafia et al. 2023). The orange dot is GRB170817A. We assumed for both populations the same jet half-opening angle θj = 3.4 deg. In the BHNS panel, we also show with orange solid lines the 50% contour for radio peaks when assuming θj = 15 deg.
In the text
thumbnail Fig. A.1.

Mass distribution of our NSNS population in the M1, M2 plane. The filled blue coloured regions contain 50%, 90%, and 99% of the binaries. The black dashed lines and the grey lines represent respectively the contours for the predicted dynamical ejecta and disc mass, assuming the SFHo EoS (left panel) and the DD2 EoS (right panel). The pink line indicates the condition for a HMNS remnant (Mrem > 1.2MTOV). Red stars and contours show the best fit and 90% credible regions for the known Galactic NSNS (Özel & Freire 2016; Farrow et al. 2019) systems that merge within a Hubble time. Yellow and aquamarine lines represent the 50% confidence regions for the component masses in GW190425 (Abbott et al. 2020a) and GW170817 (Abbott et al. 2019), both constructed using the publicly available low-spin-prior posterior samples.
In the text
thumbnail Fig. A.2.

Mass distribution of our BHNS population at redshift z = 0 (fiducial model in Broekgaarden et al. (2021)) in the MNS, MBH plane. The filled blue coloured regions contain 50%, 90%, and 99% of the binaries. The black dashed lines and the dark red lines represent respectively the contours for the predicted dynamical ejecta and disc mass, assuming the SFHo EoS (upper panel) and the DD2 EoS (lower panel). Violet and green lines represent the 50% and 90% confidence regions for the component masses in GW200115, GW200105 (Abbott et al. 2021) and GW230529, both constructed using the publicly available low-spin-prior posterior samples.
In the text
thumbnail Fig. B.1.

Gamma-ray burst afterglow detection rate as a function of the EM detection limit threshold for our fiducial NSNS (upper panels) and BHNS (lower panels) populations. The red, green, and blue colours indicate the radio, optical, and X bands, respectively, assuming the ETΔ configuration. In each panel we also report in grey and black the ET2L and ETΔ+2CE configurations. The solid line indicates all the detectable binaries, the dashed and dotted lines the detectable binaries with ΔΩ90% < 100deg2 and the ones with ΔΩ90% < 10deg2, respectively.
In the text
thumbnail Fig. B.2.

Gamma-ray burst prompt detection rate as a function of the EM detection limit threshold for our fiducial NSNS (left panel) and BHNS (right panel) populations. We account for the Fermi/GBM duty cycle. In each panel we report in orange, grey and black the ETΔ, ET2L and ETΔ+2CE configurations, respectively. The solid line indicates all the detectable binaries, the dashed and dotted lines the detectable binaries with ΔΩ90% < 100deg2 and the ones with ΔΩ90% < 10deg2, respectively.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.