© S. Ronchini et al. 2022

1. Introduction

The first detection of gravitational waves (GW) from the merger of two neutron stars, GW 170817, by the Advanced LIGO (The LIGO Scientific Collaboration 2015) and Virgo (The Virgo Collaboration 2015) interferometers and the observed electromagnetic (EM) counterparts represented a milestone for multimessenger astronomy, in particular for probing the nature of EM transients that originated from BNS mergers (Abbott et al. 2017a,b,c). The extensive study of the EM emission allowed us to obtain an exhaustive picture of the outflow properties of both the relativistic structured jet (Troja et al. 2017; Hallinan et al. 2017; Mooley et al. 2018; Ghirlanda et al. 2019) and the isotropic subrelativistic kilonova outflow (Pian et al. 2017; Smartt et al. 2017; Kasen et al. 2017). Considering the distance range expected to be reached by LIGO, Virgo, and KAGRA in the next runs (Abbott et al. 2020) for the detection of BNS mergers and the astrophysical rate of this type of system R_BNS = [10 − 1700] Gpc⁻³ yr⁻¹ (The LIGO Scientific Collaboration 2021a) derived in the first three runs of observations by the Advanced GW detectors, the detection rate is expected to range from tens to hundreds per year in the next observational runs. Therefore, taking into account that only a fraction of them produce a jet and that only a fraction of them point along the line of sight, the probability of joint GW and SGRB detection with current instruments is quite limited (Metzger & Berger 2012; Chen & Holz 2013; Clark et al. 2015; Patricelli et al. 2016, 2022; Ghirlanda et al. 2016; Beniamini et al. 2019; Saleem et al. 2018a,b; Saleem 2020; Perna et al. 2022; Colombo et al. 2022). The advent of third-generation (3G) GW observatories, such as the Einstein Telescope (ET; Punturo et al. 2010) and Cosmic Explorer (CE; Abbott et al. 2017d), will be a substantial leap forward compared to the current instruments for the cosmological census of GW sources. ET is an underground observatory, designed as three low and three high frequency interferometers nested in a triangular shape forming three 10 km arms. CE, instead, is an L-shaped surface interferometer 40 km in length. ET and CE are expected to detect ∼10⁵ BNS per year, reaching redshifts well above the star formation peak (Maggiore et al. 2020, Evans et al. 2021); for optimally located and oriented systems, ET will be able to detect BNS up to z ∼ 5, while CE will do so up to z ∼ 10. For comparison, the current GW detectors are expected to reach z ∼ 0.2, even considering new upgrades planned for the fifth run of observations in 2025 (Abbott et al. 2020). Furthermore, accessing low frequencies makes ET capable of detecting BNSs well before the merger and follow their inspiral not only for minutes but also for several hours in the case of closeby events (Chan et al. 2018; Grimm & Harms 2020; Nitz & Dal Canton 2021; Li et al. 2022). This enables us to use the Earth rotation to determine the sky localization even when ET is operating as a single detector (Chan et al. 2018; Grimm & Harms 2020; Nitz & Dal Canton 2021; Li et al. 2022). These authors found that a fraction of a few percents of BNS detected by ET will have a sky localization uncertainity ΔΩ < 100 deg². This fraction increases to a few tens of percents when ET is operating in network with CE. Considering a network of three 3G detectors located in Europe, USA, and Australia a few tens of percents of detected BNS will be localized within 10 deg² and almost the totality within 100 deg² (Li et al. 2022, Borhanian & Sathyaprakash 2022, Mills et al. 2018). A detailed study of performance metrics of several configurations of network including CE, ET, and the second generation detectors and their upgrade, are given in Borhanian & Sathyaprakash (2022). Another proposed GW observatory to be mentioned is the Neutron star Extreme Matter Observatory (NEMO, Ackley et al. 2020), a 4-km detector whose sensitivity at high frequencies will be relevant for the investigation of neutron star physics and for multimessenger purposes (Sarin & Lasky 2022).

Wide field-of-view (FOV) γ-ray and X-ray observatories will become more and more crucial to work in synergy with GW detectors which will make it possible to reach even larger distances. They represent a unique way to detect high-z counterparts. Being intrinsically fainter than the high energy emission, the kilonova emission will be hardly detected at distances larger than a redshift of 0.2–0.3 even with next generation optical wide FOV observatories, such as the Vera Rubin Observatory (see e.g., Maggiore et al. 2020). Larger distances can be reached with extremely large telescopes, such as ELT, but only arcmin-arcsec localized GW sources can be effectively pointed by these telescopes. The next generation of GW detectors observing with sensitive wide FOV high-energy satellites will enable (i) to build large samples of joint detections, (ii) to have joint detections at high redshift, and (iii) to observe, in the local Universe, counterparts from larger viewing angles with respect to the jet axis which will enable precise parameter estimation of the source progenitor properties. This will allow us to deeply investigate the connection between BNS mergers and short GRBs, to unveil the jet structure (Hayes et al. 2020; Beniamini et al. 2019; Biscoveanu et al. 2020), to estimate how many BNS mergers produce jets (Salafia et al. 2022), how the production of the jet is connected to the progenitor properties, and what is the relation between merger remnants and the properties of the high-energy emission. Accessing a large sample of joint detections at high-redshift will enable us to evaluate cosmological parameters, including the Hubble constant and the dark energy equation of state (Zhao & Wen 2018; de Souza et al. 2022; Yu et al. 2021; Chen et al. 2021), and test modified gravity at cosmological distances by comparing the GW luminosity distance and the one derived from the electromagnetic side (Belgacem et al. 2019; Mancarella et al. 2022). Identifying the BNS host galaxies and having a redshift will be of primary importance for these science cases. γ-ray and X-ray telescopes can potentially give precise (arcmin-arcsec) sky localization to drive the ground-based follow-up.

To maximize the scientific return of the next generation GW detectors in the multimessenger context it is of paramount importance to evaluate the expected joint detections, to identify instrument scientific requirements and optimal observation strategy. In this work, we focus on the synergy of 3G GW detectors with high-energy (γ-ray and X-ray) telescopes. The novelty of our study consists in the development of a theoretical framework that starting from an astrophysically-motivated population of BNSs and the modeling of the associated high-energy EM signal is able to reproduce the statistical properties of current SGRB catalogs. This approach, based on calibrating the BNS population and the jet properties using observed SGRBs, allows us to be less dependent on the uncertainties related to the BNS population synthesis models about the local merger rate and its redshift distribution. On the GW side, we evaluate the detection efficiency of 3G GW detectors, as well as their capabilities in localizing the source. On the EM side, we provide a comprehensive overview of the detectability of prompt and afterglow emission associated with the SGRB.

The population synthesis model for the cosmic BNS merger rate density, the GW detector configurations and the modeling of the GW signal, the modeling and calibration of the prompt and afterglow emission are described in Sect. 2. In Sect. 3 we show the predictions for joint detection of GW along with prompt and afterglow emission exploring several observational strategies. We consider ET as single observatory and in different networks of 3G GW detectors observing together with future high-energy facilities, which can operate in survey and pointing modes. We conclude highlighting the perspectives for joint GW and high-energy detections, discussing the optimal instruments and observational strategies to completely exploit the synergy between 3G GW detectors and γ- and X-ray satellites, and the impact on the multimessenger science case.

2. Modeling and data analysis

2.1. The BNS population model

The population of BNS mergers is produced as detailed in Santoliquido et al. (2021). The catalog of isolated compact binaries is generated with MOBSE (Mapelli et al. 2017; Giacobbo et al. 2018; Giacobbo & Mapelli 2018), which is a population synthesis-code taking into account the metallicity dependence of mass-loss rate of massive stars, the uncertainty related to the common envelope evolution, natal kicks and mass transfer efficiency. The BNS merger rate density as a function of redshift is obtained with the COSMOℛATE code (Santoliquido et al. 2020). Here, COSMOℛATE adopts the star formation rate density and average metallicity evolution of the Universe from Madau & Fragos (2017). We assume a metallicity spread σ_Z = 0.3. Redshift evolution and overall normalization of BNS merger density depend on several assumptions regarding the common envelope prescription (parametrized by the common envelope ejection efficiency α_CE), the natal kick model, the supernova explosion, the mass transfer via accretion, and the shape of the initial mass function. Here, we describe electron-capture supernovae as in Giacobbo & Mapelli (2019) and assume the delayed supernova model (Fryer et al. 2012) to decide whether a core-collapse supernova produces a black hole or a neutron star. When a neutron star forms from either a core-collapse or an electron-capture supernova, we randomly draw its mass according to an uniform distribution between 1 and 2.5 M_⊙. The assumption of a uniform and broad mass distribution is consistent with the last results from GW observations (The LIGO Scientific Collaboration 2021b). In the appendix we discuss how using a mass distribution sharply peaked around 1.33 M_⊙ as observed in Galactic double neutron-star systems (Özel & Freire 2016) can impact our results. We generate the natal kicks according to v_kick ∝ m_ej/m_rem v_H05, where m_ej is the mass of the ejecta, m_rem is the mass of the compact remnant (neutron star or black hole) and v_H05 is a random number extracted from a Maxwellian distribution with one-dimensional root-mean square σ = 265 km s⁻¹ (Hobbs et al. 2005). This formalism ensures low kicks for most electron-capture and ultra-stripped supernovae (Giacobbo & Mapelli 2020). We assume a value of α_CE = 3 for common envelope, we calculate the concentration parameter λ as in Claeys et al. (2014) and the mass transfer formalism from Hurley et al. (2002). We draw the mass of the progenitor primary star from a Kroupa mass function (Kroupa 2001). For the initial mass ratios, orbital periods and eccentricities we use the distributions inferred by Sana et al. (2012). With such combination of assumptions the entire population consists of ∼9 × 10⁵ BNS mergers per year (observed at Earth by an ideal instrument), and the local BNS merger rate (R_BNS = 365 Gpc⁻³ yr⁻¹) is consistent with the most updated constraints provided by GW observations during the first, second and third run of observations by the Advanced LIGO and Virgo detectors (R_BNS ∈ [10 − 1700] Gpc⁻³ yr⁻¹, The LIGO Scientific Collaboration 2021b; Abbott et al. 2021a,b). Moreover, the local rate assumed here is perfectly consistent with the collection of local BNS merger rates given by Mandel & Broekgaarden (2022) and including both observed (from GW, kilonovae, SGRBs, and Galactic pulsar binaries’ observations) and theoretically predicted rates. The BNS systems are isotropically distributed in the sky and randomly oriented.

2.2. GW detector configuration, GW signal modeling and parameter estimation

