Search for 10–1000 GeV Neutrinos from Gamma-Ray Bursts with IceCube

Gamma Ray Bursts, or GRBs, (Piran 2004) are among the most powerful objects in the Universe. During the milliseconds to hundreds of seconds "prompt phase," copious amounts of kilo electronvolt ro mega electronvolt photons are released. The fireball scenario (Piran 1999) assumes a radiation-dominated electron-positron plasma ejected in a jet at a high Lorentz factor, Γ ≳ 300. The prompt phase is followed by a decaying multiwavelength afterglow that can be observed from radio to X-rays for days to years (Kangas & Fruchter 2021). GRBs are empirically classified as either long for a duration of the prompt phase T₉₀, larger than 2 s, or as short otherwise. Long GRBs are associated with supernovae, see e.g., Hjorth (2013) and Cano et al. (2017) for reviews. Short GRBs are associated with the merger of compact objects as corroborated with the observation of gravitational waves for GRB/GW170817A (Abbott et al. 2017).

GRBs have been proposed as cosmic ray accelerators (Vietri 1995; Waxman 1995). The interaction of these cosmic rays with local environment radiation and/or matter would result in peta electronvolt neutrinos during the prompt phase (Waxman & Bahcall 1997; Zhang & Kumar 2013), tera electronvolt precursor neutrinos prior to the prompt phase (Mészáros & Waxman 2001; Razzaque et al. 2003; Murase & Ioka 2013), and exa electronvolt neutrinos during the afterglow (Waxman & Bahcall 2000). IceCube has discovered an all-sky, all-flavor, flux of neutrinos from ∼10 TeV to ∼10 PeV (Aartsen et al. 2020a; Abbasi et al. 2021, 2022a). GRBs are a long-standing candidate to explain, at least in part, this flux.

Neutrons and protons entrained on a GRB fireball can lead to gig-electronvolt-scale neutrinos (Bahcall & Mészáros 2000). Sub-photospheric neutron-proton decoupling in the fireball leads to 1–10 GeV neutrinos. These neutrinos are too low energy for the study reported here. Decoupled neutrons collide below the photosphere with fireball-entrained protons (Murase et al. 2013). These collisions produce charged and neutral pions. Charged pions decay, directly and indirectly, into neutrinos and charged leptons. In the collision scenario, typical neutrino energy is tens of giga electronvolts, an energy range that is tested in this work. For both decoupling and collision scenarios, neutrinos have a quasi-thermal spectrum and precede the prompt phase by tens of seconds.

To date, there is no evidence of neutrino emission from GRBs. IceCube has conducted searches for tera electronvolt to peta electronvolt neutrinos in coincidence with GRBs during the prompt phase (Abbasi et al. 2012a; Aartsen et al. 2015, 2016, 2017a). Recently, IceCube has searched for neutrino-GRB coincidences in the tera electronvolt to peta electronvolt energy range during the prompt, precursor, and afterglow phases (Abbasi et al. 2022b). During the prompt phase, IceCube results demonstrated that GRBs cannot be responsible for more than ∼1% of IceCube's extragalactic neutrino flux. For time windows around GRBs of 10⁴ s, centered on the prompt phase, GRBs cannot contribute to more than 24% of the extragalactic neutrino flux. ANTARES, which operated in the tera electronvolt to peta electronvolt energy range, conducted searches for neutrinos in coincidence with GRBs without finding a coincidence (Albert et al. 2017, 2021a, 2021b). Super-Kamiokande has conducted a search for neutrinos from GRBs above 8 MeV also with null results (Orii et al. 2021).

The brightest GRB of all time, GRB 221009A (Burns et al. 2023), was detected while this publication was being prepared and is not included in the list of 2268 GRBs we analyze here. Murase et al. (2022) have calculated neutrino fluxes for GRB 221009A in the sub-photospheric neutron-proton collision model. GRB 221009A was studied in neutrinos from mega electronvolts to peta electronvolts by IceCube (Abbasi et al. 2023a, 2023b) and no evidence for neutrino emission was found. The search for 10–1000 GeV neutrinos from GRB 221009A by IceCube is methodologically identical to the work presented here, differing only in the duration of time windows used. GRB 221009A was also studied by KM3NeT and no evidence for neutrino emission was found either (KM3NeT Collaboration 2022).

