FERMI LAT STACKING ANALYSIS OF SWIFT LOCALIZED GRBs

The Fermi Gamma-ray Space Telescope consists of two primary instruments, the GBM and the LAT. The GBM has 14 scintillation detectors that together view the entire unocculted sky. Triggering and localization are performed using 12 sodium iodide (NaI) detectors with different orientations placed around the spacecraft in four clusters of three. Two bismuth germanate (BGO) scintillators are placed on opposite sides of the spacecraft so that at least one detector can see any triggered event. GBM spectroscopy uses both the NaI and BGO detectors, sensitive between 8 keV and 1 MeV and between 150 keV and 40 MeV, respectively. Together they provide nearly 4 decades of energy for unprecedented sensitive spectroscopic studies of GRBs (Meegan et al. 2009).

The LAT is a pair conversion telescope comprising a 4 × 4 array of silicon strip trackers and cesium iodide (CsI) calorimeters covered by a segmented anti-coincidence detector to reject charged-particle background events. The LAT covers the energy range from 20 MeV to more than 300 GeV with an FOV of ∼2.4 sr. The deadtime of the LAT is nominally 26 μs, which is crucial for observations of high-intensity transient events such as GRBs. The LAT triggers on many more background events than celestial γ rays. Onboard background rejection is supplemented on the ground using event class selections that accommodate the broad range of sources of interest (Atwood et al. 2009).

The Swift spacecraft consists of the Burst Alert Telescope (BAT; Barthelmy et al. 2005), XRT (Roming et al. 2005), and Optical and Ultraviolet Telescope (UVOT; Roming et al. 2005). The BAT is a wide-field, solid-state γ-ray detector, covering an FOV of 1.4 sr and an imaging energy range of 15–150 keV. The instrument's coded mask allows for positional accuracy of 14 within seconds of the burst trigger. The XRT is a focusing XRT covering an energy range of 0.2–10 keV and providing a typical localization accuracy of ∼3''. The UVOT is a clear-aperture Ritchey-Chretien telescope that provides optical and ultraviolet photometry and subarcsecond positional accuracy of the long-lived afterglow counterparts to the prompt emission from GRBs.

The analysis presented here focuses on two distinct techniques with which to search for a signal in stacked data. The first method takes the sum of the observed counts collected over a specific duration, energy range, and region of interest (ROI) on the sky centered on each burst's best known position and compares it to the counts collected over an interval deemed to adequately represent the expected background at the time of trigger. The significance of the stacked counts compared to the stacked background is then estimated through Gaussian statistics.

The size of the ROI is typically chosen to reflect the 95% containment radius of the LAT energy-dependent PSF at a particular energy. For this analysis, we examine both a fixed 10° and 12° energy-independent ROI for each burst, as well as an energy-dependent ROI with a size ranging from 12° at 100 MeV to 22 at 10 GeV. The size of the energy-dependent ROI reflects the 95% containment of the LAT PSF when considering the P7REP_SOURCE_V15 instrument response functions. From these ROIs, we select "Source" class events that occurred within the first 2700 s after each GRB's trigger time (T₀) with an energy between 75 MeV and 30 GeV. The "Source" data class was specifically optimized for the study of point-like sources, with stricter cuts against non-photon-background contamination in comparison to the "Transient" data class that is typically used to study GRBs on very short timescales (Ackermann et al. 2012a). We also excluded events with an estimated angle from the Earth's zenith greater than 105° to guard against possible contamination from photons that originate from the Earth's limb. The selection of 75 MeV as the minimum energy is motivated by the possibility of detecting emission from bursts with significant attenuation emission below 100 MeV (e.g., Ackermann et al. 2013b). The choice of 2700 s reflects the longest amount of time that most sources on the sky can continuously stay in the LAT FOV before being occulted by the Earth (the exception being those located near the north orbital pole).

Each of our ROI selections has its own set of advantages and disadvantages. The energy-dependent ROI has the advantage of collecting photons within a radius that is commensurate with the LAT 95% containment radius at a given energy, potentially reducing background contamination at high energies, where the LAT PSF is narrowest. The 12° energy-independent ROI, on the other hand, makes no assumptions regarding the exact shape of the LAT PSF, at the cost of collecting more background at higher energies. The 10° energy-independent ROI is an intermediate solution, encompassing less of the PSF tail at low energies, but still collecting higher levels of background at high energies compared to the energy-dependent ROI. The use of a 10° ROI also reflects the better PSF of the Source class selection used in this analysis, compared to the Transient class selection, for which a 12° ROI is typically used.