In this work, we consider three GW detector configurations: (1) ET as a single detector located in Sardinia (Italy), (2) ET in a network with CE (40 km arms) located in the LIGO-Livingston site, (3) ET in a network with 2CE of 40 km arms, with the second CE located in Australia. We assume the full sensitivity configuration for ET (referred as ET-D configuration, Hild et al. 2011). For CE we use the sensitivity of the 40 km arms detector given by Evans et al. (2021). For CE as well as for each of the three combinations of high and low frequency interferometers of ET we assume a duty cycle of 0.85. Starting from the BNS population described in Sect. 2.1, we inject GW signals constructed using a post-Newtonian formalism, in particular assuming a TaylorF2 waveform (Buonanno et al. 2009). Since the spin is expected to be small in binaries of neutron stars that merge within a Hubble time (Burgay et al. 2003), the effects due to spins are neglected. The signal-to-noise ratio (S/N) of each inspiral is then computed through a matched-filter. A network S/N threshold equal to 8 is used to select GW detections. The parameter estimation is obtained using the GWFish code (Harms et al. 2022) and computing the elements of the inverse of the Fisher matrix (Grimm & Harms 2020). The code is publicly available at this repository¹. In GWFish the waveforms, modeled in frequency domain, are implemented from scratch and the code includes an interface with LALSimulation². Detector networks consider Earth’s rotation, which has an impact on the localization capabilities. The most important computational aspects (waveform derivatives and Fisher matrix inversion) are carried out using an hybrid approach for derivatives (analytic + numeric) and a combination of matrix normalization and singular value decomposition for matrix inversion (which takes care of the quasi singularity of the Fisher matrices). Such an approach approximates the true likelihood with a Gaussian profile, which is valid in the limit of a high information content in the data. In this regime the likelihood is highly peaked and the role of priors is negligible. Hereafter, the uncertainty on the sky localization ΔΩ is given as 90% credible region, while the uncertainty on other GW parameters, such as θ_v and luminosity distance, is given at 1σ confidence level.

2.3. Prompt emission modeling

For the estimation of the high-energy EM emission from the BNS merger we assume that only a fraction f_j of the mergers produces a jet which is able to penetrate the post-merger ejecta and to successfully break-out. In principle, f_j could depend on the properties of the binary system, but here we do not assume any dependency. The EM counterpart is evaluated in the assumption of a universal jet structure (Lipunov et al. 2001; Salafia et al. 2015, 2019, 2020), which is the same of GRB 170817A. In particular, we adopt:

$$\begin{array}{c}\hfill ϵ(\theta )=\frac{\mathrm{d}E}{\mathrm{d}\mathrm{\Omega }}\propto \frac{1}{1+{(\theta /{\theta }_c)}^{s_{ϵ}}}\end{array}$$(1)

and

$$\begin{array}{c}\hfill \mathrm{\Gamma }(\theta )=1+\frac{{\mathrm{\Gamma }}_0-1}{1+{(\theta /{\theta }_c)}^{s_{\mathrm{\Gamma }}}},\end{array}$$(2)

where ϵ(θ) is the angular structure of the local emissivity and Γ(θ) the bulk Lorentz factor profile (θ is the polar angle measured from the jet symmetry axis). The choice of θ_c and Γ₀ influences the number of prompt emission detections, while the off-core specific shape of the jet has more impact on the brightness of the afterglow component and the prompt tail. From the modeling of radio afterglow of GRB 170817A (Ghirlanda et al. 2019), the best fit parameters of the jet structure are the following:

$$\begin{array}{c}\hfill s_{ϵ}=5.5_{-1.4}^{+1.3},s_{\mathrm{\Gamma }}=3.5_{-1.7}^{+2.1},{\theta }_c=(3.4\pm 1.0)\text{deg}.\end{array}$$(3)

For our modeling, we fix for simplicity s_ϵ = s_Γ = s and we adopt:

$$\begin{array}{c}\hfill s=4,{\theta }_c=3.4\text{deg}.\end{array}$$(4)

Moreover, we fix Γ₀ = 500, which is in agreement with usual values inferred for GRBs (e.g., Ghirlanda et al. 2018).

These parameters are the same for the modeling of both prompt and afterglow emission. A broader structure with s = 2 and same θ_c is also tested (Rossi et al. 2002). Hereafter we refer to Stru1 for the case s = 4 and θ_c = 3.4 deg and Stru2 for the case s = 2 and θ_c = 3.4 deg.

The dependency of our results on the choice of the structure parameters is discussed later in the paper. The aperture angle assumed here is consistent with the values found in SGRBs with well constrained jet break (typical values of θ_c are in the interval ∼1 ° −8°, see Fong et al. 2014, Jin et al. 2018 and references within). Similar values are found in hydrodynamical simulations of jet propagation in NS mergers, for example Lazzati et al. (2018) predicts a jet aperture angle of ∼5° (see also Geng et al. 2019; Wu & MacFadyen 2018). Also Lamb et al. (2022) find a rather small jet core (θ_c ∼ 2.2°) and a jet profile compatible with Eqs. (1) and (2), with best fit values s_ϵ ∼ 3.1 and s_Γ ∼ 1.8.

We define E_t the beaming corrected energy budget of the jet, which is connected to the isotropic energy measured by an observer aligned to the jet (on-axis) in the following way:

$$\begin{array}{c}\hfill E_{\mathrm{iso}}=\frac{4\pi }{\mathrm{\Delta }\omega }{E}_t,\end{array}$$(5)

where Δω ∼ 4π(1 − cos(θ_c)). E_t is related to ϵ(θ) through the following relation:

$$\begin{array}{c}\hfill E_t={\displaystyle {\int }_{\mathrm{\Omega }}ϵ(\theta )\mathrm{d}\mathrm{\Omega }.}\end{array}$$

We note that E_t does not represent the energy content of the jet, but rather the energy which is converted into radiation through dissipation mechanisms (such as internal shocks or magnetic reconnection). We extract the value of E_t from the following probability distribution:

$$\begin{array}{c}\hfill P(E_t/E^{\ast })\propto {(E_t/E^{\ast })}^{-{\lambda }_E}{e}^{-E^{\ast }/E_t},\end{array}$$(6)

which is a power law with a low-energy exponential cut-off. The prompt emission spectrum is assumed to be angle-independent in the jet comoving frame and follows a smoothly broken power law, with low- and high-energy photon indices of 2/3 and 2.5, respectively. These are typical values measured for SGRBs (Nava et al. 2011). In the frame of an on-axis observer, we assume that the spectral peak energy has a log-normal distribution, that is to say:

$$\begin{array}{c}\hfill P(log(E_p))\propto exp(-\frac{{(log(E_p)-log({\mu }_E))}^2}{{\sigma }_E^2}),\end{array}$$(7)

where μ_E and σ_E are the mean value and the standard deviation of the distribution, respectively. The assumed distribution function corresponds to that adopted by Ghirlanda et al. (2016). Such distribution in the observer frame is equivalent to assume a monochromatic distribution of $E_p^{\prime }$ in the jet comoving frame and a log-normal distribution of Γ₀ (indeed $E_p={\mathrm{\Gamma }}_0{E}_p^{\prime }$). For the computation of the burst duration, we compute the time t₉₀ such that the energy emitted from 0 to t₉₀ is equal to the 90% of the total energy released by the GRB. The light curve is approximated with a linear rising phase from t = 0 to the peak time t = t_p, followed by a decreasing phase which evolves according to high latitude emission (HLE). Therefore we assume that a single pulse dominates the SGRB emission. For t_p we consider a log-normal distribution as well, as measured in the rest frame of the burst:

$$\begin{array}{c}\hfill P(log(t_p))\propto exp(-\frac{{(log(t_p)-log({\mu }_{\tau }))}^2}{{\sigma }_{\tau }^2}).\end{array}$$(8)

The angle-dependent prompt emission light curve is computed adopting the approach of Ascenzi et al. (2020). Specifically, we work in the approximation of an infinitesimal duration emission at a fixed radius R₀, namely the burst duration in the rest frame is negligible compared to the time needed by the jet to propagate up to R₀. The jet moves with a bulk Lorentz factor Γ₀. As shown by Oganesyan et al. (2020), the presence of a jet structure like the one adopted here produces an X-ray light curve with a steep decay followed by a shallower phase, known as plateau. The timescale ${\tau }_p=R_0/2c{\mathrm{\Gamma }}_0^2$, which corresponds to the typical prompt emission timescale, determines the duration of the plateau and larger values of τ_p produce a longer plateau. In order to find a fiducial combination of R₀ and Γ₀, we considered all the SGRBs from the Swift-XRT catalog with a plateau phase. We verified that the end time of the plateau typically spans in the range 10² − 10³ s. Having fixed Γ₀ = 500, we therefore adopted R₀ = 5 × 10¹⁴ cm in order to obtain a typical plateau duration around few hundred seconds. With this approach we are able to predict the prompt emission output, as well as the contribution at later time given by the jet wings due to HLE effects, which can be relevant for the detectability of the SGRB afterglow, especially in the X-rays and for off-axis observers.

2.4. Calibration of the prompt emission model

As illustrated in the previous section, our prompt emission model depends on seven free parameters which determine the observational properties of our SGRB sample. Here, we show how our prompt emission model is calibrated such that the BNS population is able to reproduce the number of SGRBs detected by Fermi-GBM and their statistical properties reported in the Fermi-GBM catalog³.

In particular, the choice of our fiducial parameters is done comparing our simulated observables with the distributions of observed peak energy, duration, peak flux and fluence reported by Ghirlanda et al. (2016). This calibration sample is obtained by selecting SGRBs in the Fermi catalog with peak flux F_p > 5 ph cm⁻² s⁻¹ to ensure the completeness of the sample. For the parameter estimation we use Goodman & Weare’s Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler⁴. The likelihood function is $log\mathcal{L}=log(P_C(N_{\mathrm{pre}}|N_{\det }))+\sum _{i=1}^{4}log(P_{\mathrm{KS},i})$, where P_C(N_pre|N_det) is the Poissonian probability that Fermi-GBM detects N_pre SGRBs per year given that the average observed number is N_det, while P_KS, i is, for each observable i, the Kolmogorov–Smirnov probability that the predicted and observed distributions come from the same population. The adopted likelihood is valid in the approximation that each observable is independent from the others.

The detection rate of SGRBs R_SGRB, det can be related to the BNS rate R_BNS as:

$$\begin{array}{c}\hfill R_{\mathrm{SGRB},\det }={\displaystyle \int }\left[{\frac{\mathrm{d}{N}_{\mathrm{SGRB}}}{\mathrm{d}z\mathrm{d}t}|}_{\det }\right]\mathrm{d}z,\end{array}$$(9)

where

