ON THE FERMI-GBM EVENT 0.4 s AFTER GW150914

Due to the nature of GBM backgrounds, it is difficult to define a spectral background model. The background contains a superposition of astrophysical, terrestrial, solar, and instrumental sources, each with its own time variability. Therefore, an assumption is made that in each GBM energy channel, the background is smoothly varying with time. With this assumption, one selects a time interval before and after the source time interval and fits a polynomial in time to the count evolution to these time intervals but excluding the source time interval. This polynomial is then interpolated through the source time interval to estimate the background counts in each energy channel during the time of source activity. A major assumption of this approach is that the background does not vary during the source time interval more than the smooth nature apparent in the surrounding background time interval. This assumption is most likely valid if the source time interval is short compared to the timescale of the variations of the fitting polynomial.

To determine the polynomial shape of the temporal background count evolution, a fit to statistically fluctuating data must be performed. Under the assumption that these counts should be Poisson distributed, the statistic to choose for maximum-likelihood estimation (MLE) is the Poisson likelihood (e.g., the so-called Cash statistic; Cash 1979):

$$\begin{eqnarray}&&L=\displaystyle \prod _{i=1}^{N}\displaystyle \frac{{M}_{i}^{{S}_{i}}{e}^{-{M}_{i}}}{{S}_{i}!},\end{eqnarray} \tag{ 1 }$$

where M_i and S_i are the model and detected counts in the ith of N bins, respectively. Taking the $-2\mathrm{log}L$ and approximating the factorial term with Stirling's approximation (Abramowitz & Stegun 2002),⁶ we have

$$\begin{eqnarray}&&-2\mathrm{log}L\approx 2\displaystyle \sum _{i=1}^{N}\,{M}_{i}-{S}_{i}\mathrm{log}{M}_{i}+{S}_{i}\mathrm{log}{S}_{i}-{S}_{i}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&=2\displaystyle \sum _{i=1}^{N}\,{M}_{i}-{S}_{i}+{S}_{i}(\mathrm{log}{S}_{i}-\mathrm{log}{M}_{i}).\end{eqnarray} \tag{ 3 }$$

The polynomial coefficients can be estimated via minimizing the log-likelihood of the background data and the polynomial integrated over each time bin.

In the RMFIT⁷ software, generally used for GBM GRB data, the approach for determining the background is similar to the above method, except a weighted least-squares (${\chi }^{2}$) statistic is chosen to fit the polynomial coefficients via a two-pass method. First, the polynomial coefficients are estimated by solving the least-squares linear equations assuming that the data variance is equal to the number of counts. Then, the fit is performed again, but using the variance of the model determined in the first pass. This is intended to compensate for the fact that at low counts, the Gaussian approximation of the Poisson distribution is poor.

In the high-count regime, this approach may be valid. However, due to the fact that a polynomial is fit in each energy channel and that the background has fewer counts at high energies, the ${\chi }^{2}$ statistic can be very biased (Baker & Cousins 1984; Humphrey et al. 2009; Andrae et al. 2010). To estimate this effect, we simulate synthetic light curves obeying a zero-order polynomial evolution in time over a period of 20 s. We choose a set of levels for the background mean rate (θ) ranging from 1 to 100 cts s⁻¹. For each θ, we simulate 500 light curves, and fit a polynomial to them with two different statistics criteria: Cash and ${\chi }^{2}$. We find that the Cash statistic well recovers the true θ even at very low counts, while both ${\chi }^{2}$ methods have a large bias at low counts (see Figure 2). Therefore, we will use the Poisson likelihood to fit the background in our spectral analysis, but will compare the results to what is found using RMFIT and ${\chi }^{2}$.

In the low source count regime, the choice of spectral fitting statistic is crucial. The total count rate in a source interval is distributed as a Poisson process, but we have more information via the background fit to better constrain the spectral fit parameters. The total counts in a spectral bin (S_i), is composed of source and background (B_i). Since we have estimated the background via a polynomial model, we can include this information in the Poisson likelihood. Let $P(t;{\boldsymbol{\theta }})$ be the estimated polynomial background rate model in a given energy channel where t is time and ${\boldsymbol{\theta }}$ are the polynomial coefficients. Then, for a given source interval from time T₁ to T₂, we have

$$\begin{eqnarray}&&{B}_{i}={\int }_{{T}_{1}}^{{T}_{2}}P(t;{\boldsymbol{\theta }}){dt}.\end{eqnarray} \tag{ 4 }$$

Now, ${\boldsymbol{\theta }}$ are fitted coefficients with associated statistical errors, and therefore B_i also has statistical error (${\sigma }_{B}$). Using the covariance matrix for ${\boldsymbol{\theta }}$, we can estimate ${\sigma }_{B}$ via standard Gaussian error propagation⁸ (see Barlow 1993 for a review of standard Gaussian error propagation). Error propagation assuming that the resulting error distribution is Gaussian leads to Gaussian uncertainties on the model estimated B_i quite naturally.⁹

