THE FIRST FERMI-LAT GAMMA-RAY BURST CATALOG

We select Transient class data with reconstructed energies in the 100 MeV–100 GeV range. The lower limit is chosen to reject events with poorly reconstructed directions and energies. Moreover, for Pass 6, the LAT response is not adequately verified at E < 100 MeV energies and the contamination from cosmic rays misclassified as gamma-rays is also significantly increased. The upper limit (UL) was chosen at 100 GeV since we do not expect to detect GRB photons at such high energies due to the opacity of the universe and the limited effective area of the LAT. We select events in a circular region of interest (ROI) that is centered on the best available GRB localization. The LAT PSF depends on the event energy and off-axis angle and has been studied using Monte Carlo simulations. We use the resulting description of the PSF to increase the sensitivity of our analyses. For the event-counting and joint spectral-fitting analyses, we select a variable ROI radius that depends on the event energy and the off-axis angle of the GRB in such a way as to select almost all the events consistent with the position of the GRB given our PSF while rejecting much of the residual cosmic-ray background, increasing the signal-to-noise ratio of the selected data. To accomplish this, we split the events in logarithmically spaced bins in energy and for each bin we select only the events contained in an ROI around the source having a radius corresponding to the 95% containment radius of the PSF evaluated at an energy equal to the geometric mean of the bin's energy range. For the duration estimation using Transient data, we deal with longer time periods, thus we dynamically adjust the radii of the energy-dependent ROIs to follow the variation of the off-axis angle with time. On the other hand, for the LLE duration estimations and the joint GBM-LAT spectral analyses, we use a single set of radii calculated using the PSF corresponding to the GRB off-axis angle at trigger time. The exact dependence of the LLE PSF on the off-axis angle is not available yet. Instead, only two possible LLE PSFs are available for setting the ROI radii: one for observations with off-axis angles greater than 40° and the other for observations closer to the center of the FoV. Finally, for cases for which the GRB localization error is non-negligible (i.e., for GBM or LAT localizations), we increase the radius of each ROI by setting it equal to the sum in quadrature of the localization error and the 95% containment radius of the PSF. For GRBs localized by the Fermi GBM, we also added in quadrature a 3° systematic error. The maximum-likelihood analysis utilizes the PSF information internally while calculating the probability of each event being associated with the GRB, thus no optimization of the ROI radius, as above, is necessary. For the maximum-likelihood analyses, we use a fixed-radius ROI set at 12°, a value larger than the 99% containment radius of the Transient LAT PSF evaluated for a 100 MeV event on axis.

We apply a cut to limit the contamination from γ-rays produced by interactions of cosmic rays with the Earth's upper atmosphere. For our maximum-likelihood analysis, we use the Fermi Science Tool gtmktime⁷⁴ to select only the time intervals (the "Good Time Intervals" or GTIs) in which no portion of the ROI is too close to the Earth's limb. Because the Earth's limb lies at a zenith angle of 113° and we wish to take into account the finite angular resolution of the detector, we exclude any events taken when the ROI is closer than 8° to the Earth's limb or equivalently when it intersects the fiducial line at 105° from the local zenith. For special cases, when the position of the GRB is very close to the Earth's limb, we compensate for the loss of exposure due to this cut by reducing the size of the ROI and simultaneously increasing the maximum zenith angle to 110°. This increases the duration of the GTI significantly, allowing deeper exposures for searches of late γ-ray activity. For all the other analyses (namely, event-counting analyses and joint spectral fitting), we do not apply a cut to select GTIs as above, but rather we process the whole observation and instead reject individual events reconstructed farther than 105° from the local zenith.

The response of a GBM detector depends on the continuously varying position of the GRB in its FoV, with its effective area decreasing as the angular distance between the detector boresight and the source (θ_GBM) increases. Because of this, when θ_GBM is large, any systematic effects due to imperfect modeling of the spacecraft or the individual detectors become relatively important (Goldstein et al. 2012). For this reason, we use the data from the GBM NaI detectors that have angles θ_GBM < 50° at the time of the trigger and the BGO detector facing the GRB at the time of the trigger.

We also exclude any detector occulted by other detectors or the spacecraft during any part of the analyzed time interval, as advised in Goldstein et al. (2012).

Since θ_GBM usually changes with time, the GBM Collaboration released RSP2 files that contain several response matrices corresponding to short consecutive time intervals (every 2° of slew of the detector about the source). With a suitable weighting scheme, as described in Section 3.4.1, these files provide an adequate description of the GRB detector responses.

Finally, in some cases, bright GRBs trigger an ARR, causing rapid variations of θ_GBM with time for some of the GBM detectors. These variations create further variations in those detector responses and background rates. In fact, due to its orbital and angular dependence, the background of those detectors can be very hard to predict. Also, since the RSP2 files might not be binned finely enough in time to cover these rapid variations, we excluded data from detectors that have such rapid variations.

The background in the LAT data is composed of charged cosmic rays (CRs) misclassified as γ-rays, astrophysical γ-rays coming from Galactic and extragalactic diffuse and point sources, and γ-rays from the Earth's limb produced by interactions of CRs in the upper atmosphere. The backgrounds for the Transient-class and LLE data are dominated by the CR component, while for the cleaner Diffuse class the backgrounds are dominated by astrophysical γ-rays. The CR component of the background depends primarily on the geomagnetic coordinates of the spacecraft and on the direction of the GRB in instrument coordinates (since the LAT's effective area varies strongly with inclination angle). The component from the Earth's atmosphere depends on the angle between the GRB and the limb (i.e., on the zenith angle of the GRB) and is strongest toward the limb. Finally, the astrophysical background γ-ray component depends on the GRB direction and is typically stronger at low Galactic latitudes.

