SEARCHING FOR NEEDLES IN HAYSTACKS—LOOKING FOR GAMMA-RAY BURST γ-RAYS WITH THE FERMI/LAT DETECTOR

One of the more surprising results of the Compton Gamma-Ray Observatory was the EGRET discovery of an 18 GeV photon associated with GRB 940217 (Hurley et al. 1994). This was not an entirely unique event; about a half-dozen bursts were seen over the course of the mission with photons above 100 MeV (Catelli et al. 1998; Dingus 2003). Since the GRB spectral energy distribution at lower energies has been well characterized by a modified power law with peak fluxes at energies of the order of 200 keV, the existence of photons at energies 10⁴ times higher puts a significant constraint on any viable model of the GRB phenomenon (e.g., Band et al. 2009). The obvious question is whether high-energy photons are associated with all GRBs or only with a small sub-class. We are now in a unique period when two space missions are providing GRB coordinates with arcsecond accuracy (Swift; Gehrels et al. 2004) and the ability to measure photon energies to 300 GeV (Fermi/LAT; Atwood et al. 2009). Once these two missions cease functioning, the resources to investigate GRBs over a broad energy range are not likely to be available for decades. Thus, it is a matter of some urgency to extract the maximum amount of information from the assets now at hand. Prior to the launch of the Fermi Gamma-ray Space Telescope, it was possible to speculate that the LAT instrument would detect more than 200 GRB events per year (Dingus 2003). In the two-year period since the launch of the Fermi Gamma-ray Space Telescope, the Gamma-ray Burst Monitor (GBM; Meegan et al. 2009) has reported approximately 475 GRBs, i.e., a rate of about 250 per year. Over essentially the same period, only 17 bursts³ have been identified by the Fermi/LAT. Given that roughly 50% of all GBM triggers point to directions within 66° of the LAT boresight, the fraction of bursts with photons above 100 MeV is apparently not much larger than 7%. Thus, it is an important task to apply the best possible techniques to explore whether this fraction would be significantly bigger if the detection threshold could be lowered. Such studies might significantly inform the design of future missions while extracting additional information from rather expensive facilities.

From launch through 2010 March 16, 55 GRBs have been detected (mainly by Swift, as well as INTEGRAL: Winkler et al. 2003; AGILE: Giuliani et al. 2008) at celestial coordinates that were simultaneously viewed by the Fermi/LAT. The angular resolutions of the Swift BAT and XRT are arcminutes and arcseconds, respectively. Compared with the mean point-spread width for the Fermi/LAT of a few degrees, these GRB coordinates are known with zero error. This sample is an excellent target for statistical methods that can take advantage of the precisely determined source direction. The basic technique that is employed here is the matched filter, most familiar to those detecting signals in the time domain. The underlying assumption of this method is that the characteristics of both the signal and background are a priori known functions of one or more variables. Since the matched filter maximizes the signal-to-noise ratio, moderate departures from optimality degrade the filter performance relatively slowly, making this a valuable tool for investigating the possible existence of faint signals. We present the sample and data processing in the following section, and details of the signal detection technique in Section 3. Results are given in Section 4 followed by a brief discussion in Section 5.

A selection of 148 precisely located GRBs (<18) was obtained from the Swift Web page⁴ devoted to cataloging all bursts. The data correspond to the period from launch until the end of 2010 March. Most of these events were observed by the Swift/XRT, augmented by a few from INTEGRAL and AGILE. The resulting list was cross-correlated with Fermi spacecraft attitude data obtained from the Fermi Science Support Center⁵ (FSSC) to identify only those events within 66° of the LAT boresight.

Having identified the principal focus of this investigation, two sets of LAT photon data were obtained from the FSSC corresponding to the 55 precisely located GRB fields and to 464 additional fields taken at random on the sky to study the background behavior. Throughout the following, they will be identified as "GRB" or "random" data. The FSSC data selection criteria corresponded to any of the three event classes, "transient," "source," or "diffuse," photon energies above 100 MeV but below 300 GeV, a zenith angle less than 105° and time intervals when the DATA_QUAL parameter was 1 and the South Atlantic Anomaly parameter, IN_SAA, was 0. The 105° zenith-angle limit might be considered a bit risky but with an orbiting altitude of 565 km, the Earth's limb appears at a zenith angle of 11283 which is easily resolvable. For the random LAT fields, fictitious GRB directions were chosen randomly over the sky but constrained by the less than 66° boresight requirement. Further cuts were imposed of less than 105 for the angle between the photon and GRB directions and a time window spanning from t₀ − 100 s to t₀ + 150 s. With LAT point-spread width maxima of the order of 4°, the 105 cut was imposed to allow up to 2.6σ errors without worrying about the effects of significant variations in the background flux over the effective solid angle of 0.105 sr. All data fields were additionally constrained to lie above or below the Galactic plane by more than 10°. Such additional criteria reduced the number of GRB fields from 55 to 41. Six of these have been previously associated with significant fluxes of energetic photons (>100 MeV) and thus provided guidance about the characteristics of the signals we were seeking.

