SEARCH FOR MUON NEUTRINOS FROM GAMMA-RAY BURSTS WITH THE IceCube NEUTRINO TELESCOPE

Gamma-ray bursts (GRBs) are among the most violent events in the universe and among the few plausible candidates for sources of the ultra-high-energy cosmic rays. The so-called long-duration GRBs (≳2 s) are thought to originate from the collapse of a massive star into a black hole (Woosley 1993), whereas short-duration GRBs (≲2 s) are believed to be the result of the merger of two compact objects (e.g., neutron stars) into a black hole (Eichler et al. 1989). Though quite different in nature both scenarios are consistent with the currently leading model for GRBs, the fireball model (Meszaros & Rees 1993), with the energy source (central engine) being the rapid accretion of a large mass onto the newly formed black hole. In this model, a highly relativistic outflow (fireball) dissipates its energy via synchrotron or inverse Compton radiation of electrons accelerated in internal shock fronts (Narayan et al. 1992; Rees & Meszaros 1994; Sari & Piran 1997). This radiation in the keV–MeV range is observed as the γ-ray signal. In case of long GRBs, the energy in gamma rays is typically of ${\cal O}(10^{51}$–10⁵⁴ erg× Ω/4π) where Ω is the opening angle for the γ-ray emission. Short GRBs are observed to release about a factor 100 less energy.

In addition to electrons, protons are thought to be accelerated via the Fermi mechanism, resulting in an E⁻² power-law spectrum with energies up to 10²⁰ eV (Waxman 1995; Vietri 1995). The normalization of the proton spectrum is usually given in relation to the energy in electrons. The latter is linked to the energy in γ-ray photons through the synchrotron and inverse Compton energy-loss mechanisms. Protons of ${\cal O}(10^{15}\,\mathrm{eV})$ interact with the keV–MeV photons forming a Δ⁺ resonance which decays into pions (Waxman & Bahcall 1997). In the decay of the charged pions, neutrinos of energy ${\cal O}(10^{14}\,\mathrm{eV})$ are produced with the approximate ratios (ν_e:ν_μ:ν_τ) = (1:2:0),⁴⁰ changing to about (1:1:1) at the Earth due to oscillations (Learned & Pakvasa 1995; Athar et al. 2006). First calculations of this prompt neutrino flux (Waxman & Bahcall 1997; Alvarez-Muniz & Halzen 1999) used average GRB parameters and the GRB rate measured by BATSE to determine an all-sky neutrino flux from the GRB population. The AMANDA-II neutrino telescope (Achterberg et al. 2007, 2008) performed searches for this so-called Waxman–Bahcall GRB flux or similar GRB fluxes (Achterberg et al. 2007, 2008) with negative results.

In a similar way, the so-called precursor neutrinos can be generated when the expanding fireball is still inside the progenitor star (Razzaque et al. 2003). In this case, the accelerated protons interact with matter of the progenitor star or synchrotron photons. However, due to the large optical depth the synchrotron photons cannot escape the fireball and, hence, no γ-ray signal is observed. The time delay between the start of this neutrino emission and the prompt γ-ray signal is expected to be about 100 s.

Observations of the early and late afterglow phases reveal that a large fraction of GRBs show X-ray flares superposed on the decaying light curve. Sometimes these flares are interpreted as a restart of the central engine that already generated the prompt emission (Burrows et al. 2007). If this is true, neutrino production with a similar spectrum as the prompt emission can be expected in the afterglow phase up to 10⁴ s after the γ-ray signal (Murase & Nagataki 2006). Furthermore, production of neutrinos with energies around 10¹⁸ GeV is expected when the shock fronts collide with the interstellar medium or the progenitor wind (Waxman & Bahcall 2000).

In our analysis, we search for muon neutrinos from GRBs recorded by satellites between 2007 June 1 and 2008 April 4 in all three phases. For the prompt phase, we utilize both an unbinned likelihood and a binned method. We find that the unbinned likelihood method has a significantly better discovery potential and is therefore used to obtain the limits presented in this paper. For searches in other emission phases, we perform only an unbinned search. The paper is structured as follows. In Section 2, we define the neutrino spectra used for the different phases, followed by a description of the IceCube detector in Section 3. Afterward, in Sections 4 and 5 the data sets and simulations are discussed, respectively. In Section 6, the unbinned likelihood and binned methods are described and their performance is compared. Section 7 then presents the results followed by a discussion of systematic uncertainties. Finally, Section 8 sets the results into context with other observations.

The searches in this paper rely on the directional, temporal, and spectral information obtained from satellite-based γ-ray observations which are distributed via the Gamma-ray burst Coordinate Network (GCN⁴¹). Primarily, this information comes from Swift (Burrows et al. 2005; also X-ray and UV observations), but also from Konus-Wind⁴², SuperAGILE (Tavani et al. 2008), Integral (Mereghetti 2004), and other satellites of the Third Interplanetary Network⁴³. In our analyses, we only consider bursts which occurred at a declination above −5°. The southern sky is dominated by downgoing muons created by cosmic-ray interactions with the atmosphere. By restricting our searches to the northern sky, the background from downward going muons is drastically reduced. The resolution of the GRB position from the satellites is better than 01, well below the resolution of the IceCube detector. It is therefore neglected in these analyses.

2.1. Prompt Emission

We calculate the expected prompt neutrino spectrum in the internal shock scenario of the fireball model following Guetta et al. (2004) which is based on Waxman & Bahcall (1997). For reference, we list all formulae used in our calculations in the Appendix together with a definition of the various parameters. Prompt neutrino emission from GRBs is the result of meson production in collisions of accelerated protons and the observed γ-rays in the keV–MeV range. It is therefore expected to occur during the same time frame as the γ-ray emission and to track the photon energy spectrum. This is reflected in the similar functional form of F_γ (Equation (A1)) and F_ν (Equation (A3)). Due to the Δ⁺ resonance condition, the neutrino energy is predicted to be inversely proportional to the photon energy, which is illustrated in the definition of α_ν and β_ν in Equation (A6). The further break in the neutrino spectrum above an energy ₂ is due to synchrotron cooling of high-energy pions and muons before producing neutrinos. The expected energy fluence in neutrinos is directly proportional to the measured energy fluence in photons (Equation (A8)). Here, some of the measured photon indices lead to diverging γ-ray spectral integrals if integrated from zero to infinity in energy. As the photon spectrum will not follow a broken power-law spectrum to arbitrarily high or low energies, we limit the integration range for all GRBs from 1 keV to 10 MeV for which broken power-law spectra have been observed by γ-ray satellites.

