AN ALL-SKY SEARCH FOR THREE FLAVORS OF NEUTRINOS FROM GAMMA-RAY BURSTS WITH THE ICECUBE NEUTRINO OBSERVATORY

Ultra high energy cosmic rays (UHECRs), defined by energy exceeding 10¹⁸ eV, have been observed for decades (Zatsepin & Kuzmin 1966; Abbasi et al. 2008; Abraham et al. 2010), but their sources remain unknown. Gamma-ray bursts (GRBs) are among the most plausible candidates for these particles' origins. If this hypothesis is true, neutrinos would also be produced in $p\gamma$ interactions at these sources. While the charged cosmic rays are deflected by galactic and intergalactic magnetic fields, neutrinos travel unimpeded through the universe and disclose their source directions. The detection of high energy neutrinos correlated with gamma-ray photons from a GRB would provide evidence of hadronic interaction in these powerful phenomena and confirm their role in UHECR production.

GRBs are the brightest electromagnetic explosions in the universe and are observed by dedicated spacecraft detectors at an average rate of about one per day (Meszaros 2006; von Kienlin et al. 2014). GRB locations are distributed isotropically at cosmic distances. Their prompt gamma-ray emissions exhibit diverse light curves and durations, lasting from milliseconds to hours. The origins of these extremely energetic phenomena remain unknown, but their spectra have been studied extensively over a wide range of energies (Sakamoto et al. 2011; Ackermann et al. 2013; Gruber et al. 2014). The Fermi detector covers seven decades of energy in gamma-rays with its two instruments and has observed photons with source-frame-corrected energies above 10 GeV early in the prompt phase of some bursts (Bissaldi et al. 2015). These measurements support efficient particle acceleration in GRBs.

The prevailing phenomenology that successfully describes GRB observations is that of a relativistically expanding fireball of electrons, photons, and protons (Piran 2004; Fox & Mészáros 2006; Meszaros 2006). The initially opaque fireball plasma expands by radiation pressure until it becomes optically thin and emits the observed gamma-rays. During this expansion, kinetic energy is dissipated through internal shock fronts that accelerate electrons and protons via the Fermi mechanism to the observed gamma-ray and UHECR energies (Waxman 1995; Vietri 1997). These high energy protons will interact with gamma-rays radiated by electrons and create neutrinos, for example, through the delta-resonance:

$$\begin{eqnarray}&&p+\gamma \to {{\rm{\Delta }}}^{+}\to {\pi }^{+}+n\to {e}^{+}+{\nu }_{e}+{\bar{\nu }}_{\mu }+{\nu }_{\mu }+n.\end{eqnarray} \tag{ 1 }$$

To date, no neutrino signal has been detected in searches for muon neutrinos from GRBs in multiple years of data from AMANDA, the partially instrumented IceCube, and the completed IceCube detector (Achterberg et al. 2007; Abbasi et al. 2010, 2011b, 2012a; Aartsen et al. 2015d), nor in four years of data by the ANTARES collaboration (Presani 2011; Adrián-Martínez et al. 2013a, 2013b). High energy ${\nu }_{\mu }$ charged-current interactions produce high energy muons that manifest as extended Cherenkov light patterns in the South Pole glacial ice, referred to as "tracks"; and Southern Hemisphere bursts were often excluded from searches for this signal in order to remove the dominant cosmic-ray-induced muon background. Adding the low-background "shower" channel from interactions other than charged-current ${\nu }_{\mu }$ gives enhanced sensitivity to such bursts, improving the sensitivity of the overall correlation analysis. Both the shower and track analyses are sensitive to the combined flux of neutrinos and anti-neutrinos and are unable to distinguish between them.

This paper is ordered as follows. In Section 2, we describe the neutrino spectra predicted by the different fireball models on which we place limits. In Section 3, we describe the IceCube detector and data acquisition system. We discuss the simulation and reconstruction of events in IceCube in Section 4. We detail the event selection techniques and likelihood analysis in Sections 5 and 6. Finally, in Section 7 we present the results of this all-sky three-flavor shower search in combination with those from the ${\nu }_{\mu }$ track searches, and conclude in Section 8.

In this search for GRB neutrinos, we examine data during the time of gamma-ray emission reported by any satellite for each burst. We do not consider possible precursor (Razzaque et al. 2003) or afterglow (Waxman & Bahcall 2000; Murase & Nagataki 2006; Burrows et al. 2007) neutrino emission far outside of this prompt window. In Section 7 we place limits on two classes of GRB prompt neutrino flux predictions: models normalized to the observed UHECR flux (Katz et al. 2009) and models normalized to the observed gamma-ray flux for each burst.

The cosmic-ray-normalized models (Waxman & Bahcall 1997; Waxman 2003; Ahlers et al. 2011) assume that protons emitted by GRBs are the dominant sources of the highest energy cosmic rays observed, and with these models we place limits on this assumption. The neutrino spectral break and normalization depend on the bulk Lorentz boost factor Γ and gamma-ray break energy in the GRB fireball. Typical values required for pion production lead to the emission of ∼100 TeV neutrinos. In Section 7, we place limits on the expected neutrino flux at different break energies.

The gamma-ray-normalized models (Hümmer et al. 2012; Zhang & Kumar 2013) do not relate the observed cosmic ray flux to neutrinos produced in GRB fireballs, and with these models we place limits on internal fireball parameters. We consider three types of gamma-ray-spectrum-normalized fireball models, calculated on a burst-by-burst basis, that differ in their neutrino emission sites. The internal shock model relates the neutrino production radius to the variability timescale of the gamma-ray light curves (Waxman & Bahcall 1997; Hümmer et al. 2012; Zhang & Kumar 2013). The photospheric model places the radius at the photosphere through combinations of processes such as internal shocks, magnetic reconnection, and neutron-proton collisions (Rees & Mészáros 2005; Murase 2008; Zhang & Yan 2011). The internal collision-induced magnetic reconnection and turbulence (ICMART) model favors a neutrino production radius ∼10 times larger than the standard internal shock model due to a Poynting-flux-dominated outflow that remains undissipated until internal shocks destroy the ordered magnetic fields (Zhang & Yan 2011; Zhang & Kumar 2013). These models assert that gamma-ray emission and proton acceleration occur at the same radius. This equivalence is not necessarily true for scenarios other than the single-zone internal shock model (Bustamante et al. 2015a), but allows the predicted neutrino flux to scale linearly with the proton-to-electron energy ratio in the fireball. Additionally, the models addressed in this work do not account for a possible enhancement to the high energy neutrino flux due to acceleration of secondary particles (Klein et al. 2013; Winter et al. 2014).