In this work, we present a study of 2268 GRBs detected by satellite-borne instruments with 8 yr of IceCube-DeepCore data in the 10–1000 GeV energy range. The present study covers the precursor, prompt, and early afterglow phases of GRBs for a total duration of up to 500 s centered on the prompt phase. We find no evidence for 10–1000 GeV neutrino emission by GRBs. We use two analysis methods, one to search for the most significant GRB-neutrino coincidence and another to search for an ensemble of GRBs that may be too weak to be detectable individually with neutrinos but may be significant as a group. We find no evidence of neutrino emission for these 2268 GRBs in either analysis and we set limits on the time-integrated, all-neutrino-flavor, neutrino flux. The results presented here are compared to limits for approximately tera electronvolts or greater that IceCube has published. Because GRB 221009A has an extremely high energy fluence, we compare predicted signal expectations for the sub-photospheric neutron-proton collision model for 2264 GRBs, for which energy fluence has been reported, to GRB 221009A. We find that, under the model assumptions, the signal expectation for GRB 221009A is a factor of 6–8 higher than for the other 2264 GRBs combined. Also, the total background for GRB 221009A is significantly lower. Thus, we derive the best possible limit on jet baryon loading using neutrino limits on GRB 221009A by IceCube (Abbasi et al. 2023b). Assuming a jet Lorentz factor of 300 (800), the baryon loading on GRB 221009A is lower than 3.85 (2.13) at a 90% confidence level.

The IceCube Neutrino Observatory (Aartsen et al. 2017b) consists of an array of 5160 digital optical modules (DOMs) on a total of 86 strings embedded within 1 km³ of Antarctic ice at the South Pole. All DOMs include a downward-facing photomultiplier tube (PMT; Aartsen et al. 2020b), and associated electronics (Abbasi et al. 2009) enclosed within a glass vessel and are deployed from 1.45 to 2.45 km below the surface. IceCube is optimized for approximately tera electronvolt or greater observations, matching an inter-string spacing of ∼125 m.

Six DeepCore strings were installed on the vertices of a hexagon with a side of 42 m. At the hexagon center is the central standard IceCube string. Two additional DeepCore-infill strings have been placed inside the hexagon with even smaller horizontal separation. The combination of these eight strings and seven nearby standard IceCube strings form the DeepCore sub-detector (Abbasi et al. 2012b).

The physics region of DeepCore is a cylinder of 125 m in radius and 350 m in height. DeepCore is optimized for ≳10 GeV neutrino observations. For DeepCore and DeepCore-infill strings, 50 DOMs are installed with 7 m vertical spacing, between 2.1 and 2.45 km of depth. The other 10 DOMs are installed with a 10 m vertical separation between 1.8 and 1.9 km of depth. These DOMs are used for enhanced down-going cosmic ray muon rejection. The region between 2000 and 2100 m is not instrumented in DeepCore and DeepCore-infill strings as glacial ice in this region has worse optical properties (Aartsen et al. 2013). All the DOMs in the DeepCore strings are equipped with high quantum-efficiency PMTs. The DOMs on the DeepCore-infill strings have a mixture of standard and high quantum-efficiency PMTs.

We use IceCube-DeepCore data collected between 2012 April 26, and 2020 May 29, with a total live time of 7.68 yr, which corresponds to an uptime of 95%. The data sample used in this publication, called GRECO-Astronomy (GeV Reconstructed Events with Containment for Oscillations), is described in detail in Abbasi et al. (2022c). GRECO-Astronomy has sensitivity over the entire sky (4π sr) that is only weakly dependent on decl. As they produce virtually identical signatures, IceCube-DeepCore cannot distinguish ν from $\bar{\nu }$. The GRECO-Astronomy data set is sensitive to all neutrino flavors via cascades and starting tracks.

Cascades, or showers, are the product of neutral current interactions for all neutrino flavors, as well as charged-current interactions of ν_e and the majority of ν_τ . For cascades, all energy is deposited in a small volume that is contained in DeepCore. This leads to good energy resolution, but poor angular resolution (Abbasi et al. 2022c). Starting tracks are predominantly muons originating from ν_μ charged-current interactions. In this case, part of the energy is deposited in a shower at the neutrino interaction vertex, plus an outgoing, relatively long-range muon. If sufficiently long, tracks can provide smaller angular uncertainty estimates than cascades. The GRECO-Astronomy data set includes event reconstructions for both cascades and starting tracks. The median angular resolution ranges from ∼40° for cascades with reconstructed visible energy of 10 GeV to a few degrees for starting tracks with ≳100 GeV reconstructed energy. In this work, we keep only events with reconstructed energy above 10 GeV, which reduces the average GRECO-Astronomy event rate from 4.6 mHz (16.6 events per hour) to 4.17 mHz (15 events per hour). This choice is taken because the angular uncertainty of lower reconstructed-energy events is judged too poor to be used here. There are 1,010,151 events in the final data sample used in this work.

