HIGH-ENERGY ELECTROMAGNETIC OFFLINE FOLLOW-UP OF LIGO-VIRGO GRAVITATIONAL-WAVE BINARY COALESCENCE CANDIDATE EVENTS

The initial LIGO (Abbott et al. 2009) and Virgo (Accadia et al. 2012) gravitational-wave (GW) detectors took their last science data between 2009 July and 2010 October before going offline for several years for major upgrades to advanced detector configurations (Acernese 2009; Harry 2010). For its sixth science run (S6), LIGO consisted of two 4 km baseline interferometric detectors—one instrument H1 at Hanford, WA, and another instrument L1 at Livingston, LA. The Virgo 3 km interferometer V1 in Cascina, Italy took concurrent data during their second (VSR2 covering 2009 July–2010 January) and third (VSR3 covering 2010 August–October) science runs. Together, they represent the most sensitive worldwide network of GW detectors to date. Figure 1 shows typical strain sensitivity of each instrument as a function of frequency, as well as the distance at which an optimally oriented merger of two compact objects (black holes or neutron stars) would produce a nominal signal-to-noise ratio (S/N) of eight in each detector.

Zoom In Zoom Out Reset image size

Compact binary coalescence (CBC) is the most anticipated GW source for first and second-generation ground-based detectors. These systems are very strong emitters of GWs, and we are confident of their existence through the discovery of a handful of galactic NS/NS binary systems where one object is a radio pulsar which is modulated by the orbit—the most famous of these systems is the Hulse–Taylor pulsar PSR1316+16 (Weisberg & Taylor 2005). The number and lifetime of these systems can be used to obtain an estimate of about one NS/NS merger event per Mpc³ per million years (Abadie et al. 2010), with up to two orders of magnitude uncertainty due largely to our limited knowledge of the pulsar luminosity function and limited statistics. The merger rate translates into an estimate of ∼0.02 detectable NS/NS merger events per year for the initial detector network, and ∼40 per year for the advanced detectors once they reach design sensitivity. NS/BH systems are also of interest for electromagnetic (EM) follow-up and are stronger GW emitters. However we have not yet observed any NS/BH binary systems, and have generally poor knowledge of the black hole mass distribution.

While no NS/NS or NS/BH binary coalescence events have been detected during S6/VSR2+3 using GW data alone (Abadie et al. 2012), it is conceivable that an EM counterpart might allow us to resolve a rare event otherwise too weak to distinguish from background in the first generation detectors. In the rapidly approaching advanced detector era, searching for EM counterparts can play an important role in increasing our confidence of detection for otherwise marginal events and will provide astrophysical context for GW detections.

Advanced LIGO is likely to begin taking its first science data in 2015, with Virgo following a year later (Aasi et al. 2013). As they reach design sensitivity, the advanced detectors are expected to begin an era of regular detection of GW events, making the search for EM counterparts triggered by GW detections an enticing possibility. In this study, we demonstrate strategies for searching high-energy archival EM data for counterparts using representative GW background events from real data taken during the most recent LIGO-Virgo science runs and high-energy EM survey data recorded at the same time. We then measure EM detection efficiencies under various plausible emission models and the probability of accidental coincidence given realistic GW background and sky-localization accuracy.

The two instruments chosen for this EM follow-up study include the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) on board the Fermi spacecraft, and the soft X-ray All-sky Monitor (ASM; Levine et al. 1996) on board RXTE. Both are particularly well suited to an offline follow-up of GW events because of their large, regular coverage of the entire sky (Table 1), and because they save a large amount of archival survey data which allows for sensitive offline searches.

The high-energy sky itself offers the advantage of being relatively clean compared to optical wavelengths, which is important for maintaining a low probability of false-coincidence given the large GW sky-location uncertainty. The use of offline survey data for follow-up has many practical differences from other triggered searches in the EM spectrum for GW counterparts. (Evans et al. 2012; Kanner et al. 2012; Aasi et al. 2014). The offline search does not rely on rapid GW data analysis and continuous coordination with EM observational facilities. An important consequence is that the offline coincidence can generally be applied to many more events, allowing the follow-up of weaker, more marginal candidates.

Gamma-ray bursts (GRBs) are flashes of gamma-rays observed approximately once per day. Their isotropic distribution in the sky was the first evidence of an extra-galactic origin, and indicated that they were extremely energetic events. The duration of prompt gamma-ray emission shows a bi-modal distribution which naturally groups GRBs into two categories (Kouveliotou et al. 1993). Long GRBs emit 90% of their prompt radiation on timescales longer than 2 s. They have been associated with the collapse of rapidly rotating massive stars (Hjorth & Bloom 2012). sGRBs with prompt emission <2 s and a generally harder spectrum are thought to arise from the merger of two neutron stars, or a neutron star and black hole (Nakar 2007). It is this favored progenitor model which makes sGRBs and associated afterglow emission a promising counterpart to GW observations.

The Swift satellite has revolutionized our understanding of short GRBs over the last several years by the rapid observation of X-ray afterglows, providing the first localization, host identification, and redshift information. The beaming angle for short GRBs is highly uncertain, although limited observations of jet breaks in some afterglows imply half-opening angles of ${{\theta }_{j}}\;\sim$ 3°–14° (Liang et al. 2008; Fong et al. 2012). The absence of an observable jet break sets a lower limit on the opening angle which is generally weak (due to limits in sensitivity), though in the case of GRB 050724A, late-time Chandra observations were able to constrain ${{\theta }_{j}}\gtrsim 25{}^\circ$ (Grupe et al. 2006). The observed spatial density of short GRBs and limits on beaming angle result in an NS/NS merger event rate roughly consistent with that derived from galactic binary pulsar measurements.

Although the beaming factor of $\sim \theta _{j}^{2}/2$ means most merger events seen by the advanced GW detectors will not be accompanied by a GRB, this is somewhat compensated by the fact that the ones that are beamed toward us have stronger GW emission. Current estimates for coincident GW-sGRB observation for advanced LIGO-Virgo are a few per year assuming an NS/NS progenitor model (Coward et al. 2012; Metzger & Berger 2012). The rate increases by a factor of eight if all observed short GRBs are instead due to NS/BH (10 ${{M}_{\odot }}$) mergers which are detectable in GWs to twice the distance. Coincident rates for the initial detectors go down by a factor a thousand due to their factor of ten worse sensitivity.

About 80% of short GRBs seen by Swift are accompanied by some kind of X-ray afterglow (Gehrels & Razzaque 2013), which are often observable for about a day. The observed X-ray afterglows are thought to occur from synchrotron emission at the shock front where the outgoing jet meets the local inter-sellar medium. Simulations show that such afterglow emission becomes very weak off-axis, with little possibility of detection at twice the opening angle (van Eerten & MacFadyen 2011). Searches for orphan afterglow signals, without the presence of detected prompt GRB emission, are quite difficult to confirm due to various sources of transient background. However with incomplete coverage of the gamma-ray sky, it is quite possible that the first GW-EM association will be with an afterglow signal.

In addition to the jet-driven burst and afterglow, other EM emission associated with a compact merger can be a promising channel for GW-EM coincidence, particularly if the EM radiation is less-beamed or even isotropic. A few short GRBs (∼10%) have shown clear evidence of high-energy flares which preceed the primary burst by 1–10 s, and possibly up to 100 s (Troja et al. 2010). The precursors can be interpreted as evidence of some activity during or before merger, such as the resonant shattering of NS crusts (Tsang et al. 2012), which depends on the NS equation-of-state and could radiate isotropically. A class of short GRBs (∼25%) also contain a period of extended soft X-ray emission on a timescale of ∼10–100 safter the initial spike, which is difficult to explain with the standard jet scenario. It has been proposed that the extended emission could arise from a relativistic wind powered by a short-lived rapidly rotating protomagetar star surviving post-merger (Metzger et al. 2008; Gao et al. 2013; Nakamura et al. 2013; Zhang 2013), which could also be considerably more isotropic than the jet afterglow. Since short GRBs themselves may be subject to a large beaming factor, the more numerous nearby NS/NS mergers seen in GWs will provide a unique opportunity to search for these potentially weak exotic EM phenomena, and isolate them from jet behavior.

