Linking the rates of neutron star binaries and short gamma-ray bursts.

Numerical simulations suggest only BNS mergers whose remnants promptly collapse into a black hole or produce a hypermassive neutron star (see Sarin and Lasky, 2021, for a review on BNS post-merger remnants), can launch jets that power short gamma-ray bursts (e.g., Margalit et al., 2015; Murguia-Berthier et al., 2017; Ciolfi, 2018). This constraint implies that any binary resulting in a remnant of mass $\lesssim 1.2M_{\rm{TOV}}{}$ does not produce a short gamma-ray burst. Here, $M_{\rm{TOV}}{}$ is the maximum, non-rotating neutron star mass known as the Tolman-Oppenheimer-Volkoff mass; a property of the unknown nuclear equation of state. Explicitly, this implies that $f_{\rm{s,BNS}}{}$ can be expressed as

Here $p_{\textrm{BNS}}$ is the probability of a successful jet launch, $\pi(m_{1},m_{2})$ is the astrophysical distribution of primary and secondary mass.

This relatively straightforward functional dependence has two caveats. First, X-ray plateaus of short gamma-ray bursts suggest that many gamma-ray bursts are powered by neutron star central engines (e.g., Rowlinson et al., 2013; Sarin et al., 2020), implying the $1.2\times M_{\rm{TOV}}{}$ constraint may be incorrect. Second, whether a jet is successful in producing prompt gamma-ray emission depends on the Lorentz factor of the jet and its threshold for producing prompt gamma-ray emission, both of which are unknown (Rhoads, 2003; Lamb and Kobayashi, 2016). Therefore, even when jets are successfully launched, they may not all produce prompt gamma-ray emission. We return to these two caveats in Sec. IV and V.

The question of whether an NSBH merger produces a short gamma-ray burst depends on whether enough matter is disrupted to create an accretion disk of sufficient mass around the remnant black hole. The mass in the disk must be large enough to feed a relativistic jet, and power it long enough so that it can break through the merger ejecta (e.g., Foucart, 2020). In general, there are two outcomes of an NSBH merger that depend primarily on the progenitor binary’s mass ratio and the black hole spin: 1) the neutron star plunges into the black hole without being tidally disrupted, or 2) the neutron star is disrupted. To good approximation, disruption occurs when the disruption radius $r_{\rm{disrupt}}{}$ is greater than the black hole’s innermost stable circular orbit (ISCO) $r_{\rm{ISCO}}{}$, i.e. $r_{\rm{disrupt}}{}\geq r_{\rm{ISCO}}{}$.

The disruption radius can be approximated by (e.g. Lee and Kluzniak (1999); Wiggins and Lai (2000))

We assume $k=1$, consistent with results found in the literature Lee and Kluzniak (1999); Wiggins and Lai (2000).

If the black hole spin is aligned with the orbital angular momentum, the ISCO radius can be approximated in terms of the black hole’s mass $M_{\rm BH}$ and dimensionless spin $\chi_{\rm BH}$ as

Here, $6\leq f(\chi_{\text{BH}})\leq 9$ for a co-rotating orbit, and $1\leq f(\chi_{\text{BH}})\leq 6$ for a counter-rotating orbit Bardeen et al. (1972). From Eqs. (3) and (4), we see that disruption is favoured for low mass ratios $q=M_{\text{NS}}/M_{\text{BH}}$ and high, co-rotating black hole spins. The disruption fraction $f_{\rm{s,NSBH}}{}$ can be calculated if we know the population distribution of NSBH progenitor masses and spins, as well as the nuclear equation of state that relates the neutron star progenitor mass to its radius. In other words, this disruption fraction can be written functionally as

where

is the probability of a successful jet launch in an NSBH merger (i.e., a Heaviside step function that evaluates to $0$ when $r_{\text{disrupt}}>r_{\text{ISCO}}$), and $\pi(m_{\text{NS}},m_{\text{BH}},\chi_{\text{BH}})$ is the astrophysical distribution of neutron star and black hole masses and black hole spins of the components that participate in NSBH mergers.

