© The Authors 2023

1 Introduction

Nearly five years after the first joint observation of a gravitational wave (GW) signal and a short gamma-ray burst (sGRB) from a binary neutron star (BNS) merger (Abbott et al. 2017a,b; Goldstein et al. 2017), multi-messenger astronomy has established itself as a central part of modern astrophysics. The ability to observe astrophysical events both gravitationally and electro-magnetically allows for unique insights into (astro)physical phenomena, potentially probing new fundamental understandings in physics. The implications of multi-messenger astrophysics stretch beyond GW physics, as this concept can be applied to many different areas within physics and astrophysics. For example, one can use joint detections to compute the Hubble constant (for instance Chen et al. 2021a,b), rule out or confirm certain binary progenitor models (for example Eichler et al. 1989), improve our understanding on neutron star (NS) interiors and equation of state (EOS; for example Metzger 2019), or carry out tests of general relativity (for instance Kim 2021; Clark et al. 2015). The number of joint observations will only increase in the near future; BNS and black hole – neutron star (BHNS) mergers, the main candidates for joint detections due to their production of GW signals as well as electromagnetically observable sGRBs, will be observed more frequently. Current detectors such as aLIGO (Harry & the LIGO Scientific Collaboration 2010) and Virgo (Acernese et al. 2015) will eventually operate at their design sensitivities, and the planned Einstein Telescope (ET, Punturo et al. 2010) and Cosmic Explorer (CE, Reitze et al. 2019) will observe GW sources at an unprecedented rate in this frequency band (~ 10–1000 Hz). Additionally, next to already established high-energy (HE) detectors, such as Fermi/GBM (Meegan et al. 2009) and Swift/BAT (Gehrels et al. 2004), newer instruments such as Insight-HXMT (Song et al. 2022) and GECAM-GRD (Xiao et al. 2022; Zhang et al. 2019; launched in 2017 and 2020, respectively) will prove extremely useful in the search of sGRBs as electromagnetic (EM) counterparts to GW observations.

Previously, we built open-source Python-based software, the Gravitational Wave Universe Toolbox (Yi et al. 2022a,b, the GWToolbox or Toolbox hereafter for short) to make GW astrophysics easily accessible. The GWToolbox is available at its own website, or can alternatively be downloaded as a Python package¹. The Toolbox can be used to simulate observations with different kinds of GW detectors on various common GW source populations, among which are BNS and BHNS mergers. It can return a synthetic catalogue of BNSs and BHNSs, which could be detected with a user-specified GW detector, signal-to-noise threshold (S/N), and observation duration. For each event in the catalogue, the Toolbox returns the masses of the binary, luminosity distance, orbital inclination, and effective spin. We have now added² a sGRB model to this that allows the Toolbox to return the properties of a potential sGRB from these sources, as observed by a user-selected HE detector. This enables the user to study the detectability of these events as a GRB with respect to a certain HE detector, and investigate joint GW-EM detection rates. This new functionality makes the GWToolbox a multi-messenger simulator, and paves the way towards more applications in studies such as prospects of joint detections with different GW-HE detectors, sub-threshold strate-gies³, GRB physics, and compact binary merger history. The first study was conducted in this work, while we discuss several other planned projects in Sect. 4 in more detail. While this research is in a sense similar to a recent study by Ronchini et al. (2022), it should be noted that our project touches upon several slightly different aspects of the prospects of multi-messenger astronomy, meaning that these two studies can be viewed as complementary to one another.

In this paper, we describe the sGRB model that has been added to the Toolbox, and investigate joint detection rates predicted by the upgraded GWToolbox. The paper is organised as follows: in Sect. 2, we first give a brief introduction on how the GWToolbox gives the synthetic catalogue of BNSs and BHNSs. Then we describe how GRB emission from a BNS or BHNS progenitor is simulated in our model, and how we determined the detectability of the GRB for a certain HE detector. In Sect. 3, we first reproduced individual and population properties of historical GRBs, and by doing that we fixed the parameters in the underlying population and GRB models. Subsequently, we used the simulation to predict future detection rates, with different combinations of GW and HE detectors, with sub-threshold joint-observation strategies, and specifically for off-axis sGRBs. In Sect. 4 we summarise our findings and provide future prospects.

2 Methods

2.1 The GWToolbox: producing a synthetic GW catalogue

Here we give a brief explanation on how a synthetic GW catalogue is generated in the GWToolbox. For a detailed description, we refer to Yi et al. (2022a). For a certain GW detector, its noise power spectrum and the antenna pattern are denoted as S_n(f) and F_+,× respectively. The S_n(f) is either from literature if the GW detector is from a default list, or calculated with FINESSE (Brown et al. 2020) if it has a user-customised configuration. The detector response to the GW from a compact binary merger is: $$h\left(f\right)=F_{+}{h}_{+}\left(f\right)+F_{\times }{h}_{\times }\left(f\right),$$(1)

where h_+,× are the ‘plus’ and ‘cross’ polarisations of the GW waveform, which are calculated with a certain waveform approximant (IMRPhenomD, for example Husa et al. 2016; Khan et al. 2016). The S/N (also denoted as ρ) of the GW is calculated as: $${\rho }^2=4{\displaystyle {\int }_{f_{\text{low}}}^{f_{\text{high}}}\frac{\left|h^2\left(f\right)\right|}{S_{\text{n}}\left(f\right)}\text{d}f}.$$(2)

For a certain S/N, we can define a probability of detection of a source as: $$\mathcal{D}\left(\Theta \right)={\displaystyle ∯\text{dΩd}{\text{Ω}}^{\prime }\mathscr{H}}\left({\rho }^2-{\rho }_{\star }^2\right)/{\left(4\pi \right)}^2,$$(3)

where Θ denotes the intrinsic parameters (for example the masses and effective spin, depending on the waveform approximant) and the luminosity distance of the source, Ω and Ω′ denote the sky location angles (θ, ϕ), the inclination angle (t) and the polarisation angles (ψ). ℋ is the Heaviside step function. With a merger rate model of the population, ṅ(Θ), the expected distribution of detectable sources is: $$N_{\text{D}}\left(\Theta \right)=\frac{T}{1+z}\frac{\text{d}{V}_{\text{c}}}{\text{d}z}\dot{n}\left(\Theta \right)\mathcal{D}\left(\Theta \right),$$(4)

and the expected number of detections is an integration of N_D(Θ) over the possible parameter space. A synthetic catalogue is obtained by a MCMC sampling from N_D(Θ). Figure 1 is an example of the simulated catalogues from the GWToolbox. The GWToolbox has one default parameterised BNS model (pop-A) and two BHNS population models (pop-A and pop-B). For a detailed description of these population models, we refer to Yi et al. (2022a). The common hyper-parameters for BNS, BHNS-pop-A and pop-B are listed in Table 1 (including hyper-parameters for the GRB model, which are introduced in Sect. 2.3). The synthetic catalogues shown in Fig. 1 were generated under the population model which is consistent with Table 1. For a BNS, there are five extra hyper-parameters describing the mass distribution of the NS, and its effective spin. In our treatment, we assume no dependence of the GRB properties on these parameters. For BHNS-pop-A, there are nine extra parameters describing the mass distributions of the BH and NS, and their effective spin. The values of these hyper-parameters influence the fraction of BHNS system which passes the R_tid > R_ISCO criterion (see Sect. 2.2). We list these hyper-parameters and their default values in Table 2. The difference between BHNS-pop-A and pop-B is that the latter includes another Gaussian mass peak centred at ~40 M_⊙ on the mass function of pop-A to take into account the scenario of the BH originating from pair-instability supernovae.

