A Joint Fermi-GBM and Swift-BAT Analysis of Gravitational-wave Candidates from the Third Gravitational-wave Observing Run

We present Fermi Gamma-ray Burst Monitor ( Fermi-GBM ) and Swift Burst Alert Telescope ( Swift-BAT ) searches for gamma-ray / X-ray counterparts to gravitational-wave ( GW ) candidate events identi ﬁ ed during the third observing run of the Advanced LIGO and Advanced Virgo detectors. Using Fermi-GBM onboard triggers and subthreshold gamma-ray burst ( GRB ) candidates found in the Fermi-GBM ground analyses, the Targeted Search and the Untargeted Search, we investigate whether there are any coincident GRBs associated with the GWs. We also search the Swift-BAT rate data around the GW times to determine whether a GRB counterpart is present. No counterparts are found. Using both the Fermi-GBM Targeted Search and the Swift-BAT search, we calculate ﬂ ux upper limits and present joint upper limits on the gamma-ray luminosity of each GW. Given these limits, we constrain theoretical models for the emission of gamma rays from binary black hole mergers.


INTRODUCTION
The detection of GW170817 (Abbott et al. 2017) coincident with the short gamma-ray burst GRB 170817A (Goldstein et al. 2017;Savchenko et al. 2017) was a ground-breaking discovery for the multimessenger era.Not only was it the first binary neutron star (BNS) merger detected by the gravitational-wave (GW) instruments Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2014), it was also the first, and to date only, GW detection with a confirmed electromagnetic (EM) counterpart.Since then, the search for EM emission from more of these extreme events has been at the forefront of multimessenger astronomy, particularly in the gamma-ray energy band since GRB 170817A demonstrated that BNS mergers are a progenitor of short gamma-ray bursts (GRBs) (Abbott et al. 2017).GWs have also been observed from the mergers of other compact objects, such as binary black hole (BBH) and neutron star-black hole (NSBH) systems (Abbott et al. 2019;Abbott et al. 2021;Abbott et al. 2021;Abbott et al. 2021); however, no additional EM counterparts have been confirmed as they have been inconclusive (Connaughton et al. 2016;LSC and Virgo and Fermi-GBM Team 2019a,b) or are still under debate (Graham et al. 2020;Ashton et al. 2021;Bustillo et al. 2021;De Paolis et al. 2020;Palmese et al. 2021).
GRB 170817A was first reported by the Fermi Gamma-ray Burst Monitor (GBM; Meegan et al. 2009), a space-based gamma-ray instrument sensitive from 8 keV to 40 MeV.This wide energy range of Fermi -GBM combined with its large field-of-view (FoV) and rapid alert abilities make it an ideal platform to search * Deceased, August 2020.
for gamma-ray counterparts to GWs in real time.
Fermi -GBM also provides continuous time tagged event (CTTE) data with a 6-hour latency that enables sensitive searches for short GRBs on the ground.Two of these searches are the Untargeted Search, a blind search of Fermi -GBM data for short GRBs, and the Targeted Search, which uses an external time to search for a short GRB (Blackburn et al. 2015;Goldstein et al. 2019).
Both were previously used to look for sub-threshold GRBs coincident with GWs from the first two LIGO-Virgo observing runs.Additionally, the Burst Alert Telescope (BAT; Barthelmy et al. 2005) on-board the Neil Gehrels Swift Observatory (hereafter referred to as Swift) provides excellent sensitivity to detecting hard X-ray and gammaray transients (Gehrels et al. 2004).Swift-BAT primarily runs in a survey mode that continuously evaluates photon rate increases and potential GRB triggers.An increase in the observed photon rate can trigger the on-board image-processing algorithms which can yield ∼arcminute GRB localizations.
Ideally, Swift-BAT would detect and localize a GRB produced by a binary merger independently of the GW detection.If a GRB does not trigger an on-board detection, continuous count rate lightcurves are still available for offline ground searches.Although Swift-BAT has been used to search for public and sub-threshold GWs during the LIGO-Virgo observing runs, this work presents the first systematic search of Swift-BAT data from a LIGO-Virgo observing run.
The first observing run (O1) operated from September 2015 to January 2016, producing the first detection of GWs from a BBH merger (GW150914; Abbott et al. 2016).Burns et al. (2019) used the Fermi -GBM Targeted Search to identify both triggered and sub-threshold GRB candidates in coincidence with GW candidates from O1.The most significant gammaray candidate found by the search was within 0.4 s of GW150914; however, it could not be confirmed as a counterpart due to its weak signal and poor localization (Connaughton et al. 2016;Greiner et al. 2016;Connaughton et al. 2018).
The second observing run (O2) took place from November 2016 to August 2017, resulting in the detections of GW170817, GRB 170817A, and the kilonova AT2017gfo (Chornock et al. 2017;Cowperthwaite et al. 2017;Nicholl et al. 2017;Soares-Santos et al. 2017;Tanvir et al. 2017;Margutti & Chornock 2021).Following O2, the LIGO Scientific and Virgo Collaboration published its first catalog of GW signals called the Gravitational-Wave Transient Catalog 1 (GWTC-1; Abbott et al. 2019) using a re-analysis of data from both O1 and O2.Hamburg et al. (2020) searched for GRBs coincident to the GWs reported in GWTC-1, using Fermi -GBM triggers as well as sub-threshold GRB candidates from the Untargeted and Targeted Searches, but found no additional counterparts beyond GRB 170817A.
The third observing run (O3) occurred from April 2019 to March 2020 with a month-long commissioning break during October 2019.It benefited from improvements to the sensitivity and duty cycle of the GW detectors made after O2 (Buikema et al. 2020;Abbott et al. 2021;Acernese et al. 2019).This observing run provided 56 public GW candidates in real-time with information from their preliminary analysis.More detailed analyses were published by the LIGO, Virgo, and KAGRA (LVK) Collaboration in a series of GWTCs (GWTC-2; Abbott et al. 2021, GWTC-2.1;Abbott et al. 2021, GWTC-3;Abbott et al. 2021) with GWTC-3 providing a cumulative list of 79 GW signals from O3 with a probability of astrophysical origin (p astro ) > 0.5 -an 8-fold increase relative to O2.Among these candidates was the detection of a second confident signal classified as a BNS merger, GW190425, whose total mass is larger than that known from Galactic neutron star binaries (Abbott et al. 2020).
Additionally, GW191219 163120 and GW200115 042309 provided the first detections of NSBH systems with p astro > 0.5.Another possible NSBH, GW200105 162426, fell just outside the p astro > 0.5 criterion in the GWTC-3 analysis (Abbott et al. 2021).There were also two confident detections with ambiguous classifications, GW190814 (Abbott et al. 2020b) and GW200210 092254, that represent a black hole merging with either a light black hole or a heavy neutron star.An overwhelming majority of the remaining candidates are most likely BBH in origin.
In this paper, we search Fermi -GBM and Swift-BAT data for short GRB counterparts to GW candidates from O3, discussed in Section 2.1.Section 2.2 provides an overview of the Fermi -GBM Untargeted Search as well as improvements made to the Fermi -GBM Targeted Search.Section 2.3 describes the Swift-BAT subthreshold search.We present the results of the search with Fermi -GBM triggers and the Untargeted Search in Section 3.1; with the Fermi -GBM Targeted Search, including a new joint ranking statistic that takes the spatial coincidence into account, in Section 3.2; and with the Swift-BAT sub-threshold search in 3.3.Section 3.4 presents the results from both Fermi -GBM and Swift-BAT for the marginal GWs identified in Section 2.1.Furthermore, Section 4 divides the discussion of GWs with p astro > 0.5 into two groups depending on their estimated secondary component mass m 2 .For mergers with a possible neutron star component, we present the flux and isotropic equivalent luminosity upper limits from both Fermi -GBM and Swift-BAT (Section 4.1).For the BBH mergers, we compare the lack of observed gamma-ray emission to that predicted by theoretical models (Section 4.2).We discuss upper limits to the marginal GWs in Section 4.3.Finally, in Section 5 we summarize our results and discuss future plans for using the sub-threshold searches for GWs.