We calculate $f_{\rm{s,BNS}}{}$ and $f_{\rm{s,NSBH}}{}$ using the above prescription and a suite of results from Broekgaarden et al. (2021) that implements different population-synthesis prescriptions for the binary evolution physics—see that Paper, in particular their Tables 1 and 2, for a description of each of the models. From these various models, we take the mass distribution of neutron stars and black holes that participate in NSBH and BNS mergers and the orbital period before the second binary component explodes as a supernova, the latter informing the spin distribution of black holes that participate in NSBH mergers (Qin et al., 2018; Broekgaarden et al., 2021). We also marginalize over the unknown equation of state. In particular, we use piecewise polytropic equations of state, with polytropic indices uniformly distributed between 1.1 and 4.5, and neutron star central pressure between $10^{32.6}$ and $10^{33.6}\,\mathrm{dyne/cm^{2}}$ (see Ref. Hernandez Vivanco et al. (2020) for an explanation and implementation of these parameter choices), and the additional constraint that $M_{\rm TOV}$ is between 2.2 and $2.4\,\mathrm{M_{\odot}}$ consistent with combined constraints from multimessenger observations of GW170817, GW190425, and NICER observations (Raaijmakers et al., 2021). We also assume $0.05M_{\odot}$ ejecta in every BNS merger, consistent with observations of AT2017gfo (e.g., Smartt et al., 2017). For $f_{\rm{s,NSBH}}{}$ we also explore the affect of an additional constraint; that the disk mass from an NSBH merger must be above $0.075\,\mathrm{M_{\odot}}$ to ensure the jet can successfully break through the ejecta as hypothesised in e.g., Zappa et al. (2019). We calculate $f_{\rm{s,NSBH}}{}$ enforcing this constraint using relations for the ejecta and disk mass in Bhattacharya et al. (2019). We find that the discrepancy from the additional constraint of a minimum disk mass creates a change in $f_{\rm{s,NSBH}}{}$ smaller than due to the unknown equation of state and we therefore ignore this additional constraint in subsequent sections. The results where only disruption is necessary for $f_{\rm{s,NSBH}}{}$ are shown in Fig. 1.

Figure 1 shows the jet-launching fractions $f_{\rm{s,NSBH}}{}$ (orange) of NSBH mergers respectively for a diverse set of compas simulations (Broekgaarden et al., 2021) with various modifications. In particular, the rapid supernovae model (Rapid SNe), fixed supernovae kick (fixed SNe kick), no pulsational pair instability supernovae (No PISN), fixed common envelope efficiency (Fixed CE efficiency), unstable case BB mass transfer, and fixed mass transfer efficiency. For a full description of these simulations we refer the Reader to (Broekgaarden et al., 2021). In general, the different models predict that anywhere between $\approx 0.5-10\%$ of NSBH mergers disrupt sufficient matter to launch a jet accounting for the unknown equation of state. The smallest disruption values are found for the rapid supernova model (Fryer et al., 2012), which enforces a mass gap between neutron stars and black holes (cf. Özel et al., 2010; Farr et al., 2011). This mass gap combined with the relative lack of highly spinning black holes in our model (Qin et al., 2018; Chattopadhyay et al., 2021; Broekgaarden et al., 2021) makes it much more difficult to tidally disrupt the neutron star leading to lower values of $f_{\rm{s,NSBH}}$.

We also calculate the BNS jet-launching fraction $f_{\rm{s,BNS}}{}$ (blue) given the constraint that the remnant mass $M_{\rm{rem}}\gtrsim 1.2M_{\rm{TOV}}{}$ for the same set of compas simulations and two phenomenological distributions based on observations of neutron stars in our Galaxy using radio and gravitational waves (Alsing et al., 2018; Sarin et al., 2020; Galaudage et al., 2021). In particular, we use the mass distributions from Sarin et al. (2020) with a mixing fraction of GW170817-like and GW190425-like mergers of $50\%$, and the mass distribution inferred from the population of double neutron stars seen in radio and gravitational waves by Galaudage et al. (2021).

Most models predict $f_{\rm{s,BNS}}{}$ between $20-60\%$ apart from the compas models with fixed mass transfer efficiency and rapid supernovae (Fryer et al., 2012), which predict $f_{\rm{s,BNS}}{}={5.1}_{-0.2}^{+0.2}\%$ and $f_{\rm{s,BNS}}{}={4.9}_{-1.3}^{+2.1}\%$, respectively. That is, in these models only $5\%$ of BNS mergers launch a jet. The former model fixes the mass transfer efficiency to $0.5$, which leads to significantly lower quantities of heavy double neutron star systems, while the latter as mentioned above enforces a mass gap between neutron stars and black holes. In both models, mergers are more likely to form long-lived neutron star remnants which, within the jet-launching model we use (motivated by numerical results (e.g., Murguia-Berthier et al., 2017; Ciolfi, 2018)), are unable to produce short gamma-ray bursts. Such a low success rate is unlikely to be correct (e.g., Margalit and Metzger, 2019; Beniamini et al., 2019), which may hint at either of these models being incorrect or that the assumption of requiring a black hole to launch a jet capable of producing a short gamma-ray burst is flawed. We return to this point in Sec. IV.1 and Sec. V. Apart from these two models, all models predict a similar $f_{\rm{s,BNS}}{}$ between $20-60\%$. However, even a value of $60\%$ is potentially problematic in light of constraints from GW170817, which suggests that most BNS mergers should produce short a gamma-ray bursts (e.g., Beniamini et al., 2019).