The catalogue simulated in this way does not include the orbital inclination angles, since the angles are marginalised when calculating the probability of detection. Since the inclination angle t is a crucial property of GRBs, we want to recover t from 𝒟(Θ). We note that 𝒟(Θ)is equivalent to the fraction in a uniform sample of (θ, ϕ, t, ψ) that results in ρ ≥ ρ_⋆. We know: $$\rho \propto \sqrt{{\left(\frac{1+{\cos }^2\iota }{2}\right)}^2{F}_{+}^2+\cos {\iota }^2{F}_{\times }^2,}$$(5)

where we define the right-hand side as $\overline{\rho }$. For a GW event in the catalogue, with a detection probability 𝒟: we draw a sample of uniformly distributed (θ, ϕ, t, ψ), and calculate the corresponding $\overline{\rho }$. We randomly select a $\tilde{\rho }$ which is larger than the (1 – p)th quantile. The corresponding t is assigned as the inclination angle of the source (Schutz 2011). In Fig. 2, we plot the distribution of t of the simulated detected sample of BNSs. We can find in the histogram that for ρ_⋆ = 0, ι follows an isotropic distribution (sin ι). When ρ_⋆ = 8, the distribution of ι peaks at around ~30°. This is a result of the increasing GW signal towards lower inclinations combined with the sin ι of the solid angle.

Fig. 1 Redshift distribution of a simulated BNS merger catalogue from the GWToolbox for a LIGO observation time of 10 yr (blue) and one week of ET observations (orange). The underlying BNS population model is consistent with that in Table 1.

Fig. 2 Inclination angle distribution of GW detection of different S/N threshold.

2.2 BNS and BHNS mergers that result in a sGRB

In general, the merger remnant after a BNS or BHNS is only able to produce a sGRB if it is surrounded by a sizeable accretion disk powering a relativistic jet (for instance Bernuzzi 2020; Fryer et al. 2015; Zhang 2018). For BNS systems, in cases when the binary promptly collapses to a black hole (BH), it is known that the remnant satisfies this criterion. However, depending on the initial binary mass, different scenarios are possible for the final fate of the BNS merger remnant: apart from a BH, a BNS may form a short-lived hyper- or supermassive NS collapsing to a BH, or result in an infinitely stable supermassive NS (for example Ciolfi 2018). There still exists uncertainty regarding the ability of a NS remnant to contain an accretion disk and relativistic jet. However, recent observations as well as numerical evidence suggest that NS remnants may indeed launch jets resulting in a sGRB (Sarin & Lasky 2021). As such, in the Toolbox, we assume that all of the simulated BNSs launch a relativistic jet after merger, triggering a sGRB.

For a BHNS system, only in the case that the NS is swallowed by the BH before it is tidally disrupted, there will be no associated EM counterpart (Fryer et al. 2015). Therefore, we impose a condition for the association of a prompt GRB with a BHNS as: $$R_{\text{tid}}>R_{\text{ISCO}},$$(6)

where R_tid is the tidal radius, within which the NS is tidally disrupted by the BH. It can be calculated with: $$R_{\text{tid}}={\left(\frac{M_{•}}{M_{\text{NS}}}\right)}^{1/3}{R}_{\text{NS}},$$(7)

where M_• and M_NS are masses of the BH and NS respectively, (Fig. 3) and R_NS is the radius of the NS. R_ISCO is the radius of the Innermost Stable Circular Orbit (ISCO) around the BH, within which the NS is thought to be swallowed by the BH. The component masses of the BHNS are outputs from the GW catalogue of the Toolbox. However, the R_NS depends on the mass-radius-relation of a NS. Therefore, in the GWToolbox, we build in four possible mass-radius-relations (MRRs): one relativistic Brueckner-Hartree-Fock EoS (MPA1; Müther et al. 1987), two relativistic mean-field theory EOSs (MS0, MS2), (Müller & Serot 1996), and one nonrelativistic potential model (PAL1; Prakash et al. 1988). We show all MRRs in Fig. 3. In general, any tabulated MRR can be incorporated.

Fig. 3 Built-in MRRs in the GWToolbox. We use the Brueckner-Hartree-Fock EoS (MPA1), two relativistic mean-field theory EOSs (MS0, MS2), and one nonrelativistic potential model (PAL1).

2.3 The sGRB model

With a catalogue of GRB progenitors, we still need a GRB emission model so as to simulate the GRB from those progenitors. Here we employ and adapt a generalised geometrical model for GRBs from Ioka & Nakamura (2001), which is independent of the details of the physics, and is therefore compatible with various physical GRB models. Its mathematical structure allows for potential future adaptations for an increased complexity and accuracy to be easy to implement. Additionally, it serves the purpose of the Toolbox: while it should recover the GRB with a reasonable accuracy, its simplicity keeps the computational cost low allowing for quick computations. The core quantity is the observed flux at the observed instance T: $$F_v\left(T\right)=\frac{vD}{\gamma \beta }{\displaystyle {\int }_0^{2\pi }\varphi {\displaystyle {\int }_0^{{\alpha }_m}{\alpha }^2\text{d}\alpha {\displaystyle {\int }_{v\gamma \left(1-\beta \right)}^{v\gamma \left(1+\beta \right)}\frac{\text{d}{v}^{\prime }}{v^{\prime }}\frac{{j^{\prime }}_{v^{\prime }}}{{\left(1-{\mu }^2\right)}^{3/2}}},}}$$(8)

where γ is the bulk Lorentz factor of the jet, $\beta =\sqrt{1-{\gamma }^2}$; α, α_m and µ are geometry parameters describing the relative position between the emitting region and the observer. ${j^{\prime }}_{v^{\prime }}$ is the emissivity in the jet-comoving frame, which is assumed to be isotropic. For more details, see Ioka & Nakamura (2001) as well as Woods & Loeb (1999) & Piran (1999).

The emissivity ${j^{\prime }}_{v^{\prime }}$ is constructed as follows: $$j_{v^{\prime }}^{\prime }=A_{0}f\left(v^{\prime }\right)\pi \left(t\right)\psi \left(r\right){\displaystyle \prod \left(\text{Ω}\right)},$$(9)

where π(t) and ψ(r) are, respectively, the temporal and radial distribution of the emission, in the comoving frame. Π(Ω) is the angular structure of the jet, A₀ is a constant, and f(v′) denotes the spectrum in the comoving frame. Here, we employ the most simplified case, where a top-hat jet with instantaneous emission at time t₀ and radius r₀ is assumed. We get: $$\pi \left(t\right)=\delta \left(t-t_0\right),$$(10) $$\psi \left(r\right)=\delta \left(r-r_0\right),$$(11) $$\begin{array}{c}{\displaystyle \prod \left(\text{Ω}\right)=H\left(\Delta \theta -\left|\theta -{\theta }_v\right|\right)}\\ H\left[\cos \varphi -\left(\frac{\cos \Delta \theta -\cos {\theta }_v\cos \theta }{\sin {\theta }_v\sin \theta }\right)\right]\end{array}.$$(12)

Δθ and θ_υ are the half-opening angle and viewing angle of the jet, respectively. For the spectrum f(v′), we use the shape of the Band function (Band et al. 1993): $$f\left(v^{\prime }\right)={\left(\frac{v^{\prime }}{v_0^{\prime }}\right)}^{1+{\alpha }_B}{\left(1+\frac{v^{\prime }}{v_0^{\prime }}\right)}^{\left({\alpha }_B-{\beta }_B\right)/s}.$$(13)

We set α_B = –1, β_B = –2.2, and s = 1, as these are typical values resulting in a spectrum whose shape matches the Band spectrum (Preece et al. 2000; Ioka & Nakamura 2018). With the assumptions of Eqs. (10)–(12), Eq. (8) becomes $$F_v\left(T\right)=\frac{2c{A}_0{r}_0{\gamma }^2}{D^2}\frac{\text{Δ}\varphi \left(T\right)f\left(v\gamma \left(1-\beta \cos \theta \left(T\right)\right)\right)}{{\left({\gamma }^2\left(1-\beta \cos \theta \left(T\right)\right)\right)}^2}.$$(14)

