Fermi GBM follow-up of LIGO-Virgo binary black hole mergers -- detection prospects

TL;DR

This paper assesses the likelihood of detecting gamma-ray counterparts to LIGO-Virgo binary black hole mergers with Fermi-GBM, using simulations to estimate detection ratios and constrain emission models.

Contribution

It introduces a simulation-based framework to estimate the gamma-ray detection prospects for binary black hole mergers and constrains emission models based on expected detection ratios.

Findings

Most models predict a BBH-to-GRB ratio between 5 and 20.

Optimistic models suggest ratios as low as 2.

Pessimistic models could have ratios up to 700.

Abstract

Fermi-Gamma-ray Burst Monitor observed a 1 s long gamma-ray signal (GW150914-GBM) starting 0.4 s after the first gravitational wave detection from the binary black hole merger GW150914. GW150914-GBM is consistent with a short gamma-ray burst origin; however, no unambiguous claims can be made as to the physical association of the two signals due to a combination of low gamma-ray flux and unfavorable location for Fermi-GBM. Here we answer the following question: if GW150914 and GW150914-GBM were associated, how many LIGO-Virgo binary black hole mergers would Fermi-GBM have to follow up to detect a second source? To answer this question, we perform simulated observations of binary black hole mergers with LIGO-Virgo and adopt different scenarios for gamma-ray emission from the literature. We calculate the ratio of simulated binary black hole mergers detected by LIGO-Virgo to the number of gamma-ray counterpart detections by Fermi-GBM, BBH-to-GRB ratio. A large majority of the models considered here predict a BBH-to-GRB ratio in the range of 5 to 20, but for optimistic cases can be as low as 2 or for pessimistic assumptions as high as 700. Hence we expect that the third observing run, with its high rate of binary black hole detections and assuming the absence of a joint detection, will provide strong constraints on the presented models.

Figures

Tables

Table 1. Table 1: Central values for the observed BBH merger parameters, taken from The LIGO Scientific Collaboration et al. ( 2018 ) . Masses are in units of Solar masses. Numbers next to GW events are for referencing them in figures.

BBH name	M1	M2	Mf	$a_f$	DL/Mpc	$z$	E${}_{\mathrm{GW}}{}^{}/M_{\odot }{c}^2$	SNR
GW150914	35.6	30.6	63.1	0.69	430	0.09	3.1	23.7
GW151012 (1)	23.3	13.6	35.7	0.67	1060	0.21	1.5	9.7
GW151226 (2)	13.7	7.7	20.5	0.74	440	0.09	1.0	13.0
GW170104 (3)	31.0	20.1	49.1	0.66	960	0.19	2.2	13
GW170608 (4)	10.9	7.6	17.8	0.69	320	0.07	0.9	13
GW170729 (5)	50.6	34.3	80.3	0.81	2750	0.48	4.8	13
GW170809 (6)	35.2	23.8	56.4	0.70	990	0.20	2.7	13
GW170814 (7)	30.7	25.3	53.4	0.72	580	0.12	2.7	15
GW170818 (8)	35.5	26.8	59.8	0.67	1020	0.20	2.7	15
GW170823 (9)	39.6	29.4	65.6	0.71	1850	0.34	3.3	15

Table 2. Table 2: The BBH-to-GRB ratio in different scenarios. We consider three types of opening angles: fixed narrow ( θ j = 20 ∘ subscript 𝜃 𝑗 superscript 20 \theta_{j}=20^{\circ} ), distributed randomly between 10 and 40 degrees and isotropic ( θ j = 90 ∘ subscript 𝜃 𝑗 superscript 90 \theta_{j}=90^{\circ} ). We show two limits in GW SNR (8 and 12). The energy extraction scenarios are described in Section 3 , while the 3 cases are presented in Section 3.6 . Column 5 is highlights the benchmark case, indicating the most similarity to short GRB opening angles and GW observations.

Scenario	Section	${\theta }_j={20}^{\circ }$		${\theta }_j=$Unif.${[10,40]}^{\circ }$		${\theta }_j={90}^{\circ }$
		SNR$>$8	SNR$>$12	SNR$>$8	SNR$>$12	SNR$>$8	SNR$>$12