Binary population-synthesis models typically predict a significant number of BNS with higher masses than those seen in the Galactic population (Vigna-Gómez et al., 2018; Chattopadhyay et al., 2020; Mandel and Smith, 2021). This discrepancy may be due to selection effects that prevent heavier neutron stars (as in GW190425 (Abbott et al., 2020)) from being observable in radio (e.g., Galaudage et al., 2021; Vigna-Gómez et al., 2021) or may hint at a systematic issue with population-synthesis results with regards to determining neutron star masses (Mandel and Smith, 2021). Our constraint that immediate formation of a black hole is required to launch a jet (i.e., $M_{\rm{rem}}\lesssim 1.2\times M_{\rm{TOV}}{}$), implies that population-synthesis results predict a higher $f_{\rm{s,BNS}}{}$ than one predicted by the distribution of BNS just in our Galaxy. For example, taking the distributions derived in Kiziltan et al. (2013) and Farrow et al. (2019), we find $f_{\rm{s,BNS}}{}={0.01}_{-0.01}^{+0.03}$ and $f_{\rm{s,BNS}}{}={0.02}_{-0.01}^{+0.03}$, respectively; such low values are unlikely to be correct. It is worth noting that these distributions were derived before the observation of GW190425, and are inconsistent with the masses inferred in that BNS.

The two terms $f_{\rm{s,BNS}}{}$ and $f_{\rm{s,NSBH}}{}$ dictate the fractions of BNS and NSBH mergers that can produce a short gamma-ray burst. A significant (but unknown) fraction of these jets are pointed away from Earth and are therefore not observable or are too dim to be detected by current gamma-ray detectors. We denote the observable fractions of short gamma-ray bursts from BNS and NSBH mergers as $\eta_{\rm{BNS}}$ and $\eta_{\rm{NSBH}}$, respectively.

Historically, the beaming fraction $f_{b}$ has been estimated through the relation $f_{b}=1-\cos(\theta_{\rm j}{})$ Frail et al. (2001), where $\theta_{\rm j}{}$ is the opening angle of the gamma-ray burst jet²²2This is sometimes defined as $f_{b}=(1-\cos{\theta_{\rm j}{}})^{-1}$, however we use the original definition as presented in Ref. (Frail et al., 2001). This implicitly assumes that the jet structure is top-hat with a large Lorentz factor such that all gamma-ray bursts observed within the opening angle $\theta_{\rm j}{}$ are observable and those outside are not. Constraints on the beaming fraction range from $\sim 10^{-3}-0.015$ (e.g., Berger, 2014), which implicitly assumes that short gamma-ray bursts produced by BNS and NSBH mergers have the same opening angles so that $\eta_{\rm{BNS}}=\eta_{\rm{NSBH}}$.

In general, we do not expect the beaming fractions of BNS and NSBH jets to be the same given the different jet-launching environment. The beaming fraction is a function of the initial jet structure, jet-launching mechanism, and the propagation of the jet through the surrounding ejecta. Typically, BNS mergers have more ejecta than NSBH mergers along the poles (e.g., Foucart, 2020; Kyutoku et al., 2021), which is the region where jets are launched. Numerical simulations already show that propagation through ejecta helps collimate gamma-ray burst jets Barbieri et al. (2019); Urrutia et al. (2021). We therefore believe that gamma-ray burst jets launched in NSBH mergers are likely to be less collimated than ones launched in BNS mergers. Similarly, depending on binary properties, some NSBH mergers are unable to form a disk of sufficient mass to power a jet long enough to make its way out of the merger ejecta (e.g., Zappa et al., 2019), which will likely lead to a cocoon-like outflow. It is worth noting that the lower amount of overall mass in the polar region may alleviate this concern (Just et al., 2016; Kyutoku et al., 2021).

The above arguments implicitly assume gamma-ray emission is not produced outside the ultra-relativistic core, something we now know to be incorrect following observations of GRB170817A Matsumoto et al. (2019). These observations also suggest that other gamma-ray bursts may be viewed from outside the opening angle of the jet, a suggestion that has led to reclassification of some faint short gamma-ray bursts such as GRB150101B as ones viewed off axis Troja et al. (2018). The presence of off-axis gamma-ray bursts makes it difficult to parameterise the beaming fraction as simply being related to the opening angle of the jet. To further complicate matters, there is no robust model for calculating the flux of the prompt emission that can robustly predict the vast range of gamma-ray burst energies. In light of GRB170817A, one recent approach to estimate the flux of the prompt emission is to compute the peak flux at a given observer viewing angle using the same structured jet profile used to describe the afterglow for a given gamma-ray detector band-pass and threshold (e.g., Howell et al., 2019; Salafia et al., 2019). For our purposes, we use this aforementioned prescription, accompanied by a choice of jet structure, to compute the flux for an isotropic distribution of observers and calculate the fraction of short gamma-ray bursts that are observable with the Fermi telescope.

To understand the beaming fractions of BNS $\eta_{\rm{BNS}}$ and NSBH $\eta_{\rm{NSBH}}$, we use two different viewing-angle dependent energy distributions for the jet structures: a Gaussian and a power-law structured jet respectively given by