We visualise the geometry of our setup in Fig. 4 for an off-axis (top) and on-axis GRB (bottom).

In Eq. (14), the constant A₀ is related to the total energy of the GRB (E_GRB). In order to find the relation between A₀ and E_GRB, we conduct the following derivation. In the comoving frame, the radiating energy of the ejecta is: $$E_{\text{GRB}}^{\prime }={\displaystyle \int j^{\prime }\text{d}V{\displaystyle \int \text{d}}{v}^{\prime }}{\displaystyle \int \text{dΩ\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\displaystyle \int \text{d}t,}}$$(15)

where j′ is the emissivity in the comoving frame, ∫ dV, ∫ dv′, ∫ dΩ and ∫ dt are integrals over volume, frequencies, emission directions and time respectively. Assuming isotropy of emission in the comoving frame, the direction integral gives 4π, and with the assumption of instantaneous emission at an infinitely thin layer with a top-hat jet, the volume and time integrals give $${\displaystyle \int \text{d}V}{\displaystyle \int \text{d}t=2\pi \left(1-\cos \text{\hspace{0.17em}\hspace{0.17em}Δ}\theta \right)r_0^2.}$$

The frequency integral of the Band spectrum gives $${\displaystyle \int f\left(v^{\prime }\right)\text{d}{v}^{\prime }=-A_0{v}_0^{\prime }\frac{s^2\Gamma \left(\frac{\alpha }{2}+1+\frac{2}{s}\right)\Gamma \left(-\frac{\beta }{s}+1-\frac{2}{s}\right)}{\left(\alpha +2\right)\left(\beta +2\right)\left|s\right|\Gamma \left(-\frac{\beta -\alpha }{2}\right)}.}$$

Therefore, Eq. (15) becomes: $$E_{\text{GRB}}^{\prime }=8{\pi }^2\left(1-\cos \text{\hspace{0.17em}\hspace{0.17em}Δ}\theta \right)r_0^2{A}_0ℐv_0^{\prime },$$(16)

where $$ℐ=-\frac{s^2\Gamma \left(\frac{\alpha }{2}+1+\frac{2}{s}\right)\Gamma \left(-\frac{\beta }{s}+1-\frac{2}{s}\right)}{\left(\alpha +2\right)\left(\beta +2\right)\left|s\right|\Gamma \left(-\frac{\beta -\alpha }{2}\right)},$$

which is 6.64 when s = 1, α = −1, and β = −2.2. Note that energy is not Lorentz invariant, and in the observer frame, the energy transformed to: $$E_{\text{GRB}}=\gamma E_{\text{GRB}}^{\prime }.$$(17)

So, we have $$E_{\text{GRB}}=A_0\gamma 8{\pi }^2\left(1-\cos \text{\hspace{0.17em}\hspace{0.17em}Δ}\theta \right)r_0^2{A}_0ℐv_0^{\prime }.$$(18)

Equivalently: $$A_0=\frac{E_{\text{GRB}}}{\gamma 8{\pi }^2\left(1-\cos \text{\hspace{0.17em}\hspace{0.17em}Δ}\theta \right)r_0^2{A}_0ℐv_0^{\prime }}.$$(19)

It should be noted that E_GRB depends on the mass of the accretion disk surrounding the merger remnant as well as the accretion efficiency, both of which depend on the binary parameters (Bernuzzi 2020; Fryer et al. 2015). Where efforts have been made to connect these quantities (for example Radice et al. 2018; Bernuzzi et al. 2020; Nedora et al. 2022), as of yet the exact relationship between the binary properties and GRB energy is not known. As such, taking E_GRB as a direct free parameter instead of indirectly through, for instance, the binary mass is currently the safest option. More so, the total GRB energy is typically well-approximated by assuming a standard reservoir clustered around some central energy value (Frail et al. 2001; Zhang 2018).

Note that Eq. (14) assumes an instantaneous emission from one location r = r₀ and t = r₀/c. The corresponding lightcurve is a single pulse which rises sharply at T = 0 for on-axis cases. In reality the lightcurve is a superposition of multiple pulses, resulting in a more gradual rise. We recover this feature by implementing a temporally extended jet launch: every GRB pulse is launched with a different delay time t_d (in the source frame) relative to the merger, and t_d follows an extended probability distribution. The lightcurve can be obtained by convolving the original lightcurve with the probability distribution of (the red-shifted) t_d, for which we use a Gaussian with a width t_d,std(1 + z), centred at 3t_d,std(1 + z).

As a summary, our GRB requires the following parameters. Firstly, we list the intrinsic parameters, consisting of γ (the bulk Lorentz factor of the jet), r₀ (the radius from the central engine, where the gamma-ray emission happens), E_GRB (the total energy released into prompt gamma-ray, in the observer frame, which is related with the normalisation parameter A₀), α_B, β_B, s and ${v^{\prime }}_0$ (the intrinsic spectrum parameters), Δθ (the half opening angle of the jet cone, and t_d,Std (the standard deviation of the intrinsic jet launch delay time distribution). Secondly, there are also several extrinsic parameters that are required for the model: ι (the viewing angle of the observer, or equivalently, the inclination angle of the binary orbital plane), and D or z (the luminosity distance from the source to the observer, or equivalently, the redshift of the GRB).

For each GRB in the Toolbox, the Δθ, log r₀, log E_grb, ${v^{\prime }}_0$, γ are assigned as Gaussian randoms. The corresponding means and standard deviations are Δθ_mean and Δθ_std, log r_0,mean and log r_0,std, log E_grb,mean and log E_grb,std, ${v^{\prime }}_{0,\text{mean}}$ and ${v^{\prime }}_{0,\text{std}}$, γ_mean and γ_std respectively. We summarise these hyper-parameters for the underlying GRB population in Table 1. The inclination angle ι and the distance D are taken directly from the GW part of the GWToolbox.

To illustrate and visualise the model, in Figs. 5 and 6 we show the lightcurve of an on- and off-axis sGRB respectively. The lightcurve, either in units of phs s⁻¹ cm⁻² or ergs s⁻¹ cm⁻², can be calculated with: $$L_{\text{ph}}\left(t\right)={\displaystyle \int \frac{F_v\left(t\right)}{hv}\text{d}v,}$$(20)

or $$L_{\text{ergs}}\left(t\right)={\displaystyle \int F_v\left(t\right)\text{d}v},$$(21)

respectively, where the integral over the photon frequencies corresponds to the energy range of the detector. In Fig. 5, we attempt to recover the properties of the on-axis sGRB GRB131004A, observed by Fermi/GBM (von Kienlin et al. 2014, 2020; Gruber et al. 2014; Bhat et al. 2016). With the hyper-parameters listed in Fig. 5, our model finds: T₉₀ = 0.94 s (compared to an observed 1.152 s), a fluence (10–1000 keV) of 5.35 × 10⁻⁷ (observed 5.099 × 10⁻⁷) ergs cm⁻², a 1-s and 64-ms peak flux of respectively 5.72 (6.8) and 7.96 (9.8) ph cm⁻² s⁻¹, and E_peak = 122 (observed 118) keV. Overall, these simulated observables match their detected counterparts well for the purposes of the GWToolbox. Where the general shape of the lightcurve in Fig. 5 is recovered well, an estimation of the background photon flux in the observation highlights some discrepancies; for the observations, we find a typical background of ~8 ph cm⁻² s⁻¹, while our modelled background is 0.456ph cm⁻² s⁻¹ for Fermi/GBM. This has a direct influence on the GRB detection rate, but since we calibrate the GWToolbox detection rate to observations, this is not problematic. However, future revisions of the Toolbox might need a more sophisticated background model.

We also show the lightcurve according to our model in the same energy range, for the off-axis GRB170817A observed by Fermi/GBM (Goldstein et al. 2017) in Figs. 5 and 6. Here, we find T₉₀ = 2.41s (compared to an observed 2.0 ± 0.5s), a fluence of 3.69 × 10⁻⁷ (observed 2.8 ± 0.2 × 10⁻⁷) ergs cm⁻², a 64-ms peak flux of 3.53 (compared to 3.7 ± 0.9) ph cm⁻² s⁻¹, and an E_peak of 192 keV (observed 185 ± 62) keV. These findings all match the observed values quite well. Also here, the background signal is slightly off, which is unproblematic for the same reasons as the GRB131004A spectrum. In the case of GRB170817A, it is important to also compare with the observed delay time between the peak of the GW signal and the start of the GRB. The observed time between the starting time⁴) of the sGRB (defined as the time corresponding to 10% of the peak) is 1.424±0.102s, whereas we find 0.85 s. It appears that our prediction comes close, but still slightly underestimates the delay time with these choices for the hyper-parameters. This discrepancy stems from our efforts to obtain a somewhat realistic E_p, which is a direct result of our choice for ${v^{\prime }}_0$, Δθ-θ_υ. Where the E_p that we found corresponds to observations, it came with the cost of an atypical ${v^{\prime }}_0$ and a slightly too low delay time, the latter of which is caused by the difference between Δθ and θ_υ. This strongly suggests that our top-hat jet model may contain several limitations which become exposed for off-axis GRBs, specifically regarding the angular dependence of the GRB spectrum; for a more detailed discussion, see Sect. 4 and Ioka & Nakamura (2019). Nonetheless, keeping in mind that we expect the majority of GRB observations to be on-axis, if its main limitation is an inaccurate computation of off-axis spectral peaks, our model is accurate enough to serve the purpose of the GWToolbox. The aim of the Toolbox is not to perfectly match all observations made by HE detectors; it should simply be able to recover the majority of the observations reasonably well, with low computational effort. This allows the user to make quick computations of joint detectabilities.