For the Transient event class analyses, we use the Background Estimation tool (Vasileiou 2013; "BKGE," hereafter), which was developed by the LAT collaboration and takes into account all of these dependencies. It can estimate the total expected backgrounds for any given ROI and period of time with an accuracy of ∼10–15% (Abdo et al. 2009e). It also provides separate estimates for the Galactic diffuse emission and for everything else, namely the sum of CRs and extragalactic diffuse emission (the "isotropic component"). Note that the BKGE cannot estimate the backgrounds from the Earth's limb. However, the zenith angle cut described in Section 2 is very effective at reducing this component to negligible levels, thus this limitation does not generally constitute an obstacle.

Our maximum likelihood analysis of Transient-class data uses a background model calculated by a combination of the isotropic component provided by the BKGE tool and the Galactic diffuse emission template provided by the LAT Collaboration.⁷⁶

The maximum likelihood analysis using the cleaner Diffuse class data, which was performed for validation studies (see Appendix A), uses the Galactic diffuse emission template plus the public template describing the "isotropic background" (extragalactic diffuse emission and CR background) as a single spectrum of the intensity averaged over the whole sky. The BKGE does not produce estimates for Diffuse-class events. For the time scales analyzed in this study, the contribution from point sources is typically negligible, so we do not take them into account in the background models.

For the joint GBM-LAT spectral analysis, we used the estimates provided by BKGE of the total background in the energy-dependent ROI for the background for the LAT. For technical reasons related to the broad PSF of the LLE class, we cannot use the BKGE to estimate the LLE background. Instead, we evaluate it directly from the LLE data associated with each individual observation. First, in order to ensure enough events in every time bin, we bin the LLE data in time with a coarse binning of 5 s, from well before the trigger time to well after the end of the burst, as measured by the GBM. We then fit the background rate as a function of time b(t) by taking into account the variation of the exposure due to the changing orientation of the LAT. Phenomenologically, we adopt the function b(t) = p₀ + p₁C(t) + p₂C(t)², where C(t) = cos [θ(t)] and θ is the off-axis angle. The parameters p₀, p₁, and p₂ are obtained by fitting the "pre-burst" and "post-burst" time windows simultaneously. We use a conservative definition of these time windows based on the burst duration as measured by the GBM. In particular, the "post-burst" data start well after the end of the low-energy emission as seen by the GBM. Finally, the fit parameters allow us to compute the background rate at any time during the burst and we use the covariance matrix from the fit to evaluate the uncertainty of this prediction. We compared this simple model with an alternative prescription b(t) = pol(t)*C(t), where the degree of the polynomial function pol(t) is increased until a good fit to the data is obtained. Typically, a polynomial of degree 1 or 2 was sufficient, although in a few cases a higher degree (3 or 4) was necessary. The expression above is motivated by the fact that, as a first approximation, the effective area of the LAT to CRs scales as cos (θ) and that we can model the CR contribution on-axis (θ = 0) with a polynomial. The two prescriptions gave very similar results in all cases. An example of the standard prescription is shown in Figure 1. The excess visible at ∼200 s after the trigger time T₀ is significantly above the estimation of the background and we claim this is due to the GRB emission. Note also that in this case, the background estimation was particularly difficult due to the contamination of the Earth's limb caused by the ARR while in several other cases the background is smoother.

Zoom In Zoom Out Reset image size

We use the GBM CSPEC event data from before and after the GRB prompt phase to obtain a model for the background, similar to the procedure followed for the LLE data above. For each selected detector, we integrate the CSPEC spectra over all the energy channels to obtain a light curve and then select two off-pulse time intervals: one before and one after the GRB prompt emission (see the left panel in Figure 2). We fit polynomial functions f(t) of increasing degree D to the data from these two time intervals, minimizing the χ² statistic, until we reach a good fit (i.e., with a reduced χ² ≃ 1). Then, we consider the light curves corresponding to each of the 128 channels separately, again with data from the off-pulse intervals, and we fit them with a polynomial of degree D by minimizing the Poisson log-likelihood function.⁷⁷ After each fit, we check by eye that the residuals are consistent with the statistical fluctuations. If this is not the case, we repeat the procedure from the beginning, changing our choice for the off-pulse intervals, until a good fit is achieved. The set of 128 polynomial functions constitutes our background model and the predicted number of background events b_i in the ith channel of the background spectrum is the integral of the corresponding polynomial function f_i (describing the rate) between t₁ and t₂:

$$\begin{equation} b_{i} = \frac{\int _{t_{1}}^{t_{2}}f_{i}(t)dt}{t_{2}-t_{1}}. \end{equation} \tag{ 1 }$$