The second technique we employ to search for a subthreshold signal in the stacked LAT data consists of a joint likelihood analysis using the analysis tools developed by the LAT team (ScienceTools version v9r30p1).⁶³ In standard unbinned likelihood fitting of individual sources, the expected distribution of counts for each burst is modeled as a point source using an energy-dependent LAT PSF and a power-law source spectrum with a normalization and photon index that are left to vary as free parameters. Galactic and isotropic background components are also included in this model. The Galactic component, gll_iem_v05, is a spatial and spectral template that accounts for interstellar diffuse γ-ray emission from the Milky Way. The normalization of the Galactic component is kept fixed during the fit. The isotropic component, iso_source_v05, provides a spectral template to account for all remaining isotropic emission that is not represented in the Galactic diffuse component and therefore accounts for contributions from both residual charged-particle backgrounds and the isotropic time-averaged celestial γ-ray emission. The normalization of the isotropic component is allowed to vary during the fit. Both the Galactic and isotropic templates are publicly available.⁶⁴

We can then derive the probability density of the observed data given this model and create a "likelihood" function by treating this probability density as a function of the model parameters. The likelihood function is essentially the probability of observing the data given the chosen model for a range of model parameters. By maximizing the likelihood function with respect to our parameters of interest, we can then estimate the model parameters that make the observed data the most probable.

The maximum likelihood technique can be expanded to apply to multiple sources to constrain or estimate a set of parameters thought to be common to those sources. In this case, the backgrounds and properties of each source can be modeled individually in source-specific likelihoods, as described above. The source-specific likelihoods can then be multiplied to yield a joint likelihood function that is used for inference on the common parameter of interest. In this way, the joint likelihood becomes a source stacking technique, which has been previously applied to LAT data in the search for dark matter (Ackermann et al. 2011) and to study the extragalactic background light (Ackermann et al. 2012b). This method is implemented in the Fermi Science Tools as the Composite2 routine.

To quantify the significance of a potential excess above background, we employ the common likelihood-ratio test (Neyman & Pearson 1928). For this test, we form a test statistic (TS) that is the ratio of the likelihood evaluated at the best-fit parameters under a background-only, null hypothesis, i.e., a model that does not include a point-source component, to the likelihood evaluated at the best-fit model parameters when including a candidate point source at the center of the ROI. According to Wilks's theorem, this ratio is distributed approximately as χ², so we choose to reject the null hypothesis when the TS is greater than TS = 25, roughly equivalent to a 5σ rejection criterion.

The data selection for this method is identical to that performed for the counting analysis in the case of the 10° and 12° ROIs described in Section 4.1. We note that the Fermi Science Tools do not currently include a means of performing a likelihood analysis on data selected using an energy-dependent ROI, as is possible with the counting analysis.

The proper selection of background regions with which to compare the stacked signal is a crucial component to the counting analysis and provides a control sample for the joint likelihood analysis. Since the spacecraft's geomagnetic coordinates and its celestial pointing can vary significantly over the 2700 s duration under consideration, the use of an off-and-on source method of background determination may not always adequately represent the true background. In addition, the use of intervals immediately prior to the burst trigger or after the end of the prompt emission, as observed by GBM or BAT, to estimate a background interval could serve to bias our investigation, as LAT-detected GRBs have exhibited both delayed and extended high-energy emission on timescales that exceed the durations traditionally defined by observations in the keV−MeV energy range (Ackermann et al. 2013b) and subthreshold emission prior to the burst trigger cannot be ruled out. Therefore, we avoid using intervals immediately preceding or succeeding the burst activity for background estimation.