Fig. 4 Schematic overview of the structure and geometry of our sGRB model, for a jet with opening angle Δθ, viewing angle θυ and distance to earth D. We depict an off-axis GRB in the top panel, while the bottom shows an on-axis one.

Fig. 5 Lightcurve of a GRB131004A-like on-axis GRB, according to the GWToolbox. We use the following hyper-parameters: γ = 120, log(r0/cm) = 12.95, log(EGRB/ergs) = 50.1, ${v^{\prime }}_0=0.7$ keV, Δθ = 15°, θυ = 10°, and z = 0.717. In blue, we show the simulated lightcurve of the on-axis GRB131004A, within an energy range of 10–1000 keV. The red data points correspond to the measured photon count rate by Fermi/GBM.

Fig. 6 Lightcurve of a GRB 170817A-like off-axis GRB, according to the GWToolbox. We use γ = 120, log(r0/cm) = 13.03, log(EGRB/ergs) = 49.65, ${v^{\prime }}_0=8$ keV, Δθ = 15°, θυ = 18.60°, and z = 0.0099. In blue, we show the simulated lightcurve of the on-axis GRB170817A, within an energy range of 10–1000 keV. The red data points correspond to the measured photon count rate by Fermi/GBM.

Fig. 7 Flowchart showing the sGRB detection test pipeline of BNSs.

Fig. 8 Flowchart showing the sGRB detection test pipeline of BHNSs.

2.4 GRB detectability against HE detectors

To decide whether a simulated sGRB is detectable by a selected HE instrument, the Toolbox computes the peak of the lightcurve (using Eqs. (20) or (21)), and adds to this peak some typical background flux within the specified detector energy range. A GRB is detected if the sum of the GRB flux and background flux is above the detection threshold of the HE instrument. For Swift/BAT sGRBs, a typical background is 0.3 ph s⁻¹ cm⁻² (Barthelmy et al. 2005), which, with a flux threshold of 1.5 ph s⁻¹ cm⁻² (Coward et al. 2012), leads to a detection if the total flux is >5 times that of the background. For simplicity, we assume the same factor of 5 for other detectors, and since we know the effective threshold of these other HE instruments, we can set the background accordingly. Even if a GRB is bright enough to be observed by a certain detector, it may still remain undetected, for either it is not in the detector’s field-of-view (FoV), or it does not pass the detector’s complex trigger algorithm. We treat this by introducing two parameters for each detector: (i) the FoV as a percentage of the full sky, and (ii) a factor f_c to account for the correction due to the trigger algorithm. For each GRB that passes the sensitivity threshold, we compare a random variable u in [0,1] with FoV/f_c. If u < FoV/f_C, we consider the GRB detected, otherwise we flag it undetected. In Figs. 7 and 8, we summarise the detection algorithm for a BNS and BHNS respectively. In Table 3, we list the parameters for our built-in HE detectors.

3 Results and discussion

3.1 Comparing the GRB model with observations

Before investigating the GWToolbox GRB observations, we comment on our population parameters. In the GWToolbox we use R₀ = 50 Gpc⁻³ yr⁻¹ and τ = 10 Gyr for the BNS population, yielding a local (z = 0) merger rate density of ~125 Gpc⁻³ yr⁻¹. This is in accordance with estimates from GW observations, of which the most recent findings give 13–1900 Gpc⁻³ yr⁻¹ (The LIGO Scientific Collaboration 2022). While our estimate sits in the lower end of the predicted range and more representative choices for R₀ and τ may exist, we note that these choices (in combination with our choice for Δθ which is based on observations in Fong et al. 2015) yield an expected on-axis sGRB rate (that is whose jet is directly in our line-of-sight) of ~11 GRBs Gpc⁻³ yr⁻¹. This agrees with a previously predicted rate of $8_{-3}^{+5}$ Gpc⁻³ yr⁻¹ (Coward et al. 2012). In this study, we assume that all generated sGRBs originate from BNSs. This assumption is justified as follows: the Toolbox finds that the GRB-launching BHNS population (that is satisfying Eq. (6)), for which we use R₀ = 10 Gpc⁻³ yr⁻¹ and τ = 10 Gyr to match observational estimates (Abbott et al. 2021) in the universe (~2300 in 1 yr) is about 17 times smaller than its BNS equivalent (~40 000 per year). It is therefore negligible for our purposes. For a LIGO-O3 configuration, the GWToolbox finds a GRB detection rate of ~0.68 per year, which we deem reasonably accurate as compared to the one BNS that was detected during O3 (Abbott et al. 2020a).

To verify how realistically the GWToolbox can recover the observed GRB population, we let it observe BNSs for 1 yr with ET, using a S/N threshold of 0. This is an approximation to all BNS mergers in the universe within 1 yr. Due to computational cost, we choose not to use a longer observing time. While this could lead to a potentially low number simulated sGRBs to compare with observations, we regard this as unproblematic as our goal is to only recover observations reasonably well. We constructed a detected GRB population from this by observing the BNSs with Swift/BAT. In Fig. 9, we show in blue several observables of this population. We compare these with the full observed Swift/BAT GRB catalogue for GRBs with T₉₀ ≤ 2 s (Lien et al. 2016).