We used GRBweb (Coppin 2021) to select GRBs for this study. This catalog brings together information from a variety of public GRB databases, such as the General Coordinate Network (GCN) Notices and Circulars, and compiles electromagnetic observational data from a large number of instruments including Fermi-GBM (Meegan et al. 2009), Fermi-LAT (Atwood et al. 2009), Swift-BAT (Barthelmy et al. 2005), the IPN network (Hurley et al. 2013), etc. Within the time period studied, 2297 GRBs were recorded in GRBweb. Of these, we have selected 2268 GRBs for this work. We have excluded 29 GRBs that do not have a sky localization or localization uncertainty. For each GRB, we obtain four pieces of information from GRBweb: the sky localization, the localization uncertainty, the duration, and the energy fluence. When multiple detectors observe a given burst, we use the localization information for the detector with the smallest localization uncertainty. In GRBweb, localization uncertainty is assumed to be the 1σ uncertainty for a 2D normal distribution, implying that the true position of the GRB lies in the uncertainty circle 39% of the time. Some instruments, e.g., Fermi-GBM provide more detailed information than the 1σ localization uncertainty. For these cases, GRBweb calculates the equivalent 1σ circularized uncertainty that provides the correct coverage. For GRBs localized by IPN, it is not possible to calculate 1σ circularized uncertainty since these GRBs are located in a box in the sky and only 3σ bounds are provided. For IPN, we use the larger of the two uncertainty box dimensions to generate a circularized uncertainty region. Each instrument that observes a GRB will typically provide a start time, T₀, and duration of the prompt phase, T₉₀. GRBweb determines the largest duration that covers all the T₉₀ observations. We call this duration T_prompt. In the majority of bursts used in this work, T_prompt matches the T₉₀ of a single observing instrument, which is most frequently Fermi-GBM. The energy fluence f_γ , during the prompt phase, is chosen by GRBweb to be that of the instrument with the widest energy observation band. Often this will be Fermi-GBM, as it has sensitivity from 10 keV to 10 MeV. The energy fluence is not used in the search for neutrino-GRB coincidences, but it is used in the interpretation of the results.

Each burst location and localization uncertainty is considered as a prior described by a probability density function. When the best localization of a given burst is by Fermi-GBM, we use HEALPix maps (Górski et al. 2005) produced by Fermi-GBM. Fermi-GBM maps starting in early 2018 were released publicly (Goldstein et al. 2020). Maps for bursts prior to 2018 were processed similarly using the same GBM Data Tools. However, the metadata in the files have not been fully qualified and the files have not been uploaded to HEASARC ⁶⁹ . Therefore, we use preliminary maps by Goldstein & Wood (2022). If a different instrument provides a better localization than Fermi-GBM, then synthetic HEALpix maps were created. These synthetic maps have uniform (top-hat) probability over the circularized 1σ localization uncertainty obtained via GRBweb. If the GRBweb localization uncertainty is smaller than 1°, then the synthetic map is created with a radius of 1°. As can be seen on the left panel of Figure 1, this lower bound of 1° is used in the vast majority of synthetic maps. The lower bound of 1° is used for computational reasons and does not significantly affect this analysis, as the typical GRECO-Astronomy event has an angular uncertainty much larger than 1°.

Zoom In Zoom Out Reset image size

Of the 2268 GRBs studied in this work, 1390 have Fermi-GBM HEALPix skymaps and 878 have synthetic skymaps. All HEALPix maps have been rebinned to a setting of Nside = 64, corresponding to 49,152 pixels over the entire sky. Each pixel has an angular area of ∼0.84 square degrees, which is usually smaller than the angular resolution for the best-reconstructed events in the GRECO-Astronomy data set.

We have conducted two analyses in this work. First, we search for the statistically most significant GRB in a temporal and directional coincidence with GRECO-Astronomy neutrino events. Second, we search for the most statistically significant group of GRBs in temporal and directional correlation with GRECO-Astronomy neutrinos. This latter search can potentially find a subset of the 2268 GRBs as neutrino emitters, even though none of the GRBs were statistically significant on their own. This second analysis uses a binomial test to statistically combine the results of the first analysis.

3.1. Search for the Most Significant Neutrino-GRB Correlation

Each GRB is studied using six predefined time windows with durations of 10, 25, 50, 100, 250, and 500 s centered in the middle of the T_prompt window. These overlapping time windows cover or over-cover the vast majority of prompt phases of the 2268 GRBs in our list. Longer-than-prompt time windows are used to identify potential precursor and early afterglow neutrino-GRB correlations. The right panel of Figure 1 shows the distribution of T_prompt for all the GRBs used in this work, as well as the duration of the time windows that we study. To characterize GRB-neutrino emission, we use the time-integrated neutrino number flux:

$\begin{eqnarray}&&F(E)={\int }_{{t}_{a}}^{{t}_{b}}\displaystyle \frac{{{dN}}_{\nu +\bar{\nu }}}{{dEdAdt}}{dt},\end{eqnarray} \tag{ 1 }$

where t_a and t_b set the time window being considered and ${{dN}}_{\nu +\bar{\nu }}/{dEdAdt}$ is the number of neutrinos per unit of area, energy, and time. Equal values of F(E) can be obtained for arbitrary neutrino time profiles within t_a and t_b .