neutrino wind, pessimistic		31.1	29.1	20.7	15.3	5.6	3.9
neutrino wind, normalized	3.1	15.7	14.6	10.4	8.6	2.8	2.5
neutrino wind, optimistic		10.7	10.5	7.1	6.5	2.0	1.9
Blandford Znajek, pessimistic		304.7	181.9	146.7	94.9	33.4	24.0
Blandford Znajek, normalized	3.2	44.5	32.6	34.9	21.8	7.1	5.0
Blandford Znajek, optimistic		16.4	14.1	10.8	9.0	2.9	2.4
Charged BH, pessimistic		158.4	121.3	101.6	66.2	22.3	13.9
Charged BH, normalized	3.3	22.6	21.0	14.6	11.7	3.9	2.9
Charged BH, optimistic		10.0	9.7	6.5	5.9	1.9	1.8
$L_{\gamma }$=E${}_{\mathrm{GW}}{}^{}\times {ϵ}_{\mathrm{GW}}$, pessimistic		720.2	218.3	565.9	121.3	55.4	21.6
$L_{\gamma }$=E${}_{\mathrm{GW}}{}^{}\times {ϵ}_{\mathrm{GW}}$, normalized	3.4	176.0	62.4	99.0	29.1	11.7	5.1
$L_{\gamma }$=E${}_{\mathrm{GW}}{}^{}\times {ϵ}_{\mathrm{GW}}$, optimistic		31.7	17.1	18.3	9.2	3.7	2.2
$L_{\gamma }\propto M_f^{-2.0}$, pessimistic		16.9	15.0	11.0	8.7	3.0	2.5
$L_{\gamma }\propto M_f^{-2.0}$, normalized	3.5	13.4	12.3	9.1	7.4	2.5	2.2
$L_{\gamma }\propto M_f^{-2.0}$, optimistic		11.9	10.9	7.8	6.9	2.2	2.0
$L_{\gamma }$=$5.7\times {10}^{48}$ erg s-1(const.)		114.8	54.6	64.9	28.7	10.7	5.0
$L_{\gamma }$=$1.8\times {10}^{49}$ erg s-1(const.)	3.5	31.1	18.3	19.6	10.2	4.2	2.7
$L_{\gamma }$=$5.7\times {10}^{49}$ erg s-1(const.)		14.6	11.8	9.6	7.2	2.5	2.0
$L_{\gamma }\propto$ random, pessimistic		79.2	42.8	49.2	22.1	8.8	4.8
$L_{\gamma }\propto$ random, normalized	3.7	29.8	20.6	18.0	11.2	4.1	2.8
$L_{\gamma }\propto$ random, optimistic		15.4	12.5	9.6	7.6	2.6	2.1

Equations

$${\mathrm{SNR}}^2=4\int \frac{|h_{+}{F}_{+}+h_{\times }{F}_{\times }{|}^2}{S_h(f)}df.$$

$$\mathrm{SN}{\mathrm{R}}_{\mathrm{S}}^2=4\frac{{\mathcal{A}}^2}{D_L^2}\left(⟨F_{+}^2⟩(1+{\cos }^2\iota {)}^2+4⟨F_{\times }^2⟩{\cos }^2\iota \right)\underset{f_{\mathrm{low}}}{\overset{f_{\mathrm{ISCO}}}{\int }}\frac{f^{-7/3}}{S_h(f)}df,$$

$$P(\iota )\propto (1+6{\cos }^2\iota +{\cos }^4\iota {)}^{3/2}\sin \iota .$$

$$L_{\nu \bar{\nu }}\approx 2.9\times 10^{47}{\left(\frac{f_3(a_f)}{f_3(0.68)}\right)}^{-4.8}{\left(\frac{M_f}{62M_{\odot }}\right)}^{-3/2}{\left(\frac{\dot{M}}{M_{\odot }{s}^{-1}}\right)}^{9/4}erg{s}^{-1}$$

$$L_{\mathrm{BZ}}=\frac{1}{3c}{\left(\frac{{\Psi }_h{\Omega }_h}{4\pi }\right)}^2=\frac{1}{12c^3}{G}^2{M}_f^2{B}^{2}f(a_f),$$

$$L_Q=C\frac{c^5}{G}{\hat{q}}^2{\hat{a}}^{-15},$$

$$P\left(x=0|N_{\mathrm{BBH}},p=\frac{1}{N_{\mathrm{BGR}}}\right)={\left(1-\frac{1}{N_{\mathrm{BGR}}}\right)}^{N_{\mathrm{BBH}}}.$$

$$N_{\mathrm{BBH}}=\frac{\log P_{\mathrm{conf}}}{\log \left(1-\frac{1}{N_{\mathrm{BGR}}}\right)}\approx -N_{\mathrm{BGR}}\log P_{\mathrm{conf}}.$$

Full text

gravitational waves, gamma rays: general, gamma-rays: GW150914-GBM

1 Introduction

Gravitational waves (GWs) were directly observed for the first time (GW150914) from the coalescence of two black holes (BH; Abbott et al., 2016b), opening a new window on the universe. There is an intense, ongoing effort to detect electromagnetic (EM) counterparts for GW observations (Abbott et al., 2016a; Racusin et al., 2017; Abbott et al., 2017a; Burns et al., 2019; Abbott et al., 2019). This effort resulted in the detection of the first unambiguous gamma-ray counterpart to a GW signal (Goldstein et al., 2017; GBM-LVC, 2017; Abbott et al., 2017a), from a binary neutron star merger.

Fermi Gamma-ray Burst Monitor (GBM) observations around GW150914 uncovered a 1 s long gamma-ray signal (designated GW150914-GBM), 0.4 s after GW150914, having a broadly consistent sky location Connaughton et al. (2016, 2018). This event did not trigger GBM, but was found in an off-line search (e.g. Blackburn et al., 2015; Goldstein et al., 2016; Kocevski et al., 2018; Goldstein et al., 2019). Having been identified using a hard spectral template and lasting $\sim$ 1 s, GW150914-GBM is consistent with a short GRB. The chance association between GW150914 and GW150914-GBM (P=$2.2\times 10^{-3}$, Connaughton et al. (2016)) is not low enough to claim common origin for the two signals. Fermi-GBM routinely follows up BBH merger events. So far there has been no unambiguous gamma-ray counterpart detection (e.g. Goldstein, 2017; Burns, 2015; Hamburg, 2017).