First of all, we find that apart from the intrinsic GRB rate matching with Coward et al. (2012), we are also able to successfully reproduce actual observation rates for all detectors. In the case of Swift/BAT depicted in Fig. 9, we find ~8 sGRBs per year, which is in line with the observed rate of 8 sGRB yr⁻¹ (Lien et al. 2016). It should be noted that we manually chose our correction factor f_c such that the resulting detection rate would (roughly) coincide with observations. Therefore, this comparison should be seen as nothing more than a check that this calibration was successful. Concerning other observables, it is especially important that we are able to recover the observed redshift distribution to at least some degree of accuracy, as the predictions for joint detection rates that we make further along in this paper heavily rely on an accurate sGRB redshift distribution. Looking at Fig. 9, we observe that in fact most histograms (redshift, fluence, and E_p) show a good match. Only the histograms for the duration T₉₀ and 1-second peak flux seem to indicate that a part of the observed population cannot be reproduced by the GWToolbox. In the case of T₉₀, sGRBs shorter than ~0.3 seconds are not generated by the Toolbox. The cause of this is our manually implemented time delay between merger and jet launch, quantified by the parameter t_d,std. Since we set this standard deviation as well as the corresponding mean to a fixed number (0.4 × (1 + z)s and 1.2 × (1 + z)s respectively), there cannot be significant variations in T₉₀. Still, since the duration is not nearly as decisive for the detectability of the GRB as, for example, the fluence or the redshift, this is not problematic for our purposes. Somewhat similarly, the 1-s peak photon flux shows no detections below ~1 ph cm⁻² s⁻¹. This is a direct result of our choice to approximate the complex Swift/BAT triggering algorithm by setting a relatively high threshold (a 20-ms peak photon flux of 1.5 phs cm⁻² s⁻¹) and introducing a correction factor f_c. Where a more accurate representation would result in the inclusion of more low-flux (and low-fluence) sGRBs, this would be too complex for the (current version of) the GWToolbox. Nonetheless, the accuracy of our recovered redshift distribution clearly does not visibly suffer from the lack of low-energy sGRBs so for the purposes of this paper, we also deem this limitation acceptable.

3.2 Joint detection rates

In Table 4, we show the number of GW detections per year for LIGO (design sensitivity) and ET, as well as joint GW-sGRB detections for four different HE detectors. We give numbers for three different S/N thresholds for the GW detection. In this section, we only discuss observations with the standard S/N of 8, whereas sub-threshold detections with S/N 7 or 6 are addressed in Sect. 3.3. While expected joint rates for ET are in line with other predictions (for example Ronchini et al. 2022), our predictions for LIGO are lower than those in other studies (for instance Howell et al. 2019; Clark et al. 2015). This is due to the low number of off-axis GRB observations according to our model. We further discuss this in Sect. 3.4. In Figs. 10 and 11, we visualise the redshift distribution of our predicted GW and sGRB observations, for LIGO (design sensitivity) and ET respectively. These figures do not only generally provide insights into joint detection rates, but also give information on sub-threshold and off-axis joint observations. Figures 10–13 all use the same Fermi/GBM GRB distribution. In this section, we discuss the former topic of these three, whereas he latter two are addressed in Sects. 3.3 and 3.4.

A feat that immediately stands out in Table 4 is the substantial difference in yearly joint detections between LIGO and ET. This is visualised in Figs. 10 and 11. To accompany these, in Fig. 12 we also show the redshift distribution of Fermi/GBM with LIGO-O3. In these figures, joint detections can only be made in the overlapping region. We see that the LIGO horizon is a main cause of the small number of joint GW-EM observations in our current era. Typical observable sGRBs by Fermi/GBM lie at redshifts ~0.6–0.8, while even at design sensitivity LIGO only reaches up to redshift ~0.l. For LIGO-O3, the overlap is very minimal and hardly visible. For this configuration, we expect ~0.027 joint detections per year out of ~0.681 GW detections per year. This is a surprisingly low number of joint observations, given that we in fact had one joint detection in O2 (Abbott et al. 2017a). We refer to Sect. 3.4 for an explanation. Even though LIGO-design shows an increase in joint detections with Fermi/GBM (~0.141 yr⁻¹) over LIGO-O3, drastic instrument changes would have to be made in order to significantly extend the horizon out to where more BNSs can be observed gravitationally.

Our predictions for ET show the effects of such changes, depicted in Fig. 11. Not only does ET observe significantly more sources than LIGO up to a redshift of 0.1, its horizon stretches out to much larger redshifts, clearly increasing the number of joint observations. We find that ET makes sure that the GW coverage for each HE detector in our model is close to 100%, that is almost every observed GRB will also have a detected GW counterpart, as the joint observation rates found in Table 4 are very close to the yearly detection rates by the HE instruments themselves. This is an important finding, as it will allow us to place constraints on sGRB progenitors: if all sGRBs are accompanied by a GW signal, we will know that compact binaries are the sole origins of these bursts. It is important to realise that this is a theoretical prediction, as the GW detectors do not have a 100% duty cycle, so in reality the GW coverage for every HE instument will be lower. Additionally, it should be noted that the while HE detectors that we use in the GWToolbox will not be operational during the ET era, we can still use them as reference instruments to extract useful information, some of which will apply to future instruments as well. With future HE detectors such as AMEGO (McEnery et al. 2019), HERMES (Fiore et al. 2020), TAP-GTM (Camp et al. 2019) and THESEUS (Amati et al. 2018, 2021; Ciolfi et al. 2021; Rosati et al. 2021), we will be able to make joint observations at redshifts far beyond the current Fermi/GBM limit of ~2.

For illustrative purposes, we also show a ‘half-ET’ configuration in of Fig. 13. This detector has a sensitivity that is halfway in-between LIGO (design sensitivity) and ET, which we obtained by multiplying the ET noise curve with a factor 10.89 (see Yi et al. 2022b), representing an early stage of the ET capability or GW detectors in-between LIGO-design and ET-design. This configuration yields ~408 GW observations per year, ~6 of which would be observed by Fermi/GBM as well. This shows that even an ET operating at a fraction of its capacity would be a major improvement over current-era detectors. As it is expected that ET will not immediately operate at its designed sensitivity, this is a very bright prospect. This is especially the case for the number of joint detections, as the ‘half-ET’ horizon lies at a red-shift of ~ 1, which is much further than the typical LIGO horizon, leading to more and higher-redshift joint detections.

Fig. 9 Fluence, 1-second peak photon flux, T90 and redshift histograms. We compare the GWToolbox simulation of Swift/BAT observing sGRBs for 1 yr (blue) with all sGRBs in the Swift/BAT GRB catalogue (red; Lien et al. 2016).

Fig. 10 Redshift distribution of BNSs through detection of GWs and sGRBs. In blue, we show the distribution for LIGO (design sensitivity), where the shaded region with the dotted line represents the distribution for those GW sources with S/N between 6 and 8. In orange, we plot the distribution of Fermi/GBM GRBs (the same distribution as in Figs. 11–13), where the red region corresponds to the off-axis GRBs, generated with the addition to our model as described in Sect. 3.4.

Fig. 11 Redshift distribution of BNSs through detection of GWs and sGRBs. In blue, we show the distribution for ET. Again, in orange, we plot the distribution of Fermi/GBM GRBs (same distribution as Figs. 10, 12 and 13), where the red region corresponds to the off-axis GRBs, generated with the addition to our model as described in Sect. 3.4. The orange region with dotted line represents those GRBs with a peak flux between 50 and 100% of the Fermi/GBM detection threshold.

Fig. 12 Redshift distribution of LIGO-O3 BNSs (blue), and Fermi/GBM BNSs measured through sGRBs (orange). Once again, the Fermi/GBM distribution is the same as in Figs. 10, 11 and 13.