The unbinned maximum likelihood method is a common approach to search for time-dependent neutrino sources (Braun et al. 2008). In this work, we use an extended likelihood, which is a function of the number of neutrino signal events, n_s , the location of the sky that is being studied, Ω, conditional on the expected number of background events, n_b , and per-neutrino information x _i . The neutrino information x _i is reconstructed direction (R.A., α_i and decl., δ_i ), neutrino angular uncertainty estimator σ_i , neutrino arrival time t_i , and reconstructed energy E_i . The extended likelihood function is:

$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L }({\boldsymbol{\Omega }},{n}_{s}| {n}_{b},\left\{{{\boldsymbol{x}}}_{i}\right\}) & = & \displaystyle \frac{{\left({n}_{s}+{n}_{b}\right)}^{N}{e}^{-({n}_{s}+{n}_{b})}}{N!}\\ & & \times {\prod }_{i=1}^{N}\left(\displaystyle \frac{{n}_{s}}{{n}_{s}+{n}_{b}}{ \mathcal S }({{\boldsymbol{x}}}_{i},{\boldsymbol{\Omega }})+\displaystyle \frac{{n}_{b}}{{n}_{s}+{n}_{b}}{ \mathcal B }({{\boldsymbol{x}}}_{i},{\boldsymbol{\Omega }})\right),\end{array}\end{eqnarray} \tag{ 2 }$

with the index i iterating over all N candidate neutrino events in a given time window. The Poisson term is frequently used by IceCube in transient neutrino analyses, including GRB studies (Aartsen et al. 2015, 2016, 2017a; Abbasi et al. 2022b), in which n_s + n_b is a relatively small number. ${ \mathcal S }(\cdot )$ and ${ \mathcal B }(\cdot )$ are the signal and background probability density functions (PDFs) respectively. Both ${ \mathcal S }(\cdot )$ and ${ \mathcal B }(\cdot )$ are described as the product of a directional term and an energy term.

The directional term of ${ \mathcal S }(\cdot )$ is represented by a Kent function (Kent 1982)

$\begin{eqnarray}&&{S}_{\mathrm{space}}({\rm{\Delta }}{\psi }_{i})=\displaystyle \frac{\kappa }{4\pi \sinh \kappa }{e}^{\kappa \cos {\rm{\Delta }}{\psi }_{i}},\end{eqnarray} \tag{ 3 }$

where Δψ_i is the angular difference between each reconstructed event direction and the direction, Ω, being studied and $\kappa =1/{\sigma }_{i}^{2}$. The choice of a Kent function, instead of a normal distribution is made because the angular uncertainty of GRECO-Astronomy events is relatively large. The background space PDF, B_i,space, is a function of zenith only due to the approximate azimuthal symmetry of the IceCube-DeepCore detector, and is determined from 7.68 yr of data, but excluding events from the corresponding time window around all of the GRBs. Background events in GRECO-Astronomy are ∼60% due to atmospheric neutrinos, and ∼40% due to down-going atmospheric muons (Abbasi et al. 2022c).

The energy term of the signal PDF is determined from simulations that assume an energy spectrum of E^−2.5. This choice of spectrum was made so that the most sensitive energy range of the study would correspond to approximately 20–500 GeV. The energy term of the background PDF is determined from 7.68 yr of data, again, excluding events from the time window around all of the GRBs.

For each of the six time windows of a GRB, we perform a maximum likelihood fit over the entire sky, placing Ω at the center of each bin of the Nside = 64 HEALPix map. The maximization is done for n_s only, and we obtain a best-fit value ${\hat{n}}_{s}$ for each Ω and time window. For each GRB time window we define the sky-map Λ_ν (Ω):

$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{\nu }({\boldsymbol{\Omega }}) & = & 2\cdot \mathrm{log}\left[\displaystyle \frac{{ \mathcal L }({\hat{n}}_{s},{\boldsymbol{\Omega }})}{{ \mathcal L }({n}_{s}=0,{\boldsymbol{\Omega }})}\right]\\ & = & -2{\hat{n}}_{s}+2{\sum }_{i=1}^{N}\mathrm{log}\left[\displaystyle \frac{{\hat{n}}_{s}S({{\boldsymbol{x}}}_{i},{\boldsymbol{\Omega }})}{{n}_{b}B({{\boldsymbol{x}}}_{i},{\boldsymbol{\Omega }})}+1\right].\end{array}\end{eqnarray} \tag{ 4 }$

Note that the sky-map Λ_ν (Ω) depends only on neutrino data. Up to this point, we have not used GRB positions or localization uncertainties. To include this information, and similarly to Abbasi et al. (2022b), we use the previously described Fermi-GBM HEALPix maps, or synthetic maps, as a prior on the final test statistic. We define a weight $w({\boldsymbol{\Omega }})=P({\boldsymbol{\Omega }})/{P}_{\max }$ for each GRB map. Here P(Ω) is the GRB localization probability evaluated at each sky location Ω, and P_max is the maximum probability for the map of this GRB. The combined sky map of neutrino and GRB information, Λ_final(Ω) is defined as:

$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Lambda }}}_{\mathrm{final}}({\boldsymbol{\Omega }})=2\cdot \mathrm{log}\left[w({\boldsymbol{\Omega }})\displaystyle \frac{{ \mathcal L }({\hat{n}}_{s},{\boldsymbol{\Omega }})}{{{ \mathcal L }}_{\nu }({n}_{s}=0,{\boldsymbol{\Omega }})}\right]\\ \,=\,{{\rm{\Lambda }}}_{\nu }({\boldsymbol{\Omega }})+2\mathrm{log}(w({\boldsymbol{\Omega }})).\end{array}\end{eqnarray} \tag{ 5 }$

This all-sky map has a maximum value at $\hat{{\boldsymbol{\Omega }}}$, which for each GRB and each time window, defines the test statistic (TS) used in this work,

$\begin{eqnarray}&&{TS}={{\rm{\Lambda }}}_{\mathrm{final}}(\hat{{\boldsymbol{\Omega }}}).\end{eqnarray} \tag{ 6 }$

To estimate the significance of TS, we use a large number of scramblings. This is a standard technique in IceCube, in which the time of each event is randomized within the live time of IceCube being studied. This procedure keeps the reconstructed neutrino direction detector coordinates constant but randomizes the reconstructed neutrino R.A. Given the very large duration of the live time compared to the total duration of all the time windows that we study, each scrambling data set is a good approximation of the background expectation. In this work, we have used 1000,000 scramblings.

The significance of correlations between neutrino events and each GRB for a given time window is quantified as a p-value by comparing the measured TS to the background-only TS distribution. The background-only TS distribution is obtained via scramblings.

At this stage, we have a total of six p-values for each GRB, corresponding to the six time windows that have been studied. For each GRB we choose the most significant of the six p-values. We calculate a trial correction for the look-elsewhere effect for this choice by using the scrambled data cumulative distribution function of the p-value of the most significant time window. For each GRB we now have the most significant time window for correlation and the best-fit ${\hat{n}}_{s}$ for this time window and a single per-GRB p-value.

Finally, the most significant per-GRB p-value is selected and trial-corrected to account for 2268 GRBs in the search. This results in the most significant GRB, its most significant time window, the best-fit ${\hat{n}}_{s}$ for this GRB time window, and the posttrial p-value.

To characterize the performance of this analysis, we use sensitivity and 5σ discovery potential. With repeated simulated signal injected over background scramblings, we can calculate the signal TS distribution. The sensitivity is the signal strength needed so that 90% of the signal-injected TS is above the median TS for background scramblings. The 5σ discovery potential is defined as the signal strength required for the injected median TS to be greater than the 5σ threshold of the background scrambling TS distribution.

Figure 2 shows the 90% sensitivity and 5σ discovery potential for all six time windows for the most significant GRB identified in this analysis, bn 140807500. Neutrino acceptance and background rates vary slowly with decl. in GRECO-Astronomy. This results in sensitivity and 5σdiscovery potential that vary by, at most, a factor of ∼2 from the northern sky (best) to the southern sky (worst) (Abbasi et al. 2022c). This is unlike GRB studies with TeV–PeV neutrinos in which sensitivity and 5σ discovery potential for northern declinations is a factor of ∼20 better than for southern declinations (Abbasi et al. 2022b). At the South Pole, the intense down-going cosmic ray muon background worsens southern sky sensitivity for tera electronvolt to peta electronvolt neutrinos compared to the northern sky. As the background rate in the work presented here is very low, sensitivity and 5σ discovery potential do not depend strongly on the GRB localization uncertainty. This is because the angular uncertainty of GRECO-Astronomy events is typically larger than GRB localization uncertainty. For example, a GRB localized to 15° has ∼15% higher background rate than one localized to 1°. GRB localization uncertainty has a smaller influence on sensitivity and 5σ discovery potential than declination.

Zoom In Zoom Out Reset image size

3.2. Binomial Test

We statistically combine per-GRB p-values for the most significant time window, to search for a subgroup of GRBs that may be significant neutrino emitters. This case can be interesting when each individual GRB is not statistically significant by itself. The binomial probability, which has been used by IceCube in other works, e.g., Aartsen et al. (2020c), is given by:

$\begin{eqnarray}&&P(k)={\sum }_{m=k}^{N}\displaystyle \frac{N!}{(N-m)!m!}{{p}_{k}}^{m}{\left(1-{p}_{k}\right)}^{N-m}.\end{eqnarray} \tag{ 7 }$

Here the pretrial binomial p-value P(k) denotes the probability of k or more GRBs with p-values smaller than p_k appearing among the background, where p_k is the kth smallest pretrial p-value in the final GRB list. Taking into account the number of independent GRBs (2268), we adjust the probability of observing a given result by chance with the cumulative density function (CDF) made out of the most significant binomial p-values obtained from 100 million null hypothesis binomial tests. The smallest pretrial binomial p-value P(k) obtained at a size of subgroup k is corrected for trials and reported as posttrial binomial p-value P_binom.