$$\begin{array}{c}\hfill {\frac{\mathrm{d}{N}_{\mathrm{SGRB}}}{\mathrm{d}z\mathrm{d}t}|}_{\det }=⟨{\mathrm{\Omega }}_{\mathrm{obs}}⟩\times {\varphi }_{\det }(z)\times f_j\frac{\mathrm{d}{N}_{\mathrm{BNS}}}{\mathrm{d}z\mathrm{d}t}.\end{array}$$(10)

We note that ⟨Ω_obs⟩ is the time-averaged sky coverage of the instrument and f_j is the fraction of BNS able to form a jet. The quantity $\frac{\mathrm{d}{N}_{\mathrm{BNS}}}{\mathrm{d}z\mathrm{d}t}$ is the number of BNS mergers observed at Earth per unit time and redshift and it is provided by the population synthesis model. ϕ_det(z) is the variable that combines the detection efficiency of the instrument (which takes into account the k-correction) and the beaming of the SGRB emission. Namely, ϕ_det(z) represents the probability that a SGRB at redshift z is detectable by a given instrument at Earth, that is to say:

$$\begin{array}{c}\hfill {\varphi }_{\det }(z)=\frac{N_{\mathrm{SGRB}}(F>F_{\lim },z)}{N_{\mathrm{SGRB}}(z)},\end{array}$$

where $N_{\mathrm{SGRB}}(z)={\displaystyle \frac{\mathrm{d}{N}_{\mathrm{SGRB}}}{\mathrm{d}z\mathrm{d}t}}$, $N_{\mathrm{SGRB}}(F>F_{\lim },z)={\displaystyle \frac{\mathrm{d}{N}_{\mathrm{SGRB}}}{\mathrm{d}z\mathrm{d}t}}(F>F_{\lim })$ and F_lim is the limiting flux in the band of the instrument. In turn, N_SGRB(F > F_lim, z) depends on the luminosity function of SGRBs, their average spectral properties, as well as on the jet structure.

The inclusion of f_j as free parameter of our model makes our estimates of absolute numbers of detections independent on the overall normalization of the BNS population model which is still subject to large uncertainties (Santoliquido et al. 2021). Indeed the likelihood function is written in such a way that, whatever the BNS merger rate normalization is, the value of f_j is optimized to reproduce the current average rate of Fermi-GBM detections. Namely, in Eq. (10) the quantity $\frac{\mathrm{d}{N}_{\mathrm{BNS}}}{\mathrm{d}z\mathrm{d}t}$ is uniquely determined by the BNS population model, ϕ_det(z) depends on the prompt emission model and therefore the MCMC converges to that value of f_j that reproduces R_SGRB, det, which is known.

The MCMC sampling has been performed for two different jet structures Stru1 and Stru2. For both structures, we obtain the posterior distribution of the model parameters. The confidence intervals of the parameters are reported in Table 1, while the corner plots of the posterior distributions are reported in Figs. B.5 and B.6 for Stru1 and Stru2, respectively. The modification of the off-core slope of the jet structure does not influence the best fit parameters, except for f_j, which tends to be smaller for s = 2. Such result is justified by the fact that a broader structure implies a larger detectability of the prompt emission at larger viewing angles, which corresponds to a larger number of off-axis detections. Therefore, for a given BNS population and an average SGRB detection rate, a broader structure requires a smaller fraction of BNS to be able to produce SGRB. The confidence interval we obtain for f_j is not particularly tight and it is compatible with other works which combine EM observations of SGRBs and BNS rates from GW observations (Ghirlanda et al. 2019; Salafia et al. 2022; Sarin et al. 2022; Beniamini et al. 2019; Lamb & Kobayashi 2017).

2.5. Forward shock modeling

The forward shock (FS) light curves are produced using the Python package afterglowpy (Ryan et al. 2020). For the jet kinetic structure we assume the same profile adopted for the structure of the comoving frame emissivity reported in the prompt emission modeling (see Eq. (1)). This implies that we are assuming that the conversion efficiency from kinetic energy to radiated energy has no angular dependence. The detectability of X-ray emission is evaluated for both Stru1 and Stru2. Afterglowpy includes the parameter θ_w which is the angular extension of the jet wings. For Stru1 we assume θ_w = 15 deg, while for Stru2 we assume θ_w = 30 deg. The impact of the choice of θ_w on the detectability of X-ray emission is discussed later in the text.

The model depends also on the micro-physical parameters n₀, p, ϵ_e and ϵ_B which are the ISM density, the slope of the electron energy distribution, the fraction of energy carried by electrons and the fraction of energy carried by the magnetic field, respectively. In the following, we perform our simulation using two setups of parameters: 1) n₀ = 2.5 × 10⁻⁴ cm⁻³, p = 2.2, ϵ_e = 0.1 and ϵ_B = 1.3 × 10⁻⁴, which are the best fit values obtained from the multiwavelength (X-ray, optical, radio) modeling of the X-ray afterglow of GRB 170817A (Ghirlanda et al. 2019). These values are also consistent with the confidence intervals reported by Wu & MacFadyen (2019). We call this configuration FS-GW17. 2) For the ISM density we take the fiducial interval derived by Fong et al. (2015)n₀ ∈ (3 − 15)×10⁻³ cm⁻³. The other parameters are fixed to p = 2.2, ϵ_e = 0.1 and ϵ_B ∈ [0.01 − 0.1]. n₀ and ϵ_B are uniformly extracted from the confidence intervals reported above. This configuration is more representative for SGRB population with respect to FS-GW17, and we call it FS-SGRB.

If we call η the fraction of kinetic energy which is transformed into radiation, then $E_{\mathrm{rad}}/E_{\mathrm{kin}}=(1-\eta )/\eta =\stackrel{\sim }{\eta }$ and we assume $\stackrel{\sim }{\eta }$ randomly distributed in the interval [0.01 − 0.1]. Figure 1 shows a collection of X-ray light curves produced with our model at different viewing angles, including both HLE and FS (configuration FS-SGRB). The curves are obtained with a random extraction from the probability distributions of the prompt parameters described in Sect. 2.3 and the parameter of each distribution is extracted from the Monte Carlo posterior described in Sect. 2.4.

Fig. 1. X-ray light curves of SGRBs with viewing angle 0 ° < θv < 15° predicted by our model. The lines are produced with a random extraction from the posterior distribution of prompt parameters. In the computation we assume Stru1 and we include both HLE and FS (configuration FS-SGRB). The color of each line indicates the viewing angle. The X-ray luminosity is computed at 5 keV.

3. Results

3.1. Joint detection of GWs and the prompt emission

For the detection of the prompt emission from our SGRBs population, we consider the following γ-ray instruments: the Gamma-ray Burst Monitor on board of Fermi (Fermi-GBM, Meegan et al. 2009), the Swift-Burst Alert Telescope (Swift-BAT, Barthelmy et al. 2005), the Gravitational wave high-energy Electromagnetic Counterpart All-sky Monitor (GECAM, Li et al. 2020), the ECLAIRs telescope and the Gamma-Ray burst Monitor (GRM) on board of the Space-based multi-band astronomical Variable Objects Monitor (SVOM, Cordier et al. 2015; Götz & SVOM Collaboration 2012), the X/Gamma-ray Imaging Spectrometer (XGIS) on board of the Transient High-Energy Sky and Early Universe Surveyor (THESEUS, Amati et al. 2018, 2021; Ciolfi et al. 2021; Rosati et al. 2021), the Gamma-ray Transient Monitor (GTM) on board of the Transient Astrophysical Probe (TAP, Camp et al. 2019), and the High Energy Rapid Modular Ensemble of Satellites (HERMES, Fiore et al. 2020). The first three are currently operating missions; Swift and Fermi have been observing for more than 10 years, and GECAM was launched in 2020. SVOM is expected to be launched in 2023. Nominally they are not expected to be operative in mid ’30s and ’40s when ET and CE will start observations, but we take them as reference instruments. THESEUS, TAP, and HERMES are mission concepts. The HERMES-Pathfinder project consisting in the deployment of six cube-satellites is expected to be operative in the next few years (Fiore et al. 2021). In this work we consider the final HERMES as a full constellation of cube-satellites. In order to establish whether the prompt emission is detected, we compute the peak flux during 1 second exposure and we compare this value with the detection threshold of the instrument. The viewing angle of each injection is assumed to be exactly the same used for the evaluation of the GW signal, and thus extracted from the randomly oriented BNS population (see Sect. 2.1).

In order to estimate the number of joint detections, we simulate the high-energy emission (keV-MeV) according to the model described in Sect. 2.3. The results relative to one year of joint GW+γ-ray detections are reported in Tables 2 and 3 for Stru1 and Stru2, respectively, where we considered ET and ET+CE as GW networks. The uncertainty intervals are computed simulating 1000 years of observations and for each year the prompt emission parameters are randomly extracted from the Monte Carlo posterior distribution. In the computation of the absolute numbers, we took into account the FOV of each instrument, assuming that the probability of detection is FOV/4π, which implicitly assumes an optimistic duty cycle of 100%. In the case of Fermi-GBM, Swift-BAT and SVOM, we adopt a more realistic value of sky-averaged FOV, which takes into account Sun and Moon occultations (Burns et al. 2016). In the case of THESEUS-XGIS the FOV we assumed (∼2 sr) is relative to the energy band 2–150 keV. Since the FOV in the energy band above 150 keV is 2π sr, the reported number of joint detections is a conservative estimate. Our results show that, already with ET alone, ∼60 − 70% of all the SGRBs will have a detectable GW counterpart, and this fraction approaches 100% if we consider ET operating with CE. The less steep structure profile of Stru2 increases the number of joint detections for the majority of the satellites. Depending on the properties of each satellite, Tables 2 and 3 show instruments giving a few tens of detections per year and other hundreds of detections per year. While instruments such as GECAM are optimal for statistical studies by giving a large number of detections, instruments with a smaller number of detections but able to localize the source, such as THESEUS-XGIS, are crucial for science cases requiring the identification of the host galaxy, the knowledge of the source redshift, and the complete multiwavelength characterization of the source emission. Furthermore, the sensitivity of the instruments considered for the detection of the prompt emission maximizes in different energy bands. While instruments such as GECAM are appropriate for detecting GRBs with a harder spectrum, GRBs less energetic (intrinsically or because viewed off-axis) will peak at lower energies (soft/hard X-rays) and for them instruments such as XGIS are more suitable. The presence of multiple instruments will be extremely valuable to cover the entire energy range typical of the prompt emission of short GRBs.