Instead, we attempt to locate an interval at least one orbit prior to the GRB trigger that best matches the spacecraft's observing conditions at the time of the trigger. This includes matching the same off-axis angle between the GRB sky coordinates and the LAT boresight, the spacecraft's geomagnetic coordinates in orbit, and the angle of the GRB location to the Earth's zenith at the time of trigger. These three criteria serve to match the charged-particle background, the Galactic and isotropic backgrounds, and possible contamination from Earth limb photons, respectively, to that observed by the LAT at the time of the GRB trigger. A period of 30 sidereal orbits (171,915 s) prior to trigger adequately matches these geomagnetic and pointing criteria and provides an interval with which to estimate the background during the GRB. An example of the spacecraft's orbit and orientation with respect to the position of a GRB in our sample and its associated background orbit can be seen in Figure 2.

Zoom In Zoom Out Reset image size

We note that this selection does not adequately represent the expected background for the 15 GRBs that either triggered ARRs of the spacecraft or occurred while the spacecraft was performing a ToO, since maneuvers take the spacecraft out of survey mode and initiate a custom pointing that keeps the GRB in the LAT FOV. As a result, there is no previous interval that matches the spacecraft's orbit and orientation during the repoint, and as such, these bursts are excluded from the counting analysis. Since the background for each individual GRB is being modeled in the likelihood analysis, this method can adequately take into account the spacecraft motion for these bursts. This leaves a total of 64 bursts for which we can apply the counting analysis. Therefore, the application of the likelihood technique to the bursts for which we have the most comprehensive observations is a significant advantage of this method over the counting analysis.

In order to validate the effectiveness of this background selection for bursts that did not trigger an ARR or occurred during a ToO, we extract the observed photons over a 2700 s interval from 1000 random locations on the sky covering an energy range from 75 MeV to 300 GeV, using a 10°, 12°, and energy-dependent ROI. We then extract the observed photons at the same location, but 30 sidereal orbits prior to the original selection interval. We can then test whether the observed photons during the signal and background intervals are consistent with being drawn from the same distribution by examining the resulting significance distribution. We define the counting significance as $(S-B)/\sqrt{S+B},$ where S and B are the counts in the signal and background intervals, respectively (Li & Ma 1983).

The resulting distribution for the 12° ROI has a mean value of μ₁₂ = 0.09 and variance of σ₁₂ = 1.01. Among our 1000 trials, we obtain three false positives above 3σ, consistent with expectations, as we expect 99.7% of values drawn from our significance distribution to lie within 3σ of the mean. The results for the 10° and energy-dependent ROIs were consistent with these values. Therefore, our background selection is found to be robust for any single source and background interval selection.

Using the data selection criteria described in Section 4.1, we can directly compare the extracted signal to the counts accumulated using the same selection criteria, but collected during the background orbits described in Section 4.3. When considering the entire interval from T₀ to 2700 s and an energy range from 75 MeV to 30 GeV, we obtain a total accumulated signal and background of S₁₀ = 1711 and B₁₀ = 1476, S₁₂ = 2365 and B₁₂ = 2161, and S_EROI = 750 and B_EROI = 626 counts for a significance of 4.16σ, 3.03σ, and 3.34σ for the 10°, 12°, and energy-dependent ROIs, respectively.

Introducing cuts on the data, as we have here with the extraction radius, in order to maximize the observed significance introduces a bias in the analysis due to the so-called look-elsewhere effect. The chance that the observed significance could have arisen at random owing to the size of the parameter space that was searched can be accounted for by applying trial factor corrections to the final significance. Here we searched the data using three different extraction radii, and therefore the final significance needs to be attenuated as 1−(1 − α)^N, where α is the probability of observing such a value by chance⁶⁵ of the detection and N represents the number of trials. For the analysis presented above, our 4.16σ detection therefore becomes an erf ${\left(\tfrac{4.16}{\sqrt{2}}\right)}^{3}$ ∼ 3.9σ detection. We note that this correction factor assumes that each data set is statistically independent, which is overly conservative in our case, since the three different ROI selections result in data sets that are subsets of each other. Therefore, we regard the 3.9σ detection significance as a conservative lower limit to the true significance of the signal over the background.

Focusing on the 10° ROI analysis, we can create a stacked light curve that is co-aligned to the trigger time of each GRB. This light curve of stacked signal and background counts, binned to 100 s intervals, is shown in the top panel of Figure 3, with the middle panel showing the Gaussian significance of their difference and the bottom panel showing the summed LAT effective area during the observations. The overall shape of the light curve reflects the evolution in the total effective area. The observed count rate rises as additional fields containing GRBs enter the LAT FOV and falls as they exit or as the spacecraft enters the SAA, whereby data taking is disabled. The total effective area, and hence the stacked counts light curve, peaks near the co-aligned trigger, as all fields are predicated to be in the LAT FOV at this time owing to our sample selection criteria. The stacked light curve shows a ∼3σ excess roughly 30 s after the co-aligned trigger, followed by additional periods of excess signal over background, although a consideration of trial factors lessens the significance of any one peak in such a time-resolved analysis.