Detection of the characteristic signature of a CBC in GWs prior to an associated sGRB will provide unambiguous support for the compact merger progenitor model. Along with this critical piece of the sGRB puzzle, GWs provide a largely complimentary set of information: component masses (whether black holes are involved) and spins, system inclination, luminosity distance; while the EM counterpart provides information about EM energetics, a precise location, local and host environment, and redshift. This makes joint GW-EM detections particularly valuable and a key goal in GW astronomy.

GW-EM coincidence can be approached from a variety of angles. In this analysis, we use a collection of all-sky, all-time GW events to conduct a targeted search in high-energy EM archival survey data. Another strategy demands real-time analysis and localization of the GW data (Abadie et al. 2012; Cannon et al. 2012; Singer et al. 2014), suitable for targeted observation with narrow-field instruments with the goal of catching a short-lived afterglow signal. This rapid EM-follow-up approach was succesfully tested in initial LIGO/Virgo for a handful of online GW events using both a collection of ground-based optical telescopes (Aasi et al. 2014) and Swift-XRT (Evans et al. 2012). The online strategy places heavy demands on both the rapid processing and interpretation of GW data, as well as the capability of pointed telescopes to effectively cover the large GW-determined sky location (Abbott et al. 2012; Kanner et al. 2012; Singer et al. 2013). In return they provide the deepest observations, with sufficient sensitivity to observe, for example, the faint but promising optical/IR kilonova merger counterpart (Metzger & Berger 2012; Kasliwal & Nissanke 2014).

Known EM transients can also trigger specialized searches in the GW data itself, where the precise time, sky-location, and potentially model information from the EM event may greatly constrain the space of GW signals to look for. The reduced search can be performed at a lower threshold in the GW data due to both computational and statistical (false-alarm probability) considerations (Was et al. 2012). In the case of known GRBs, such searches have been used to place minimum distance limits dependent on GW emission model (Abadie et al. 2012a; Aasi et al. 2014), and for two short GRBs have been able to uniquely exlude the possibility of a binary merger progenitor from a plausible nearby host (Abbott et al. 2008; Abadie et al. 2012b).

For this anaysis, the offline coincidence between GW and EM observation is done in time and sky location, providing an ability to exclude non-coincident background from either search. Timing for a GW signal is very good, and model spiral waveforms can be placed to fractions of a millisecond (Fairhurst 2009). The coincidence search window in time is determined by the expected delay between coalescence and the emission of EM radiation, which ranges from seconds for prompt gamma-rays to hours for afterglow emission. Sky localization using GWs is based primarily on triangulation across a light-travel baseline of 10–30 ms (Fairhurst 2011), but can be improved by including the signal amplitude and phase (Sidery et al. 2014; Singer et al. 2014). For the initial detector network, sky position for a merger event can be determined to several tens of square degrees (Fairhurst 2009).

Figure 2 shows an outline of the GW-EM coincidence pipeline. We begin with GW events found by a standard matched-filter analysis for coincident CBC signals across two or all three of the GW detectors (Abadie et al. 2012). The original results of this analysis identified a loud signal that had been injected into the data as a blind test of the pipeline (as seen in Figure 3), but did not find any other outliers that stood out above the coincident background of the GW instruments. The matched filter analysis gives some information about sky location (from timing) and distance (from amplitude), but better parameters are obtained through a coherent Bayesian follow-up (Veitch & Vecchio 2010; Veitch et al. 2014) which calculates the posterior probability of all possible physical signals (with varying component masses, location, distance, inclination, etc) using their coherent overlap with the data from all detectors. This is particularly useful for the relatively common case of a signal detected in H1 and L1, but too weak to be detected in Virgo, thus missing V1 timing information. The Bayesian analysis is still able to use the Virgo data to exclude regions of the sky where, for example, Virgo is particularly sensitive relative to the LIGO instruments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The sky location and distance estimates can be further refined by assuming that compact mergers arise from a known galaxy, for which we have reasonably complete catalogs out to the LIGO horizon (White et al. 2011). Each galaxy in the catalog can be assigned a host probability which is proportional to its blue-light luminosity (as a proxy for its mass, which we assume is proportional to the rate of compact binary mergers), as well as its overlap with the GW-derived posterior distributions in distance and sky-location. This reduces the search from several hundred square degrees to a handful of essentially point-source locations (given the angular resolution of ASM), at the risk of the true host being absent from the catalog.

The offline search for an EM counterpart from GBM or ASM is then triggered by the time of coalescence, and the GW-derived skyamp or the list of possible host galaxies with corresponding probabilities. We search over the parameter space of expected EM signals, and set nominal thresholds above which we consider GW-EM candidate coincident events.

As two compact objects orbit, they lose orbital energy and angular momentum to GWs. The orbit responds by shrinking and the orbital period decreases, causing even more rapid loss of orbital energy to gravitational radiation. This process ultimately leads in a runaway process to merger. At any moment during the inspiral, the GW emission is well characterized by high-order analytic post-Newtonian approximations to general-relativity. By following the orbit adiabadically through the detector bandwidth, a characteristic chirp waveform of increasing amplitude is produced that sweeps through the sensitive band (∼40–2000 Hz) of ground-based GW detectors.

LIGO-Virgo GW data from S6/VSR2+3 has been searched for CBC events with total system mass < 25 ${{M}_{\odot }}$, and the results have been reported in Abadie et al. (2012). The search made use of matched filtering, correlating the GW data in each detector against a bank of theoretical templates composed of model inspiral chirp signals for non-spinning systems with various component masses and orientations. Matched filtering describes the optimal linear filter for maximizing S/N in the presence of stationary Gaussian noise. For the filter corresponding to a frequency domain signal $\tilde{h}(f)$, the matched filter produces an expected single-detector S/N of,

$$\begin{eqnarray}&&{{\rho }^{2}}=\int _{0}^{\infty }\frac{4|\tilde{h}(f){{|}^{2}}}{{{S}_{n}}(f)}df\end{eqnarray} \tag{ 1 }$$

where ${{S}_{n}}(f)$ is the one-sided power spectral density of the stationary noise.

Due to the presence of non-Gaussian noise in the detector data, a ${{\chi }^{2}}$ consistency test (Allen 2005) is applied to events which are initially identified as local maxima in matched-filter S/N. A re-weighted S/N $\hat{\rho }$ is obtained by downgrading events using an empirically determined relationship that takes the level of inconsistency with the model signal into account (Babak et al. 2013),