METHOD
In this section, we summarize the set of GW signals that we analyze from O3.We also present the search methods used to find coincident gamma-ray and hard X-ray emission with Fermi -GBM and Swift-BAT.

GW Trigger Selection
The analysis reported here focuses on GW candidates identified during O3.These were selected by four separate analysis pipelines (i.e., GstLAL, Multi-Band Template Analysis (MBTA), PyCBC, and cWB) and published in GWTC-3 (Abbott et al. 2021).Each pipeline calculates both a false alarm rate (FAR) from a background noise hypothesis and a p astro for each candidate assuming a compact binary coalescence source.Candidate signals with p astro > 0.5 in any pipeline are selected for detailed analysis with a full estimation of the potential astrophysical source parameters.The one exception is GW candidates identified by the minimally modeled cWB pipeline, which requires a time-matched confirmation with p astro > 0.1 in one of the other pipelines in order to ensure they originated from a compact binary coalescence.In total, there were 79 GWs identified with p astro > 0.5 during O3.Table 1 shows the candidate identifier, date, time, and p astro for these GWs.
The remaining subset of GW signals with a FAR below 2 yr −1 and p astro ≤ 0.5 in a given pipeline are considered marginal GW candidates.As of GWTC-3, there are 6 marginal candidates which cannot be attributed to instrumental or environmental causes (Table 2).We exclude these candidates from our main analysis; however, since the existence of a gamma-ray counterpart could potentially prove an astrophysical origin, we perform separate searches around each marginal candidate.

Fermi-GBM Searches
Fermi -GBM has 12 sodium iodide (NaI) and 2 bismuth germanate (BGO) detectors that are strategically positioned to cover the full sky, unocculted by the Earth (Meegan et al. 2009).The flight software on-board Fermi -GBM triggers on an event when there is an influx of gamma rays at a level greater than 4.5σ above the background rate in at least two NaI detectors (Paciesas et al. 2012).Additionally, the downlink of CTTE data enables searches for GRBs below Fermi -GBM's onboard triggering threshold using ground-based computing resources.With 2 µs timing resolution and full coverage of the unocculted sky over the energy range from 8 keV to 40 MeV, CTTE data has significantly expanded the sensitivity of the Fermi -GBM instrument and its sub-threshold searches.

Untargeted Search
The Fermi -GBM Untargeted Search is a blind search that automatically scans the CTTE data for significant count rate increases in at least two NaI detectors.The algorithm was originally developed for detecting terrestrial gamma-ray flashes (Briggs et al. 2013) and has since been adapted to search for short GRBs with fluxes below the on-board triggering threshold.The Untargeted Search runs through eighteen timescales ranging from 64 ms to 31 s and five energy bins from 27 keV to 985 keV, and short GRB candidates are identified when at least two detectors exceed 2.5σ and 1.25σ above the background rate.Each candidate is given a reliability score based on whether the geometry of the detectors with significant flux is consistent with the observation of a distant astrophysical source.Currently, short GRB candidates with durations less than 2.8 s and reliability classifications of low, medium, and high are publicly distributed via GCN. 1  In this work, we combine short GRB candidates detected by the Untargeted Search with GBM-triggered GRBs and examine their temporal offsets from the GWs listed in Table 1.Theoretical models predict the tem-1 https://gcn.gsfc.nasa.gov/fermigbm subthresh archive.htmlporal offset between merger time and the production of gamma-rays to range from 0.01 s to 10 s depending on the conditions producing the gamma-ray emission (Zhang 2019).For GRB 170817A, the only known short GRB associated with a GW, the temporal offset was 1.7 s with a duration (T 90 ) over which 5-95% of the GRB flux (50-300 keV) was detected of 2 s (Abbott et al. 2017).This is consistent with a range of physically viable scenarios (e.g., Lin et al. 2018;Salafia et al. 2018;Zhang et al. 2018) where the temporal offset is correlated with burst duration.We therefore choose to subtract the burst duration timescale from the temporal offset when performing our analysis.Doing so increases the observed significance of simulated short GRB counterparts and yields no loss in detection sensitivity at the 3σ level in alternative scenarios where the temporal offset is the same for all simulated GRBs.
After calculating the time offsets for each GW-GRB pair minus the burst duration, the smallest resulting time offset for each GW is taken.For GBM-triggered GRBs, we use the T 90 as a measure of the duration.For GRB candidates from the Untargeted Search, we use the most significant timescale over which the GRB candidate was detected, which scales linearly with T 90 for on-board triggered GRBs.A background distribution is produced in the same way by replacing the observed GW times with random times during which there are no reported GW signals.This yields a distribution of temporal offsets minus the burst duration between unrelated GWs and the GRB sample.In both the search and background samples, positive and negative time offsets are allowed, with no maxima imposed.GW triggers occurring during Fermi passage through South Atlantic Anomaly (SAA) are also included.See the results presented in Section 3.1 for a comparison of the cumulative signal and background distributions.