We must use a likelihood for Poisson data with a Gaussian background to fit the spectral data. This is embodied in the PG-statistic (Arnaud et al. 2015):

$$\begin{eqnarray}-2\mathrm{log}L & = & 2\displaystyle \sum _{i=1}^{N}\,{M}_{i}+{t}_{s}{f}_{i}-{S}_{i}\mathrm{log}({M}_{i}{t}_{s}{f}_{i})\\ & & +\displaystyle \frac{1}{2{\sigma }_{B,i}^{2}}{({B}_{i}-{t}_{s}{f}_{i})}^{2}-{S}_{i}(1-\mathrm{log}{S}_{i}).\end{eqnarray} \tag{ 5 }$$

Here, t_s is the source interval duration, and f_i is the profiled-out background model.¹⁰ This is the statistic we will minimize in our spectral analysis.

In the RMFIT software, the background is assumed to be perfect with no variance. Therefore, only the total counts are assumed to have statistical variability. To deal with this, RMFIT modifies the Poisson likelihood by adding the estimated background from Equation (4) to the estimated spectral model counts:

$$\begin{eqnarray}&&-2\mathrm{log}L=2\displaystyle \sum _{i=1}^{N}({M}_{i}+{B}_{i}))-{S}_{i}\\ &&+{S}_{i}({ln}({S}_{i})-{ln}({M}_{i}+{B}_{i}))).\end{eqnarray} \tag{ 6 }$$

Since it is possible that the method of Blackburn et al. (2015) recovers a significant fluctuation above the typical background, it is important to understand if the source has the spectral parameters reported in Connaughton et al. (2016). Therefore, we simulate the problem. We fit the background of the data before the trigger over a duration of 10 s to obtain the average background rate in each channel. We then simulate 10 s background from this sampled spectrum assuming each channel has Poisson fluctuations with their means derived from the sampled spectrum. Next, we simulate a 1 s source with a power-law spectrum. We choose the power-law index to be –1.4 and three different amplitudes: A = 2 × ${10}^{-3},\times {10}^{-4},0.0$. The first amplitude models the claim of Connaughton et al. (2016), the second simulates an order of magnitude weaker source, and finally the last amplitude implies no source. For each set of spectral parameters, we generate 100 sets of source and background simulations and create TTE data via the method described in Burgess (2014).

For each simulation, we fit the simulated source spectra with three different methods. First, we fit the background and spectrum in RMFIT. Next, we use the Poisson method to fit the background, and fit the source spectra with the standard PGStat statistic. To minimize the likelihood, we use MINUIT and our own RSP matrix convolution software for computational simplicity (hereafter MLEfit), but verified our results by producing PHA files for source and background and fitting them with XSPEC.¹¹ We effectively compare three different methods for fitting spectra. Figure 3 shows the results of the simulations where we assume the same source spectrum reported in Connaughton et al. (2016).

Zoom In Zoom Out Reset image size

MLEfit using PGStat recovers the true simulated amplitude more accurately than any other method. The same has been shown to be true for a simulation of different spectral slopes at identical amplitudes (see Figure 7.3 in Arnaud et al. 2011). RMFIT is biased to slightly higher amplitude values. When the source amplitude is reduced by an order of magnitude, MLEfit is about to find the correct amplitude with either statistic though the modified C-stat is slightly biased toward higher values. On the other hand, RMFIT recovers a source amplitude resembling what is found in Connaughton et al. (2016; see Figure 3). Similarly, when there is no source, MLEfit finds a source amplitude near zero, while RMFIT still finds a source amplitude of $\sim 3\times {10}^{-3}$. This has two implications. MLEfit is able to distinguish between source and background while RMFIT is including background in its source spectrum. The fact that RMFIT finds a source spectrum similar to what is reported in Connaughton et al. (2016) even when there is no source indicates that the background spectrum can be modeled with a power law having an amplitude and overall spectrum very similar to the reported source spectrum in Connaughton et al. (2016). We further demonstrate this in the following section.

Combining the low-count bias of ${\chi }^{2}$, which can underestimate the background level and the assumption of no variance on this underestimated background via the use of RMFIT's modified C-stat, it is possible to find a source spectrum and flux that is unrealistically too large as we will demonstrate using our approach to spectral fitting.

Following Connaughton et al. (2016), we generate 10 sets of detector response matrices (RSPs) corresponding to the 10 locations listed in the same work, using the standard stand-alone response generator. We use the GBM team-created software RMFIT to bin the data, make energy selections equivalent with those in Connaughton et al. (2016), as well as fit the background to the surrounding count light curve. We fit a power-law spectrum to the data, using the 10 different RSPs discussed above. Table 1 details these fits. Within the statistical errors, these are compatible with the spectra derived by Connaughton et al. (2016), in particular the details given in their Section 3.2 and Figure 5.