$$\begin{eqnarray}\hat{\rho }=\left\{ \begin{array}{ccccccccccccccc} \frac{\rho }{{{[(1+{{(\chi _{r}^{2})}^{3}})/2]}^{1/6}}} & {\rm for}\;\;\chi _{r}^{2}\gt 1, \\ \rho & {\rm for}\;\;\chi _{r}^{2}\leqslant 1. \\ \end{array} \right.\end{eqnarray} \tag{ 2 }$$

Events from two or more detectors are searched for coincidence in time and mass parameters. Coincident events are ranked by the quadrature sum, ${{\rho }_{c}}$, of their re-weighted S/N. This combined S/N for each event is compared against an estimated background distribution derived from running the same coincidence analysis many times while applying unphysical relative time-shifts between the detector data (as in Figure 3). When the event has a ${{\rho }_{c}}$ value above a predetermined threshold corresponding to some low false-alarm rate (e.g., one false alarm from background per several 1000 yr), the event is considered a GW candidate. Background estimation is performed separately for various two-week epochs, detector combinations, and mass bins in order to isolate configurations with varying background levels.

The S6/VSR2+3 analysis identified one outlier event with an estimated false-alarm-rate of ∼1/7000 yr⁻¹. The event was ultimately revealed as a blind simulated signal injected in the data, which served as an end-to-end test of the detection strategy.⁸ Otherwise, no events stood out clearly from the expected background.

Bayesian integration is a natural way to determine whether or not a coherent GW signal is present in the data from a network of instruments, and if so, what are the likely parameters. While the coincident matched-filter analysis described in the previous section imposes some consistency on time and mass, a coherent analysis requires from the start a common physical gravitational waveform projected onto each instrument, making it particularly useful for a network of detectors.

If we let ${{\mathcal{H}}_{1}}$ represent the hypothesis that a CBC signal is present in the data, and let ${{\mathcal{H}}_{0}}$ represent only noise, the likelihood ratio provides the optimal statistic to distinguish the two,

$$\begin{eqnarray}&&{\Lambda }=\frac{P(d|{{\mathcal{H}}_{1}},I)}{P(d|{{\mathcal{H}}_{0}},I)}\end{eqnarray} \tag{ 3 }$$

where I represents any prior information we have about the system. Given that our signal hypothesis ${{\mathcal{H}}_{1}}$ represents a large population of possible signals with different waveforms determined by their set of binary parameters $\vec{\theta }$, proper evaluation of the likelihoood involves marginalizing over the parameter space,

$$\begin{eqnarray}&&P(d|{{\mathcal{H}}_{1}},I)=\int P(\vec{\theta }|{{\mathcal{H}}_{1}},I)P(d|{{\mathcal{H}}_{1}},\vec{\theta },I)d\vec{\theta }\end{eqnarray} \tag{ 4 }$$

where $P(\vec{\theta })$ represents the prior probability distribution of parameters $\vec{\theta }$ in our signal population—for example it may reflect that we expect systems to be uniformly distributed in volume with random orientations with respect to the detectors.

In addition to the question of whether or not a signal is present, we are also interested in determining the physical parameters of the binary system that are implied by a set of data. This is represented by the posterior probability distribution over the parameters $\vec{\theta }$,

$$\begin{eqnarray}&&P(\vec{\theta }|d,{{\mathcal{H}}_{1}},I)=\frac{P(\vec{\theta }|{{\mathcal{H}}_{1}},I)P(d|\vec{\theta },{{\mathcal{H}}_{1}},I)}{P(d|{{\mathcal{H}}_{1}},I)}.\end{eqnarray} \tag{ 5 }$$

If we are interested in the probability distribution only over a subset ${{\vec{\theta }}_{A}}$ of the parameters, where $\vec{\theta }\equiv \{{{\vec{\theta }}_{A}},{{\vec{\theta }}_{B}}\}$, we marginalize over those that remain,

$$\begin{eqnarray}&&P({{\vec{\theta }}_{A}}|d,{{\mathcal{H}}_{1}},I)=\int P(\vec{\theta }|d,{{\mathcal{H}}_{1}},I)d{{\vec{\theta }}_{B}}.\end{eqnarray} \tag{ 6 }$$

Nested sampling (Skilling 2004) is a computationally efficient alternative to MCMC techniques (Christensen & Meyer 2001), both of which are designed to intelligently sample and integrate the probability distribution of a high-dimensional parameter space of possible signals. For the case of CBC, we assume the spin of each neutron star is small and there are nine relevant physical parameters: two for the component masses, sky location, and orientation, as well as distance, time, and phase. We apply a Bayesian coherent follow-up based on nested sampling (Veitch & Vecchio 2010) to a large number of coincident events identified from the matched filter analysis. The output of the nested sampling routine is a set of sample vectors which trace the estimated 9-dimensional posterior probability distribution of system parameters (Equation (5)), as well as the integrated likelihood obtained by summing over them all (Equation (4)).

Extra-galactic events detectable by the initial LIGO-Virgo network are close enough so that a galaxy catalog can be used to identify probable hosts. While this is only marginally useful for a search using GW data alone due to the poor angular resolution of current GW detector networks, it becomes quite important for EM follow-ups where the EM sky resolution is significantly better (Nissanke et al. 2013; Hanna et al. 2014). Inclusion of a galaxy catalog for this offline EM search greatly reduces the computational cost and false-positive rate of the ASM analysis by reducing an area of ∼150 square degrees to just tens of individual points.

The Gravitational-wave Galaxy Catalog (GWGC; White et al. 2011) was designed specifically to aid in GW searches with the initial detector network. It contains distance, type, location, geometry, and blue-light luminosity information for ∼50 k galaxies within 100 Mpc, taken from the HyperLEDA database and other sources. To choose probable host galaxies from the catalog, we compare their distance and sky location against the probability distribution derived from the Bayesian analysis on GW data. Galaxies which occur in regions of high probability are considered for EM follow-up. At advanced LIGO/Virgo distances, galaxy catalogs suffer from incompleteness, with only the brighter galaxies making the flux limits of surveys, and deep surveys often covering only fractions of the sky. Hanna et al. (2014) outline some of the considerations necessary to effectively use incomplete catalogs in that regime. For this initial LIGO/Virgo study, we assume the GWGC is 100% complete out to the distance of detectable NS/NS mergers.

The GW posterior probability distribution must be estimated from the discrete sampling provided by the nested sampling procedure. We estimate the density of posterior samples at each galaxy location and distance using Gaussian kernel density estimation (KDE). We treat distance and sky location independently due to limited sampling of the posterior, typically thousands of points.

For a galaxy with blue-light luminosity ${{L}_{10}}({{10}^{10}}\;{{L}_{\odot }})$ at distance r < 100 Mpc and location ${\Omega }$, the estimated distribution from the N posterior samples is

$$\begin{eqnarray}&&{{P}_{{\rm GW}}}(r)=\frac{1}{N}\mathop{\sum }\limits_{i}^{N}\frac{1}{{{r}^{2}}}w({{r}_{i}}-r,{{\sigma }_{r}}=10-20\%)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{P}_{{\rm GW}}}({\Omega })=\frac{1}{N}\mathop{\sum }\limits_{i}^{N}w({{\vec{x}}_{i}}-\vec{x},{{\sigma }_{{\Omega }}}=3{}^\circ )\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&P({{L}_{10}},r,{\Omega })={{L}_{10}}{{P}_{{\rm GW}}}(r){{P}_{{\rm GW}}}({\Omega })\end{eqnarray} \tag{ 9 }$$

where $w(\mu ,\sigma )$ represents the 1D and 2D Gaussian kernels, normalized to give P = 1 in the case of a uniform posterior. Figure 4 shows an example of the kernel density estimation applied to the distance parameter, as well as a sample of matching galaxies, weighted by all three factors in Equation (9), against the GW skymap. The distance KDE bandwidth is chosen according to the estimated individual galaxy distance errors, which are around 10–20% depending on the original source catalog. The 3° sky bandwidth is chosen to balance statistical errors (resulting Poisson errors are at the level of a few %) with resolution of the GW network. Finally, we have undone the r² distance prior from the assumption of homogeneously distributed sources during the Bayesian integration as the galaxies themselves already include this volume effect.

Zoom In Zoom Out Reset image size