GBM Targeted Search
The Fermi -GBM Targeted Search was developed for multimessenger follow-up observations (Blackburn et al. 2015).It uses CTTE data to scan around an external trigger time for gamma-ray emission typical of a short GRB.For follow-up of the GWs in Table 1, we search from −1 s to +30 s around the GW time to ensure we do not miss unexpectedly delayed gamma-ray emission from a counterpart short GRB, even after accounting for temporal offsets up to 10 s relative to the GW time.Starting 1 s before the GW time provides a comfortable buffer to account for the fact that the trigger times can vary by a few milliseconds for GW signals that are identified by multiple pipelines.The scan is repeated for eight characteristic emission timescales which increase The Targeted Search achieves greater sensitivity than the on-board triggering algorithm by processing the data from all 14 detectors coherently rather than focusing on significant signals present in detector pairs.This allows for the detection of weaker signals below the Fermi -GBM on-board triggering threshold (Kocevski et al. 2018).To do this, three spectral templates representing spectrally hard, normal, and soft GRBs (Table 3) are folded through the GBM detector responses to produce an expected count rate for a given astrophysical source location and flux.This expected count rate is then compared to the observed counts through a loglikelihood ratio, where di represents the background-subtracted measurements in each detector, σ n is the standard deviation of the background measurement, σ di is the standard deviation of the expected data (background+signal), r i,j is the location-dependent instrumental response for the spectrum denoted by index j, and s is the intrinsic source photon flux at the Earth.See Blackburn et al. (2015) for a full derivation.
The log-likelihood ratio quantifies the probability that an astrophysical source is present versus a backgroundonly hypothesis.It is first computed separately for each point on the sky and spectral template at a given time and emission duration.During this process we estimate the best-fit photon flux for each spectral template by finding the value s best that maximizes the log-likelihood ratio.Since the best-fit photon flux maximizes the likelihood, which is effectively a product of Gaussian dis-tributions, the variance on this photon flux equals the variance of the likelihood: where σ di includes both background and source contributions, with the latter evaluated at s best .Signal injection studies using the normal spectral template demonstrated that this formulation yields the expected error coverage levels for true source fluxes near 1 × 10 −7 erg cm −2 s −1 and below, which is the relevant flux range for this sub-threshold analysis in Fermi -GBM.We marginalize the log-likelihood ratio over all possible source amplitudes using a modified power law prior designed to both avoid divergence and to produce a luminosity distribution for the observed source flux that is invariant with respect to source distance (Blackburn et al. 2015): The net result is a hypothesis test formulation following Bayes' theorem.
The amplitude-marginalized log-likelihood ratios for individual spectral templates, L ′ j (d), at each time and duration are then averaged over all sky positions and templates using a uniform prior to formulate the full marginal log-likelihood ratio, where the sum over j covers the hard, normal, and soft spectral templates.The marginal results from all scanned times and durations are then sorted according to the largest value of Λ after filtering out known detector effects.
Localization maps estimating the probability of finding the true source location at each point on the sky are produced for the top ranking candidates using the log-likelihood ratio of the best-fitting spectrum for each candidate.This is done by noting that the log-likelihood Table 3. Spectral templates used by the Fermi-GBM Targeted Search.

Template
Type Parameters hard Cut-off Power-law (Goldstein et al. 2016) E peak = 1500 keV, α = −1.5 normal Band (Band et al. 1993) Band (Band et al. 1993) ratio asymptotically approaches the behavior of a χ 2 distribution according to Wilks' theorem (Wilks 1938) with a statistical probability given by The statistical probability is then convolved with Gaussian kernels to account for systematic errors, which are predominantly induced by the difference between the true source spectrum and the three spectral templates, imperfect knowledge of the detector response, and whether atmospheric scattering is taken into account for a given spacecraft rocking angle.As a final step, we set the region blocked by the Earth in Fermi -GBM to zero and re-normalize the map to account for the fact that gamma-ray sources are not visible through the Earth and an implicit assumption that the signal has a non-terrestrial origin.The Targeted Search method was previously used to search for sub-threshold counterparts to GWs identified during the O1 and O2 observing runs (Hamburg et al. 2020).A number of improvements were made to it in preparation for O3 (Goldstein et al. 2019): 1. Removal of the lowest 4-12 keV energy channel in the NaI detector data helped remove detector noise as well as Galactic transients.
2. Better background fitting during approach and exit from the SAA.This reduces local particle background triggers that were present in about 1% of searches and formed the dominant non-GRB background in the high log-likelihood ratio parameter space.
3. Better detector response models with a more complete treatment of the effects from the backscattering of high energy gamma-ray photons off the Earth's atmosphere.The atmospheric scattering effects are currently applied when the zenith of Fermi -GBM is within ±5 • of its nominal rocking angle of 130 • with respect to the Earth's geocenter, which occurs for ∼70% of measurements.
4. Decreasing the resolution from 1 • to 5 • for the grid of sky positions analyzed during the search.This provided an order of magnitude improvement in execution time with no notable loss in sensitivity or degradation of localization capability.
These changes necessitated a recalculation of the estimated systematic uncertainty applied to the localization maps generated for the top-ranking search candidates.
An initial study of 34 sub-threshold short GRB detections modeled this uncertainty as a 2.7 • Gaussian systematic (Goldstein et al. 2019).A more detailed model was developed for this work using a larger sample of 3,000 simulated short GRB detections.It consists of a weighted pair of Gaussian shapes normalized over the sky with the standard deviation σ 1 of the first Gaussian always smaller than that of the second Gaussian, σ 2 .The parameters of each Gaussian were determined as functions of the most probable zenith angle for each candidate and the spacecraft rocking angle relative to the Earth.They range from 1.6 • -6.0 • for σ 1 and 6.4 • -60.4 • for σ 2 , with the fractional contribution of the first Gaussian spanning 0.42-0.77.

GBM Targeted Search Ranking Statistic
We use a ranking statistic R to characterize the significance of a coincidence between the GW candidates from the catalog and the short GRB candidates found by the Targeted Search.Following the formulation in Hamburg et al. (2020), the statistic takes into account the probability of astronomical origin of the GW, p astro ; the fraction of the GW localization not occulted by the Earth for Fermi -GBM, p visible ; the time offset of the GRB candidate from the GW time, ∆t; and the FAR from the best-fitting spectral template of the GRB candidate, FAR GBM .We update the formulation to include the spatial association probability p assoc and the duration D of the gamma-ray emission: The spatial association probability p assoc quantifies whether the localizations of a sub-threshold gamma-ray candidate and GW are consistent with being produced by the same source.It is computed according to where S represents a signal hypothesis with both localizations produced by the same source and B denotes a background hypothesis where the localizations are unrelated.Both S and B are constructed from integrals over all sky positions.In this context, ρ GBM is the probability density per unit area reported by the localization maps produced for gamma-ray candidates identified by the Targeted Search.Likewise, ρ GW is the probability density of the localization maps produced for each GW and ρ uniform = 1/4π is the unit density of a uniform spatial distribution on the sky.The duration of gamma-ray emission D is incorporated into the temporal weight, where ∆t is the temporal offset between the GW time and the start of the candidate gamma-ray emission identified by the Targeted Search and D is the candidate emission timescale.As discussed in Section 2.2.1, this is designed to account for scenarios where the observed temporal offset scales with burst duration, which is expected from a broad range of models describing the observations of GW170817/GRB 170817A.The best-fit value of D, given by the candidate with the largest value of Λ, is a good proxy for burst duration because it scales proportionally with T 90 when the Targeted Search is applied to confirmed short GRBs in Fermi -GBM.We enforce a minimum value of |∆t − D| = 1 ms to avoid divergence and account for the millisecond scale uncertainty between the GW merger times of signals identified by multiple pipelines.We also apply a minimum value of FAR GBM = 6.43×10 −6 Hz.This is equal to observing a single GRB candidate over the length of the background sample used to compute FAR GBM .
We tested the impact of these updates to the ranking statistic by using the Fermi -GBM response generator2 to inject short GRBs into CTTE data from the locations of modeled BNS mergers in Abbott et al. (2020a).The start time of each GRB was offset from the GW time using the duration of each burst, as given by T 90 .We then applied the Targeted Search to this dataset and ranked the candidates according to Equation 6 as well as the older method from Hamburg et al. (2020).Doing so resulted in a factor of 1.7 increase in the number of joint detections relative to the ranking statistic formulation from our older method.
Since the true time offset model is not known, we repeated the GRB injection study using the following alternative models for the start time of the injected GRB relative to the GW: 1. Offset of half T 90 to test a scaling factor less than the total burst duration.
2. No time offset to bound emission scenarios where the GRB occurs a few ms after the GW (Zhang 2019).
3. Fixed offset of 0.5 s assuming most gamma-ray counterparts have a characteristic time delay which is half the median T 90 of short GRBs observed in Fermi -GBM.
These models were chosen with a bias towards testing time offsets shorter than the typical duration of a short GRB since the inclusion of D in the updated temporal weight naturally performs better at longer time offsets than the 1/|∆t| weight used in Hamburg et al. (2020).
The updated ranking statistic outperformed the older method in all scenarios, albeit with a smaller increase in the relative number of joint detections compared to the scenario where temporal offset scales with burst duration.