An EM counterpart to a stellar mass binary BH (BBH) merger is not commonly expected, however there are a wealth of proposed mechanisms that could at least in principle provide detectable emission. Further observations in GWs and gamma-rays are needed to confirm or render such an association unlikely. In this paper we calculate the ratio of BBH mergers to gamma-ray counterparts for Fermi-GBM in different scenarios, assuming the GW150914 association is real. We calculate concrete estimates for future BBH observations during the ongoing LIGO-Virgo third observing run (O3).

A large number of scenarios have been outlined for producing observable gamma-rays from BBH mergers. Loeb (2016) sketches a scenario where a rapidly rotating core of a massive star collapses into two BHs that merge and accrete the remaining stellar material and produce GW150914-GBM. Dai et al. (2017) argues that the heat associated with dynamical friction for the BBH would remove the star as a potential accretion source, unless the BBH forms close to the center. In this case the dynamical friction speeds up the merger which has observable imprints on the waveform (Fedrow et al., 2017). Perna et al. (2016) propose a disk around one of the BHs that survives up to the time of the merger and it is subsequently accreted. This scenario was also discussed by Woosley (2016) and challenged by Kimura et al. (2017). Zhang (2016) presents a scenario, where charged BHs orbit each other before merging, Lyutikov (2016) argues however that the required charge is prohibitive. BBHs associated with GRBs can also be used to constrain GRB emission mechanisms (Veres et al., 2016).

In considering the possible reasons why other observed BBH mergers did not produce gamma-rays, aside from the uncertainty on the physical mechanism, we note that GW150914 is one of the most massive BBH merger observed. Together with the unfavorable position for Fermi-GBM, GW150914-GBM was observed close to the sensitivity limit. If there is a positive correlation between the total mass and the EM output, it is not surprising that there are no EM counterparts observed for lower mass or more distant BBH mergers.

Here we attempt to answer the following question: if GW150914 and GW150914-GBM were indeed associated, how many LIGO-Virgo BBH mergers will Fermi-GBM have to follow-up to detect a second source. We perform simulated observations of BBH mergers with LIGO-Virgo, adopt different scenarios for gamma-ray emission and calculate the ratio of BBH mergers detected by LIGO-Virgo to the number of gamma-ray counterpart detections by Fermi-GBM. An alternative formulation of the above question is: given an emission model, how many BBH follow-up observations (that are non-detections) are needed to rule out the specific model for GW150914-GBM.

The exposure function of GBM depends on a number of factors. Any given point on the sky will be observed roughly 60-70% of the time. Here we take 65% as an average. The two main reasons of the loss of coverage are: the source being occulted by Earth and the fact that Fermi is turned off during passage through the South Atlantic Anomaly.

We structure the paper in the following way: in Section 2 we describe the simulation of BBHs, in Section 3 we describe possible scenarios by which gamma-rays and GWs can be related. We discuss our results in Section 4 and conclude in Section 5. Physical constants have their usual meaning. For cosmological calculations we use $\Omega_{m}=1-\Omega_{\Lambda}=0.286$ and $H_{0}=69.6~{}{\rm km}{\hbox{~{}s}}^{-1}\hbox{~{}Mpc}^{-1}$ (Bennett et al., 2014).

2 Binary black hole merger simulations

2.1 Mass and spin assumptions

We simulate the primary mass of the BBH systems ($M_{1}$) following a power law distribution (Kovetz et al., 2017) $P(M_{1})\propto M_{1}^{\alpha_{M}}$ between $M_{\rm min}=5M_{\odot}$ and $M_{\rm max}=100M_{\odot}$ and $\alpha_{M}=-2.35$ from estimates in Abbott et al. (2017b). At high masses the distribution is cut by an exponential function, so the differential distribution is $dN/dM_{1}=C\left(M_{1}/5M_{\odot}\right)^{\alpha_{M}}\exp{(-M_{1}/M_{\rm cut})}$ where C is a constant, and $M_{\rm cut}=60M_{\odot}$. The secondary mass ($M_{2}$) follows a uniform distribution between $M_{\rm min}$ and M1.

The spins of individual BHs in the first Gravitaitonal Wave Transient Catalog (GWTC-1, The LIGO Scientific Collaboration et al., 2018) are not extreme, and thus they will not have a dominant effect in determining the spin parameter of the final BH. Based on this observation, we use non-spinning BH components in our simulations. In addition, the distribution of final spins in from our simulations is consistent with the final spins presented in GWTC-1.

2.2 Spatial distribution

We generate a population of binary black hole systems and analyze their LIGO-Virgo signal. The spatial distribution of BBHs is uniform in volume up to a distance Dmax. Dmax is determined from the requirement that all systems up to this distance with optimal viewing properties (face-on or inclination angle, $\iota=0$) should be detectable at signal to noise ratio (SNR, see Section 2.3) levels corresponding to at least borderline detection (e.g. SNR$>$8). This guarantees that we account for every BBH system potentially detectable by LIGO-Virgo for a given sensitivity. Within 4 Gpc, a BBH at the high end of the mass distribution with a favorable inclination, will be detected marginally. We thus fix D${}_{\rm max}=$ 4 Gpc.

We assign a total angular momentum vector to every BBH pointing at random directions. This will correspond to the direction the jet is launched, and also defines the inclination angle, $\iota$ as seen from Earth. Any non-zero individual spins would have an effect on the direction total angular momentum.