The central theme of this paper is the use of a matched filter to enhance the signal-to-noise ratio for the detection of GRB photons in a diffuse background. There are four independent variables that are involved: energy (E), photon angle with respect to the GRB direction (θ), photon time of arrival (t), and photon event class (c ∈ {1, 2, 3}). The most obvious consideration is that the photon direction should be correlated with the GRB coordinates within an angle related to the LAT point-spread function (PSF) for that specific photon energy. Photons that fall far outside this criterion are almost certainly background. The energy criterion is less intuitive. The energy spectra for both the GRB signal and the background are approximately represented by decreasing power laws with exponents Γ_GRB and Γ_back, respectively. Unfortunately, the background has a somewhat harder spectrum with Γ_back ≃ 1.5 while Γ_GRB ≃ 2.2. Thus, higher energy photons are more likely to represent background than signal. Photon directional accuracy is also critical, so photons with small uncertainty (i.e., high energy) should be weighted more heavily. The time filter is the most problematic because GRB light curves lack a uniform shape. Particularly for the most energetic burst photons, the time delay relative to the low-energy trigger can be substantial, up to thousands of seconds (Hurley et al. 1994). However, uniquely identifying such late quanta with a GRB can be problematic, especially in view of the lack of a distinguishing spectral characteristic. One feature that is generally observed is a prompt flux within a few seconds of the low-energy trigger. This can be seen in Figure 1 of Ghisellini et al. (2010). Most of the depicted LAT light curves exhibit a significant flux within 10 s of the burst. There are a few exceptions such as GRB 090328, so it must be recognized that not every fish can be caught with the same net. This led us to adopt a standard GRB light curve described by a broken power law with a break at t_b = 5 s. Prior to t_b, the flux rises as (t/t_b)$^{\alpha _0}$ after which it falls as (t/t_b)^−α. The value for α₀ was set rather arbitrarily at 1/2, while α was picked to be 7/5 in concert with the general behavior of a large number of GRB afterglows. We emphasize that all the filter parameters were determined prior to applying them to either the GRB or random data sets.

When we began this analysis, we were unaware of the nuanced distinctions between the three LAT photon event classes, "transient," "source," and "diffuse," identified as 1, 2, and 3. Since we were looking for a transient point source, it seemed reasonable to include the "transient" quanta on equal footing with the other two. Only after going through the complete development of the matched filter search did we discover that the signal content of the "transient" photons is extremely small, ∼3.5%. This was determined by looking at the event class photon counts for the most intense GRBs, 080916C, 090510, 090902B, and 091003. The ratio, r_GRB, of "transient" and "source" GRB photons to "diffuse" is 0.34 and 0.44, respectively, roughly independent of the particular burst. On the other hand, the same ratios for the selected random LAT fields, r_back, are 9.65 and 0.85. Consequently, we realized that we needed to include a weighting factor, w_c, based on event class that reflected this wide dispersion of information content. Serendipitously, the simple existence of a strong correlation of GRBs with "diffuse" class photons provided an independent verification of the validity of the matched filter technique.

From these considerations, the mathematical form of the matched filter weight functions is given by

$$\begin{equation} w_E = \frac{1}{4\pi {\sigma }^2(E)}\left(\frac{E}{E_{\rm{th}}}\right)^{\Gamma _{\rm back}-\Gamma _{_{\rm GRB}}} \end{equation} \tag{ 1 }$$

$$\begin{equation} w_\theta = 2 e^{-\frac{{\theta }^2}{2{\sigma }^2(E)}} \end{equation} \tag{ 2 }$$

$$\begin{equation} w_t = \left \lbrace \begin{array}{@{}ll@{}} c \cdot \left(\dsty\frac{t}{t_b}\right)^{\alpha _0} & 0 \le t \le t_b \\ \ms c \cdot \left(\dsty\frac{t}{t_b}\right)^{-\alpha } & t_b < t \le t_c \\ \ms c=2.056382 & t_c = \frac{19}{2}t_b, \end{array}\right. \end{equation} \tag{ 3 }$$