Fermi-GBM Flux Upper Limits
For GW signals without a significant counterpart detection in Fermi -GBM we compute the gamma-ray flux upper limits as a function of sky position using the Targeted Search because it is the most sensitive analysis method employed by Fermi -GBM.To do this, we use the normal spectral template from Table 3 and the 1 s gamma-ray emission duration from the Targeted Search since they are characteristic of typical short GRBs (von Kienlin et al. 2020;Poolakkil et al. 2021).This results in a set of upper limits for each sky position at times ranging from −1 s to +30 s around the GW time.We then choose the maximum observed upper limit measurement for each sky position, guaranteeing that the specified confidence level of the upper limit applies over the entire search period.
We construct the upper limits from the best-fit photon flux amplitude s best and its Gaussian error σ L discussed in Section 2.2.2 according to where N is the significance level of the upper limit.We use a 3σ upper limit level for reporting upper limits over the full 10-1000 keV energy range of standard GRB flux measurements in Fermi -GBM following the convention established in Goldstein et al. (2019).We also compute a second 5σ upper limit over a 15-350 keV range to match the convention used by Swift-BAT (see Section 2.3.3)when combining the upper limits from both instruments.

Swift-BAT Searches
Swift-BAT is a coded-aperture, large FoV (2.2 sr at 10% coding fraction), hard X-ray instrument on-board Swift.Its detector plane contains 32,768 CZT detector elements, positioned under a coded aperture mask and a graded-Z fringe shield that helps lower the background rate (Barthelmy et al. 2005).The BAT covers an energy range from 15 keV to 350 keV and monitors large portions of the sky with the goal of detecting GRBs.Once triggered, the BAT can localize a GRB to 1-3 arcmin accuracy, prompting the Swift spacecraft to slew and point its two narrow-field instruments-the X-ray Telescope (XRT; Burrows et al. 2005) and the Ultraviolet/Optical Telescope (UVOT; Roming et al. 2005)-for follow-up observations.The BAT's localization accuracy is quantified by the instrument's partial coding fraction, i.e. the fraction of detectors exposed to an event at a given time and sky position.If the coding fraction for a given trigger is 0%, then BAT will not be able to localize the event.The BAT averages ∼90 GRB on-board triggers per year, among which ∼10% are short in duration (Gehrels et al. 2009).The on-board GRB triggers are complemented by subsequent on-ground rates and image data processing, in turn allowing for dedicated searches for GRB emission.With no GWs from Table 1 triggering an on-board BAT detection, we conduct an offline follow-up analysis from the ground to search for the corresponding hard X-ray counterpart emission.

Swift-BAT Rates Data Search
The BAT flight software continuously assesses the signal-to-noise ratio (SNR) of the observed count rates.If an SNR exceeds the given threshold value determined by a number of rate-trigger criteria (Fenimore et al. 2004), the triggering algorithm subsequently checks the corresponding image data for the final confirmation and the localization of the potential burst.The detection is confirmed only if the image SNR threshold is surpassed (≳ 6.5) and no other sources have been previously reported at the event localization.For every confirmed detection, BAT records event data containing counts' arrival times, location on the detector plane, and energy.With its large effective area (∼2600 cm 2 for 100 keV photon detection at launch), the event data volume collected by BAT is too big to be stored on-board and, due to the limitations of the Swift downlink bandwidth, it is not possible to transfer all the event data to the ground.As such, until recently, the only way to conduct an offline, on-ground follow-up analysis of untriggered and sub-threshold events relied upon the rates light curves in four energy channels (15-25 keV, 25-50 keV, 50-100 keV, and 100-350 keV) with three time binnings (64 ms, 1 s, and 1.6 s) and their corresponding 64 s images in a single energy bin (15-50 keV).The recently developed Gamma-ray Urgent Archiver for Novel Opportunities (GUANO) technique circumvents this issue by retrieving BAT event data extending to ∼200 s long windows surrounding the trigger times from various astrophysical events (e.g., GWs, GRBs, fast radio bursts, neutrinos, etc.; Tohuvavohu et al. 2020).However, in this paper we do not use the GUANO technique since a significant number of the considered GW triggers were detected prior to the GUANO deployment in December 2019.Instead, we conduct the analysis using the regular rates data from BAT and leave the analogous study using the GUANO data for the future GW observing runs.
To conduct the untriggered and sub-threshold search for hard X-ray counterparts coincident with the LVK triggers in Table 1, we developed a code analogous to Lien et al. (2014).The search process begins by extracting the raw light curves from within the central region of the BAT FoV binned in 64 ms, 1 s, and 1.6 s time intervals.We opt to use the 1 s binned data to calculate the average background rate and standard deviation, σ bg , starting at −1 s before the GW trigger time, and extending to +30 s after.Using the raw light curves, we compute the average background rate, r bg , spanning a time window outside the signal interval, and spanning ∼ 800 s (excluding the instrument slews or SAA).The signal significance, S, is then computed from σ bg , using where r sig is the threshold signal rate.The background uncertainty is estimated as , where N is the number of data points in the considered portion of the lightcurve, r bg,i is the i th background rate measurement, and r bg is the mean background rate over the considered time interval.Furthermore, we visually inspect each light curve to ensure that no peaks originate from detector noise.Once this is done, we check whether there are any potential counterparts to the GW, defined as a ≥ 5σ detection above background.

NITRATES Response Functions
To produce the BAT instrument response functions appropriate for converting from photon counts to a source flux in the rate data domain, we use the Non-Imaging Transient Reconstruction And TEmporal Search (NITRATES; DeLaunay & Tohuvavohu 2021).The NITRATES response modeling takes into account both coded and uncoded parts of the detector, and thus includes responses appropriate for all counts recorded in the rates data.In addition, these responses allow for potential GRB detection from outside the BAT's coded FoV, as well as higher sensitivity across the entire FoV.The instrument responses were created by simulating photon beams onto the Swift Mass Model (SwiMM) using Geant4, a particle-interaction simulator software toolkit (Allison et al. 2016).We produce Detector Response Matrices (DRMs) for 31 different incident directions, covering the ∼ 2.2 sr sky area (corresponding to the 10-% coding fraction) where the responses are well calibrated.For a complete description of the BAT instrument response modeling, see Sec. 5 and Appendix A in DeLaunay & Tohuvavohu (2021).