The statistical error of the integral is computed using the covariance matrix from the fit.⁷⁸ Since the background for GBM detectors is much less predictable than for LLE data, we determine the off-pulse regions manually. In order to minimize the statistical and systematic errors (and hence ensure a reliable background estimate), the off-pulse time intervals must be close to the GRB's signal, have a long enough duration, and also possibly have a smooth part of the light curve without bumps or other structures. Moreover, the number of counts in each channel is much smaller than the total number of counts used to determine D. Thus, the larger the value of D, the more f_i can pick up statistical fluctuations in some channels, giving a slightly wrong interpolation for those channels in the pulse region. Thus, we try to find off-pulse intervals well described by low-order polynomials (ideally, D = 1). Unfortunately, this is not always possible. For example, for GRBs triggering ARRs, the background can vary quickly in response to the change of pointing, requiring higher order polynomials to describe it. This effect introduces some additional noise in the spectrum, but it is unlikely to introduce any bias in the fit results, given its random nature. Note that it is not possible to fix the shape of the polynomial, since the background shows spectral evolution and thus every channel needs to be considered independently. In some cases, even with high-order polynomials, fitting the model to the background can be difficult and even impossible without being completely arbitrary (see the right panel in Figure 2 for an example). In those cases, we opt for excluding the problematic detector from the analysis. These issues are not solvable at present given our current understanding of the detectors and their backgrounds. More advanced techniques to deal with the backgrounds are currently under investigation by the Fermi-GBM Collaboration (Fitzpatrick et al. 2012).

Zoom In Zoom Out Reset image size

To determine the significance of the detections of sources using the maximum likelihood analysis, we consider the "Test Statistic" (TS) equal to twice the logarithm of the ratio of the maximum likelihood value produced with a model including the GRB over the maximum likelihood value of the null hypothesis, i.e., a model that does not include the GRB. The probability distribution function (PDF) of the TS under the null hypothesis represents the probability that a measured signal is consistent with the statistical fluctuations. The PDF in such a source-over-background model cannot, in general, be described by the usual asymptotic distributions expected from Wilks' theorem (Wilks 1938; Protassov et al. 2002). However, it has been verified by dedicated Monte Carlo simulations (Mattox et al. 1996) that the cumulative PDF of the TS in the null hypothesis (i.e., the integral of the TS PDF from some TS value to infinity) is approximately equal to a $\chi ^2_{n_{{\rm dof}}}/2$ distribution, where n_dof is the number of degrees of freedom (dof) associated with the GRB. The factor of 1/2 in front of the TS PDF formula results from allowing only positive source fluxes.

Since we model the GRB spectrum as a PL with two dof and we fix the localization, the TS distribution should nominally follow $(1/2)\chi ^2_{2}$. This is formally correct if the localization of the GRB is provided by an independent dataset (i.e., from another instrument). However, when the input localization is not sufficiently precise, we optimize it using the same dataset used for detecting the source, thereby introducing two additional free parameters (R.A. and Decl.). In this case, the TS distribution should follow $(1/2)\chi ^2_{4}$. In practice, the steps of detection and localization are iterated many times and a detection step is performed using an ROI centered on the position found by a prior localization step. Therefore, the datasets used in each step are not exactly overlapping. For this reason, we expect some deviation from $(1/2)\chi ^2_4$ distribution. For simplicity, we set a unique threshold of TS_min = 20 for our analysis independent of the origin of the localization. This formally corresponds to two slightly different one-sided Gaussian equivalent thresholds, 4.1σ for $\chi ^2_{2}$ and 3.5σ for $\chi ^2_{4}$. Additionally, we check the calibration of the detection algorithm on a sample of "fake GBM triggers" generated as described in Section 2.2. With the aforementioned value of TS_min, we obtain zero false detections on the "fake GBM triggers" sample (see Section 4.1 for more details).

We compute the localizations with the LAT in two steps. The first step provides a coarse estimation of the GRB position and is performed using the Fermi ScienceTool gtfindsrc. At this stage, we look for an excess consistent with the LAT PSF and we do not assume a particular background model. Although this method is quick and robust, it assumes that the likelihood function is parabolic and symmetric in azimuth around the found position and so the provided localization error can be slightly underestimated. Therefore, this step is only used to obtain an initial seed for the follow-up analysis.

For a more accurate localization, we use the Fermi ScienceTool gttsmap, which starts from the best-fit background model obtained by the likelihood fit and builds a map of the TS in a grid around the best available localization of the source. The GRB spectral parameters are fit at each position in the grid, along with all free parameters of the background model. The grid size and spacing are set based on the localization error obtained in the first step. The final LAT localization corresponds to the position of the maximum of the TS map. Its statistical uncertainty is derived by examining the distribution of TS values around it. Following Mattox et al. (1996), we interpret changes in the TS values in terms of a χ² distribution with two dof to account for the flux and spectral index of the GRB. Specifically, a confidence-level (CL) uncertainty is given by the TS map contour that corresponds to a decrease from the maximum value by a value equal to the CL quantile of the $\chi ^2_2$ distribution. For example, the 90% (68%) CL corresponds to a decrease of the TS from its maximum value by 4.61 (2.32).

We estimate the probability of each γ-ray being associated with the GRB by using the Fermi ScienceTool gtsrcprob. The probability computation takes into account the spectral, spatial (extent), and temporal (flux) information of all the components in the source model and the response of the LAT (PSF and effective area) to the particular event. The probabilities are assigned via likelihood analysis and are computed starting from the best-fit model. In particular, the probability that a photon is produced by a component i is proportional to M_i, given by

$$\begin{equation} M_{i}(\epsilon ^{\prime },p^{\prime }, t)=\int d\epsilon dp\,S_{i}(\epsilon ,p,t)\,R(\epsilon ,p;\epsilon ^{\prime },p^{\prime },t), \end{equation} \tag{ 2 }$$