The time window of the prompt emission is determined from the information published in the GCN circulars and reports and checked against the measured γ-ray emission curves (available from GCN⁴¹). In case of early or late emission outside the window specified by the satellite experiments the window is extended accordingly. The exact window definitions for the bursts are listed in Table 1. In previous publications on stacked searches for neutrinos from GRBs (Achterberg et al. 2007, 2008) it was assumed that the sum of all GRB spectra follows a Waxman–Bahcall GRB spectrum (Waxman & Bahcall 1997). However, the burst parameters can vary significantly from burst to burst (Guetta et al. 2004; Becker et al. 2006) and the GRB population used here (mostly observed by the Swift satellite) is different from the BATSE population on which the calculations from Waxman and Bahcall were based. Therefore, we use the measured parameters to calculate the neutrino spectra for each GRB individually as has already been done in the case of single, bright bursts (Stamatikos et al. 2005; Abbasi et al. 2009b). The parameters (see the Appendix) for each GRB are listed in Table 1. In case a parameter has not been measured for a particular burst the average value listed in Table 2 is used. The table also contains assumed values of parameters that have not been measured for any of the GRBs. The resulting neutrino spectra for all 41 GRBs are shown in Figure 1 together with the standard Waxman–Bahcall spectrum.⁴⁴ Clear differences in the shapes of the two summed spectra are observed together with an ${\cal O}(10)$ times lower overall fluence for the individual spectra. The latter is caused by the much higher sensitivity of Swift compared to BATSE. The observed differences stress the importance of using individual fluences in analyses.

Zoom In Zoom Out Reset image size

Table 1. Burst Parameters, γ-Ray, and Neutrino Spectra of All 41 GRBs (for Definitions of Parameters See the Appendix)

	General Burst Parameters						γ-Ray Spectrum				ν Spectrum
	T0	R.A.	Decl.	T1 − T0	T2 − T0	z	fγ	γ	αγ	βγ	fν	ν,1	ν,2	αν	βν	γν
GRB 070610	20:52:26	298.8	26.2	−0.8	+4.4	2.00a	2.3e+00	0.20a	1.76	2.76a	2.2e−17	0.35	8.54	0.24	1.24	3.24
GRB 070612A	02:38:45	121.4	37.3	−4.7	+418.0	2.00a	1.1e+02	0.20a	1.69	2.69a	1.2e−15	0.35	8.54	0.31	1.31	3.31
GRB 070616	16:29:33	32.2	56.9	−2.6	+602.2	2.00a	2.2e+02	0.20a	1.61	2.61a	2.8e−15	0.35	8.54	0.39	1.39	3.39
GRB 070704	20:05:57	354.7	66.3	−57.3	+400.8	2.00a	5.4e+01	0.20a	1.79	2.79a	4.9e−16	0.35	8.54	0.21	1.21	3.21
GRB 070714B	04:59:29	57.8	28.3	−0.8	+65.6	0.92	1.1e+01	0.20a	1.36	2.36a	5.1e−17	0.85	13.34	0.64	1.64	3.64
GRB 070724B	23:25:09	17.6	57.7	−2.0	+120.0	2.00a	1.1e+03	0.08	1.15	3.15	1.8e−15	0.85	8.54	−0.15	1.85	3.85
GRB 070808	18:28:00	6.8	1.2	−0.7	+41.4	2.00a	1.6e+01	0.20a	1.47	2.47a	2.8e−16	0.35	8.54	0.53	1.53	3.53
GRB 070810B	15:19:17	9.0	8.8	+0.0	+0.1	2.00a	1.4e−01	0.20a	1.44	2.44a	2.6e−18	0.35	8.54	0.56	1.56	3.56
GRB 070917	07:33:57	293.9	2.4	−0.1	+11.4	2.00a	3.7e+01	0.21	1.36	3.36	6.0e−16	0.33	8.54	−0.36	1.64	3.64
GRB 070920A	04:00:13	101.0	72.3	+15.1	+75.0	2.00a	5.3e+00	0.20a	1.69	2.69a	5.7e−17	0.35	8.54	0.31	1.31	3.31
GRB 071003	07:40:55	301.9	10.9	−7.6	+167.4	2.00a	3.0e+01	0.80	0.97	2.97	3.8e−14	0.09	8.54	0.03	2.03	4.03
GRB 071008	21:55:56	151.6	44.3	−11.0	+14.0	2.00a	1.2e+00	0.20a	2.23	3.23a	6.4e−18	0.35	8.54	−0.23	0.77	2.77
GRB 071010B	20:45:47	150.5	45.7	−35.7	+24.1	0.95	2.1e+03	0.03	1.25	2.65	5.0e−17	5.71	13.16	0.35	1.75	3.75
GRB 071010C	22:20:22	338.1	66.2	−2.0	+20.0	2.00a	3.2e+01a	0.20a	1.00a	2.00a	1.8e−15	0.35	8.54	1.00	2.00	4.00
GRB 071011	12:40:13	8.4	61.1	−9.5	+63.8	2.00a	3.2e+01	0.20a	1.41	2.41a	6.3e−16	0.35	8.54	0.59	1.59	3.59
GRB 071013	12:09:19	279.5	33.9	−5.9	+23.4	2.00a	3.7e+00	0.20a	1.60	2.60a	4.8e−17	0.35	8.54	0.40	1.40	3.40
GRB 071018	08:37:41	164.7	53.8	−50.0	+417.7	2.00a	1.1e+01	0.20a	1.63	2.63a	1.4e−16	0.35	8.54	0.37	1.37	3.37
GRB 071020	07:02:27	119.7	32.9	−3.0	+7.4	2.15	2.6e+01	0.32	0.65	2.65	5.3e−15	0.20	8.13	0.35	2.35	4.35
GRB 071021	09:41:33	340.6	23.7	−31.4	+252.2	2.00a	1.3e+01	0.20a	1.70	2.70a	1.4e−16	0.35	8.54	0.30	1.30	3.30
GRB 071025	04:08:54	355.1	31.8	+38.5	+193.8	2.00a	5.9e+01	0.20a	1.79	2.79a	5.4e−16	0.35	8.54	0.21	1.21	3.21
GRB 071028A	17:41:01	119.8	21.5	+0.0	+48.9	2.00a	2.4e+00	0.20a	1.87	2.87a	2.0e−17	0.35	8.54	0.13	1.13	3.13
GRB 071101	17:53:46	48.2	62.5	−1.9	+10.0	2.00a	3.7e−01	0.20a	2.25	3.25a	2.1e−18	0.35	8.54	−0.25	0.75	2.75
GRB 071104	11:41:23	295.6	14.6	−5.0	+17.0	2.00a	3.2e+01a	0.20a	1.00a	2.00a	1.8e−15	0.35	8.54	1.00	2.00	4.00
GRB 071109	20:36:05	289.9	2.0	−5.0	+35.0	2.00a	3.2e+01a	0.20a	1.00a	2.00a	1.8e−15	0.35	8.54	1.00	2.00	4.00
GRB 071112C	18:32:57	39.2	28.4	−5.0	+30.0	0.82	6.3e+01	0.20a	1.09	2.09a	4.3e−16	0.95	14.07	0.91	1.91	3.91
GRB 071118	08:57:17	299.7	70.1	−25.0	+110.0	2.00a	5.6e+00	0.20a	1.63	2.63a	6.8e−17	0.35	8.54	0.37	1.37	3.37
GRB 071122	01:23:25	276.6	47.1	−29.4	+47.3	1.14	5.4e+00	0.20a	1.77	2.77a	1.7e−17	0.69	11.97	0.23	1.23	3.23
GRB 071125	13:56:42	251.2	4.5	−0.5	+8.5	2.00a	3.4e+02	0.30	0.62	3.10	3.8e−14	0.24	8.54	−0.10	2.38	4.38
GRB 080121	21:29:55	137.2	41.8	−0.4	+0.4	2.00a	7.9e−02	0.20a	2.60	3.60a	3.6e−19	0.35	8.54	−0.60	0.40	2.40
GRB 080205	07:55:51	98.3	62.8	−10.1	+105.3	2.00a	1.3e+01	0.20a	2.08	3.08a	8.0e−17	0.35	8.54	−0.08	0.92	2.92
GRB 080211	07:23:39	44.0	60.0	−10.0	+50.0	2.00a	1.3e+02	0.35	0.61	2.62	3.1e−14	0.20	8.54	0.38	2.39	4.39
GRB 080218A	20:08:43	355.9	12.2	−12.8	+18.6	2.00a	2.6e+00	0.20a	2.34	3.34a	1.3e−17	0.35	8.54	−0.34	0.66	2.66
GRB 080307	11:23:30	136.6	35.1	+1.7	+146.1	2.00a	8.0e+00	0.20a	1.78	2.78a	7.4e−17	0.35	8.54	0.22	1.22	3.22
GRB 080310	08:37:58	220.1	−0.2	−71.8	+318.7	2.43	9.6e+00	0.20a	2.32	3.32a	7.2e−17	0.27	7.47	−0.32	0.68	2.68
GRB 080315	02:25:01	155.1	41.7	−5.0	+65.0	2.00a	4.3e−01	0.20a	2.51	3.51a	2.0e−18	0.35	8.54	−0.51	0.49	2.49
GRB 080319C	12:25:56	259.0	55.4	−0.3	+51.2	1.95	1.5e+02	0.11	1.01	1.87	1.8e−15	0.68	8.68	1.13	1.99	3.99
GRB 080319D	17:05:09	99.5	23.9	+0.0	+50.0	2.00a	2.4e+00	0.20a	1.92	2.92a	1.8e−17	0.35	8.54	0.08	1.08	3.08
GRB 080320	04:37:38	177.7	57.2	−60.0	+40.0	2.00a	2.8e+00	0.20a	1.70	2.70a	2.9e−17	0.35	8.54	0.30	1.30	3.30
GRB 080325	04:09:17	277.9	36.5	−29.3	+170.5	2.00a	5.2e+01	0.20a	1.68	2.68a	5.7e−16	0.35	8.54	0.32	1.32	3.32
GRB 080328	08:03:04	80.5	47.5	−2.2	+117.5	2.00a	1.0e+02	0.28	1.13	3.13	5.4e−15	0.25	8.54	−0.13	1.87	3.87
GRB 080330	03:41:16	169.3	30.6	−0.5	+71.9	1.51	1.0e+00	0.20a	2.53	3.53a	3.1e−18	0.50	10.20	−0.53	0.47	2.47