The GBM (Meegan et al. 2009) on board the Fermi spacecraft measures photon rates from 8 keV–40 MeV. The instrument consists of 12 semi-directional NaI scintillation detectors and 2 BGO scintillation detectors which cover the entire sky not occluded by the Earth (about 65%). The lower-energy NaI detectors have an approximately ${\rm cos} \theta$ response relative to angle of incidence, and relative rates across detectors are used to reconstruct the source location to a few degrees. The BGO detectors are much less directional, and can be used to detect and resolve the higher energy spectrum above ∼200 keV.

GBM produces on-board triggers for GRB events by looking for multi-detector rate excess over background across various energy bands and timescales. In the case of a trigger, individual photon information is sent to the ground and the event is publicly reported. Those events which have been confirmed as GRBs have already been studied in coincidence with LIGO-Virgo data (Abadie et al. 2012). In addition to the triggered events, survey data is available which records binned photon counts over all time. For this offline analysis of short transients, we use the CTIME (time-resolved) daily data, which is binned at 0.256 s over 8 energy channels for each detector. A new GBM data product for continuous time-tagged events was introduced in 2012 that provides complete individual photon information with 2 μs and 128 energy channel resolution. This data product was not available for un-triggered times during the initial LIGO-Virgo science runs, but will enhance offline sensitivity to short bursts when the advanced GW detectors begin operation.

In this section, we develop a procedure to coherently search GBM detector data for modeled events (Blackburn et al. 2013). The idea is that by processing multiple detector data coherently, we can obtain a greater sensitivity than when considering one detector at a time. Greater computational resources available offline (versus on-board) also allow for more careful background estimation to be done. For this analysis, we can relax to some extent the strict 2-detector coincidence requirement used to veto spurious events on board as the GW trigger means much less time and sky area is considered. These advantages help to balance out the coarse time resolution of CTIME data, which reduces offline sensitivity to very short bursts prior to 2013.

Each detector is subject to a substantial time-varying background from bright high-energy sources that come in and out of the wide field of view (FOV), as well as location-dependent particle and Earth atmospheric effects. This background must be estimated and subtracted out to look for any prompt excess. In this analysis where we are interested in the background estimate for a short foreground interval $[-T/2,+T/2]$ where $T{\mkern 1mu} \sim$1 s, we estimate the background using a polynomial fit to local data from [$-10\;T,+10\;T$] (minimum $\pm 5\;{\rm s}$), excluding time [$-3\;T/2,+5\;T/2$] around the foreground interval to avoid bias from an on-source excess. The polynomial degree is determined by the interval length to account for more complicated background variability over longer intervals. It ranges from 2 (minimum) to $1+0.5{{{\rm log} }_{2}}T$. A separate fit is done for each detector/channel combination, and determines the background rate estimation over the foreground interval, as well as its systematic uncertainty from fitting error. Data from channels with poor fits (large ${{\chi }^{2}}$) are excluded from the analysis.

High-energy cosmic rays striking a NaI crystal can result in long-lived phosphorescent light emission. The detector may interpret this is a rapid series of events, creating a short-lived jump in rates for one or multiple channels, and severely distorting the background fit if not accounted for. They are identified with a simple procedure that scans for rapid 1-bin spikes in the background interval. The affected bins are removed from the background fitting. Cosmic rays affecting the foreground interval are handled differently, as described in Section 5.1.

A likelihood ratio combines information about sources and noise into a single variable. It is defined as the probability of measuring the observed data, d, in the presence of a particular true signal H₁ (with some source amplitude s > 0) divided by the probability of measuring the observed data in noise alone ${{H}_{0}}(s=0)$,

$$\begin{eqnarray}&&{\Lambda }(d)=\frac{P(d|{{H}_{1}})}{P(d|{{H}_{0}})}.\end{eqnarray} \tag{ 10 }$$

When signal parameters such as light curve, spectrum, amplitude s, and sky-location $\alpha ,\delta$ are unknown, one can either marginalize over the unknown parameters, or take the maximum likelihood over the range to obtain best-fit values.

For binned, uncorrelated Gaussian noise,

$$\begin{eqnarray}&&P({{d}_{i}}|{{H}_{1}})=\mathop{\prod }\limits_{i}\frac{1}{\sqrt{2\pi }{{\sigma }_{{{d}_{i}}}}}{\rm exp} (-\frac{{{(\tilde{{{d}_{i}}}-{{r}_{i}}s)}^{2}}}{2\sigma _{{{d}_{i}}}^{2}})\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&P({{d}_{i}}|{{H}_{0}})=\mathop{\prod }\limits_{i}\frac{1}{\sqrt{2\pi }{{\sigma }_{{{n}_{i}}}}}{\rm exp} (-\frac{{{\tilde{{{d}_{i}}}}^{2}}}{2\sigma _{{{n}_{i}}}^{2}})\end{eqnarray} \tag{ 12 }$$

where we have used $\tilde{{{d}_{i}}}={{d}_{i}}-\langle {{n}_{i}}\rangle$ to represent the background-subtracted measurements in each detector-time-energy bin, ${{\sigma }_{{{n}_{i}}}}\;\;{\rm and}\;\;{{\sigma }_{{{d}_{i}}}}$ for the standard deviation of the background and expected data (background+signal), r_i for the location/spectrum-dependent instrumental response, and s for the intrinsic source amplitude at the Earth. Maximizing the likelihood ratio is the same as maximizing the log-likelihood ratio $\mathcal{L}={\rm ln} {\Lambda }$,

$$\begin{eqnarray}&&\mathcal{L}=\mathop{\sum }\limits_{i}\left[ {\rm ln} \frac{{{\sigma }_{{{n}_{i}}}}}{{{\sigma }_{{{d}_{i}}}}}+\frac{{{\tilde{{{d}_{i}}}}^{2}}}{2\sigma _{{{n}_{i}}}^{2}}-\frac{{{(\tilde{{{d}_{i}}}-{{r}_{i}}s)}^{2}}}{2\sigma _{{{d}_{i}}}^{2}} \right].\end{eqnarray} \tag{ 13 }$$

The dependence of the response factors r_i on sky location is complicated, so the likelihood ratio is calculated over a sample grid of all possible locations. Assuming a single location, the remaining free parameter is the source amplitude s. The variance in the background-subtracted detector data includes both background and source contributions,

$$\begin{eqnarray}&&\sigma _{{{d}_{i}}}^{2}=\sigma _{{{n}_{i}}}^{2}+{{r}_{i}}s+\sigma _{{{r}_{i}}}^{2}\;{{s}^{2}}\quad (s\geqslant 0)\end{eqnarray} \tag{ 14 }$$

with $\sigma _{{{r}_{i}}}^{2}$ representing Gaussian-modeled systematic uncertainty in the instrumental response. Source terms are only included for physical s ≥ 0, else their contribution is zero. The background contributes Poisson error, as well as any systematic variance $\sigma _{{{b}_{i}}}^{2}$ from poor background fitting which is also assumed to be Gaussian,

$$\begin{eqnarray}&&\sigma _{{{n}_{i}}}^{2}=\langle {{n}_{i}}\rangle +\sigma _{{{b}_{i}}}^{2}.\end{eqnarray} \tag{ 15 }$$

We find ${{s}_{{\rm best}}}$ which maximizes $\mathcal{L}$ by setting the derivative $d\mathcal{L}/ds$ to zero using Newton's method.

To consider all possible source amplitudes we need to integrate the likelihood P(d ∣ s) (Equation (11)) over a prior on s,

$$\begin{eqnarray}&&P(d)=\int P(d|s)P(s)ds.\end{eqnarray} \tag{ 16 }$$