2.3 SNR calculation

We use the LALsim-inspiral tool111https://www.lsc-group.phys.uwm.edu/daswg/docs/howto/lal-install.html and generate frequency domain waveforms for the simulated BBH mergers (Jiménez-Forteza et al., 2017; Husa et al., 2016; Khan et al., 2016). BBHs are uniformly distributed in a $D_{\rm max}=$4 Gpc radius sphere. The masses that enter LALsim-inspiral are $M_{\rm det}=M_{\rm source}(1+z)$, where $M_{\rm source}$ are the simulated masses and $z$ is the redshift. We take the generated $h_{+}$ and $h_{\times}$ waveforms, and using equations 17-18 of Schutz (2011), we convolve them with the antenna responses $F_{+}$ and $F_{\times}$ for Hanford and Livingston to obtain the simulated SNR:

[TABLE]

We integrate the waveform from $f_{\rm low}=20\hbox{~{}Hz}$ and take the lower of the two instruments’ SNR. If it is above 8, we consider it a detection. Alternatively we consider SNR of 12 as a more conservative detection threshold. Past and future estimated sensitivity ($S_{h}(f)$) of LIGO are compiled in (Abbott et al., 2018), we use the ’Mid high/Late low’ sensitivity for both Hanford and Livingston as a proxy for O3 sensitivity.

We need a large sample of high SNR GW detections to have meaningful statistics. E.g. for $\approx 1000$ BBHs with SNR $>$ 8 with the mass distribution given in Section 2.1 out to Dmax= 4 Gpc, we need $\sim$400,000 simulated BBHs. Because most of these systems will not contribute above SNR = 8, in order to speed up the calculation, we do a preliminary SNR calculation (SNRS) using average antenna functions and a simplified waveform. For a BBH with total mass $M_{t}$, and inclination angle $\iota$ to our line of sight, the signal to noise ratio measured by LIGO-Virgo can be approximated as (Dalal et al., 2006; O’Leary et al., 2016):

[TABLE]

where $\mathcal{A}=\sqrt{\frac{5}{96}}\pi^{-2/3}\left(\frac{GM_{t}}{c^{3}}\right)^{5/6}c$, $D_{L}$ is the luminosity distance, and $F_{+}$ and $F_{\times}$ are taken to be their square-averaged value 1/5. $f_{\rm ISCO}(M_{t})=c^{3}6^{-3/2}\pi^{-1}G^{-1}M_{t}^{-1}\approx$ $438(M/10M_{\odot})^{-1}\hbox{~{}Hz}$ is the GW frequency at the innermost stable circular orbit around the final BH.

First we calculate the simple waveform and if it indicates an SNR${}_{\rm S}>$5, we subsequently calculate the SNR using the advanced waveform. By this method we only need to calculate 5% of the total cases using the advanced waveform. The chance of having an SNR${}_{\rm S}<5$ but SNR$>8$ is less than $10^{-3}$.

The probability that an observed GW from a compact binary merger will have an inclination $\iota$ is (Schutz, 2011):

[TABLE]

We plot this distribution in Figure 1 (red curve), overlaid with the distribution of simulated GW signals as a function of inclination angle (green histogram). The figure illustrates that our simulations capture the expected inclination dependence of GW detections.

2.4 GBM sensitivity

The flux measured by Fermi-GBM for GW150914-GBM   in the 10-1000 keV interval is $F\approx 2.4\times 10^{-7}\hbox{~{}erg}\hbox{~{}cm}^{-2}{\hbox{~{}s}}^{-1}$. The sensitivity of Fermi-GBM varies as a function of the spacecraft’s geographic location, the source position in spacecraft coordinates (Meegan et al., 2009) and other unmodeled factors.

To estimate the sensitivity of GBM for a GRB with a particular flux, we use the following method: We consider the fluxes of triggered GRBs in the Fermi-GBM catalog (Gruber et al., 2014; Bhat et al., 2016) and analyze the low-end tail of the distribution. This means fluxes typically below $3\times 10^{-7}\hbox{~{}erg}\hbox{~{}cm}^{-2}{\hbox{~{}s}}^{-1}$. GRBs below this flux level can potentially go undetected for GBM even if they occur within the field of view. The lowest flux where a GRB was triggered is $\sim 7\times 10^{-8}\hbox{~{}erg}\hbox{~{}cm}^{-2}{\hbox{~{}s}}^{-1}$, while the flux above which essentially all GRBs that are within GBMs field of view are detected is $\sim 3\times 10^{-7}\hbox{~{}erg}\hbox{~{}cm}^{-2}{\hbox{~{}s}}^{-1}$. In this flux interval we fit the distribution with a one-sided Gaussian (i.e defined only below the peak) and find that such a function gives a good description of the distribution of fluxes. Using this function we employ an acceptance-rejection method to decide if a GRB with a simulated flux is detected. We illustrate this function as a gradient below the histogram in Figures 2 through 5. Simulated flux values falling on darker hues have a larger chance for acceptance. The blue dashed curve marks the 50% acceptance rate.

2.5 Opening angles