Notes. Columns: T₀—trigger time of satellite (UT), R.A.—right ascention of GRB (°), decl.—declination of GRB (°), T₁ − T₀—start of prompt window (s), T₂ − T₀—end of prompt window (s), f_γ (MeV⁻¹ cm⁻²), _γ (MeV), f_ν (GeV⁻¹ cm⁻²), _ν,1 (PeV), _ν,2 (PeV). The parameters f_γ and f_ν are the fluxes at _γ and _ν,1 of the gamma-ray and neutrino spectrum, respectively (see also the Appendix). ^aParameter has not been measured. Instead, an average value is used (see Table 2).

Download table as:  ASCIITypeset image

Table 2. Average Values of GRB Parameters Taken from Becker (2008)

Parameter		Average Value
fγ		1.3 MeV−1 cm−2 a
z		2
γ		0.2 MeV
αγ		1
βγ		αγ + 1
Lisoγb		1051 erg cm−2
Γjetb		300
tvarb		0.01 s
eb		0.1
Bb		0.1
feb		0.1

Notes. For parameter definitions see the Appendix. ^aCorresponds to ${\cal F_\gamma } = 10^{-5}\,\mathrm{erg}\,\mathrm{cm}^{-2}$ between 10 keV and 10 MeV. ^bNot measured for any GRB.

Download table as:  ASCIITypeset image

2.2. Precursor Emission

2.3. Wide Window Emission

IceCube (Achterberg et al. 2006), the successor of the AMANDA experiment and the first next-generation neutrino telescope, is currently being installed in the deep ice at the geographic South Pole. Its final configuration will instrument a volume of about 1 km³ of clear ice in depths between 1450 m and 2450 m. Neutrinos are reconstructed by detecting the Cherenkov light from charged secondary particles, which are produced in interactions of the neutrinos with the nuclei in the ice or the bedrock below. The optical sensors, known as Digital Optical Modules (DOMs), consist of a 25 cm Hamamatsu photomultiplier tube (PMT) housed in a pressure-resistant glass sphere and associated electronics (Abbasi et al. 2009a). They are mounted on vertical strings where each string carries 60 DOMs. The final detector will contain 86 such strings spaced horizontally at approximately 125 m intervals.⁴⁵ Physics data taking with IceCube started in 2006 with nine strings installed. The completion of the detector construction is planned for the year 2011. The analyses described here use data taken with the 22 string configuration of the detector, which operated between 2007 May 31 and 2008 April 5.

The data acquisition (DAQ) system of IceCube (Abbasi et al. 2009a) is based on local coincidences of photon signals (hits, threshold 0.25 photoelectrons) in neighboring or next-to-neighboring DOMs on a string within 1 μs. All data from DOMs belonging to a local coincidence are read out and the digitized waveforms are sent to a computer farm at the surface. In order to pass the trigger a minimum number of eight DOMs in local coincidences within a time window of 5 μs is required. If this condition is fulfilled the waveforms are combined to an event and the number and arrival times of the Cherenkov photons are extracted. Here, the relative timing resolution of photons within an event is about 2 ns. The absolute time of an event is determined by a GPS clock to a precision of ${\cal O}({\mu }\mathrm{s})$, which is more than sufficient for our analyses.

Data are transferred from the South Pole to a computer center in the North via satellite. For the analyses described in this paper, we consider only muons produced in charged current interactions

$$\begin{equation} \nu _{\mu } + N \rightarrow \mu + X \end{equation} \tag{ 1 }$$

as only the track-like hit pattern of muons allows for a good angular resolution. In order to fit satellite bandwidth restrictions, a filter removes events which do not qualify as good upgoing neutrino candidates.⁴⁶ For this purpose, an initial track is reconstructed for each event using the line-fit algorithm (Ahrens et al. 2004). This is a simple but fast analytic track reconstruction based on the measured hit times in the DOMs. The transferred data form the basis for our analyses.

For our analyses, we use data taken with the IceCube detector in its 22 string configuration from 2007 June 1 to 2008 April 4. To prevent bias in our analyses, the data within the −1 hr to +3 hr windows (on-time data) are initially only checked for detector stability until all parameters of the analysis have been fixed, i.e., they are not used in the optimization of the analyses. The remaining off-time data that pass basic quality criteria (95% of all the data collected in the 22 string configuration) amount to 268.9 days of live-time. This long off-time window allows for a precise experimental determination of the background rate in the on-time windows.