Zoom In Zoom Out Reset image size

The median ratio of the exposure during the signal and background orbit is found to be μ = 1.05, with a variance of σ = 0.072, highlighting that the exposures of the signal and background orbits are well matched, but not exactly equal. Small mismatches in exposure such as these may be a fundamental limitation of the stacked counting analysis performed here. Such differences are small for the source/background comparison for any one burst, but when adding 79 such comparisons, these small mismatches contribute to a non-negligible difference in the total summed exposure, which can be seen in the bottom panel of Figure 3. Ultimately, normalizing the observed counts by the estimated exposure for each burst and performing a summed rate comparison can account for these differences.

Next, we examine the difference between the cumulative signal and background over much wider time intervals. In the top panel of Figure 4, we show the cumulative signal and background counts as a function of time since T₀ – 1000 s for a 10° ROI selection, with the bottom panel again showing the significance of their difference. The effective area of the stacked observations drops to zero above T₀ + 2700 s as all fields exit the LAT FOV, and as such, the cumulative light curve levels off accordingly. The excesses seen in Figure 3 can be clearly traced as local maxima in Figure 4. Although the difference between the signal and background varies above and below zero significance prior to the co-aligned trigger, the signal clearly begins to diverge from the background at T₀, above which the significance climbs to approximately ∼3σ at T₀ + 2700 s. We note that this cumulatative significance differs from that quoted above because the integration period here begins at T₀ – 1000 s.

Zoom In Zoom Out Reset image size

We also examine the dependence of the signal significance on the minimum energy E_min used in our selection criteria. In Figure 5, we plot the resulting signal significance as a function of E_min when considering the interval from T₀ to 2700 s for a 10° ROI selection, showing that the detection significance does not improve when only considering higher-energy photons. The dashed line, representing the number of bursts with photons contributing to the modified selection criteria, falls steadily with increasing E_min.

Zoom In Zoom Out Reset image size

Finally, an identical analysis of the preburst observations, covering an interval from T₀ – 2700 s to the co-aligned trigger, reveals no excess emission above the background. For this control sample, using a 10° ROI selection, we obtain a total accumulated signal and background of S = 1704 and B = 1681 counts, respectively, for a significance of 0.40σ.

Using the data selection criteria described in Section 4.1 and the analysis method outlined in Section 4.2, we first performed a joint likelihood analysis on the ensemble of GRB locations, including those that triggered ARRs of the spacecraft or occurred during ToOs, for an integration period covering the entire 2700 s post trigger and an energy range from 75 MeV to 30 GeV, using both a 10° and 12° ROI. The resulting TS of the joint likelihood fit is 59 and 58.6, respectively, roughly representing a ∼7σ detection significance of a point source being present in excess to the expected background using either ROI size. The TS of an identical analysis performed on the background orbits defined in Section 4.3 is consistent with zero.

Focusing on the 10° ROI analysis, we examine the dependence of the TS value on the minimum energy E_min used in our selection criteria. In Figure 6, we plot the resulting TS for a range of E_min values, showing that the detection significance rises as low-energy photons are excluded from the fit, peaking at E_min = 300 MeV, before eventually leveling out above 400 MeV. The dashed line represents the number of analysis intervals for which there we observed photons matching the selection criteria. As E_min rises, the number of intervals with photons contributing to the joint likelihood analysis falls, with nearly one-third of the burst positions being removed from the sample when E_min = 1000 MeV. Because the detection significance peaks at E_min = 300 MeV, we will focus the remainder of our joint likelihood analysis on an energy selection criterion covering 300 MeV to 30 GeV.

Zoom In Zoom Out Reset image size