$$\begin{equation} w_c =\begin{array}{ll} \left({\dsty\frac{r_{\rm GRB}(i)}{r_{\rm GRB}(3)}}\right)\!\! \Big/ \!\! \left({\dsty\frac{r_{\rm back}(i)}{r_{\rm back}(3)}}\right), & 1 \le i \le 3 \end{array} \end{equation} \tag{ 4 }$$

$$\begin{equation} w(E,\theta,t,c) = w_E \cdot w_\theta \cdot w_t \cdot w_c. \end{equation} \tag{ 5 }$$

The PSF for the LAT has been specified by the FSSC team in terms of a half-angle that includes 68% of a random sample. We have fit this to a power law in energy (Akerlof & Yuan 2007) to obtain separate parameter sets corresponding to conversion in either half of the LAT tracker using the form

$$\begin{equation} {\sigma _{68\%}(E)}=\left(\frac{E}{E_{\rm{th}}}\right)^{-\delta }{\sigma }_0+{\sigma }_1. \end{equation} \tag{ 6 }$$

For this analysis, we have assumed that the LAT PSF is Gaussian as described above for w_θ with a characteristic width determined by the containment condition:

$$\begin{equation} \sigma _{\rm PSF}(E) = \frac{2}{3}{\sigma _{68\%}}(E). \end{equation} \tag{ 7 }$$

Although the actual LAT angular error distribution function exhibits non-Gaussian tails, this approximation is unlikely to affect subsequent conclusions. As indicated above, the weight value for each photon is computed as a product of the four elements, w_E, w_θ, w_t, and w_c. The time filtering algorithm substantially ignores all photons that occur more than a few tens of seconds after the GRB trigger. For the parameters selected, this can best be described by an effective T₉₀ window time for the signal of 16.864 s starting 1.624 s after the burst. The sum over all photons, ∑w_i, is the statistic used to identify events with significant fluxes of high-energy photons. In assessing this technique, the mean contribution from background photons turns out to be negligibly small. The identification problem is driven entirely by infrequent but large fluctuations of ∑w_i contributed by background quanta. As an additional safeguard, a weight-sharing measure was introduced to ensure that the entire event weight was not generated by a single extreme photon. This parameter is defined by

$$\begin{equation} \zeta = \frac{3(w_{0}\cdot w_{1}\cdot w_{2})^{1/3}}{w_{0} + w_{1} + w_{2}}. \end{equation} \tag{ 8 }$$

In this formula, w₀, w₁, and w₂ are the maximum three photon weights for a particular event. The combined statistic, ζ∑w_i, incorporates both the sum of photon weights and the additional requirement that no single photon contributes an overwhelming fraction of the value. The downside is that at least three photons must be associated with an event, an issue that will arise later.

The significances of the ∑w_i and ζ∑w_i statistics were determined by Monte Carlo simulations of the background. For each field of the GRB data set, 2000 events were generated using the measured background spectral distribution to determine energy and uniform distributions over solid angle and time to determine θ and t. The number of background photons for each field was set by the measured background rate. By looking at the weight statistics for simulated events with and without injected GRBs, we gained reasonable confidence that the matched filter algorithm would correctly identify a significant fraction of real bursts. This work also demonstrated that the modified sharing-compensated statistic, ζ∑w_i, was the stronger measure.

As described earlier, 41 GRB fields were selected for examination. These are listed in Table 2. Six of these events have been previously identified with >100 MeV LAT-detected photons. These are GRB 080916C (Abdo et al. 2009a), GRB 090323 (Ghisellini et al. 2010), GRB 090328A (Ghisellini et al. 2010), GRB 090510 (Abdo et al. 2009b), GRB 090902B (Abdo et al. 2009c), and GRB 091003 (Ghisellini et al. 2010).