Fig. 13 Redshift distribution of BNSs detected with a ‘half-ET’ configuration (blue) which has a sensitivity halfway between LIGO and ET, and Fermi/GBM BNSs measured through sGRBs (orange, same distribution as in Figs. 10–12). Also here, the shaded region with the dotted line corresponds to those sources with a S/N between 6 and 8. S/N8: 2038, S/N6: 6019 in 5 yr.

3.3 Sub-threshold detections

In the multi-messenger astronomy era, we have the possibility of using GW observations to search for possible EM counterparts in data from HE instruments, as well as the other way around. Since this might lead to the discovery of sGRBs or GW signals that would otherwise not have been found, a quantitative investigation could provide us valuable insights. Here, we address the potential gains of a targeted search for sub-threshold observations of BNSs.

Table 4 depicts the number of GW and joint observations for different GW-HE combinations as predicted by the Toolbox, for S/N thresholds of 8, 7 and 6 respectively. With the Toolbox, we predict that for a threshold of 6, about 35 of 100 detected events is of astrophysical origin⁵. With HE data at hand, it would be quite easy to filter out the real GW signals whose properties match with those of the detected GRBs. We find that the fraction of astrophysical events drops significantly beyond 6; at 5.5, we already find ourselves at about 1 out of 100 detections that has an astrophysical origin. However, we might be able to still pinpoint this single GW detection with sGRB data; it is generally up to the GW-EM community which lowest S/N value maximises the scientific output of such studies Lynch et al. (2018). Table 5 shows a type of results similar to Table 4, but for sub-threshold sGRB detections. Here, we do not go down further than a fraction 0.3 of the current treshold, as the Toolbox models the typical flux of the background as a fraction 0.2 of the threshold, that is we need to stay above the background.

Table 4 clearly shows that the LIGO sensitivity is an important limiting factor in the current era. Figure 10 visualises this, as we observe that many of the sources between S/N = 8 and S/N = 6 are at higher redshift (≥0.1, up to ~0.15), beyond the LIGO horizon with S/N = 8. In this case, decreasing the threshold shifts the LIGO horizon more towards the peak of the z-distribution of the HE detectors, thus increasing the number of joint detections. Decreasing the S/N to 6 would lead to an increase in joint detections of, for example, 179 and 185% for Fermi/GBM and GECAM with LIGO, respectively. As such, targeted searches for such sub-threshold GW observations would be extremely beneficial not only for the number of joint detections, but also for the BNS merger history. As the weaker GW signals typically originate from those BNSs further away, this would allow us to investigate the BNS redshift distribution and merger history in more detail. Additionally, with a larger sample extending to higher redshifts, we would be able to test gravitational waveform models in a larger part of their parameter space. Interestingly, lowering the threshold for HE detectors does not yield substantially more joint observations, as we notice in Table 5. While the number of sGRB detections increases for lower thresholds, the majority of them is beyond the LIGO horizon and therefore does not contribute to the joint detection rate. Only a few intrinsically faint sGRBs below redshift ~0.1 are added to the number of joint detections. Therefore, we can conclude that searches for sub-threshold GRBs using LIGO detections should not be a priority. The yield of such campaigns would be very limited, and focus should be on sub-threshold GWs with information from sGRB data at hand instead.

Tables 4 and 5 show different patterns for the ET era. First of all, in Table 4 we see barely any increase in the number of joint detections as we lower the threshold and increase the number of GW observations by ET. This ties in well with our previous remarks on the fact that nearly 100% of all GRBs in the ET era will also be detected gravitationally. The depletion of current-era HE detectors during the ET era that is clearly visible here can be interpreted as a motivation for the development of new HE detectors that can keep up with the ET sensitivity. Future detectors should cover at least part of redshifts ~2–5 to make a significant contribution to the number of joint detections with ET. Still, without new detectors reaching to much larger redshifts, a search for sub-threshold observations in GRB data of instuments like the current ones during the ET era could yield many new joint detections. In Table 5, we can see increases up to 270% for Fermi/GBM, and joint observation numbers as high as 227 per year for GECAM. Similarly to before, lowering the HE detection threshold leads to more GRBs with higher redshift, most of which have a GW counterpart. Sub-threshold searches should not only be focused on higher redshifts (beyond ~2), but also on those BNSs with low inclination. Those GRBs that are further away are much more likely to be observable if their jet points in our line-of-sight; off-axis GRBs are too weak to be observed at such redshifts. BNSs with inclination angles below ~15º should be prioritised, as this is the typical width of a jet opening angle (Sarin et al. 2022; Fong et al. 2015; Salafia & Ghirlanda 2022).

In Fig. 13, the shaded region around the main blue population visualises the BNS detections by the ‘half-ET’ configuration that have a S/N between 6 and 8. This population contains ~956 BNSs per year, which results in a total number of 1364 when added to the original 408 sources with S/N ≥ 8. Most sub-threshold events are contained at redshifts beyond 0.2–0.3, reaching out to ≥1. The total population of BNSs with an S/N of 6 or higher yields ~21 joint detections with half-ET and Fermi/GBM. This is clearly higher than the ~6 joint detections for those with S/N ≥ 8. Where the redshift distribution of half-ET BNSs does not overlap completely with Fermi/GBM sGRBs, including sub-threshold GW detections can partially take care of this limitation as Fig. 13 shows us that many BNSs are added at the high-redshift end. This results in about 50% of Fermi/GBM detections that have a GW counterpart which is detectable by half-ET. Where this is clearly not close to the near 100% predicted for ET, it shows that even a half-ET configuration can be extremely beneficial and valuable for the GW-EM community, especially when compared with current-era GW instruments.

3.4 Off-axis detections

The sole joint GW-GRB detection thus far contained an off-axis GRB (Abbott et al. 2017a; Goldstein et al. 2017; Ioka & Nakamura 2018). With the Toolbox we can predict the expected detection rate of such GRBs. We report the following results: we find ~0.0035 joint observations with LIGO-O3 and Fermi/GBM per year that are off-axis (with a general joint detection rate of ~0.027 yr⁻¹). Assuming this number is realistic, chances that the very first joint GW-sGRB observation contains an off-axis GRB are very slim, making GW170817/GRB170817A an extremely unlikely event. A more logical explanation is the fact that these numbers expose the limitations of the top-hat jet model. In reality, the energy of the GRB leaks to a wider angle than we describe in this model, leading to observable (off-axis) GRBs whose inclination is well outside the jet cone. As such, with our current model, we treat many GRBs with higher inclination angles as undetected, where in reality they are visible to HE instruments.

In this Section, we improve our off-axis detection rate by adding a second component to the model. Typically, a jet is surrounded by a wide-angle, mildly relativistic cocoon surrounding the relativistic jet (Zhang 2018; Ciolfi 2018). We model this cocoon as the main source of off-axis detections. Still, since the vast majority of independent HE detections of sGRBs consists of on-axis GRBs, the cocoon should be several orders of magnitude less energetic than the main jet, not significantly contributing to the HE instruments’ detection rate. We model the wide-angle cocoon emission by running a second population of GRBs, next to our main population described in Table 1. For this second population, we use log E_{GRB, MEAN} = 48, γ = 150 and a typical Δθ_mean = 60º with Δθ_std = 1º (see for example Ramirez-Ruiz et al. 2002). All other parameters remain unchanged. With this population, we can still recover the properties of GRB170817A-like GRBs, this time without having to increase ${v^{\prime }}_0$ to atypical values.