3.3. Systematic Uncertainties

The most significant GRB-neutrino correlation identified among 2268 GRBs is the Fermi-GBM GRB bn 140807500. The p-value, corrected for the six time windows, is 4.6 × 10⁻⁵ and is found for a time window of 100 s with the best-fit number of events ${\hat{n}}_{s}=1.08$. After correcting for trials for 2268 GRBs, we obtain a posttrial p-value of 0.097 (1.3σ). Burst bn 140807500 triggered Fermi-GBM on 2014 August 7. This is a short GRB with T₉₀ = 0.512 ± 0.202 s and an energy fluence of (1.289 ± 0.014) × 10⁻⁶ erg cm⁻² (Coppin 2021). Because only Fermi-GBM identified this burst, the T₉₀ matches T_prompt.

This analysis uses the Fermi-GBM maps whenever they are available, where the most likely location of this GRB is R.A. α = 1983 and decl. δ = 317. Two neutrino candidates are identified in IceCube-DeepCore in the 100 s time window around the center of T_prompt for bn 140807500. Event one is identified as a cascade detected 33.67 s before the center of the T_prompt window. This cascade has a reconstructed energy of 50 GeV. The best-fit reconstructed direction is α = 768 and δ = − 515, which is 1335 away from the most likely GRB location. Event one has a circularized directional uncertainty of 354. Event two is a starting track detected 38.58 s after the center of the T_prompt window. The reconstructed energy is 221 GeV. The best-fit reconstructed direction is α = 2002 and δ = 333, which is 23 away from the most likely GRB location. Event two has a circularized directional uncertainty of 10. Figure 3 shows the localization probability sky map of bn 140807500 provided by Fermi-GBM, the Λ_ν all-sky scan from neutrino events within the 100 s time window, and the Λ_final sky map. Figure 3 shows that the second neutrino candidate is responsible for a relatively high test statistic. Table 1 shows the three most significant GRB-neutrino coincidences identified.

Zoom In Zoom Out Reset image size

Table 1. Top Three Most Significant GRBs in the Individual Search Analysis

GRB	Start Time	R.A.	Decl.	Loc.	Loc.	TW	${\hat{n}}_{s}$	E2 F(E) Upper Limit		Synth.	Pretrial	Posttrial
Name				Unc.	Satellite			γ = 2.5∣100GeV	γ = 2.0	Prior	p-value	p-value
	(UTC)	(deg)	(deg)	(deg)		(s)		(10−5erg cm−2 s−1)
bn 140807500	2014-08-07 11:59:33	198.3	31.7	4.4	GBM	100	1.08	6.42	2.62	No	4.6 × 10−6	0.097
bn 160804968	2016-08-04 23:13:34	82.3	−23.0	9.1	GBM	100	2.81	11.52	5.97	No	9.6 × 10−4	⋯
GRB 160802A	2016-08-02 06:13:30	28.1	71.4	1.0	IPN	50	2.74	8.46	4.13	Yes	1.3 × 10−3	⋯

Note. The table shows the GRB name; start time of T_prompt; R.A.; decl.; circularized localization uncertainty; most significant time window duration; best-fit ${\hat{n}}_{s}$ number of neutrino events in coincidence with the GRB; all-flavor neutrino number flux F(E) upper limit times E² for the most significant time window TW for two spectral indices; whether or not a synthetic prior was used; the GRB pretrial p-value, corrected for the six time windows used in this analysis; and the posttrial p-value for the most significant burst. The R.A. and decl. for bn 140708500 and bn 160804968are for P_max of the GBM map.

Download table as:  ASCIITypeset image

The smallest binomial p-value is P(k = 1) = 0.099 with index k = 1 with threshold p-value p_k = 4.6 × 10⁻⁵. After correcting for trials due to testing multiple thresholds, the posttrial binomial p-value of 0.65 is obtained. So no additional neutrino event-GRB correlations, besides bn 140807500, are identified with the binomial test.

We have followed calculations by Murase et al. (2013) and Murase et al. (2022), for the neutron-proton collision scenario, to obtain an all-flavor neutrino signal expectation for the GRBs used in this work. The authors of these works provided us with neutrino spectra that allow the calculation of E² F(E) (Murase 2022). Given a GRB, the model depends on the gamma-ray energy fluence of the prompt phase f_γ , redshift z, and the product of baryon loading times the neutron-proton opacity ξ_N τ_np. Following Murase et al. (2013), we adopt the definition of the baryon loading as the ratio of the equivalent isotropic energies in baryons to gamma-rays, ξ_N = E_N,iso/E_γ,iso, with ξ_N = 5 and τ_np = 1 as default model values.

