Searches for Neutrinos from Gamma-Ray Bursts Using the IceCube Neutrino Observatory

An unbinned likelihood ratio method (Achterberg et al. 2006; Braun et al. 2008; Aartsen et al. 2017a) is combined with frequentist statistics to assign a probability that a subset of neutrino candidate events is consistent with background. For a sample of N candidate neutrino events with characteristics x_i , the likelihood can be written as

$$\begin{eqnarray}&&{ \mathcal L }({n}_{s}| {n}_{b},{x}_{i})={P}_{N}\prod _{i=1}^{N}\left[{p}_{s}S({x}_{i})+{p}_{b}B({x}_{i})\right],\end{eqnarray} \tag{ 2 }$$

where p_s = n_s /(n_s + n_b ), p_b = n_b /(n_s + n_b ), and P_N is the Poisson probability to observe N events, assuming n_s signal events and n_b expected background events,

$$\begin{eqnarray}&&{P}_{N}=\displaystyle \frac{{({n}_{s}+{n}_{b})}^{N}}{N!}{e}^{-({n}_{s}+{n}_{b})}.\end{eqnarray} \tag{ 3 }$$

S and B denote the probability distribution functions (PDFs) describing the spatial and energy distribution of signal events and background events, respectively. For the signal energy PDF, an E^−γ power-law spectrum is assumed. The use of an energy spectrum is to improve the sensitivity to identify GRB neutrinos from the background, without assuming a specific model of GRB neutrino emission. In some of the analyses, the energy spectrum of the search is allowed to float for each GRB; in other analyses it is fixed to a generic E⁻² spectrum based on first-order Fermi acceleration. The signal space PDF, shown in Equation (4), uses a 2D Gaussian to test the compatibility of the neutrino candidate's reconstructed position, x _ν , with the source position, x _GRB,