Swift-BAT Flux Upper Limits
For GWs without a 5σ detection above background in Swift-BAT, we estimate the flux upper limit from the observed photon counts.We compute these limits for each time bin in the search by calculating the necessary number of counts that would result in a 5σ detection at that time from the estimated background uncertainty σ bg , assuming a 1 s emission duration.We then select the largest counts value obtained in the search bins for each GW since this is guaranteed to satisfy the 5σ criteria over the full search window.We convert the photon counts to flux upper limits within the partially coded BAT FoV, over a 15-350 keV energy range as a function of sky position by applying the NITRATES instrument response functions using the normal spectral template (Table 3) employed for upper limits computed with the Fermi -GBM Targeted Search in Section 2.2.4.This is done over 31 locations in a grid covering the ∼ 2.2 sr BAT FoV.

Combined & Marginal Flux Upper Limits
For GWs without a detected counterpart in Fermi -GBM or Swift-BAT, we combine the 5σ confidence level flux upper limits described in Sections 2.2.4 & 2.3.3 to produce joint flux upper limits as a function of sky position using both instruments.We do this by selecting the more constraining upper limit at each position since the individual limits result from independent measurements.This allows us to provide a single upper limit map for each GW that simultaneously leverages the wide FoV provided by Fermi -GBM as well as the additional coverage and enhanced sensitivity of Swift-BAT.
We also provide marginalized flux upper limits (S UL,marg ) that we compute by integrating the upper limits over the sky using the probability density of the GW localization as a prior where S UL is the position-dependent upper limit at a given confidence level and ρGW is the probability density of the GW localization normalized over the visible portion of the sky.This reduces the set of upper limits for each GW to a single, characteristic upper limit that accounts for the most likely location of the GW source.

Isotropic-Equivalent Luminosity Upper Limits
We compute upper limits on the isotropic-equivalent gamma-ray luminosity L iso in the cosmological restframe energy range of 1 keV-10 MeV for GWs without a detected counterpart in Fermi -GBM or Swift-BAT according to where D L is the median luminosity distance of the GW, S UL,marg is the marginalized confidence level flux upper limit described in Section 2.4, and k is the standard bolometric correction factor given by k ≡ 10 MeV/1+z 1 keV/1+z where z is the redshift inferred from D L .In this case, dN dE (E) represents the normal spectral shape from Table 3, which is used to generate the marginalized flux upper limit.We chose to use the median D L and the marginalized flux upper limit for each GW rather than marginalizing L iso from the values estimated at each sky position in order to exclude the lowlikelihood modes of the distance luminosity posteriors for GW200308 173609 and GW200322 091133, as discussed in Abbott et al. (2021).All other GWs yield similar values for L iso regardless of whether we use the individual or median values for D L in the calculation.

RESULTS
This section presents the results for the searches from Fermi -GBM and Swift-BAT for coincident gamma-ray emission to the GW candidates presented in Tables 1  and 2.

Triggered and Untargeted Search Results
We compare the distribution of temporal offsets between the GWs listed in Table 1 and the closest gammaray signal in Fermi -GBM, minus its burst duration, shown in Figure 1.The GRB sample here comprises 214 GRBs triggered by Fermi -GBM during O3 and 479 short GRB candidates detected by the Untargeted Search.The background distribution was composed by choosing ∼ 10 4 random times in O3 during which there were no reported GWs.Confidence intervals for the search sample were determined empirically by Monte Carlo sampling of the background distribution.
There are no significant deviations from the background, with the search sample lying largely within the 68% confidence interval.The shortest temporal offset given by the sample distribution is approximately 10 minutes and is within the 95% confidence interval.Such large offsets are not expected for on-axis prompt emission from GRBs associated with binary mergers (Vedrenne & Atteia 2009; Zhang 2019); therefore, we find no further evidence for a GW/gamma-ray association.

Targeted Search Results
We ran the Targeted Search on the 79 GWs shown in Table 1.Fermi -GBM was in the SAA for 15 of those times, therefore the detectors were turned off and no data were obtained.Figure 2 shows the cumulative rate above a given value of the marginalized log-likelihood ratio Λ separated according to the best-fitting spectral template.The background distribution is determined by randomly selecting times during the O3 livetime without known GW triggers.This represents the FAR of the search and describes the frequency of expected false positives as a function of Λ.No significant gamma-ray signals were found in coincidence with GWs.
Figure 3 shows the ranking statistic R from Section 2.2.3 mapped to a p-value, defined as the number of more highly ranked background candidates as compared to the total number of background candidates or p i = N (R > R i )/N , where N is the number of background gamma-ray candidates and i is the candidate index within the search sample.The plots show no signif-icant deviations from background, yielding no evidence for a GRB counterpart to the GW signals.

Swift Results
We ran the pipeline described in Section 2.3 on the 1 s binned light curves from Swift-BAT for the 79 candidates listed in Table 1.The goal of the pipeline is to examine whether any emission surpasses the 5σ threshold above the observed background rate level.We visually inspect each light curve to ensure that no detector noise or malfunction affects the reported results.We identify detector noise as a fast duration peak seen only in some of the energy channels.Once identified, the detector noise is subtracted down to the level of the average background rate.There are 14 GWs for which the data are either unavailable or background dominated because they occurred during either the Swift-BAT SAA passage or a slew.Separately, 13 GWs were within the BAT FoV (>10% partial coding) and had image survey data available at the time of interest.We report no significant hard X-ray detection in the Swift-BAT rate data coincident with the reported GW triggers at the 5σ level.We also ran the standard BAT analysis on the survey data (longer timescales, ∼ 300 s) for the 13 inside-BAT-FoV candidates and also report no excess X-ray emission.3

Results for Marginal GW Canidates
The Fermi -GBM and Swift-BAT searches were applied separately to the 6 marginal GW signals from Table 2 in an effort to identify EM counterparts, that could prove an astrophysical origin.The closest time offset from the GBM Triggered and Untargeted Searches was observed for GW191118 212859, which occurred 42 minutes before the on-board trigger of GRB 191118A and corresponds to a p-value of 0.1.The most significant Targeted Search candidate was found for GW200105 162426 using the hard spectral template.It has a pre-trial p-value of 0.1 as estimated by the ranking statistic distribution for background described in Section 2.2.3.Applying a trials factor of 3 to account for the 3 spectral templates used by the Targeted Search increases the p-value to 0.3.No 5σ detections were found for any marginal candidates using the Swift-BAT rates data search.

SCIENCE DISCUSSION
Compact binary mergers containing a neutron star component are likely candidates for gamma-ray emis- The cut m 2 ≤ 3 M ⊙ was chosen to include all systems with at least one neutron star component up to the maximal allowed neutron star mass of 2.16-3.0M ⊙ (Bombaci 1996;Kalogera & Baym 1996;Margalit & Metzger 2017;Rezzolla et al. 2018).It may include a few ambiguous BBH mergers due to the uncertainty on the maximum allowed neutron star mass.We favor this approach over a stricter cut due to the limited number of systems in O3 with light secondary component masses.Additionally, the discussion of possible BBH mergers does not suffer from the loss of a few ambiguous candidates, particularly given the large number of systems with m 2 > 3 M ⊙ .
In Section 4.1 we discuss the absence of GRB detections in coincidence with the 6 GWs classified with a possible neutron star component and present flux upper limits from Fermi -GBM and Swift-BAT.In Section 4.2 we discuss the BBH mergers, providing flux upper limits and exploring how these limits may rule out certain theoretical models.−2.8 M ⊙ and the lowest secondary component mass m 2 = 1.17 +0.07 −0.06 M ⊙ of all the GWs with a possible neutron star (Abbott et al. 2021); it represents a potential NSBH merger.GW200115 04230 has a primary mass of m 1 = 5.9 +2.0 −2.5 M ⊙ suggesting a low-mass black hole, and a secondary mass of m 2 = 1.44 +0.85  −0.29 M ⊙ which is consistent with a neutron star (Abbott et al. 2021).GW200210 092254 possesses a primary component mass of m 1 = 24.1 +7.5 −4.6 M ⊙ and a secondary component mass of m 2 = 2.83 +0.47 −0.42 M ⊙ that could either be a light black hole or a heavy neutron star (Abbott et al. 2021).It is unclear if this source is a BBH or NSBH merger.