In addition to the basic quality criteria, the on-time data are tested for stability by fitting a Gaussian to the distribution of low level event counts in each 1 s bin. In order to check for unexpected periods of dead time, an exponential is fitted to the distribution of time delays between events. Data with deviations from the expected shapes are inspected more closely for possible causes which make the data unsuitable for our analysis. Here, most observed deviations can be attributed to transitions between runs and are therefore not critical. Out of 48 northern hemisphere bursts in the time period under investigation, the data for seven do not pass the stability/quality criteria or have gaps during the prompt/precursor emission windows. For all remaining 41 GRBs, both tests show excellent agreement with no indications of abnormal behavior of the detector during the on-time periods. The on-time data cover 100% of the prompt and precursor windows for the 41 selected bursts. For the extended window seven out of the 41 bursts exhibit gaps at the beginning and/or end of the window. For these seven bursts, the extended window is shortened accordingly (see Table 3). With this correction, the data cover 94% of the extended time windows of all GRBs. The missing 6% are due to larger gaps in data taking.

After the data have been transferred to the North, a more precise determination of the direction of an event is achieved by fitting a muon-track hypothesis to the hit pattern of the recorded Cherenkov light in the detector using a log-likelihood reconstruction method (Ahrens et al. 2004). A fit of a paraboloid to the region around the minimum in the log-likelihood function yields an estimate of the uncertainty on the reconstructed direction (Neunhöffer 2006). At this point, the data sample with an event rate of 3.3 Hz is still dominated by several orders of magnitude by misreconstructed downgoing atmospheric muons as demonstrated in Figure 2, which shows a comparison between data and Monte Carlo (see Section 5). The quantities shown are later used to reject these misreconstructed atmospheric muons and improve the sensitivity of the analyses.

• 1.  
  θ_rec. Reconstructed zenith angle.⁴⁷
• 2.  
  σ_dir. The uncertainty on the reconstructed track direction (quadratic average of the minor and major axes of the 1σ error ellipse).
• 3.  
  L_red. log  of the likelihood value of the reconstructed track divided by the number of degrees of freedom (number of hit DOMs minus number of fit parameters). This has proven to be a powerful variable for separating signal and background as visible from Figure 2.
• 4.  
  L_U/D. Difference in log-likelihood value between the reconstructed track and one containing a bias to be reconstructed as downgoing. The bias is zenith angle dependent and follows the rate of downgoing atmospheric muons. The rationale behind this is that a track is much more likely to originate from an atmospheric muon than from a muon generated in a neutrino interaction. Only high-quality upgoing tracks have high L_U/D values.
• 5.  
  N_dir. The number of photons detected within a −15 to +75 ns time window with respect to the expected arrival time for unscattered photons from the muon-track hypothesis.
• 6.  
  θ_min. Minimum zenith angle from a fit of a two-track hypothesis to the light pattern. For this, the light pattern is divided into two separate sets of hits based on the mean hit time of the event, and a track hypothesis is fitted to each hit set separately. A cut on the smaller zenith angle, θ_min, of the two tracks is very effective against the so-called coincident muons, where two muons from different atmospheric showers pass the detector in fast succession mimicking an upgoing track. In this case, often at least one of the two tracks is reconstructed as downgoing whereas for good-quality upgoing neutrinos both tracks appear most often as upgoing.
• 7.  
  E_rec. Reconstructed muon energy at point of closest approach to the center of gravity of hits in an event (Zornoza et al. 2008). This is the calibrated output of a likelihood reconstruction method evaluating the measured number of photons in each DOM with respect to the corresponding probability density function (PDF) of a given track-energy hypothesis.Muons carry a significant fraction of the original neutrino energy, and at the energies of highest acceptance in our analyses (∼100 TeV), have a range of about 10 km. The dominant energy-loss mechanisms for these muons are Bremsstrahlung and pair production which grow with increasing energy, thereby increasing the amount of Cherenkov light emitted and allowing one to estimate the muon energy.

Zoom In Zoom Out Reset image size

Signal neutrinos are generated from the direction of the GRB with the corresponding spectrum using a port of the anis code (Gazizov & Kowalski 2005) called neutrino-generator. Here, the change in position of the source in the detector coordinate system during the respective time window is taken into account. The neutrino emission is assumed to be constant during this time. In addition, three types of background are simulated. At the beginning of the analysis chain, the data sample is dominated by downgoing atmospheric muons, produced by interactions of cosmic rays in the atmosphere, which are reconstructed as upgoing. As more cuts are applied the data sample starts to be dominated by coincident muon events. Both event classes are simulated with the corsika air shower simulation package (Heck et al. 1998). Finally, we consider the irreducible background of atmospheric neutrinos generated in the same interactions as the atmospheric muons. These neutrinos are simulated as an all-sky flux with neutrino-generator and weighted according to the Bartol spectrum (Barr et al. 2004).

Neutrinos are tracked from the surface through the Earth taking into account absorption, scattering, and neutral current regeneration of neutrinos (Gazizov & Kowalski 2005). Information on the structure of the Earth is taken from the Preliminary Reference Earth Model (Dziewonski & Anderson 1981). Muons originating from neutrino interactions near the detector and atmospheric muons are traced through rock and ice taking into account continuous and stochastic energy losses (Chirkin & Rhode 2004). The photon signal in the DOMs is determined from a detailed simulation of the propagation of Cherenkov light from muons and showers through the ice (Lundberg et al. 2007) which includes the modeling of the changes in absorption and scattering length with depth due to dust layers (Ackermann et al. 2006). This is followed by a simulation of the DOM electronics and the trigger. The simulated DOM signals are then processed in the same way as the data.

Muon neutrinos from GRBs show up as an excess of tracks above the background from the direction of the GRB within a certain time window. Background from misreconstructed atmospheric muons can be suppressed by applying quality cuts on reconstructed quantities. Background from atmospheric neutrinos on the other hand is indistinguishable from cosmic neutrinos. Hence, once a high-purity (atmospheric) neutrino sample has been selected, quality cuts which aim at rejecting misreconstructed tracks cannot further improve the signal to background ratio. However, as neutrinos from GRBs are expected to exhibit a harder energy spectrum than that of atmospheric neutrinos, information on the muon energy allows to further increase the sensitivity of the analysis to neutrinos from GRBs.

We have analyzed the data both with an unbinned likelihood and a binned method. In order to enhance the chances for a discovery, we do not analyze the 41 GRBs individually but as a population (stacked analysis).⁴⁸ This allows us to set the most stringent limit on the tested models but might not be optimal in other cases, e.g., if one burst has a much higher neutrino fluence relative to the rest of the GRBs than expected. Also, in case of a discovery a stacked analysis only allows one to calculate the neutrino fluence from the whole burst population but not from individual GRBs. In the following, we describe the unbinned likelihood and binned methods and compare their performances in the case of the prompt emission scenario. The more sensitive method is then used to obtain the results presented in Section 7.

6.1. Unbinned log-likelihood Method