In order to investigate the time dependence of the excess signal, we also performed a joint likelihood analysis over a range of integration times before and after the trigger time, going from from T₀ – 1000 s to T₀ and from T₀ to T₀ + 2700 s, in 50 s intervals. The TS as a function of integration time is shown in Figure 7. The red dotted line represents the TS of an identical analysis performed on each burst's associated background orbits. The TS of the joint likelihood fit to the data collected over the 1000 s prior to the co-aligned trigger is consistent with zero. The TS rises slightly for shorter integration duration approaching T₀, but is below TS = 9, roughly equivalent to 3σ, for all integration durations prior to T₀, except for the 50 s interval covering T₀ – 50 to T₀, at which point it approaches TS ∼ 10. At all subsequent integration times after T₀, the significance of the signal excess above background is greater than 3σ, reaching 5σ within 200 s after T₀. The TS values exhibit local maxima as a function of time, reflecting the arrival of photons in excess to our background model, before leveling off at TS = 59 as the accumulated effective area drops to zero. The data collected from the same fields during the background intervals show no such excess when fit to our background models, being consistent with TS ∼ 0 for roughly all integration periods under consideration.

Zoom In Zoom Out Reset image size

Before we can assess the final significance of our likelihood analysis, we again need to take trial factors into account. Since we considered two different ROI sizes and 10 different minimum energy boundaries, this introduces N = 12 trials. Therefore, a TS = 59, or $\sqrt{59}\sim 7.68\sigma ,$ detection becomes an erf ${\left(\tfrac{7.68}{\sqrt{2}}\right)}^{12}$ ∼ 7.46σ detection. Again, since the ROI and minimum boundary selections did not result in statistically independent data sets, we regard this attenuated significance as a conservative lower limit to the true significance of the detection.

By scanning the joint likelihood function over a range of photon indices and flux normalizations for the point source in our model, we can obtain a joint likelihood profile that is a function of these two common parameters of interest. Finding the maximum of this profile allows us to estimate the stacked source flux and characteristic photon index of the sample. A contour plot of the resulting joint likelihood profile for the analysis covering the entire 2700 s interval post trigger and an energy range from 300 MeV to 30 GeV is shown in Figure 8. The best-fit photon flux and photon index of the combined data are F_ph = 6.4 × 10⁻⁸ photons cm⁻² s⁻¹ and Γ = −1.92, respectively. The 68%, 95%, and 99% confidence levels (C.L.) for the photon flux estimate are shown as the solid, dashed, and dot-dashed lines, respectively.

Zoom In Zoom Out Reset image size

For comparison, a typical photon flux upper limit for an individual well-observed GRB that remained in the LAT FOV for an entire 2700 s duration is on the order of F_ph,UL ∼ 10⁻⁶ (Ackermann et al. 2012b), making it clear why none of the bursts in our sample were individually detected. The best-fit photon index is representative of indices measured at late times in previously detected LAT bursts (e.g., GRB 110731A) and consistent with photon indices measured by XRT of GRB afterglows at late times.

In order to ensure that one or more GRBs are not dominating the observed excess seen in Figure 7, we examine the TS distribution for the 2700 s integration period succeeding T₀ in Figure 9. The distribution peaks at TS = 0 for 41 GRBs, nearly half the sample under consideration. The remaining bursts are evenly distributed between 1 ≲ TS ≲ 10, with one burst (GRB 110903A) at TS = 23. This burst is just under the fiducial 5σ threshold for a burst to be considered detected by the LAT. Removing this burst from the sample, we find that the TS of the joint likelihood analysis for the 2700 s second integration interval drops from TS = 55 to TS = 41, still yielding a strong detection of the remaining bursts. Moreover, the best-fit flux and index remain relatively unchanged, at F_ph = 6.0 × 10⁻⁸ photons cm⁻² s⁻¹ and Γ = −2.09.

Zoom In Zoom Out Reset image size

The counting and joint likelihood analyses employed here have both revealed the existence of a signal in excess of the expected background for the GRB locations in our sample. Here we compare the two methods to determine whether there are any inherent benefits to using one method over the other. The significance of the detected signal as a function of integration times, covering an energy range from 300 MeV to 30 GeV using a 10° ROI, for each method is displayed in Figure 10. In order to ensure a proper comparison, we have removed the bursts that triggered ARRs of the spacecraft or occurred during ToOs from the likelihood analysis.

Zoom In Zoom Out Reset image size