Figures 2 and 3 show the distribution in redshift of the joint detections, in the specific case of Fermi-GBM in synergy with ET and ET+CE, respectively. These figures are produced injecting 10⁵ BNS mergers (extracted from the BNS population, and thus following the astrophysical merger rate evolution with z), with the viewing angle uniformly distributed in the range 0 ° −15°. Even if we distribute the injections over a wider range of angles, the number of joint detections (given as fraction of γ-ray detections having an associated GW counterpart) does not change, since, while there are GW sources detectable at θ_v > 15°, there are no γ-ray detections for θ_v > 15° (see Figs. 4 and B.2). The γ-ray detections are computed considering the best fit configuration of parameters derived from MCMC posterior distribution, except for f_j, which has been assumed to be one, for better visibility. The value of f_j shifts the green and red histograms vertically but does not change their relative ratio. In the case of ET alone, the GW detector is so sensitive that up to z ≃ 0.8 − 1.0 the probability that a SGRB has a detectable GW counterpart is close to 100%. Adding CE to the network, the GW detection efficiency remains close to 100% for redshifts above the BNS merger peak. In this section we do not include the results for ET+2CE, since the addition of another CE would not further increase the joint detection efficiency in the redshift range accessible to γ-ray instruments for SGRBs. The γ-ray missions that are more suited to maximize the joint detection rates are those with large FOV and best sensitivity around MeV energies. However, another parameter to take into account is the localization accuracy by the γ-ray detectors (given in Tables 2 and 3), which defines what are the instruments that are able to drive the follow-up observations by ground-based telescopes, crucial for obtaining the source redshift and to completely characterize the multiwavelength emission of the source.

Fig. 2. Histogram of the joint ET+Fermi-GBM detections (red), together with the distribution of injected BNS (blue), the mergers detected by ET (orange), the mergers detected by Fermi-GBM (green). The histogram is normalized to Ninj = 105, which is the number of injected BNS mergers in the angle range 0 ° < θv < 15° and in the redshift range 0 < z < 4. For visualization purposes, the fraction of BNS producing a jet, fj, has been assumed to be one.

Fig. 3. Same as Fig. 2, but considering Fermi-GBM+(ET+CE).

Figure 4 shows the distribution of the ET+γ-ray joint detections as a function of the viewing angle, assuming Stru1, for three redshift bins: z < 0.5, 0.5 < z < 1.5, z > 1.5. Figure B.2 shows the same, but assuming Stru2. Again, as example, we consider Fermi-GBM as γ-ray detector, but we obtain consistent results for the other γ-ray satellites considered in the present paper. On the y axis we report the ratio N_JD/N_inj between the joint detections over the total injected BNS (assuming f_j = 1) per redshift bin and viewing angle bin. For the viewing angle we consider a linearly spaced grid, where the width of the single bin is θ_c/3. The uncertainty bands are obtained with 15 random extractions of the prompt parameters and each realization considers 5 × 10⁴ BNS injections with viewing angle 0 ° < θ_v < 15°. As before, a different assumed value of f_j just shifts the bands vertically. The plot shows a mild decrease of the ratio N_JD/N_inj for very small viewing angles (θ_vv ≲ θ_c/3). Such effect is due to the fact that SGRB viewed on-axis have very high values of peak energy, meaning that the bulk of the flux is above the Fermi-GBM band. Therefore, considering two identical SGRBs at the same redshift, the one viewed at θ_v ∼ θ_c appears slightly brighter than the one viewed at θ_v ∼ 0. Moreover, the number of sources per unit angle scales as sin(θ_v), therefore the small decrease of N_JD/N_inj around θ_v ∼ 0 is related to the paucity of sources contained in the first angle bin.

Fig. 4. Angle dependency of the joint GW+γ-ray detections for three redshift bins, considering ET+Fermi-GBM. The bands are obtained with random realizations extracted from the MCMC posterior samples. For each redshift bin, the y-axis gives the number of joint detections normalized to 5 × 104 BNS injections with viewing angle 0 ° < θv < 15°. We assume that all BNS produce a jet (fj = 1).

3.2. Joint detection of GWs and γ-rays from cocoon shock break-out

In the derivation of the joint GW+γ-ray detections, we have assumed that all the SGRBs share a common jet structure and that the γ-ray emission is given by the dissipation of the internal energy of the jet. However, for GRBs observed at large viewing angles, namely θ_v ≫ θ_c, the γ-ray emission from the shock break-out (SBO) of the jet should be taken into consideration. Indeed, even if the formation and successful break-out of a relativistic jet through the post-merger ejecta is not guaranteed, the formation of a cocoon at very wide angles is commonly expected in the BNS mergers (Ramirez-Ruiz et al. 2002; Nakar & Piran 2017; Gottlieb et al. 2018a). The shocks driven by the mildly relativistic expansion of the cocoon can produce γ-ray emission which is potentially detectable at least in the local Universe. Moreover, several studies found a remarkable agreement among the properties of the gamma-ray emission of GW 170817 and the shock-breakout model (Kasliwal et al. 2017; Gottlieb et al. 2018b; Nakar et al. 2018; Bromberg et al. 2018; Pozanenko et al. 2018). The γ-ray emission detected is not emitted by the wings of a structured jet, but rather by the shock breakout of the cocoon produced by the interaction of the jet with the NS merger ejecta (Gottlieb et al. 2018b; Bromberg et al. 2018). In order to investigate this scenario, as a reference, we take the spectrum (cut-off power law with photon index α = 0.62 ± 0.40 and cut-off energy E_c = 185 ± 62 keV) and the Fermi-GBM peak flux ((3.1 ± 0.7)×10⁻⁷ erg cm⁻² s⁻¹) of the γ-ray flash associated to GW 170817 (Goldstein et al. 2017) and we compute the maximum distance D_L, max at which this emission can be detected by the γ-ray instruments considered in the present work. In a first order approximation, the SBO emission can be considered isotropic, therefore we do not include any angle dependency in our treatment. In order to evaluate the number of γ-ray detections associated to SBO, we compute the number of BNS at a luminosity distance D_L < D_L, max. In our derivation, we assume that each BNS produces a SBO. Since the detection efficiency of any 3G GW detector is ∼100% for distances D_L < D_L, max, the number of γ-ray detections per year that we report corresponds also to the number of joint GW+γ-ray detections associated to SBO from BNS mergers. The results are reported in Table 4.

3.3. Joint detection of GWs and X-ray emission

In this section, we evaluate the detectability of the afterglow emission with future X-ray telescopes. Specifically, we show the expected rate of BNS mergers which will have both a detectable GW signal and an X-ray emission associated to the afterglow phase of the SGRB. For the X-ray emission, we consider two components: (1) the standard forward shock (FS) emission, (2) the HLE associated to the last photons emitted during the prompt phase. Under the assumption of a structured jet, the HLE produces an X-ray flux which can be comparable or even dominant with respect to the FS emission (Oganesyan et al. 2020; Panaitescu 2020). We evaluate the joint X-ray and GW detections considering satellites in survey mode as well as satellites pointing to the GW source. For the pointing strategy, we use the localization capabilities of ET alone and in the network of 3G observatories. Since the number of sources suitable for pointing observations result to be high (especially for network of 3G detectors), we also evaluate other source parameters estimated from the GW signals, such as viewing angle and distance, which can be used to down-select sources to be observed.

3.3.1. GW sky localization

Figure 5 shows the sky localization of ET, ET+CE and ET+2CE. The figures show the error on sky position, expressed in deg², as a function of redshift for one year of BNS injections. The injections are extracted from the BNS population described in Sect. 2.1, thus following the same redshift evolution of the merger rate. We show the sky localization both for sources with 0 < θ_v < 15° and for sources with no selection on θ_v. We chose θ_v maximum equal to 15° as it is consistent with the largest viewing angle up to which the X-ray emission from a Stru1 jet is observable by the analyzed WFX-ray satellites (see section 3.3.3). For the combination ET+CE (ET+2CE), a considerable fraction of sources detected at z < 1 (z < 3) has a sky localization smaller than 10 deg², which is small enough to enable prompt and efficient multiwavelength search for EM counterparts.

Fig. 5. Redshift distribution of the sky-localization uncertainty (given as 90% credible region) for three detector configurations: ET, ET+CE and ET+2CE. The absolute numbers are relative to one year of observation and assuming a duty cycle of 0.85 as described in the text. Plots on the right show BNS systems with a viewing angle 0 < θv < 15°. (a) ET, all angles. (b) ET, θv < 15°. (c) ET+CE, all angles. (d) ET+CE, θv < 15°. (e) ET+2CE, all angles. (f) ET+2CE, θv < 15°.

In Tables 5 and 6 we summarize the sky localization capabilities of ET, ET+CE and ET+2CE for systems observed with θ_v < 15° and for systems observed at all angles, respectively. The tables list the number of detections per year with sky localization uncertainty ΔΩ < 1, 10, 100, 1000 deg². Figure B.1 shows the redshift distribution of ΔΩ as scatter plot with a color code indicating the S/N for ET, ET+CE and ET+2CE, considering cases with θ_v < 15° and cases with no selection on θ_v.

Figure B.3 in appendix shows the cumulative distribution of the S/N for BNS detected with ET, ET+CE and ET+2CE with a sky uncertainty ΔΩ < 100 deg², for θ_v < 15° (panel a) and without selection on θ_v (panel b). The S/N distribution shows that, even adopting a larger S/N threshold, for instance S/N > 12, the number of sources with 8 < S/N < 12 is negligible and therefore not relevant for the evaluation of the joint GW+X-ray detection rate at least for well localized sources.

3.3.2. Viewing angle and luminosity distance from GWs

Here, we explore how the GW observations can give constraints about the viewing angle with respect to the axis perpendicular to the BNS orbital plane (assumed to be parallel to the jet axis) and the luminosity distance of the system. Figure 6 shows how the detected sources are distributed in the plane Δθ_v − θ_v, considering ET, ET+CE and ET+2CE. Δθ_v is the 1σ error on the viewing angle. Only sources at z < 1 and with Δθ_v < 40 deg are shown. For all the networks, a considerable fraction of sources have Δθ_v/θ_v < 0.1. Moreover, since the gradient of the GW amplitude is maximum at θ_v = 90 deg, the θ_v parameter is better constrained when the source is seen edge-on. This result is particularly important to prioritize the BNS events to be followed-up. Indeed, in the era of 3G GW detectors we expect hundreds of BNS detections per day. Due to the limited observational time at each observatory, with a pointing strategy only a limited fraction of GW-detected signals can be followed up. For example, aiming at detecting the high-energy emission, one can exclude all the sources edge-on for which the viewing angle estimate is more precise from the GW parameter estimation.

Fig. 6. Distribution of Δθv vs. θv relative to one year of observation, in the case of ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the redshift. We report only BNS mergers detected at z < 1 and with a Δθv < 40°. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.