The 35 GRB fields with no previous claims for LAT detections were the target of this investigation. We recognized that the most convincing argument for a true LAT identification should rely on the statistical distributions for the matched filter weights in LAT fields with similar characteristics. To increase that number as much as possible, the LAT fields of view (FoVs) were segmented into 24 circular tiles embedded on a spherical surface as shown in Figure 1. Each tile subtends a cone with a half-angle of 105. This tiling scheme was applied to both the GRB and random field data sets to realize 697 and 8802 independent directions in space satisfying all the criteria described previously. Taking advantage of the fact that each field observation was blocked into a 250 s segment, the number of independent observations was multiplied by 5 by regarding each 50 s time slice as a separate sample. Thus, there are 3485 background measurements taken from LAT observations obtained simultaneously with the candidate GRB fields and an additional 44,010 samples taken under similar but not identical conditions. One particular concern for an analysis of this type is that false positives will selectively occur as the sample photon rate rises substantially above the mean. Evidence that this is not the case here is shown in Figure 2 which plots the cumulative distribution of the total number of photons within the LAT FoV over a 250 s interval. These rates explicitly exclude contamination from photons beyond the 105° zenith-angle cut. As shown in the plot, the two low-fluence GRB events reported here are not associated specifically with fields with high-ambient background rates. The similarity of the distributions for GRB and random fields also shows that the GRB data are not pathological as far as rates are concerned.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The statistical weighting scheme with the additional fractional sharing requirement clearly identified four of the six events previously associated with high-energy gamma-ray emission. The exceptions were GRB 090323 and GRB 0903028A for which photon counts were very sparse during the first 50 s following the burst trigger. Of the two, 090323 was marginally detected by the matched filter technique but GRB 0903028A provided no sign of a signal whatsoever. This gives some order of magnitude indication of the likely efficiency of the method. Three of the six were extremely intense with 100 or more energetic photons (GRB 080916C, GRB 090510, and GRB 090902B). As will be discussed later, the existence of such a substantial bright fraction is a rather surprising situation.

Most significantly, our analysis identified two additional bursts, GRB 080905A and GRB 091208B, as probably associated with energetic photon detections in the LAT. The best estimate for the probability of these occurring by chance alone was obtained by performing identical searches on random LAT fields with the same criteria. Thus, 1 out of 44,010 random fields generated matched filter weights exceeding the value for GRB 080905A and 11 similar fields generated weights exceeding the value for GRB 091208B. Multiplying by a trial factor of 35 for the number of localized fields considered, the associated probabilities are 8.0 × 10⁻⁴ and 8.7 × 10⁻³. To check that these correlations were simply not due to preferentially higher background rates for the GRB exposures, we also performed similar calculations for the LAT data confined to ∼20 uncorrelated directions and five independent time intervals from the same GRB data sets. All of these results are plotted as cumulative distributions in Figure 3. The statistical similarity of the GRB off-axis and random field data demonstrates that the GRB data set is not correlated with anomalous environmental conditions such as higher cosmic ray background rates. A list of photons associated with these two events is provided in Table 1.

Zoom In Zoom Out Reset image size

Two independent observations at significance levels of 0.1% and 1% already make this a reasonably credible result. However, as we discovered the importance of LAT "diffuse" photons for identifying point sources such as GRBs, we realized that independent of the matched filter weight analysis, a simple count of the number of such quanta was statistically quite significant. Thus, for the 35 well-localized GRB fields that do not include the previously reported LAT-detected GRBs, the total number of "diffuse" photons is 25 for a rate of 0.714 ± 0.143. For the GRB off-axis fields and the random fields, the comparable numbers are 1384 for 3485 fields and 17,414 for 44,010 to yield per field rates of 0.397 ± 0.011 and 0.396 ± 0.003. The random probability of such differences is less than 2.7%. As can be seen from the lists of photon counts in Tables 1 and 2, the two new GRB identifications are strongly correlated with these "diffuse" photons as well. Thus, the combined probability that these claims might be randomly generated must be less than 1 × 10⁻⁵.

GRB 080905A is a short burst with optical counterparts and a redshift of 0.128 (Rowlinson et al. 2010), while GRB 091208B is a long burst also with optical counterparts (Pagani et al. 2010) and a redshift of 1.063 (Wiersema et al. 2010; Perley et al. 2010). As listed in Table 1, the estimated number of >100 MeV photons for each event is about 3. This is substantially lower than any previous LAT GRB identifications. Estimating the corresponding fluence of >100 MeV photons raised some interesting statistical issues that are perhaps not well appreciated. To first order, one would anticipate simply adding up all the energies of the detected quanta and dividing by the effective detector area. Although this produces a number, it bears only a casual relation to the intrinsic properties of the GRB in question. The reason is that the spectral number distribution, dN/dE follows a power law with exponent between −3 and −2. Thus, the mathematical mean is defined but the variance is infinite. This leads to the conclusion that sums of energies drawn from the identical parent population will vary by factors of the order of unity, independent of the number of photons in the sample. In our case, the statistical uncertainty is already substantial since only three photons are detected for each event. Thus, we opted for a solution that invoked the least number of new assumptions: we adopted the LAT fluence values for bright bursts computed by Ghisellini et al. (2010) and listed in their Table 1. We next assumed that the spectral distribution for all bursts was an approximately universal function of energy so that fluences would scale with photon number. This led to estimates of the LAT fluence value of 3 × 10⁻⁶ erg cm⁻² for both events. Unfortunately, the proportionality of fluence to photon counts is not very evident. As can be seen in Table 2 for the two bursts, GRB 090510 and GRB 090902B, for essentially the same number of photons, the fluences differ by a factor of 16. This problem may require a better realization of how fluences should be estimated. For future purposes when dealing with events with more numerous quanta, there is a robust measure of fluence: multiply the number of photons by the median energy. In the continuum limit, a power-law distribution generates the following relation between E_mean and E_median over the interval from threshold to infinity:

$$\begin{equation} {E_{\rm mean}}=\frac{\Gamma - 1}{\Gamma - 2}\frac{1}{2^{\Gamma - 1}}E_{\rm median}. \end{equation} \tag{ 9 }$$

Thus, a robust proxy for gamma-ray burst fluence can be obtained by multiplying this estimate of E_mean by the number of detected photons.

The fact that >100 MeV photons have been observed with both long and short bursts suggests that the production mechanism is generic to relativistic jets, independent of the detailed structure of the progenitor system. So far, most of the high-energy photon observations have been associated with large fluences at the lower energies measured by the Swift/BAT and the Fermi/GBM detectors. To explore such fluence correlations, Figure 4 shows a scatterplot of LAT fluences versus GBM fluences for the six bright bursts previously identified by the LAT collaboration, the two bursts uncovered by the statistical analysis described in this paper and the co-added sum of nine bursts with no strong perceptible radiation above 100 MeV. For these purposes, GBM fluences were defined by photons in the range of 8–1000 keV within the T₉₀ time interval. The LAT fluences were obtained from Ghisellini et al. (2010) and the scaling method described previously.

Zoom In Zoom Out Reset image size

The salient feature of Figure 4 is the overall linear correlation of the LAT and GBM fluences with a ratio of unity between the two. The short burst identified by our analysis, GRB 080905A, lies above the trend line, a feature previously noted for this class of events (see the discussion of GRB 081024B in Abdo et al. 2010 and GRB 090510 in Ackermann et al. 2010). GRB 091208B is completely consistent with the higher fluence events. An interesting aspect of this graph is the estimate for the total co-added LAT fluence for nine bursts with GBM fluences roughly equal to or greater than GRB 091208B, but which lie a factor of 10 lower but with rather large errors. If one believes that all GRBs have a similar generic behavior, it is possible that a large fraction of these events have delayed high-energy emission that evades our estimates. It is more likely that there is a real population difference between events such as 080916C, 090926A, and 091208B and the majority of all other GRBs. Examining this question with larger data sets may shed some light on the mechanisms that drive high-energy emission processes.

Although the statistics are admittedly crude, we can use this analysis to set some limits on the fraction of high-energy photon GRBs at the present levels of sensitivity. Starting with a sample of 41 well-localized events, we have identified two with high statistical significance. The remaining number of LAT "source" and "diffuse" photons above background is no more than approximately three. The efficiency of this technique is about 50% for a total, so as many as 16 GRBs out of 41 fields could be associated with >100 MeV photons. Thus, instead of 17 LAT-detected GRBs over a two-year period, there might be 45 events or 19% of all viewable GBM triggers with a lower limit of 25 events or 11% of the total. That is lower than any estimates prior to the Fermi launch.

By using optimal signal detection techniques, we have demonstrated that the existence of GRB high-energy photon emission can be extended downward in brightness by a significant factor. Increasing the dynamic range over which such phenomena can be measured puts a number of constraints on physical models that can describe these events. As we have shown, these statistical methods can set interesting limits on the behavior of ensembles of events at lower levels than is possible for single isolated bursts, raising serious questions about the homogeneity of the GRB phenomenon. After another two years in orbit, the Fermi data will increase enough to allow more careful re-examination of the initial conclusions presented here.

The most surprising aspect of this work is the very small number of GRBs that could be positively identified with high-energy emission as a result of the substantially lower fluence thresholds. Almost any model would anticipate an increase in number proportional to the inverse of the ratio of detection thresholds raised to the power of unity or greater. Monte Carlo estimates indicate that our detection efficiency was of the order of 50% for the two events reported. This does not come close to explaining the apparent deficit of low-intensity events compared to the number of bursts similar to GRBs 080916C, 090510, and 090902B. We look forward to pursuing this problem during this precious and unique period when both Fermi and Swift are working magnificently.

We thank Chris Shrader, director of the Fermi Science Support Center for his considerable help in obtaining and interpreting the Fermi mission data products. Fang Yuan provided valuable assistance in the initial process of learning how to access and manipulate the Fermi data. This research is supported by NASA grant NNX08AV63G and NSF grant PHY-0801007.