For a given set of detector data d, the likelihood $P(d|s)\;{\rm over}\;s$ is almost the product of individual Gaussian distributions (not quite Gaussian because ${{\sigma }_{d}}$ depends on s). The product of Gaussian distributions with mean values ${{\mu }_{i}}$ and standard deviations ${{\sigma }_{i}}$ is itself Gaussian with mean and variance,

$$\begin{eqnarray}&&{{\mu }_{{\rm prod}}}=\frac{\sum {{\mu }_{i}}/\sigma _{i}^{2}}{\sum 1/\sigma _{i}^{2}},\quad \sigma _{{\rm prod}}^{2}=\frac{1}{\sum 1/\sigma _{i}^{2}}.\end{eqnarray} \tag{ 17 }$$

In this case, ${{\mu }_{i}}=\tilde{{{d}_{i}}}/{{r}_{i}}\;{\rm and}\;{{\sigma }_{i}}={{\sigma }_{{{d}_{i}}}}/{{r}_{i}}$. We estimate the variance of $\mathcal{L}$ over s as,

$$\begin{eqnarray}&&\sigma _{\mathcal{L}}^{2}=\frac{1}{\sum r_{i}^{2}/\sigma _{{{d}_{i}}}^{2}},\quad \sigma _{{{d}_{i}}}^{2}\;{\rm evaluated}\;{\rm at}\;{{s}_{{\rm best}}}\end{eqnarray} \tag{ 18 }$$

and choose a scale-free prior with fixed $\beta =1$, so that our choice of form for an amplitude prior at the Earth does not translate back into a luminosity distribution that varies with distance.

One difficulty with any power-law prior is that it diverges for $s\to 0$. We enforce a finite and well-behaved prior by multiplying by a prefactor,

$$\begin{eqnarray}&&P(s)=\left[ 1-{{e}^{-{{(s/\gamma {{\sigma }_{\mathcal{L}}})}^{\beta }}}} \right]{{s}^{-\beta }}\end{eqnarray} \tag{ 19 }$$