Figure 7 shows the distribution in redshift of ΔD_L/D_L, where ΔD_L is 1σ error on the luminosity distance D_L. We consider ET, ET+CE and ET+2CE. Only sources with ΔD_L/D_L < 1 are shown. For ET operating alone the relative error can be below 10% only for sources with z ≲ 0.2, while for ET+CE and ET+2CE the same is true for z ≲ 1.5 and z ≲ 3, respectively. The plots also show that the constraint on distance is tighter for sources seen edge on. This is due to the strong degeneracy between luminosity distance and viewing angle and therefore the error on the first is significantly reduced when the second is well constrained.

Fig. 7. Distribution of the relative error on the luminosity distance DL as a function of redshift relative to one year of observation, in the case of ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the viewing angle. Only cases with ΔDL/DL < 1 and at z < 3 are shown. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.

3.3.3. Wide-FOV X-ray telescopes: survey and pointing observations

For the detection of the afterglow emission from our SGRBs population, we consider the following X-ray instruments: the Einstein Probe (EP, Yuan et al. 2018) which is scheduled for launch by the end of 2022, the ECLAIRs instrument on board of SVOM (Cordier et al. 2015) which is scheduled for launch in 2023, and three mission concepts, THESEUS (Amati et al. 2018, 2021), the Wide Field Imager on board of TAP (Camp et al. 2019), and the Gamow Explorer (White et al. 2021). The launch date of this last is expected to be in 2028–2032. A summary of the properties of these telescopes is shown in Table 7.

If the instrument operates in survey mode, the probability of detecting the source (if the flux is above the sensitivity limit) is given by ∼FOV/4π. Instead, if the instrument operates only in pointing mode, the light curve is monitored only after a time t_resp from the trigger time, which is the sum of the time to respond to a trigger from ground and the time to repoint the instrument. In order to evaluate how each instrument will sample the afterglow temporal evolution, we use the sensitivity curve (the sensitivity curve of EP is taken here⁵, while for all the other X-ray missions, we use the sensitivity curves reported in the references cited at the beginning of this section), which relates the minimum detectable flux as a function of the exposure time. The temporal sampling is computed using an iterative approach. Starting from the beginning of the light curve, we determine the length ΔT of each bin at time t_i such that

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Delta }T}{\displaystyle {\int }_{t_i}^{t_i+\mathrm{\Delta }T}F(t)\mathrm{d}t>F_{\lim }(\mathrm{\Delta }T),}\end{array}$$(11)

where F_lim(ΔT) is derived from the sensitivity curve. The afterglow light curve is classified as undetectable if, even increasing ΔT, the Eq. (11) is never satisfied. For Gamow, THESEUS and EP we consider sensitivity curves including a median Galactic absorption N_H = 5 × 10²⁰ cm⁻². For TAP the published sensitivity curve does not include this information.

For all the results reported in this section, we simulate 25 realizations of one year of BNS mergers, each with parameters extracted randomly from the MCMC posterior distributions. For each realization we consider a number of injections N_inj = min(N_req, 10⁴), where N_req is the number of BNS in one year that satisfy the selection criteria (e.g., about S/N, θ_v and ΔΩ). The choice N_inj = min(N_req, 10⁴) is made in order to have enough statistics and such that the full parameter space is correctly covered. Then the number of detections correspondent to one year is obtained rescaling by the factor N_req/N_inj. The parameters of the prompt model influence the afterglow light curve because both the FS and HLE brightness depend on the E_iso which is predicted by our model. First we evaluate the expected rate of joint GW+X-ray detections, considering the telescopes operating in survey mode. The results are reported in Table 8. Due to the good sensitivity down to 2 keV and the large FOV, THESEUS-XGIS can contribute substantially to the total amount of X-ray detections. Thanks to the larger FOV (1.1 sr) and the better sensitivity, EP gives a larger number of detections compared to THESEUS-SXI and TAP.

The probability of detection can be increased considering the possibility to point GW sources localized with enough precision through the GW signal. In such case, wide field X-ray (WFX-ray) telescopes can point to the sky error box and start the search of the X-ray counterpart. We explore this scenario simulating one year of BNS mergers and selecting only those that are detected with an error on sky position better than 100 deg². The typical FOV of WFX-ray telescopes is larger than 1 sr, so an error region of 100 deg² can be well contained inside the FOV of the instrument. A sky region of 100 deg² is around 1/10 of a typical FOV of a WFX-ray telescope. If, instead, we select sources localized better than 1000 deg² and considering that the bulk of sources detectable in the X-rays are at z ≲ 1, we expect, with respect to the cases with ΔΩ < 100 deg²: (1) a factor ∼10 more joint GW+X-ray detections with ET; (2) no substantial variation for ET+CE, since sources with 100  deg² < ΔΩ < 1000  deg² are located at z ≳ 1; 3) no substantial variation for ET+2CE, since N(ΔΩ < 100 deg²)∼N(ΔΩ < 1000 deg²) for z ≲ 1.

For the selected subsample of well localized sources, we predict the X-ray emission and we check if it is detectable with the telescopes listed before. The numbers of THESEUS-SXI and Gamow are identical since we consider for them the same sensitivity, and the slightly different FOV does not influence pointing observations. In Table 9 we show the number of expected joint GW+X-ray detections for one year of observation with ET, ET+CE and ET+2CE. We consider the configuration FS-SGRB for the FS parameters and we assume that the X-ray telescopes are able to point toward the sky position provided by GW detectors within 100 s from the BNS merger. A response time of 100 s is a short amount of time to communicate the trigger from the ground to the satellites, but the possibility of accessing lower frequencies by the next generation GW detectors will make it possible premerger alerts. We evaluated that, in the case of ET, a fraction of ∼30% and ∼60% of sources localized better than 100 deg² at the merger are above the detection threshold (S/N = 8) and with a sky-localization better than 1000 deg², which is well within the FOV of the satellites analyzed here, 20 and 10 min before the merger, respectively (see Table B.2). For ET+CE (ET+2CE) this fraction is ∼5% (6%) 10 min before the merger, leaving the absolute number of triggers to be followed up of order of several hundreds. The typical exposure time for detection in X-rays is $t_{\exp }={66}_{-42}^{+249}$ s for SXI and Gamow, $t_{\exp }={96}_{-29}^{+219}$ s for TAP, $t_{\exp }={239}_{-173}^{+171}$ s for EP, where the uncertainties are reported at 1σ level of confidence and they are computed considering a random sampling of the MCMC posterior distribution.

Considering ET alone and FS configuration (FS-SGRB), Tables 8 and 9 show that the survey mode is more promising than the pointing mode in terms of number of GW+X-ray detections per year. For ET+CE and ET+2CE, due to the larger number of well localized sources, the number of joint detections is considerably larger than the survey mode case. Hence, the exploitation of sky localization provided by GW detectors can enhance substantially the probability of identifying the EM counterpart of BNS mergers. We point out that the reported numbers for pointing mode do not take into account the duty factor of the instruments, which is assumed to be 100%. Moreover, each X-ray mission, depending on the scientific objectives, can dedicate only a limited amount of time for the follow-up of GW triggers. Therefore, if we call ϵ_GW the fraction of observational time dedicated to GW follow-up and ϵ_DC the duty cycle, a more realistic estimate of the number of joint detections would be N_JD, real = ϵ_GWϵ_DCN_JD, where N_JD are the numbers reported in this section about the pointing mode. Regarding the survey mode, instead we just have N_JD, real = ϵ_DCN_JD. Another pivotal point of the future pointing strategy, which can strongly influence the real joint detection efficiency, is the trigger selection. As shown in Table 5, the number of sources localized better than 100 deg² are of order of 100 for ET and of 10⁵ for ET+CE and ET+2CE. In order to not lose observational time, the triggers to be observed will require to be prioritized. A strategy based on the knowledge of the satellite detection efficiency, the distance and viewing angle from the GW signals (see Sect. 3.3.2) will be mandatory to select the triggers with a higher chance of joint detection. Therefore the reported numbers for pointing mode should be considered as potential detections, which can be achieved only provided that an efficient selection of the GW triggers is performed to restrict the follow-up to cases with highest probability of detection.

In Table 10 we consider longer response time (relaxing the assumption of t_resp = 100 s) and we show how the detection rate significantly decreases if the telescope needs more time to respond to the trigger and point the GW localization. The significant decrease of the chance to detect the afterglow emission pointing with one hour delay is due to a combination of the rapid decline of the GRB X-ray flux and the sensitivity of wide field X-ray instruments. Among the population of detected BNS, nearby events similar to GW170817 are a tiny fraction, while the large majority of joint detections comes from on-axis (or slightly off-axis) short GRBs distributed at cosmological distances, with an afterglow light curve which peaks at early times and followed by rapid drop of the flux. For these cases a rapid response is fundamental to catch the first decaying part of the X-ray emission. In Table 11, we show how the prediction of joint GW+X-ray detections depends on the assumption of FS parameters and on the inclusion of HLE. In this case, we assume ET+2CE as GW network and t_resp = 100 s. In the specific, we test three possible setups, where we compute the X-ray light curve: (1) including HLE + FS emission, with FS-SGRB configuration; (2) including HLE + FS emission, with FS-GW17 configuration; (3) including only FS emission, with FS-SGRB configuration. Our results show that the assumption of FS parameters like GW 170817 leads to less bright FS emission, mainly due to lower value of ϵ_B. Moreover, the tables show that the inclusion of HLE is non negligible in the estimation of detectability of the X-ray afterglow. The reasons are mainly twofold: (1) the HLE produces an initial steep decay phase which is well visible in X-rays and exceeds the FS contribution at least in the first tens of seconds, (2) the HLE associated to a structured jet typically produces a plateau phase for sources observed at θ_v ≲ θ_c and this contribution can be as relevant as the FS emission itself.