Of the 2268 GRBs studied here, four do not have a measured energy fluence in GRBWeb so we exclude these GRBs from the calculation. Only 198 GRBs have a measured redshift, of which 187 are long and 11 are short GRBs. Among these, the mean redshift for long GRBs is 2.0 and, for short GRBs is 0.7. In our calculations for GRBs without a measured redshift we adopt a default value of z = 0.7 for short GRBs and z = 2.0 for long GRBs. With a Lorentz factor of Γ = 300 (800), and default ξ_N τ_np = 5 we obtain an all-flavor neutrino signal, using GRECO-Astronomy simulation with E_reco > 10 GeV, in IceCube-DeepCore of 0.64 (1.39) events for all 2264 GRBs combined. As a systematic check, we have redone the calculation for 2264 GRBs by assuming z = 0.25 and z = 1.0 for short and long GRBs respectively, that do not have a measured redshift. We find for this check a signal expectation of 0.64 (1.24). Therefore we find that the lack of redshift information for most GRBs does not significantly affect the neutrino signal prediction for the sub-photospheric model.

The all-flavor signal expectation for 2264 GRBs can be compared to that for GRB 221009A. For the latter, we use z = 0.151 (de Ugarte Postigo et al. 2022), an energy fluence measured by Konus-Wind of f_γ = 1.2 × 10⁵⁵ erg (Frederiks et al. 2023), ξ_N τ_np = 5, and Lorentz factor of 300 (800). This results in an expectation of 4.72 (8.56) events. Alternatively, using the energy fluence measured by Fermi-GBM, f_γ = 1.0 × 10⁵⁵ erg (Lesage et al. 2023), results in a factor of 1.2× lower signal expectation, as F(E) is proportional to f_γ .

Remarkably, GRB 221009A has a larger neutrino signal expectation in IceCube-DeepCore than the other 2264 GRBs combined, by a factor of 6–8. Even if we assume redshift values for bursts without a redshift measurement that are significantly smaller than typical, we still find that GRB 221009A has a larger signal expectation over the other 2264 GRBs combined by a similar factor. A comparison of E² F(E) for GRB 221009A and 2264 GRBs combined is shown in Figure 4. There are two reasons why GRB 221009A has a larger neutrino signal expectation than the combination of the GRBs studied here. First, neutrinos of higher energy are easier to identify and correlate with GRBs because they typically have lower angular uncertainty and because the detector has a higher efficiency for detecting them. The peak energy of E² F(E), seen in Figure 4, is inversely proportional to (1 + z) and the redshift of GRB 221009A is small. So GRB 221009A has higher energy neutrinos than a GRB at z = 2. Second, GRB 221009A has an extremely high energy fluence, and the height of the peak of E² F(E) seen in Figure 4 is proportional to f_γ .

Zoom In Zoom Out Reset image size

The parameter product ξ_N τ_np is unknown for all GRBs. We assume τ_np = 1, which leaves the baryon loading unknown. Knowledge of the baryon loading can provide information about the environment for the formation of the jet. The best constraint on ξ_N τ_np can be derived from GRB 221009A for two reasons: GRB 221009A's high signal expectation compared to the set of 2264 GRBs and the lack of correlated neutrino events reported in Abbasi et al. (2023a) and because, for 2264 GRBs, with a per-GRB window of 500 s, background accumulates over 1,132,000 s (13 days), while for GRB 221009A, the background is only over a single time window of 2200 s. For GRB 221009A, we calculate the constraint on ξ_N using the 2200 s time window, as this includes observations prior to the prompt phase. Stating that the constraint from GRB 221009A is the most restrictive, assumes that GRB 221009A has similar characteristics as all other GRBs. This assumption may not be valid, e.g., because of the unusually high isotropic equivalent energy of GRB 221009A. Therefore there is still important, complementary information provided by the neutrino flux constraints we have placed on the other 2268 GRBs.

With a Lorentz factor of 300 (800) and a time window of 2200 s, the upper limit on E² F(E) for GRB 221009A is a factor of 1.65 (2.98), lower than predicted by the model in Murase et al. (2022), but updated to use the Konus-Wind energy fluence. This upper limit directly translates into a constraint of ξ_N τ_np < 3.85 (2.13) for a Lorentz factor of 300 (800). In Figure 5, the left panel shows the canonical model prediction for GRB 221009A and GRECO-Astronomy limits from Abbasi et al. (2023b) for a Lorentz factor of Γ = 400. The right panel shows the constraint on ξ_N τ_np as a function of the Lorentz factor.

Zoom In Zoom Out Reset image size

Lesage et al. (2023) provide lower bounds on the Lorentz factor of GRB 221009A under various modeling scenarios. Requiring no photon–photon e⁺ e⁻ pair production for the highest energy photons observed in Fermi-GBM and using a single zone model, they set a constraint of Γ ≳ 1560. However, under more realistic scenarios, the lower bound is a factor of 2 lower or Γ ≳ 780. Figure 5 is agnostic to the choice of Lorentz factor, but taking a value of Γ = 780, results in a constraint on the baryon loading that is significantly lower than the canonical theoretical value.