After filtering at the South Pole and transfer to the north, the data sample is still dominated by downgoing muons which are reconstructed as upgoing. These muons are rejected by applying quality cuts on reconstructed quantities (for a description see Section 4).

$$\begin{eqnarray} && \theta _\mathrm{rec} > 85^\circ ; \sigma _\mathrm{dir} < 3^\circ ; \theta _\mathrm{min} > 70^\circ ; L_{U/D} > 30;\nonumber\\ && L_\mathrm{red} \le \left\lbrace \begin{array}{@{}l@{\quad }l} 7.8 \quad {\rm for }\quad \ N_\mathrm{dir} < 7\\ 8.5 \quad {\rm for }\quad \ N_\mathrm{dir} = 7\\ 9.5 \quad {\rm for }\quad \ N_\mathrm{dir} > 7 \end{array}\right.. \end{eqnarray} \tag{ 2 }$$

After these cuts, a high-purity upgoing (atmospheric) neutrino sample remains with an event rate of 2.1 × 10⁻⁴ Hz. Table 4 lists the cut efficiency for data and signal. A comparison of data with simulations is displayed in Figure 3. In general, good agreement between the atmospheric neutrino Monte Carlo and data is observed. Small deviations for example visible at low N_dir are most likely due to atmospheric muon background which is not accounted for by the Monte Carlo due to its limited statistics at this cut level. The cumulative point-spread function is shown in Figure 4 (left). The median angular resolution is about 15. Figure 4 (right) shows the muon neutrino effective area for different declination bands.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Number of Signal and Background (Off-time Data) Events for the Unbinned Method at Different Cut Levels

Cut Level	Prompt Window		Precursor Window		Wide Window
	No. Events	a (%)	No. Events	a (%)	No. Eventsb	a (%)
Signal
Filter	0.062	100	1.8	100	7.6	100
Final	0.033	53	0.53	29	2.8	37
Background
Filterc	7.6 × 108	100	7.6 × 108	100	7.6 × 108	100
Finalc	4846	6.5 × 10−4	4846	6.5 × 10−4	4846	6.5 × 10−4
Final (window)d	6.1 × 10−4	8.2 × 10−11	6.1 × 10−4	8.2 × 10−11	8.8 × 10−2	1.2 × 10−8

Notes. ^aEfficiency relative to filter level. ^bFor a fluence equal to the computed upper limit in Section 7. ^cAll events after cuts. ^dExpected events in cones with radii of 23 (contain 70% of signal events; see Figure 4) around GRBs within respective time window (prompt: T₂ − T₁ from Table 1; precursor: 100 s; wide: 4 hr). Note, that these cones are not used in the evaluation of the data with the unbinned likelihood method.

Download table as:  ASCIITypeset image

The data sets after quality cuts are the starting point for the unbinned likelihood method. In contrast to binned methods where the event is rejected if it lies outside the cut region (binary selection), unbinned likelihood methods do not discard events but use PDFs to evaluate the probability of an event belonging to signal or background population. The unbinned likelihood method used here is similar to that described in Braun et al. (2008). The signal, $S(\vec{x}_i)$, and background, $B(\vec{x}_i)$, PDFs are each the product of a time PDF, a directional PDF, and an energy PDF, where $\vec{x}_i$ denotes the directional, time, and energy variables.

The directional signal PDF is a two-dimensional Gaussian distribution with the two widths being the major and minor axes of the 1σ error ellipse of the paraboloid fit described in Section 4. The time PDF is flat over the respective time window and falls off on both sides with a Gaussian distribution. The width σ of the Gaussian is determined by the length of the time window with a maximum of σ = 25 s and a minimum of σ = 2 s. The Gaussian accounts for possible small shifts in the neutrino emission time with respect to that of the γ-rays and prevents discontinuities in the likelihood function. The sensitivity of the method depends only weakly on the exact choice of σ. The energy PDF is determined for each GRB individually. It is derived from the energy-estimator distribution of the tracks of the corresponding signal Monte Carlo data set (weighted to an E⁻² spectrum⁴⁹) after final cuts (see Equation (2)). The signal PDFs of the GRBs are combined using a weighted sum (Abbasi et al. 2006)

$$\begin{equation} S_\mathrm{tot}(\vec{x}_i) = \frac{\sum _{j=1}^{N_\mathrm{GRBs}} w_j\,S_j(\vec{x}_i)}{\sum _{j=1}^{N_\mathrm{GRBs}} w_j}, \, \end{equation} \tag{ 3 }$$

where $S_j(\vec{x}_i)$ is the signal PDF of the jth GRB and w_j is a weight that in the case of the prompt and precursor window is proportional to the expected number of events in the detector according to the fluences described in Section 2. In the case of the extended window, we use w_j = 1 for all GRBs in order to make the search as general as possible.

For the directional background PDF, the detector asymmetries in zenith and azimuth must be taken into account. This is accomplished by evaluating the data in the detector coordinate system. The directional background PDF is hence derived from the distribution of all off-time events after final event selection in the zenith–azimuth plane of the detector. The time distribution of the background during a GRB can be assumed to be constant resulting in a flat time PDF. The energy PDF is determined in the same way as for the signal PDF with weights corresponding to the Bartol atmospheric neutrino flux.

All PDFs are combined in an extended log-likelihood function (Barlow 1989)

$$\begin{equation} \ln \left({\cal L}(\langle n_s \rangle)\right) = -\langle n_s \rangle - \langle n_b \rangle + \sum _{i=1}^{N} \ln \left(\langle n_s \rangle \,S_\mathrm{tot}(\vec{x}_i) + \langle n_b \rangle \,B(\vec{x}_i) \right), \, \end{equation} \tag{ 4 }$$

where the sum runs over all reconstructed tracks in the final sample. The variable 〈n_b〉 is the expected mean number of background events, which is determined from the off-time data set. The mean number of signal events, 〈n_s〉, is a free parameter which is varied to maximize the expression

$$\begin{eqnarray} \ln \left({\cal R}(\langle n_s \rangle)\right) &=& \ln \left(\frac{{\cal L}(\langle n_s \rangle)}{{\cal L}(0)}\right)\nonumber\\ &=& {-}\langle n_s \rangle + \sum _{i=1}^N \ln \left(\frac{\langle n_s \rangle \,S_\mathrm{tot}(\vec{x}_i)}{\langle n_b \rangle \,B(\vec{x}_i)} + 1 \right) \end{eqnarray} \tag{ 5 }$$

in order to obtain the best estimate for the mean number of signal events, $\widehat{\langle n_s \rangle }$.

Zoom In Zoom Out Reset image size

To determine whether a given data set is compatible with the background-only hypothesis 10⁸ background data sets for the on-time windows are generated from off-time data by randomizing the track times while taking into account the downtime of the detector. For each of these data sets the $\ln ({\cal R})$ value is calculated, yielding the distribution shown in Figure 5. The probability for a data set to be compatible with background is given by the fraction of background data sets with a larger $\ln ({\cal R})$ value. For comparison, the plot also displays the $\ln ({\cal R})$ distributions for background data sets with one and two injected Monte Carlo signal events, respectively. The signal events are randomly distributed among the GRBs, where the assignment probability to a specific GRB is proportional to the expected number of events from that burst. The energy of the neutrinos is generated according to the spectra calculated in Section 2.1.