Finally, we test how the expected joint detection rate of GW+X-rays from BNS mergers depends on the assumption of the jet structure. The specific structure Stru1 assumed so far, namely the one derived from GW 170817, has a very steep off-core profile, meaning that the probability of detecting prompt and/or afterglow emission from SGBRs viewed at θ_v > θ_c is lower than the one of detecting a shallower off-core jet profile, such as Stru2. Table 12 shows the difference in joint GW+X-ray detections in pointing mode assuming Stru1 and Stru2. The overall rate of joint GW+X-ray detections for Stru2 is, within the uncertainties, comparable with the one of Stru1. This is due to the fact that, though the X-ray emission is detectable at larger viewing angles for Stru2, the best fit value for the fraction of BNS producing a jet, f_j, derived from the Monte Carlo is smaller for Stru2 than the one of Stru1. Figures 8 and 9 show the distribution of viewing angle of GW+X-rays joint detections as a function of redshift, adopting the jet structure Stru1 and Stru2, respectively. For each detection, we use a color code for the sky-localization uncertainty. The plot is realized simulating 1000 GW+X-rays detections (ET+2CE + SXI) which satisfy the requirement ΔΩ < 100 deg², assuming t_resp = 100 s and FS-SGRB configuration. Below z ∼ 0.5 − 1, for both the structures, the X-ray emission is detectable up to viewing angles much larger than θ_c = 3.4 deg. In the bottom panel of Fig. 8, we show how the GW+X-rays joint detections are distributed in redshift and we distinguish between cases detected at θ_v < 5° and θ_v > 5°. We adopt θ_v = 5° as a threshold because, according to Fig. 4, this is the maximum angle up to which prompt emission is detectable by Fermi-GBM (similarly for XGIS onboard of THESEUS), considering an average on redshift. This implies that all the sources detected in X-rays at θ_v > 5° have a non detectable prompt emission in γ-ray. In the case of Stru1, sources with θ_v > 5° represent ∼40% of all the GW+X-rays joint detections. Similarly, for Stru2 the threshold angle for detectability of prompt emission is θ ∼ 10° and ∼50% of all the GW+X-rays joint detections is above this limit. Independently on the assumption on the jet structure, this result means that there is a significant number of off-axis events concentrated below z ∼ 0.5 − 1 which are only detectable with WFX-ray instruments. In terms of absolute numbers of off-axis GW+X-rays joint detections, a combination of ET+2CE + SXI would observe few a tens of such events per year. For TAP and Einstein Probe a similar fraction (∼40 − 50%) of events are expected to be observed off-axis. This result demonstrates the fundamental role of WFX-ray telescopes in synergy with GW detectors for the detection of off-axis emission from BNS mergers that would be otherwise undetactable for γ-ray instruments. We finally evaluated how the choice of a larger value of the angular extension of the jet wings, θ_w (used in Afterglowpy), could increase the number of detections. Increasing θ_w also the maximum viewing angle at which the source is detectable increases. We estimated that increasing θ_w from 15 deg (30 deg) to 45 deg, the number of joint detections increases of a factor < 5% (∼10%) for Stru1 (Stru2).

Fig. 8. Redshift distribution of the joint GW+X-ray detections. Panel a: shows the distribution of the viewing angle as a function of redshift for a sample of 1000 joint GW+X-ray detections, considering ET+2CE and SXI and under the assumption Stru1. Only BNS detections with GW sky localization ΔΩ < 100 deg2 are selected. The assumed jet core angle is θc = 3.4 deg. The color bar indicates the GW sky localization uncertainty of each detection. Panel b: the distribution in redshift of the same sample, where we distinguish between detections with θv < 5° and θv > 5°. The histogram is normalized to the number of total joint detections with ΔΩ < 100 deg2.

Fig. 9. Same as in Fig. 8, but for Stru2.

3.3.4. Medium and small-FOV sensitive X-ray observatories

When the GW sources are localized at the level of arcminutes the X-ray afterglow can be directly detected by pointing small-FOV instruments. However, there are no BNSs localized at arcminute precision through the GW signal when ET is operating as a single detector, and the number is also negligible when ET is in a network of 3G detectors. When the GW localization is of order of 1–10 deg², medium-FOV X-ray instruments can scan the entire localization region by performing mosaic observations. Considering the GW sources localized better than 10 deg² from Table 6, the sources that can be followed-up by medium-FOV instruments are around 10 when ET is operating as a single detector, but several thousands (hundred of thousands) when the ET+CE (ET+2CE) network is observing. However, in the case of ET+CE, the distances of these sources are limited to z < 1.2. For larger z the way to detect the X-ray counterpart remains through the use of wide-FOV γ-ray and/or X-ray satellites. In this case, the more sensitive small-FOV satellites can operate at a later time by pointing the precise localization provided by the wide-FOV satellites, following-up and characterizing the detected counterparts.

During the activity period of 3G GW detectors, the mission Athena is expected to be operative (Nandra et al. 2013). The exceptional sensitivity of Athena X-ray telescope and its focal plane instruments will allow us to detect the faint afterglow emission from SGRBs and GW counterparts hours to years (for closer events) after the merger and to make high resolution spectra (Piro et al. 2021). Athena will host the Wide Field Imager (WFI) with a FOV of 0.4 deg² and the X-IFU high resolution spectrometer with a FOV of 5 arcmin equivalent diameter. The WFI is able to directly point GW signals with sky localization uncertainties smaller than its FOV (a few tens per year for ET+CE and thousands for ET+2CE considering all the orientations of the binary systems), and in addition is capable to carry out a mosaic of sky regions extended up to 10 deg². On the other hand, X-IFU requires a precise localization (≃2 arcmin) to directly point the source. Even considering a GW network with ET+2CE, there are no sources with such a small uncertainty on the GW sky location. Therefore, instead of being directly triggered by the GW detection, X-IFU can be used either for following up precise localizations delivered by WFX-ray telescopes or by the Athena WFI itself. Taking into account that Athena will likely have a Time of Opportunity (TOO) response time of 4 h, we find that the totality of the sources jointly detected by ET+2CE and the WFX-ray telescopes can be detected and followed-up by Athena X-IFU.

In order to evaluate the Athena-WFI capabilities to detect GW counterparts, we consider the sources localized by the GW detectors within a region of 10 deg². We assume an average number of exposures to fully cover the GW sky region equal to N_exp ∼ ΔΩ_GW/FOV_WFI, where FOV_WFI = 0.4 deg², and an exposure time T_exp = 10⁴ s for each mosaic observation. For simplicity, we work in the approximation that the source is randomly located within the 90% GW error region and that each position is equi-probable. This implies that, if the GW region is tiled in N_exp subregions, the probability that the source is located inside the kth subregion is P(k) = 1/N_exp and is the same for all the subregions. We consider a TOO response time T_TOO of 4 h and we include a dead time interval T_d of 30 min corresponding to the time necessary to slew from one subregion of the mosaic to the adjacent one⁶. Hence, the time necessary to point and detect the source is T(k) = T_TOO + k × (T_exp + T_d). Operatively, for each source we extract a random number k in the range [0, N_exp] and we check if at time T(k) after the GW trigger the X-ray afterglow emission is above the sensitivity of Athena-WFI corresponding to 10⁴ s of exposure.

Given the several science goals of Athena, only a fraction of its observational time can be realistically dedicated to the follow-up of GWs; specifically we assume one month over one year of observations. While for ET alone the number of sources with ΔΩ < 10 deg² is around 10, and all of them can in principle be followed up by Athena in this total amount of time, the number of the sources detected by ET+CE with ΔΩ < 10 deg² is of several thousands. In order to maximize the effectiveness of the search and not lose time on sources without any chance of detection, a further selection of the sources to be followed-up results to be necessary. We select the sources by excluding all those with θ_v > 50° (for which the error on θ_v is relatively small, see Fig. 6, and thus can be reliably considered off-axis) and z > 0.5. We select randomly among these sources (∼130) up to reach a total amount of time of one month of observation, and find that Athena-WFI is able to detect $N_{\mathrm{GW}+\mathrm{Athena}}={ϵ}_{\mathrm{FOR}}\times (6_{-3}^{+2})$ events per year, where ϵ_FOR is a factor that takes into account the TOO efficiency and it corresponds to a field of regard of ∼50% (i.e., only a fraction ϵ_FOR of the events can be successfully followed-up by Athena). Including another CE in the GW network does not increase significantly the number of Athena detections ($N_{\mathrm{GW}+\mathrm{Athena}}={ϵ}_{\mathrm{FOR}}\times 8_{-2}^{+3}$ detections per year), because even increasing the number of sources well localized, the maximum number of Athena detections is anyway limited by the time dedicated to TOO. If we exclude more sources by lowering the threshold on the viewing angle to θ_v < 30° instead of θ_v < 50°, the number of joint detection increases to ${ϵ}_{\mathrm{FOR}}\times {14}_{-3}^{+4}$ yr⁻¹ (${ϵ}_{\mathrm{FOR}}\times {27}_{-3}^{+3}$ yr⁻¹) for ET+CE (ET+2CE). The knowledge about the luminosity distance and the viewing angle coming from GW analysis can considerably help for the selection of a golden sample of events that have a larger probability to be detected, giving the possibility to maximize the GW+X-ray joint detection efficiency.

4. Discussion and conclusions

Our study evaluates the perspectives for multimessenger astronomy in the ET era focusing on the search of high-energy signals, which is a unique way to detect EM counterparts of GW signals at high redshift. We estimate the expected rate of joint detections of GWs and γ-ray and/or X-ray signals from BNS mergers considering prompt and afterglow emissions, different GW detector configurations (ET, ET+CE, ET+2CE), and different observational strategies (survey and pointing mode) for several satellites. Our theoretical framework starts from an astrophysically-motivated population of BNS mergers and predicts the high-energy signals associated with these mergers making the BNS population able to reproduce all the statistical properties of currently observed SGRBs. This approach makes it possibile to reproduce not only the observational features of SGRBs, such as luminosity, duration and spectral properties, but also the average detection rate with current γ-ray instruments, such as Fermi-GBM. This enables us to normalize the number of joint detections, evaluate the fraction of BNS producing a jet, and be less dependent of the BNS merger absolute number, which is still largely uncertain.

Our main results are summarized as follows:

• Regarding the joint detection of GWs and γ-ray emission, we find that already considering ET alone, more than 60% of all the SGRBs having a detectable prompt emission will also have a detectable GW counterpart. This percentage approaches 100% if we include CE as an additional interferometer in the network. From a few tens to hundred joint detections are expected per year, mainly depending on the FOV of the high-energy satellites. However, although increasing the number of joint detections is important for statistical studies, characterizing the sources and identifying the host galaxies is of primary importance. As shown by Swift in many years of observations, this is possible and effective when the source is localized at arcmin level enabling to drive the follow-up of ground-based telescopes. Mission concepts such as THESEUS-XGIS, will be essential to detect and localize high-energy counterpart of BNS mergers with enough precision (arcmin uncertainty). This is particularly important when the sources are located at large redshift, namely z > 1 where also in the optimistic scenario of a network of 3G GW detectors, the GW sky localization will be larger than 1 deg² from the GW signals. Instruments such as HERMES (full constellation of cube-satellites) can also be a good compromise of number of detections and sky localization uncertainties which are however order of a few deg².

• Considering the possible shock breakout of the cocoon produced by the interaction of the jet with the NS merger ejecta, we show that sensitive (more than the current ones) and wide FOV γ-ray instruments are required.

• Regarding the joint detection of GWs and X-ray emission, we demonstrate the important role of WFX-ray missions, such as SVOM, Einstein Probe, Gamow, THESEUS, and TAP. We predict tens of detections per year when these instruments operate in survey mode. The joint detection rate in survey mode is limited by the chance of having the source inside the FOV at the moment of the merger. We propose an additional observational strategy, the pointing strategy, which could in principle enhance the probability of joint detection exploiting the information about sky localization provided by GW instruments (see below).