Here $E_{0}$ is the on-axis isotropic equivalent energy, $\theta$ is the angle from the jet axis, $\theta_{\rm core}{}$ is the half opening angle of the ultra-relativistic core of the jet, and $b$ is the power-law slope. Often these models are parameterised with an additional angle $\theta_{\rm{wing}}$, beyond which the energy is zero, implying no gamma rays can be produced beyond this angle. This angle is typically taken to be some integer multiple of $\theta_{\rm core}{}$ (e.g., Cunningham et al., 2020). The parameters of the Gaussian and power-law jet structured models are different for BNS and NSBH as we expect BNS jets to be more collimated, as described in Sec. II.

In Fig. 2 we plot the beaming fractions $\eta_{\rm{BNS}}$ and $\eta_{\rm{NSBH}}$ for a realistic population of BNS and NSBH out to 330 Mpc and 590 Mpc (corresponding to the horizon distances for aLIGO operating at design sensitivity (Abbott et al., 2018), where the NSBH horizon distance assumes $1.4\,\mathrm{M_{\odot}}$ + $10\,\mathrm{M_{\odot}}$ systems), and viewing angles that are uniform in $\cos(\iota)$. The red violins correspond to the Gaussian structured jet while the blue violins correspond to the power-law structured jet. The priors for the BNS jets are motivated by fits to the afterglow of GRB170817A (Ghirlanda et al., 2019) and expectations of structured jets produced in NSBH mergers (Barbieri et al., 2019).

Figure 2 shows two main trends: 1) power-law structured jets tend to have higher beaming fractions than gaussian jets and 2) NSBH jets have higher beaming fractions than BNS jets. The former is due to the energy distribution of power-law jets as opposed to Gaussian jets, with power-law jets having more energy at higher angles away from the jet axis than Gaussian. The latter is a realisation of our prior assumption that NSBH jets are less collimated. It is worth noting that the range of beaming fractions predicted by either a Gaussian or power-law structured jet is broad for both BNS and NSBH jets and is an encapsulation of both our uncertainty in the jet-structure parameters and their priors and the unknown prompt emission generation mechanism. This uncertainty will improve with future multi-messenger events like GW170817 (Beniamini et al., 2019; Biscoveanu et al., 2020; Farah et al., 2020). In general, we find that $\eta_{\rm{BNS}}\lesssim 0.1$ and $\eta_{\rm{NSBH}}\lesssim 0.3$ for the less collimated power-law jet model at $90\%$.

In the following sections, we use these models to estimate constraints on each of these terms given current short gamma-ray burst, BNS and NSBH merger-rate measurements. We also forecast each of these rates assuming future gravitational-wave and electromagnetic observations.

As described in Sec. II, each term in Eq. 1 is a function of several physical parameters. These parameters can be informed through gamma-ray burst modelling and population synthesis studies, as well as inferences on the observed population. As a first attempt, we consider constraints imposed by a model with minimal assumptions and directly comparable to historical measurements. Our minimal working model assumes the beaming fraction of both BNS and NSBH is the same ($\eta_{\rm{BNS}}=\eta_{\rm{NSBH}}$), and is parameterised only by the jet opening angle i.e., that all gamma-ray bursts are observable as long as they are observed within some opening angle $\theta_{\rm j}$. To wit

We place priors directly on $f_{\rm{s,BNS}}{}$ and $f_{\rm{s,NSBH}}{}$, which uses two pieces of information. First, that we expect a large majority of BNS mergers to produce short gamma-ray bursts (e.g., Margalit and Metzger, 2019; Beniamini et al., 2019) and second, that the two bona fide NSBH mergers detected with gravitational waves should not have produced short gamma-ray bursts (e.g., Abbott et al., 2021a). These two statements can be translated into priors as

Equation 10 is a uniform distribution motivated by observations of GW170817 (Beniamini et al., 2019). Equation 11 is a binomial probability distribution for zero successes from two independent trials, where $p$ is the success probability for each trial. We also enforce that $f_{\rm{s,BNS}}{}\gtrsim f_{\rm{s,NSBH}}{}$ due to theoretical motivations outlined in Sec. II.