6.2. Binned Method

For the binned method, a machine learning algorithm was trained to separate signal and background. The algorithm used was a Support Vector Machine (SVM; Cortes & Vapnik 1995) with a radial basis function kernel. It was provided with the best reconstructed track direction in detector coordinates as well as many quality parameters including σ_dir, L_red, L_U/D, N_dir, and θ_min which are described in Section 4. The SVM was trained using the off-time filtered data as background and all-sky neutrino simulation weighted to the sum of the individual burst spectra as signal. The optimum SVM parameters (kernel parameter, cost factor, margin) were determined using a coarse, and then fine, grid search with a fivefold cross-validation technique at each node, as described in Hsu et al. (2003).

The resulting SVM classification of events is shown in Figure 6. The final cut on this classifier is optimized to detect a signal fluence with at least 5σ (significance) in 50% of cases (power) by minimizing the Model Discovery Factor (MDF) according to Hill et al. (2006). The MDF is the ratio between the signal fluence required for a detection with the specified significance and power and the predicted fluence. The angular cut around each GRB is then calculated to keep 3/4 of the remaining signal after the cut on the SVM classifier. In this way, there is one cut on the SVM classifier for all GRBs, but different angular cuts around each GRB according to the angular resolution of the detector in that direction. The optimum SVM cut is determined to be at a value of 0.25. Table 5 displays the signal and background event rates at different cut levels.

Zoom In Zoom Out Reset image size

Table 5. Number of Signal and Background (Off-time Data) Events for the Binned Method in the Prompt Window at Different Cut Levels

Cut Level	Signal		Background
	No. Events	Efficiencya (%)	No. Events	Efficiencya (%)
Filter	0.062	100	77 × 107	100
Final	0.023	37	4.7	6.1 × 10−9

Note. ^aRelative to filter level.

Download table as:  ASCIITypeset image

6.3. Comparison of the Two Methods

We compare the performance of the unbinned likelihood and binned method by means of their discovery potential for the prompt neutrino emission scenario. Figure 7 displays for a significance of 5σ the power as a function of the MDF for the two methods. For a power of 50% the unbinned method shows an improvement in the MDF of about a factor 1.8 compared to the binned method. It is therefore used to derive the results presented in the following section. The large gain in sensitivity is partly due to the explicit use of energy information in the likelihood. The discovery potential for the expected fluxes is further improved by weighting the bursts according to the expected number of signal events in the detector in the unbinned method (Equation (3)). In the binned method, all bursts are treated equally.

Zoom In Zoom Out Reset image size

We apply the unbinned likelihood method to the on-time data sets after neutrino candidate event selection with the final cuts (Equation (2)). For all three emission scenarios, the values of $\ln ({\cal R})$ and $\widehat{\langle n_s \rangle }$ are zero and hence consistent with the null hypothesis. Therefore, we derive 90% CL upper limits⁵⁰ on the fluence from the 41 GRBs in the prompt phase of 3.7 × 10⁻³ erg cm⁻² (72 TeV–6.5 PeV) and on the fluence from the precursor phase of 2.3 × 10⁻³ erg cm⁻² (2.2–55 TeV), where the quoted energy ranges contain 90% of the expected signal events in the detector. Further information is listed in Table 6. The limits, which are displayed in Figure 8, are not strong enough to constrain the models. The 90% CL upper limit for the wide time window is 2.7 × 10⁻³ erg cm⁻² (3 TeV–2.8 PeV) assuming an E⁻² flux.

Zoom In Zoom Out Reset image size

To illustrate, we counted the number of events after final cuts in cones with radii 23 around the GRB positions (contain 70% of signal events; see Figure 4) within the corresponding time windows (prompt: T₂ − T₁ from Table 1; precursor: 100 s; wide: 4 hr). In addition, we analyzed the data in the prompt window with the binned method after final cuts. In all cases, zero events remain which is consistent with the results of the unbinned likelihood method.

As described previously, we use the off-time data to determine the background rate in the on-time windows. This technique removes many potential sources of uncertainties in the calculation of the significance of a possible signal that are introduced when using a simulation of the background. However, this method makes the assumption that the rates of data during the off-time and on-time windows are the same. Furthermore, we use Monte Carlo for the signal simulation and the derivation of upper limits, which involves the propagation of particles through the Earth and ice, and the simulation of the detector response. The most important sources for systematic uncertainties are discussed below in detail. Their effects on the upper limits are summarized in Table 7.

Ice simulation. Inaccuracies in the ice simulation can lead to a wrong estimate of the efficiency of the detector to neutrinos from GRBs. Data-Monte Carlo comparisons and variation of simulation parameters indicate that the systematic uncertainty from this aspect of the simulation is about ±15%.

DOM efficiency. A ±10% uncertainty in the efficiency of the optical modules in the detection of photons leads to a corresponding uncertainty in the number of expected events from a GRB. This effect is nonlinear and spectrally dependent, so simulation was generated spanning the range of uncertainties to determine the resulting change in signal event rates.

Neutrino and muon propagation. Theoretical uncertainties on muon energy losses and the neutrino–nucleon cross section, determined from the uncertainty on the CTEQ6 PDFs (Pumplin et al. 2002), contribute a 5% uncertainty on the neutrino event rate in the detector.

Background rate. After final cuts, the variation of the event rate over the data-taking period is about ±5%. In order to account for potential differences at the time of the bursts, the background data rate is varied by this amount. This results in a shift of the upper limits of less than ±1% and is therefore negligible.

A search with the AMANDA detector for muon neutrinos in the prompt and precursor phase was conducted for GRBs detected between 1997 and 2003 (Achterberg et al. 2008) with null result. The analysis of prompt neutrinos contained 419 bursts observed by the BATSE experiment (Paciesas et al. 1999) as well as by others with similar characteristics. The large number of bursts allowed that analysis to set an upper limit only a factor 1.4 above the prediction of the Waxman–Bahcall prompt emission model. A further search with the AMANDA detector for muon neutrinos was conducted for 85 GRBs detected between 2005 and 2006 (Strahler 2009), primarily by the Swift satellite, also with null result. Due to the smaller number of bursts, the upper limit from this analysis is much less restrictive. Converting the upper limit obtained from the 419 bursts to a fluence limit from 41 standard Waxman–Bahcall bursts (for definition see footnote 41) yields the dotted line in Figure 8. Due to the 10 times smaller number of bursts available, the limit presented in this paper is about a factor 3 worse. For the precursor emission model, the AMANDA limit (Achterberg et al. 2008; Kuehn 2007) is much less restrictive as only 60 bursts, detected between 2001 and 2003, were used. It is shown in Figure 8 as a dash-dotted line. Here, our analysis improves on the AMANDA upper limit despite the fact that 30% fewer bursts were investigated.