• ET operating as a single observatory can detect several hundreds of BNS mergers per year with ΔΩ < 100 deg², among them about ∼70 BNS/yr have a viewing angle θ_v < 15°, namely sources with potentially detectable high-energy emission. The inclusion of CE in the network significantly improves the localization capabilities by reaching several thousands (hundreds with θ_v < 15°) of BNS detections per year with ΔΩ < 10 deg². The network of CE+ET is able to localize a few hundreds (a few tens with θ_v < 15°) with ΔΩ < 1 deg² up to redshift 0.3. This number increases to thousands (hundreds on-axis) per year up to redshift 1 with the network of ET+2CE.

• Considering WFX-ray telescopes slewing to the sky position provided by GW detectors for sources with a good sky localization (better than 100 deg²), we evaluate that the X-ray afterglow emission from ten to hundreds SGRBs could be detected during one year of observation assuming a response time t_resp = 100 s.

• A rapid communication of the GW sky location and a rapid repointing of the high-energy satellites is crucial to maximize the efficiency of the poiting strategy. Indeed, due to the typical rapid decline of the X-ray flux of SGRBs, we showed that a longer delay time between the merger and the beginning of X-ray monitoring would decrease the probability of detection by a factor ∼10 or more going from t_resp = 100 s to t_resp = 1 h. The low-frequency sensitivity of ET will make it possible to provide a good sky localization even several minutes before the merger (see Chan et al. 2018; Li et al. 2022; Nitz & Dal Canton 2021, Banarjee et al., in prep.). This will allow us to send information about the sky location of the source well before the merger and X-ray telescopes can realistically start earlier the slewing procedure. We evaluate that a not negligible fraction of GW events, localized within 100 deg² at merger, are detectable tens of minutes before the BNS merger with a sky localization well within the FOV of the WFX-ray satellites.

• While the number of detections for WFX-ray satellites in survey mode is reliable, the number of detections in pointing is optimistic and relies on a perfect prioritization of the trigger to be followed (see below).

• We show that using WFX-ray instruments makes it possible to detect a significant fraction of BNS mergers which are too off-axis to have a detectable γ-ray emission and that would be otherwise missed in absence of these instruments. The presence of WFX-ray monitors is also fundamental to provide the spectral coverage at smaller energies than γ-rays and to precisely localize the source. The WFX-ray telescopes can provide sky location accuracy of ∼ arcmin, which enables to trigger follow-up observations by ground-based optical telescopes, radio arrays and exceptionally sensitive X-ray instruments, such as the X-IFU instrument onboard of Athena. We evaluate that the totality of the sources jointly detected by the GW detectors and the WFX-ray telescopes can be detected and followed-up by the Athena X-IFU.

• The networks of ET+CE and ET+2CE are expected to localize a large number of BNS signals with ΔΩ < 10 deg². An instrument such as the WFI (FOV of 0.4 deg²) on-board of Athena is able to carry out a mosaic of these sky regions by detecting from a few to a few dozen signals by using 1 month of Athena observations. The joint detection number comes from an observational strategy which maximizes the chance of detections by removing the BNS signals with smaller probability to have observable jet. In particular, this is done by removing BNS with larger viewing angles and distances.

• Due to the high number of BNS events expected to be released by the next generation GW detectors, for the pointing observational strategy, it will be crucial to select and prioritize the ones to be followed-up. The selection based on sky-localization will not be enough to avoid to lose observational time of the high-energy satellites when the network ET+CE and ET+2CE will operate. On the basis of the different scientific goals, specific selection based on the GW source parameters such as distance and viewing angles, but also premerger sky localization, will be mandatory to maximize the joint detection efficiency. This shows the importance in the ET era to send in low-latency (updating) information about GW parameter estimation.

Our work shows that to maximize the science return of 3G GW detectors in the multimessenger context, it is necessary to develop sensitive wide-FOV high-energy instruments able to detect but also localize a large number of short GRBs. Instruments, such as THESEUS, combining γ-ray and WFX-ray telescopes, are of primary importance to guarantee the observation of γ-ray and X-ray counterparts, by detecting ∼10 − 100 well localized SGRBs up to high redshifts, for which also the host galaxy can be identified. This is of primary importance for evaluating the cosmological parameters and testing general relativity at cosmological scales. WFX-ray telescopes are required to detect GRB off-axis and make it possible to explore the jet structure, its interconnection with the ejecta, and the nature of low-luminosity SGRBs with a unique level of detail. Medium-size X-ray instruments, such as Athena WFI, will become effective (operating mosaic of the GW sky localization) for the X-ray counterpart search when a network of GW detector will observe. Mosaic observations with deep exposure by Athena WFI can enhance the possibility to detect off-axis X-ray emission. Sensitive instruments such as Athena X-IFU will be crucial to follow-up and characterize the X-ray emission at later times with respect to the merger.

Acknowledgments