We show our predicted off-axis detection rates in Table 6. This population satisfies the condition of not significantly increasing the number of independent HE detections as compared to the main population: we find, for example, one to two extra detections with Fermi/GBM per year. The joint detection rates we obtain with LIGO are more in line with predictions made by previous studies (for instance Howell et al. 2019; Clark et al. 2015). Additionally, we highlight a recent study by Colombo et al. (2022) due to its similarities with ours; Colombo et al. (2022) report an O4 (see also design sensitivity, Abbott et al. 2020b), and a joint detection rate for aLIGO-Fermi/GBM and aLIGO-Swif/BAT of ${0.17}_{-0.13}^{+0.26}$ and ${0.03}_{-0.02}^{+0.04}$ yr⁻¹ respectively. For both detector combinations we find similar estimates, albeit ours are slightly higher. Our sub-threshold predictions show the same behaviour, although here the two studies match one another better: we find, for example, (including the off-axis population) ~0.953 joint aLIGO-Fermi/GBM observations per year with S/N ≥6, against ${0.75}_{-0.55}^{+1.16}$ yr⁻¹ in Colombo et al. (2022). Differences can largely be attributed to different choices for the local BNS merger density in the population model, and to the model setup for the GRB emission.

Interestingly, we find that for all HE detectors, the joint off-axis detection rate with LIGO exceeds the number on-axis GRB observations listed in Table 4. Here, we can clearly see the impact of the inclination angle and redshift. With a typical cocoon opening angle about four times larger than the GRB jet opening angle, this second population generally contains more sources whose emission is in the line-of-sight of the detector. Moreover, even though their flux is relatively weak as compared to the main GRB jet itself, many of these GRBs are still detected due to their proximity. With the LIGO horizon at a red-shift of only ~0.1, a relatively high off-axis detection rate is not surprising. This shows that the properties and characteristics of GW170817/GRB170817A are not as unique as one might think at first glance; with this improvement included, an off-axis GRB as a first joint detection is in fact quite likely. For ET, the off-axis detections encompass only a small part of the total number of joint observations. When combined with the on-axis observations in Table 4, we still find a total rate that is in agreement with predictions (Ronchini et al. 2022). ET stretches out to a horizon (a redshift of ~5) that is far beyond that of HE instruments observing off-axis GRBs (only ~0.05 for Fermi/GBM). As such, the majority of ET joint observations will contain on-axis GRBs. Still, we expect more low-redshift off-axis joint detections with ET as compared to LIGO, due to its higher sensitivity. This means that the ET era is a very bright prospect, also in this regard. A potential ~20 joint off-axis detections per year with ET-GECAM and likely more with next-generation HE detectors will allow us to cover significant parts of the parameter space for off-axis GRBs. We will be able to understand the GRB emission mechanism in more detail by comparing with the high number of on-axis joint detections we will have. This may probe a higher level of understanding of GRB central engines.

4 Summary and future prospects

Based on GWToolbox, we developed an algorithm that simulates the potential prompt GRB emission from BNS and BHNS mergers, and checked their detectability against a list of known HE detectors. Our algorithm pipeline of sGRB detectability determination is plotted as flowcharts in Figs. 7 and 8. We employed the algorithm to simulate observations of sGRBs and joint observations with GW and HE detectors of BNS and BHNS, and studied the prospect of sub-threshold observation strategies in both GW and HE observation.

In the list below, we summarise our main findings in this study. Afterwards, we discuss several future prospects for additions to the model as well as possible topics to explore using the GWToolbox.

• We have constructed a sGRB model for the GWToolbox using a top-hat jet, operating with an accuracy sufficient for the purposes of the Toolbox. The addition of this model increases the number of possible uses for the Toolbox making it a valuable tool to investigate the detectability and properties of GW sources in the multi-messenger astronomy era.

• Our findings stress the significance of the impact that ET will have on multi-messenger astronomy. Where currently the LIGO sensitivity is the main limiting factor to the relatively low number of observations, we can quantitatively see that ET completely turns this around and will detect BNSs at an unprecedented rate. We have showed that a ‘half-ET’ instrument which is more than ten times less sensitive than ET already increases the number of joint detections significantly.

• We find that in the ET era, the GW coverage will be close to 100% for each HE detector. Almost every detected GRB will have a GW counterpart that is observed by ET. In reality, this percentage will be lower due to the fact that GW detectors do not have a 100% duty cycle. Still, reaching 100% coverage would allow us to confirm that compact binaries are the sole origin of sGRBs.

• A targeted search campaign in LIGO data for sub-threshold GW detections has the potential of being successful in increasing the number of joint observations. LIGO is not able to detect most independent GRB detections by HE instruments. With information about the source from GRB data at hand, it will be relatively simple to pinpoint GW signals within the LIGO datastream in terms of localisation, chirp time and parameter space. We show an increase in joint detections of, for example, ~185% for GECAM if we would be able to find GW detections down to an S/N of 6. This highlights the value of multi-messenger observations, as a sub-threshold search campaign will clearly increase the number of observed BNSs, as well as the parameter space in which these BNSs find themselves.

• In the ET era, this will be the other way around as we will have many more GW observations than GRB detections. Those GWs further away measured by ET can be used to search for sub-threshold GRBs. Since ET stretches quite far out, we will be able to find more (for instance ~170% more if we search down to a threshold half of the current Fermi/GBM threshold) joint detections if we perform a targeted search to sub-threshold GRBs. Such a search should be focused on high-redshift (for example beyond ~2 for Fermi/GBM) GRBs, as closeby BNSs will be covered by both HE and GW detectors. Additionally, such campaigns are more likely to be successful when focusing on BNSs with inclinations lower than ~15º, as this is the typical opening angle of a GRB jet.

• The top-hat jet model in the GWToolbox generally predicts joint detection rates that are quite low compared to other predictions, especially with LIGO. This is a direct result of the low number of off-axis detections predicted by the model. We achieve more realistic rates by including a second population, modelling the wide-angle cocoon around the jet, showing that the majority of joint detections in the LIGO era originates from off-axis GRBs. This feat makes the characteristics of GRB170817A not as unlikely as one may think based on the top-hat model. With ET, ~ 100% of electromagnetically detected off-axis GRBs will have an GW counterpart. However, the ET horizon stretches far beyond the off-axis detection horizon of current HE instruments, so on-axis GRBs will dominate the ET joint detection rate.

Our model has a number of limitations that can be addressed in future work. Here, we list potential improvements that serve as future prospects for continuation of this work. Additionally, we address general applications of the newly upgraded GWToolbox now functioning as a multi-messenger simulation programme.

4.1 Future improvements

• We employed the top-hat jet model in this study, and a nested wider and less energetic top-hat jet to mimic the off-axis emissions. With increasing evidences, people now believe the jet of GRB should be structured (Salafia & Ghirlanda 2022). Our methodology is compatible with a structured jet, by substituting the top-hat function Π(Ω) in the emis-sivity with a more general continuous angular distribution. The time-evolving spectrum in this case will not take the form in Eq. (14), but should be integrated through Eq. (8) respectively.

• The comoving spectrum indices α_B, β_B and s are not sampled from a distribution like γ, ${v^{\prime }}_0$ etc., but share the values for all sGRBs in the catalogue. In reality these parameters can differ from burst to burst, depending on the local emission condition and mechanism. We argue that the influence on detectability of bursts from the diversity on these spectrum parameters is largely degenerate with that from γ and ${v^{\prime }}_0$. Our current treatment is a trade-off between the generality and degrees-of-freedom of the model, although extension to include six more hyper-parameters to describing distributions of these spectrum parameters would be straightforward. Other possible improvements to the spectrum could be more along the lines of Ioka & Nakamura (2019), where we would implement an intrinsic angular dependence in both γ and ${v^{\prime }}_0$.