GRBs involve jetted emission (Rhoads, 1999; Harrison et al., 1999). The energetics of the GRBs strongly depend on the value of the opening angle. Jet opening angles are measured from the achromatic break in the afterglow lightcurve of GRBs. The break occurs when the increasing angular size observed ($\approx 1/\Gamma$, where $\Gamma$ is the bulk Lorentz factor of the outflow) of the emitting surface grows larger the angular size ($\theta_{j}$). This results in a dearth of emitting surface and thus a steeper fading. There is a better handle on the opening angles in long GRBs, but short GRBs have an increasing number of measured jets (Fong et al., 2015).

In our simulations, the GRB that accompanies the GW signal is jetted with a half opening angle $\theta_{j}$. If our line of sight to the jet axis ($\theta_{v}$) is less than $\theta_{j}$, we may in principle observe the GRB, subject only to the instrument’s sensitivity and observing conditions.

It is unclear how the opening angles of BBH counterparts may be distributed. Here we discuss 3 cases: a fixed narrow 20∘ opening angle; a wide jet with $\theta_{j}=90^{\circ}$ opening angle corresponding to isotropic emission; a case where we associate uniformly random opening angles between 10∘ and 40∘. For simplicity, we only consider top-hat jets, without any angular structure. In other words, we expect to detect GRBs for which the observer lies within the jet and no flux is expected for observers outside of the jet.

For the observed jet opening angles, the uniform 10-40*∘* distribution is the most realistic. In addition, looking at the reported SNR in Table 1, the SNR$>$12 requirement is more in line with GW detections. We thus consider this our benchmark case and highlight column 5 in Table 2.

3 Possible paths to link GW and EM energy for binary black hole mergers

It is unclear if BBH mergers are accompanied by significant EM radiation, and indeed the most likely scenario is that they occur in very low density environments where it is difficult to extract e.g. the rotational energy of the BHs and channel it to gamma-rays.

We note that all the outlined scenarios to extract gamma-rays from a BBH merger suffer from non-trivial shortcomings, or critiques (Lyutikov, 2016). We consider the enumerated processes only as a guide to give an order of magnitude answer to the fundamental question we are seeking to answer: after how many BBH merger detections during O3, should we expect the next gamma-ray counterpart, assuming a particular scenario for GW150914-GBM.

We use the isotropic-equivalent gamma-ray luminosity of GW150914-GBM, $L_{\gamma,{\rm iso}}=1.8\times 10^{49}\hbox{~{}erg}{\hbox{~{}s}}^{-1}$ (Connaughton et al., 2016) as a scaling factor to guide our models. Some scenarios involve the final spin and the radiated energy in form of GWs ($E_{\rm GW}$). We use the formalism of Jiménez-Forteza et al. (2017) to calculate $E_{\rm GW}$ and the final spin parameter $a_{f}$. In this approach, these quantities are determined via polynomial functions based on the results of detailed general-relativistic simulations, e.g. Husa et al. (2016) (see also Buonanno et al., 2008, for an alternative approach). For our simulations, we use the formalism of Jiménez-Forteza et al. (2017).

As indicated before, the gamma-ray luminosity of GW150914-GBM is an isotropic-equivalent luminosity. The physical models considered in this section however yield the total available luminosity. We account for the beaming correction between the two quantities by assuming the opening angle of GW150914-GBM is representative of the simulated BBH population. Thus when normalizing any given model to the flux of GW150914-GBM, the correction between the isotropic-equivalent (measured) luminosity and the total (model) luminosity is included in the scaling factor (the gamma-ray production efficiency is accounted for in a similar way). The fact that GW150914-GBM’s opening angle may not be representative for the BBH population is one of the reasons we consider pessimistic and optimistic cases (Section 3.6). This is a possible factor that can drive a dimmer or brighter expected population, given GW150914-GBM and a model.

3.1 Neutrino driven wind

Neutrino driven winds are routinely invoked as mechanisms launching the GRB jet (e.g. Ruffert & Janka, 1998). It implies a dense accretion disk that emits neutrinos towards the axial region of the system where $\nu\bar{\nu}$ collisions result in e± pairs that in turn drive a relativistic jet. Due to the high accretion rate required, Li et al. (2016) find that this mechanism is unlikely to be at work in the case of GW150914-GBM. Nonetheless, since this is one of the leading jet launching mechanisms, we consider the implication of this scenario here.

Zalamea & Beloborodov (2011) finds an empirical relation for the neutrino driven wind luminosity:

[TABLE]

$f_{3}(a)=3+Z_{2}-[(3-Z_{1})(3+Z_{1}+2Z_{2})]^{1/2}$ and $Z_{1}=1+(2-a^{2})^{1/3}[(1+a)^{1/3}+(1-a)^{1/3}]$ and $Z_{2}=(3a^{2}+Z_{1}^{2})^{1/2}$. In this scenario, we vary the uncertain accretion rate $\dot{M}$ and scale it to match the luminosity of GW150914-GBM (additional details in Section 3.6). The results of simulations using this model are illustrated in Figure 2

3.2 Blandford Znajek (BZ) mechanism

The BZ mechanism (Blandford & Znajek, 1977) extracts the rotational energy ($E_{\rm rot}=M_{f}c^{2}\{1-\left[(1+(1-a_{f}^{2})^{1/2})/2\right]^{1/2}\}$) of the BH with the help of magnetic fields, where $a_{f}$ is the spin parameter of the merged BH. Here we use the formulae by Komissarov & Barkov (2010); Reynolds et al. (2006) for the extracted EM luminosity:

[TABLE]

where $\Omega_{h}\approx c^{3}/GM_{f}$ is the angular velocity of the BH, and $\Psi_{h}\approx 2\pi R_{g}^{2}B$ the magnetic flux threading the BH horizon at one of its hemispheres, $R_{g}=GM_{f}/c^{2}$ is the gravitational radius, B the magnetic field strength and $f(a)=f_{1}^{4}(a)f_{2}^{2}(a)$ ($f_{1}(a)=2-a+2(1-a)^{1/2}$ and $f_{2}(a)=a/2(1+(1-a^{2})^{1/2})$) accounts for the approximations in the expressions of $\Psi_{h}$ and $\Omega_{h}$. The main unknown here is the value of the magnetic field. We set this value so that the simulations scale with the GW150914-GBM   observations (see Section 3.6 for a description of the scalings). Figure 3 shows the flux distribution obtained from this model.

3.3 Charged black holes

Zhang (2016) proposed a mechanism, where at least one of the BHs carries a significant charge and Poynting flux is extracted from the system (see also Liu et al., 2016, for a similar approach). To approximate the luminosity in this scenario, we consider equation (7) of Zhang (2016):

[TABLE]

where C=49/120000 a numerical factor, $\hat{q}$ is a dimensionless charge, in units of critical charge, $Q_{c}$ ($Q_{c}=2\sqrt{G}M_{f}$) and $\hat{a}$ is the normalized Newtonian distance of the two BHs, defined as $\hat{a}=(1+\sqrt{1-a^{2}})/2$. We scale the luminosity to GW150914-GBM and vary the dimensionless charge to obtain the gamma-ray flux of the simulated population.

3.4 Gamma-rays as a fraction of GW energy release

During the BBH merger a fraction of the mass is radiated in form of GWs. While the coupling between GW and matter is extremely weak, here we explore a scenario where the gamma-ray output of the system depends on the GW energy output. We do not imply that the GW energy is channeled into gamma-rays, but the gamma-ray generation correlates with the emitted GW energy. We can parametrize the gamma-ray luminosity as a fraction of $E_{\rm GW}$, $L_{\gamma}=\epsilon_{\rm GW}E_{\rm GW}$. Abbott et al. (2016c) estimates the total energy emitted in GWs: $E_{\rm GW}=5.3\times 10^{54}\hbox{~{}erg}$. To normalize the GW energy to e.g. GW150914-GBM, we have $\epsilon_{\rm GW}=3.4\times 10^{-6}{\hbox{~{}s}}^{-1}$. We refer to this model as ’gamma-GW fraction’ and show the associated gamma-ray flux distribution in Figure 5.

3.5 Power law dependence on the final mass

For this scenario, we assume the EM luminosity, $L_{\gamma}$ has a power law dependence on the final BH mass, $M_{f}$: $L_{\gamma}=L_{*}(M_{f}/M_{*})^{p}$. We may use e.g. $L_{*}=1.8\times 10^{49}\hbox{~{}erg}{\hbox{~{}s}}^{-1}$ and $M_{*}=62M_{\odot}$ for normalization.

This scenario is best suited to illustrate the scaling with final mass for the gamma-ray flux (see Figure 6). The final mass of GW150914 was at the high end of the detected mass distribution and the source was relatively nearby. If $L_{\gamma}$ correlates positively with $M_{f}$ or $p>0$, we expect gamma-ray counterparts to be a rare occurrence. If however the gamma-ray emission is negatively correlated with the final BH mass ($p<0$), that would suggest that GW150914-GBM   was a common event.

3.6 Input for simulations

For each emission model described above, we consider the following subcases:

*Pessimistic: * In this case we simulate the unmodeled parameters (magnetic field, $B$ for BZ; accretion rate, $\dot{M}$ for neutrino annihilation; normalized charge $\hat{q}$ for the charged BHs, etc) for each scenario according to a log-normal distribution with a mean that is half order of magnitude ($10^{1/2}$) lower than what is required to produce a GW150914-GBM-like signal, and has a standard deviation of half an order of magnitude ($10^{1/2}$). 2. 2.

Normalized to GW150914-GBM: The parameters of the scenario are chosen so that they reproduce on average, the Fermi-GBM luminosity of GW150914-GBM, with a standard deviation of half an order of magnitude. 3. 3.

*Optimistic: * This is similar to the pessimistic case, but the mean of the log-normal distribution is half an order of magnitude higher than the one needed to reproduce GW150914-GBM. 4. 4.

In addition to the above three cases, we also carry out simulations where the gamma-ray luminosity is required to match exactly GW150914-GBM. The purpose such a scenario is to gauge if currently observed BBHs would be expected to produce a detectable gamma-ray signal (for the BBH-to-GRB ratio ratio in these cases see Figures 2, through 4, and Section 3.8).

3.7 Gamma-rays from a random, unknown process

It is possible that a hereto unknown process is responsible for GW150914-GBM, that is unrelated to the GW observables (see e.g. Bartos et al., 2017, for BBH mergers occurring in active galactic nuclei). We explore this scenario by randomly assigning gamma-ray luminosity to simulated BBH mergers. We constrain the mean and standard deviation of the flux to the 3 scenarios outlined in Section 3.6 (normalized, pessimistic and optimistic), but independent of any physical property of the systems.