Possible Neutron Star in the System
All GWs except GW191219 163120 were observable by both Fermi -GBM and Swift-BAT (Table 4).GW191219 163120 was observed by Fermi -GBM but no Swift-BAT data are available due to a slewing behavior of the spacecraft at the GW time.Neither instrument  Table 4. Flux upper limits for the 6 GWs from O3 with pastro > 0.5 that are classified with a possible neutron star component.The 3σ upper limits are computed for the 10-1000 keV energy range over the FoV of Fermi-GBM.The 5σ upper limits are computed for the combined coverage of Fermi-GBM and Swift-BAT with both instruments matched to the 15-350 keV energy range of Swift-BAT.The columns labeled Min and Max correspond, respectively, to the minimum and maximum upper limits obtained for points within the 90% credible region of the GW localization.The Marginal upper limit is computed by integrating the upper limits produced at individual locations over the full sky using the GW localization as a weighted prior, normalized to the visible portion of the sky.Also shown is the visible coverage percentage of the full GW localization for Fermi-GBM alone, Swift-BAT alone, and the combined FoV from both instruments.detected an EM counterpart to these GWs.Therefore, we compute flux upper limits for each GW using the methods described in Section 2. Table 4 presents the minimum and maximum flux upper limits from Fermi -GBM and Swift-BAT over the 90% credible regions of the GW localizations, as well as sky-marginalized flux upper limits.Incorporating the combined measurements from Fermi -GBM and Swift-BAT, we also use the sky-marginalized 5σ flux limits to generate upper limits on the isotropic-equivalent luminosity (Figure 5).The combined 5σ flux upper limits over the 15-350 keV energy range can be seen in Figure 6.
The lack of a gamma-ray counterpart to BNS merger GW190425 has three plausible explanations.First, the combined coverage of ∼60% of the total GW localization implies that the GW source may not have been visible to Fermi -GBM and Swift-BAT.Second, GW190425 has an estimated luminosity distance of D L = 0.15 +0.08 −0.06 Gpc which is 4 times larger than that to GW170817 (Abbott et al. 2019).At this distance, the luminosity of GW170817 would fall well below the upper limit for GW1904254 (Fletcher et al. 2019).(Figure 5), indicat-ing that a counterpart similar to the one for GW170817 would have been unobservable to Fermi -GBM or Swift-BAT.Finally, the inclination angle of this GW is poorly constrained, with the 90% credible level extending to a viewing angle of 70 • with respect to the jet axis.This encompasses scenarios where the observed off-axis flux would be below the detection limits, even if the central engine of GW190425 was powerful enough to be detected on-axis by Fermi -GBM and Swift-BAT at its measured distance.

Probable BBH Mergers
There are a total of 73 GWs with the criterion of m 2 > 3 M ⊙ in > 95% of the posterior probability.All of these have estimated primary and secondary component masses much larger than the maximum expected neutron star mass of 3 M ⊙ .Therefore, they are most likely GW signals from BBH mergers.
Of these GWs, 10 occurred during SAA for Fermi -GBM, but had data from Swift-BAT.Likewise, Swift-BAT does not have data for 9 GWs, either because Swift-BAT was in the SAA or slewing, but data are available from Fermi -GBM.Finally, there are 5 GWs that do not have data from either Fermi -GBM or Swift-BAT due to being in the SAA and/or slewing.Neither instrument identified an EM counterpart for the GWs with data coverage.As a result, we compute flux upper limits for each GW according to the methods described in Section 2. Table 5 presents the minimum and maximum upper limits over the 90% credible region of the GW localization as well as the marginalized flux upper limits.Joint flux upper limits as a function of sky position and the corresponding isotropic-equivalent luminosity limits for the GWs that have data coverage are provided in a separate data release (Wood et al. 2023).For the GWs with Fermi -GBM data, we look into constraining theoretical models of gamma-ray emission from BBHs using the 3σ flux upper limits, as they provide broad spectral coverage over the 10-1000 keV energy range.