We note that we have used the full 128-channel TTE data instead of the 8-channel TTE data used in Connaughton et al. (2016) because a proper MLE statistic is accurate regardless of the binning as pointed out in Cash (1979).

The source spectrum reported in Connaughton et al. (2016) is more similar to the background (cosmic X-ray background plus Earth albedo spectrum; Ajello et al. 2008) than of short GRBs. In Section 2.3, we demonstrated that RMFIT will effectively find a source spectrum even when no source is present and is therefore fitting strong, positive fluctuations of the background. To further assess this claim, we fit 1.024 s intervals before and after the tentative detection. We iteratively select the intervals as source and fit the surrounding light curve as background. We use the Cash statistic to fit the background and MLEfit with PGStat to fit the source since we have demonstrated that the RMFIT methods are insufficient for low-count analysis. As an additional check, we also fit the spectra with RSPs generated from each of the 10 locations used in Connaughton et al. (2016).

Figure 4 shows that the spectra of the surrounding time intervals are very similar to the spectral parameters of the putative event, regardless of position on the sky used to generate the RSPs. Since the location is poorly constrained, then the range of parameters covered by these plots could be considered valid. Therefore, all recovered spectra are essentially the same indicating that the source is most likely a strong positive background fluctuation.

Zoom In Zoom Out Reset image size

Here, we demonstrate the spectrum of the source candidate using Cash statistics to fit the background, and PGStat to fit the spectra. Table 1 shows the best fit for each of the 10 locations across the LIGO localization arc. As can be seen, the amplitudes are much lower and the fluences reduced correspondingly. Also, the spectral slope is substantially softer. It is still most likely that the source is a strong background fluctuation because the spectrum is not inconsistent with zero and consistent with what is expected from the combination of extragalactic background and Earth's albedo. While the larger error in spectral slope still allows a hard spectrum, the reality of these solutions is excluded by the INTEGRAL/ACS limits (see end of this section and the next one).

Next, we examine the predicted counts from our spectral fits with MLEfit, RMFIT, and those reported in Connaughton et al. (2016). Figure 5 shows the count light curve of NaI 5 integrated over 11.5–930 keV. The background model coming from the polynomial fit is shown with its 1σ Gaussian errors. The total Poisson error of the event is not too far from the background, and indeed there are several peaks of similar amplitude within this time slice. Using the RSP corresponding to R.A. = 155.3, decl. = −43.2 and the spectral parameters corresponding to the fits at that location (see, e.g., row 2 of Table 1), we convolve the power-law photon model to produce the predicted number of counts during the event. The counts predicted from MLEfit are within the observed range of background subtracted counts, but the results from our own RMFIT fits and those of Connaughton et al. (2016; using RMFIT) overpredict the number of counts significantly.

Zoom In Zoom Out Reset image size

As a secondary check, we attempt a Bayesian fit of the data where both the background polynomials in each channel and the spectra are fit simultaneously. This avoids the ambiguity of background error propagation and allows us to directly assess the background's effect on the spectral fit. The method will be detailed in J. M. Burgess (2016, in preparation); however, we briefly detail the approach here. The full temporal and spectral model is defined as a piecewise function in time:

$$\begin{eqnarray*}f({t}^{i},\varepsilon ;{\theta }_{n}^{j},{\boldsymbol{\phi }})=\left\{\begin{array}{l}B({t}^{i};{\theta }_{n}^{j})\& {t}^{j}\lt {t}_{{\rm{start}}}\\ B({t}^{i};{\theta }_{n}^{j})+S({t}^{i},\varepsilon ;{\boldsymbol{\phi }})\& {t}_{{\rm{start}}}\leqslant {t}^{j}\leqslant {t}_{{\rm{stop}}}\\ B({t}^{i};{\theta }_{n}^{j})\& {t}^{j}\gt {t}_{{\rm{stop}}},\end{array}\right.\end{eqnarray*}$$

where tⁱ is the time in the $i{\rm{th}}$ interval, ${\theta }_{n}^{j}$ is the $j{\rm{th}}$-order polynomial coefficient in the $n{\rm{th}}$ energy channel, B is the background polynomial, and S is the source photon model with ${\boldsymbol{\phi }}$ parameters. The ${t}_{{\rm{start}}}$ and ${t}_{{\rm{stop}}}$ correspond to the source interval times. For the data at hand, we use a zero-order polynomial, B in each channel, and again choose a power-law function for S. A log-uniform prior is chosen for ${\theta }_{n}^{0}$ and the power-law amplitude and a uniform prior for the power-law spectral index. Unlike the previous analysis, we are fitting both the background and source simultaneously, so we can take both the background and source as distributed via the Poisson process.

Figure 6 shows the resulting parameter posterior contours indicating that the parameters are mostly unconstrained and the amplitude is not inconsistent with zero. Moreover, Figure 7 shows the light curves from the data with the background and source+background posteriors for the NaI and BGO detectors. The source+background histogram is clearly consistent with the background only, suggesting that the source is a background fluctuation and not an astrophysical source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We again simulate the spectral parameters in Connaughton et al. (2016) to test if this method could recover the parameters of a real source. Figures 8 and 9 show the results indicating that a true source with the reported spectral parameters would be recovered by the method.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To access the presence of an actual source, we use the Bayesian tool to fit the data using the background polynomial only. Using the Monte Carlo samples from the source+background fit and the background-only fit, we compute the deviance information criteria (DIC; Spiegelhalter et al. 2002) for each model, which allows us to perform model comparison between them. The DIC compares model fit and penalizes for model complexity, which can be a function of both data and model. We find a DIC of 18341 and 18347 for background-only and source+background, respectively. Choosing the lowest DIC implies that the background-only model is the best description of the data. This proposes an interesting contrast to the claimed detection in Connaughton et al. (2016) via the method of Blackburn et al. (2015). It is likely that the included information in the simultaneous fit of the full resolution data is a better characterization, and this method correctly finds no significant source in agreement with theory and the non-detection in all other observatories.

Zoom In Zoom Out Reset image size

Based on the previous analysis, we conjecture that the mismatch in the resulting parameters depends on the statistical significance of the residual source counts per event. In order to test this conjecture, we reanalyzed with PGStat all 58 short GRBs from the GBM time-integrated 4 year catalog (Gruber et al. 2012), which had the power-law model as BEST-fit model. Figure 10 shows the comparison of the resulting parameters with those obtained using RMFIT (and listed in Gruber et al. 2012). We find that PGStat reproduces the fit parameters of RMFIT for strong GRBs, while for faint GRBs the spectra are systematically softer and the fluences smaller, exactly as shown in the above analysis of the GBM event.

Zoom In Zoom Out Reset image size

The above analysis demonstrates the systematic effect of the differences between the RMFIT package and the PGStat analysis. We can go one step further and ask, how do the different best-fit parameters affect the prediction of the detectability by INTEGRAL/ACS, and how do these predictions compare with the actual detections by INTEGRAL/ACS? Taking two short GRBs from the Gruber et al. (2012) catalog that belong to the 20% lowest fluence group, we fold the spectral parameters (with their catalog sky position) through the INTEGRAL/ACS response and obtain predictions in the 3σ–4σ range. Checking the INTEGRAL/ACS, we indeed find these GRBs at these confidence levels, demonstrating that our PGStat analysis provides results that are verified by detections with an independent instrument. These two short GRBs are added in Figures 10 and 11. As previously, RMFIT overpredicts the counts and, correspondingly, the fluence in INTEGRAL/ACS.

Zoom In Zoom Out Reset image size

Typical short GRBs have a spectral index of –1.4 ± 0.2 (Gruber et al. 2012), and our PGStat analysis of the faintest triggered short GRBs finds them at indices as soft as −1.6 (with large errors; see Figure 11). The best-fit spectral index of the GBM event close to GW150914 GRB 120212353 is −1.6 (for a position outside the LIGO arc), and the spectral indices for positions along the LIGO arc cover the range −2.0... −1.5 (also with large errors; see Figure 11). While there is some overlap between the two distributions, the GBM–GW150914 event is systematically softer than even the softest/faintest short GRBs.

We can take it yet another step further and test the validity of the PGStat analysis using GBM subthreshold short GRBs, i.e., GRBs that were detected by Swift/BAT, but not by GBM. We pick two GRBs (one faint and one stronger) out of the seven for which Burns et al. (2016) determined visibility by GBM and which were recovered in ground-based searches of the untriggered GBM data (see Table 4 of Burns et al. 2016). For these bursts, the PGStat analysis results in slightly softer and dimmer spectra than the analysis with RMFIT, consistent with the above analysis. The uncertainties in the fluence are small, and exclude zero fluence, i.e., the result recovers proper flux from these two GRBs. This shows that the PGStat analysis is capable of recovering real GBM subthreshold events. In contrast, PGStat analysis of GW150914-GBM indicates an amplitude consistent with zero.

To further validate the Bayesian method, we apply it to one of these weaker GRBs (GRB 140516A). The method is able to clearly distinguish between source and background (Figure 12) with a DIC preferring a source even though the total count rate above background is relatively low (even lower then the GBM-GW150914!). This further strengthens our claim that were GBM–GW150914 a real astrophysical short GRB, the method would be able to identify it. Caution must be taken when assessing the presence of a source based on total count rate alone when spectral information is available. The additional spectral information, as well as the variability of the background, place severe constraints on the presence of weak events.