We evaluate Eq. 9 using the above priors and rate measurements $\mathcal{R}_{\rm{BNS}}{}=320^{+490}_{-240}$ Gpc^-3yr^-1 (Abbott et al., 2021b), $\mathcal{R}_{\rm{NSBH}}{}=130^{+112}_{-69}$ Gpc^-3yr^-1 (Abbott et al., 2021a), and $\mathcal{R}_{\rm{SGRB}}{}=8^{+5}_{-3}$ Gpc^-3yr^-1 (e.g., Coward et al., 2012; Wanderman and Piran, 2015). With the assumption that both NSBH and BNS merger jets have the same opening angle, measurements of the merger and short gamma-ray burst rates place strong constraints on the opening angle of the jet and the BNS merger rate. We show this two-dimensional marginalised posterior distribution in Fig. 3. We find that the opening angle of gamma-ray burst jets is constrained to $\theta_{\rm j}{}={15^{\circ}}^{+4^{\circ}}_{-5^{\circ}}{}$ ($90\%$ credible interval), consistent with other values in the literature (e.g., Berger, 2014; Williams et al., 2018). Note that we place a hard cut off at $\theta_{\rm j}{}=20^{\circ}$ as larger opening angles would be inconsistent with constraints on short gamma-ray burst energetics (Frail et al., 2001; Fong et al., 2012). Similarly, the BNS merger rate is constrained to $\mathcal{R}_{\rm{BNS}}{}=375^{+436}_{-211}\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}}$, which is approximately $15\%$ more informative than the second gravitational-wave catalog merger rate estimate (Abbott et al., 2021b). We measure $f_{\rm{s,BNS}}{}=0.69^{+0.27}_{-0.37}$ ($90\%$ credible interval) and that $f_{\rm{s,BNS}}\gtrsim 0.3$ with $95\%$ confidence ruling out population synthesis models with rapid supernovae and fixed mass transfer efficiency (see Fig. 1). Within these assumptions, the NSBH merger rate and jet launching fraction is not constrained from the prior. We also find an anti-correlation between $\theta_{\rm j}{}$ and $\mathcal{R}_{\rm{BNS}}{}$ which agrees with physical intuition: small jet opening angles must be compensated by a larger BNS merger rate to explain the measured short gamma-ray burst rate.

We now relax the assumption that jets launched in NSBH mergers have the same opening angle as those in BNS mergers; as described in Sec. II, we expect the former to be less collimated than the latter. We again evaluate Eq. 9 using current rate measurements, but now assume $\theta_{\rm j}{}$ is different for NSBH and BNS, parameterising both as uniform distributions between $1-20\,\mathrm{{}^{\circ}}$ motivated by estimates of the beaming fraction of short gamma-ray bursts and their overall energetics (e.g., Fong et al., 2012; Berger, 2014). The other priors are the same as used above. We also enforce $\theta_{\rm j,NSBH}{}\gtrsim\theta_{\rm j,BNS}{}$ due to theoretical motivations outlined in Sec. II. The posterior distribution on a subset of model parameters are shown in Fig. 4.

Under this model with fewer assumptions, we measure BNS and NSBH opening angles to $\theta_{\rm j,BNS}{}={14.8^{\circ}}^{+4.0^{\circ}}_{-5.1^{\circ}}{}$ and $\theta_{\rm j,NSBH}{}={15.3^{\circ}}^{+4.0^{\circ}}_{-5.1^{\circ}}{}$ respectively. These opening angles are effectively the same, and may suggest that within the model where the beaming fractions are given by $1-\cos{\theta_{\rm j}{}}$, there is no significant difference between BNS and NSBH jets. This may imply that there is a maximum angle away from the jet axis at which prompt gamma-ray radiation cannot be produced, and that this is the same for both BNS and NSBH jets or that the difference is smaller than we can currently probe.

The analysis presented hitherto also provides updated rate measurements that are slightly improved over current estimates in the literature. We measure the BNS merger rate as $\mathcal{R}_{\rm{BNS}}{}=384^{+431}_{-213}\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}}{}$ and the NSBH merger rate as $\mathcal{R}_{\rm{NSBH}}{}=132^{+109}_{-70}\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}}{}$ ($90\%$ credible interval), which are approximately $16\%$ and $3\%$ more informative than the BNS and NSBH merger rates from the second gravitational-wave catalog (Abbott et al., 2021b), respectively.

The previous Section used the $f_{b}=1-\cos{\theta_{\rm j}{}}$ relation as an estimate of the “beaming fraction” to allow direct comparison with historical measurements. However, in reality, we do not expect this relation to be correct for two reasons. 1) It is unlikely all gamma-ray bursts have the same jet opening angle and 2) prompt gamma-ray radiation can be produced outside the ultra-relativistic core as has been argued for GRB170817A (Matsumoto et al., 2019). Here we address the impact of ignoring these two issues on our analysis.

As described in Sec. II, interaction with the ejecta around the polar region helps collimate jets (Barbieri et al., 2019; Urrutia et al., 2021). We do not expect every BNS merger or NSBH merger to eject the same amount of matter so one cannot reasonably expect there to be a universal opening angle $\theta_{\rm j}{}$ across multiple BNS mergers, let alone between BNS and NSBH mergers. Moreover, both the propagation of a jet and the overall jet structure are strongly connected to the initial jet structure and energetics (Nativi et al., 2021). Both of these are sensitive to binary parameters such as the masses, spins, nuclear equation of state, etc, although it is worth noting that no robust physical relationship exists in the literature linking these variables.