We calculate the per-GRB predictions normalized to the measured gamma-ray spectra numerically with a wrapper of the Monte-Carlo generator SOPHIA (Mücke Reimer et al. 2014), taking into account the full particle production chain and synchrotron losses in inelastic $p\gamma$ interactions. For these calculations, we parametrize the reported gamma-ray spectrum of each GRB as a broken power-law approximation of the Band function (Band et al. 1993; Zhang & Kumar 2013). We retrieve GRB parameters from the circulars published by satellite detectors on the Gamma-ray Coordinates Network⁵³ and the Fermi GBM database (Gruber et al. 2014; von Kienlin et al. 2014). We compile the relevant temporal, spatial, and spectral parameters used in this analysis on our GRBweb database, presented on a publicly available website⁵⁴ (Aguilar 2011). The prompt photon emission time (T₁₀₀) is defined by the most inclusive start and end times (T₁ and T₂) reported by any satellite. We use the most precise localization available. Following the same prescription of our previous model limit calculations (Abbasi et al. 2011b, 2012a; Aartsen et al. 2015d), if the fluence is unmeasured, we use an average value of 10⁻⁵ erg cm⁻²; if the gamma-ray break energy is unmeasured, we use 200 keV for ${T}_{100}\gt 2\;{\rm{s}}$ bursts and 1000 keV for ${T}_{100}\leqslant 2\;{\rm{s}}$ bursts; and if the redshift is unmeasured, we use 2.15 for ${T}_{100}\gt 2\;{\rm{s}}$ bursts and 0.5 for ${T}_{100}\leqslant 2\;{\rm{s}}$ bursts.

The neutrino flux predictions depend on several unmeasured quantities; we use a variability timescale of 0.01 s and an isotropic luminosity of 10⁵² erg cm⁻² for long bursts and a variability timescale of 0.001 s and an isotropic luminosity of 10⁵¹ erg cm⁻² for short bursts, which are consistent with the literature (Hümmer et al. 2012; Zhang & Kumar 2013; Baerwald et al. 2014). If the redshift is known for a particular burst, we calculate the approximate isotropic luminosity from the redshift, photon fluence, and T₁₀₀ (Hümmer et al. 2012).

Figure 1 illustrates neutrino spectra from the three models with benchmark fireball parameters. These benchmark parameters are bulk Lorentz boost factor ${\rm{\Gamma }}=300$ and proton-to-electron energy ratio, or baryonic loading, f_p = 10. These spectra are presented as per-flavor quasi-diffuse fluxes, in which we divide the total fluence from all GRBs in the sample by the full sky $4\pi$ steradians and one year in seconds, and scale the total number of bursts to a predicted average 667 observable bursts per year, which has been used in our previous IceCube publications (Abbasi et al. 2010, 2011b, 2012a; Aartsen et al. 2015d). The actual number of bursts observed by satellite detectors in each year is less than the predicted average because of detector field of view limitations and obstruction by the Sun, Earth, and Moon. For all GRBs we assume that neutrino oscillations over cosmic baselines result in equal fluxes for all flavors (Aartsen et al. 2015a, 2015c; Bustamante et al. 2015b; Palladino et al. 2015; Palomares-Ruiz et al. 2015).

Zoom In Zoom Out Reset image size

The IceCube detector (Achterberg et al. 2006) consists of 5160 digital optical modules (DOMs) instrumented over 1 km³ of clear glacial ice 1450 to 2450 m below the surface at the geographic South Pole. IceCube is the largest neutrino detector in operation. Each DOM incorporates a 10 in. diameter photomultiplier tube (PMT) (Abbasi et al. 2010). Signal and power connections between the DOMs and the computer farm at the surface are provided by 86 vertical in-ice cables or "strings" that each connect 60 DOMs spaced uniformly. Adjacent strings are separated by about 125 m. The DeepCore array (Abbasi et al. 2012b) is made up of a more densely spaced subset of strings that are located in the clearest ice at depths below 2100 m and contain higher quantum efficiency PMTs.

Sensor deployment began during the 2004–2005 austral summer. Physics data collection began in 2006 with the 9-string configuration and continued with partial detector configurations through completion of the 86 strings in 2010 December. We conducted the analysis presented in this paper using data taken from 2010 May through 2013 May, with one year using the 79-string configuration and two years using the completed detector.

The DOM PMTs detect the Cherenkov radiation of relativistic charged particles produced in deep inelastic neutrino-nucleon scattering in or near the detector volume. Data acquisition (Abbasi et al. 2009) begins once the output current exceeds the threshold of 1/4 of the mean peak current of the pulse amplified from a single photo-electron. If another such "hit" is recorded in a neighboring or next-to-neighboring DOM within 1 μs, the full waveform information is recorded. DOMs with hits failing this local coincidence condition report a short summary of their recorded waveform for inclusion in data records. The digitized waveforms are sent to the computers at the surface, where they are assembled into "events" by a software trigger (Abbasi et al. 2009; Kelley et al. 2014). The relative timing resolution of photons within an event is about 2 ns (Achterberg et al. 2006). Data are sent daily from the detector to computers in the Northern Hemisphere via satellite.