We thank Yann Bouffanais, Boris Goncharov, Andrea Maselli, Om Sharan Salafia, Giulia Stratta for valuable discussions. We thank Michele Maggiore, Francesco Iacovelli, Michele Mancarella, and Stefano Foffa for checking the gravitional-wave results and for very helpful discussions. We acknowledge Stefano Bagnasco, Federica Legger, Sara Vallero, and the INFN Computing Center of Turin for providing support and computational resources. We thank Ssohrab Borhanian, Bangalore Sathyaprakash, Man Leong Chan, Yufeng Li, Ik Siong Heng for useful information and cross-check on the gravitational wave data analysis. We thank Lorenzo Amati, Fabrizio Fiore, Paul O’Brien and Nicholas White for useful information about some of the X-ray instruments. We thank Maria Grazia Bernardini for the information about SVOM. MB and SR acknowledge financial support from MIUR (PRIN 2020 grant 2020KB33TP_001). MB and GO acknowledge financial support from the AHEAD2020 project (grant agreement n. 871158). BB, GG and MB acknowledge financial support from MIUR (PRIN 2017 grant 20179ZF5KS). GG acknowledges ASI-INAF agreement n. 2018-29-HH. MM and FS acknowledge financial support from the European Research Council for the ERC Consolidator grant DEMOBLACK, under contract no. 770017. The research leading to these results has been conceived and developed within the Einstein Telescope Observational Science Board. The latest public version of MOBSE can be downloaded from this repository. GWFish is publicly available at this repository (https://github.com/janosch314/GWFish)

References

Appendix A: Computation of the observables in the structured jet scenario

In this section we show how the peak flux, the fluence and the isotropic energy are derived in the structured jet framework. The derivation is done assuming a unique jet structure, with fixed values of R₀ and Γ, defined in the main text. For a given viewing angle θ, we write the rest frame photon spectrum as

$$\begin{array}{c}\hfill N(E,\theta ,t)=N_0(\theta ,t)f(E/E_0(\theta ,t)),\end{array}$$(A.1)

where we included the time dependency only in the normalization, assuming that there is not a strong spectral evolution during the GRB pulse. f(E/E₀(θ)) is the spectral shape and it is assumed to be a smoothly broken power law, whose spectral slopes are assumed to be angle independent. f is normalized in such a way that f(1) = 1. E₀(θ, t) is the peak energy and we express it as

$$\begin{array}{c}\hfill E_0(\theta ,t)=E_0(\theta ,t=t_p)P_E(t/t_p,\theta ),\end{array}$$(A.2)

where t_p is the peak time of the pulse and P(t/t_pθ) describes the temporal profile, in principle angle-dependent. Analogously, we express the time behaviour of the photon spectrum normalization as:

$$\begin{array}{c}\hfill N_0(\theta ,t)=N_0^p(\theta )P_N(t/t_p,\theta ).\end{array}$$(A.3)

In both cases, we approximate the temporal profiles of peak energy and normalization as

$$\begin{array}{c}\hfill P_X(t)=\{\begin{array}{cc}t/t_p\hfill & 0<t<t_p\hfill \\ {(t/t_p)}^{-{\alpha }_X(\theta )}\hfill & t>t_p\hfill \end{array}X=E,N.\end{array}$$(A.4)

At the peak time we write the photon spectrum as

$$\begin{array}{c}\hfill N^p(E,\theta )=N_0^p(\theta )f(E/E_0(\theta )).\end{array}$$(A.5)

We define also

$$\begin{array}{c}\hfill R_E(\theta )=E_0(\theta )/E_0(\theta =0)=E_0(\theta )/{\widehat{E}}_0\end{array}$$(A.6)

and

$$\begin{array}{c}\hfill R_F(\theta )=\frac{N^p{(\theta )|}_{E=E_0(\theta )}}{{\widehat{N}}^p{|}_{E={\widehat{E}}_0}}.\end{array}$$(A.7)

The functions α_E(θ), α_N(θ), R_E(θ) and R_F(θ) are computed numerically using the structured jet model. Hereafter the notation $\widehat{x}$ indicates that the quantity x is evaluated at viewing angle θ = 0. In the observer frame, the peak photon spectrum is:

$$\begin{array}{c}\hfill N^p(E,\theta )={(1+z)}^2{N}_0^p(\theta )f\left(\frac{1}{R_E}\frac{(1+z)E}{{\widehat{E}}_0}\right)\end{array}$$(A.8)

and using eq. A.7

$$\begin{array}{c}\hfill N^p(E,\theta )={(1+z)}^2{\widehat{N}}_0^p{R}_{F}f\left(\frac{1}{R_E}\frac{(1+z)E}{{\widehat{E}}_0}\right).\end{array}$$(A.9)

The value of ${\widehat{N}}_0^p$ is derived from the rest frame on-axis isotropic energy ${\widehat{E}}_{\mathrm{iso}}=E_t/(1-cos{\theta }_c)$, where E_t is extracted from the assumed probability distribution defined in the main text. In the specific,

$$\begin{array}{c}\hfill {\widehat{E}}_{\mathrm{iso}}=4\pi D_L^2{\displaystyle \int dt{\int }_{1\text{keV}}^{{10}^4\text{keV}}{F}_{\nu }(t)dE,}\end{array}$$(A.10)

where D_L is the luminosity distance and F_ν is the flux density. Using eq. A.1:

$$\begin{array}{c}\hfill {\widehat{E}}_{\mathrm{iso}}=4\pi D_L^2{\displaystyle \int dt{\widehat{N}}_0(t){\widehat{E}}_0^2(t){\int }_{1\text{keV}/{\widehat{E}}_0(t)}^{{10}^4\text{keV}/{\widehat{E}}_0(t)}xf(x)dx,}\end{array}$$(A.11)

which can be written as

$$\begin{array}{c}\hfill {\widehat{E}}_{\mathrm{iso}}=4\pi D_L^2{\widehat{N}}_0^p{\displaystyle \int {\widehat{P}}_N(t/t_p){\widehat{E}}_0^2(t)I_E(t)dt,}\end{array}$$(A.12)

where $I_E(t)={\int }_{1\text{keV}/{\widehat{E}}_0(t)}^{{10}^4\text{keV}/{\widehat{E}}_0(t)}xf(x)dx$.

This procedure allows us to write ${\widehat{N}}_0^p$ in terms of known quantities, namely

$$\begin{array}{c}\hfill {\widehat{N}}_0^p=\frac{{\widehat{E}}_{\mathrm{iso}}}{4\pi D_L^2\int {\widehat{P}}_N(t/t_p){\widehat{E}}_0^2(t)I_E(t)dt}.\end{array}$$(A.13)

Once ${\widehat{N}}_0^p$ is determined, the peak photon flux integrated in the energy band of a given instrument is

$$F_p(\theta )={\displaystyle {\int }_{(1+z)E_1}^{(1+z)E_2}{N}^p}(E,\theta )dE$$(A.14)

$$=(1+z){\widehat{N}}_0^p{R}_F(R_E{\widehat{E}}_0){\displaystyle {\int }_{(1+z)E_1/R_E{\widehat{E}}_0}^{(1+z)E_2/R_E{\widehat{E}}_0}f}(x)dx.$$(A.15)

The duration of the burst t₉₀ is computed imposing that the 10-1000 keV fluence accumulated up to t₉₀ is 90% of the total fluence, namely

$${\displaystyle {\int }_0^{t_{90}}d}t{\displaystyle {\int }_{(1+z)\text{ keV}}^{(1+z){10}^4\text{ keV}}E}N(E,\theta ,t)dE$$(A.16)

$$=90\%{\displaystyle {\int }_0^{\infty }d}t{\displaystyle {\int }_{(1+z)\text{ keV}}^{(1+z){10}^4\text{ keV}/}E}N(E,\theta ,t)dE.$$(A.17)

For computational reasons, the last integral in time is performed up to a time t_max, which is the maximum time for which the γ-ray emission is still detectable by the instrument. Finally the isotropic luminosity is defined as E_iso/t₉₀.

Appendix B: Detection efficiency for different neutron-star mass distribution

In Fig. B.4 we show the detection efficiency (detected signals over injected ones) of ET as a function of redshift for a Gaussian mass distribution as observed in Galactic double neutron-star systems in radio (Özel & Freire 2016) and a uniform distribution of NS mass consistent with gravitational waves observations (The LIGO Scientific Collaboration 2021b). For the Gaussian distribution we fixed the mean value to ⟨M⟩ = 1.33M_⊙ and the standard deviation to σ_m = 0.1M_⊙. The detection efficiency obtained with the two distributions is consistent up to a redshift equal to 0.5, and then the Gaussian distribution leads to a lower number of GW detections accountable by the fact that the GW signal amplitude scales with the mass.

Fig. B.1. Sky localization (90% credible region in deg2) as a function of the redshift for one year of BNS injections with viewing angle in the range 0 ° < θv < 15° and considering ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the S/N of each detection. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.

Fig. B.2. Same as Fig. 4, but for Stru2.

Fig. B.3. Cumulative distribution of the S/N for BNS mergers localized better than 100 deg2 for the three GW networks considered in this work. Panel (a) shows cases with θv < 15°, while in panel (b) all angles are considered.

Fig. B.4. Detection efficiency of ET as a function of redshift assuming a Gaussian distribution (⟨M⟩ = 1.33M⊙, σM = 0.1M⊙) and a uniform distribution (M ∈ [1.0 − 2.5]M⊙). Only cases with θv < 15° are considered. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text. The ratio between the orange bars and the blue bars corresponds to the function ϕ(z) defined in the text.

In order to evaluate the impact on joint detections, we estimate, as example, the loss of joint GW+γ-rays detections with ET+Fermi-GBM. We assume the probability of having a joint detection as

$$\begin{array}{c}\hfill P_{\mathit{JD}}(z)=P_{\mathit{EM}}(z)P_{\mathit{GW}}(z),\end{array}$$

where P_EM(z) is the probability of detecting the EM counterpart, and P_GW(z) the probability of detecting the associated GW signal. Each probability corresponds to the ratio between the number of detected and injected sources.

$$\begin{array}{c}\hfill P_i(z)=\frac{N_{det,i}(z)}{N_{\mathit{inj}}(z)},i=GW,EM.\end{array}$$

The assumption of NS mass distribution only affects P_GW(z), since we specified in the main text that we neglect any dependence on the NS mass for the formation of a jet. We call P_GW, flat(z) and P_GW, Gauss(z) the GW detection efficiency for flat and Gaussian distributions, respectively. The variation on the rate of joint detections can be computed as

$$\begin{array}{c}\hfill \frac{d{N}_{JD,Gauss}}{\mathit{dz}}(z)=\frac{d{N}_{\mathit{inj}}(z)}{\mathit{dz}}{P}_{JD,Gauss}(z)=\end{array}$$

$$\begin{array}{c}\hfill =\varphi (z)\frac{d{N}_{\mathit{inj}}(z)}{\mathit{dz}}{P}_{\mathit{EM}}(z)P_{GW,flat}(z)=\varphi (z)\frac{N_{JD,flat}}{\mathit{dz}},\end{array}$$

where we wrote P_GW, Gauss(z) = ϕ(z)P_GW, flat(z). The total variation, integrated in redshift, is

$$\begin{array}{c}\hfill N_{JD,Gauss}={\displaystyle \int }\varphi (z)\frac{N_{JD,flat}}{\mathit{dz}}dz.\end{array}$$

In Fig.B.4 the function ϕ(z) corresponds to the ratio between the two histograms. Finally, knowing the function ϕ(z) and the redshift distribution of joint detections, we derive the variation in joint detection rate for a Gaussian mass distribution. We obtain that the joint GW+γ-rays detections with ET+Fermi-GBM decreases of a factor ∼27% (from 33/yr to 24/yr), assuming a Gaussian mass distribution.

Fig. B.5. Corner plot of the posterior distribution of the parameters adopted for the prompt emission model. The assumed structure is Stru1. The number of MCMC steps is chosen such that the auto-correlation time reaches a plateau, as described in Foreman-Mackey et al. (2013).

Fig. B.6. Same as Fig. B.5, adopting Stru2.

All Tables

All Figures

	Fig. 1. X-ray light curves of SGRBs with viewing angle 0 ° < θv < 15° predicted by our model. The lines are produced with a random extraction from the posterior distribution of prompt parameters. In the computation we assume Stru1 and we include both HLE and FS (configuration FS-SGRB). The color of each line indicates the viewing angle. The X-ray luminosity is computed at 5 keV.
In the text

Fig. 2. Histogram of the joint ET+Fermi-GBM detections (red), together with the distribution of injected BNS (blue), the mergers detected by ET (orange), the mergers detected by Fermi-GBM (green). The histogram is normalized to Ninj = 105, which is the number of injected BNS mergers in the angle range 0 ° < θv < 15° and in the redshift range 0 < z < 4. For visualization purposes, the fraction of BNS producing a jet, fj, has been assumed to be one.

In the text

	Fig. 3. Same as Fig. 2, but considering Fermi-GBM+(ET+CE).
In the text

	Fig. 4. Angle dependency of the joint GW+γ-ray detections for three redshift bins, considering ET+Fermi-GBM. The bands are obtained with random realizations extracted from the MCMC posterior samples. For each redshift bin, the y-axis gives the number of joint detections normalized to 5 × 104 BNS injections with viewing angle 0 ° < θv < 15°. We assume that all BNS produce a jet (fj = 1).
In the text

Fig. 5. Redshift distribution of the sky-localization uncertainty (given as 90% credible region) for three detector configurations: ET, ET+CE and ET+2CE. The absolute numbers are relative to one year of observation and assuming a duty cycle of 0.85 as described in the text. Plots on the right show BNS systems with a viewing angle 0 < θv < 15°. (a) ET, all angles. (b) ET, θv < 15°. (c) ET+CE, all angles. (d) ET+CE, θv < 15°. (e) ET+2CE, all angles. (f) ET+2CE, θv < 15°.

In the text

	Fig. 6. Distribution of Δθv vs. θv relative to one year of observation, in the case of ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the redshift. We report only BNS mergers detected at z < 1 and with a Δθv < 40°. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.
In the text

	Fig. 7. Distribution of the relative error on the luminosity distance DL as a function of redshift relative to one year of observation, in the case of ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the viewing angle. Only cases with ΔDL/DL < 1 and at z < 3 are shown. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.
In the text

Fig. 8. Redshift distribution of the joint GW+X-ray detections. Panel a: shows the distribution of the viewing angle as a function of redshift for a sample of 1000 joint GW+X-ray detections, considering ET+2CE and SXI and under the assumption Stru1. Only BNS detections with GW sky localization ΔΩ < 100 deg2 are selected. The assumed jet core angle is θc = 3.4 deg. The color bar indicates the GW sky localization uncertainty of each detection. Panel b: the distribution in redshift of the same sample, where we distinguish between detections with θv < 5° and θv > 5°. The histogram is normalized to the number of total joint detections with ΔΩ < 100 deg2.

In the text

	Fig. 9. Same as in Fig. 8, but for Stru2.
In the text

	Fig. B.1. Sky localization (90% credible region in deg2) as a function of the redshift for one year of BNS injections with viewing angle in the range 0 ° < θv < 15° and considering ET (a), ET+CE (b) and ET+2CE (c). The color bar indicates the S/N of each detection. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text.
In the text

	Fig. B.2. Same as Fig. 4, but for Stru2.
In the text

	Fig. B.3. Cumulative distribution of the S/N for BNS mergers localized better than 100 deg2 for the three GW networks considered in this work. Panel (a) shows cases with θv < 15°, while in panel (b) all angles are considered.
In the text

Fig. B.4. Detection efficiency of ET as a function of redshift assuming a Gaussian distribution (⟨M⟩ = 1.33M⊙, σM = 0.1M⊙) and a uniform distribution (M ∈ [1.0 − 2.5]M⊙). Only cases with θv < 15° are considered. A duty cycle of 0.85 has been assumed for the GW detectors as described in the text. The ratio between the orange bars and the blue bars corresponds to the function ϕ(z) defined in the text.

In the text

	Fig. B.5. Corner plot of the posterior distribution of the parameters adopted for the prompt emission model. The assumed structure is Stru1. The number of MCMC steps is chosen such that the auto-correlation time reaches a plateau, as described in Foreman-Mackey et al. (2013).
In the text

	Fig. B.6. Same as Fig. B.5, adopting Stru2.
In the text