3.8 Application to currently observed BBH mergers and GBM upper limits

There are 10 BBH mergers detected prior to the start of O3 (The LIGO Scientific Collaboration et al., 2018), with only one putative gamma-ray counterpart. We note that in this work we used the O3 sensitivity for the simulations that is superior to the previous two observing run’s sensitivity. The BBH-to-GRB ratio is expected to be lower for BBHs detected in a less sensitive configuration, because of the lower average distance of the sources imply higher gamma-ray flux. At the same time the expected number of events is also lower.

Gamma-ray emission will likely be beamed, while the GW emission is close to isotropic. Thus for individual events it is difficult to rule out any particular model. Even if a model suggests a large flux, it is possible that the jet was beamed away from Earth, thus undetectable for GBM. This could be the case e.g. for GW170814 (7), GW151226 (2) and GW170608 (4) in Figure 2. The neutrino driven wind model predicts a detectable flux for these events. Interestingly these three BBH mergers imply fluxes within the sensitivity of GBM for the gamma-GW fraction model (Figure 5), and the charged BH scenario (Figure 4) as well. For the BZ model, only GW170814 (7) produces detectable gamma-ray flux (Figure 3).

Fermi GBM reported on a weak source following GW170814 (7) (Goldstein, 2017), but it was deemed unassociated with GW170814, based on the small overlap between source locations. Part of the location region of GW151226 (2) was occulted for Fermi and no significant sources were detected (Burns, 2015). The location region of GW170608 (4) was covered well by GBM and no sources were reported (Hamburg, 2017).

4 Discussion

We have calculated the expected LIGO-Virgo detected BBH mergers to gamma-ray counterparts (BBH-to-GRB ratio) that could be observed by Fermi-GBM in different scenarios. The results of this study are summarized in Table 2.

The magnitude of the gamma-ray flux depends on the adopted scenario. Among the BBHs detected by LIGO and Virgo GW150914 has a relatively large final mass and it is located close compared to other GW events. It is thus in line with expectations, that in scenarios where the gamma-ray luminosity positively correlates with the mass of the final BH we expect a large number of BBH observations before another counterpart is detected. Conversely, for the scenarios where the EM luminosity is inversely proportional with the final mass, the required number of BBH observations before another gamma-ray counterpart is observed is low ($\lesssim$ 10). This is displayed in Figure 6: for negative $p$ values, essentially all the generated GRB flux that reaches the observer is detected. The BBH-to-GRB ratio is governed by the opening angle distribution. As we move to positive $p$ values, the ratio increases because we scale our calculations to GW150914, involving a nearby, massive BH. For $p>0$ the bulk of the simulated BHs will be less massive thus generating lower gamma-ray flux.

Figures 2 through 5 illustrate the calculation of the BBH-to-GRB ratio in the cases where the assumed mechanism reproduces GW150914-GBM exactly. The green histogram shows the gamma-ray flux for all the detectable BBHs with SNR 12 or 8 for LIGO-Virgo. The initial cut for the jets that point elsewhere reduce the numbers to the gray histogram. Considering the Fermi-GBM live-time, sky exposure and detection threshold further reduces the detected GRBs to the red histogram. For each scenario, we show the expected fluxes from the observed BBH population as well, indicated by numbers on the histograms.

4.1 Favored and unfavored mechanisms

We find that the neutrino driven wind mechanism gives the lowest BBH-to-GRB ratio (comparable to the power low scenario with $p=-2$). This can be understood as we scale the gamma-ray flux to the relatively nearby and high mass GW150914 event. In this scenario there is a strong dependence on the final mass (Equation 4). Even in the pessimistic scenario we can constrain this mechanism (rule it out) after $\sim$15 BBH observations.

If the EM output depends on the final BH mass as a power law, our simulations show that for positive power law indices ($p$, where $L_{\gamma}\propto M_{f}^{p}$) the BBH-to-GRB ratio ranges from $\sim$ 5 to a few times 10 (for $p=2$ and SNR$>$12) or to $\lesssim$ 200 (for $p=2$ and SNR$>$8) as shown in Figure 6.

4.2 Observed BBHs

The neutrino driven wind scenario results in 3 (out of 9) of the observed GWs with confidently detectable flux for GBM (Figure 2) while the simulations give a BBH-to-GRB ratio of 7.7. The charged BH scenario gives similarly 3 GWs in the detectable flux range for GBM, albeit with lower flux (see Figure 4). The gamma-GW fraction scenario places 3 GW in the marginally detectable flux regime, but well below the 50 % detection threshold (see Figure 5). The BZ scenario results in detectable counterpart for a single known GW (see Figure 3). Taken at face value these results indicate that with an increasing number of non-detections, the neutrino driven wind and the charged BH model can be ruled out first, then the gamma-GW fraction scenario. Finally the BZ mechanism requires the most non-detections for it to be ruled out.

4.3 Role of the opening angle and signal to noise limit

The assumed opening angle significantly change the resulting BBH-to-GRB ratio: wider jets are more likely to include the detector in their aperture thus decreasing the BBH-to-GRB ratio. We focused the discussion so far on the benchmark case of uniformly distributed opening angles between 10 and 40*∘*. We have also calculated a narrow (20 degree) and wide, 90 degrees case as well, to get a sense of how the number of GRBs will change.

