Searching for electromagnetic counterpart of LIGO gravitational waves in the Fermi GBM data with ADWO
Key Words.:gamma rays: general – gravitational waves – (stars:) gamma rays bursts: general – (stars:) gamma rays bursts: individual
Aims:The Fermi collaboration identified a possible electromagnetic counterpart of the gravitational wave event of September 14, 2015. Our goal is to provide an unsupervised data analysis algorithm to identify similar events in Fermi’s Gamma-ray Burst Monitor CTTE data stream.
Methods:We are looking for signals that are typically weak. Therefore, they can only be found by a careful analysis of count rates of all detectors and energy channels simultaneously. Our Automatized Detector Weight Optimization (ADWO) method consists of a search for the signal, and a test of its significance.
Results:We developed ADWO, a virtual detector analysis tool for multi-channel multi-detector signals, and performed successful searches for short transients in the data-streams. We have identified GRB150522B, as well as possible electromagnetic candidates of the transients GW150914 and LVT151012.
Conclusions:ADWO is an independently developed, unsupervised data analysis tool that only relies on the raw data of the Fermi satellite. It can therefore provide a strong, independent test to any electromagnetic signal accompanying future gravitational wave observations.
We present a new method to search for non-triggered, short-duration transients in the data-set of the Fermi Gamma-ray Space Telescope’s Gamma-ray Burst Monitor (GBM) (Carson, 2007; Meegan et al., 2009). The method, called Automatized Detector Weight Optimization (ADWO), combines the data of all available detectors and energy channels, identifying those with the strongest signal. This way, we are able to separate potential events from the background noise and present the statistical probability of a false alarm. Although it is possible to apply our ADWO method to look for non-triggered short gamma-ray bursts (SGRBs), ADWO works the best if a potential event at a given time (and, if available, a given celestial position) is provided as an input. Thus, ADWO is ideal to search for electromagnetic (EM) counterparts of gravitational wave (GW) events, when the time of the event is well known from the GW-detectors’ observation.
On September 14, 2015 at 09:50:45.391 UTC the two detectors of the advanced Laser Interferometer Gravitational-Wave Observatory (LIGO) simultaneously observed a transient gravitational-wave signal GW150914 (Abbott et al., 2016c). The signal is originated from the merger of a binary black hole (BBH) system at low redshift () (Abbott et al., 2016a).
GBM observations revealed a weak transient source above 50 keV, 0.4 s after the GW event, with a false alarm probability of 0.0022 (Connaughton et al., 2016). This weak transient, with a duration of s, does not appear to be connected with any other previously known astrophysical, solar, terrestrial, or magnetospheric activity. Its localization is ill-constrained but consistent with the direction of GW150914. The duration and spectrum of the Fermi transient event suggest that the radiation was arriving at a large angle relative to the direction where the Fermi Large Area Telescope (LAT) was pointing.
The electromagnetic transient was a result of a custom pipeline looking for prompt gamma-ray counterparts in GBM (Blackburn et al., 2015; Kelley et al., 2013), optimized for LIGO/Virgo GW candidate events. The automatic GBM pipelines (looking for GRBs) did not find any transients.
Neither the Fermi LAT observation (Ackermann et al., 2016) above 100 MeV nor the partial Swift follow-up (Evans et al., 2016) in the X-ray, optical and UV bands, nor the INTEGRAL observations (Savchenko et al., 2016) in the gamma-ray and hard X-ray bands found any potential counterparts to GW150914, they only provide limits on the transient counterpart activity.
However, from a theoretical point of view, EM counterparts such as short duration GRBs associated with GW events are not excluded. Recently, Perna et al. (2016) proposed a scenario where a double black hole merger is accompanied by a SGRB. The evolution of the system starts with two low-metallicity massive stars that are orbiting around each other (de Mink et al., 2009; Marchant et al., 2016; Mandel & de Mink, 2016). Their orbit is so tight initially that their rotational periods are synchronized with the orbital period. Due to the fast rotation, these stars evolve homogeneously and never expand (as described by Szécsi et al., 2015, for single, homogeneously evolving stars). This way, the stars avoid the supergiant phase and thus a common envelope evolution, which reduces the theoretical uncertainties involved. Assuming that (at least) one of the supernova explosions leaves a long-lived disk behind, Perna et al. (2016) predict that this scenario leads to a relativistic jet to be launched during the merger of the black holes. The burst-duration timescale they derive from their models is in the order of 5 ms. In light of these theoretical models that predict not only the existence of black hole mergers but even the subsequent production of a SGRB, it is quite reasonable to look for EM transients of any possible gravitational wave detection.
LVT151012, the second GW candidate transient event occurred on October 12, 2015 at 09:54:43.555 UTC (Abbott et al., 2016b; The LIGO Scientific Collaboration et al., 2016). They report a false alarm probability of 0.02, and consider it not to be low enough to confidently claim this event as a real GW signal.
This paper is organized as follows: in Section 2 we describe our method, in Section 3 we test our ADWO method with the short-duration GRB150522B and in Section 4 with the SGRB-like signal that accompanied the GW150914 event. We find that our analysis of these signals are in accordance with the results of Connaughton et al. (2016). In Section 5 we apply ADWO to look for a potential EM counterpart of the event LVT151012.
2 Input data and methods
2.1 Fermi GBM overview
The Fermi GBM includes two sets of detectors: 12 thallium activated Sodium Iodide (NaI(Tl)) and two Bismuth Germanate (BGO) scintillation detectors (Meegan et al., 2009). The NaI(Tl) detectors measure the low-energy spectrum (8 keV to MeV) while the BGO detectors have an energy range of keV to MeV. The effective area of the detectors varies with the photon energy and the angle of incidence, with a maximum of cm (NaI(Tl)) and cm (BGO).
Signals from the photomultipliers are analyzed on-board, and the pulse height analysis (PHA) converts the peak heights into 128 PHA channels. The signal distribution in these PHA channels as a function of the incoming photon energy and geometry is described by the detector response matrix (DRM). The DRMs contain the effective detection area as the function of the angular dependence of the efficiency, energy deposition and dispersion, detector non-linearity, as well as the atmospheric and spacecraft scattering. The PHA distribution is usually wide for high-energy photons (especially above MeV), as some photons will scatter prior to detection. The DRMs are provided as a standard data product for each GBM trigger, but neither the program nor the data are public.
It is important to note that the 128 PHA channels have different energy ranges from detector to detector, according to the detector’s actual setup. The PHA channels are aggregated into different data products, e.g. CTIME data, which consist of accumulated spectra from each detector with a 8-channel energy and 64/265 ms time resolution.
A GBM trigger occurs when the count rates of two or more detectors exceed the background with a given threshold (). The trigger algorithms include four energy ranges ( keV, keV, keV, and keV) and ten timescales (from 16 ms to 8.192 s). A total of 120 different trigger algorithms can be specified, from which usually 75 operate simultaneously.
2.2 Automatized Detector Weight Optimization (ADWO)
The basic problem of the event analysis is to find the parameters of an event in multi-detector multi-channel time series when the approximate time and direction of the expected signal are given. To calculate the significance of such an event as described by PHA counts, one should take the typical background noise and the spectral model into account. To obtain the background-induced PHA counts, the assumed synthetic spectrum is multiplied by the DRM and binned. This is then compared to the PHA counts derived from the combination of the signal and the background with tools such as XSPEC for fitting Gaussian signals using , and C-Stat for Poisson signals (Arnaud, 1996).
Contrary to this detection method, here we do not assume the event time, only a possible time interval is given. Our goal is to find the strongest weights and the best time position in this interval using a weighted signal from the multi-detector multi-channel continuous data. The simplest method would be to compare the sum of the count rates within and outside the given time interval. This approach, however, is not the most effective one in a multi-channel multi-detector environment, since for a maximum signal-to-noise ratio usually only those detectors should be summed (selected for the analysis) which produce the strongest signals. Noisy energy channels and not illuminated detectors with very low DRM should either not be taken into account, or only with a low weight. A further complication arises from the fact that we know neither the direction of the event (and, therefore, if a given detector is illuminated or not), nor the spectra.
Our solution for these problems is the following: we give different weights to different energy channels () and detectors (), and optimize the Signal’s Peak to Background’s Peak Ratio (SPBPR). The weights are positive and normalized as . We do not restrict these weights any further, i.e. we do not include any information about the DRM (which we do not know anyway, without any spectral and directional information).
If the background subtracted intensity in the th detector’s th energy channel is , we define our composite signal as . The Signal’s Peak is the maximum of within the given time interval, and the Background’s Peak is the maximum of outside this interval. Our goal is to maximize the ratio of these two maximums. The best weights will be built up by iteration, maximizing SPBPR as a function of and . These and weights create an optimal filter among the spectra and detectors. The algorithm will provide not only maximum value of SPBPR, but will search the best weights and the exact time of this maximum, within the pre-defined interval.
We call this algorithm the Automatized Detector Weight Optimization. ADWO is similar to the GRB satellites’ triggering mechanism, but includes several improvements. For example, while Fermi’s trigger algorithm selects the and factors to be or , here we allow intermediate values too. Additionally, the condition that at least two detectors exceed a threshold simultaneously, is not required anymore, since the ADWO algorithm will automatically produce the best weights. For a signal with time-evolving spectrum, ADWO will determine the best trigger time position.
We applied Matlab’s/Octave’s fminsearch routine to find the maximum via the Nelder & Mead Simplex algorithm. The algorithm always started from an equal weights position. The analysis of signal data points against the background data points with 80 dimensional data took several minutes on a 4-core Intel i7 processor with 8GB memory, depending on the linear algebra packages used by the programs. After the search converged, the differences of the weights on the final simplex were below (the sum of the weights is 1). The sample Matlab/Octave code is available on GitHub (https://github.com/zbagoly/ADWO).
2.3 Analysis of the Fermi GBM data
Since November 2012, the Fermi GBM continuous time-tagged event (CTTE) data is present for each detector with a time precision of s, in all the 128 PHA energy channels (Meegan et al., 2009). Here we use the same CTIME energy channels of Connaughton et al. (2016), with limits of 4.4, 12, 27, 50, 100, 290, 540, 980 and 2000 keV (denoted with , resp.). Since we look for spectrally hard events, we use only the upper 6 energy channels in the 27-2000 keV range (). The exclusion of the low energy channels also reduces the background contamination from soft particle events, such as Cygnus X-1 and other weak variable X-ray sources, since their flux is usually small above 27 keV.
All the 12 NaI(Tl) and both BGO detectors were included in the analysis. Since the BGO detectors’ low energy PHA channels start above keV, the corresponding 27-100 keV energy channels are empty. Overall, we have non-zero time series.
For each detector and for each channel, the CTTE s event data is filtered with a 64 ms wide moving average filter at 1 ms steps, producing the light curve. This filtering is important as the photon event data are quite sparse (the intensity is quite low; for the GW150914 event there is, on average ms between photons in a given detector and energy channel). Our 64 ms window contains photons on average, while this window size corresponds to the typical triggered CTIME light curve resolution.
Without filtering, the photon-photon correlation in time that we search for would disappear. Very narrow filters are worthless because of the sparsity constraint, while much wider filters will smooth and filter out short transients, lowering ADWO’s sensitivity. As a byproduct, the smoothing also acts as a low-pass filter which reduces the Poisson noise. The moving average filter is the simplest choice here: e.g. using some prior knowledge about the signal’s shape, a matching filter tuned to the signal would improve the sensitivity at the expense of generality.
Fermi operates in survey mode most of the time, with slewing at 4 degrees per minute. This creates a continuously changing background, which should be accounted for, since ADWO would be optimal without directional changes (as it uses the correlation between the detectors and channels). One possibility would be to take the detailed satellite positional information into account and create a physical model to determine the background for a hundreds of seconds (Szécsi et al., 2013). However, we expect that the slow slew will not suppress the sensitivity to the kind of short (sec) transients that we are looking for. Therefore, a much simpler, 6th order polynomial background fit was subtracted for each channel and detector, similar to the method of Connaughton et al. (2016). We chose the typical background window to be s around the search window, depending on the CTTE data availability: this window can contain the majority of GBM’s GRB ligtcurves and covers approx. of Fermi’s orbit.
To test ADWO, we analyze the short GRB150522B gamma-ray burst, with Ts and erg/cm fluence. These parameters are comparable to the EM companion values of GW150914, as reported by Connaughton et al. (2016). Fermi triggered on May 22, 2015 at 22:38:44.068 UTC, and full CTTE data of s interval relative to the trigger is analyzed, using a 6 s long signal window centered on the trigger. The ADWO obtains a maximal SPBPR of 3.12, and reveals the double pulse shown in the Fermi GBM quicklook data (Fig. 1). The detector and energy channel weights are given in Tables 1-2.
To determine the significance we generated a Poisson-distributed synthetic signal, using the background photon data of the interval, and repeated ADWO for Monte-Carlo (MC) simulations with the same window width. There was no simulation with bigger SPBPR value than 3.12, therefore we estimate the false alarm rate to be below , and the false alarm probability to be below , analogously to Connaughton et al. (2016).
4 The GW150914 event
We apply the ADWO method on the Fermi CTTE data set covering the event of GW150914: the 6 s long signal window was centered on September 14, 2015 09:50:45 UTC (391ms before trigger). Here we investigate a s time background interval that adds up as s before and s after the time of the possible event. The ADWO has converged (Fig. 2) and the obtained maximal SPBPR is 1.911, 474 ms after the GW trigger.
We repeated ADWO for MC simulations using this data: 86 cases had bigger SPBPR than 1.911. The false alarm rate is , giving a false alarm probability of , which is higher than 0.0022, the value given by Connaughton et al. (2016).
Considering the positional errors of GW150914 on the sky, it can be easily shown that there’s a high () probability that a similar error ring will intersect with Fermi GBM’s field of view in the case of LVT151012 too. Therefore, we apply ADWO on the Fermi GBM CTTE data around the event of LVT151012, covering s, centered on October 12, 2015 at 09:54:43 UTC. We find a relatively strong signal at 09:54:44.207 UTC in the 6 s signal window, with a SPBPR of 1.805 (Fig. 3). The sum of the keV weights is higher than in the case of GW150914, i.e. here the signal is softer than GW150914 at the peak ( MeV), but harder than GRB15522 at the peak ( keV).
We made MC simulations: 308 cases had bigger SPBPR than the original observation. hence the false alarm rate is , and the false alarm probability is estimated to be .
When cross-checking the lightning detections made by WWLLN (Rodger et al., 2009) with the Fermi’s positions and times, we find no TGF candidates (storm activity) within 500km of the spacecraft position and s around the peak.
Although here we applied our ADWO method to look for particular events, we point out that it is entirely possible to use this unsupervised data analysis method for a general search for non-triggered, short-duration Fermi GBM events. Automatized search processes are important, as the total data-set collected by the Fermi’s 8-years operation is significantly larger than the triggered data-set. It is likely that there are several potential EM events observed but not triggered. For example, based on the CTIME 256ms data product, Gruber & Fermi/GBM Collaboration (2012) estimates untriggered SGRB/month in the Fermi observations. It is thus a worthwhile future task to identify potential SGRB candidates in the non-triggered Fermi GBM data-set using ADWO. Alternatively, we can cross-check those already found by other algorithms.
As our ADWO method is independently developed, and only relies on the raw data of the satellite, it can provide a strong, independent test to any future signal. In regard of the current expectation that LIGO will detect several GW events in the near future, many of which may have a weak EM transient counterpart such as a SGRB, it is of crucial importance to identify those potential EM signals. We therefore expect that ADWO will be successfully applied in the future to find SGRB counterparts of the GW events observed by LIGO. The analysis of the GW151226 event as well as the improvement of ADWO is the topic of a forthcoming paper.
Acknowledgements.This work was supported by OTKA grants NN111016 and NN114560. The authors wish to thank the World Wide Lightning Location Network (http://wwlln.net), a collaboration among over 50 universities and institutions, for providing the lightning location data used in this paper. Thanks for the computational resources of the Wigner GPU Laboratory of the Wigner RCP of the H.A.S, A.C. Elbeze for pointing out that GRB150522B was the second GRB of the day and the anonymous referee for valuable comments on the paper.
- Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJ, 818, L22
- Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, Phys. Rev. D, 93, 122003
- Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, Physical Review Letters, 116, 061102
- Ackermann, M., Ajello, M., Albert, A., et al. 2016, ApJ, 823, L2
- Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
- Blackburn, L., Briggs, M. S., Camp, J., et al. 2015, ApJS, 217, 8
- Carson, J. 2007, Journal of Physics Conference Series, 60, 115
- Connaughton, V., Burns, E., Goldstein, A., et al. 2016, Accepted for publication in ApJL, ArXiv e-prints [\eprint[arXiv]1602.03920]
- de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243
- Evans, P. A., Kennea, J. A., Barthelmy, S. D., et al. 2016, MNRAS[\eprint[arXiv]1602.03868]
- Gruber, D. & Fermi/GBM Collaboration. 2012, in Proc. of the Gamma-Ray Bursts 2012 Conference (GRB 2012). May 7-11, 2012. Munich, Germany., 36
- Kelley, L. Z., Mandel, I., & Ramirez-Ruiz, E. 2013, Phys. Rev. D, 87, 123004
- Mandel, I. & de Mink, S. E. 2016, MNRAS, 458, 2634
- Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50
- Meegan, C., Lichti, G., Bhat, P. N., et al. 2009, ApJ, 702, 791
- Perna, R., Lazzati, D., & Giacomazzo, B. 2016, ApJ, 821, L18
- Rodger, C., Brundell, J., Holzworth, R., & Lay, E. 2009, in Am. Inst. Phys. Conf. Proc., Coupling of thunderstorms and lightning discharges to near-Earth space: Proceedings of the Workshop, Corte (France), 23-27 June 2008, Vol. 1118, 15–20
- Savchenko, V., Ferrigno, C., Mereghetti, S., et al. 2016, ApJ, 820, L36
- Szécsi, D., Bagoly, Z., Kóbori, J., Horváth, I., & Balázs, L. G. 2013, A&A, 557, A8
- Szécsi, D., Langer, N., Yoon, S.-C., et al. 2015, A&A, 581, A15
- The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, B. P., et al. 2016, ArXiv e-prints [\eprint[arXiv]1606.04856]