Constraining gamma-ray emission models from BBH mergers
EM radiation from BBH mergers is not expected due to the challenges associated with forming an accretion disk during the merger process.Nevertheless, Connaughton et al. (2016) reported GW150914-GBM, a weak gamma-ray signal following the first LIGO-Virgo detection of BBH GW150914 (Abbott et al. 2016), and, preliminary Fermi-GBM flux upper limits reported during the O3 online analysis more recently, the Zwicky Transient Facility identified a potential EM counterpart to GW190521 (Graham et al. 2020).While the associations between these detections and the corresponding GWs remain nebulous (Connaughton et al. 2018;Ashton et al. 2021;Bustillo et al. 2021;De Paolis et al. 2020;Palmese et al. 2021), a wide spectrum of models have been developed to invoke EM emission from BBH-mergers, all with non-negligible difficulties (e.g., Loeb 2016b; Perna et al. 2016;Zhang 2016;Perna et al. 2018).
To test the association between GW150914 and GW150914-GBM, Veres et al. (2019) assumed some of these BBH emission models and derived a modeldependent BBH-to-GRB ratio which represents the expected number of BBH mergers to be detected by LVK before a gamma-ray counterpart might be observed by Fermi -GBM.Since the number of LVK BBH merger detections has reached the BBH-to-GRB ratio for a few models reported by Veres et al. (2019), we attempt to constrain them by computing gamma-ray flux upper limits for each model and examining the implications with respect to individual BBH mergers.
We consider four models for relating potential gammaray emission to the energy present in BBH mergers: a neutrino-antineutrino annihilation powered jet mechanism (ν ν; Ruffert & Janka 1998), a charged BH (Q; Zhang 2016), the Blandford-Znajek mechanism (BZ; Blandford & Znajek 1977), and a model where the gamma-ray energy is proportional to the emitted GW energy (E GW ).A detailed summary of these models and their parameters can be found in Veres et al. (2019).We note that all of the above scenarios suffer from non-trivial critiques but are used here to be widelyinclusive of the broad spectrum of proposed mechanisms for gamma-ray production.The intrinsic properties of each model (e.g., magnetic field strength, charge of the black hole, etc.) are determined by setting them to values consistent with the observed luminosity of GW150914-GBM (Connaughton et al. 2016).
For each GW, we use the posterior distributions of BBH parameters (e.g., final mass, distance, inclination, rotation parameter, etc.) from GWTC-3 and derive a distribution of gamma-ray fluxes for the different models.We then compare the distribution of fluxes to the 3σ marginalized flux upper limit from Fermi -GBM, shown in Table 5.This is performed for three GRB jet geometries: an isotropic emitter (i.e., an opening angle of 90 • ), an opening angle distributed uniformly between 10-40 • , and a fixed 20 • opening angle.All jets are assumed to have a top-hat angular structure and are assigned pointing directions by sampling from the inclination posterior of each GW.Due to relativistic beaming, emission is   Example of gamma-ray flux expected for the 4 different models: charged BH model (Q), neutrino-antineutrino annihilation powered jet model (ν ν), Blandford-Znajek (BZ), and gamma-ray energy as a fraction of GW energy (EGW) in the case of GW191216 213338.The 3σ (10-1000 keV) marginalized flux upper limit is indicated by the vertical line.The legend contains the fraction of cases above the 3σ upper limit; note that for the jetted emission models (uniform, fixed), cases with zero expected flux are not shown.Here the limit is violated in 100% and >99% of cases for the ν ν and Q models, respectively (assuming isotropic emission).
strongly suppressed for GRBs with jet opening angles smaller than the viewing angle.In order to simplify the treatment of such cases, we assign zero flux to jets whose inclination is larger than the opening angle.
Figure 7 shows an example of the flux distribution for the four different models and the three jet geometries compared to the GBM upper limits for GW191216 213338.There is a dearth of jetted cases (green, red) compared to isotropic emission (light blue) at higher fluxes.This is explained by the inclination angle-distance degeneracy of the GW parameter estimation.Point estimates with smaller distances and thus higher flux will preferentially have jets pointed away from our line of sight.When we impose a jet opening angle ≲ 40 • on such systems, they will not include the observer within their aperture in most of the cases.Conversely, the point estimates with largest distances point preferentially towards the observer, thus there will be no strong differences between the jetted and isotropic cases at low flux values.
We classify GW191216 213338 as noteworthy because the predicted gamma-ray flux distribution from at least one model violates the GBM flux upper limit by more than 10%.The 10% limit in the rest of this section refers to the isotropic emission model.Of the 58 probable BBH mergers with Fermi -GBM data coverage described in Section 4.2, 18 are considered noteworthy according to this criterion.The remaining 40 did not yield the necessary number of cases above the GBM flux upper limit in any of the models.
Out of the four models considered here, the ν ν model violates the Fermi -GBM flux upper limit in most of the cases.Of the noteworthy GWs (denoted by * in Table 7, Appendix A), 15 exceed the GBM limit in more than 10% of cases for this model.In particular, GW190924 021846, GW191216 213338 (Figure 7) and GW200202 154313 exceed the GBM limit in {100, 30, 20}%, {100, 31, 19}% and {100, 36, 24}% of the cases respectively (the 3 numbers represent the 3 different opening angle choices).Interestingly, for these three events, in the isotropic emission scenario the ν ν can be ruled out, and the non-detection in gamma-rays can constrain the jet geometry.This is due to GW191216 213338 and GW200202 154313 being the two closest BBH signals observed during O3 with luminosity distances of 0.34 +0.12 −0.13 Gpc and 0.41 +0.15 −0.16 Gpc, respectively (Abbott et al. 2021).In addition to being nearby (luminosity distance of 0.55 ± 0.22 Gpc), GW190924 021846 has a relatively low final mass (13.9 +2.8 −0.9 M ⊙ ), which leads to higher flux in the ν ν model.
For the Q model, the most constraining GWs are also GW191216 213338 and GW200202 154313.They violate the GBM flux upper limit for the different jet geometries in {100, 31, 19}% and {96, 33, 22}% of cases, respectively.In total, 12 of the probable BBH mergers have larger than 10% of their flux estimates above the upper limit for this model.
GW191109 010717 is the most constraining for the BZ scenario.It violates the gamma-ray upper limit in {26, <0.1, <0.1}% of cases for the 3 different opening angle choices.It has the fourth-highest total mass M = 112 +20 −16 M ⊙ in O3 and is reasonably close at D L = 1.29 +1.13  −0.65 Gpc (Abbott et al. 2021).The only other GW with more than 10% of the flux estimates above the upper limit for the BZ mechanism is GW190521 074359.
For the E GW scenario, GW191216 213338 is again the most constrained.It violates the gamma-ray upper limit in {27, 1.8, 0.7}% of the cases for the three jet opening angle choices.There are 3 GWs with 10% of the flux estimates from this model above the GBM flux upper limit (Table 7 continued, Appendix A).
In summary, we provide constraints on theoretical models of gamma-ray emission from BBH mergers using the flux upper limits from Fermi -GBM.We find that for most BBH mergers the models considered here do not predict gamma-ray flux over the upper limit.Under our model assumptions, this can be understood as a consequence of the larger average distance of BBH mergers during O3 compared to that of GW150914.We also find that the ν ν model is the most constrained.This model has the lowest BBH-to-GRB ratio in Veres et al. (2016), and indeed, observations reveal that 18 out of 58 cases for the ν ν model yield an appreciable flux above the upper limit.The expected flux in this model is inversely proportional with the final mass; the average BBH merger in O3 was less massive than GW150914, resulting in larger predicted gamma-ray flux.

Marginal GWs
Although all 6 marginal GWs from Table 2 have p astro < 0.5 they are of interest for EM follow-up.This is because GW200311 103121 may have a possible BNS origin and the remaining 5 candidates have possible NSBH origins (Abbott et al. 2021;Abbott et al. 2021).In particular, the possible NSBH merger GW200105 162426 was noted as a clear outlier from experimental backgrounds despite not satisfying the p astro > 0.5 criteria used to identify GW signals with a likely astrophysical origin.It also has the highest observed p astro of all the marginal GWs.
The 5 marginal GWs with possible NSBH origins were visible to Fermi -GBM while the remaining one, GW200311 103121, occurred when Fermi -GBM was in the SAA.None of the marginal GWs have appreciable coverage in Swift-BAT.No significant counterparts were found.As with GW190425, this may be due to unfavorable viewing angles with respect to the jet axis, larger observational distances such as the 0.27 +0.12 −0.11 Gpc distance to GW200105 162426 (Abbott et al. 2021), and partial sky coverage for candidates other than GW190426 152155.It therefore remains ambiguous as to whether these signals are real compact binary coalescences.Nevertheless, we provide in Table 6 the flux upper limits for each marginal GW calculated according to the same methods described in Section 4.1 since they may provide emission model constraints if future analyses can identify an astrophysical progeni-tor with a favorable viewing angle with respect to the jet axis.Figure 8 displays the 5σ confidence level flux upper limit map for GW200105 162426 since it is the marginal GW with the highest probability of having an astrophysical origin.The marginalized 5σ flux upper limit of GW200105 162426 yields an isotropic-equivalent luminosity upper limit of L iso = 2.1 × 10 48 erg s −1 when combined with its 0.27 Gpc distance.The data release (Wood et al. 2023) associated with this work provides flux upper-limits as a function of sky position for the remaining marginal GWs.

