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

Gamma-ray bursts ( GRBs ) are considered as promising sources of ultra-high-energy cosmic rays ( UHECRs ) due to their large power output. Observing a neutrino ﬂ ux from GRBs would offer evidence that GRBs are hadronic accelerators of UHECRs. Previous IceCube analyses, which primarily focused on neutrinos arriving in temporal coincidence with the prompt gamma-rays, found no signi ﬁ cant neutrino excess. The four analyses presented in this paper extend the region of interest to 14 days before and after the prompt phase, including generic extended time windows and targeted precursor searches. GRBs were selected between 2011 May and 2018 October to align with the data set of candidate muon-neutrino events observed by IceCube. No evidence of correlation between neutrino events and GRBs was found in these analyses. Limits are set to constrain the contribution of the cosmic GRB population to the diffuse astrophysical neutrino ﬂ ux observed by IceCube. Prompt neutrino emission from GRBs is limited to  1% of the observed diffuse neutrino ﬂ ux, and emission on timescales up to 10 4 s is constrained to 24% of the total diffuse ﬂ ux.


INTRODUCTION
Gamma-ray bursts (GRBs) are short bursts of gamma radiation and are among the most energetic events in the Universe (Zhang & Mészáros 2004;Mészáros 2006). The primary burst of gamma rays, called the prompt emission, lasts for about 10 −3 s to 10 3 s. GRBs are broadly classified into two categories based on the duration of their prompt emission: short GRBs (for bursts shorter than 2 s) and long GRBs (for bursts longer than 2 s) (Mazets et al. 1981;Norris et al. 1984;Kouveliotou et al. 1993). Any particle emission observed prior to and after the prompt emission is referred to as precursor and afterglow emission, respectively. Short GRBs are generally observed to have a harder energy spectrum than long GRBs (Zhang et al. 2012). Although their exact emission mechanism is not well understood, the predominant model for GRB phenomenology includes the emission of a relativistic fireball triggered by the interaction of accreting matter onto a compact central object (Kumar & Zhang 2015;Mészáros 2006). The recent observation of gravitational wave emission from a neutron star (NS) binary merger, GW170817, in coincidence with the short GRB 170817A (Abbott et al. 2017), confirmed that short GRBs can be produced by mergers of compact objects. Long GRBs have been previously linked to the core collapse of super-massive stars by the observation of coincident supernovae (Hjorth & Bloom 2012;MacFadyen & Woosley 1999).
Central engines of GRBs drive a highly relativistic jet beamed into a narrow opening angle (Zhang & Mészáros 2004;Rhoads 1997). The jet fireball is hypothesized to be a plasma arising from a quasi-thermal equilibrium between radiation and e − e + pairs. Multiple shells of plasma can be emitted, which propagate outward from the central engine into the interstellar region with a varying Lorentz factor (Paczynski 1986;Goodman 1986). When two shells collide, a shock wave will develop that can accelerate charged particles to higher energies via first order Fermi acceleration (Krymskii 1977;Bell 1978). At a later stage, the relativistic outflow from the fireball will interact with the interstellar medium, leading to external shocks (Mészáros 2006). In some models, protons and ions are accelerated at the sites of internal and external shocks to energies in excess of 10 20 eV leading to emission of Ultra-High-Energy Cosmic Rays (UHECRs) (Waxman 1995;Vietri 1995;Murase & Nagataki 2006). UHECRs observed in coincidence with GRBs would offer direct evidence of this hadronic acceleration. However, cosmic rays get deflected by the intergalactic magnetic fields as they propagate through space. Thus they neither point back to their sources, nor reach us at the same time as the GRB gamma rays. Fortunately, neutrinos offer an alternative approach to identify the progenitors of cosmic rays. Fermi accelerated protons can interact with the gamma rays produced in the fireball and lead to the photo-meson production of pions, which can create an accompanying burst of neutrinos (Waxman & Bahcall 1997). These interactions can take place through the following channels: p + γ → n + π + ; π + → µ + + ν µ ; µ + → e + +ν µ + ν e . (1) Assuming sufficient pion production, GRBs can collectively produce a diffuse neutrino flux observable at Earth above 0.1-1 PeV (Waxman & Bahcall 1997;Gandhi et al. 1998;Winter et al. 2014;Bustamante et al. 2015;Baerwald et al. 2015;Biehl et al. 2018). Neutrinos effectively only interact through the weak force, propagating through the Universe without deflection and thus point back to their sources. Detecting high-energy neutrinos correlated with GRBs would establish them as cosmic-ray acceleration sites.
The IceCube Neutrino Observatory is currently the most sensitive instrument for the detection of astrophysical neutrinos (Ahlers et al. 2018). In 2013 the IceCube Collaboration first reported the discovery of an astrophysical neutrino flux (Aartsen et al. 2013). This was later corroborated by the discovery of a hard spectrum of muon events in the Northern Hemisphere (Aartsen et al. 2016a). While candidate neutrino sources have been identified by IceCube (Aartsen et al. 2018) and others (Stein et al. 2021), the origin of the diffuse astrophysical flux is not fully understood and may have several classes of progenitors. Because searches for neutrinos in coincidence with brief, transient phenomena are nearly background free, IceCube has pursued several different types of analyses designed to search for neutrino correlations with the prompt phase of GRB observations (Abbasi et al. 2012;Aartsen et al. 2015Aartsen et al. , 2016bAartsen et al. , 2017a, but found no associations. These results are consistent with non-detections in analyses performed by AMANDA (Achterberg et al. 2007(Achterberg et al. , 2008 and ANTARES (Albert et al. 2017a(Albert et al. , 2021a. The most recent IceCube results (Aartsen et al. 2017a;Abbasi et al. 2021a) have put constraints on the single-zone fireball models of GRB neutrino and UHECR production during the prompt phase. IceCube has also performed dedicated searches for neutrinos coincident with gravitational wave candidates, such as GW170817, detected by LIGO and Virgo (Adrián-Martínez et al. 2016;Albert et al. 2017bAlbert et al. ,c, 2019Aartsen et al. 2020). To date, no correlation has been found on 10 3 s time scales, but additional studies are ongoing (Keivani et al. 2021).
Recent observations by imaging air Cherenkov telescopes (IACTs) have shown that TeV particles can be produced during the afterglow phase, up to several days after the prompt emission (Acciari et al. 2019;Abdalla et al. 2019Abdalla et al. , 2021. Additionally, it has been shown that ∼10% of GRBs have an observed gamma-ray precursor that precedes the main burst by a few tens of seconds, but in extreme cases, up to 10 minutes (e.g. Burlon et al. (2009);Charisi et al. (2015); Coppin et al. (2020)). The complementary studies presented in this paper extend the search for neutrino correlations to precursor time windows, as well as extended precursor and afterglow time windows of up to −14 to +14 days around GRB gamma-ray triggers. The low background rates allow for highly sensitive searches on the scale of days to weeks. All the analyses use 7.16 years of IceCube muon neutrino candidate events. Section 2 describes the IceCube detector and the event selection for the neutrino dataset used in the four analyses. In Section 3, we describe the catalog of GRBs that was used in our analyses, as well as the different selection cuts on the GRB sample that were considered for the respective analyses. In Section 4, we describe the methods used for evaluating the statistical significance of results, and in Section 5 we describe the analysis approach and results for each study. Our limits and interpretation on neutrino emissions from cosmic GRB populations are presented in Section 6. Section 7 then provides the concluding remarks and outlook.

DETECTOR & EVENT SELECTION
The IceCube Neutrino Observatory is a cubic-kilometer-scale Cherenkov detector buried deep in the South Pole ice. IceCube consists of 5,160 digital optical modules (DOMs) arranged in an array of 86 strings deployed 1,450 to 2,450 meters below the ice surface (Aartsen et al. 2017b). The strings are arranged in a hexagonal grid with 125 m spacing between adjacent strings and with each string containing 60 DOMs, spaced 17 m apart along the string. Each DOM houses a downward facing photo-multiplier tube (PMT) inside a spherical transparent glass capsule. The DOMs are designed to be sensitive to the Cherenkov radiation produced by the secondary particles that result from interactions of neutrinos with the ice. Cherenkov radiation is produced when a charged particle travels faster than the speed of light in a dielectric medium, resulting in conical emission of photons along the path of the charged particle. DOMs record the waveforms of Cherenkov photons, from which the number of photo-electrons and their arrival time can be extracted (Aartsen et al. 2017b). This information is combined from all DOMs to reconstruct the Cherenkov light cones and infer the energy and direction of the particles that produced them.
IceCube is sensitive to all three flavors of (anti-)neutrinos; however, the data set used in these analyses is optimized to select charged-current interactions from muon (anti-)neutrinos, as they offer the best pointing resolution. These interactions result in the production of a muon that will propagate through the detector in a straight line, depositing light along its track. The typical angular resolution for these tracks is 1°for muons with energies 1 TeV.
The sample used for this paper is the IceCube gamma-ray follow-up (GFU) data consisting of well-reconstructed muon tracks collected from 2011-05-13 through 2018-10-14 (Aartsen et al. 2017c). The vast majority of events that trigger the IceCube detector are not astrophysical neutrinos, but muons produced in cosmic-ray air showers. In the Southern hemisphere, atmospheric muons are the dominant background and are observed at a rate of 2.7 kHz (Aartsen et al. 2017b). Since only neutrinos can propagate through the Earth without being absorbed, this background vanishes in the Northern hemisphere, where atmospheric neutrinos dominate the background. A selection with different data quality cuts for the Northern and Southern hemisphere is therefore used, which reduces these backgrounds to 6.6 mHz integrated over the full sky (Aartsen et al. 2017c). A detailed account of this event selection is given in (Aartsen et al. 2017c).

GRB CATALOG
Space-based gamma-ray observatories, such as Swift and Fermi, as well as a variety of ground-based observatories, continuously monitor the sky for high-energy gamma-ray activity (Gehrels et al. 2004;Meegan et al. 2009). These observatories provide hundreds of GRB measurements per year, which we used to construct a GRB catalogue to enable our coincidence study. GRBweb (Coppin 2021) is an IceCube project designed to combine the observational data from all major GRB observatories into a single database. Input data to GRBweb primarily originates from online GRB catalogs, such as those by Fermi (von Kienlin et al. 2020; Fermi-LAT collaboration 2019), Swift (Lien et al. 2016), and IPN (Hurley et al. 2013), and from the automated parsing of GCN circulars (Barthelmy 1998). An all-inclusive GRB catalog is thus constructed. Internally, GRBweb makes use of a set of automated Python scripts that process and save the data into a SQL database. GRBweb currently contains over 7,500 GRBs and is updated on a weekly basis.
GRBweb uses a set of predefined conditions 1 to determine which data will be used if the burst was observed by multiple detectors. For instance, the burst direction is set equal to the localization with the smallest angular uncertainty. Burst times, T 0 , correspond to the earliest time at which gamma-ray activity was reported. For the GRB duration, two variables are considered: the conventional T 90 and new, composite variable called the T 100 . The T 100 is defined as the time difference between the last and earliest reported time of gamma-ray activity. Each analysis selected a subset of GRBs between 2011 and 2018. The selection was motivated by the effect of localization uncertainties on the sensitivity of the analysis, as well as requirements on timing or precursor information. The differences between the analyses is summarized in Table 1 and explored in detail in Section 5. As the data set ends in October 2018, the only IACT-detected burst in our sample is GRB 180720B (Abdalla et al. 2019), which is discussed in Section 6.

STATISTICAL METHOD
Each of the analyses makes use of an unbinned likelihood ratio method to quantify the potential correlation between GRB observations and IceCube events, using a blind analysis technique. In this section, we present details on the likelihood ratio method used to determine the p-value of individual GRBs, as well as the binomial test used to assess the group of p-values from all GRBs.

Likelihood Ratio & Test Statistic
An unbinned likelihood ratio method (Braun et al. 2008;Achterberg et al. 2006;Aartsen et al. 2017a) is combined with frequentist statistics to assign a probability that a subset of neutrino candidate events is consistent with background. For a sample of N candidate neutrino events with characteristics x i , the likelihood can be written as where p s = n s /(n s + n b ), p b = n b /(n s + n b ), and P N is the Poisson probability to observe N events, assuming n s signal events and n b expected background events, S and B denote the probability distribution functions (PDFs) describing the spatial and energy distribution of signal events and background events, respectively. For the signal energy PDF, an E −γ power-law spectrum is assumed, where the spectral index γ is either a fit parameter or is fixed to γ = 2, depending on the specific analysis (details in Section 5). The signal space PDF, shown in Eq. (4), uses a 2D Gaussian to test the compatibility of the neutrino candidate's reconstructed position, x ν , with the source position, x GRB , where σ 2 is the quadratic sum of the uncertainty on the reconstructed neutrino and GRB direction. Since the contribution of signal events to the total data set is expected to be extremely small, the data can be used to construct the background PDFs. In particular, the background energy PDF is assumed to follow the energy distribution of data. Due to the detector geometry, the approximation of azimuthal symmetry can be used to describe the background space PDF solely as a function of zenith. Any neutrino emission that occurs within a given analyzed time window is fitted as constant emission during the time window. The specific time windows used in each analysis are described in Section 5.
The likelihood Eq. (2) is evaluated for different values of n s using the PDFs described above, withn s designating the value of n s that maximizes the likelihood. A likelihood ratio with respect to the null hypothesis L(n s = 0) is then used to obtain the following test statistic (TS): High-energy events in temporal and close spatial coincidence with a GRB will result in a large TS value. The p-value of an observation is determined by comparing the observed test statistic to a test statistic distribution created from scrambled data (see Section 4.4). Two of the four analyses use localizations provided by Fermi -GBM for some GRBs and therefore have an additional step to determine the test statistic. Starting in early 2018, the Fermi-GBM collaboration began releasing a HEALPix skymap (Górski et al. 2005) with the localization probability as a function of sky position for each GRB localized by GBM (Goldstein et al. 2020). Maps prior to 2018 were processed in a similar way using the GBM Data Tools, but the metadata in the files have not been fully qualified and the files have not yet been uploaded to the final HEASARC archive, therefore we use the preliminary files from Goldstein & Wood (2022). These maps contain a per-pixel probability that a given GRB originates from that direction. An all-sky scan of neutrino data is performed, in which TS orignial (equal to TS from Eq. (5)) is calculated at every pixel of the skymap and then penalized by the probability, P GBM , of that pixel where P GBM,max is the maximum probability on the entire skymap. The position of the GRB on the sky and the number of signal events are thus both fitted to find the combination which maximizes TS final .

Stacked Likelihood Analysis
Equation (2), which describes the likelihood for a single well-localized GRB, can be easily modified to describe the likelihood of N GRB well-localized GRBs. When considering multiple sources, n s corresponds to the total number of signal events summed over all GRBs. Each GRB is assumed to have the same time-integrated neutrino flux at Earth. The expected number of neutrinos from each GRB will therefore be proportional to the effective area at the declination of the burst, A ef f (δ), which is calculated assuming an E −2 neutrino energy spectrum. The signal PDF is then replaced by the sum of the N GRB signal PDFs, weighted by their relative contribution to the number of signal events where S j and δ j are the signal PDF and declination of the j'th burst, respectively. This method, known as a stacked likelihood analysis, provides a single measure of the total neutrino emission from a set of GRBs.

Cumulative Binomial Test
When not performing a stacked search, analyzing a selection of N GRBs will result in a p-value for each individual burst. A trial-correction method is thus needed to determine if one or more of the obtained p-values provides a statistically significant result. Arranging the p-values from smallest to largest, their values are denoted as p 1 , p 2 , ..., p N . Under the null hypothesis of no neutrino emission, the correlations of neutrino events with GRBs will only occur randomly, and these N p-values are expected to follow a uniform distribution between 0 and 1. The probability that k or more p-values are smaller than or equal to p k is thus given by the binomial probability: Evaluating Eq. (8) over all potential values of k, the smallest P (k) is selected. An empirical trial-correction factor is then applied to account for the fact that N potential p-values were scanned. This trial correction is found by determining the fraction of background-only realizations that produce a more significant result. A visualization of this In this example, k = 10 (i.e. the ten GRBs with the smallest ppre values) was found to be most significant subset, with P (k) = 0.043 (red line). Note that the accumulation of p = 1 at large k is typical for low-background searches (GRBs withns = 0 are assigned p = 1). Right: The trial correction process is shown, with the most significant P (k) from the left plot corresponding to a post-trial p-value of 0.4. The cumulative distribution in this plot is created by performing the analysis on scrambled data sets.
trial-correction procedure is shown in Figure 1. This procedure also ensures that the rare overlap of a single neutrino contributing to two analyzed GRBs (thus creating a correlation between two p-values) is also accounted for in the final post-trial p-value. On the other hand, given that the successive probabilities P (k) are strongly correlated with one another, the overall trial correction factor is modest in the end.

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

ANALYSIS
IceCube has previously reported limits on neutrino production during the prompt phase of GRB observations and found no association, yielding limits that have constrained several leading models of GRBs (Aartsen et al. 2017a). However, recent observations of gamma-ray emission outside of the prompt phase motivate a more comprehensive search. The four analyses described below allow for the possibility of neutrino emission outside of the prompt phase. As with previous GRB searches, all analyses are primarily based on the spatial and temporal correlation between neutrino events and GRB observations, while they differ in their specific assumptions about the neutrino emission. The first two analyses are the most general, with the "Extended TW" analysis opening the observation window to up to one day prior and two weeks following the prompt emission. The "Precursor/Afterglow" analysis focuses on a smaller sample of well-localized GRBs in our catalogue, while treating precursor emission separately from emission during and after the prompt phase. The final two analyses focus specifically on the precursor phase, with the "GBM Precursor" analysis examining GRBs for which Fermi-GBM detected gamma-ray emission prior to the prompt phase, The time windows used in the "Extended TW" analysis. These 10 time windows are evaluated for all GRBs, regardless of their duration. The first 9 time windows are centered on the T100, while the longest time window is asymmetric around the T100. The T100 of a given GRB is the duration between the earliest and latest measured time of gamma-ray activity. The shortest time window that completely contains the T100 is used to determine the prompt limit. In the example shown above, the prompt time window is 10 3 seconds. (b) A schematic representation of how the best-fit time windows are obtained for the precursor (left) and the prompt+afterglow searches (right) in the "Precursor/Afterglow" analysis. (c) Illustration of a gamma-ray light curve in which the GRB prompt phase is preceded by a precursor. The time range of the precursor, extended by 2 s on either side, is used to search for coincident neutrinos in the "GBM Precursor" analysis.
while the "Stacked Precursor" analysis examines well-localized bursts for which no precursor emission was observed. Together, these analyses provide both a comprehensive and a model-independent approach to the search for neutrino production in GRBs. Each analysis and its results are described in more detail below. Limits and interpretation of the results are then presented in Section 6.

Extended TW
This analysis searches for neutrinos in a range of time windows (TW) and uses the largest GRB catalog of the four analyses. All GRBs with known duration were included, given they were observed during the detector livetime. 163 GRBs do not have a reported T 100 in GRBweb at the time of writing, and so these GRBs were excluded from the selection. The final number of GRBs in this analysis is 2,091.
The p-values are calculated for these 2,091 GRBs in 10 pre-determined time windows and with a fixed E −2 neutrino energy spectrum in the signal hypothesis. The first 9 time windows range from 10 seconds to 2 days, centered on the T 100 of the GRB, and the final time window is asymmetric with a 1 day precursor and 14 day afterglow time window (see Figure 2a). The shortest pre-defined time window that fully envelopes the T 100 interval will be used to describe the prompt emission. The seven shortest time windows range up to 10 3 seconds, which includes most GRB prompt phases. The 10 second time window is sufficiently small to study the short GRBs, because the neutrino background is effectively zero at this timescale. For the longer time windows the background becomes non-negligible, which is what motivates the decision to search up to ∼ two weeks but not longer. For a well-localized GRB (σ 1 • ), the 15-day time window has an expected background on the order of one neutrino candidate event. For poorly-localized GRBs this longer window is less sensitive, which is why the emphasis was placed on the afterglow region where higher-energy neutrinos are predicted (Asano & Murase 2015). The test statistic in each time window is calculated using Eq. (6). From those 10 p-values, the most significant one is selected to represent the GRB.
Since each GRB is studied in 10 time windows, a correction is required to compensate for the look-elsewhere effect. Because the time windows are correlated, an effective trial correction is used. For every given GRB, the smallest p-value of the 10 time windows is selected. Background scrambles are then used to determine the probability of obtaining a Table 2. Post-trial p-value of the binomial test for the "Extended TW" study of 2,091 GRBs. The binomial test was run on four subsets of GRBs split by hemisphere and prompt gamma-ray duration. The number of GRBs in each sub-population is indicated in parentheses.
Northern Long (960) Northern Short (183) Southern Long (814) Southern Short (134) 0.038 0.799 0.898 0.849 Table 3. The top GRB result found in each sub-population in the "Extended TW" study of 2,091 GRBs. The number of GRBs in the sub-population is indicated in parentheses. The column titled ppre gives the p-value for the GRB, without a correction for the size of the sub-population. This ppre p-value has been corrected for searching 10 time windows. The column ppost provides the corresponding p-value including the look-elsewhere correction for both the population size and for searching ten time windows. value that is smaller than or equal to this result. Next, the p-values that have been corrected for searching 10 time windows are evaluated in a binomial test (Section 4.3) to search for evidence of a sub-set of GRBs with significant neutrino emission. Four binomial tests are evaluated, where the total GRB sample has been divided into four sub-populations by duration and hemisphere. Following von Kienlin et al. (2021), GRBs are separated into short and long classes according to whether the measured T90 is less than or greater than 2 s, respectively. This is intended to account for the different progenitor classes of merger events (Abbott et al. 2017) and core collapse supernovae (Hjorth & Bloom 2012), albeit in a simplistic way since the two classes are known to overlap. This overlap explains why GRB 170817A falls into the long GRB class with a T90 of 2.048 s despite being a known short GRB from a binary neutron star merger. In this case, its misclassification does not notably affect the results as GRB 170817A would not have appeared in the top 5 short GRBs from the Southern Hemisphere (listed in Appendix A.2) even if it were included in that class. Future studies will explore grouping methods which account for the fact that GRB 170817A is a short GRB. The split by hemisphere allows for the increased sensitivity of the analysis to GRBs in the Northern sky.
The final p-value for each binomial test is determined by comparing the result found for data with the results found for scrambled data sets. The final p-values are summarized in Table 2 and are consistent with background. The most significant GRB (pre-trial) from each sub-population is listed in Table 3. All GRBs with a p-value less than 1% are listed in Appendix A.2.

Precursor/Afterglow
In the "Precursor/Afterglow" analysis, the time window for the signal region is treated as a parameter which can be estimated from the data. The time window can be fitted from a minimum duration of 0.5 s up to a maximum of 2 weeks, long enough to cover a wide range of possible neutrino emission time scales. Only GRBs with a positional angular uncertainty of less than 0.2°are considered in this analysis (thus they can be approximated as point sources), which keeps the background low within the full time window range. GRBs within the first 14 days and the last 14 days in the GFU data are further removed for this analysis. This results in a selection of 733 GRBs. This includes 53 GRBs that did not have a reported T 100 in GRBweb and were not included in the Extended TW analysis. Two separate searches are performed for every GRB, denoted as the "precursor" search and "prompt+afterglow" search, each with fit parameters corresponding to the number of signal events (n s ), the spectral index (γ), and the width of the emission time window (T w ). The only difference between the two searches is how the time window T w is defined with respect to the start of the prompt phase T 0 (see Figure 2b): • For prompt+afterglow searches, the time window extends forwards from T 0 to the time (T 0 + T w ) after the start of the prompt phase.
For every GRB, the analysis returns the best-fit parameters for the respective search and the p-value. This results in two lists of 733 p-values, one list for each search. A binomial test (Section 4.3) is performed on each list to search for a subset of GRBs with significant neutrino emission. The final p-value after each binomial test is determined by comparing the result found for data with the results found for scrambled data sets.

GBM Precursor
The two analyses above searched for neutrino correlations with GRBs in an agnostic manner. It is possible to use additional gamma-ray observations to perform more sensitive dedicated analyses. A previous study by Coppin et al. (2020) analyzed the light curves of all Fermi -GBM bursts up to 2020, leading to a sample of 217 GRBs that exhibit signs of gamma-ray precursors. Of those 217 bursts, there are 133 GRBs that overlap with the IceCube data set examined here. A dedicated search is performed to investigate potential neutrino production coincident with the emission phase of the observed gamma-ray precursors using these 133 GRBs. The time windows of the analysis are set to those of the identified gamma-ray precursors, extended by 2 s on either side to obtain a restrictive yet conservative range. Summed over all 133 GRBs, a total time window of 3.3 · 10 3 s is examined. The signal hypothesis in the likelihood uses a fixed E −2 neutrino energy spectrum. Out of 133 bursts, 100 GRBs were localized solely by Fermi -GBM. For those 100 bursts, a TS is defined following Eq. (6). Otherwise, the TS is based on Eq. (5).
Performing the analysis on all GRBs results in 133 individual p-values. To trial correct the result, a procedure similar to that of the binomial correction is applied. However, instead of considering the kth most significant p-value, the product of the k most significant p-values is considered. Comparing this result to the distribution of the same statistic (product of the k-most significant p-values) found in scrambled data sets yields the final p-value for the analysis.
As this analysis only targets GRBs with observed gamma-ray precursors, it has the smallest background of all four analyses. The neutrino flux to which the analysis is sensitive (Sec. 4.4) corresponds to only 1.7 neutrino events on average. Considering events within a relative combined neutrino and GRB angular uncertainty of 5σ, no neutrino events were observed in temporal coincidence with any of the 133 precursors. A final p-value of 1 was thus obtained.

Stacked Precursor
One of the results by Coppin et al. (2020) was that almost all (>95%) gamma-ray precursors were found to occur within 250 s of the prompt emission. A fourth analysis was therefore performed to enable a larger systematic search for precursor neutrinos, not limited to Fermi -GBM bursts for which a gamma-ray precursor could be resolved. For all well-localized 2 GRBs within the IceCube GFU data period, corresponding to 872 bursts, a generic time window of 250 s prior to the T 0 time was searched for excess neutrinos. All 33 well-localized GRBs from the GBM Precursor analysis are included in this search. For the set of 872 GRBs, a single stacked test statistic is constructed according to Eq. (5) and assuming an E −2 signal spectrum. Due to the stacking procedure, this analysis automatically leads to a single p-value, not requiring a trial correction procedure.
When the stacked precursor analysis is applied to these 872 GRBs, no excess of neutrino events is found (n s = 0). Overall, five low-energy events within the given time window are in loose spatial coincidence with the examined GRBs, listed in Tab. 4. Given the length of the time window and the number of bursts, these coincident events are fully consistent with the background expectation. Similar to the GBM precursor analysis, a final p-value of 1 is thus obtained. Since none of the analyses found evidence for neutrino emission from GRBs, limits can be placed on the GRB neutrino flux. These limits can be compared to model predictions and can constrain the fraction that GRBs contribute to the diffuse neutrino flux measured by IceCube (Stettner 2019). It is worth noting that upper limits on the diffuse flux presented in previous IceCube publications (Abbasi et al. 2012;Aartsen et al. 2015Aartsen et al. , 2016bAartsen et al. , 2017a correspond to the flux from GRBs observable by current gamma-ray satellites. In contrast, for the results presented here the limits are placed on the total contribution to the diffuse flux from all GRBs in the Universe. For the "Extended TW" and "Stacked precursor" analyses, the upper limits on the time integrated flux from all analysed GRBs, dN dEdA , are converted to a diffuse flux where the factor ∆Ω normalises the flux by the solid angle subtended by the GRBs, i.e. ∆Ω = 4π for an all-sky sample, ∆t ∼ 7.16 yr is the livetime of the IceCube data used in this work, d corrects for the field-of-view (FOV) and dead-time of the GRB telescopes, and z accounts for the contribution from GRBs that are too dim to be observable with current gamma-ray satellites. Determining the value of d and z requires specifying concrete characteristics of the GRB satellites. In the "Extended TW" analysis, which uses GRBs detected by a variety of satellites, the canonical estimate is made that with no dead-time and a all-sky FOV, a total of 667 GRBs would be observed every year (Aartsen et al. 2017a). In the "Stacked precursor" analysis, only well-localized GRBs are used, most of which were detected by Swift. Limits are therefore set on the time-integrated flux of the subset of Swift bursts, and then corrected using Eq. (9) based on the characteristics of the Swift telescope. In the "Precursor/Afterglow" analysis, instead of finding a limit on the stacked flux and using it in Eq. (9), a simulation of the neutrinos from individual GRBs following the cosmological distribution is used to provide simulated data sets for the analysis, allowing to set limits directly on the total flux by varying the injected signal strength from the population.
6.1. Extended TW 2 A less restrictive cut of 1.5 • is applied on the localization compared to the "Precursor/Afterglow" analysis, as the significantly smaller time window and reduced number of fit parameters lead to a smaller effective background. The stacked test statistic presented in Section 4.1 is used to place an upper limit on contributions of GRBs to the quasi-diffuse flux measured by IceCube. Flux is injected using an E −2.28 (Stettner 2019) spectrum until 90% of trials yield a stacked test statistic above the unblinded value. This injected flux is converted to a diffuse flux using Eq. (9). This procedure is repeated for all ten time windows and the prompt. For the prompt, the shortest time window that includes the entire reported T 100 is used (see Figure 2). These 90% confidence level limits are set for each time window for various subsets of GRBs. The four sub-populations analyzed with a binomial test are presented in Figure 4. Limits are also placed on all GRBs observed in the Northern and Southern Sky, as well as all short and long GRBs ( Figure 5). In each plot, the stacked limit from all 2,091 GRBs is shown for reference.
Previous IceCube studies (Aartsen et al. 2017a) constrained the prompt contribution of GRBs observable by current gamma-ray satellites to ∼1% of the diffuse flux observed by IceCube. The prompt limit presented here applies to all GRBs in the Universe and corresponds to 1%. The limits are similar despite analyzing nearly twice as many GRBs in this analysis. The difference is the inclusion of z term to account for GRBs that are too far away to be observable with current gamma-ray satellites. Given its fluence trigger threshold of ∼ 10 −8 erg cm −2 s −1 , the Swift-BAT detector can be assumed to view all canonical GRBs with an isotropic equivalent luminosity L iso ≥ 10 50 erg cm −2 up to a redshift of z ∼ 1.3. Assuming that all GRBs have identical properties in terms of neutrino emission and that they follow the redshift evolution described by Lien et al. (2014), the contribution from GRBs outside the observable redshift threshold can be calculated using the procedure outlined in (Ahlers & Halzen 2014). This results in z value of ∼ 2.1. The inclusion of this z term leads to a similarly constraining prompt limit compared to previous studies.

Precursor/Afterglow
These results were used to calculate limits for the total contribution of all long GRBs to the observed diffuse neutrino flux (see Figure 6). Based on the study by Lien et al. (2014) using Swift observations of GRBs to extrapolate to the whole Universe, a mean rate of 4,571 long GRBs per year is estimated to be beamed at Earth. We use the software package FIRESONG (Tung et al. 2021) to simulate neutrino emission from this cosmic GRB population, with different emission periods. Emission windows ranging from 100 s to 14 days were considered in the simulations. The redshift distribution of the cosmic GRB population from the study by Lien et al. (2014) was used to simulate the GRB rate density in FIRESONG. We assume the case where no luminosity evolution is required and for simplicity assume that the GRBs are standard candle neutrino sources.
The Precursor/Afterglow analysis only considers 733 GRBs across 7.16 years, out of which 546 GRBs are long GRBs which were observed by Swift. Hence we use the GRB detection efficiency as a function of redshift from the study by Lien et al. (2014) as well as the sky coverage and the survey time of Swift/BAT to down-select from the 4,571 Figure 5. Left: The stacked limit on the quasi-diffuse flux (at 1 TeV) for all Northern GRBs (blue) and all Southern GRBs (red). Right: The stacked limit on the quasi-diffuse flux for all long GRBs (dark blue) and the long GRBs split by northern and southern sky (red and light blue). In both figures, the right axis presents this limit as a fraction of the quasi-diffuse flux (Stettner 2019). The limits for all 2,091 GRBs are shown in green for reference. Each dot indicates the 90% confidence limit for the given time window, and the dashed lines show the limit for the prompt. sources created by the FIRESONG simulations to a sample of 546 GRBs to recreate the observation biases of the Swift sample. These 546 flux values down-selected from the simulation are used to inject signal with an E −2.28 (Stettner 2019) spectrum into a simulated neutrino dataset. To repeat the original analysis under the same conditions, an additional 187 sources are added with no signal (i.e. represent background only), bringing the sample to 733 sources again. The likelihood analysis and binomial test are then performed as before on this simulated sample. This sequence of steps for every emission window considered is repeated using different simulated neutrino datasets to produce 1000 trials. When GRBs are assumed to produce the entire diffuse flux (Stettner 2019), 100% of these injected trials result in a binomial test result more significant than the unblinded result. The fraction of the diffuse flux is then reduced to identify the flux where 90% of trials produce a binomial test result more significant than the unblinded result. This is summarized in Figure 6, where the points show this upper limit (90% confidence level) for a range of neutrino emission time windows. The respective fraction for the different emission durations considered represents the total allowed contribution of all long GRBs to the observed astrophysical diffuse neutrino flux.

Stacked Precursor
Precursor neutrinos have been predicted in models in which a jet initially has to burrow its way through remnant layers of the progenitor star. A prediction for the diffuse precursor neutrino flux from such sources was made by (Razzaque et al. 2003) and is shown in Figure 7. The red and green lines correspond to progenitors that have a remnant outer hydrogen (H) or helium (He) shell, respectively. This analysis is able to exclude the H-shell model by He shell model prediction H shell model prediction H shell model limit HESE best fit (7.5 yr) ν µ +ν µ best fit (9.5 yr) Analysis limit (E −2 ) Analysis limit (differential) Figure 7. Model predictions by (Razzaque et al. 2003) are shown for two progenitor scenarios. The red and green line corresponds to central engines enveloped by an outer hydrogen and helium shell, respectively. The "Stacked precursor" analysis can exclude the H-shell model by a factor 10, but does not constrain the He-shell model. Aside from these model comparisons, generic flux upper limits are also shown. Limits for an E −2 spectrum are displayed by the black line, where the x-range corresponds to the energy band that contributes 90% of all signal events. The solid grey line shows the differential flux, normalized per decade of energy, to which the analysis is sensitive. These generic upper limits extrapolate the total GRB flux in a different way than the model prediction/limits, as explained in subsection 6.3. For reference, the astrophysical neutrino flux observed by IceCube is also shown (Abbasi et al. 2021b;Stettner 2019).
a factor 10, but cannot constrain the He-shell model. To be consistent with the model prediction, the H-shell upper limit shown in Figure 7 assumes that the diffuse GRB neutrino flux results from 10 3 GRBs per year (Razzaque et al. 2003). This is in contrast to the model-independent upper limits in Figure 7, which rely on the conversion outlined in Eq. (9). This latter approach incorporates updated information about the redshift distribution of GRBs that was unavailable when the model was released. As a result, these generic limits are slightly more conservative.

CONCLUSION
The results from the four analyses presented in this paper cover 2209 unique GRBs. These GRBs are investigated for neutrino correlations from the precursor, prompt, and afterglow emission regions in a comprehensive manner and all the four analyses report observations consistent with background expectations. We obtain a range of upper limits to the astrophysical diffuse neutrino flux. We show that prompt emission from all GRBs in the Universe is limited to 1% of the diffuse astrophysical neutrino flux. Neutrino emission limits range from 1%-2% in timescales up to 10 3 s using the historic assumption of 667 GRBs observable by satellites per year. Neutrino emission is constrained to less than 24% for timescales up to 10 4 s by simulating GRB populations using the FIRESONG (Tung et al. 2021) module. These constraints are shown for additional neutrino emission timescales in Figures 4, 5 and 6. By looking for neutrinos on time scales motivated by the observation of gamma-ray precursors, we were able to constrain physical models such as those presented by Razzaque et al. (2003) (see Figure 7). Table 5 highlights the result for the GRB 180720B which was observed by H.E.S.S. to have a high-energy afterglow counterpart 10 hours after the start of prompt emission (Abdalla et al. 2019). The exceptionally bright GRB 130427A (Ackermann et al. 2014) was also analyzed, but none of the four analyses reported neutrino excesses in correlation with it. The results for the most significant GRBs from the individual analyses are presented in Appendix A and B. A comprehensive summary of all the results from each analysis can be accessed from the supplementary materials. Table 5. Summary of GRB information and best-fit results for the GRB 180720B. The positional coordinates are defined using the equatorial coordinates: right ascension (α) and declination (δ). These two quantities are expressed in degrees. T100 represents the total period of observation of prompt γ fluence from the GRB and is expressed in seconds.ns represents the best-fit number of signal events,γ represents best-fit spectral index andTw is the best-fit time window (expressed in days). E 2 F is the 90% confidence level upper limit on the time integrated flux normalization at 1 TeV in GeV/cm 2 . For Precursor results, the limit is evaluated for the best fit values ofγ andTw, and for the Afterglow results γ=2.0 and Tw=12.1 hrs is used to match with the H.E.S.S observation period (Abdalla et al. 2019). For GRB 130427A, γ=2.0 was considered for all the E 2 F calculations. For Extended TW, the limit is evaluated for the prompt phase, for Precursor/Afterglow Tw=14.0 days is used.

GRB information
Extended TW Precursor Despite non-detections in these analyses, GRBs cannot be fully eliminated as potential sources of high-energy neutrinos. A larger class of neutrino detector, such as the proposed IceCube-Gen2 detector, will offer increased sensitivity to observe potential neutrino-emitting sub-populations of GRBs. APPENDIX A. EXTENDED TW: MOST SIGNIFICANT GRBS This appendix provides details for the most significant GRBs in the Extended TW analysis. The threshold for inclusion in these tables is a p-value below 1%. For the case of short GRBs detected in the southern hemisphere, a threshold of 5% was chosen to include five GRBs in Table A.2. This was done to match the results of the binomial test, which indicated the most significant subset of the southern sky short GRBs was k=5. The other sub-populations had a much larger value of k and would have increased the p-value threshold beyond a reasonable cutoff. Although the results of the binomial test were consistent with background, these additional GRBs are included in Table A.2 for completeness.

A.1. Variables
The GRB name is based on GCN notices, with a preference for the name provided by Fermi-GBM. If the GRB is not observed by Fermi-GBM, then the name is based on the date of observation with the format YYMMDD. The final letter indicates the order of detection if more than one GRB is observed in a single day. The right ascension (α) and declination (δ) in J2000.0 equatorial coordinates, as well as the localization uncertainty (σ) are all listed in degrees. T 0 indicates the MJD trigger time of the GRB. The T 100 is provided in seconds.
The most significant time window selected by this analysis is indicated by "TW" and is given in either seconds or days with units provided in the column. The test statistic (TS) and best fit number of signal events (n s ) are provided as well as the p-value in the final column. The p-value has an effective trial correction for searching 10 time windows, but it is not corrected for the size of the given sub-population.

A.2. Results by Sub-Population
The GRBs are split into sub-populations by duration and hemisphere. Short GRBs are defined as T 90 < 2 s and long GRBs are T 90 ≥ 2 s. In this analysis, the northern sky is defined as a declination greater than −5 • , while the southern sky includes all declinations less than −5 • . The most significant GRBs for each sub-population are shown in the following tables. The information is split into two sections: GRB information based on GCN notices and Fit results based on the analysis.    Table 10. Summary of best fit results and information for the GRBs in the top 20 results for the Precursor search. Each GRB is named based on the date when it was observed with the standard formatting of YYMMDD and a letter denoting the order in which the bursts were detected on the same day (A, B etc.). The positional coordinates for the bursts are defined using the equatorial coordinates: right ascension (α) and declination (δ). These two quantities are expressed in degrees. The burst timing (T0) is expressed using the dating convention of Modified Julian Date (MJD). z represents the redshift. T100 represents the total period of the GRB observation and is expressed in seconds.ns represents the best-fit number of signal events ns,γ represents best-fit spectral index andTw is the best-fit time window (expressed in seconds). 'TS' represents the test statistic and 'Significance' represents the significance of the pre-trial local p-value for one-sided Gaussian distributions.  -Table 11. Summary of best fit results and information for the GRBs in the top 20 results for the Prompt+Afterglow search. This is similar to Table 10 but shows results from the prompt+afterglow search instead.