so that $P(s)$ reaches a maximum constant value of ${{(1/\gamma {{\sigma }_{\mathcal{L}}})}^{\beta }}$ for small s. The tunable parameter $\gamma$ sets the number of standard deviations at which the prior begins to plateau, and we use $\gamma =2.5$. We then approximate P(s) as constant over a range of ${{\sigma }_{\mathcal{L}}}$ for any $s\gt 0$, and include a correction to account for clipping of the Gaussian for non-physical $s\lt 0$, which can be represented by the error function. The final approximation for the amplitude-marginalized log-likelihood becomes,

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{L}(d) & = & {\rm ln} {{\sigma }_{\mathcal{L}}}+{\rm ln} \left[ 1+{\rm Erf}(\frac{{{s}_{{\rm best}}}}{\sqrt{2}{{\sigma }_{\mathcal{L}}}}) \right]+\mathcal{L}(d|{{s}_{{\rm best}}}) \\ {} & {} & +\left\{ \begin{array}{ccccccccccccccc} {\rm ln} \left[ 1-{{e}^{-{{({{s}_{{\rm best}}}/\gamma {{\sigma }_{\mathcal{L}}})}^{\beta }}}} \right]-\beta {\rm ln} {{s}_{{\rm best}}} & {{s}_{{\rm best}}}\gt 0 \\ -\beta {\rm ln} (\gamma {{\sigma }_{\mathcal{L}}}) & {{s}_{{\rm best}}}\leqslant 0 \\ \end{array} \right. \\ \end{array}\end{eqnarray} \tag{ 20 }$$

which contains factors from the Gaussian width, fractional overlap with $s\gt 0$, maximum likelihood at ${{s}_{{\rm best}}}$, and scaling from P(s), respectively. Finally we are free to calibrate the log-likelihood by subtracting the expected $\mathcal{L}(d)$ calculated for no signal at a reference sensitivity: ${{\mathcal{L}}_{{\rm ref}}}=-\beta {\rm ln} \gamma +(1-\beta ){\rm ln} {{\sigma }_{{\rm ref}}}.{{\sigma }_{\mathcal{L}}}$ represents the source amplitude required for a $1\sigma$ excess in the combined data, and is around 0.05 photons/s/cm$^{2}\times {{(T/1{\rm s})}^{-1/2}}$ [50–300 keV] for a typical source spectrum and background level. Figure 5 shows the coherent S/N expected from all detectors for a 0.512 s-long event with normal GRB spectrum and constant amplitude of 1.0 photons/s/cm$^{2}$, and compares it to the S/N expected from individual detectors alone in the 50–300 keV band.

Zoom In Zoom Out Reset image size

A GW trigger provides an accurate time of coalescence t_c as well as an approximate sky location to tens or hundreds of square degrees. This matches well the time and sky location constraints provided by a prompt observation by GBM, so that coincidence between the two instruments is particularly effective. To make use of the coherent analysis described in the previous section, we must first obtain an expected all-sky GBM response for a given counterpart. We make use of three representative source spectra using a set of precomputed response tables. The representative spectra are Band functions with the parameters listed in Table 2. For the soft spectrum, ${{E}_{{\rm peak}}}\;{\rm and}\;\beta$ are undefined. The tables provide channel-dependent expected counts for the direct and spacecraft-scattered response as a function of source location to ∼1°resolution. Contributions from atmospheric scattering are available as counts summed over the two CTIME channels covering 50–300 keV as a function of source and Earth relative position for the two most common spacecraft rocking angles.

Table 2.  Model Band Spectral Parameters Used to generate GBM Detector Response Tables Used in the Coherent Analysis

Spectrum	${{E}_{{\rm peak}}}$	α	β
normal	230 keV	1	2.3
hard	1 MeV	0	1.5
soft	⋯	2	⋯

Note. The soft and hard spectra were originally designed to produce a representative response over the primary GBM-GRB localization band 50–300 keV, and a number of more realistic response models are currently being developed which will be more approriate for full-spectral all-sky analysis. The construction of these all-sky, all-time models is difficult computationally due to the many degrees of freedom in the detector response.

Download table as:  ASCIITypeset image

For each GW candidate event, we search a window [−30, 30 s] relative to t_c for prompt excess between 0.256 and 8 s long. Emission outside of the standard accretion timescale covering the first seconds following merger is speculative. However, any events detectable in GWs will be much closer than known sGRBs to date, and therefore they provide an opportunity to search below threshold for weak and possibly less-beamed precursor or extended emission. To appropriately tile the search in time and duration T, we use rectangular windows with T spaced by powers of two (0.256, 0.512, 1 s, etc). Their central times are sampled along the search interval in units of T/4 to provide an even mismatch in S/N across search windows. For each emission model tested, the likelihood ratio is then marginalized over all windows and spectra.

The GW data also provides a rough sky-location, which can be represented as a prior probability distribution over the sky ${{P}_{{\rm GW}}}({\Omega })$. This GW prior is multiplied by the GBM likelihood ratio, which is also a function of sky position, before marginalization over sky location is done,

$$\begin{eqnarray}&&{{{\Lambda }}_{{\rm GW}-{\rm GBM}}}=\int d{\Omega }{{P}_{{\rm GW}}}({\Omega }){{{\Lambda }}_{{\rm GBM}}}({\Omega })\end{eqnarray} \tag{ 21 }$$

where the GBM likelihood ratio depends on source spectrum, on-source time window, and location. While the sky prior can in principle be improved by incorporating knowledge of the anisotropic local mass distribution (using a galaxy catalog), we note that as neither the GW network nor GBM have precise (sub-degree) sky localization, the gains from sharpening the prior in this way are small, and catalog incompleteness and luminosity-rate relationships introduce unnecessary complications. The situation is different when considering better-localizated X-ray (and optical/ratio) counterparts.

The background for GW-GBM analysis is characterized by the probability of observing a given ${{{\Lambda }}_{{\rm GW}-{\rm EM}}}$ from random coincidence between GW and GBM data. We calculate the expected distribution of ${{{\Lambda }}_{{\rm GW}-{\rm EM}}}$ by running the coincidence search defined by Equation (21) on GW background events derived from time-shift analysis. The time-shift sample represents 100 artificial time-shifts (∼seconds) applied to data between the LIGO Hanford and Livingston sites, as well as a further 10 time-shifts ($\sim$day) applied between the GW and GBM data. Thus we obtain a background sample representing characteristic noise from $\times 1000$ the original foreground live-time.

GRB simulated signals are injected into the data by using the instrument response model to predict the signal contribution from a source with a given spectrum and fluence. In our analysis, we simulate a standard candle signal following the normal-type Band spectrum with a fluence at the Earth of 1 photon/cm²/s in 50–300 keV assuming a source at 30 Mpc. The signal assumes normal GRB spectral parameters according to Table 2 and lasts for 1 s, corresponding to an isotropically emitted energy of about 6.6 × 10⁴⁶ erg. When simulating a signal at distance D, the fluence is reduced by a factor of ${{(30\;{\rm Mpc}/D)}^{2}}$.

The ASM (Levine et al. 1996) on the RXTE surveyed the X-ray sky between 1.5–12 keV from 1996–2011. The instrument consisted of three long-and-narrow shadow-mask X-ray cameras rigidly attached to a common motorized axis of rotation. Each camera saw a 6 by 90 degree FOV, and together covered about 3% of the sky at any one time. As the cameras scanned across the sky while tiling 90 s dwells, they were able to localize a steady source to about 01. The duty cycle of the ASM along with constraints from the observing schedule of other RXTE instruments resulted in a randomly chosen location being scanned by an ASM camera a few times per day.

An ASM camera observation consists of a superposition of the shadow pattern of all X-ray sources within the FOV. The data can be modeled as the contributions from active sources with known locations along with a diffuse background. Amplitudes for sources are treated as free parameters, and estimated using a linear least-squares fit to the data (Levine et al. 1996), with errors dominated by photon counting statistics. The amplitude from any additional position of interest within the FOV can also be estimated by adding the location as a source during the fit. For a weak source, the error on estimated amplitude is dominated by the diffuse background which typically contributes a 3-σ error of ∼20 mCrab (4.8 × 10⁻¹⁰ erg/cm²/s between 2–10 keV).

The ASM measurements consist of irregularly spaced flux observations with varying uncertainties from specific regions in the sky, which makes a search strategy for an arbitrary flux excess complicated. We use limited knowledge of the sGRB X-ray afterglow signal to narrow the search to signals that decay like broken power laws in flux beginning from the initial time-of-merger t₀ as determined by the GW data.

Following Zhang et al. (2006), we model a canonical sGRB X-ray afterglow by a double-broken power law with a short region of rapid decline (power law index $\alpha \approx 3$) lasting 1e2–1e3 seconds, followed by a standard afterglow decay ($\alpha \approx -1.2$), and finally a break and raid decay ($\alpha \approx -2$) after 10⁴ to 10⁵ seconds. Additionally there may be a plateau or extended emission somewhere after the rapid decay phase around 10² to 10³ seconds. We attempt to create a generic light-curve with minimal free parameters which covers possible observed light curve scenarios, in context of the limited sampling of ASM data, and also allows for the possibility of delayed or non-standard X-ray emission due to off-axis observation or some other unobserved phenomenon.

For this we sample double-broken model power-law light curves with break times sampled on a grid (Figure 6) and freely varying power-law indices. We label N local measurements of flux excess d_i by ASM at times t_i with standard deviations ${{\sigma }_{i}}$. A given light curve parameterization defines expected relative counts at the measurement times ${{L}_{i}}=L({{t}_{i}})$. We can then define the weighted sum of measurements c_i which maximizes the S/N of the summed data,

$$\begin{eqnarray}&&D=\mathop{\sum }\limits_{i}\frac{{{L}_{i}}{{d}_{i}}}{\sigma _{i}^{2}}.\end{eqnarray} \tag{ 22 }$$

The measurements d_i are already assumed to be zero-mean, giving a measured S/N for the coherent sum D of,

$$\begin{eqnarray}&&\rho =\frac{D}{\sqrt{{{\sum }_{i}}L_{i}^{2}/\sigma _{i}^{2}}}.\end{eqnarray} \tag{ 23 }$$

The expected S/N of each measurement is added in quadrature for the total expected S/N of the coherent sum,

$$\begin{eqnarray}&&\langle \rho \rangle =\sqrt{\mathop{\sum }\limits_{i}\frac{L_{i}^{2}}{\sigma _{i}^{2}}}.\end{eqnarray} \tag{ 24 }$$

If we are fortunate to have an observation soon after the beginning of the afterglow, the first measurement is likely to dominate the sum due to the rapid decay. If only late-time measurements are available, the coherent sum naturally bins the data in order to gain in sensitivity. A positive power-law index for the first segment of the light curve can handle situations where the afterglow does not begin immediately after the burst, for example due to beaming effects.

Zoom In Zoom Out Reset image size

For each pair of power-law breakpoints (${{t}_{1}},{{t}_{2}}$), we fit the three power-law indices (${{\alpha }_{1}},{{\alpha }_{2}},{{\alpha }_{3}}$) using a minimization routine under the constraints $-3\lt {{\alpha }_{1}}\lt 2,-2\lt {{\alpha }_{2}}\lt 1$, and $-3\lt {{\alpha }_{3}}\lt 0$. This allows for a wide variety of afterglow waveforms including the standard canonical double-broken power law, an initially rising light curve which then decays normally, or an isolated burst of extended emission that might occur minutes to days after merger. For each set of ASM measurements from a location of interest [0, 4 d] about t₀, we record the fit parameters which give the maximum S/N, as well as a variety of auxiliary parameters useful for characterizing non-Gaussian noise.

Table 1.  High-energy Photon Survey Instruments Used to Search for EM Counterparts

Instrument	Energy	FOV	Resolution	Cadence
Fermi GBM	8 keV–40 MeV	65%	>5°	1.5 h
RXTE ASM	1 keV–10 keV	3%	01	1.5 h

Note. Field of View (FOV) is represented as a percentage of the entire sky.

Download table as:  ASCIITypeset image

The ASM point source flux estimation is well-suited to follow up locations of likely host galaxies for GW events. The shadow-mask flux reconstruction isolates contributions to about 01 on the sky. This is small enough that a large fraction of the sky can be excluded by requiring known galaxy coincidence, but also large enough so that all but the most extended nearby galaxies can be targeted by a single location. For each host location, we generate ASM flux measurements from available camera dwells by adding that single location as a test source along with other known active X-ray sources during the fit of each dwell's shadow-mask data. Derived flux measurements from that location are then searched numerically for the five parameter light curve (two breakpoints and three decay factors) which maximize S/N of the coherent sum. Because the data are sparse, it is difficult to maximize S/N numerically for freely varying breakpoints, so they are sampled on a fixed grid and only the three power-law indices are searched with standard numerical minimization routines as described above.

We take the maximum S/N $\rho$ observed across all model light curves to weigh each candidate host galaxy according to likelihood. The conversion from $\rho$ to likelihood is determined by a global S/N distribution of X-ray background derived from the same methods. The distribution is fit to a function $f(\rho )$ parameterized by a Gaussian (maximum value 1) with power-law tail (Figure 7), and this empirical fit is used to calculate the likelihood=$1/f(\rho )$ at each S/N measurement. For a given GW event with a set of possible hosts, we obtain a final rank for the coincident observation,

$$\begin{eqnarray}&&{{{\Lambda }}_{{\rm GW}-{\rm ASM}}}=\mathop{\sum }\limits_{{\rm host}}\left[ {{P}_{{\rm GW}}}({\rm host})\times f{{({{\rho }_{{\rm host}}})}^{-1}} \right]\end{eqnarray} \tag{ 25 }$$

where the likelihood ${{{\Lambda }}_{{\rm GW}-{\rm ASM}}}$ represents the probability of an afterglow signal in the ASM data, marginalized over all possible hosts. The hosts must maintain a ${{P}_{{\rm GW}}}\gt 1$ (a threshold representing the case for a galaxy of typical $\sim$ MW mass and no location or distance prior information), and at maximum 200 hosts are scanned to limit computational expense. ${{P}_{{\rm GW}}}({\rm host})$ is derived, as in Equation (9), from the sky and distance overlap of each host with the prior distributions determined by the GW data and host mass (luminosity). To convert the likelihood into a more S/N-like quantity, we define rank $r=\sqrt{2{{{\Lambda }}_{{\rm GW}-{\rm ASM}}}}$.

Zoom In Zoom Out Reset image size

We simulate a family of X-ray afterglow signals by adding model light curves directly to the derived ASM flux measurements for a given location. We assume that the error estimate for each flux measurement remains dominated by contributions from the diffuse background, so that Poisson or other systematic error from the source can be neglected. The model light curve is then injected into the data from a reconstructed location at a time of interest (Figure 8). For a given light curve shape, we recover the distribution of maximum S/N measurements from the template search as a function of injected amplitude, or equivalently injected distance for a standard candle source.

Zoom In Zoom Out Reset image size

For the case of a GW-ASM coincidence, we begin with a set of simulated NS/NS coalescence events detected with the GW analysis procedure. For S6/VSR2+3, a large number of simulations were done at random sky locations and distances in order to evaluate the overall detection efficiency of the pipeline (Abadie et al. 2012). To characterize the sensitivity of the joint GW-ASM coincidence analysis which assumes signals originating from a set of host galaxies, we first create an artificial galaxy at the location and distance of the GW simulation, with a luminosity taken from a (luminosity-weighted) random draw of galaxies from the catalog within some compatible volume. We then synthesize ASM data from this fake location and distance by appropriating the ASM flux measurements from the randomly chosen galaxy. Reconstructed ASM flux measurements from the remaining true galaxies within the GW-derived sky region are included in the coincidence analysis as before, as well as the fake galaxy with the simulated lightcurve added to the ASM data. We then obtain a distribution of the joint GW-ASM detection statistic for a given CBC system as a function of distance and intrinsic EM amplitude.

We test performance of the method against a standard X-ray afterglow typical of Swift XRT observations of those from short GRBs (Figure 6). This light curve has power-law decay indices of −2, −1.2, and −2 between breakpoints at $3\times {{10}^{2}}\;\;{\rm and}\;\;5\times {{10}^{4}}\;{\rm s}$, and a 1–10 keV luminosity of 10⁴⁶ erg s⁻¹ at the first break-point. We also assume a Crab-like spectrum, for convenient conversion into ASM band counts which are Crab calibrated.

In addition, we inject a representative sample of X-ray burst signals which occur minutes to days after the merger, motivated by theoretical magnetar wind-driven scenarios (Gao et al. 2013; Zhang 2013). We test a short, medium, and late-time burst beginning with a rise $\propto \;{{t}^{2}}\;{\rm until}\;{{t}_{1}}={{10}^{2}}$, 10³, and 10⁴ s, decaying as ${{t}^{-1}}\;{\rm until}\;{{t}_{2}}=2{{t}_{1}}$, and then falling as t⁻² afterward. The simple wind-driven emission scenarios consider a total amount of energy released as the magnetar spins down, and then either dissipating via internal or external shocks, releasing some fraction of the dissipated energy in X-rays. We set the total energy released in X-rays over the entire period of the burst to be 10⁴⁹ erg in 2–10 keV, resulting in bursts with instantaneous luminosities slightly greater than that for the standard afterglow at peak. In this way, they can be thought of representing the extended emission/flares seen in some sGRB afterglows. As in the standard afterglow model, we assume a Crab-like spectrum for convenient translation into ASM counts.

The GW-ASM background distribution is found by running the coincidence search using GW background events derived from time-shift analysis, similar to that for GBM. The same 10 additional time-shifts (spaced in multiples of 4 days) are applied between the GW and ASM data to increase statistics. We then obtain the distribution of the GW-ASM rank statistic expected from accidental coincidence between GW and ASM background. This distribution, as well as that for the family of simulated afterglow signals, are shown in Figure 9.

Zoom In Zoom Out Reset image size

We have tested the end-to-end GW-EM pipelines for both GBM and ASM on a selection of real GW background noise events from the LIGO-Virgo S6/VSR3 joint science run, as well as simulated GW-EM events injected in coincidence in both the GW and EM data streams. By comparing the properties of GW background and simulations, we identify a number of selection cuts that help reduce the non-Gaussian background present in the EM data. The ASM selection cut is based on the reduced ${{\chi }^{2}}$ of the ASM data after the best-fit double-broken power-law light curve has been subtracted (Figure 10). However, due to the sporadic coverage of the ASM measurements, it is difficult to select on any morphological properties of the light curves themselves, beyond the weak requirement that the measurements be consistent with one or more of the possible templates.

Zoom In Zoom Out Reset image size

We place a nominal threshold of 18 on ${{{\Lambda }}_{{\rm GW}-{\rm ASM}}}$, the GW-ASM coincident likelihood ratio. Because the likelihood is largely determined by the Gaussian-like shape of the ASM background distribution, we represent this as a threshold of 6 on the alternative rank parameter $r=\sqrt{2{\Lambda }}$. The distribution of r for the GW background events, as well as a variety of simulations is shown in Figure 9. With the ASM coincidence requirement applied, the loudest GW background event drops from a combined S/N of 11.0 to 9.13 (Figure 14), a reduction of ∼17%. ASM is generally sensitive to afterglow signals within the LIGO horizon, however detectability varies widely since it is highly dependent on exactly when, or if ever, the ASM observation(s) occur relative to the peak in the light curves (Figure 11). Similarly, detections may fail if the GW localization is not sufficient to pick out the correct galaxy host among the top 200 candidates, which happens about 7% of the time at 10 Mpc, and 32% of the time at 30 Mpc.

Zoom In Zoom Out Reset image size

In the GBM detector, most short background events are due to phosphorescent events in the NaI detectors, which can arise from cosmic rays. While we implemented a rough cut to exclude such events from contaminating the background fits during the analysis, we implement a further selection to remove remaining particle events from the foreground interval. The particle events are unique in that they show up in a single detector, and are generally soft in reconstructed energy. Thus we implement two cuts based on the ratio of detected channel 0 S/N in the loudest and second-loudest NaI detectors, as well as between channel 0 and channel 1 in the loudest detector. With these cuts we are able to remove most of the remaining soft particle background.

We also implement a sky-coincidence cut between the GW sky location and GBM by comparing the GBM likelihood ratio with and without the GW sky prior folded in. Since the sky prior typically covers a fractionally small region of the sky for a well-localized event (100 s of square degrees), we expect an appropriately normalized coincident observation to have a correspondingly higher likelihood ratio than one that assumes an isotropic all-sky prior. We expect a factor that can be roughly in line with the fractional reduction in sky area (∼400), but since this can be moderated by many effects (averaging windows, systematic errors in reconstructed location), we choose an empirical factor of e² as a coincidence requirement. This cut is still able to reject a large fraction of the background where we do not expect any systematic increase in likelihood ratio due to incorporation of a GW sky prior (Figure 12).

Zoom In Zoom Out Reset image size

Finally we place a nominal threshold of 30 on the coincident GBM likelihood ratio, which corresponds roughly to that of the weakest un-triggered GRBs detectable by this method (Figure 13). After this threshold and the previously mentioned selection cuts, there is only one coincident GW-GBM event remaining, which has a combined S/N ${{\rho }_{c}}$ of 8.1 (Figure 14), ∼26% lower than the previous loudest event.

Zoom In Zoom Out Reset image size

The cumulative histograms of S/N of the time-shift GW background events before and after EM coincidence are shown in Figure 14. This can be qualitatively compared to the published background distribution shown in Figure 3.

The original low-threshold GW events used in this analysis were subject to a single detector S/N threshold of 5.5. This means that for a double-coincident H1L1 event, the absolute minimum network S/N is $\sqrt{2\times {{5.5}^{2}}}=7.78$, before any reduction by the ${{\chi }^{2}}$ consistency check (Equation (2)). Near this threshold, one must be careful when assuming simple range scaling relationships between combined S/N thresholds and effective range of the analysis because of the upstream cuts.

Past the knee in the original GW background distribution ${{\rho }_{c}}\gtrsim 9$, the inverse relationship between distance and combined S/N should be robust, so we can estimate a factor of ≳15% increase in range for this sample at equivalent false-alarm rate given events where EM coincidence can provide at least a factor of 1/1000 background rejection. While the background sample used in this study reliably probes about three orders of magnitude in rate, the extended background studies from Abadie et al. (2012) and represented in Figure 3 show that for H1L1 events near the end of the LIGO S6 with chirp mass $3.48\leqslant {{\mathcal{M}}_{c}}\lt 7.40$, the background distribution scales roughly as ∼${{100}^{-{{\rho }_{c}}}}$ in combined S/N over eight orders of magnitude in rate. For the special case in which both a prompt gamma-ray and X-ray afterglow counterpart are observable, the accidental coincidence factors measured here give a large joint rate of background rejection around ∼10⁷, which implies a factor of ∼1.5 increased range to such events if maintaining a fixed false-alarm-rate given an original threshold of ${{\rho }_{c}}\gt 10$ (Camp et al. 2013) before EM coincidence.

The efficiency of the EM pipelines to our standard-candle signal injections is shown in Figure 15. The efficiency fractions are calculated from the end-to-end simulation of a joint NS/NS inspiral signal with corresponding prompt and afterglow EM counterpart, and the process is triggered by a low-threshold (${{\rho }_{c}}\gtrsim 8.5$) GW candidate. The presence of a GW trigger is required, and the represented fraction does not include efficiency factors from the analysis of GW data itself. However, the EM follow-up efficiency is still generally influenced by the quality of the GW sky localization.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

GBM views the entire unocculted sky (65%) when not in the South Atlantic Anomaly (∼15%), and this duty-cycle dominates the efficiency factor out to 40 Mpc. This is not surprising as our injeciton amplitude was chosen to be moderately detectable at 30 Mpc (and still several hundred times weaker than a typical sGRB).

The chance of detecting an X-ray afterglow signal with ASM is much more variable (Figure 11) due primarily to the large variability in delay between onset of the afterglow and the first available measurement. The ASM follow-up is also more sensitive to the sky localization accuracy from the GW trigger due to the choice to follow-up only the most probable 200 individual galaxy host locations. Increased distance both increases the sky area uncertainty (due to decreased GW S/N), and increases the area-density of galaxies on the sky. We do not observe ASM counterparts above threshold beyond ∼30 Mpc. For the X-ray burst model waveforms, longer duration bursts (at equivalent total fluence) were relatively easier to detect as they better matched the ASM cadence.

In this paper, we introduced a strategy to search for high-energy EM counterparts to GW binary colescense events in archival satellite survey data. We designed an end-to-end GW candidate follow-up pipeline and tested it on a large number of background binary colescence GW events from time-shifted initial LIGO/Virgo data during their most recent science runs (S6/VSR2+3). The representative noise events were subject to a fully coherent Bayesian parameter estimation analysis in order to obtain sky and distance posterior probability distributions. These distributions were used to obtain a set of probable hosts from a catalog of nearby galaxies.

Two custom follow-up methods were designed to search for both a prompt gamma-ray counterpart in offline data from the Fermi GBM within ±30 s of the GW coalesence time and consistent with the GW sky location, as well as a generic family of X-ray afterglow lightcurves in data from the RXTE ASM at the locations of possible host galaxies, and parametrized by a generic broken power law with characteristic timescales of minutes to days. The requirement of a GBM or ASM coincident counterpart reduced the number of background events by factors of ${{10}^{-4}}\;{\rm and}\;{{10}^{-3}}$, respectively, reducing the GW ampltiude of the loudest suriving background event by ∼15–20%.

Both EM pipelines were tested on a set of joint GW-EM simulated signals, where GW simulations corresponded to NS/NS mergers at random times and sky locations/orientations. The GBM follow-up pipeline maintained sensitivity to weak coincident prompt bursts with a normal GRB specturm and a fluence corresponding to an isotropic-equivalent energy release of ${{E}_{{\rm iso}}}\sim 6.6\times {{10}^{46}}\;{\rm erg}$ over 1 s, detecting about 55% of the brightest counterparts (due to Earth occultation) and about 20% of counterparts placed at 40 Mpc (the most distant simulations due to GW sensitivity). The simulated signal was several hundred times weaker than a typical sGRB, implying a wide range of possible GRB luminosities that could be probed using GW triggers up to advanced LIGO distances. The ASM follow-up pipeline was able to detect a typical sGRB afterglow signal with about 40% efficiency at 10 Mpc.

The strategy of targeted offline follow-up of GW candidates in EM data presented here is one of a couple possibilities for joint GW-EM observation—others being rapid online GW analysis to alert pointed telescopes, and using known EM events like GRB alerts to do a targeted specialized search in GW data. All of these strategies will likely be in use in the advanced detector era. The main benefits of having archival EM data at hand is that survey data is continuously available, and not rate or schedule limited. Thousands or even millions of GW candidates can be followed-up offline, provided there is a low-enough probability of accidental coincidence. In the case of a confident GW detection, automated archival EM searches will also be able to look for subthreshold/untriggered (Gruber 2012; Blackburn et al. 2013), potentially exotic EM emission scenarios (e.g., soft flares which are difficult to target on-board due to large background variation), and in the case of no counterpart, place upper limits on the EM emission thorough software simulations.

While Fermi will continue operating in the foreseeable future, the RXTE mission ended in early 2012 after 16 yr of operation. In the advanced LIGO era, the MAXI mission (Matsuoka et al. 2009) on board the ISS has a similar wide-field soft X-ray camera and continuously surveys the X-ray sky at a sensitivity several times better than RXTE-ASM. Another promising candidate for X-ray follow-up is ISS-Lobster (Camp et al. 2013)—a NASA-proposed wide field ($30{}^\circ \times 30{}^\circ$) soft X-ray imaging telescope which would also be able to survey the sky (as well as do directed observations). The wide field imager on ISS-Lobster would be ∼100 times as sensitive as RXTE-ASM, matching well the expected factor of 10 increase in GW horizon distance for advanced LIGO. Moreover the ability of ISS-Lobster to repoint for an online GW trigger would allow it to begin observing a rapidly fading afterglow signal as soon as it appears on the right side of the Earth. Both Fermi-GBM and ISS-Lobster would be able to search for counterparts with or without the help of a galaxy catalog due to their extremely wide FOV.