where S_i(, p, t) is the predicted counts density from the component at energy , position p, and (observed) time t, and the integral is the convolution over the instrument response R(, p; ', p', t). In general, the predicted count density is the sum of the different contributions S_i(, p, t), including the extended backgrounds (such as the isotropic component and the Galactic diffuse emission), background point sources (nearby bright sources), and the GRB under study. Each contribution is described by a model, the parameters of which are optimized during the maximum likelihood fit. We simplify the calculation by not including nearby bright sources since, in these short time scales, they do not contribute significantly to the total number of counts. Once we compute the maximum likelihood model for the observed number of counts, we assign to each event the probability of it being associated with a particular component i.

Because the flux varies with time, we perform the calculation in several time bins so that the flux is never averaged over long time intervals. We tested schemes for defining the time intervals including linear, logarithmic, and Bayesian-blocks (Scargle et al. 2013) binnings and the results were stable among the different choices. For consistency with the other parts of the analysis, we chose the same logarithmically spaced time bins used in the time-resolved spectral analysis described in Section 3.5 below.

Consider a GRB as an impulse f(t) superimposed on a background signal b(t). Depending on the unknown shape of f(t), there will be a particular time scale δt and a particular start time t₀ maximizing the quantity

$$\begin{equation} S = \frac{\int ^{t_{0}+\delta t}_{t_{0}}f(t)dt}{\sqrt{\int ^{t_{0}+\delta t}_{t_{0}}b(t)dt}}, \end{equation} \tag{ 3 }$$