To investigate the effect of assuming a universal jet opening angle on our results above, we perform a more realistic simulation. In particular, we rewrite Eq. 9 as

where $\pi(\theta_{\rm j}{})$ is the astrophysical distribution of jet opening angles relaxing the assumption implicit in Eq. 9 of a universal opening angle. We simulate the short gamma-ray burst rate using Eq. 12 with $\mathcal{R}_{\rm{BNS}}{}=320\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}{}}$, $\mathcal{R}_{\rm{NSBH}}{}=130\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}{}}$, $f_{\rm{s,BNS}}{}=0.3$, $f_{\rm{s,NSBH}}{}=0.02$ and $\theta_{\rm j}{}$ drawn from a normal distribution with mean $\mu=15^{\circ}$, and $\sigma=4^{\circ}$. We then perform our analysis by estimating the average opening angle $\theta_{j,avg}$ with uniform priors on $f_{\rm{s,BNS}}{}$, $f_{\rm{s,NSBH}}{}$, and $\theta_{j,avg}$ between $0.01-0.2$, $0.1-1$, and $1-30^{\circ}$, respectively. We use the measured merger rates from the second gravitational-wave transient catalog as our priors for $\mathcal{R}_{\rm{BNS}}{}$ and $\mathcal{R}_{\rm{NSBH}}{}$.

In Figure 5, we show the one and two-dimensional posterior distributions of all parameters obtained by fitting the simulated data. The orange lines refer to the injected values while the green curves represent the prior. All parameters are recovered correctly, with the injected parameter within the $1\sigma$ credible interval, suggesting that given our current range of uncertainty on the rate of NSBH and BNS mergers and short gamma-ray bursts ignoring the distribution of $\theta_{\rm j}{}$ does not bias our results. However, we find that the uncertainty on $\theta_{j,avg}$ is larger by $\approx 40\%$ compared to analysis where a single universal jet opening angle is assumed implying that the uncertainties on $\theta_{\rm j}{}$ derived in Sec. III.1 are underestimated by at least $40\%$.

In the previous Section, we implicitly assumed a “top-hat” jet structure i.e., that there is a maximum opening angle $\theta_{\rm j}{}$ outside of which we observe no gamma-ray emission. As discussed in Sec. II, we know from observations of GRB170817A, that this is incorrect and that the jet is likely structured.

In this Section, we link the rates of short gamma-ray bursts and neutron star binaries together using Eq. 1 with priors on $\eta_{\rm BNS}$ and $\eta_{\rm NSBH}$ informed by the Gaussian structured jet model i.e., the beaming fractions shown in Fig. 2 in red. We use the same priors for all other parameters as in Sec. III.1 and enforce $\eta_{\rm NSBH}\gtrsim\eta_{\rm NSBH}$ due to the theoretical reasons outlined in Sec. II.

With current measurements, our constraints with a structured-jet model are similar to the ones obtained with a top-hat model in Sec. III.1. In particular, we measure $f_{\rm{s,BNS}}{}=0.69^{+0.28}_{-0.37}$ and $\mathcal{R}_{\rm{BNS}}{}=360^{+407}_{-238}\,\mathrm{\rm{Gpc}^{-3}\rm{yr}^{-1}}$ with $\mathcal{R}_{\rm{NSBH}}{}$ and $f_{\rm{s,NSBH}}{}$ practically indistinguishable from the prior. This suggests that given current measurements, and the priors used on $f_{\rm{s,BNS}}{}$ and $f_{\rm{s,NSBH}}{}$, this framework can not probe the difference between a top-hat jet and a Gaussian jet structure. The analysis does provide an improved constraint on the beaming fraction of $\eta_{\rm{BNS}}=0.03^{+0.05}_{-0.02}$ and $\eta_{\rm{NSBH}}=0.08^{+0.07}_{-0.05}$ for BNS and NSBH jets, approximately $60\%$ more informative than the prior. However, with structured-jet models the “beaming fractions” themselves are not particularly enlightening, and one would instead want constraints on the parameters such as the distribution of on-axis isotropic equivalent energy, $E_{0}$, or the half-opening angle of the ultra-relativistic core, $\theta_{\rm core}{}$ that determine the beaming fraction. These constraints require Bayesian hierarchical inference, a development we leave to future work. We note that given current measurements, the constraints on these parameters will not be particularly informative but could become a powerful complementary way to probe the jet structure in the future.