SUMMARY AND FUTURE DIRECTIONS
Using the 79 GW candidates with p astro > 0.5 from O3 that were reported in GWTC-3, we searched for coincident EM counterparts with Fermi -GBM and Swift-BAT.This represents the most comprehensive follow-up to date of the O3 run in the hard X-ray and gamma-ray regime.We found no significant counterparts in either instrument.For the one BNS merger, GW190425, with p astro > 0.5 there are several possible reasons for the non-detection of a counterpart: • The combined Fermi -GBM and Swift-BAT coverage of the GW localization area was ∼60%, meaning that the GW source may have been outside the FoV for both instruments.
• The distance to GW190425 was four times larger than the estimated distance for GW170817, causing the observed flux to be below the detection threshold in both instruments if it had the same intrinsic luminosity and viewing angle as GW170817.
• The viewing angle may have been too far away from the jet axis to detect emission under scenarios with a clean or structured jet.
GW190425 is therefore unconstrained by our observations.
In contrast to GW190425, the large number of BBH detections in this sample allowed us to begin placing constraints on certain models of gamma-ray emission from BBH mergers, despite the larger average distance for this class compared to the BNS mergers.The most constrained model was the ν ν model where two GWs, GW191216 213338 and GW200202 154313, were predicted to produce an observable flux in Fermi -GBM.With the number of GWs from BBH mergers to increase in the fourth LVK observing run (O4), we expect this model to become more constrained and ruled out as a potential explanation of EM emission from BBH mergers.With O4 having begun on May 24, 2023, we expect an increase in GW detections by a factor of 5 (Abbott et al. 2018;Petrov et al. 2022) and a bountiful regime of EM follow-up data.This will greatly increase the need for further instantaneous, wide FoV gamma-ray/X-ray observations in order to detect the EM counterparts to these GWs and localize them, especially given the absence of a counterpart detection during O3.Towards this end, the Fermi -GBM Targeted Search updates presented in this paper will be used in future LVK observing runs.In the absence of detections, flux upper limits, both marginalized and as a function of sky position, will be provided to the community.Additionally, Swift-BAT GUANO and NITRATES will be used during the next observing run.

Figure 1 .
Figure 1.The cumulative distribution for the minimum time offsets between the O3 GW triggers and GRBs found by either the GBM on-board triggering algorithms or the Untargeted Search.Confidence intervals for the search sample were determined by Monte Carlo sampling of the background.

Figure 2 .
Figure2.The cumulative rate above a given value of the marginalized log-likelihood ratio Λ, separated into three plots according to the best fitting spectral template, found in the Targeted Search.The orange line is the foreground distribution of GRB candidates found with the Targeted Search around the given GW merger time.The black dotted line is the randomly selected background sample and the green shading represents the 68%, 95% and 99.7% confidence intervals around it.sion,particularly if the inspiral process results in tidal disruption of the neutron star(Burns 2020).In contrast, BBH mergers are not expected to produce gamma-rays outside of a few exotic scenarios (e.g., Loeb 2016a;Perna et al. 2016;Zhang 2016;Dai et al. 2017).Therefore, using the standard convention of m 1 > m 2 , we divide our discussion into two sections based on the secondary component mass m 2 : 1. Mergers with a possible neutron star: m 2 ≤ 3 M ⊙ (5% credible level)

Figure 3 .
Figure 3. Cumulative fraction versus p-value of the updated ranking statistic R. The solid black line is the foreground distribution of GRB candidates found with the Targeted Search around the given GW merger time.The dashed black line is the randomly selected background sample and the blue, green and orange lines represent the 68%, 95% and 99.7% confidence intervals for background, respectively.There are 6 GWs (i.e., GW190425, GW190814, GW190917 114630, GW191219 163120, GW200115 042309, and GW200210 092254) where ≥ 5% of posterior probability lies below the dashed red line in Figure4.GW190425 is the least massive system from O3 and the second BNS merger detected by LIGO-Virgo.It has a primary mass m 1 = 2.1 +0.5 −0.4 M ⊙ and a secondary mass of m 2 = 1.3 +0.3 −0.2 M ⊙(Abbott et al. 2021(Abbott et al. , 2020).GW190814 has a low mass secondary component, estimated at m 2 = 2.6 +0.1 −0.1 M ⊙ , while its primary component has an estimated mass of m 1 = 23.3+1.4−1.4 M ⊙ .It is unclear whether this source is a BBH or a NSBH merger, since the secondary component could either be a light black hole or a heavy neutron star(Abbott et al. 2021).GW190917 114630 was identified as a BBH merger by the GstLAL pipeline, but its secondary component mass of m 2 = 2.1 +1.1 −0.4 M ⊙ is a strong indicator for a NSBH origin(Abbott et al. 2021).GW191219 163120 has a large

Figure 4 .
Figure 4.The inferred 90% credible regions of the component masses for all GWs with pastro > 0.5 from O3 are shown in grey (Abbott et al. 2021; Abbott et al. 2021; Abbott et al. 2021).The red dashed line marks the upper bound of m2 allowed for our classification of systems with a possible neutron star component, which are marked by colored contours.

Figure 5 .
Figure5.The 5σ upper-limits on isotropic-equivalent luminosity Liso for the 6 GWs in O3 identified with a possible neutron star component and pastro > 0.5.The black data point is the measured Liso from GRB 170817A and the black dashed line is the approximate lower bound for Liso of GRBs detected on-board Fermi-GBM(Abbott et al. 2017).

Figure 6 .
Figure 6.The 5σ flux upper-limit as a function of sky position for the 6 GWs from O3 identified with a possible neutron star component and pastro > 0.5.The purple gradient represents the combined Fermi-GBM and Swift-BAT flux upper limits for source positions at each point on the sky.The star symbol represents the zenith direction of Fermi-GBM, the square symbol represents the center of the Swift-BAT FoV, and the green contour represents the 90% credible area of the LVK localization.The blue region is the non-visible portion of the sky which is occulted by the Earth for Fermi-GBM and outside the Swift-BAT FoV.

Figure
Figure 5 continued.
Figure 7.Example of gamma-ray flux expected for the 4 different models: charged BH model (Q), neutrino-antineutrino annihilation powered jet model (ν ν), Blandford-Znajek (BZ), and gamma-ray energy as a fraction of GW energy (EGW) in the case of GW191216 213338.The 3σ (10-1000 keV) marginalized flux upper limit is indicated by the vertical line.The legend contains the fraction of cases above the 3σ upper limit; note that for the jetted emission models (uniform, fixed), cases with zero expected flux are not shown.Here the limit is violated in 100% and >99% of cases for the ν ν and Q models, respectively (assuming isotropic emission).

Figure 8 .
Figure 8.The 5σ upper-limits as a function of sky position for GW200105 162426, the marginal GW with the highest pastro.The purple gradient represents the combined Fermi-GBM and Swift-BAT flux upper limits for source positions at each point on the sky.The star symbol represents the zenith direction of Fermi-GBM, the square symbol represents the center of the Swift-BAT field-of-view, and the green contour represents the 90% credible area of the LVK localization.The blue region is the non-visible portion of the sky which is occulted by the Earth for Fermi-GBM and outside the Swift-BAT FoV.

Table 2 .
Marginal GWs from O3 without clear instrumental or environmental causes(Abbott et al. 2021; Abbott et al. 2021)