which is the significance of the signal in the Gaussian regime. The pair δt, t₀ corresponds to the highest sensitivity to the signal of this particular GRB. Our source-detection method searches for the closest pair to δt, t₀ by resizing and shifting the time bins and selecting the light curve that contains the single bin with the highest significance. Since the typical event rate inside the LLE ROI is not particularly large (∼10–20 Hz for the background), the Gaussian approximation implicit in Equation (3) is not always justified. The significance S in each bin is thus derived from the Poisson probability of obtaining the observed number of counts given the expectation from the background, by converting this probability to an equivalent sigma level for a one-sided standard normal distribution. Our algorithm starts by defining a conservative window around the trigger time, with a total duration depending on the GBM burst duration T₉₀. Then, a set of 10 bin sizes δt is defined depending on T₉₀. For each of these bin sizes, the algorithm computes 11 light curves with shifted bins, i.e., with bins centered on t₀ + (i/20) δt (i = 0...10). For each of these 10 × 11 light curves, the background function b(t) is fit to the data outside the GRB window (as described in Section 3.1) and the algorithm seeks the bin with the largest significance S inside the GRB window. This value is then corrected for the number of trials, i.e., by the number of bins N in the current light curve. If p is the probability corresponding to S, then the corrected-for-trials probability is p' = 1 − (1 − p)^N. This new probability is converted to a Gaussian-equivalent significance S' and the pre-trials significance for the detection of the GRB is defined as S_pre = max (S'), where the maximum is computed over the 110 light curves. Since the data have been rebinned multiple times, a post-trial probability is finally computed to account for these dependent trials. For this purpose, we performed 3 × 10⁶ Monte Carlo simulations of the background, running our algorithm and recording S_pre for each realization. The resulting distribution of S_pre is well described by a Lorentzian function 1 + [(x − x₀)/r_c]²]^−β, with x₀ = 1.36, r_c = 7.38, and β = 41.8 (χ² = 43.2 with 38 dof). We use this function to convert the pre-trials significance S_pre into a post-trials significance S_post.

We consider as LLE-detected the GRBs that have post-trial significances S_post > 4σ, which correspond to chance probabilities P < 3 × 10⁻⁵. We ran our algorithm on the 733 GRBs of the GBM sample (see Section 2.2) and so we expect no false-positive detections using this arbitrary threshold.

We describe the duration of a GRB detected by the LAT using the parameter T₉₀ (Kouveliotou et al. 1993). A simple measurement of T₉₀ starts with the construction of the integral distribution of the number of background-subtracted events accumulated since the trigger time. As the GRB becomes progressively fainter, the distribution flattens and eventually reaches a plateau.

The calculation of the duration of the emission consists of finding the times where the integral distribution reaches the 5% and 95% levels of its total height (called T₀₅ and T₉₅, respectively), and calculating their difference T₉₀ ≡ T₉₅ − T₀₅. Our duration estimation method is based on the above simple prescription, but is also extended to estimate the statistical uncertainty of the results and accounts for the effects of effective area variations over time (for its application to the Transient-class events).

Because of the unavoidable statistical fluctuations involved in the process of detecting incoming GRB flux, a GRB observed under identical conditions by a number of identical detectors will in general produce different detected light curves and hence different duration estimates. Our method quantifies the uncertainties on the duration estimates associated with these statistical fluctuations. In short, it accomplishes this by treating the detected light curve as the true one (i.e., that of the incoming γ-ray flux), producing a set of simulated light curves by applying Poissonian fluctuations on the detected one, estimating the durations of the simulated light curves, and calculating a single duration estimate and its uncertainty from the distribution of simulated duration estimates.

Our method starts by constructing the integral distribution of the accumulated background-subtracted events curve in small steps in time. For each step, the number of expected background events is estimated and the number of detected events is counted. At the end of each step, an algorithm checks for the presence of a plateau by searching for statistically significant increases in the average value of the points added last to the curve. If a certain number of steps does not increase the integral distribution, a plateau is reached and the construction stops. A set of simulated light curves are then produced by adding Poisson noise to the observed light curve and the corresponding integral distributions are produced. A duration estimation is made for each of the simulated light curves and the results (T₀₅, T₉₅, and T₉₀) are recorded. After the durations of all the simulated light curves have been measured, the median and a (minimum-width) 68% containment interval are calculated for each distribution and used as our measurements and ±1σ errors. In case the light curve contains multiple peaks separated by quiescent periods, the algorithm, depending on the intensity of each peak and the duration of the intermittent quiescent periods, might set the beginning of the plateau at the end of the last peak or during a quiescent period. In the latter case, some of the late emission might not be fully accounted for by the produced duration. However, the returned statistical errors would be appropriately increased in both cases, indicating the uncertainty of identifying the end of the emission.

Any changes in the off-axis angle of the GRB during an observation will change the effective area of the LAT, affecting the light curve. For example, a GRB observation that involves an ARR will in general start with a moderate to large off-axis angle that will then rapidly decrease and stay small for most of the rest of the observation. Because the effective area of this observation will be small before the ARR starts, the count rate will be artificially decreased and this would cause a bias in the measurement of T₀₅ if it were simply based on counts. To account for this effect, we weight the simulated light curves by the inverse of the exposure.

To illustrate this method, we present in Figure 3 the case of GRB 080916C and the duration measurement using the Transient-class data. These curves are used as the basis for the simulations. Figure 4 shows the distribution of T₀₅, T₉₅, and T₉₀, as measured from the simulations. These distributions are used to define the duration and associated error. In this particular case, some excess events were observed at late times (about ∼400 s), as can be seen in Figure 3. Consequently, a small fraction of the simulated light curves gave T₉₀ and T₉₅ values that were very close to ∼400 s, which caused a small increase in the duration estimates and the errors for the positive fluctuations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In some cases, a GRB observation can be interrupted before the GRB emission becomes too weak to be detectable (i.e., before reaching a plateau in the integral distribution).

Such interruptions can happen if the GRB exits the FoV of the LAT, it becomes occulted by the Earth, or the LAT enters the South Atlantic Anomaly (SAA), suspending observations. In these cases, only a lower limit on the duration can be obtained (with no errors), equal to the time interval between T₀₅ and the interruption of the observation.

We apply this method to both Transient-class data and LLE data. In the former case, we use the BKGE to estimate the background, while for the latter case we use the polynomial fit, as described in Section 3.1.1. Note however that in the calculation of the duration the exposure weighting is performed only for Transient-class data, since the effective area for the LLE class has not been characterized yet.

As a cross check, we also apply a different algorithm to the LLE data. We consider the light curve with the binning that gives the highest significance, as obtained by the algorithm explained in Section 3.3.1, and we measure T₀₅, T₉₅, and T₉₀ on the integral distribution obtained from that light curve. We verified that the numbers obtained with this simple method are always within the errors obtained with the other method. Thus, we will only provide the set of results related to the first algorithm.

We start by selecting the GBM detectors as described in Section 2 and estimate the expected backgrounds as described in Section 3.1.2. We then use the Fermi Science Tool gtbin to extract the observed spectrum (source + background) from the GBM TTE data. We obtain the response of a GBM detector in the interval to be analyzed (t₁–t₂) using the RSP2 file for the detector for the time interval. Because the RSP2 file contains several response matrices corresponding to consecutive time intervals that in general are shorter than t₁–t₂, we sum the matrices of all the sub-intervals included in t₁–t₂ using an appropriate weighting scheme. Specifically, if c_i is the counts detected in the sub-interval covered by the ith matrix and C = ∑_jc_j is the number of counts detected between t₁ and t₂, then the weight for the ith response matrix is

$$\begin{equation} w_{i}=\frac{c_{i}}{\sum _{j}c_{j}}. \end{equation} \tag{ 4 }$$

To sum the matrices, we use the tool addrmf, part of NASA HEASARC's FTOOLS.⁸¹

For the analysis of LAT observations of all GRBs detected inside the LAT FoV, we use the Transient-class events as described in Section 2. We bin the LAT data in 10 logarithmically spaced energy bins between 100 MeV and 250 GeV and use an energy-dependent ROI as described in Section 2.1.1. We derive the observed spectrum and the response matrix using the Fermi Science Tools gtbin and gtrspgen. We also use the BKGE to obtain a background spectrum containing the contributions from all the background sources, as described in Section 3.1.1.

Note that for GRBs detected by the LLE photon counting analysis outside the LAT FoV, we used only GBM data for the spectral analysis.

We load the spectra and response matrices in XSPEC v.12.7.⁸² For GBM data, we exclude from the fit all of the NaI channels between 33 keV and 36 keV (corresponding to the iodine K-edge; see Meegan et al. 2009) and ignore the channels at the extremes of the spectra (channels below 8 keV and channels 127 and 128 for NaI; channels 1, 2, 127, and 128 for BGO). We do not exclude any energy bin in the LAT spectrum, since we already selected the data before binning them. We jointly fit the GBM and LAT data with several models (described below), minimizing the negative log-likelihood. This likelihood function is derived from a joint probability distribution, obtained by modeling the spectral counts as a Poisson process and the background counts as a Gaussian process. For the latter, the Gaussian standard deviation for the ith channel is given by $\sigma _{i}= \sqrt{\vphantom{A^A}\smash{{{\sigma _{{\rm stat},i}^2+\sigma _{{\rm sys},i}^2}}}}$, where $\sigma _{{\rm stat},i}^2$ and $\sigma _{{\rm sys},i}^2$ are the statistical and the systematic variances, respectively. The maximum likelihood principle assures that the derivatives of the likelihood function with respect to the parameters are null for the best-fitting set of parameters. Exploiting this, one can treat the means of the Gaussian functions describing the background counts as nuisance parameters and remove them from the fitting procedure by expressing them as functions of the other parameters. This is a rather standard statistical procedure and leads to the formulation of a so-called profile likelihood function. PG-stat is defined as the natural logarithm of this function (see the XSPEC website⁸³ for more details). The fitting algorithms implemented in XSPEC find local minima for the statistic, but they can fail to converge to the global minimum. This is a known issue with gradient-descent algorithms (Arnaud et al. 2011). To mitigate this problem, we perform multiple fits (from 10 to 40) for each model, each time starting from a different set of values for the parameters, and we keep as the putative best fit the set giving the lowest overall value for the statistic. If the fitting algorithm finds an even better minimum for the statistic while computing error contours for this set of parameters, we adopt that as the new putative best fit and restart the error computation, iterating the procedure until no new minimum is found.

Traditionally, GRB spectra have been described using the phenomenological "Band function" (Band et al. 1993) or a model consisting of a PL with an exponential cutoff (also called a "Comptonized model"). Another common choice is the smoothly broken power law (SBPL; Ryde 1999). Recently, the logarithmic parabola has been shown to be a good description the spectra of some GRBs, especially in time resolved analyses (Massaro et al. 2010; Massaro & Grindlay 2011). We call these four spectral models main components. One of the first results by Fermi was the need for multi-component spectral models for some GRBs, showing a high-energy excess over the main component that has been modeled with an additional PL (Ackermann et al. 2010b; Abdo et al. 2009a). In one case, Fermi observed a high-energy cutoff that required the addition of an exponential cutoff to the PL component in the spectral model (Ackermann et al. 2011), for a total of three components (Band, power law, and exponential cutoff). In the following, we will call the PL and the exponential cutoff functions additional components, to emphasize the fact that we add them to the main components when needed. Some authors have claimed the presence of a thermal component, modeled by a blackbody emission spectrum (see e.g., Guiriec et al. 2011; Zhang et al. 2011, and references therein). However, a careful time-resolved analysis is needed in order to investigate and characterize such a component, which is beyond the scope of the present analysis. Thus, we did not include a blackbody component in our spectral fits. Hereafter, N(E) is the differential photon flux (in units of cm⁻² s⁻¹ keV⁻¹) expected from a model at a given energy E (in keV) and k is a normalization constant whose units depend on the model. We have four main model components:

• 1.  
  Comptonized model (a PL with an exponential cutoff):

  $$\begin{equation} N(E)\equiv kE^{-\alpha }e^{-\frac{E}{E_{0}}}, \end{equation} \tag{ 5 }$$

  where α is the photon index and E₀ is the cutoff energy.
• 2.  
  Logarithmic parabola, defined following Equation (9) in Massaro et al. (2010):

  $$\begin{equation} N(E) \equiv \frac{S_{p}}{E^{2}}10^{-b \log {E/E_{p}}^{2}}, \end{equation} \tag{ 6 }$$

  where S_p is the height of the spectral energy distribution at the peak frequency, E_p is the peak energy, and b represents the curvature of the spectrum.
• 3.  
  Band model (Band et al. 1993): two PLs joined by an exponential cutoff:

  $$\begin{eqnarray} B(E) &=& N(E) \nonumber\\ &\equiv & k \left\lbrace \begin{array}{@{}ll@{}}E^{\alpha }e^{-E/E_{0}} &{when } E < (\alpha -\beta)E_{0}\\ \left[(\alpha -\beta)\,E_{0}\right]^{\alpha -\beta }&\\ \quad \times E^{\beta }e^{-(\alpha -\beta)} &{\rm when } E > (\alpha -\beta)E_{0}\\ \end{array}. \right.\nonumber\\ \end{eqnarray} \tag{ 7 }$$

  Note that this is the representation that uses the e-folding energy E₀ (keV) instead of the peak energy E_p, where E_p = (2 + α)E₀. α and β are the (asymptotic) photon index at low energy and the photon index at high energy, respectively.
• 4.  
  Smoothly broken power law (Ryde 1999): two PLs joined by a hyperbolic tangent function with an adjustable transition length:

  $$\begin{equation} N(E)\equiv k \left(\frac{E}{E_{{\rm piv}}}\right)^{\frac{\alpha +\beta }{2}} \left[ \frac{\cosh {\left(\frac{\log {(E/E_{0})}}{\delta }\right)}}{\cosh {\left(\frac{\log {(E_{{\rm piv}}/E_{0})}}{\delta }\right)}} \right]^{\frac{\alpha +\beta }{2}\delta \log _{e}{(10)}}, \end{equation} \tag{ 8 }$$

  where E_piv is a fixed pivot energy, α and β are, respectively, the photon index of the low-energy and the high-energy PLs, E₀ is the e-folding energy, and δ is the energy range over which the spectrum changes from one PL to the other.

Here are the definitions of our additional model components:

• 1.  
  Power law: N(E) ≡ kE^−α, where α is the photon index.
• 2.  
  Exponential cutoff: $e^{-\frac{E}{E_{0}}}$.

Because of the variety of spectral models, we have considered a number of functions composed of one main component and one or more additional components: Band, Band + power law, Band + power law with exponential cutoff ($\equiv B(E)+kE^{-\alpha }e^{-E/E_{0}}$), Band with exponential cutoff ($\equiv B(E)e^{-E/E_{0}}$), Comptonized, Comptonized + power law, Comptonized + power law with exponential cutoff, logarithmic parabola, SBPL, and SBPL + power law.

To take into account the relative unknown uncertainties in the inter-calibration between the different detectors, for bright bursts we also apply an effective area correction (Bissaldi et al. 2011): we scale the model under examination by a multiplicative constant, with the constant being fixed to 1 for the LAT (taken as reference detector), but free to assume different values for all the other detectors. For GRBs for which we do not use LAT data, we choose one of the NaI detectors as the reference. While for bright bursts adding such a correction changes the best-fit parameters and the value of the statistic, for the other bursts it is essentially inconsequential, since in the latter cases statistical errors dominate over the inter-calibration uncertainties. For such spectra, the multiplicative factors are unconstrained during the fit and therefore we removed them. After the best fit is found, we fix all the factors to their best-fit values and we proceed with the error computation. The correction factors typically have values between 0.95 and 1.05 for the NaI detectors and between 0.75 and 1.25 for the BGO detectors.

The main focus of the spectral analysis performed here is to characterize the GRB spectrum, which requires selecting the most appropriate spectral model. We define the best model for a given GRB as the simplest one giving a "good" value for the test statistic (PG-stat, S in the following) and no evident structures in the residuals. Since S is based on a Poisson likelihood, we do not have a simple goodness-of-fit test comparable with the χ² test when minimizing the χ² statistic. The actual expected value S* for the statistic S is a function of the number of counts N in the spectrum and the background model and its uncertainties and can be estimated using Monte Carlo simulations. We assume a model $m_{0}(\boldsymbol {p})$ (for example, the Band model) with the best-fitting set of parameters $\boldsymbol {p}_{0}$ as the null hypothesis and we generate 1 million realizations of $m_{0}(\boldsymbol {p}_{0})$ and the corresponding background spectrum using the fakeit command of XSPEC. Each realization $r^{i}_{p_{0}}$ is obtained by adding Poisson noise to the count spectrum obtained by summing the observed background spectra and $m_{0}(\boldsymbol {p}_{0})$. Correspondingly, each realization of the background spectrum is obtained by adding Gaussian noise to the observed background spectrum, using a total variance composed of the statistical and the systematic variance of the observed background. Then, we fit $m_{0}(\boldsymbol {p})$ to $r^{i}_{p_{0}}$ and we record the value for the statistic S_i. In Figure 5, we show an example of a distribution for S obtained using the Band model and a χ² distribution for the same number of dof as reference. Note that depending on the case, the two distributions can be very different. We can now use the distribution for S resulting from the simulations to compute the probability of obtaining the observed value for S under the null hypothesis m₀. This approach requires a large number of simulations, so we applied it just for the subsample of GRBs for which we claim the detection of an extra component (see below and Section 4.4.1).

Zoom In Zoom Out Reset image size

In order to compare different models, we considered them in pairs. Let us consider the model m₀ with n_{0, dof} and m₁. If S₀ < S₁ and n_{0, dof} ⩽ n_{1, dof}, then m₀ describes the data better using fewer or the same number of parameters and we consider it a better fit following the definition given at the beginning of this section. If S₀ ≃ S₁ and n_{0, dof} = n_{1, dof}, the two models are equivalent and we should report the results for both models. Anyway, this never happened in our analysis. On the other hand, if S₀ > S₁, then m₁ better fits the data and we have to decide if the improvement is significant enough to justify the added complexity. In the literature, there are different ways to quantify this improvement, sometimes incorrectly (see, for example, the discussion in Protassov et al. 2002). One of the standard methods is the likelihood-ratio test, which uses as a test statistic the difference in S (ΔS) between the two models. In the case of nested models m₀ and m₁, Wilks' theorem (Wilks 1938) assures under certain hypotheses that the quantity ΔS asymptotically follows a χ² distribution with n = n_{0, dof} − n_{1, dof} dof. Unfortunately, in all the cases of interest here, the theorem's hypotheses are not satisfied and the reference distribution for ΔS is not known. In general, one should perform dedicated Monte Carlo simulations to obtain the reference distribution. Performing such simulations for each pair of models is not practical. Thus, we select three cases of interest (i.e., Band versus Band + power law, Band versus Band with exponential cutoff, and Band + power law versus Band + power law with exponential cutoff) and we perform several million simulations to evaluate the reference distributions. We use the same procedure as above, using the simplest model as the null hypothesis, but we fit both m₀ and m₁ to each simulated dataset, recording ΔS. At the end of the simulation, the distribution for ΔS is used to compute the probability P of obtaining a ΔS greater than the observed value, which corresponds to the complement of the cumulative distribution function. In Figure 6, we plot this function for the three cases. We fix an arbitrary threshold at P_th(> ΔS) = 1 × 10⁻⁵, where the statistical error on the simulated distribution, visible toward the tail, is still low. P_th corresponds to a significance level of ∼4.2σ and defines a threshold for ΔS above which we claim a significant detection of an extra component. Specifically, P_th corresponds to ΔS = 25 for Band versus Band + power law, ΔS = 28 for Band versus Band with exponential cutoff, and ΔS = 20 for Band + power law versus Band + power law with exponential cutoff.

Zoom In Zoom Out Reset image size

Figure 11 shows the flux and fluence measured by the LAT in the "GBM" (top two panels) and "LAT" (bottom two panels) time intervals as a function of the durations of these time intervals (i.e., GBM and LAT T₉₀, respectively). The fluxes and fluences presented in these figures are for the 100 MeV–10 GeV energy range. As can be interestingly seen in the bottom right panel of Figure 11, within the first 3 years of operations the LAT detected four very high fluence bursts (GRBs 080916C, 090510, 090902B, and 090926A) that are outliers with respect to the main distribution of the LAT-detected GRBs. We will revisit these hyper-fluent bursts in Section 5.2, where we discuss the energetics of Fermi-LAT detected GRBs.

Zoom In Zoom Out Reset image size

We evaluate localizations from the LAT for all GRBs detected by the maximum likelihood analysis by searching for the maximum of the TS map according to the procedure described in Section 3.2.2. We present our results in Table 5, in which we report the position of the maximum of the TS map (i.e., the LAT localization) along with its 68%, 90%, and 95% statistical errors.

We report the energies and arrival times of a set of interesting high-energy photons that, according to our likelihood analysis (as described in Section 3.2.3), have a high probability (P > 0.9) of being associated with the GRBs. Specifically, we give information for the following events:

• 1.  
  The highest energy Transient-class LAT γ-ray in the "GBM" time window (Table 6);
• 2.  
  The highest energy Transient-class LAT γ-ray in the interval starting from GBM T₉₅ and ending at the end of the "EXT" window (i.e., from the end of the measured duration in the GBM data up to the end of the LAT measured duration; Table 7);
• 3.  
  The highest energy Transient-class LAT γ-ray detected in the time-resolved likelihood analysis (Table 8).

The results are shown in Tables 6–8. These results show that the detection of high-energy events with GRB point source probabilities P > 0.9 is not strongly correlated with features in the GBM light curve. In a few cases, such as GRB 090510, such events are coincident with bright pulses in the GBM light curve, but more often the most energetic event is detected after the intense low-energy emission, as with the 33.39 GeV event detected at T₀+81.75 s from GRB 090902B, which is the highest energy ever observed from a burst. GRB 100728A is particularly interesting since a 13.54 GeV event was detected ∼90 minutes after the trigger time. This is the only case in which we observe such a late event and it can potentially confirm that high-energy γ-rays can arise very late in time, as observed from GRB 940217 by EGRET (Hurley et al. 1994). On the other hand, GRB 100728A is not significantly detected at the time that the highest energy event is observed (similar to GRBs 090217 and 100116A, reported in Table 8), thus the probability P = 0.987 that the 13.54 GeV event is associated with the burst must be taken with caution. Considering the trials factors, this probability would be further reduced, weakening the case for high-energy emission that persists for hours. A detailed analysis of the probability corrected by the trials factors would be non-trivial as the background strongly varies as a function of the location in orbit and is beyond the scope of this paper. Figure 12 shows the energies and arrival times for the highest energy γ-rays associated with LAT GRBs. The estimated errors are computed from the energy dispersion in the IRFs and are of the order of 10% for energies >1 GeV. When possible, we also indicate the source-frame energy.

Zoom In Zoom Out Reset image size

To study the temporal decay of the extended emission detected by the LAT, we utilized the time-resolved analysis described in Section 3.5. We first visualized any detected extended emission using flux light curves (shown in Appendix B) and then calculated the peak-flux value F_p and the time of the peak flux t_p, quantities shown in the two top panels of Figure 14. In the time-resolved analysis, we adaptively changed the size of the time bin width in order to significantly detect the source, so F_p corresponds to the average flux in the time bin of the maximum. As a result, it is more precise (i.e., with a smaller uncertainty) for bright GRBs than for faint GRBs.

The four most luminous bursts detected by the LAT have some of the highest peak fluxes in the ensemble, all exceeding 10⁻³ cm⁻² s⁻¹. Among the rest of the bursts, GRBs 081024B and 110721A also have notably high peak fluxes. GRB 100728A was at the edge of the FoV at the time of the GBM trigger and was detected only at later times. It has by far the lowest peak flux of all GRBs, at least an order of magnitude lower than the rest of the population; however, its value is possibly affected by large systematic errors.

We also applied the methods described in Section 3.5 to the subsample of GRBs with detected extended emission. We detected temporal breaks in the decay of the extended emission of three bright GRBs: 090510, 090902B, and 090926A. In the top panel of Figure 13, we show their luminosities as a function of rest-frame time, as well as the best-fitting broken PL models. The later points in the light curves are very important to constrain the break, but they also would be the most affected by any unaccounted-for systematic uncertainties arising, for example, from the background estimation or the exposure calculation. In the bottom panel of Figure 13, we again report the luminosity as a function of the rest-frame time, but for all the GRBs in the subsample. In Table 9, we report the results of this analysis. For the three GRBs with temporal breaks, we report the decay index starting from the peak flux and before the break α₁, the decay index after the break α₂, and the break time t_b. For all other GRBs, we report the decay index for the whole extended emission starting from the peak flux and the decay index for the light curve starting from the end of the low-energy (GBM) emission.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Referring to Table 9, we also define the "late-time decay index" α_L, which corresponds to the decay index measured after the GBM T₉₅ (α_L = α) for all GRBs except the three for which we detect temporal breaks, for which it corresponds to the decay index after the break (α_L = α₂). In the third panel of Figure 14, we report α_L for all of the GRBs of the subsample.