There are two main interaction topologies that can manifest in hit DOMs during an event: a thin track or a near-spherical shower. Track-like light emission is generated by muons created either through cosmic ray interactions in the atmosphere or ${\nu }_{\mu }$ interacting in the ice. Shower-like light emission, which is the signal of this analysis, is generated by electromagnetic cascades from charged current interactions of ${\nu }_{e}$ and ${\nu }_{\tau }$ and hadronic cascades from neutral current interactions of ${\nu }_{e}$, ${\nu }_{\tau }$, and ${\nu }_{\mu }$ in the ice. ${\nu }_{\alpha }$ stands for ${\nu }_{\alpha }+{\bar{\nu }}_{\alpha }$ unless otherwise specified.

Electromagnetic cascades of photons, electrons, and positrons develop through radiative processes, such as bremsstrahlung and pair production. Once the energy of each particle in the cascade is below the critical energy, non-radiative processes, such as ionization and Compton scattering dominate and the shower stops growing. Hadronic cascades result from showers of baryons and mesons produced by deep inelastic scattering interactions of all neutrino types. Interactions of ${\bar{\nu }}_{e}$ with electrons via the Glashow resonance at 6.3 PeV would also produce cascade signal for this analysis; however, this resonance has not yet been observed.

We use Monte Carlo simulations of neutrinos interacting in the IceCube detector for the signal hypothesis in our event selection and optimization for this search. Neutrinos are simulated with the NEUTRINO-GENERATOR program, a port of the ANIS code (Gazizov & Kowalski 2005). We use NEUTRINO-GENERATOR to distribute neutrinos with a power-law spectrum uniformly over the entire sky and then we propagate them through the Earth and ice. The simulated neutrino-nucleon interactions take cross sections from CTEQ5 (Lai et al. 2000). The Earth's density profile is modeled with the Preliminary Reference Earth Model (Dziewonski & Anderson 1981). The propagation code takes into account absorption, scattering, and neutral-current regeneration.

Although we use data outside of GRB gamma-ray emission time windows for our background, simulated neutrinos and muons generated in cosmic ray air showers are useful checks to characterize our background and estimate the signal purity of our final data sample. We apply the Honda et al. spectrum (Honda et al. 2007) for the atmospheric neutrino background. We use the CORSIKA simulation package (Heck et al. 1998) to simulate cosmic ray air showers. Muons are traced through the ice and bedrock incorporating continuous and stochastic energy losses (Chirkin & Rhode 2004). The PMT detection of Cherenkov light from muon tracks and showers is simulated using ice and dust layer properties determined in detailed studies and simulations (Lundberg et al. 2007) (Ackermann et al. 2006). Finally, we simulate the DOM triggering and signal from all interactions. We process these signals in the same way as described in Section 3.

We use the geometric pattern, timing, and amount of recorded photons in an interaction to reconstruct and identify the incident particle. We run simple analytic reconstructions automatically in data filters at the Pole and more complex reconstructions on the reduced data in the north. These reconstructions provide powerful discriminating features that we employ in our event selection methods, described in the next section.

One of the analytic calculations used to identify shower-like events is the "tensor-of-inertia" algorithm of Ahrens et al. (2004) that determines an interaction vertex by calculating a "center of mass," where the mass terms correspond to the number of photoelectrons recorded by each DOM. The eigenvector along the longest principle axis provides a reasonable guess for the incident particle trajectory. Using the rigid body analogy, a spherical shower should provide nearly even tensor-of-inertia eigenvalues, while an elongated track should have one eigenvalue much smaller than the other two.

A second analytic calculation is the "line-fit" algorithm of Ahrens et al. (2004) that ignores the geometry of the Cherenkov cone and the optical properties of the medium and assumes light travels along a one-dimensional path at some constant velocity through the detector. This algorithm determines a best-fit track, but provides a good indicator for spherical showers, which exhibit slower best-fit velocities than tracks. We recalculate this fit during the final level of event selection using an improved version of the algorithm that discards photons due to uncorrelated noise (Aartsen et al. 2014a).

We also perform fast reconstructions based on an analytic likelihood function. In general, we determine the set of parameters ${\boldsymbol{a}}=\{{a}_{1},{a}_{2},\ldots ,{a}_{m}\}$ of an event that maximizes the likelihood function

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{x}}| {\boldsymbol{a}})=\displaystyle \prod _{i=1}^{N}p({{\boldsymbol{x}}}_{i}| {\boldsymbol{a}})\end{eqnarray} \tag{ 2 }$$