Constraining the fraction of BNS that successfully launch a jet ($f_{\rm{s,BNS}}{}$) has several important implications. As mentioned in Sec. II, numerical simulations suggest a gamma-ray burst jet is only produced when the remnant mass of a BNS merger is $\gtrsim 1.2\times M_{\rm{TOV}}{}$. Depending on the BNS mass distribution, this could imply that a significant fraction of BNS mergers do not produce short gamma-ray bursts, which must be represented in the respective rates. However, X-ray afterglow observations of several short gamma-ray bursts have features that are difficult to explain without requiring a neutron star engine (e.g., Rowlinson et al., 2013; Sarin et al., 2020). These X-ray observations indicate that $\gtrsim 30\%$ of short gamma-ray bursts produce long-lived neutron stars. Placing the constraint that a jet is only launched if $M_{\rm{rem}}\gtrsim 1.2\times M_{\rm{TOV}}{}$, our different suite of models combined predict a BNS jet-launching fraction $f_{\rm{s,BNS}}{}={0.30}_{-0.24}^{+0.31}$ ($90\%$ credible interval). This fraction is lower than expectations in light of GW170817 (Margalit and Metzger, 2019; Beniamini et al., 2019). For example, Beniamini et al. (2019) predicted $f_{\rm{s,BNS}}{}\approx 0.6-1$ using the short gamma-ray burst luminosity function (Wanderman and Piran, 2015) and various structured jet models, hinting at a potential inconsistency.

Here, we explore how many observations of BNS and NSBH mergers are necessary to determine whether a black hole central engine is necessary to launch a jet i.e., how long to confidently determine $f_{\rm{s,BNS}}{}\neq 1$? We create simulated data following Eq. 1 with $\eta_{\rm{BNS}}=0.02$ and $\eta_{\rm{NSBH}}=0.04$ for BNS and NSBH beaming fractions, respectively, and $f_{\rm{s,BNS}}{}=0.3$ and $f_{\rm{s,NSBH}}{}=0.02$ for the jet-launching fractions of BNS and NSBH mergers. We estimate the constraints for varying number of observations, simulating new constraints on the merger rates using the current distribution and decreasing the uncertainty with $1/\sqrt{N}$ where $N$ is the number of observations. We use broad uninformative priors on the jet-launching fractions and priors informed by the Gaussian jet structure for the beaming fraction. We also enforce the constraint that $f_{\rm{s,BNS}}{}\gtrsim f_{\rm{s,NSBH}}{}$ and $\eta_{\rm{BNS}}\gtrsim\eta_{\rm{NSBH}}$.

Within a single year of observations at design sensitivity, we expect $48$ and $111$ BNS and NSBH mergers allowing us to confirm if most BNS mergers can successfully launch a jet i.e., rule out $f_{\rm{s,BNS}}{}\gtrsim 0.8$ with $90\%$ confidence. We note that this may not necessarily solve the problem in our understanding of whether short gamma-ray bursts can be produced by mergers producing long-lived neutron stars, as the missing fraction could be gamma-ray burst jets that have low Lorentz factors such that they can not produce prompt gamma-ray emission. However, by this time we may likely have a good understanding of the BNS mass distribution (Chatziioannou and Farr, 2020) such that we could isolate this fraction of ‘failed’ gamma-ray bursts (e.g., Lamb and Kobayashi, 2016) from the fraction of remnants with mass $M_{\rm rem}\lesssim 1.2\times M_{\rm{TOV}}{}$, providing an answer to one of the most fundamental questions in gamma-ray burst physics.

The afterglow observations of GRB170817A confirmed that gamma-ray burst jets are structured, but the exact jet structure is unclear. However, more pertinent to our discussion is whether the jet structures of gamma-ray bursts launched in NSBH and BNS mergers are different and ultimately whether they have the same beaming fraction i.e., whether $\eta_{\rm{NSBH}}=\eta_{\rm{BNS}}$. As described in Sec. II, we expect the structures and therefore beaming fractions to be different. This is borne out from numerical simulations which suggest that gamma-ray burst jet structure is a shaped largely by the jet’s interaction with the ejecta (Urrutia et al., 2021; Nativi et al., 2021) and the expectation that the amount and location of ejecta in BNS and NSBH mergers are different (Kyutoku et al., 2021).

We use the same procedure and parameters as outlined above to estimate when we will be able to rule out both beaming fractions being equal. However, unlike previously, we use priors motivated by the fixed supernovae kick population synthesis model on the jet-launching fraction and broad uninformative priors on $\eta_{\rm{BNS}}$ and $\eta_{\rm{NSBH}}$. We again enforce the constraint that $f_{\rm{s,BNS}}{}\gtrsim f_{\rm{s,NSBH}}{}$ and $\eta_{\rm{BNS}}\gtrsim\eta_{\rm{NSBH}}$.

Within $6$ months of observation at design sensitivity, we can confirm if the beaming fraction of short gamma-ray bursts launched in NSBH and BNS are different and in what way. In particular, with $6$ months of observation for this simulation, we can constrain $\eta_{\rm{BNS}}={0.02}_{-0.02}^{+0.02}$ and $\eta_{\rm{NSBH}}={0.16}_{-0.12}^{+0.11}$ (90% credible interval).