$$\begin{eqnarray}&&{S}_{\mathrm{sp}}({{\boldsymbol{x}}}_{\nu },\sigma | {{\boldsymbol{x}}}_{\mathrm{GRB}})=\displaystyle \frac{1}{2\pi {\sigma }^{2}}\exp \left(-\displaystyle \frac{| {{\boldsymbol{x}}}_{\nu }-{{\boldsymbol{x}}}_{\mathrm{GRB}}{| }^{2}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 4 }$$

where σ² is the quadratic sum of the uncertainty on the reconstructed neutrino and GRB direction. Since the contribution of signal events to the total data set is expected to be extremely small, the data can be used to construct the background PDFs. In particular, the background energy PDF is assumed to follow the energy distribution of data. Due to the detector geometry, the approximation of azimuthal symmetry can be used to describe the background space PDF solely as a function of zenith. Any neutrino emission that occurs within a given analyzed TW is fitted as constant emission during the TW. The specific TWs used in each analysis are described in Section 5.

The likelihood Equation (2) is evaluated for different values of n_s using the PDFs described above, with ${\hat{n}}_{s}$ designating the value of n_s that maximizes the likelihood. A likelihood ratio with respect to the null hypothesis L(n_s = 0) is then used to obtain the following test statistic (TS):

$$\begin{eqnarray}&&\mathrm{TS}=2\cdot \mathrm{ln}\left[\displaystyle \frac{L({\hat{n}}_{s})}{L({n}_{s}=0)}\right]=-2{\hat{n}}_{s}+2\sum _{i=1}^{N}\mathrm{ln}\left[\displaystyle \frac{{\hat{n}}_{s}S({x}_{i})}{{n}_{b}B({x}_{i})}+1\right].\end{eqnarray} \tag{ 5 }$$

High-energy events in temporal and close spatial coincidence with a GRB will result in a large TS value. The p-value of an observation is determined by comparing the observed test statistic to a test statistic distribution created from scrambled data (see Section 4.4).

Two of the four analyses use localizations provided by Fermi-GBM for some GRBs and therefore have an additional step to determine the TS. Starting in early 2018, the Fermi-GBM collaboration began releasing a HEALPix skymap (Górski et al. 2005) with the localization probability as a function of sky position for each GRB localized by GBM (Goldstein et al. 2020). Maps prior to 2018 were processed in a similar way using the GBM Data Tools, but the metadata in the files have not been fully qualified and the files have not yet been uploaded to the final HEASARC archive, therefore we use the preliminary files from Goldstein & Wood (2022). These maps contain a per-pixel probability that a given GRB originates from that direction. An all-sky scan of neutrino data is performed, in which TS_orignial (equal to TS from Equation (5)) is calculated at every pixel of the skymap and then penalized by the probability, P_GBM, of that pixel:

$$\begin{eqnarray}&&{\mathrm{TS}}_{\mathrm{final}}={\mathrm{TS}}_{\mathrm{original}}+2\times \left[\mathrm{ln}({P}_{\mathrm{GBM}})\ -\mathrm{ln}({P}_{\mathrm{GBM},\max })\right],\end{eqnarray} \tag{ 6 }$$

where ${P}_{\mathrm{GBM},\max }$ is the maximum probability on the entire skymap. The position of the GRB on the sky and the number of signal events are thus both fitted to find the combination which maximizes TS_final.

Equation (2), which describes the likelihood for a single well-localized GRB, can be easily modified to describe the likelihood of N_GRB well-localized GRBs. When considering multiple sources, n_s corresponds to the total number of signal events summed over all GRBs. Each GRB is assumed to have the same time-integrated neutrino flux at Earth. The expected number of neutrinos from each GRB will therefore be proportional to the effective area at the decl. of the burst, ${{ \mathcal A }}_{\mathrm{eff}}(\delta )$, which is calculated assuming an E⁻² neutrino energy spectrum. The signal PDF is then replaced by the sum of the N_GRB signal PDFs, weighted by their relative contribution to the number of signal events:

$$\begin{eqnarray}&&S=\displaystyle \frac{{\sum }_{j=1}^{{N}_{\mathrm{GRB}}}{{ \mathcal A }}_{\mathrm{eff}}({\delta }_{j})\cdot {S}_{j}}{{\sum }_{j=1}^{{N}_{\mathrm{GRB}}}{{ \mathcal A }}_{\mathrm{eff}}({\delta }_{j})},\end{eqnarray} \tag{ 7 }$$

where S_j and δ_j are the signal PDF and decl. of the jth burst, respectively. This method, known as a stacked likelihood analysis, provides a single measure of the total neutrino emission from a set of GRBs.

When not performing a stacked search, analyzing a selection of N GRBs will result in a p-value for each individual burst. A trial-correction method is thus needed to determine if one or more of the obtained p-values provides a statistically significant result. Arranging the p-values from smallest to largest, their values are denoted as p₁, p₂,...,p_N . Under the null hypothesis of no neutrino emission, the correlations of neutrino events with GRBs will only occur randomly, and these N p-values are expected to follow a uniform distribution between 0 and 1. The probability that k or more p-values are smaller than or equal to p_k is thus given by the binomial probability:

$$\begin{eqnarray}&&P{(k)\equiv P(n\geqslant k| N,{p}_{k})=\sum _{m=k}^{N}\displaystyle \frac{N!}{(N-m)!m!}{p}_{k}^{m}(1-{p}_{k})}^{N-m}.\end{eqnarray} \tag{ 8 }$$

Evaluating Equation (8) over all potential values of k, the smallest P(k) is selected. An empirical trial-correction factor is then applied to account for the fact that N potential p-values were scanned. This trial correction is found by determining the fraction of background-only realizations that produce a more significant result. A visualization of this trial-correction procedure is shown in Figure 1. This procedure also ensures that the rare overlap of a single neutrino contributing to two analyzed GRBs (thus creating a correlation between two p-values) is also accounted for in the final post-trial p-value. On the other hand, given that the successive probabilities P(k) are strongly correlated with one another, the overall trial-correction factor is modest in the end.

Zoom In Zoom Out Reset image size

To test the neutrino flux to which the analyses are sensitive, simulated muon-neutrino and muon-anti-neutrino events based on full detector Monte Carlo (MC) can be injected into the data sample. The energy spectrum of the MC-generated neutrinos can be fixed to E⁻² or other spectra as desired. To describe the background, we use data with event times that are randomly selected from a uniform distribution over the livetime of the data set while accounting for the detector downtime. Keeping the event coordinates fixed in the detector frame, scrambling background event times correspondingly randomizes the R.A. Such a pseudo-data set is called scrambled data. Signal can then be "injected" by adding the signal-like events to it. In each realization of scrambled data plus signal, the number of injected MC events is drawn from a Poisson distribution with a fixed mean μ. By varying the value of μ, the threshold can be determined at which 90% of all realizations result in a TS value that is larger than the median of the background TS distribution of the analysis. We define this threshold, and the neutrino flux to which it corresponds, as the sensitivity of the analysis. Once the analysis has been performed, upper limits can be calculated in a similar way, by comparing to the observed TS in data rather than the median TS of the background-only realizations.

This analysis searches for neutrinos in a range of TWs and uses the largest GRB catalog of the four analyses. All GRBs with known duration were included, given they were observed during the detector livetime. 163 GRBs do not have a reported T₁₀₀ in GRBweb at the time of writing, and so these were excluded from the selection. The final number of GRBs in this analysis is 2091.

The p-values are calculated for these 2091 GRBs in 10 pre-determined TWs and with a fixed E⁻² neutrino energy spectrum in the signal hypothesis. The first nine TWs range from 10 s to 2 days, centered on the T₁₀₀ of the GRB, and the final TW is asymmetric with a 1 day precursor and 14 day afterglow TW (see Figure 2(a)). The shortest predefined TW that fully envelopes the T₁₀₀ interval will be used to describe the prompt emission. The seven shortest TWs range up to 10³ s, which includes most GRB prompt phases. The 10 s TW is sufficiently small to study the short GRBs, because the neutrino background is effectively zero at this timescale. For the longer TWs the background becomes non-negligible, which is what motivates the decision to search up to around two weeks but not longer. For a well-localized GRB (σ ≲ 1°), the 15 day TW has an expected background on the order of one neutrino candidate event. For poorly localized GRBs this longer window is less sensitive, which is why the emphasis was placed on the afterglow region where higher-energy neutrinos are predicted (Asano & Murase 2015). The TS in each TW is calculated using Equation (6). From those 10 p-values, the most significant one is selected to represent the GRB.

Zoom In Zoom Out Reset image size

Since each GRB is studied in 10 TWs, a correction is required to compensate for the look-elsewhere effect. Because the TWs are correlated, an effective trial correction is used. For every given GRB, the smallest p-value of the 10 TWs is selected. Background scrambles are then used to determine the probability of obtaining a value that is smaller than or equal to this result. Next, the p-values that have been corrected for searching 10 TWs are evaluated in a binomial test (Section 4.3) to search for evidence of a sub-set of GRBs with significant neutrino emission.

Four binomial tests are evaluated, where the total GRB sample has been divided into four sub-populations by duration and hemisphere. GRBs are separated into short and long classes according to whether the measured T₉₀ is less than or greater than 2 s, respectively. This is intended to account for the different progenitor classes of merger events (Abbott et al. 2017) and core-collapse supernovae (Hjorth & Bloom 2012), albeit in a simplistic way since the two classes are known to overlap. This overlap explains why GRB 170817A falls into the long GRB class with a T₉₀ of 2.048 s despite being a known short GRB from a binary neutron star merger. In this case, its misclassification does not notably affect the results as GRB 170817A would not have appeared in the top five short GRBs from the Southern Hemisphere (listed in Appendix A.2) even if it were included in that class. Future studies will explore grouping methods which account for the fact that GRB 170817A is a short GRB. The split by hemisphere allows for the increased sensitivity of the analysis to GRBs in the northern sky.

The final p-value for each binomial test is determined by comparing the result found for data with the results found for scrambled data sets. The final p-values are summarized in Table 2 and are consistent with background. The most significant GRB (pre-trial) from each sub-population is listed in Table 3. All GRBs with a p-value less than 1% are listed in Appendix A.2.

In the "Precursor/Afterglow" analysis, the TW for the signal region is treated as a parameter which can be estimated from the data. The TW can be fitted from a minimum duration of 0.5 s up to a maximum of two weeks, long enough to cover a wide range of possible neutrino emission timescales. Only GRBs with a positional angular uncertainty of less than 02 are considered in this analysis (thus they can be approximated as point sources), which keeps the background low within the full TW range. GRBs within the first 14 days and the last 14 days in the GFU data are further removed for this analysis. This results in a selection of 733 GRBs. This includes 53 GRBs that did not have a reported T₁₀₀ in GRBweb and were not included in the Extended TW analysis. Two separate searches are performed for every GRB, denoted as the "precursor" search and "prompt+afterglow" search, each with fit parameters corresponding to the number of signal events (n_s ), the spectral index (γ), and the width of the emission TW (T_w ). The only difference between the two searches is how the TW T_w is defined with respect to the start of the prompt phase T₀ (see Figure 2(b)):

• 1.  
  For precursor searches, the TW extends backwards from T₀ to the time (T₀ − T_w ) before the start of the prompt phase.
• 2.  
  For prompt+afterglow searches, the TW extends forwards from T₀ to the time (T₀ + T_w ) after the start of the prompt phase.

For every GRB, the analysis returns the best-fit parameters for the respective search and the p-value. This results in two lists of 733 p-values, one list for each search. A binomial test (Section 4.3) is performed on each list to search for a subset of GRBs with significant neutrino emission. The final p-value after each binomial test is determined by comparing the result found for data with the results found for scrambled data sets. Figure 3 demonstrates different numbers of individual GRBs with 2σ, 3σ, and 4σ p-values that can typically result in a final post-trial p-value of 3σ significance.

Zoom In Zoom Out Reset image size

The final post-trial p-values of the binomial tests are 0.495 for the precursor search and 0.486 for the prompt+afterglow search. The results are comparable with the median background expectation (i.e., p-value of 0.5). The results for the top 20 GRBs in each search are summarized in Appendix B.

The two analyses above searched for neutrino correlations with GRBs in an agnostic manner. It is possible to use additional gamma-ray observations to perform more sensitive dedicated analyses. A previous study by Coppin et al. (2020) analyzed the light curves of all Fermi-GBM bursts up to 2020, leading to a sample of 217 GRBs that exhibit signs of gamma-ray precursors. Of those 217 bursts, there are 133 GRBs that overlap with the IceCube data set examined here. A dedicated search is performed to investigate potential neutrino production coincident with the emission phase of the observed gamma-ray precursors using these 133 GRBs. The TWs of the analysis are set to those of the identified gamma-ray precursors, extended by 2 s on either side to obtain a restrictive yet conservative range. Summed over all 133 GRBs, a total TW of 3.3 × 10³ s is examined. The signal hypothesis in the likelihood uses a fixed E⁻² neutrino energy spectrum. Out of 133 bursts, 100 GRBs were localized solely by Fermi-GBM. For those 100 bursts, a TS is defined following Equation (6). Otherwise, the TS is based on Equation (5).

Performing the analysis on all GRBs results in 133 individual p-values. To trial-correct the result, a procedure similar to that of the binomial correction is applied. However, instead of considering the kth most significant p-value, the product of the k most significant p-values is considered. Comparing this result to the distribution of the same statistic (product of the k-most significant p-values) found in scrambled data sets yields the final p-value for the analysis.

As this analysis only targets GRBs with observed gamma-ray precursors, it has the smallest background of all four analyses. The neutrino flux to which the analysis is sensitive (Section 4.4) corresponds to only 1.7 neutrino events on average. Considering events within a relative combined neutrino and GRB angular uncertainty of 5σ, no neutrino events were observed in temporal coincidence with any of the 133 precursors. A final p-value of 1 was thus obtained.

One of the results by Coppin et al. (2020) was that almost all (>95%) gamma-ray precursors were found to occur within 250 s of the prompt emission. A fourth analysis was therefore performed to enable a larger systematic search for precursor neutrinos, not limited to Fermi-GBM bursts for which a gamma-ray precursor could be resolved. For all well-localized ⁶⁷ GRBs within the IceCube GFU data period, corresponding to 872 bursts, a generic TW of 250 s prior to the T₀ time was searched for excess neutrinos. All 33 well-localized GRBs from the "GBM Precursor" analysis are included in this search. For the set of 872 GRBs, a single stacked test statistic is constructed according to Equation (5) and assuming an E⁻² signal spectrum. Due to the stacking procedure, this analysis automatically leads to a single p-value, not requiring a trial-correction procedure.

When the stacked precursor analysis is applied to these 872 GRBs, no excess of neutrino events is found (n_s = 0). Overall, five low-energy events within the given TW are in loose spatial coincidence with the examined GRBs, listed in Table 4. Given the length of the TW and the number of bursts, these coincident events are fully consistent with the background expectation. Similar to the GBM precursor analysis, a final p-value of 1 is thus obtained.

Table 4. Properties of the Neutrinos Arriving in Coincidence with GRBs in the Precursor Stacking Search

GRB Name	Angular Separation	Loc. Uncertainty	Energy Proxy	Time Delay (s)
GRB 130131B	103	26	676 GeV	54.0 s
GRB 141220A	20	22	47 GeV	247.3 s
GRB 160314B	58	12	1023 GeV	158.4 s
GRB 160705B	52	15	794 GeV	91.5 s
GRB 160912A	61	23	525 GeV	106.4 s

Note. For each neutrino, the angular separation from the GRB, angular uncertainty of the neutrino direction, a proxy for the event energy, and the arrival time before the GRB is shown.

Download table as:  ASCIITypeset image

The stacked TS presented in Section 4.1 is used to place an upper limit on contributions of GRBs to the quasi-diffuse flux measured by IceCube. Flux is injected using an E^−2.28 (Stettner 2019) spectrum until 90% of trials yield a stacked TS above the unblinded value. This injected flux is converted to a diffuse flux using Equation (9). This procedure is repeated for all 10 TWs and the prompt. For the prompt, the shortest TW that includes the entire reported T₁₀₀ is used (see Figure 2).

These 90% confidence level limits are set for each TW for various subsets of GRBs. The four sub-populations analyzed with a binomial test are presented in Figure 4. Limits are also placed on all GRBs observed in the northern and southern sky, as well as all short and long GRBs (Figure 5). In each plot, the stacked limit from all 2091 GRBs is shown for reference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Previous IceCube studies (Aartsen et al. 2017a) constrained the prompt contribution of GRBs observable by current gamma-ray satellites to ∼1% of the diffuse flux observed by IceCube. The prompt limit presented here applies to all GRBs in the universe and corresponds to ≲1%. The limits are similar despite analyzing nearly twice as many GRBs in this analysis. The difference is the inclusion of _z term to account for GRBs that are too far away to be observable with current gamma-ray satellites. Given its fluence trigger threshold of ∼10⁻⁸ erg cm⁻² s⁻¹, the Swift-BAT detector can be assumed to view all canonical GRBs with an isotropic equivalent luminosity L_iso ≥ 10⁵⁰ erg cm⁻² up to a redshift of z ∼ 1.3. Assuming that all GRBs have identical properties in terms of neutrino emission and that they follow the redshift evolution described by Lien et al. (2014), the contribution from GRBs outside the observable redshift threshold can be calculated using the procedure outlined in Ahlers & Halzen (2014). This results in _z value of ∼2.1. The inclusion of this _z term leads to a similarly constraining prompt limit compared to previous studies.

These results were used to calculate limits for the total contribution of all long GRBs to the observed diffuse neutrino flux (see Figure 6). Based on the study by Lien et al. (2014) using Swift observations of GRBs to extrapolate to the whole universe, a mean rate of 4571 long GRBs per year is estimated to be beamed at Earth. We use the software package FIRESONG (Tung et al. 2021) to simulate neutrino emission from this cosmic GRB population, with different emission periods. Emission windows ranging from 100 s to 14 days were considered in the simulations. The redshift distribution of the cosmic GRB population from the study by Lien et al. (2014) was used to simulate the GRB rate density in FIRESONG. We assume the case where no luminosity evolution is required and for simplicity assume that the GRBs are standard candle neutrino sources.

Zoom In Zoom Out Reset image size

The "Precursor/Afterglow" analysis only considers 733 GRBs across 7.16 yr, out of which 556 GRBs are long GRBs which were observed by Swift. Hence we use the GRB detection efficiency as a function of redshift from the study by Lien et al. (2014) as well as the sky coverage and the survey time of Swift/BAT to down-select from the 4571 sources created by the FIRESONG simulations to a sample of 556 GRBs to recreate the observation biases of the Swift sample. These 556 flux values down-selected from the simulation are used to inject signal with an E^−2.28 (Stettner 2019) spectrum into a simulated neutrino data set. To repeat the original analysis under the same conditions, an additional 187 sources are added with no signal (i.e., represent background only), bringing the sample to 733 sources again. The likelihood analysis and binomial test are then performed as before on this simulated sample. This sequence of steps for every emission window considered is repeated using different simulated neutrino data sets to produce 1000 trials. When GRBs are assumed to produce the entire diffuse flux (Stettner 2019), 100% of these injected trials result in a binomial test result more significant than the unblinded result. The fraction of the diffuse flux is then reduced to identify the flux where 90% of trials produce a binomial test result more significant than the unblinded result. This is summarized in Figure 6, where the points show this upper limit (90% confidence level) for a range of neutrino emission time windows. The respective fraction for the different emission durations considered represents the total allowed contribution of all long GRBs to the observed astrophysical diffuse neutrino flux.

Precursor neutrinos have been predicted in models in which a jet initially has to burrow its way through remnant layers of the progenitor star. A prediction for the diffuse precursor neutrino flux from such sources was made by Razzaque et al. (2003) and is shown in Figure 7. The red and green lines correspond to progenitors that have a remnant outer hydrogen (H) or helium (He) shell, respectively. This analysis is able to exclude the H-shell model by a factor 10, but cannot constrain the He-shell model. To be consistent with the model prediction, the H-shell upper limit shown in Figure 7 assumes that the diffuse GRB neutrino flux results from 10³ GRBs per year (Razzaque et al. 2003). This is in contrast to the model-independent upper limits in Figure 7, which rely on the conversion outlined in Equation (9). This latter approach incorporates updated information about the redshift distribution of GRBs that was unavailable when the model was released. As a result, these generic limits are slightly more conservative.

Zoom In Zoom Out Reset image size

The GRB name is based on GCN notices, with a preference for the name provided by Fermi-GBM. If the GRB is not observed by Fermi-GBM, then the name is based on the date of observation with the format YYMMDD. The final letter indicates the order of detection if more than one GRB is observed in a single day. The R.A. (α) and decl. (δ) in J2000.0 equatorial coordinates, as well as the localization uncertainty (σ) are all listed in degrees. T₀ indicates the MJD trigger time of the GRB. The T₁₀₀ is provided in seconds.

Table 9. Most Significant Short GRBs from the Southern Hemisphere

GRB Information							Fit Results
GRB Name	α	δ	σ	T0	z	T100	TW	TS	${\hat{n}}_{s}$	p-value
GRB 140511095	329.76	−30.06	8.50e+00	56788.095	⋯	1.41	2 days	20.27	2.93	9.20e-03
GRB 130919173	297.35	−11.73	5.23e+00	56554.173	⋯	0.96	5000 s	16.20	1.00	1.75e-02
GRB 160411A	349.36	−40.24	2.85e-04	57489.062	⋯	1.26	2 days	8.08	0.99	2.36e-02
GRB 120212353	303.40	−48.10	1.05e+01	55969.353	⋯	0.86	500 s	17.82	1.00	2.78e-02
GRB 141102112	223.23	−17.42	9.21e+00	56963.112	⋯	0.02	15 days	18.56	2.97	2.81e-02

Note. This sub-population has a total of 134 GRBs. The variables are described in detail in Appendix A.1.

Download table as:  ASCIITypeset image

The most significant TW selected by this analysis is given in either seconds or days with units provided in the column. The TS and best-fit number of signal events (${\hat{n}}_{s}$) are provided as well as the p-value in the final column. The p-value has an effective trial correction for searching 10 TWs, but it is not corrected for the size of the given sub-population.

The GRBs are split into sub-populations by duration and hemisphere. Short GRBs are defined as T₉₀ < 2 s and long GRBs are T₉₀ ≥ 2 s. In this analysis, the northern sky is defined as a decl. greater than −5°, while the southern sky includes all decl. less than −5°. The most significant GRBs for each sub-population are shown in the following tables. The information is split into two sections: GRB information based on GCN notices and fit results based on the analysis.