The narrow opening angle increases BBH-to-GRB ratio compared to the benchmark case by anywhere between 50% to 100%. The wide, isotropic emission case decreases BBH-to-GRB ratio by a factor of 3 and up to 20.

Allowing for lower SNR (¿8) threshold has the effect of increasing the number of GWs detected and the distance limit. This in turn disproportionately increases the number of gamma-ray events with flux below the detection threshold, thus increasing the BBH-to-GRB ratio. The increase in BBH-to-GRB ratio can range from a few percent to a factor of 3 (see Table 2).

4.4 Triggers versus off-line searches

It is important to note that the GBM sensitivity considered here applies to GRBs that were triggered real-time on-board the spacecraft. GW150914-GBM was found in an off-line search. Off-line searches improve the sensitivity by a factor of few (Kocevski et al., 2018), thus improving (decreasing) the BBH-to-GRB ratio.

4.5 Expected distribution of BBH-to-GRB ratio

A model predicts on average $N_{\rm BGR}$ BBH observations for each GRB. The probability of observing [math] GRBs after $N_{\rm BBH}$ BBH detections for such a model, can be calculated from the binomial distribution:

[TABLE]

We illustrate the implications for the models by calculating this probability after exactly the same number of observed BBHs as the model’s BBH-to-GRB ratio would suggest ($N_{\rm BBH}=N_{\rm BGR}$). Using representative cases, e.g. $N_{\rm BGR}=\{5,10,20,100\}$, yields the following probabilities: $P(x=0|N_{\rm BBH}=N_{\rm BGR})=\{0.33\,,0.35,\,0.36,\,0.37\}$. In other words, after $N_{\rm BGR}$ non-detections there is still a $\sim 35$% probability that the model is correct and will lead to a BBH observation with a counterpart.

After $N_{\rm BBH}>N_{\rm BGR}$ observations with no GRB counterparts, what is the probability that the model predicting $N_{\rm BGR}$ BBHs for each GRB on average is correct? Equation 7 can be easily inverted to calculate the number of BBH observations necessary to rule out a model predicting $N_{\rm BGR}$, with a desired confidence $P_{\rm conf}$:

[TABLE]

For example, a model predicting BBH-to-GRB ratio of $N_{\rm BGR}=10$ can be ruled out with $P_{\rm conf}=10^{-2}$ confidence after $\approx 44$ non-detections. Figure 7 shows the required number of BBH observations to rule out a model with an arbitrary BBH-to-GRB ratio of $N_{\rm BGR}$ with a given confidence.

4.6 Detection prospects

Assuming a rate of 40 BBH detections during the year long observing run O3 (The LIGO Scientific Collaboration et al., 2018; Abbott et al., 2018), and with no counterpart detection, we can make the following statements for the benchmark case ($\theta_{j}$ uniformly distributed random between 10-40∘ and SNR¿12, highlighted in Table 2): The neutrino driven wind, the power law model with $p=-2,0$, and the random model can be ruled out even in the pessimistic case. The BZ model, the charged BH, the gamma-GW fraction models can be constrained in all but the pessimistic cases.

Assuming again the benchmark case, the majority of the scenarios point to 6 to 20 BBHs for every gamma-ray counterpart detected. This is encouraging for the third observing run of LIGO-Virgo, that is planned to last a year, and is predicted to observe about 40 BBH mergers. If no counterpart is detected during O3, most of the benchmark models can be ruled out. The only feasible models remaining after O3 are those that assume that GW150914-GBM was a brighter than usual event given a particular gamma-ray producing mechanism (pessimistic variants of e.g. Blandford-Znajek, charged BH or gamma-GW fraction models).

After 40 BBH observations with no counterparts, almost all isotropic emission models ($\theta_{j}=90^{\circ}$) can be ruled out with reasonable confidence. For narrowly beamed gamma-ray emission ($\theta_{j}=20^{\circ}$), most models can be ruled out, but some can only be constrained.

5 Conclusion

In this paper we investigated detection prospects of gamma-ray emission associated with BBH mergers, using different emission models. We emphasise here the caveat that EM counterparts from stellar mass binary BH mergers are unexpected and all the outlined models have considerable issues to be worked out. Nonetheless, the detection of GW150914-GBM presents and intriguing prospect that through some mechanism, not considered before, significant energy might be channeled from GW to gamma-rays in the process of a BBH merger.

Assuming GW150914-GBM was a GRB associated with GW150914, we quantitatively estimated how many BBH events Fermi-GBM should follow up in order to detect a second one with a counterpart. We find that for a majority of the models we expect 6-20 BBH mergers for every gamma-ray counterpart. Hence we expect that, after the ongoing third LIGO-Virgo observing run, either another counterpart will be found, or many of the models discussed here will be incompatible with the data.

**Acknowledgements: ** We thank Bing Zhang for valuable input, Francesco Pannarale for comments on the manuscript and the anonymous referee for a valuable comment. PV acknowledges support from Fermi grant NNM11AA01A and 80NSSC17K0750. TD and EB are supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by Universities Space Research Association under contract with NASA. NC is supported by NSF Grants PHY-1806990 and PHY-1505373.

This paper has LIGO document number LIGO-P1900147.