The AMANDA data were also analyzed for neutrinos of all flavors from GRBs (Achterberg et al. 2007). Apart from a search for neutrinos from 73 bursts detected by BATSE in 2000, another search did not rely on information from satellites but looked for a clustering of events within sliding time windows of 1 s and 100 s. Due to this more generic approach and the low number of bursts, respectively, the limits from these analyses are much less restrictive than those from Achterberg et al. (2008). However, the sliding window or similar searches are the only way to detect GRBs where either the jet does not emerge from the progenitor star (choked bursts; Meszaros & Waxman 2001) or the γ-ray signal is not observed by satellites for other reasons.

We continue with a discussion of the impact of results from high-energy γ-ray observations on expected neutrino fluences. Along with high-energy neutrinos which originate from the decay of charged pions, high-energy γ-ray photons are produced in the decay of simultaneously generated neutral pions. In addition, high-energy photons are produced in inverse Compton scattering of synchrotron photons by accelerated electrons (Falcone et al. 2008). In contrast to neutrinos, the flux of high-energy γ-ray photons at the Earth is significantly reduced due to the large optical depths for photon–photon pair production inside the source for not too large jet Lorentz factors Γ_jet ≲ 800 (Falcone et al. 2008). In addition, high-energy photons above 100 GeV are absorbed on the extragalactic background light (EBL) if they travel distances with z ≳ 0.5. Observations with air-Cherenkov telescopes like H.E.S.S. (Aharonian et al. 2009; Tam et al. 2009) or MAGIC (Albert et al. 2007; Galante et al. 2009) are also hampered by the fact that usually it takes more than 50 s (MAGIC) or 100 s (H.E.S.S.) from the observation of a GRB by a satellite to the start of data taking with these telescopes. Therefore, the prompt emission window is only partially covered or not at all. MILAGRO as an air shower array observed large parts of the sky continuously (Atkins et al. 2004, 2005; Abdo et al. 2007). However, it was mostly sensitive to energies above 100 GeV and therefore suffered significantly from γ-ray absorption on the EBL. HAWC, the successor of Milagro currently in the planning phase, will be able to detect γ rays from GRBs down to 100 GeV where about 50 (Fermi) GRBs per year will fall in its field of view (Goodman et al. 2009).

At energies above 100 GeV, there has been no definitive detection of γ-ray emission from GRBs. Milagrito (Atkins et al. 2000, 2003) and the HEGRA AIROBICC array (Padilla et al. 1998) reported evidence at the 3σ level for high-energy γ-ray emission from GRB 970417A (E_γ>650 GeV) and GRB 920925C (E_γ>20 TeV), respectively. However, subsequent searches for high-energy γ-ray emission from GRBs did not find similar signals. The limits obtained from MAGIC and H.E.S.S. are not directly comparable to our results due to their incomplete burst coverage. In Atkins et al. (2004), Atkins et al. (2005), and Abdo et al. (2007), the MILAGRO collaboration reports 99% CL upper limits for a large number of individual bursts (both long and short) between 2000 and 2006 down to 10⁻⁷ erg cm⁻² (energy range ∼100 GeV–10 TeV; the exact energy range differs from publication to publication). Extrapolating the average per burst 99% CL upper limit for the prompt window in our analysis to the energy range from 100 GeV to 10 TeV yields 3.2 × 10⁻⁷ erg cm⁻². However, the photon limits do not constrain our results as they do not account for absorption in the EBL (this would significantly worsen the limits) and include an unknown, probably dominant, contribution from inverse Compton scattering. In general, current flux predictions for high-energy gamma rays from GRBs are near or below the sensitivity of current instruments (Falcone et al. 2008), where the predicted fluxes in the energy range below ∼100 TeV are dominated by the leptonic emission component in most scenarios.

Within the internal shock (fireball) model, synchrotron self-Compton (SSC) processes between the accelerated electrons and the γ-ray photons could lead to a high-energy γ-ray peak in the GeV range detectable by Fermi in case of bright GRBs. Up to now, this has happened only for a handful of bursts (one of the photons with the highest energy from a GRB detected by Fermi came from GRB 080916C and had an energy of ∼13 GeV; Abdo et al. 2009). However, this does not necessarily disfavor the internal shock scenario. As discussed in Fan (2009), the amount of energy in high-energy photons could be suppressed by an inefficient SSC process in the extreme Klein–Nishina regime or a combination of a SSC peak at high energies and a low photon cutoff energy above which the fireball becomes optical thick.

We have performed a set of complementary searches for muon neutrinos associated in space and time with 41 GRBs that were observed in the northern sky between 2007 June and 2008 April. For the first time in searches with large GRB populations, we have calculated individual prompt neutrino spectra for all 41 GRBs using measured GRB parameters. The search results are consistent with the case of a background-only hypothesis. Therefore, we place 90% CL upper limits on the fluence from the prompt phase of 3.7 × 10⁻³ erg cm⁻² (72 TeV–6.5 PeV) and on the fluence from the precursor phase of 2.3 × 10⁻³ erg cm⁻² (2.2–55 TeV), where the quoted energy ranges contain 90% of the expected signal events in the detector. Though the number of bursts is smaller than in previous searches the larger detector allows us to improve on the limits for the precursor phase by a factor 1.4. Compared to the predictions, the limits lie a factor 72 (prompt phase) and 9.7 (precursor phase) higher. Hence, they do not allow us to constrain the models. Apart from these model-driven searches, we have also conducted for the first time a generic search for neutrino emission from GRBs in a wide window of (−1 hr to +3 hr) around each burst. Finding no evidence for a signal, we place a 90% CL upper limit on the fluence of 2.7 × 10⁻³ erg cm⁻² (3 TeV–2.8 PeV) assuming an E⁻² flux.

Launched in 2008 June, the Fermi Gamma-ray Space Telescope⁵¹ has begun to provide an expanded catalog of sources for future neutrino searches. With a much larger field of view than other satellites, Fermi has increased the GRB detection rate by more than a factor 2. At the same time, the detected bursts have on average a higher luminosity than those detected by Swift due to the lower sensitivity of the Fermi-GBM (Gamma-ray Burst Monitor) instrument. The first IceCube analysis to take advantage of the increased detection opportunities will utilize the 40 string configuration of the detector, already as large as the full IceCube along one axis. This gives it the full angular resolution power along that direction and thus provides powerful background rejection. However, this is mitigated by the comparatively poor angular resolution of the Fermi-GBM (∼3°), which is worse than the IceCube resolution. With the full 80 string detector scheduled to be completed in 2011 and an expected 100–150 detected bursts per year in the northern hemisphere by the Fermi and Swift satellites the sensitivity of IceCube to neutrinos from GRBs will soon exceed that of AMANDA. This will allow IceCube to either confirm the predicted fluxes within the next years or set stringent limits thereby disfavoring GRBs as the major sources of ultra-high-energy cosmic rays.

We acknowledge the support from the following agencies: U.S. National Science Foundation-Office of Polar Program, U.S. National Science Foundation-Physics Division, University of Wisconsin Alumni Research Foundation, U.S. Department of Energy, and National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; Swedish Research Council, Swedish Polar Research Secretariat, and Knut and Alice Wallenberg Foundation, Sweden; German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), Germany; Fund for Scientific Research (FNRS-FWO), Flanders Institute to encourage scientific and technological research in industry (IWT), Belgian Federal Science Policy Office (Belspo); the Netherlands Organisation for Scientific Research (NWO); Marsden Fund, New Zealand; M. Ribordy acknowledges the support of the SNF (Switzerland); A. Kappes and A. Groß acknowledge support by the EU Marie Curie OIF Program; J. P. Rodrigues acknowledge support by the Capes Foundation, Ministry of Education of Brazil.

