Search for Gravitational Waves Associated with Gamma-Ray Bursts Detected by Fermi and Swift During the LIGO-Virgo Run O3b

We search for gravitational-wave signals associated with gamma-ray bursts detected by the Fermi and Swift satellites during the second half of the third observing run of Advanced LIGO and Advanced Virgo (1 November 2019 15:00 UTC-27 March 2020 17:00 UTC).We conduct two independent searches: a generic gravitational-wave transients search to analyze 86 gamma-ray bursts and an analysis to target binary mergers with at least one neutron star as short gamma-ray burst progenitors for 17 events. We find no significant evidence for gravitational-wave signals associated with any of these gamma-ray bursts. A weighted binomial test of the combined results finds no evidence for sub-threshold gravitational wave signals associated with this GRB ensemble either. We use several source types and signal morphologies during the searches, resulting in lower bounds on the estimated distance to each gamma-ray burst. Finally, we constrain the population of low luminosity short gamma-ray bursts using results from the first to the third observing runs of Advanced LIGO and Advanced Virgo. The resulting population is in accordance with the local binary neutron star merger rate.


INTRODUCTION
Gamma-ray bursts (GRBs; Kumar & Zhang 2015) are intense and highly variable flashes of gamma-rays (the prompt emission), followed by a long-lasting, multiwavelength emission (the afterglow emission), typically * Deceased, August 2020. observed in X-rays, optical, radio, and sometimes in gamma-rays. They are believed to be powered by ultrarelativistic jets produced by rapid accretion onto a central compact object: a black hole (BH; Woosley 1993;Popham et al. 1999) or a magnetar (Dai & Lu 1998;Zhang & Mészáros 2001).
GRBs are divided into two classes, depending on the duration and the spectral hardness of the prompt emis-sion (Kouveliotou et al. 1993): long, soft GRBs (duration 2 s) and short, hard GRBs (duration < 2 s).
Long GRBs are thought to be associated with the core collapse of massive stars. This connection is observationally supported by the identification of supernova (SN) signatures in a number of sufficiently close long GRBs (Galama et al. 1998;Hjorth et al. 2003;Stanek et al. 2003). Core-collapsing massive stars are also expected to emit gravitational waves (GWs) if there is some asymmetry in the stellar-envelope ejection phase (Kotake et al. 2006;Ott 2009;Gossan et al. 2016). Stateof-the-art models predict that such GW radiation can be detected by current generation GW interferometers only within our Galaxy ; however, according to more extreme phenomenological models, such as long-lived bar-mode instabilities and disk fragmentation instabilities, GW radiation could be detected even for extra-galactic sources (Fryer et al. 2002;van Putten et al. 2004;Piro & Pfahl 2007;Corsi & Mészáros 2009;Gossan et al. 2016;).
Short GRBs were long believed to be associated with compact binary coalescence (CBC) composed by two neutron stars (NSs), a binary neutron star (BNS) system, or a NS and a BH, a NSBH binary (Eichler et al. 1989;Paczynski 1991;Narayan et al. 1992). The definitive proof of this association (Abbott et al. 2017a,d) came with the joint detection of the BNS merger GW signal GW170817 (Abbott et al. 2017c(Abbott et al. , 2019c and the GRB 170817A (Savchenko et al. 2017;Goldstein et al. 2017). The ground-breaking electromagnetic follow-up campaign performed after this joint detection allowed the identification of the associated kilonova emission and of the GRB afterglow emission (see Abbott et al. 2017d and references therein).
GRB 170817A was 2 to 6 orders of magnitude less energetic than other GRBs (Abbott et al. 2017a); the low luminosity of this source, together with the evolution of the X-ray and radio light curve (Troja et al. 2018;Margutti et al. 2018;D'Avanzo et al. 2018) suggested an off-axis GRB with a relativistic structured jet or a cocoon emission from the relativistic jet shocking its surrounding non-relativistic material. Subsequent very long baseline interferometry observations allowed constraints on the source size and its displacement, indicating that GW170817 produced a structured relativistic jet (Ghirlanda et al. 2019;Mooley et al. 2018).
In  we presented targeted GW follow-up of GRBs reported during the first part of the third observing run of Advanced LIGO and Advanced Virgo (O3a; 1 April 2019 15:00 UTC-1 October 2019 15:00 UTC) by Fermi 's Gamma-Ray Burst Monitor (Fermi /GBM; Meegan et al. 2009) and Swift's Burst Alert Telescope (Swift/BAT; Gehrels et al. 2004; Barthelmy et al. 2005;Tohuvavohu et al. 2020). No significant evidence for GW signals associated with the GRBs that have been followed up has been found, nor for a population of unidentified sub-threshold signals.
In this paper we present targeted GW follow-up of GRBs reported during the second part of the third observing run of Advanced LIGO and Advanced Virgo (O3b) by Fermi /GBM and Swift/BAT. O3b took place between 1 November 2019 15:00 UTC and 27 March 2020 17:00 UTC. During O3b, 35 CBC events have been identified with an inferred probability of astrophysical CBC origin of p astro > 0.5 (Abbott et al. 2021a). The majority of them are classified as mergers of binary black hole (BBH) systems; however, several events are consistent with binary systems with at least one NS (Abbott et al. 2021a). One other event with lower p astro was also published as a possible NSBH coalescence (Abbott et al. 2021b). No EM counterparts have been reported so far in association with these events; however, given their large distances ( 300 Mpc) and their large error in the sky localization (Abbott et al. 2021a), it would have been difficult to detect an EM signal in association with these GW events.
In Section 2 we discuss the sample of GRBs analyzed in this paper. In Section 3 we summarize the methods used to follow-up GRBs. In Section 4 we describe the results, and in Section 5 we present a population model analysis. Finally, in Section 6 we present our concluding remarks.

GRBS DURING O3B
Our GRB sample consists of 108 events that occurred between 1 November 2019 15:00 UTC and 27 March 2020 17:00 UTC. The vast majority of these events were identified in low-latency via notices circulated by the Gamma-ray Coordinates Network (GCN) and subsequently refined with additional data from the Swift/BAT catalog and the Fermi /GBM catalog. 1 The Vetting Automation and Literature Informed Database (VALID; Coyne 2015) is a dedicated processing system that tracks updates to the observed GRB parameters, comparing time and localization data to ensure that the latest results are used for our GW analyses, and employing an automated literature search to identify particularly noteworthy events.
We identify candidate events by classifying each GRB as long, short, or ambiguous. We classify events based on their T 90 (and its associated error δT 90 ), which is the time interval over which 90% of the total backgroundsubtracted photon counts are observed. GRBs are classified as short when T 90 +|δT 90 | < 2 s, GRBs are classified as long when T 90 − |δT 90 | > 4 s, and all remaining GRBs are labeled as ambiguous. This long/short classification based on duration is only a general trend, and is not a perfect discriminator. For more robust classification one must also consider spectral properties, most commonly the spectral hardness or peak energy of the event, but since our sample consists of observations from multiple observatories with different spectral sensitivities we do not employ such quantities when organizing our sample.
This classification process results in 7 short GRBs, 12 ambiguous GRBs, and 89 long GRBs. Of all these GRBs, only 2 have known redshifts: In keeping with previous studies of this kind (Abbott et al. 2017e, 2019d(Abbott et al. 2017e, , 2020, we apply a generic transient search to all events, regardless of classification. In order to maximize our chances at identifying potential CBC candidates, we apply our modeled search to all short and ambiguous GRBs. We also follow the same requirements on amount of data available within our network to process a given GRB. For the modeled search we select GRBs if there is a minimum amount of time in at least one detector around the time of the event. This gives us 17 events for our analysis corresponding with the observing time for the same selection criteria (96.6 % with at least one interferometer in observing mode). For the generic transient search, we perform the selection by requiring enough data in at least two interferometers. This leads to 86 GRBs to analyze and is also compatible with the network observing time of at least 2 detectors (85.3 %).

Modeled search for compact binary mergers
This analysis is carried out by a coherent matched filtering pipeline, PyGRB (Harry & Fairhurst 2011;Williamson et al. 2014), contained within the opensource PyCBC (Nitz et al. 2020) suite which also relies heavily on the LALSuite (LIGO Scientific Collaboration 2018) library. These searches seek to find candidate GW signals coincident with the GRB triggers due to the inspiral and merger of BNS or NSBH binaries. We define a window around each GRB trigger, the on-source window, which is [−5, +1] s from the GRB trigger time.
This window is based on the assumption that a GW may precede the prompt GRB emission by several seconds (Lee & Ramirez-Ruiz 2007;, and was demonstrated by GW170817 (Abbott et al. 2017b). The search also uses time surrounding the trigger, split into 6 s off-source windows, to estimate the background. In total, the search uses ∼90 min of data around each GRB trigger to assign a significance to candidate events by ranking them against the background.
The analysis requires a bank of template waveforms to carry out the matched filtering. We generated this bank using both geometric (Brown et al. 2012;Harry et al. 2014) and stochastic methods (Harry et al. 2008) for BNS and NSBH signals. The waveforms used in generating this bank are phenomenological inspiral-mergerringdown waveform models of the IMRPhenomD family Khan et al. 2016). We choose to place limits on the bank, identical to those used in the O3a template bank , such that any NS masses are limited to [1.0, 2.8]M and BH masses are within [2.8, 25]M . We conservatively set the mass cutoff between NS and BH based on an NS equation of state (Kalogera & Baym 1996). Functionally, this cutoff has no effect on the waveforms and is just used for nomenclature. The bank only contains aligned-spin BNS and NSBH binaries where the maximum dimensionless spin magnitude for NSs is 0.05 from the largest observed NS spin in a binary (Burgay et al. 2003). For BH, we limit the spin to 0.998 based on theory (Thorne 1974). Finally, we check to ensure that all potential binaries are viable GRB progenitors with the creation of an accretion disk able to power a GRB (Pannarale & Ohme 2014).
The only structural change between this bank and the bank used in the O3a modeled searches ) is the template placement for NSBH systems with total mass M < 6M . Both banks are constructed by first performing a geometric generation for a part of the parameter space. These templates are then seeded to a stochastic generation that fills the rest of the parameter space (Capano et al. 2016). The difference between the O3a and O3b banks is that the geometric generation for the O3a bank extended through the low-mass NSBH region whereas the O3b bank limits the geometric generation to the BNS region. We made this change based on a bank verification which tests a bank's ability to recover a set of signals. The result of this verification is a fitting-factor (FF) that quantitatively measures the bank's performance (Apostolatos 1995). The target for our template banks is to minimize the number of signals that have a FF less than a threshold, which we set at 0.97 for our offline searches. For the same set of signals in the low-mass NSBH region, the bank with a limited geometric generation recovers a factor of ten less signals with a fitting factor below 0.97-when compared to the extended geometric bank. These results show that the limited geometric approach creates a more sensitive template bank for our searches.
PyGRB uses this bank to rank candidate signals based on a re-weighted optimal SNR. This optimal SNR is the result of the coherent matched filter, and is re-weighted by how well the template matches the identified signal (Harry & Fairhurst 2011;Williamson et al. 2014). The search can then rank the significance of any event against the background using the off-source windows. In order to improve this ranking statistic, we artificially increase the amount of off-source data by performing time slides (Williamson et al. 2014).
To further determine the sensitivity of our searches, we inject signals into the off-source data and attempt to recover them. The signals that we choose to inject are generally in the same BNS and NSBH domains as the template bank, with a few important distinctions. Again, we replicate what was done in O3a , where the injected signals are split into three sets; a BNS set with non-aligned (precessing) spins, an aligned-spin NSBH set, and a precessing NSBH set. The NS masses in a BNS binary are selected randomly from a normal distribution with a mean of 1.4M and variance of 0.2M (Özel et al. 2012). For NSBH binaries, NS masses are selected from a normal distribution with slightly more variance (µ = 1.4M , σ = 0.4M ). The larger width reflects the greater uncertainty arising from a lack of observed NSBH systems. BH masses are randomly selected from the following normal distribution (µ = 10.0M , σ = 6.0M ). For all cases we place limits on the distributions similar to those used for the template bank. Randomly selected spin magnitudes are less than 0.4 for NSs based on the maximum observed pulsar spin (Hessels et al. 2006), and less than 0.98 for BHs (Miller & Miller 2014). For the two sets of injections that allow precessing signals, the orientations are also randomly selected. We also choose to use different waveform families than the ones used to generate the template bank to account for modeling uncertainty. We generate the BNS injections using the SpinTaylorT2 family, which are post-Newtonian approximations in the time domain (Sathyaprakash & Dhurandhar 1991;Blanchet et al. 1996;Bohé et al. 2013;Arun et al. 2009;Mikoczi et al. 2005;Bohé et al. 2015;Mishra et al. 2016). The NSBH sets both make use of the SEOBNRv3 family of waveforms. These waveforms are effective-one-body approximates that are tuned for precessing systems Taracchini et al. 2014;Babak et al. 2017).
As with the template bank, we check to ensure that generated systems are capable GRB progenitors (Pannarale & Ohme 2014). These injection sets allow us to calculate the 90% exclusion distance (D 90 ), which is the distance at which we recover 90% of the injected signals with a significant ranking statistic.

Search for generic GW transients
This analysis, carried out with the X-Pipeline software package (Sutton et al. 2010;Was et al. 2012), searches for excess power that is coherent across the GW detector network and consistent with the sky localization and time window of each GRB. Like the previous X-Pipeline analyses (Abbott et al. 2017e, 2019d(Abbott et al. 2017e, , 2020, the search time window starts 600 s before the GRB trigger time and ends at 60 s after trigger time, or T 90 after if T 90 > 60 s. This is sufficient to cover the time delay between GW emission from a progenitor and any GRB prompt emission (Koshut et al. 1995;Aloy et al. 2000;MacFadyen et al. 2001;Zhang et al. 2003;Lazzati 2005;Wang & Mészáros 2007;Burlon et al. 2008Burlon et al. , 2009Lazzati et al. 2009;). While some GW emissions, such as from core-collapse SNe, are expected to reach frequencies up to a few kilohertz (Radice et al. 2019), we restrict our search frequency range to the most sensitive band of the GW detectors, 20-500 Hz, since detecting such signals above a few hundred hertz requires extremely high GW energies (Abbott et al. 2019b, Fig. 4) and expanding the frequency range would also significantly increase the computational cost.
X-Pipeline produces time-frequency maps of the GW data coherently combined between the detectors. These maps give access to the temporal evolution of the spectral properties of the signal and enable the pipeline to search for clusters of pixels containing excess energy, referred to as events. The pipeline assigns each event a detection statistic based on energy and ranks them accordingly. A coherent consistency test, based on correlations between data in different detectors, then vetoes events that are associated with noise transients. The surviving event with the largest ranking statistic is the best candidate for a GW detection, and the search quantifies its significance as the probability of the event being produced by the background alone. This is determined by comparing the SNR of the trigger within the 660 s on-source window to the distribution of the SNRs of the loudest triggers in the 660 s off-source windows. As a requirement, the off-source data consist of at least ∼ 1.5 hours of coincident data from at least two detectors around the time of a GRB. This is small enough to select data where the detectors should be in a similar state of operation as during the GRB on-source window, and large enough so that probability estimates using artificial time-shifting of the data are at the sub-percent level.
We quantify the sensitivity of the generic transient search by injecting simulated signals into off-source data. For each waveform family injected we determine the largest significance of any surviving cluster associated with the injections. We compute the percentage of injections that have a significance higher than the best event candidate and look for the amplitude at which this percentage is above 90%, which sets the upper limit. We include O3b calibration errors (Sun et al. 2021;Acernese et al. 2021) by jittering the amplitude and arrival time according to a Gaussian distribution representative of the calibration uncertainties. As with the modeled search, these injection sets allow us to calculate 90% exclusion distances.
We choose simulated waveforms to cover the search parameter space of three distinct sets of circular waveforms: BNS and NSBH binary inspiral signals, stellar collapse, and disk instability models.
• Circular sine-Gaussian (CSG): signals representing GW emission from stellar collapses defined in Eq.
(1) of Abbott et al. (2017e) with a Q factor of 9 and varying center frequency of 70 Hz, 100 Hz, 150 Hz, and 300 Hz. In all cases, we assume an optimistic emission of energy in GWs of E GW = 10 −2 M c 2 .
• Binary inspiral: signals are characterized by a Gaussian distribution centered at 1.4M , with a width of 0.2M for NS in a BNS, and with a width of 0.4M for NS in NSBH. The distribution for GWs emitted by BNS mergers addresses the case of short GRB events as in Abbott et al. (2017e) and adopted in PyGRB search (Sec. 3.1).
• Accretion disk instability (ADI): long-duration waveforms for GWs produced by instabilities in the magnetically suspended torus around a rapidly spinning BH. The model specifics and parameters used to generate the five families of ADI signals are the same as in the previous searches (Abbott et al. 2017e, 2019d(Abbott et al. 2017e, , 2020. In the O3a search, the sensitivity to long-duration ( 10 s) signals was often limited by loud background noise transients known as glitches (Davis et al. 2021). While X-Pipeline's coherent consistency tests easily veto these glitches, many long-duration simulated signals would overlap such a glitch by chance. In these cases the simulated signal and glitch would be clustered together and subsequently vetoed together. To address this problem, we implemented an autogating procedure for O3b. For each detector, we compute the total energy in the whitened data stream over a 1 s window. If this total fluctuates by more than 50 standard deviations above the median value, then the data is zeroed out over the interval where the threshold is exceeded and we apply an inverse 1 s Tukey window at each end of the zeroed interval to transition smoothly between the whitened and zeroed data. To minimise the possibility of a loud GW transient triggering a gate, the procedure cancels a gate if there is a simultaneous energy excursion above 10 standard deviations in any other detector. The threshold of 50 standard deviations is low enough to gate the most problematic loud glitches, while being high enough that the only GWs zeroed out by the gate would have been detectable by all-sky searches. Empirically we find that this procedure is effective at reducing the impact of loud glitches without affecting the sensitivity to low-amplitude GW signals.
For both search methods, we rank each candidate by calculating a p-value, the probability of an event or a louder one in the on-source data, given the background distribution, under the null hypothesis. The p-value is calculated by counting the fraction of background trials that contain an event with a greater signal-to-noise ratio than that of the loudest on-source event.

RESULTS OF ANALYSES
We followed up 86 GRB triggers with the generic transient method and 17 GRBs (those categorized as short or ambiguous) with the modeled search. None of the analyses indicate the presence of a statistically significant GW signal associated with one or more of the GRBs. This null result is consistent with the estimated GW-GRB joint detection rate with Fermi /GBM of 0.07-1.80 per year reported previously in Abbott et al. (2019a) for the second observing run of Advanced LIGO and Advanced Virgo (O2).
We present the cumulative p-value distributions from both search methods in Figures 1 and 2. In these plots, a significant event would appear at a much lower p-value in the lower left corner of the plots, and be outside (to the left) of the 90% confidence region. Both plots show that the p-value distributions are consistent with the background.
The most significant event from the modeled search had a p-value of 1.08 × 10 −2 (GRB 200129A). Through further investigation of this candidate event, a period of excess noise in one of the detectors was discovered ∼ 20 s before the candidate time. To determine the effect of this noise on the candidate, we used BayesWave to reconstruct the glitch and then clean the data by sub-  Figure 1. The cumulative distribution of p-values for the loudest on-source events for the modeled search in O3b. If a trigger is found in the on-source the upper and lower limits are identical to the reported p-value. If no trigger is identified in the on-source window, we set an upper limit on the p-value of 1, and a lower limit equal to the fraction of off-source trials that also did not contain a trigger. The upper limits are plotted as the curve with full circles and the lower limits are plotted as the curve with empty circles. The dashed line indicates an expected uniform distribution of p-values under a no-signal hypothesis, with the corresponding 90% band as the dotted lines.
tracting the reconstruction (Cornish et al. 2021;Pankow et al. 2018). After this cleaning, we conducted a coherent matched-filtering on the cleaned data and the recovered candidate was no longer significant with respect to the background. This result suggests that much of the power of the candidate was caused by noise and not a GW. Even if there is a quiet GW at this time, it is not strong enough without the contribution from the glitch to survive ranking against the background in the analysis.
The lowest reported p-value found during O3b for the generic transient search was 7.95×10 −3 (GRB 200224B). Although this p-value is very small, it is not unexpected given the high number of GRBs analyzed.
Given that no loud GW signals were observed coincident with any of the GRBs in either of our searches, we perform a weighted binomial test to determine the probability of observing our set of p-values assuming a uniform background distribution. A small probability would suggest that there may be a population of subthreshold GW signals that our searches did not identify. This type of weighted binomial test, fully described in the Appendix of Abadie et al. (2012), uses the lowest reweighted p-values from the searches. The resulting probability for the modeled search is 0.07. If we remove GRB 200129A, for which the small p-value is the result of noise, the probability becomes 0.68, suggesting no population of weak GW signals. For the generic transient search, the test gives a probability of 0.76. These same weighted binomial tests carried out in O3a returned probabilities of 0.43 and 0.30 for the modeled and generic transient searches, respectively . In O2 (removing GW170817/GRB 170817A) and the first observing run of Advanced LIGO and Advanced Virgo (O1) the probabilities were 0.30 and 0.75, and 0.57 and 0.75, respectively (Abbott et al. 2019a. As in these previous analyses, the probabilities obtained in O3b suggest that no weak GWs can be attributed to the population of GRBs. In Fig. 3, we present the cumulative 90% exclusion distances for the 17 GRBs analyzed with the modeled search. The first of these 17 GRBs, GRB 200323A, has significantly lower exclusion distances than the rest. We can attribute this to the fact that the analysis of this GRB only used data from the Virgo interferometer. Furthermore, this GRB has a sub-optimal sky location for  the Virgo interferometer with a sensitivity, when compared to an optimal sky-location, of ∼ 30%. Both of these factors produce the relatively small exclusion distances for the first step in the histogram. Table 1 reports the median D 90 for the 17 GRBs analyzed with the mod-  . Cumulative histograms of the 90% confidence exclusion distances, D90, for accretion disk instability (ADI) signal model A (orange, thin line) and circular sine-Gaussian (CSG) 150 Hz model (green, thick line). For a given GRB and signal model, this is the distance within which 90% of simulated signals inserted into off-source data are successfully recovered with a significance greater than the loudest on-source trigger. eled search. It shows median values for all three of the injected signal types described in Sec. 3.1. For comparison, all three of these median values are 10-30% larger than those reported from the same modeled search in O3a . This difference stems from having a larger fraction of GRBs in O3b that by chance arrived with better LIGO-Virgo antenna factors on average, bringing up the median values. The individual D 90 values for each of the 17 GRBs analyzed with the modeled search can be seen in Table 2.
Similar to the modeled search, we derive a 90% confidence level lower limit on the distance for each of the 86 GRBs analyzed with the generic transient search, based on the different emission models described in Sec  Fig. 4. The limits reported depend on the sensitivity of the instruments in the network, which change with time and sky localization of the GRB events. We marginalize these limits over errors introduced by detector calibration. In Table 1, we report the median exclusion distance limits, D 90 , for the set of GRBs for the different signals described in Sec. 3.2. The limits vary by nearly an order of magnitude due to the variety of signals used in our analysis. On average the median values for the O3b generic transient search are about 50% greater than those reported in O3a . We can primarily attribute this improvement to the use of autogating in O3b: the increase in exclusion distances is highest (up to a factor of two) for the longest-duration waveforms, which are most impacted by the glitches removed by autogating (as explained in Sec. 3). The exclusion distances for the shorter-duration CSG waveforms, which are not expected to be affected by autogating, increased by about 30% on average. This is more than could be accounted for by chance differences in the LIGO-Virgo antenna factors between the two samples. Rather, the increase is likely due to improvements in the performance of the detectors themselves, such as through the reduction of noise caused by scattered light in the LIGO detectors (Soni et al. 2021) or the improvement in sensitivity of the Virgo detector (Davis et al. 2021). We report the D 90 values found for each GRB in the case of ADI model A simulated signals and CSG simulated signals with central frequency of 150 Hz in Table 2, at the end of this paper.

POPULATION STUDIES
We use the results obtained from the GW followup analysis of GRBs to put constraints on the lowluminosity short GRB population. For this purpose, we describe the short GRB population through a simple luminosity function model following (Wanderman & Piran 2015), extended at low luminosities following the procedure described in (Abbott et al. 2019d). We can then model the luminosity distribution through a power law with two breaks where L iso is the isotropic equivalent GRB luminosity and for which we have L * = 2 × 10 52 erg s −1 , L * * = 5 × 10 49 erg s −1 , α L = 0.94 and β L = 2 (Wanderman & Piran 2015). We do not take into account the measurement uncertainties for those fixed parameters as they would not significantly influence the analysis. The parameters on which we aim to put constraints us-ing the joint GW-GRB analysis are the low-luminosity power index γ L and the low-luminosity cutoff for our population L 0 . To make the dependence from these parameters clearer, we refer to the luminosity distribution as φ 0 (L iso ) ≡ φ 0 (L iso , γ L , L 0 ). A Bayesian analysis constrains the parameters γ L and L 0 using the results from the O1, O2, O3a and O3b PyGRB searches (Harry & Fairhurst 2011;Williamson et al. 2014;Abbott et al. 2019d and the results on BNS rates from Abbott et al. (2021c). Under certain conditions, NSBH mergers can also produce sGRBs (Narayan et al. 1992) and a small fraction of sGRBs can arise from local magnetar giant flares (Burns et al. 2021). For simplicity, we ignore those relatively uncommon possibilities here. We assume that BNS coalescences are the only progenitors for short GRBs, since there are restricted conditions under which an NSBH coalescence results into a short GRB (Pannarale & Ohme 2014).
First, we compute the observed cumulative rate distribution C obs R (z, γ L , L 0 ) as a function of redshift z, γ L and L 0 . To do so, we take into account the cosmic rate density for short GRB explosions ψ(z) adopting its form given in Wanderman & Piran (2015). A Band function models the energy spectrum of the short GRBs (Band et al. 1993) with power indices α Band = −0.5, β Band = −2.25 and peak energy E peak = 800 keV, and we use Eq. (1) as the luminosity distribution function for our population of short GRBs. As in Wanderman & Piran (2015), we consider short GRBs detectable in gamma-rays when their 64 ms peak photon flux is above P th 64 = 2.37 photons cm −2 s −1 in the energy window considered for Fermi /GBM, i.e. [50-300] keV. We then compute the cumulative observed rate distribution as where the differential probability of having an observed short GRB is defined as Here in Eq. (3), ψ(z) is the short GRB redshift distribution, dV /dz is the differential comoving volume and (z, γ L , L 0 ) is the efficiency curve for the Fermi /GBM detector as a function of redshift and of the lowluminosity parameters of the luminosity distribution.
Using this rate distributon, we build a prior probability distribution function (PDF) Π(γ L , L 0 ). The prior is built starting from a flat distribution in the logarithms of the local observed rate density and of L 0 , since those quantities can span over several orders of magnitude, then it is rescaled by the posterior cumulative distribution function of the BNS local rate density from Abbott et al. (2021c). This last factor formalizes the assumption that most of the short GRBs are produced in BNS coalescences. For all computations, we consider a flat ΛCDM cosmology with h 0 = 0.7, Ω m = 0.3 and Ω Λ = 0.7, in order to be consistent with the analysis done in Wanderman & Piran (2015).
We define the likelihood function L(x|γ L , L 0 ) (where x indicates our set of data) as the probability of detecting no GW transients associated with short or ambiguous GRBs during O1, O3a and O3b and of detecting one single GW transient associated to a GRB observed during the O2 run. Furthermore, we impose that the joint detection occured at the redshift measured for NGC 4993, the host galaxy of the event GW170817 (z NGC 4993 = 0.009783; Levan et al. 2017) and that the luminosity of the corresponding GRB is in the luminosity range measured for GRB 170817A L GRB 170817A = (1.6 ± 0.6) × 10 47 erg s −1 (Abbott et al. 2017a). For our purpose we use the set of GW efficiency curves computed through the PyGRB analysis of the short and ambiguous GRBs events detected during the O1, O2, O3a and O3b runs (respectively 20, 41, 32 and 17 events analyzed). 2 Given a detected GRB i during O2, we compute the probability of a joint GW detection like the one observed during this run Here η i (z) is the efficiency curve corresponding to the given GRB and dP GRB obs /dz has been defined in Eq. (3). In order to set the joint detection to have the same luminosity of GRB 170817A and the same redshift of GW170817, we choose NL(L) to be a log-normal distribution with mean L GRB170817A =L with σL being the error on the measurement ofL, and we use a Dirac delta distribution δ(z − z NGC 4993 ) because our analysis is insensitive to small variations in the assumed redshift.
Analogously, we can compute the probability of not having a joint GW detection associated to a given GRB 2 There are actually 42 efficiency curves available from the O2 PyGRB analysis, but the efficiency curve corresponding to GRB 170817A was not computed properly since the pipeline considered the GW170817 event as a background event.
detected during O1, O3a or O3b We then obtain that the probability of a single joint detection during O2 is while the probability of not having a joint detection during O1, O3a and O3b is then the obtained likelihood is Finally, we compute the posterior P (γ L , L 0 |x) ∝ L(x|γ L , L 0 )Π(γ L , L 0 ), the contour plot for which is shown on Fig. 5, with contours in blue and red corresponding respectively to the posterior 90% and 50% credible regions. The constant rate curves shape the posterior: if we fix a value of the rate, higher values for the low-luminosity cutoff L 0 favor higher values of the low-luminosity power index γ L . Each credible region's L 0 value is compatible with the luminosity value range of GRB 170817A. Finally, the 90% credible region curve does not close for low values of L 0 : this is due to the fact that we do not have any information about events down to those luminosities and for this reason we did not explore lower values for L 0 . By marginalizing the posterior PDF over L 0 , we obtain that γ L = 0.28 ± 0.45.
To present these results in the luminosity function space, we compute the rate curves dR 0 /d log L for pairs of values (γ L , L 0 ) sampled according to the posterior distribution P (γ L , L 0 |x). From this set of curves we obtain the median and credible intervals on the luminosity distribution.
The plot in the top panel of Fig. 6 shows dR 0 /d log L 50% credible intervals as functions of log L and compares them to other estimations performed in other works (Salafia et al. 2020;Tan & Yu 2020;Ghirlanda et al. 2016). It illustrates how the short GRB luminosity functions in our model peaks around L ∼ L GRB 170817A , considering this the only short GRB event observed at such a low luminosity. . The contours correspond to the 90% and 50% credible regions (respectively in blue and red) for the two parameters. The bounds regions for those two parameters are compatible with the measured luminosity from GRB 170817A (yellow dashed line with shaded area) as its value is greater than L0 for the bulk of the values of our population. The marginalized posterior for L0 peaks around L = LGRB 170817A because of the likelihood factor which requires that the joint detection happened around that value.
The plot in the bottom panel of Fig. 6 shows the inverse cumulative short GRB rate density distribution R 0 (> L) as a function of the luminosity L. The credible intervals corresponding to the sampled curve are compatible with the BNS rate density measured for Abbott et al. (2021c).
Given the present results on the low-luminosity short GRB population and the expected sensitivity for the fourth observing run of Advanced LIGO and Advanced Virgo (O4; , and only considering short GRBs detected by Fermi/GBM as onboard triggers, we estimate a joint GW-GRB detection rate of R O4 GW−GRB = 1.04 +0.26 −0.27 yr −1 during the next data collecting period.  For GRBs flagged as either short or ambiguous (see Sec. 2), we ran a template-based search for BNS and NSBH waveforms (Harry & Fairhurst 2011;Williamson et al. 2014). We also ran on all GRBs a generic transient analysis to look for GW signals (Sutton et al. 2010;Was et al. 2012). We did not find any significant GW candidate in coincidence with the GRBs we analyzed. Our results are consistent with the previously predicted detection rate of 0.07-1.8 per year for O3 (Abbott et al. 2019a). We also performed a weighted binomial test to search for a population of subthreshold GW signals in our sample. We did not find strong evidence for any such event. We used different emission models to put a lower bound on the distances of the GRB progenitors. The 90% exclusion distances are reported in Table 2 for all the GRBs in our sample, along with timing and localization information as well as information on detectors used in the analyses. Finally, we performed a population study for all GRBs analyzed with the modeled search in O1 (Abbott et al. 2017e), O2 (Abbott et al. 2019a), O3a ) and O3b. Starting from a broken power law to model our population and constraining two of its parameters through Bayesian inference, we found that our luminosity function peaks around the luminosity value measured for GRB 170817A with this model. Furthermore, the local rate density for short GRBs is compatible with that of BNS events. Based on the present population study, we provided an estimate of the joint GW-GRB detection rate for the O4 run.
This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands Organization for Scientific Research (NWO), for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific We would like to thank all of the essential workers who put their health at risk during the COVID-19 pandemic, without whom we would not have been able to complete this work. Table 2. GRB details and associated GW emission limits for each of the Fermi and Swift GRBs followed up on during O3b. The GRB Name column reports each GRB's formal designation (Barthelmy et al. 2009) Table 2 continued   Table 2 continued