These constraints will help shape our understanding of how gamma-ray burst jets get their structure and could be further decomposed into constraints on the population of parameters corresponding to the Gaussian and power-law structured jets. This increase in understanding will be critical for estimating the Hubble constant with coincident gravitational-wave and electromagnetic observations in the absence of very long baseline interferometry (Schutz, 1986; Abbott et al., 2017c; Hotokezaka et al., 2019).

All the models investigated in this work predict jet-launching fractions $f_{\rm{s,BNS}}{}\lesssim 0.6$ and $f_{\rm{s,NSBH}}{}\lesssim 0.1$. The former is enforced due to the constraint that a gamma-ray burst jet is not produced unless the remnant mass is greater than $1.2\times M_{\rm{TOV}}{}$. The latter is dependent on two critical factors. 1) Whether there is a mass gap between neutron stars and black holes and 2) the spins of black holes that participate in NSBH mergers.

Recently, Fragione (2021) estimated the fraction of NSBH mergers that tidally disrupt the neutron star for several different equations of states and for different black-hole spin models. One model in their work, which estimates black hole spins using single star models from the Geneva stellar evolution code (Eggenberger et al., 2008; Ekström et al., 2012), predicts the disruption fraction $f_{\rm{s,NSBH}}{}\gtrsim 0.5$ for all different configurations. The large $f_{\rm{s,NSBH}}{}$ is attributed to the rapidly rotating black holes produced by the Geneva code, a byproduct of inefficient angular momentum transfer in the massive stars. Such a high disruption fraction, and large black hole spins, appear to be inconsistent with the gravitational-wave data (Abbott et al., 2021a) and the lack of observed electromagnetic radiation from GW200105 and GW200115 (e.g., Zhu et al., 2021; Dichiara et al., 2021; Anand et al., 2021). We note that this may just be a selection effect, as LIGO preferentially sees more massive systems that are less likely to produce electromagnetic radiation.

To investigate whether such high disruption fractions can be ruled out with the data, we simulate data with $f_{\rm{s,BNS}}{}=0.1$, and $f_{\rm{s,NSBH}}{}=0.8$ with all other parameters kept the same as above. We again use uninformative priors on the jet launching fractions while the beaming fractions are informed by the Gaussian structured jet model.

Within $2$ years of observation at design sensitivity, we can rule out $f_{\rm{s,NSBH}}{}\lesssim 0.5$ with $95\%$ confidence. If $f_{\rm{s,NSBH}}{}\gtrsim 0.5$ then this has important implications for binary evolution, as it could either rule out a mass gap between neutron stars and black holes, or imply that the spins of black holes in NSBH are typically much larger than predicted by the tidal spin up model (Qin et al., 2018). The former would provide insights on the formation of neutron stars and black holes in supernovae, whilst the latter probes the efficiency of angular momentum transport in massive stars and the role of tides in binary evolution.

Neutron star-black hole mergers are often invoked to explain peculiar observations of some short gamma-ray bursts such as the extended prompt emission or internal plateaus (e.g., Desai et al., 2019; Gompertz et al., 2020). However, as noted above, both these features are more typically attributed to the spin-down of a nascent neutron star (e.g., Rowlinson et al., 2013), an engine that can not be born in a NSBH merger. The fraction of short gamma-ray bursts that exhibit such features is $\approx 30\%$ higher than any of our predictions.

With parameters informed by the suite of population-synthesis results described earlier and the Gaussian structured jet, we will be able to measure $f_{\rm{s,BNS}}{}$ and $f_{\rm{s,NSBH}}{}$ with $10\%$ precision within $6$ months of observation at design sensitivity. This precision may allow us to determine if the contaminated sample of short gamma-ray bursts with internal plateaus and extended emission are consistent with the disruption fraction of NSBH mergers. This would provide evidence for mechanisms such as fall-back accretion onto the black hole for explaining these peculiar observations (e.g., Desai et al., 2019).

Under the constraint that $M_{\rm{rem}}\lesssim 1.2M_{\rm{TOV}}{}$, all population synthesis models analysed in this work predict $f_{\rm{s,BNS}}{}\lesssim 0.5$. Models motivated by Galactic and extragalactic observations of double neutron stars (e.g., Sarin et al., 2020; Galaudage et al., 2021) predict slightly higher values of $f_{\rm{s,BNS}}{}\lesssim 0.6$.

As described in Sec. IV.1, we can rule out if $f_{\rm{s,BNS}}{}\gtrsim 0.8$ with $90\%$ confidence within a year of observations. Such a high value of $f_{\rm{s,BNS}}{}$ is inconsistent with any of our model predictions and may suggest that either there are significantly more higher mass mergers than those predicted by population synthesis and current observations of gravitational-wave mergers. Or, that the assumption that $M_{\rm{rem}}\gtrsim 1.2M_{\rm{TOV}}{}$ is required to produce a gamma-ray burst jet is flawed. Both hypotheses have profound implications on our understanding of binary evolution and gamma-ray burst jet launching.