<figure><img src="image_1.png" /></figure>

# Fermi-GBM Follow-up of LIGO-Virgo Binary Black Hole Mergers: Detection Prospects

P. Veres¹, T. Dal Canton²,<fnref n="7" />,⁹, E. Burns²,<fnref n="7" />, A. Goldstein³, N. Christensen⁵,⁶, and R. D. Preece¹

¹Center for Space Plasma and Aeronomic Research, University of Alabama in Huntsville, 320 Sparkman Drive, Huntsville, AL 35899, USA

²NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

³Science and Technology Institute, Universities Space Research Association, Huntsville, AL 35805, USA

⁴NASA Marshall Space Flight Center, Huntsville, AL 35811, USA

⁵Carleton College, Northfield, MN 55057, USA

⁶Artemis, Université Côte d'Azur, Observatoire Côte d'Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France

Received 2019 May 20; revised 2019 July 10; accepted 2019 July 10; published 2019 September 2

## 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 (BH) 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 the unfavorable location of Fermi-GBM. Here we answer the following question: if GW150914 and GW150914-GBM were associated, how many LIGO-Virgo binary BH mergers would Fermi-GBM have to follow up to detect a second source? To answer this question, we perform simulated observations of binary BH mergers with LIGO-Virgo and adopt different scenarios for gamma-ray emission from the literature. We calculate the ratio of simulated binary BH mergers detected by LIGO-Virgo to the number of gamma-ray counterpart detections by Fermi-GBM, the BBH-to-GRB ratio. A large majority of the models considered here predict a BBH-to-GRB ratio in the range of 5–20, but for optimistic cases it can be as low as 2, while for pessimistic assumptions it can be as high as 700. Hence, we expect that the third observing run, with its high rate of binary BH detections and assuming the absence of a joint detection, will provide strong constraints on the presented models.

**Key words:** gamma rays: general – gamma-ray burst: individual (GW150914-GBM) – gravitational waves

## 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), providing a new view of the universe. There is an intense, ongoing effort to detect electromagnetic (EM) counterparts for GW observations (Abbott et al. 2016a, 2017a, 2019; Racusin et al. 2017; Burns et al. 2019). This effort resulted in the detection of the first unambiguous gamma-ray counterpart to a GW signal (Abbott et al. 2017a, 2017c; Goldstein et al. 2017), from a binary neutron star merger.

Fermi-Gamma-ray Burst Monitor (GBM) observations around GW150914 uncovered a 1 s 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 offline search (e.g., Blackburn et al. 2015; Goldstein et al. 2016, 2019; Kocevski et al. 2018). Having been identified using a hard spectral template and lasting ~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., Burns 2015; Goldstein 2017; 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

<footnote n="7">⁷ NASA Postdoctoral Program fellow.</footnote>
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 \text{ km s}^{-1} \text{ Mpc}^{-1}$ (Bennett et al. 2014).

## 2. Binary BH 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_{\min} = 5M_\odot$ and $M_{\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(M_1/5M_\odot)^{\alpha_M} \exp(-M_1/M_{\text{cut}})$, where C is a constant, and $M_{\text{cut}} = 60M_\odot$. The secondary mass ($M_2$) follows a uniform distribution between $M_{\min}$ and $M_1$.

The spins of individual BHs in the first Gravitational 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 from our simulations is consistent with the final spins presented in GWTC-1.

### 2.2. Spatial Distribution