$$\begin{eqnarray} F_\gamma (E_\gamma) &=& \frac{{d}N(E_\gamma)}{{d}E_\gamma } \nonumber\\ &=& f_\gamma \times \left\{\begin{array}{@{}cc}\left(\frac{\epsilon _\gamma }{\mathrm{MeV}}\right)^{\alpha _\gamma } \ \left(\frac{E_\gamma }{\mathrm{MeV}}\right)^{-\alpha _\gamma } & {\rm for}\ E_\gamma < \epsilon _\gamma \\ \left(\frac{\epsilon _\gamma }{\mathrm{MeV}}\right)^{\beta _\gamma } \ \left(\frac{E_\gamma }{\mathrm{MeV}}\right)^{-\beta _\gamma } & {\rm for}\ E_\gamma \ge \epsilon _\gamma \end{array}\right.\!\!,\quad \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} &&\fl {\cal F}_\gamma = \int {d}E_\gamma \ E_\gamma F_\gamma (E_\gamma), \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} &&\hspace{-6pt}F_\nu (E_\nu) = \frac{{d}N(E_\nu)}{{d}E_\nu } = f_\nu\hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaa}\nonumber\\ &&\hspace{-6pt}\times \left\{ \begin{array}{@{}ll@{}} \left(\frac{\epsilon _{\nu,1}}{\mathrm{GeV}}\right)^{\alpha _\nu } \ \left(\frac{E_\nu }{\mathrm{GeV}}\right)^{-\alpha _\nu } & {\rm for}\ E_\nu < \epsilon _{\nu,1} \\ \left(\frac{\epsilon _{\nu,1}}{\mathrm{GeV}}\right)^{\beta _\nu } \ \left(\frac{E_\nu }{\mathrm{GeV}}\right)^{-\beta _\nu } & {\rm for}\ \epsilon _{\nu,1} \le E_\nu < \epsilon _{\nu,2} \\ \left(\frac{\epsilon _{\nu,1}}{\mathrm{GeV}}\right)^{\beta _\nu } \ \left(\frac{\epsilon _{\nu,2}}{\mathrm{GeV}}\right)^{\gamma _\nu -\beta _\nu } \ \left(\frac{E_\nu }{\mathrm{GeV}}\right)^{-\gamma _\nu } & {\rm for}\ E_\nu \ge \epsilon _{\nu,2} \end{array}\right.\!,\nonumber\\ \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} &&\fl \epsilon _1 = 7\times 10^{5}\,\mathrm{GeV} \ \frac{1}{(1+z)^2} \left(\frac{\Gamma _\mathrm{jet}}{10^{2.5}}\right)^2 \left(\frac{\mathrm{MeV}}{\epsilon _\gamma } \right), \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \fl &&\epsilon _2 = 10^7 \, \mathrm{GeV} \ \frac{1}{1+z} \ \sqrt{\frac{\epsilon _e}{\epsilon _B}} \ \left(\frac{\Gamma _\mathrm{jet}}{10^{2.5}}\right)^4 \left(\frac{t_\mathrm{var}}{0.01\,\mathrm{s}} \right) \hphantom{aaaaaaaaaa\,}\nonumber\\ \fl &&\hphantom{\epsilon _2 =}\times\sqrt{\frac{10^{52}\,\mathrm{erg\,s}^{-1}}{L_\gamma ^\mathrm{iso}}}, \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} &&\fl \alpha _\nu = 3 - \beta _\gamma, \quad \beta _\nu \ = \ 3 - \alpha _\gamma, \quad \gamma _\nu = \beta _\nu + 2, \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \fl \frac{\Delta R}{\lambda _{p\gamma }} &= \displaystyle\left(\frac{L_\gamma ^\mathrm{iso}}{10^{52}\,\mathrm{erg\,s}^{-1}} \right) \ \left(\frac{0.01\,\mathrm{s}}{t_\mathrm{var}} \right) \ \left(\frac{10^{2.5}}{\Gamma _\mathrm{jet}}\right)^4 \ \left(\frac{\mathrm{MeV}}{\epsilon _\gamma } \right),\nonumber \\ \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} \fl &&\int _{0}^\infty {d}E_\nu E_\nu F_\nu (E_\nu) = \frac{1}{8} \ \frac{1}{f_e} \ (1 - (1-\langle x_{p\rightarrow \pi } \rangle)^{\Delta R/\lambda _{p\gamma }})\hphantom{aaaaaa} \nonumber\\ &&\quad\quad \times\int _{1\,\mathrm{keV}}^{10\,\mathrm{MeV}} {d}E_\gamma \ E_\gamma F_\gamma (E_\gamma). \end{eqnarray} \tag{ A8 }$$

• 1.  
  Parameters of the γ-ray spectrum F_γ(E_γ):

  • a)  
    _γ: break energy;
  • b)  
    α_γ: spectrum index before break energy;
  • c)  
    β_γ: spectrum index after break energy;
  • d)  
    ${\cal F}_\gamma$: measured fluence in γ-rays integrated over the energy range given in the GCN circulars and reports (see footnote ⁴¹);
  • e)  
    f_γ: normalization; obtained from integral of Equation (A2).
• 2.  
  Parameters of the neutrino spectrum F_ν(E_ν):

  • a)  
    ₁: first break energy;
  • b)  
    ₂: second break energy;
  • c)  
    α_ν: spectrum index before first break energy;
  • d)  
    β_ν: spectrum index between first and second break energy;
  • e)  
    γ_ν: spectrum index after second break energy;
  • f)  
    f_ν: normalization; obtained from integral of Equation (A8).
• 3.  
  z: redshift of GRB;
• 4.  
  _e: fraction of jet energy in electrons;
• 5.  
  _B: fraction of jet energy in magnetic field;
• 6.  
  f_e: ratio between energy in electrons and protons;
• 7.  
  L^iso_γ: isotropic luminosity of the GRB;
• 8.  
  t_var: variability of the γ-ray light curve of the GRB;
• 9.  
  Γ_jet: Lorentz boost factor of the jet.

The expression $1 - (1-\langle x_{p\rightarrow \pi } \rangle)^{\Delta R/\lambda _{p\gamma }}$ in Equation (A8) estimates the overall fraction of the proton energy going into pions from the size of the shock, ΔR, and the mean free path of a proton for photomeson interactions, λ_pγ. Here, 〈x_p→π〉 = 0.2 is the average fraction of proton energy transferred to a pion in a single interaction. The expression ensures that the transferred energy fraction is ⩽1. The calculations are insensitive to the beaming effect caused by a narrow opening angle of the jet as all formulae contain the isotropic luminosity in conjunction with a 4π shell geometry, i.e., effectively use luminosity per steradian. For example, the target photon density used to calculate N_int is given by n_γ ∝ L^iso_γ/4πR², where R is the distance of the shock region from the central black hole.