• In our simulation, we assume every BNS merger or BHNS that have R_t > R_ISCO is associated with prompt GRB emission. While there are studies show that some merger remnants can not work as engine to power a GRB jet (Sarin & Lasky 2021; Shemi & Piran 1990; Ciolfi 2018; Nakar 2007; Murguia-Berthier et al. 2014, 2017; Margalit et al. 2015; Beniamini et al. 2017), and the scenario of a ’choked jet’ is also possible; here, the jet fails to break out of the surrounding ejecta, and is not able to dissipate and emit a prompt GRB (Rosswog & Davies 2002; Murguia-Berthier et al. 2014; Lazzati et al. 2018; Bromberg et al. 2011; Lazzati & Perna 2019; Salafia et al. 2020). The above-mentioned possibility will yield a fraction of failed GRBs, which is not considered in our treatment, and will cause an over-estimation of the sGRB detection rate by a factor f_fail. On the other hand, there could be sGRB whose progenitor is anything other than a BNS or BHNS merger, which we also did not include in our consideration, which will cause an underestimation of the detection rate by a factor f_other. Since we collaborate the simulator with the historical detection rates, those effect of f_otherf_fail is absorbed into the merger rate parameter R_n and the detector correcting factor f_c, and therefore will not cause under- or overestimation on the detection rate of future (joint) observations. When R_n is more accurately estimated from a larger GW population, and f_c is determined through dedicated trigger simulations, f_failf_other will be in turn constrained.

• All of the values for our hyper-parameters which are listed in Table 1 are backed up by literature. However, a recent study by Rouco Escorial et al. (2022) that came out during preparation of this manuscript suggests that slight alterations could be made. Specifically, the median for the jet opening angle of 6.1° [+3.2°, –9.3º] that the study presents differs from ours. This would yield an expected merger rate density of 360–1800 Gpc⁻³ yr⁻¹, that is a higher lower limit than our expected local rate of ~125 Gpc⁻³ yr⁻¹ with the top-hat jet approximation. These findings are especially relevant as the study is in essence a follow-up on Fong et al. (2015), on which we base our opening angle distribution. Where these do not necessarily denounce our results they do not vastly differ and, they indicate that a top-hat jet model with the current standard opening angle distribution may need to be altered in future. Any user interested in different opening angle distributions or other hyper-parameters can alter these relatively easily.

4.2 Applications on GRB physics and compact binaries merger history

As shown in Table 1, the underlying GRB population model of the simulator has hyper-parameters that describe the merger rate of the progenitor binaries and GRB physics. A comparison between the simulated catalogue with the real observation from the same detector can in turn give constraints on this parameters, and thus give information on the GRB physics and progenitor binaries merger history. More specifically, such inference can be done in a Bayesian fashion: $$p\left(ℬ|\left\{{\text{Θ}}_i\right\}\right)\propto ℒ\left(\left\{{\text{Θ}}_i\right\}|ℬ\right)p\left(ℬ\right),$$(22)

where $p\left(ℬ|\left\{{\text{Θ}}_i\right\}\right)$ is the posterior distribution of the list of the hyper-parameters of GRB population model, given a catalogue of detected GRB {Θ_i;}. Θ_i denotes a list of observables, for example T90, the fluence, redshift, E_peak et cetera, of the ith sGRB in the catalogues.

The key step is to find out the likelihood function on the right-hand side: ℒ({Θ_i}|ℬ), which has a known function form as in a inhomogeneous Poisson process: $$ℒ\left(\left\{{\text{Θ}}_i\right\}|ℬ\right)=N!{\displaystyle \prod _{i=1}^{N}p\left({\text{Θ}}_{i,\text{obs}}|ℬ\right)p\left(N|{\lambda }_{\text{GRB}}\left(ℬ\right)\right),}$$(23)

where N is the number of event in the catalogue, z_i is the red-shift of ith sGRB, p(z_i,_obs|ℬ) is probability of observing z_i given model, and λ_GRB is the expected total number in the catalogue, which is also function of model parameters ℬ. p(N|λ_GRB) is a Poisson probability with expectation λ_GRB. With simulations from GWToolbox, λ_GRB and p(Θ_i,obs|ℬ) can be obtained, and thus makes the study possible.

References

All Tables

All Figures

	Fig. 1 Redshift distribution of a simulated BNS merger catalogue from the GWToolbox for a LIGO observation time of 10 yr (blue) and one week of ET observations (orange). The underlying BNS population model is consistent with that in Table 1.
In the text

	Fig. 2 Inclination angle distribution of GW detection of different S/N threshold.
In the text

	Fig. 3 Built-in MRRs in the GWToolbox. We use the Brueckner-Hartree-Fock EoS (MPA1), two relativistic mean-field theory EOSs (MS0, MS2), and one nonrelativistic potential model (PAL1).
In the text

	Fig. 4 Schematic overview of the structure and geometry of our sGRB model, for a jet with opening angle Δθ, viewing angle θυ and distance to earth D. We depict an off-axis GRB in the top panel, while the bottom shows an on-axis one.
In the text

Fig. 5 Lightcurve of a GRB131004A-like on-axis GRB, according to the GWToolbox. We use the following hyper-parameters: γ = 120, log(r0/cm) = 12.95, log(EGRB/ergs) = 50.1, ${v^{\prime }}_0=0.7$ keV, Δθ = 15°, θυ = 10°, and z = 0.717. In blue, we show the simulated lightcurve of the on-axis GRB131004A, within an energy range of 10–1000 keV. The red data points correspond to the measured photon count rate by Fermi/GBM.

In the text

	Fig. 6 Lightcurve of a GRB 170817A-like off-axis GRB, according to the GWToolbox. We use γ = 120, log(r0/cm) = 13.03, log(EGRB/ergs) = 49.65, ${v^{\prime }}_0=8$ keV, Δθ = 15°, θυ = 18.60°, and z = 0.0099. In blue, we show the simulated lightcurve of the on-axis GRB170817A, within an energy range of 10–1000 keV. The red data points correspond to the measured photon count rate by Fermi/GBM.
In the text

	Fig. 7 Flowchart showing the sGRB detection test pipeline of BNSs.
In the text

	Fig. 8 Flowchart showing the sGRB detection test pipeline of BHNSs.
In the text

	Fig. 9 Fluence, 1-second peak photon flux, T90 and redshift histograms. We compare the GWToolbox simulation of Swift/BAT observing sGRBs for 1 yr (blue) with all sGRBs in the Swift/BAT GRB catalogue (red; Lien et al. 2016).
In the text

Fig. 10 Redshift distribution of BNSs through detection of GWs and sGRBs. In blue, we show the distribution for LIGO (design sensitivity), where the shaded region with the dotted line represents the distribution for those GW sources with S/N between 6 and 8. In orange, we plot the distribution of Fermi/GBM GRBs (the same distribution as in Figs. 11–13), where the red region corresponds to the off-axis GRBs, generated with the addition to our model as described in Sect. 3.4.

In the text

Fig. 11 Redshift distribution of BNSs through detection of GWs and sGRBs. In blue, we show the distribution for ET. Again, in orange, we plot the distribution of Fermi/GBM GRBs (same distribution as Figs. 10, 12 and 13), where the red region corresponds to the off-axis GRBs, generated with the addition to our model as described in Sect. 3.4. The orange region with dotted line represents those GRBs with a peak flux between 50 and 100% of the Fermi/GBM detection threshold.

In the text

	Fig. 12 Redshift distribution of LIGO-O3 BNSs (blue), and Fermi/GBM BNSs measured through sGRBs (orange). Once again, the Fermi/GBM distribution is the same as in Figs. 10, 11 and 13.
In the text

	Fig. 13 Redshift distribution of BNSs detected with a ‘half-ET’ configuration (blue) which has a sensitivity halfway between LIGO and ET, and Fermi/GBM BNSs measured through sGRBs (orange, same distribution as in Figs. 10–12). Also here, the shaded region with the dotted line corresponds to those sources with a S/N between 6 and 8. S/N8: 2038, S/N6: 6019 in 5 yr.
In the text