We generate a population of binary BH systems and analyze their LIGO-Virgo signal. The spatial distribution of BBHs is uniform in volume up to a distance $D_{\max}$. $D_{\max}$ 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 (S/N; see Section 2.3) levels corresponding to at least borderline detection (e.g., S/N > 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_{\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. S/N Calculation

We use the LALsim-inspiral tool<fnref n="8" /> and generate frequency domain waveforms for the simulated BBH mergers (Husa et al. 2016; Khan et al. 2016; Jiménez-Forteza et al. 2017). BBHs are uniformly distributed in a $D_{\max} = 4$ Gpc radius sphere. The masses that enter LALsim-inspiral are $M_{\text{det}} = M_{\text{source}}(1+z)$, where $M_{\text{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 S/N:

$$S/N^2 = 4 \int \frac{|h_+ F_+ + h_\times F_\times|^2}{S_h(f)} df. \quad (1)$$

We integrate the waveform from $f_{\text{low}} = 20$ Hz and take the lower of the two instruments' S/N. If it is above 8, we consider it a detection. Alternatively, we consider an S/N 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-S/N GW detections to have meaningful statistics, e.g., for $\approx 1000$ BBHs with S/N > 8 with the mass distribution given in Section 2.1 out to $D_{\max} = 4$ Gpc, we need $\sim 400,000$ simulated BBHs. Because most of these systems will not contribute above S/N = 8, in order to speed up the calculation, we do a preliminary S/N calculation (S/Ns) using average antenna functions and a simplified waveform. For a BBH with total mass $M_\iota$, and inclination angle $\iota$ to our line of sight, the S/N measured by LIGO-Virgo can be approximated as (Dalal et al. 2006; O’Leary et al. 2016)

$$S/N_S^2 = 4 \frac{A^2}{D_L^2} (\langle F_+^2 \rangle (1 + \cos^2 \iota)^2 + 4 \langle F_\times^2 \rangle \cos^2 \iota) \\ \times \int_{f_{\text{low}}}^{f_{\text{ISCO}}} \frac{f^{-7/3}}{S_h(f)} df, \quad (2)$$

where $A = \sqrt{\frac{5}{96}} \pi^{-2/3} (\frac{GM_\iota}{c^3})^{5/6} c$, $D_L$ is the luminosity distance, and $F_+$ and $F_\times$ are taken to be their square-averaged value, one-fifth. $f_{\text{ISCO}}(M_\iota) = c^3 6^{-3/2} \pi^{-1} G^{-1} M_\iota^{-1} \approx 438(M/10M_\odot)^{-1}$ 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 S/Ns > 5, we subsequently calculate the S/N 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 S/Ns < 5 but S/N > 8 is less than $10^{-3}$.

<footnote n="8">⁸ https://www.lsc-group.phys.uwm.edu/daswg/docs/howto/lal-install.html</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Distribution of inclination angles from ~8000 simulated BBHs that had a S/N > 8 out of ~164,000 calculated waveforms. The initial number of simulated BBHs was 4 × 10⁶, most of which did not contribute at S/N > 8. The red curve shows the expected distribution.</figcaption></figure>

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

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

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}$ erg cm$^{-2}$ 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}$ erg cm$^{-2}$ 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 ~7 × 10$^{-8}$ erg cm$^{-2}$ s$^{-1}$, while the flux above which essentially all GRBs that are within GBM’s field of view are detected is ~3 × 10$^{-7}$ erg cm$^{-2}$ 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.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Simulated EM counterparts based on a neutrino-driven wind model with GW S/N > 12 (see Section 3.1). The model was normalized to exactly match GW150914-GBM’s flux. The green histogram represents the simulated BBH mergers with LIGO-Virgo S/N > 12. The gray histogram shows the cases when the associated jet points toward Earth. The red histogram incorporates the sensitivity and the observing live-time of Fermi-GBM showing all the detections from the simulations. The band below the histogram illustrates the GBM detection probability, and the vertical dashed blue line indicates the 50% mark. Based on this model every 7.7 BBH observed by LIGO-Virgo should have a detectable EM counterpart by GBM. Out of the 10 observed BBH mergers, 4 (GW170814 (7), GW150914, GW151226 (2), and GW170608 (4)) yield fluxes above the GBM threshold (Section 4.2).</figcaption></figure>

## 2.5. Opening Angles

GRBs involve jetted emission (Harrison et al. 1999; Rhoads 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 light curve of GRBs. The break occurs when the increasing angular size observed (≈1/Γ, where Γ is the bulk Lorentz factor of the outflow) of the emitting surface grows larger than 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 three cases: a fixed narrow 20° opening angle; a wide jet with a $\theta_j = 90°$ 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 S/N in Table 1, the S/N > 12 requirement is more in line with GW detections. We thus consider this our benchmark case and highlight column 5 in Table 2.
<table><caption>**Table 1**
Central Values for the Observed BBH Merger Parameters, Taken from The LIGO Scientific Collaboration et al. (2018)</caption><thead><tr><th>BBH Name</th><th><i>M</i><sub>1</sub></th><th><i>M</i><sub>2</sub></th><th><i>M</i><sub>f</sub></th><th><i>a</i><sub>f</sub></th><th><i>D</i><sub>L</sub>/Mpc</th><th><i>z</i></th><th><i>E</i><sub>GW</sub>/<i>M</i><sub>☉</sub><i>c</i><sup>2</sup></th><th>S/N</th></tr></thead><tbody><tr><td>GW150914</td><td>35.6</td><td>30.6</td><td>63.1</td><td>0.69</td><td>430</td><td>0.09</td><td>3.1</td><td>23.7</td></tr><tr><td>GW151012 (1)</td><td>23.3</td><td>13.6</td><td>35.7</td><td>0.67</td><td>1060</td><td>0.21</td><td>1.5</td><td>9.7</td></tr><tr><td>GW151226 (2)</td><td>13.7</td><td>7.7</td><td>20.5</td><td>0.74</td><td>440</td><td>0.09</td><td>1.0</td><td>13.0</td></tr><tr><td>GW170104 (3)</td><td>31.0</td><td>20.1</td><td>49.1</td><td>0.66</td><td>960</td><td>0.19</td><td>2.2</td><td>13</td></tr><tr><td>GW170608 (4)</td><td>10.9</td><td>7.6</td><td>17.8</td><td>0.69</td><td>320</td><td>0.07</td><td>0.9</td><td>13</td></tr><tr><td>GW170729 (5)</td><td>50.6</td><td>34.3</td><td>80.3</td><td>0.81</td><td>2750</td><td>0.48</td><td>4.8</td><td>13</td></tr><tr><td>GW170809 (6)</td><td>35.2</td><td>23.8</td><td>56.4</td><td>0.70</td><td>990</td><td>0.20</td><td>2.7</td><td>13</td></tr><tr><td>GW170814 (7)</td><td>30.7</td><td>25.3</td><td>53.4</td><td>0.72</td><td>580</td><td>0.12</td><td>2.7</td><td>15</td></tr><tr><td>GW170818 (8)</td><td>35.5</td><td>26.8</td><td>59.8</td><td>0.67</td><td>1020</td><td>0.20</td><td>2.7</td><td>15</td></tr><tr><td>GW170823 (9)</td><td>39.6</td><td>29.4</td><td>65.6</td><td>0.71</td><td>1850</td><td>0.34</td><td>3.3</td><td>15</td></tr></tbody></table>

Note. Masses are in units of solar masses. The numbers next to GW events are for referencing the events in figures.

### 3. Possible Paths to Link GW and EM Energy for Binary BH 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_{γ,iso} = 1.8 \times 10^{49}$ erg 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 the form of GWs ($E_{GW}$). We use the formalism of Jiménez-Forteza et al. (2017) to calculate $E_{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 toward the axial region of the system where $ν\bar{ν}$ collisions result in $e^{\pm}$ 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, because this is one of the leading jet launching mechanisms, we consider the implications of this scenario here.

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

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

$$f_3(a) = 3 + Z_2 - \left[ (3 - Z_1)(3 + Z_1 + 2Z_2) \right]^{1/2} \quad \text{and}$$

$$Z_1 = 1 + (2 - a^2)^{1/3} \left[ (1 + a)^{1/3} + (1 - a)^{1/3} \right] \quad \text{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_{rot} = M_f c^2 \{1 - [(1 + (1 - a_f^2)^{1/2})/2]^{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:

$$L_{\text{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), \quad (5)$$

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$ is 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
<table><caption>**Table 2**
The BBH-to-GRB Ratio in Different Scenarios</caption><thead><tr><th rowspan="2">Scenario</th><th rowspan="2">Sections</th><th colspan="2">θ<sub>j</sub> = 20°</th><th colspan="2">θ<sub>j</sub> = Unif. [10, 40]°</th><th colspan="2">θ<sub>j</sub> = 90°</th></tr><tr><th>S/N > 8</th><th>S/N > 12</th><th>S/N > 8</th><th>S/N > 12</th><th>S/N > 8</th><th>S/N > 12</th></tr></thead><tbody><tr><td>Neutrino wind, pessimistic</td><td></td><td>31.1</td><td>29.1</td><td>20.7</td><td><strong>15.3</strong></td><td>5.6</td><td>3.9</td></tr><tr><td>Neutrino wind, normalized</td><td>3.1</td><td>15.7</td><td>14.6</td><td>10.4</td><td><strong>8.6</strong></td><td>2.8</td><td>2.5</td></tr><tr><td>Neutrino wind, optimistic</td><td></td><td>10.7</td><td>10.5</td><td>7.1</td><td><strong>6.5</strong></td><td>2.0</td><td>1.9</td></tr><tr><td>Blandford Znajek, pessimistic</td><td></td><td>304.7</td><td>181.9</td><td>146.7</td><td><strong>94.9</strong></td><td>33.4</td><td>24.0</td></tr><tr><td>Blandford Znajek, normalized</td><td>3.2</td><td>44.5</td><td>32.6</td><td>34.9</td><td><strong>21.8</strong></td><td>7.1</td><td>5.0</td></tr><tr><td>Blandford Znajek, optimistic</td><td></td><td>16.4</td><td>14.1</td><td>10.8</td><td><strong>9.0</strong></td><td>2.9</td><td>2.4</td></tr><tr><td>Charged BH, pessimistic</td><td></td><td>158.4</td><td>121.3</td><td>101.6</td><td><strong>66.2</strong></td><td>22.3</td><td>13.9</td></tr><tr><td>Charged BH, normalized</td><td>3.3</td><td>22.6</td><td>21.0</td><td>14.6</td><td><strong>11.7</strong></td><td>3.9</td><td>2.9</td></tr><tr><td>Charged BH, optimistic</td><td></td><td>10.0</td><td>9.7</td><td>6.5</td><td><strong>5.9</strong></td><td>1.9</td><td>1.8</td></tr><tr><td>L<sub>γ</sub> = E<sub>GW</sub> × ϵ<sub>GW</sub>, pessimistic</td><td></td><td>720.2</td><td>218.3</td><td>565.9</td><td><strong>121.3</strong></td><td>55.4</td><td>21.6</td></tr><tr><td>L<sub>γ</sub> = E<sub>GW</sub> × ϵ<sub>GW</sub>, normalized</td><td>3.4</td><td>176.0</td><td>62.4</td><td>99.0</td><td><strong>29.1</strong></td><td>11.7</td><td>5.1</td></tr><tr><td>L<sub>γ</sub> = E<sub>GW</sub> × ϵ<sub>GW</sub>, optimistic</td><td></td><td>31.7</td><td>17.1</td><td>18.3</td><td><strong>9.2</strong></td><td>3.7</td><td>2.2</td></tr><tr><td>L<sub>γ</sub> ∝ M<sub>f</sub><sup>-2.0</sup>, pessimistic</td><td></td><td>16.9</td><td>15.0</td><td>11.0</td><td><strong>8.7</strong></td><td>3.0</td><td>2.5</td></tr><tr><td>L<sub>γ</sub> ∝ M<sub>f</sub><sup>-2.0</sup>, normalized</td><td>3.5</td><td>13.4</td><td>12.3</td><td>9.1</td><td><strong>7.4</strong></td><td>2.5</td><td>2.2</td></tr><tr><td>L<sub>γ</sub> ∝ M<sub>f</sub><sup>-2.0</sup>, optimistic</td><td></td><td>11.9</td><td>10.9</td><td>7.8</td><td><strong>6.9</strong></td><td>2.2</td><td>2.0</td></tr><tr><td>L<sub>γ</sub> = 5.7 × 10<sup>48</sup> erg s<sup>-1</sup>(const.)</td><td></td><td>114.8</td><td>54.6</td><td>64.9</td><td><strong>28.7</strong></td><td>10.7</td><td>5.0</td></tr><tr><td>L<sub>γ</sub> = 1.8 × 10<sup>49</sup> erg s<sup>-1</sup>(const.)</td><td>3.5</td><td>31.1</td><td>18.3</td><td>19.6</td><td><strong>10.2</strong></td><td>4.2</td><td>2.7</td></tr><tr><td>L<sub>γ</sub> = 5.7 × 10<sup>49</sup> erg s<sup>-1</sup>(const.)</td><td></td><td>14.6</td><td>11.8</td><td>9.6</td><td><strong>7.2</strong></td><td>2.5</td><td>2.0</td></tr><tr><td>L<sub>γ</sub> ∝ random, pessimistic</td><td></td><td>79.2</td><td>42.8</td><td>49.2</td><td><strong>22.1</strong></td><td>8.8</td><td>4.8</td></tr><tr><td>L<sub>γ</sub> ∝ random, normalized</td><td>3.7</td><td>29.8</td><td>20.6</td><td>18.0</td><td><strong>11.2</strong></td><td>4.1</td><td>2.8</td></tr><tr><td>L<sub>γ</sub> ∝ random, optimistic</td><td></td><td>15.4</td><td>12.5</td><td>9.6</td><td><strong>7.6</strong></td><td>2.6</td><td>2.1</td></tr></tbody></table>

**Note.** We consider three types of opening angles: fixed narrow ($\theta_j = 20°$), distributed randomly between 10° and 40° and isotropic ($\theta_j = 90°$). We show two limits in GW S/N (8 and 12). The energy extraction scenarios are described in Section 3, while the three cases are presented in Section 3.6. Bolded values highlight the benchmark case, indicating the greatest similarity to short GRB opening angles and GW observations.

Section 3.6 for a description of the scalings). Figure 3 shows the flux distribution obtained from this model.

### 3.3. Charged BHs

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):

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

where $C = 49/120,000$ a numerical factor, $\hat{q}$ is a dimensionless charge, in units of critical charge, $Q_c (Q_c = 2\sqrt{GM_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 the 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 parameterize the gamma-ray luminosity as a fraction of $E_{GW}$, $L_\gamma = \epsilon_{GW}E_{GW}$. Abbott et al.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Expected gamma-ray flux distribution from BBH, arising from the Blandford Znajek mechanism (Section 3.2). The BBH-to-GRB ratio is 2183/24 = 91.0. The notations are similar to those of Figure 2.</figcaption></figure>

(2016c) estimates the total energy emitted in GWs: $E_{GW} = 5.3 \times 10^{54}$ erg. To normalize the GW energy to e.g., GW150914-GBM, we have $\epsilon_{GW} = 3.4 \times 10^{-6}$ s$^{-1}$. We refer to this model as the “gamma-GW fraction” and show the associated gamma-ray flux distribution in Figure 5.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Expected gamma-ray flux distribution from BBH, arising in the charged BH scenario (Section 3.3). Note that in this case we present GW signals with S/N > 8. The BBH-to-GRB ratio (ratio of numbers in the green to red histogram) is 7922/381 = 20.8. This is the expected number of BBH mergers GBM has to follow until the detection of the next gamma-ray counterpart. The notations are similar to those of Figure 2.</figcaption></figure>

### 3.5. Power-law Dependence on the Final Mass

For this scenario, we assume the EM luminosity, $L_γ$ has a power-law dependence on the final BH mass, $M_f$; $L_γ = L_*(M_f/M_*)^p$. We may use, e.g., $L_* = 1.8 \times 10^{49}$ erg s⁻¹ and $M_* = 62M_⊙$ 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_γ$ 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:

1. **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 an 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. **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. **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. 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 of this scenario is to gauge if currently observed BBHs could be

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Example of observed sources, where the EM energy is a fraction $\epsilon_{GW} = 3.4 \times 10^{-6}$ s⁻¹ of the emitted GW radiation ($L_γ = \epsilon_{GW}E_{GW}$, Section 3.4). In the above example the BBH-to-GRB ratio is 34.1. Notations are similar to Figure 2.</figcaption></figure>

expected to produce a detectable gamma-ray signal (for the BBH-to-GRB ratio in these cases, see Figures 2-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 an O3 sensitivity for the simulations that is superior to the previous two observing runs' sensitivity. The BBH-to-GRB ratio is expected to be lower for BBHs detected in a less sensitive configuration because the lower average distance of the sources implies 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, and thus is 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
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Effect of changing the power-law index on the BBH-to-GRB ratio. We assume the opening angles have a uniform distribution between 10° and 40°. The lines are averages of four random realizations of BBH orientations.</figcaption></figure>

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 (≤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 cases where the assumed mechanism reproduces GW150914-GBM exactly. The green histogram shows the gamma-ray flux for all the detectable BBHs with S/Ns of 12, or 8 for LIGO-Virgo. The initial cut for the jets that point elsewhere reduces 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-law 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 ~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_γ ∝ M_f^p$) the BBH-to-GRB ratio ranges from ~5 to a few times 10 (for $p = 2$ and S/N > 12) or to ≲200 (for $p = 2$ and S/N > 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 a 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 changes 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 have 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°) and wide, 90° case as well, to get a sense of how the number of GRBs will change.

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

Allowing for a lower S/N (>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 Offline Searches

It is important to note that the GBM sensitivity considered here applies to GRBs that were triggered in real time on board the spacecraft. GW150914-GBM was found in an offline search. Offline searches improve the sensitivity by a factor of
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Number of BBH observations with no GRB counterparts required to rule out a model that predicts a given BBH-to-GRB ratio. The different curves correspond to different confidence levels.</figcaption></figure>

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_{BGR}$ BBH observations for each GRB. The probability of observing 0 GRBs after $N_{BBH}$ BBH detections for such a model can be calculated from the binomial distribution:

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

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_{BBH} = N_{BGR}$). Using representative cases, e.g., $N_{BGR} = \{5, 10, 20, 100\}$, yields the following probabilities: $P(x = 0|N_{BBH} = N_{BGR}) = \{0.33, 0.35, 0.36, 0.37\}$. In other words, after $N_{BGR}$ non-detections there is still a ~35% probability that the model is correct and will lead to a BBH observation with a counterpart.

After $N_{BBH} > N_{BGR}$ observations with no GRB counterparts, what is the probability that the model predicting $N_{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_{BGR}$, with a desired confidence $P_{conf}$:

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

For example, a model predicting a BBH-to-GRB ratio of $N_{BGR} = 10$ can be ruled out with $P_{conf} = 10^{-2}$ confidence after ≈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_{BGR}$ with a given confidence.

### 4.6. Detection Prospects

Assuming a rate of 40 BBH detections during the year-long observing run O3 (Abbott et al. 2018; The LIGO Scientific Collaboration 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° and 40° and S/N > 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, and 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–20 BBHs for every gamma-ray counterpart detected. This is encouraging for the third observing run of LIGO-Virgo, which 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°$) can be ruled out with reasonable confidence. For narrowly beamed gamma-ray emission ($\theta_j = 20°$), 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 emphasize 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.

We thank Bing Zhang for valuable input, Francesco Pannarale for comments on the manuscript, and the anonymous referee for a valuable comment. P.V. acknowledges support from Fermi grant NNM11AA01A and 80NSSC17K0750. T.D. and E.B. 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. N.C. is supported by NSF grants PHY-1806990 and PHY-1505373.

This paper has LIGO document number LIGO-P1900147.

## ORCID iDs

P. Veres https://orcid.org/0000-0002-2149-9846

T. Dal Canton https://orcid.org/0000-0001-5078-9044

A. Goldstein https://orcid.org/0000-0002-0587-7042

T. B. Littenberg https://orcid.org/0000-0002-9574-578X
R. D. Preece https://orcid.org/0000-0003-1626-7335

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJL, 826, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, PhRvL, 116, 061102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, PhRvL, 116, 241102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJL, 848, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, PhRvL, 118, 221101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJL, 848, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, LRR, 21, 3

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 875, 161

Bartos, I., Kocsis, B., Haiman, Z., et al. 2017, ApJ, 835, 165

Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135

Bhat, N. P., Meegan, C. A., von Kienlin, A., et al. 2016, ApJS, 223, 28

Blackburn, L., Briggs, M. S., Camp, J., et al. 2015, ApJS, 217, 8

Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433

Buonanno, A., Kidder, L. E., & Lehner, L. 2008, PhRvD, 77, 026004

Burns, E. 2015, GCN, 18741, 1

Burns, E., Goldstein, A., Hui, C. M., et al. 2019, ApJ, 871, 90

Connaughton, V., Burns, E., Goldstein, A., et al. 2016, ApJL, 826, L6

Connaughton, V., Burns, E., Goldstein, A., et al. 2018, ApJL, 853, L9

Dai, L., McKinney, J. C., & Miller, M. C. 2017, MNRAS, 470, L92

Dalal, N., Holz, D. E., Hughes, S. A., & Jain, B. 2006, PhRvD, 74, 063006

Fedrow, J. M., Ott, C. D., Sperhake, U., et al. 2017, PhRvL, 119, 171103

Fong, W., Berger, E., Margutti, R., & Zauderer, B. A. 2015, ApJ, 815, 102

Goldstein, A. 2017, GCN, 21484, 1

Goldstein, A., Burns, E., Hamburg, R., et al. 2016, arXiv:1612.02395

Goldstein, A., Hamburg, R., Hui, C. M., et al. 2019, arXiv:1903.12597

Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJL, 848, L14

Gruber, D., Goldstein, A., Weller von Ahlefeld, V., et al. 2014, ApJS, 211, 12

Hamburg, R. 2017, GCN, 21222, 1

Harrison, F. A., Bloom, J. S., Frail, D. A., et al. 1999, ApJL, 523, L121

Husa, S., Khan, S., Hannam, M., et al. 2016, PhRvD, 93, 044006

Jiménez-Forteza, X., Keitel, D., Husa, S., et al. 2017, PhRvD, 95, 064024

Khan, S., Husa, S., Hannam, M., et al. 2016, PhRvD, 93, 44007

Kimura, S. S., Takahashi, S. Z., & Toma, K. 2017, MNRAS, 465, 4406

Kocevski, D., Burns, E., Goldstein, A., et al. 2018, ApJ, 862, 152

Komissarov, S. S., & Barkov, M. V. 2010, MNRAS, 402, L25

Kovetz, E. D., Cholis, I., Breysse, P. C., & Kamionkowski, M. 2017, PhRvD, 95, 103010

Li, X., Zhang, F.-W., Yuan, Q., et al. 2016, ApJL, 827, L16

Liu, T., Romero, G. E., Liu, M.-L., et al. 2016, ApJ, 826, 82

Loeb, A. 2016, ApJL, 819, L21

Lyutikov, M. 2016, arXiv:1602.07352

Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791

O'Leary, R. M., Meiron, Y., & Kocsis, B. 2016, ApJL, 824, L12

Perna, R., Lazzati, D., & Giacomazzo, B. 2016, ApJL, 821, L18

Racusin, J. L., Burns, E., Goldstein, A., et al. 2017, ApJ, 835, 82

Reynolds, C. S., Garofalo, D., & Begelman, M. C. 2006, ApJ, 651, 1023

Rhoads, J. E. 1999, ApJ, 525, 737

Ruffert, M., & Janka, H.-T. 1998, A&A, 338, 535

Schutz, B. F. 2011, CQGra, 28, 125023

The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, B. P., et al. 2018, arXiv:1811.12907

Veres, P., Preece, R. D., Goldstein, A., et al. 2016, ApJL, 827, L34

Woosley, S. E. 2016, ApJL, 824, L10

Zalamea, I., & Beloborodov, A. M. 2011, MNRAS, 410, 2302

Zhang, B. 2016, ApJL, 827, L31