The constraint on baryon loading we present here is complementary to that which can be set, also for GRB 221009A, using tera electronvolt to peta electronvolt neutrinos. With 10–1000 GeV neutrinos, the limit on ξ_N τ_np is more constraining as the Lorentz factor increases. This is because the peak value of E² F(E) does not depend on the Lorentz factor but the energy for which E² F(E) peaks depends linearly on the Lorentz factor. A higher value of the Lorentz factor results, on average, on higher energy neutrinos, which are easier to detect and correlate to a GRB.

The situation is reversed for tera electronvolt to peta electronvolt neutrinos in which the constraint on the baryon loading is best for low values of the Lorentz factor. For both the internal shock model and the Internal-Collision-Induced Magnetic Reconnection and Turbulence, ICMART (Zhang & Yan 2011; Zhang & Kumar 2013) model which predict tera electronvolt to peta electronvolt neutrinos, the peak value of E² F(E) is lower for larger values of the Lorentz factor. Using tera electronvolt to peta electronvolt neutrinos, the 90% confidence level upper limit on the baryon loading, derived from GRB 221009A, under the internal shock model and for Γ = 300, is ξ_N < 0.55. Under the same assumptions, but for the ICMART model, the limit is ξ_N < 2.97 (Abbasi et al. 2023b).

Figure 6 shows a comparison of limits on E² F(E) for several GRBs set by IceCube using GRECO-Astronomy and tera electronvolt to peta electronvolt searches. It is worth noticing that the recent study of correlations and GRBs for approximately tera electronvolt or greater neutrinos (Abbasi et al. 2022b) also found a correlation between bn 140807500 and event two (see Section 4). In that work, bn 140807500 was the most likely GRB-neutrino correlation among short GRBs in the northern sky (which IceCube defines as δ > − 5°). The data set used on Abbasi et al. (2022b) was collected from 2012 April to 2018 October. During this data period, 2.9% of the events in the data used in the current work are also found in Abbasi et al. (2022b). Because event two is an event that starts in DeepCore and has relatively high energy, it is not surprising that both analyses identify it.

Zoom In Zoom Out Reset image size

Using IceCube data and public GRB data, we have studied 10–1000 GeV neutrino and GRB correlations for 2268 GRBs detected over 8 yr. No evidence for neutrino emission by GRBs is found using either of the two analysis methods. In the first method, we search for the most statistically significant GRB-neutrino correlation. The most significant GRB is Fermi-GBM bn 140807500 with a posttrial p-value of 0.097. In the second method, we statistically combined the results for all 2268 GRBs to search for a set of GRBs that could be significant as a group but not individually. We do not find any additional burst, besides bn 140807500, to possibly contribute, and the p-value for this test is 0.65.

We compare sub-photospheric model predictions for a subset of 2264 GRBs, with prompt gamma-ray energy-fluence measurements, to the sub-photospheric model prediction of the Brightest of all time GRB 221009A. We find that, in the sub-photospheric model, GRB 221009A results in a neutrino signal in IceCube-DeepCore that is ≳6 larger than for the combined set of 2264 GRBs. We use previously calculated limits on neutrino emission for GRB 221009A to constrain the baryon loading of the jet. For a Lorentz factor of 300 (800), the baryon loading on GRB 221009A is lower than 3.85 (2.13) at a 90% confidence level. The set of 2268 GRBs may still be useful to constrain models besides the sub-photospheric model.

While GRECO-Astronomy cannot reach the sensitivity that the TeV–PeV searches achieve for northern sky GRBs, it covers a complementary energy range where different physical mechanisms of neutrino emission can be explored.

Future work that benefits from the use of the IceCube-Upgrade (Ishihara 2021) will enhance the sensitivity of IceCube-DeepCore to 10–1000 GeV neutrinos as the angular resolution of reconstructed events is expected to improve.

We are grateful to K. Murase for providing templates of neutrino spectra necessary to calculate neutrino fluxes for the sub-photospheric model under the collision scenario. We are grateful to J. Wood and A. Goldstein for suggestions that improved this manuscript. The authors gratefully acknowledge the support from the following agencies and institutions: USA—U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, U.S. National Science Foundation-Office of Advanced Cyberinfrastructure, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin-Madison, Open Science Grid (OSG), Partnership to Advance Throughput Computing (PATh), Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS), Frontera computing project at the Texas Advanced Computing Center, U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, Astroparticle physics computational facility at Marquette University, NVIDIA Corporation, and Google Cloud Platform; Belgium—Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany—Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden—Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; European Union—EGI Advanced Computing for Research; Australia—Australian Research Council; Canada—Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Digital Research Alliance of Canada; Denmark—Villum Fonden, Carlsberg Foundation, and European Commission; New Zealand—Marsden Fund; Japan—Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea—National Research Foundation of Korea (NRF); Switzerland—Swiss National Science Foundation (SNSF).