With the ARR and ToO bursts removed, the joint likelihood method provides stronger evidence for a signal excess when considering the entire 2700 s integration period, yielding a total signal significance of ∼5σ, whereas the counting technique provides only a marginal detection at ∼3.4σ. On smaller integration time, the two methods are roughly consistent, before diverging for longer integration intervals. The significance derived from the joint likelihood analysis rises gradually as more data are accumulated before leveling out to roughly 5σ, whereas the significance derived from the counting analysis remains largely unchanged through the inclusion of these data. Note that the final signal significance quoted here for the counting analysis differs from that reported in Section 5.1 because of the different energy ranges under consideration. Here we are examining photons over a 300 MeV to 30 GeV energy range, in order to match the analysis performed using the joint likelihood technique.

In order to examine which bursts in our sample contribute most significantly to the observed signal, we calculate the cumulative TS by rerunning the likelihood analysis for each additional burst, sorted as a function of a selected burst property. An example of this method can be seen in Figure 11, where the cumulative TS is plotted as a function of burst exposure. Here exposure is defined as an integral of the total response over the entire ROI.⁶⁶ The TS is first calculated for the burst with the lowest exposure, covering an energy range from 300 MeV to 30 GeV and a period of 2700 s post trigger, then recalculated by including the data from the burst with the next highest exposure, and so on. The result is an increasing cumulative TS that eventually peaks at the value reported in Section 5.2. The dashed line in Figure 11 represents the cumulative sample size as a function of exposure, and the dotted line demarcates the 50th percentile of the sample. If all the bursts contributed equally to the final TS value, independent of their exposure, then one would expect the cumulative TS to rise at the same rate as the cumulative sample size. This is not the case in Figure 11, where the bursts with exposures in the first 50th percentile of the sample contribute very little to the final signal significance. This would be expected, since longer observations of bursts closer to the instrument boresight would be more sensitive at detecting extended subthreshold emission. Still, Figure 11 shows that exposure alone is not entirely indicative of whether a burst will contribute to the final TS value. Several bursts with the highest exposure are seen to contribute very little to the overall signal, whereas several of the bursts with the lowest exposure do contribute to the significance of the final signal. This is also reflected by the fact that the non-ARR sample of bursts discussed in Section 5.3 is still detected at roughly ∼5σ (TS ∼ 26), compared to ∼7.7σ (TS ∼ 59) for the entire sample. Therefore, the ARR and ToO bursts contribute most, but not all, of the observed signal.

Zoom In Zoom Out Reset image size

We can test whether the observed low-energy γ-ray flux of a source plays an additional role in its detectability and hence contribution to the final TS value. To examine this, we calculated the cumulative TS as a function of the burst's peak photon flux in the 15–150 keV energy range as measured by BAT (Donato et al. 2012), and we show the results in Figure 12. The bursts with the highest observed peak flux contribute strongly to the cumulative TS, although there remain bursts with flux values above 1.0 photon cm⁻² s⁻¹ that do not appear to contribute significantly to the final TS value. Note that the sample size presented in Figure 11 differs from that in Figure 12 because not every burst in the likelihood sample was initially detected by BAT, and therefore a photon flux estimate of their prompt emission is not available.

Zoom In Zoom Out Reset image size

Given the extended nature of the signal inferred in Figure 10, we also calculate the cumulative TS as a function of the burst's X-ray brightness at 11 hr, in the 0.3–10 keV energy range, as measured by XRT. The X-ray flux at 11 hr has become a standard measure of afterglow brightness, as it samples the afterglow light curve at a period where the steep decay and plateau phases have typically ceased (Nousek et al. 2006). To calculate these values, we downloaded the XRT data for each from the Swift XRT light-curve repository (Evans et al. 2009) and fit the flux light curves with a multi-segmented afterglow model defined by Zhang et al. (2006). The fitting procedure employed is outlined in Racusin et al. (2009). The result of the X-ray brightness versus cumulative TS analysis is presented in Figure 13. There exists an even stronger trend of bursts with bright X-ray emission 11 hr post trigger contributing significantly to the final TS value, compared to the prompt γ-ray flux measured by the BAT. In fact, fewer than half of the brightest bursts in the sample contribute most of the signal, with the addition of the remaining bursts actually decreasing the final signal significance.

Zoom In Zoom Out Reset image size