given sets of observables ${\boldsymbol{x}}=\{{{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2},\ldots ,{{\boldsymbol{x}}}_{N}\}$ from N DOMs. We minimize the negative logarithm of ${ \mathcal L }$ with respect to ${\boldsymbol{a}}$. In these online likelihood reconstructions, the ${\boldsymbol{x}}$ are parametrized in terms of the difference between the recorded photon hit times and the expectation without scattering or absorption (Ahrens et al. 2004).

For a cascade hypothesis, which can be considered point-like because of the relatively small maximum shower range on the scale of IceCube instrumentation (Aartsen et al. 2014b), the time residual depends on the location and time of the interaction vertex ${\boldsymbol{a}}=\{t,x,y,z\}$. This minimization is seeded by the above center of mass vertex calculation and the earliest vertex time for which at least 4 DOMs record hits with time residuals within [0, 200] ns. For an infinite track hypothesis, the time residual depends on the distance from an arbitrary point on the track and the direction of the track ${\boldsymbol{a}}=\{x,y,z,\theta ,\phi \}$. This minimization is seeded by the above line-fit track. Both of these cascade and track hypothesis maximum likelihoods provide useful best-fit quantities for the event selection described in Section 5.

The more computationally expensive likelihood-based reconstruction that we use on the reduced data takes advantage of detailed in-ice photon propagation models (Lundberg et al. 2007; Whitehorn et al. 2013) and DOM waveform information. The detection of light follows a Poisson distribution and thus the likelihood function we maximize is

$$\begin{eqnarray}&&{ \mathcal L }(k| E,t,x,y,z,\theta ,\phi )=\displaystyle \prod _{i=1}^{{N}_{\mathrm{DOM}}}\displaystyle \frac{{\lambda }_{i}^{{k}_{i}}}{{k}_{i}!}{e}^{-{\lambda }_{i}}\end{eqnarray} \tag{ 3 }$$

for the observed and expected number of photoelectrons k and λ in each DOM from a neutrino traveling in direction $(\theta$, $\phi )$ and depositing energy E at time t and vertex $(x,y,z)$ (Aartsen et al. 2014b). This likelihood generally considers multiple light sources, e.g., for stochastic muon losses along a path through the detector. Because we define our signal to be neutrino-induced showers, it is sufficient to use a point-like hypothesis. Due to the linear relationship between Cherenkov light emission and deposited energy in ${\nu }_{e}$ charged-current interactions, we estimate the electromagnetic-equivalent energy of an event by comparing to a template electromagnetic cascade of 1 GeV (Aartsen et al. 2014b). By using Cherenkov light source-dependent spline-fit tables of photon amplitudes and time delays, we minimize (James & Roos 1975) the negative log-likelihood and extract the best-fit parameters of the particle that produced the shower.

We iterate this minimization several times and achieve improved angular reconstruction from the best-fit energy and modeled ice structure. The reconstruction is seeded by the tensor-of-inertia direction, online likelihood vertex, and best-fit energy from a simple isotropic shower and coarsely binned propagation model expectation in Equation (3). We use this reconstruction's best-fit direction and energy for each event in our likelihood analysis. We estimate a lower bound on the error in the reconstructed directions with the Cramer–Rao relation (Cramer 1945; Rao 1945) between the covariance of each fit parameter and the inverted Fisher information matrix. Figure 2 shows the angular and energy resolutions using this reconstruction at the final event selection of this cascade search. The reconstructed energy for neutral current interactions is less than the primary neutrino energy because of the dissipation to outlets without Cherenkov emission.

Zoom In Zoom Out Reset image size

Our signal in this search is one or more high energy neutrino-induced showers coincident in space and time with one or more GRBs. Before we analyze the likelihood that an event is a GRB-emitted neutrino, we must reduce the triggered data that are dominated by nearly 3 kHz of cosmic ray air shower muons to a much smaller selection of possible signal. Our event selection is optimized on a general ${E}^{-2}$ spectrum of diffuse astrophysical simulated ${\nu }_{e}$ for signal and off-time data for background; we choose "off-time" as data not within two hours of any reported GRB γ emission in order to ensure blindness to potential signal. We use the reconstructions detailed in Section 4 to remove the background and realize our final sample.

Searches for neutrino-induced showers from astrophysical and atmospheric sources have been conducted previously in IceCube (Abbasi et al. 2011a; Aartsen et al. 2014c, 2014e, 2015b). The predominant difference between these and the search presented in this paper is that we assume neutrinos come from known transient sources. The previous shower-like event selections assume a diffuse signal or constantly emitting sources, and, as a result, require much more stringent background reduction that leads to data rates nearly a factor of 100 smaller than what is needed in this search. Cascade containment constraints in the detector were imposed to reach these low backgrounds, which are achieved in this search by the effective cuts in time and space around each GRB in the unbinned likelihood analysis presented in Section 6.

The first class of background to remove is track-like events generated by muons losing their energy through continuous ionization processes. These muon tracks are relatively easily separated from our shower-like neutrino signal by use of a filter run online at the detector site. During the 79-string configuration, we use the tensor-of-inertia eigenvalue ratio and line-fit absolute speed to choose spherical events. The online cascade filter for the first two years of the 86-string configuration imposes a zenith-dependent selection to allow more events with energies below 10 TeV. We use the reduced log-likelihood calculated in the fast cascade likelihood reconstruction to remove a large fraction of the muons that are misreconstructed as upgoing. The down-going region requires more restrictive reduced log-likelihood values in combination with the same online 79-string event selection from above.

Upon transmission via satellite from the South Pole, these ∼30 Hz filtered data of spherically shaped events must be further reduced to have optimal performance with the boosted decision tree (BDT) forest algorithm described below. The background at this rate is still dominated by atmospheric muons that lose their energy through bremsstrahlung, photonuclear interactions, and pair production. These stochastic energy loss mechanisms create nearly spherical hit patterns that are difficult to differentiate from our signal when the muon track is at the edge or outside of the detector and therefore not observed. We derive variables from the more CPU-intensive reconstructions run offline to separate these muons from neutrino-induced showers at this selection level, called level 3 in this text, and for eventual machine-learning input.

Even though the background muons in the online filter event selection are shower-like in topology, we exploit the Cherenkov light hit patterns of the minimally ionizing muon track when possible to differentiate from showers produced by neutrinos. We perform this further separation of signal from muon background first by using the fast track hypothesis and cascade hypothesis analytic likelihood function reconstructions, described in Section 4. Our first level 3 selection is on events for which the ratio of the track likelihood to the cascade likelihood heavily favors the cascade hypothesis. Two other background muon features that we take advantage of are their down-going directions and typically lower energies. We thus parametrize the track likelihood reconstructed zenith in terms of reconstructed energy and remove low energy down-going events from our sample. Lastly, at this rate, we choose events based on the fraction of hit DOMs that are inside of a sphere centered on the reconstructed vertex with a radius determined by the mean hit distance. We optimize the radius and fraction filled of the sphere for reconstructed vertices inside and outside of the instrumented volume and use both volume regions for this neutrino search.

Our final event selection employs a machine learning algorithm to optimize the separation of signal and background over a space of many features. The algorithm we use is a BDT forest (Freund & Schapire 1997) that also has been used in Northern Hemisphere ${\nu }_{\mu }$ track GRB coincidence searches in IceCube (Abbasi et al. 2012a; Aartsen et al. 2015d). We built a collection of variables, many of which were influenced by past neutrino-induced shower searches conducted in IceCube (Abbasi et al. 2011a; Aartsen et al. 2014c), to train the BDT to discriminate between signal and background events. During this training, the algorithm determines the variable value that best separates signal and background at each node of each tree. Each successive tree increases, or "boosts," the weights of incorrectly classified events from the prior tree, allowing more ambiguous events to be classified correctly. The trees and forest grow until certain stopping criteria are reached. The signal/background $(+1/-1)$ score of each event that we use in our final selection is a weighted average of its scores over all trees, where this weight is related to the boost factor of each tree. For each of the three year-long detector configurations of this search, we train a separate BDT with configuration-specific signal simulation and background data.

Our BDT discrimination variables take advantage of topological and energetic differences between astrophysical neutrino and atmospheric muon spectra. Figures 3 and 4 show the distributions for simulated astrophysical neutrinos, atmospheric neutrinos, atmospheric muons, and data for three of the most powerful variables used by the BDT after level 3 and the final level, respectively. Some deviations in variable values between measured and simulated data are visible in Figures 3 and 4, and are attributed to the simulation's limited ability to produce all observed events. These deviations do not affect the analysis because the background is determined from measured data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the distributions of simulation and data with respect to BDT score. The vertical dashed line corresponds to our optimized final event selection described in the next section. The atmospheric and astrophysical ${\nu }_{\mu }$ in these plots is already preselected to be shower-like and has minimal overlap with the track search. Background is reduced by a factor of $6\times {10}^{-6}$ from the online filter level to the BDT-based final selection. Figure 6 shows the ${\nu }_{e}$ signal efficiency at this final level with respect to the online filter rate as a function of energy. The slight decrease in efficiency for the highest energy signal is due to events for which the computationally expensive reconstruction, described at the end of Section 4 and used in the likelihood analysis, fails to converge. These simulated neutrinos interact outside of the detector, account for about 1% of the ${E}^{-2}$ signal, and are discarded prior to the BDT selection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the most effective BDT variables is the same ratio of the track likelihood to the cascade likelihood used in level 3. The high energy background muons that passed in spite of the selection at level 3 are distinguished from signal with the BDT algorithm. The large likelihood ratios, seen in Figures 3 and 4, are due to the highest energy signal events with the most hit DOMs. Another potent variable is the total amount of Cherenkov light imparted in the DOMs divided by the tensor-of-inertia derived elongation of an event. The numerator separates lower energy background atmospheric muons from higher energy astrophysical signal neutrinos, while the denominator separates elongated background atmospheric muons from spherical neutrino-induced showers. As shown in Figure 3, the ratio of these two observables effectively separates the lower energy atmospheric neutrino background from our signal as well. A third useful variable is the reduced negative log-likelihood value from the analytic cascade hypothesis reconstruction.

Compared to the event samples of IceCube's searches for ${\nu }_{\mu }$-induced tracks from Northern Hemisphere GRBs, our backgrounds in this all-sky all-flavor cascade search require an event selection to a 10 times lower data rate of 0.16 mHz averaged over the three search years. The differences between neutrino-induced showers and muon-induced stochastic energy loss showers are less apparent than the differences between neutrino-induced tracks and detector-edge and coincident muon-induced tracks incorrectly reconstructed as upgoing. In the Northern Hemisphere track searches, most of these muons are able to be identified and removed, leaving only atmospheric neutrinos. As a result, the atmospheric neutrino purity with respect to muons of this all-sky all-flavor search is ∼40%, and is significantly less than the ∼90% purity of the Northern Hemisphere track searches. Nevertheless, we achieve similar sensitivity to the track search through our acceptance of ${\nu }_{e}$, ${\nu }_{\tau }$, and ${\nu }_{\mu }$ signal from GRBs over the entire sky. Figure 7 compares the effective areas for the two different GRB neutrino searches.

Zoom In Zoom Out Reset image size

Once we have selected a sample of events that resemble high energy neutrino-induced electromagnetic or hadronic showers, we must calculate the likelihood that these events are neutrino signal from observed GRBs. We use a likelihood function that incorporates the probabilities that an event is a signal neutrino from a GRB or a background atmospheric neutrino or muon. We calculate these background-like and signal-like probabilities from individually normalized probability distribution functions (PDFs) in time, space, and energy:

$$\begin{eqnarray*}&&{ \mathcal S }({{\boldsymbol{x}}}_{i})={P}_{{\rm{s}}}^{\mathrm{Time}}({t}_{i})\times {P}_{{\rm{s}}}^{\mathrm{Space}}({{\boldsymbol{r}}}_{i})\times {P}_{{\rm{s}}}^{\mathrm{Energy}}({E}_{i})\end{eqnarray*}$$

$$\begin{eqnarray*}&&{ \mathcal B }({{\boldsymbol{x}}}_{i})={P}_{{\rm{b}}}^{\mathrm{Time}}({t}_{i})\times {P}_{{\rm{b}}}^{\mathrm{Space}}({{\boldsymbol{r}}}_{i})\times {P}_{{\rm{b}}}^{\mathrm{Energy}}({E}_{i}),\end{eqnarray*}$$

where S and B are the probabilities that an event i with properties ${{\boldsymbol{x}}}_{i}$ is signal and background, respectively.

The signal time PDF is flat during the gamma-ray emission duration (T₁₀₀ defined in Section 2) and has Gaussian tails before T₁ and after T₂. The width of these Gaussian tails ${\sigma }_{{\rm{t}}}$ equals the duration of measured gamma-ray emission up to 30 s and down to 2 s to account for possible small shifts in the neutrino emission time with respect to that of the photons. We accept events out to ${T}_{100}\pm 4{\sigma }_{{\rm{t}}}$ for each GRB time window. The background time PDF is flat throughout the total period of acceptance for each GRB. Examples of the signal/background time PDF ratios for events during GRBs with different measured durations are shown in Figure 8.

Zoom In Zoom Out Reset image size

The signal space PDF is the first-order non-elliptical term in the Kent distribution (Kent 1982):

$$\begin{eqnarray}&&{P}_{{\rm{s}}}^{\mathrm{Space}}({{\boldsymbol{r}}}_{i})=\displaystyle \frac{\kappa }{4\pi \mathrm{sinh}(\kappa )}{e}^{\kappa ({\hat{r}}_{i}\cdot {\hat{r}}_{\mathrm{GRB}})},\end{eqnarray} \tag{ 4 }$$

for which the concentration parameter $\kappa =\tfrac{1}{{\sigma }_{\mathrm{GRB}}^{2}+{\sigma }_{i}^{2}}$ is the reciprocal of the uncertainty in the GRB's localization and the Cramer–Rao uncertainty in the event's reconstructed direction, ${\hat{r}}_{i}$ is the reconstructed direction of the event, and ${\hat{r}}_{\mathrm{GRB}}$ is the most precise GRB localization available. We use this spherical analog of the planar Gaussian distribution because of the large uncertainties in our reconstructed cascade event directions. The background space PDF is a spline fit to the distribution of reconstructed $\mathrm{cos}({\theta }_{\mathrm{zenith}})$ of all final analysis level off-time data events. Small variance in the background reconstructed azimuth distribution has negligible impact on this analysis, so our background space PDF only varies with zenith. Examples of the signal space PDFs for events and correlated GRBs with different directional uncertainties are shown with the background space PDF in Figure 9.

Zoom In Zoom Out Reset image size

The signal and background energy PDFs are the reconstructed energy distributions of ${E}^{-2}$ ${\nu }_{e}$ simulation and off-time data, respectively. We fit a spline to the ratio of these two PDFs. We find few background events in our final sample that have reconstructed energies above 1 PeV, so we conservatively assume a constant ratio of signal and background energy PDFs at energies above 1 PeV. The signal and background energy PDFs and their spline-fit ratio are shown in Figure 10.

Zoom In Zoom Out Reset image size

In order to choose our optimal final selection on BDT score and characterize the significance of the result, we construct a test statistic in the form of a maximum likelihood function. This function incorporates probabilities that observed events are signal and background and provides an estimator for the number of observed signal events. The likelihood function combines the above signal and background PDFs with the Poisson probability P_N of observing N events, given that the expected total number of signal + background events is ${n}_{{\rm{s}}}+{n}_{{\rm{b}}}$:

$$\begin{eqnarray}&&{ \mathcal L }(\{{{\boldsymbol{x}}}_{i}\};{n}_{{\rm{s}}}+{n}_{{\rm{b}}})={P}_{N}\displaystyle \prod _{i=1}^{N}({p}_{{\rm{s}}}{ \mathcal S }({{\boldsymbol{x}}}_{i})+{p}_{{\rm{b}}}{ \mathcal B }({{\boldsymbol{x}}}_{i})),\end{eqnarray} \tag{ 5 }$$

where

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

and ${p}_{{\rm{s}}}=\tfrac{{n}_{{\rm{s}}}}{{n}_{{\rm{s}}}+{n}_{{\rm{b}}}}$ and ${p}_{{\rm{b}}}=\tfrac{{n}_{{\rm{b}}}}{{n}_{{\rm{s}}}+{n}_{{\rm{b}}}}$ are the probabilities of observing a signal and background event.

We do not know a priori the total number of events we will observe. Therefore, we estimate n_b from the background data rate multiplied by the total search time window ${\sum }_{i}^{{N}_{\mathrm{GRBs}}}({T}_{100,i}+8{\sigma }_{{\rm{t}},i})$ in which we accept events; and this expectation is denoted by $\langle {n}_{{\rm{b}}}\rangle$. We estimate ${n}_{{\rm{s}}}$ from the value ${\hat{n}}_{{\rm{s}}}$ that maximizes ${ \mathcal L }$. Additionally, we simplify the likelihood function by dividing by its background-only hypothesis and taking the logarithm of this ratio without loss of generality. Finally, our test statistic T is the maximized value of this ratio at ${\hat{n}}_{{\rm{s}}}$:

$$\begin{eqnarray}&&T=-{\hat{n}}_{{\rm{s}}}+\displaystyle \sum _{i=0}^{N}\mathrm{ln}\left[\displaystyle \frac{{\hat{n}}_{{\rm{s}}}{ \mathcal S }({{\boldsymbol{x}}}_{i})}{\langle {n}_{{\rm{b}}}\rangle { \mathcal B }({{\boldsymbol{x}}}_{i})}+1\right].\end{eqnarray} \tag{ 7 }$$

We extend to multiple detector configurations (79 strings, first year of 86 strings, second year of 86 strings, etc.) and search channels (cascades, tracks) by adding maximized test statistics for each configuration and channel:

$$\begin{eqnarray}&&T=\displaystyle \sum _{{\rm{c}}}\left\{-{({\hat{n}}_{{\rm{s}}})}_{{\rm{c}}}+\displaystyle \sum _{i=0}^{{N}_{{\rm{c}}}}\mathrm{ln}\left[\displaystyle \frac{{({\hat{n}}_{{\rm{s}}})}_{{\rm{c}}}{{ \mathcal S }}_{{\rm{c}}}({{\boldsymbol{x}}}_{i})}{\langle {n}_{{\rm{b}}}{\rangle }_{{\rm{c}}}{{ \mathcal B }}_{{\rm{c}}}({{\boldsymbol{x}}}_{i})}+1\right]\right\},\end{eqnarray} \tag{ 8 }$$

where c represents each combination of search channel and detector configuration.

To set discovery significance thresholds, we first perform 10⁸ pseudo-search trials using only background data for a range of BDT score $\gt s$ cuts. Each background sample has its own $\langle {n}_{{\rm{b}}}{\rangle }_{\mathrm{cut}}$, with lower values for cuts on higher BDT scores. In each background-only trial, for each GRB, we choose a pseudo-random number of observed events within our ${T}_{100}\pm 4{\sigma }_{{\rm{t}}}$ search time window about the gamma-ray emission from a Poisson distribution with expectation determined by the background data rate and time window. If the number of events is 0, then T receives no contribution from the search window about that GRB. If the number of events is greater than 0, then we construct each event using the following steps: (1) choose a random time PDF value; (2) choose a random azimuth within 0 to $2\pi ;$ (3) choose a reconstructed energy by sampling from the background distribution; (4) choose a reconstructed zenith by sampling from the background distribution of events with similar energy; (5) choose an estimated error in reconstructed direction by sampling from the background distribution of events with similar energy and zenith. Finally, with the signal and background PDF values for every event, we calculate the test statistic for each trial. We obtain a distribution like that shown in Figure 11 for each possible selection on BDT score.

Zoom In Zoom Out Reset image size

We choose our optimal final selection on BDT score by injecting simulated neutrino signal along with background data and performing 10⁴ pseudo-search trials for a range of BDT score $\gt s$ cuts. The background events are selected for each GRB the same as above for the background-only trials. The simulated signal events within an 11° circle about each GRB contribute to the likelihood with a probability proportional to their simulated weight. Signal is increasingly weighted until a defined discovery or limit-setting threshold is reached. These discovery and limit-setting potentials per cut on BDT score are shown in Figure 12 for a general ${E}^{-2}$ spectrum and the Figure 1 benchmark internal shock fireball spectrum. The final selection, BDT score $\gt 0.525$, was optimized to set the best possible upper limit while suffering little loss in discovery potential. This optimization was performed for each search year's event selection.

Zoom In Zoom Out Reset image size

In three years of data, five cascade neutrino candidate events are correlated with five GRBs. These coincidences yield a combined P-value of 0.32, derived from the total test statistic value with respect to the background-only distribution. The Northern Hemisphere track searches in four years of data resulted in a single neutrino candidate event correlated with prompt GRB emission and a combined P-value of 0.46 (Aartsen et al. 2015d). Considering the atmospheric neutrino purity of each search, discussed in Section 5, the track event is almost certainly a ${\nu }_{\mu }$ while the cascade events could be high energy atmospheric muons.

Table 1 shows the time, space, and energy data for these events and GRBs. Events 1, 2, and 4 occurred on the edge or corner of the detector. Events 3 and 5 occurred within the instrumented volume. The large directional uncertainties of the second and fourth events are due to the location of their interactions at the corner and edge of the instrumented volume, with relatively few DOMs able to record the Cherenkov light. The reconstructed energy of the muon neutrino track coincidence discussed in Aartsen et al. (2015d) is an approximate lower bound on the neutrino energy because the interaction could have taken place well outside of IceCube and consequently, the muon could have lost significant energy before reaching the detector.

Table 1.  GRB and Event Properties for the Three-year Cascade and Four-year Track Search Coincidences

	Time	Angular Uncertainty	Angular Separation	Fluence/Energy
GRB 101213A	${T}_{100}=202\;{\rm{s}}$	$0\buildrel{\circ}\over{.} 0005$		$7.4\times {10}^{-6}$ erg cm−2
Event 1	${T}_{1}+109\;{\rm{s}}$	$32\buildrel{\circ}\over{.} 0$	23°	11 TeV
GRB 110101B	${T}_{100}=235\;{\rm{s}}$	$16\buildrel{\circ}\over{.} 5$		$6.6\times {10}^{-6}$ erg cm−2
Event 2	${T}_{1}+141\;{\rm{s}}$	118°	112°	34 TeV
GRB 110521B	${T}_{100}=6.14\;{\rm{s}}$	$1\buildrel{\circ}\over{.} 31$		$3.6\times {10}^{-6}$ erg cm−2
Event 3	${T}_{1}+0.26\;{\rm{s}}$	$16\buildrel{\circ}\over{.} 5$	35°	3.4 TeV
GRB 111212A	${T}_{100}=68.5\;{\rm{s}}$	$0\buildrel{\circ}\over{.} 0004$		$1.4\times {10}^{-6}$ erg cm−2
Event 4	${T}_{1}+11.7\;{\rm{s}}$	$44\buildrel{\circ}\over{.} 8$	120°	31 TeV
GRB 120114A	${T}_{100}=43.3\;{\rm{s}}$	$0\buildrel{\circ}\over{.} 04$		$2.4\times {10}^{-6}$ erg cm−2
Event 5	${T}_{1}+57.2\;{\rm{s}}$	$7\buildrel{\circ}\over{.} 9$	18°	3.8 TeV
GRB 100718Aa	${T}_{100}=39\;{\rm{s}}$	$10\buildrel{\circ}\over{.} 2$		$2.5\times {10}^{-6}$ erg cm−2
${\nu }_{\mu }$ Track Eventa	${T}_{1}+15\;{\rm{s}}$	$1\buildrel{\circ}\over{.} 3$	16°	≳ 10 TeV

Note.

^aCorresponds to the ${\nu }_{\mu }$ track search coincidence discussed in Aartsen et al. (2015d).

Download table as:  ASCIITypeset image

The benchmark internal shock, photospheric, and ICMART fireball model spectra plotted in Figure 1 yield expectations of 1.5, 2.5, and 0.07 neutrinos, respectively, for this three-year all-flavor cascade analysis. As discussed at the end of Section 5, this search achieves nearly the same sensitivity as the Northern Hemisphere track search. Thus, the combined four years of track searches and three years of shower searches yield total benchmark model expectations of 3.3, 5.4, and 0.1 neutrinos. The all-flavor shower search has an average $\langle {n}_{{\rm{b}}}\rangle$ of 11 events per year. This expectation is more concentrated at lower energies than the expected signal and is weighted by the time, space, and energy PDFs accordingly in our unbinned likelihood. The expected background during just the T₁₀₀ of each GRB is 3.5 events per year. An observation of three 1 PeV neutrinos correlated with three GRBs, with the same temporal and spatial properties as our observed Event 1 and its respective GRB in Table 1, would yield a T value over 12 and a $5\sigma$ discovery based on the background-only distribution in Figure 11.

Considering the low significance of the events correlated with GRBs, we set 90% confidence level (CL) Neyman upper limits (Neyman 1937) on models normalized to the observed flux of UHECRs as well as models normalized to the observed gamma-ray fluence of each GRB. We calculate these limits by combining the three-year cascade search results and four-year Northern Hemisphere track search results using the multiple configuration and channel test statistic given in Equation (8). Simulated neutrinos are weighted to a certain spectrum and normalization and injected over background in the pseudo-searches. The exclusion CL is the fraction of pseudo-search trials that yield $T\geqslant {T}_{\mathrm{observed}}$.

Figure 13 shows exclusion contours for double broken power law spectra of the form

$$\begin{eqnarray}{{\rm{\Phi }}}_{\nu }(E)={\phi }_{0}\times \left\{\begin{array}{ll}{E}^{-1}{\epsilon }_{b}^{-1}, & E\lt {\epsilon }_{b}\\ {E}^{-2}, & {\epsilon }_{b}\leqslant E\lt 10{\epsilon }_{b}\\ {E}^{-4}{(10{\epsilon }_{b})}^{2}, & 10{\epsilon }_{b}\leqslant E\end{array}\right.\end{eqnarray} \tag{ 9 }$$

at different first break energies ${\epsilon }_{b}$ and normalizations ${\epsilon }_{b}^{2}{\phi }_{0}$. Our combined limits largely rule out cosmic ray escape via neutron production (Ahlers et al. 2011). Mechanisms allowing for cosmic ray escape via protons (Waxman & Bahcall 1997) are disfavored as well. The Waxman & Bahcall (1997) model has been updated to account for more recent measurements of the UHECR flux (Katz et al. 2009) and typical gamma break energy (Goldstein et al. 2012).

Zoom In Zoom Out Reset image size

Figure 14 shows exclusion contours in the baryonic loading—bulk Lorentz factor parameter space for the internal shock, photospheric, and ICMART per-burst gamma-ray-normalized fireball models. The benchmark model spectra from Figure 1 are indicated by the intersection of the vertical and horizontal dashed lines. Systematic uncertainties in the ice model, DOM photon detection efficiency, and lepton propagation in the Earth and ice are added in quadrature to a total effect of $\sim 11 \%$ for each model limit. The addition of the all-sky cascade search over three years of data and GRBs improves the Northern Hemisphere track search limits by about 50%.

Zoom In Zoom Out Reset image size

We have performed a search for neutrinos of all flavors emitted by 807 GRBs during 3 years of IceCube data. This search exhibits similar sensitivity to that of previously published searches for muon neutrinos emitted by 506 Northern Hemisphere GRBs during 4 years of IceCube data. Over both search channels, we found six events that are correlated with GRBs but also consistent with background. Our limits placed on neutrino emission models normalized to the observed UHECR flux are the strongest constraints thus far on the hypothesis that GRBs are the dominant sources of this flux. Additionally, our limits placed on the latest neutrino emission models normalized to the observed gamma-ray fluence from each GRB constrain parts of the parameter space relevant for the production of UHECR protons. Models that are still allowed require increasingly lower neutrino production efficiencies through large bulk Lorentz boost factors, low baryonic loading, or large dissipation radii.

As shown in Baerwald et al. (2014), constraints on parameters involved in fireball neutrino production via internal shock collisions can be connected to the requirement that GRBs are the sources of the observed UHECRs in a self-consistent way, assuming a pure proton composition. While cosmic ray escape from GRB fireballs through neutron production is strongly constrained by our limits, Baerwald et al. (2014) show that the unexcluded parameter space in Figure 14 allows for efficient diffusive proton escape assuming a galactic-to-extragalactic source transition at the ankle of the cosmic ray spectrum. Although the allowed parameter space of this model is plausible, the average burst likely exhibits Γ and f_p values that are largely excluded for neutrino production. Nevertheless, we have only considered single-zone neutrino and cosmic ray emission models of GRBs in this work, and models of multiple emission regions predict a neutrino flux below our current sensitivity (Bustamante et al. 2015a; Globus et al. 2015).

Our acceptance of possible prompt neutrino signal can increase further through the addition of searches for ${\nu }_{\mu }$-induced tracks correlated with Southern Hemisphere GRBs that occurred during the 79-string and completed detector configurations. Moreoever, the next-generation IceCube-Gen2 detector will significantly improve the sensitivity to transient sources (Aartsen et al. 2014d; Ahlers & Halzen 2014). The continued pursuit of all neutrino flavors from observed GRBs over the entire sky will either reveal a flux that is still lower than our current sensitivity or increasingly disfavor these phenomena as sources of the highest energy cosmic rays.

We acknowledge the support from the following agencies: the U.S. National Science Foundation-Office of Polar Programs, the U.S. National Science Foundation-Physics Division, the University of Wisconsin Alumni Research Foundation, the Grid Laboratory Of Wisconsin (GLOW) grid infrastructure at the University of Wisconsin—Madison, the Open Science Grid (OSG) grid infrastructure; the U.S. Department of Energy, and the National Energy Research Scientific Computing Center, the Louisiana Optical Network Initiative (LONI) grid computing resources; the Natural Sciences and Engineering Research Council of Canada, WestGrid and Compute/Calcul Canada; the Swedish Research Council, the Swedish Polar Research Secretariat, the Swedish National Infrastructure for Computing (SNIC), and the Knut and Alice Wallenberg Foundation, Sweden; the German Ministry for Education and Research (BMBF), Deutsche Forschungsgemeinschaft (DFG), the Helmholtz Alliance for Astroparticle Physics (HAP), the Research Department of Plasmas with Complex Interactions (Bochum), Germany; the Fund for Scientific Research (FNRS-FWO), the FWO Odysseus programme, Flanders Institute to encourage scientific and technological research in industry (IWT), the Belgian Federal Science Policy Office (Belspo); University of Oxford, United Kingdom; the Marsden Fund, New Zealand; the Australian Research Council; the Japan Society for Promotion of Science (JSPS); the Swiss National Science Foundation (SNSF), Switzerland; the National Research Foundation of Korea (NRF); and the Danish National Research Foundation, Denmark (DNRF).