© The Authors 2023

1. Introduction

In the Universe, γ-ray bursts (GRBs) are among the brightest electromagnetic transient events. They consist of second- or minute-scale bursts in the soft gamma-ray band, known as the prompt phase, followed by a broadband afterglow emission in the radio, optical, and X-rays, lasting a few hours to several days. The observed photon spectra of the prompt GRB emission are typically fit by two power laws:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}{E}_{\gamma }}\propto \{\begin{array}{cc}{\left(\frac{E_{\gamma }}{E_{\gamma ,\mathrm{p}eak}}\right)}^{-\alpha }\hfill & E_{\gamma }<E_{\gamma ,\mathrm{p}eak}\hfill \\ {\left(\frac{E_{\gamma }}{E_{\gamma ,\mathrm{peak}}}\right)}^{-\beta }\hfill & E_{\gamma }>E_{\gamma ,\mathrm{peak}}\hfill \end{array}.\end{array}$$(1)

In the double power-law model, the break energy of the spectrum coincides with the peak energy E_γ,peak of the $E_{\gamma }^2\text{d}{N}_{\gamma }/\text{d}{E}_{\gamma }$ spectrum. Typical spectral indices are α ≃ 0 − 1 and β ≃ 2 − 3.

The very fact that the best fit to GRB spectra is a power-law function above E_γ,peak supports the idea that the GRB spectra are produced by nonthermal phenomena. Nonetheless, twenty years of detailed investigation of the spectra of GRBs did not result in the clear identification of the dominant radiative processes responsible for the prompt emission production. Therefore, a clear understanding of the nature of the GRB jet composition and dissipation processes therein is lacking. However, it is well established, that shocks or reconnection regions produce nonthermal populations of protons and electrons (for a review, see Zhang 2018).

Given the nonthermal nature of GRB spectra and their cosmic rate, GRBs are candidate accelerators of ultra-high energy cosmic rays (UHECRs) and are considered promising sources of high-energy neutrinos, as initially proposed in Waxman & Bahcall (1997). As a matter of fact, GRBs can provide the sufficient power to justify the spectrum of UHECRs above 10¹⁷ eV, as they are the brightest explosions in the Universe with a rate on the order of approximately one per day, as observed by the Burst And Transient Source Experiment (BATSE). Such explosions can release energies as high as ∼10⁵⁴ erg over time intervals of a few seconds (short GRBs or sGRBs) and around 100 s up to about 1000 s (long GRBs or lGRBs).

The bimodal distribution of the BATSE GRBs hints at two subclasses (Kouveliotou et al. 1993), which might be generated by two different progenitors, possibly binary mergers for the sGRBs (Eichler et al. 1989) and collapsars for lGRBs (Woosley & Bloom 2006). This separation-progenitor dichotomy might not be absolute (see Bromberg et al. 2012), as the recent detection of a lGRB in coincidence with an identified kilonova light curve proves (Rastinejad et al. 2022), followed by evidence of the Fermi Large Area Telescope (Fermi-LAT) observation (Mei et al. 2022). Recent reviews on GRBs and on models of electromagnetic and neutrino emission can be found for instance in Pitik et al. (2021) and Kimura (2022).

The first X-ray afterglow emission was observed by the Beppo-SAX satellite in 1997 (Costa et al. 1997). Before the launch of The Neil Gehrels Swift Observatory (Swift, hereafter) in 2004, the observational data could be accommodated for with an external shock model, in which the afterglow emission was induced by the interactions of the forward and reverse shock of the GRB ejecta with the circum-burst medium (Meszaros & Rees 1997). Nevertheless, the early Swift era revealed that a more complex scenario must be considered to explain the unexpected features in the X-ray light curve, such as flares and plateaus.

X-ray flares are rapid rises and falls of the intensity of the X-ray emission observed in nearly 50% of the GRBs detected by Swift (Gehrels et al. 2009). They occur tens to hundreds of seconds after the prompt trigger, and last 10 s to 10³ s. X-ray flares cannot be produced by external shocks, and they show similarities with the prompt emission (Margutti et al. 2010), which might indicate a common origin powered by a late, long-lasting activity of the central engine (Chincarini et al. 2007). In this regard, an internal-shock model is proposed at the origin of X-ray flares, in which the production is due to internal dissipation mechanisms (Fan & Wei 2005). Given the similarities with the prompt phase and the longer duration, such X-ray flares can lead to a more efficient production of high-energy neutrinos than during the prompt emission, provided that the internal dissipation mechanisms involve a population of nonthermal protons (Murase & Nagataki 2006; Kimura et al. 2017; Kimura 2022). The high-energy neutrino production is, however, suppressed if the jet is highly magnetized.

Plateaus are a shallow decay phase of the X-ray emission with time, ∝t^−a, with a typical temporal slope¹a < 0.7. Plateaus usually occur 10² s to 10⁴ s after the prompt trigger, and they are observed in roughly 60% of the Swift GRBs (Gehrels et al. 2009). The plateau can be followed by a normal (slope 1 ≲ a ≲ 2) or a steep (slope 3 ≲ a ≲ 10) decay phase in the light curve, with no spectral changes across the break. If followed by a normal decay phase, the plateau can be interpreted within the external shock model (Zhang et al. 2006). In this case, the flattening of the X-ray light curve can be attributed to a continuous energy injection supplied by a late-time activity of the central engine (Dai & Lu 1998) or to a structured jet (Oganesyan et al. 2020; Beniamini et al. 2020). The plateau followed by a steep decay phase cannot originate from external shocks (Troja et al. 2007). In this regard, many authors believe that the plateau phase is induced by internal dissipation mechanisms powered by a magnetar wind, and that the subsequent steep decay phase is associated with the abrupt collapse of the magnetar to a black hole (Lyons et al. 2010).

Neutrino telescopes have conducted many GRB searches. Prominently, IceCube has previously performed searches for spatial and temporal correlation of astrophysical neutrinos and GRBs, finding no significant association. The first exclusion was provided in Abbasi et al. (2012), in a search for neutrino emission in time windows up to ±1 d from GRBs observed from April 2008 to May 2010. In 2015, the analysis of the prompt phase of 506 GRBs in the Northern Sky from April 2008 to May 2012 yielded constraints on parameters of the single-zone fireball model, and estimated the contribution of GRB neutrinos to ∼1% of the diffuse astrophysical flux (Aartsen et al. 2015). The analysis was extended in 2017 to include three more years of IceCube data up to May 2015, and prompt GRBs observed in the whole sky (Aartsen et al. 2017a). An overall catalog of 1172 prompt GRBs was selected for this analysis. More recently, searches for generic extended time windows up to 14 d after the prompt phase and for targeted precursors have been presented in Abbasi et al. (2022). These searches, which covered more than 2000 GRBs observed from May 2011 to October 2018, place upper limits on the afterglow emission and on the contribution of GRB neutrinos to the diffuse flux, and constrain precursor GRB models. While all the aforementioned analyses adopted an IceCube sample optimized for muon selection, more suitable for point-like source searches (see Sect. 4), in Aartsen et al. (2016a) IceCube extended the search to include all-flavor high-energy neutrinos detected in three years, from May 2010 to May 2013. A total of 807 GRBs over the entire sky were analyzed.

Despite many in-depth investigations of GRB neutrinos, the emission from X-ray plateaus and/or flare afterglows has been constrained by the IceCube collaboration only as a generic afterglow emission. In the work described in this paper, we analyze the publicly available IceCube point-source data (Abbasi et al. 2021a), searching for temporal and spatial coincidence of high-energy neutrinos with the prompt phase and with the plateaus and flares of the X-ray afterglow phase to constrain a model inspired from Murase (2007), Zhang & Kumar (2013), Kumar & Zhang (2015), Kimura (2022), Guarini et al. (2022). We describe this model for clarity in Sect. 3. Along with a search from individual GRBs (referred to as “single-source search” in the following), we also perform a search for a cumulative neutrino emission (“stacking search”).

With respect to other IceCube stacking searches, our stacking method adopts physically motivated weights for the contribution of each GRB to the overall emission, thus avoiding the assumption of an equal fluence at Earth. For instance, Aartsen et al. (2017a) provides upper limits on the cumulative neutrino flux from GRBs and constraints on the parameter space of the fireball model, with the underlying assumption that the parameters determining the neutrino emission (e.g., the bulk Lorentz factor Γ and luminosity) have the same values for all GRBs. Nonetheless, it is nowadays well established that such parameters can span orders of magnitude for different GRBs. In this paper, motivated by the progress of the last decade in the understanding of the GRB phenomena, and given the increased GRB observations, we exploit the luminosity–Lorentz factor correlation (Ghirlanda et al. 2012; Lü et al. 2012) and the luminosity–peak energy correlation (Yonetoku et al. 2004) observed during the prompt phase to determine these parameters on a per-GRB case (see also Sect. 3.1). This approach improves the physical reliability of the results presented in this work compared to past searches.

The paper is structured as follows. In Sect. 2, we introduce the GRB catalogs analyzed by this work. Then, we discuss the model and we derive the stacking weights for the different GRB catalogs in Sect. 3. We describe the IceCube neutrino sample in Sect. 4. We provide details about the statistical method adopted for the single-source and the stacking searches in Sect. 5. We report the results in Sect. 6, and we discuss the related limits and physical implications in Sect. 7. Final remarks are then summarized in Sect. 8.

2. Definition of the GRB catalogs

This work analyzes three GRB catalogs. Two catalogs are based on the X-ray afterglow emission, and include GRBs with X-ray plateaus or X-ray flares, hence referred to as “plateau catalog” and “flare catalog”, respectively. A third catalog is based on the prompt phase of selected GRBs, and thus referred to as “prompt catalog”. We additionally require that the selected GRBs occur during the IceCube uptime provided with the data release (Abbasi et al. 2021a). Furthermore, the GRBs with declination δ above 80° in absolute value are excluded, as the technique used to estimate the background of IceCube neutrino events (see Sect. 5) is not efficient around the equatorial poles.

The plateau catalog comprises 260 GRBs detected by Swift/XRT, whose X-ray light curve, reconstructed by the automatic analysis of the Swift/XRT products (Evans et al. 2009), presents a power-law attenuation in time ∝t^−a, with index a < 0.7. This is the shallowest expected temporal attenuation from synchrotron radiation of external shock (Zhang et al. 2006). The flare catalog collects 200 GRBs detected by Swift/XRT, for which the online light curve analysis identifies X-ray flaring activity in the afterglow emission. A total of 85 GRBs have both X-ray plateaus and X-ray flares. They are in common to both catalogs and analyzed considering the corresponding time intervals. Examples of light curves with an X-ray plateau and an X-ray flare are provided in Fig. 1. The uncertainty on the sky localization of the GRBs that compose these two X-ray-based catalogs benefits from the very good accuracy that characterizes Swift/XRT, and that is generally smaller than the typical angular uncertainty of the IceCube track-like events (see Sect. 4).

Fig. 1. Examples of a prompt emission light curve (a) observed by Fermi/GBM, a flare (b) and a plateau (c) of the X-ray afterglow observed by Swift/XRT.

The 1751 GRBs of the prompt catalog are selected from GRBweb (Coppin 2022), an IceCube project that collects some relevant GRB parameters observed by different telescopes (Swift/XRT Lien et al. 2016, Fermi/GBM Hurley et al. 2013; von Kienlin et al. 2020, and others Ajello et al. 2019) into a single, online database. In general, Fermi/GBM is characterized by a poor angular resolution, that can be as high as ∼10° and that can exceed the typical resolutions of the track-like neutrino events. In these cases, the GRB is treated as an extended source and its spatial size is taken into account in the analysis method (see Sect. 5 for more details). Nearly 90% of the prompt catalog is composed of lGRBs.

The plateau and flare catalogs are analyzed both with the single-source and the stacking searches, described in Sect. 5, while the prompt catalog is only analyzed with the stacking search. The latter choice is motivated by the fact that IceCube has already performed extensive investigations of GRB neutrinos from the prompt phase (Abbasi et al. 2012, 2022; Aartsen et al. 2015, 2016a, 2017a). While our single-source search does not add anything to those searches, we think that a significant difference stems in using the weights defined in Sect. 3 for the stacking search. For each catalog, the stacking search is also performed on the subsamples of GRBs with measured redshift, that will be referred to as “subcatalogs”. Despite their smaller size, these subcatalogs allow us to reduce the number of assumptions on the GRBs, thus producing less strong, but more reliable results. Table 1 provides the number of GRBs in each (sub)catalog and in each hemisphere. It also anticipates the stacking weights (derived in Sect. 3), defined in terms of the prompt and X-ray isotropic equivalent fluence (S_iso and $S_{\mathrm{iso}}^{\mathrm{X}}$) and luminosity (L_iso and $L_{\mathrm{iso}}^{\mathrm{X}}$), as well as the peak energy E_γ,peak and the redshift z. Finally, it provides information about the parameters constrained in each catalog, namely the baryon loading factor ξ_p, the bulk Lorentz factor Γ, the variability timescale of the prompt phase δt_obs, or the magnetic field B (see also Sect. 3).

The catalogs are separated into the Northern (δ ≥ −5°) and Southern (δ < −5°) Hemisphere due to the different sensitivity of IceCube in these two regions of the sky. This is due to the different nature of the background, mostly downgoing muons in the Southern Hemisphere, and atmospheric neutrinos in the Northern Hemisphere.

3. The model

Protons accelerated to high energies can undergo photo-pion production with ambient photons:

$$\begin{array}{c}\hfill p+\gamma \to {\mathrm{\Delta }}^{+}\to \{\begin{array}{c}p+{\pi }^0\mathrm{BR}=2/3(\mathrm{at}\mathrm{resonance})\hfill \\ n+{\pi }^{+}\mathrm{BR}=1/3(\mathrm{at}\mathrm{resonance})\hfill \end{array}.\end{array}$$(2)

The branching ratios (BR) of the two channels are shown at the resonant energy for the production of an on-shell Δ⁺ baryon. Nonetheless, out of resonance the two channels are about equally probable, with BR ∼ 1/2 in each. While π⁰’s decay into photon pairs, the decay chain of π⁺’s leads to a copious production of high-energy GRB neutrinos, as also illustrated in Fig. 2:

$$\begin{array}{cc}& {\pi }^{+}\to {\mu }^{+}+{\nu }_{\mu }\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\mu }^{+}\to {\mathrm{e}}^{+}+{\nu }_{\mathrm{e}}+{\overline{\nu }}_{\mu }.\hfill \end{array}$$(4)

Fig. 2. Illustration of the photo-pion processes that lead to the production of high-energy neutrinos inside the dissipation region of GRB jets.

In these photo-pion interactions, each neutrino carries on average 25% of the parent pion energy. The expected ratio of each neutrino+antineutrino flavor at the source is (ν_e : ν_μ : ν_τ)_S = (1 : 2 : 0). Nevertheless, flavor oscillations of TeV–PeV neutrinos over cosmological baselines alter the proportions, and one expects on average the same flux of all flavors at Earth, (ν_e : ν_μ : ν_τ)_E = (1 : 1 : 1). Given the chain of reactions above, the emitted neutrino fluence is linked to the proton spectrum $\mathrm{d}{N}_{\mathrm{p}}/\mathrm{d}{E}_{\mathrm{p}}\propto E_{\mathrm{p}}^{-s}$, which depends on the GRB spectral shape.

The proton contribution to the GRB emission is typically expressed in terms of the baryon loading factor ${\xi }_{\mathrm{p}}=E_{\mathrm{p}}^{\mathrm{tot}}/E_{\mathrm{iso}}$, defined as the ratio of the total energy in nonthermal protons $E_{\mathrm{p}}^{\mathrm{tot}}$ to the isotropic-equivalent bolometric energy released in photons during the prompt phase E_iso. Typically, ξ_p ≳ 20 is required for the GRBs to be the sources of the observed UHECR flux (for a review, see Kimura 2022, and references therein). E_iso is defined through the observed GRB fluence during the prompt phase S_iso and the luminosity distance d_L at a given redshift z:

$$\begin{array}{c}\hfill E_{\mathrm{iso}}=\frac{4\pi d_{\mathrm{L}}^2{S}_{\mathrm{iso}}}{1+z}·\end{array}$$(5)

In turn, the GRB fluence S_iso is given by the product of the observed GRB flux F_iso (in units erg s⁻¹ cm⁻²), and the total observation time of the prompt phase T₁₀₀ as defined in GRBweb (Coppin 2022).

In terms of the GRB parameters, the proton spectrum can then be written as (Kimura 2022)

$$\begin{array}{c}\hfill E_{\mathrm{p}}^2\frac{\mathrm{d}{N}_{\mathrm{p}}}{\mathrm{d}{E}_{\mathrm{p}}}={\xi }_{\mathrm{p}}{f}_{\mathrm{p}}{E}_{\mathrm{iso}}={\xi }_{\mathrm{p}}{f}_{\mathrm{p}}\frac{4\pi d_{\mathrm{L}}^2{S}_{\mathrm{iso}}}{1+z},\end{array}$$(6)

where f_p is the fraction of energy in protons that produce neutrinos via the photo-meson interactions. It is calculated in Eq. (10) of Zhang & Kumar (2013) for an injection proton spectrum $\mathrm{d}{N}_{\mathrm{p}}/\mathrm{d}{E}_{\mathrm{p}}\propto E_{\mathrm{p}}^{-s}$ with s = 2 (that we assume through all this work), and for a range of neutrino energy from 0 to ∞. In first approximation, f_p ∼ 1/ln(E_p, max/E_p, min), where E_p, min and E_p, max are respectively the minimum and maximum energy of the nonthermal proton spectrum.

On average, one expects to observe the same flux ϕ_ν of all neutrino flavors at Earth. In the following, we focus on muon neutrinos and antineutrinos, that constitute the large majority of the events of the IceCube data analyzed by this work (see Sect. 4). The muon neutrino+antineutrino fluence $F_{\nu }=\int E_{\nu }^2{\varphi }_{\nu }(E_{\nu },t)\text{d}t$, at the neutrino energy E_ν, can be linked to the proton spectrum, and ultimately, using Eq. (6), to the observable GRB parameters (Kimura 2022):

$$\begin{array}{c}\hfill F_{\nu }=\frac{1}{8}{f}_{\mathrm{p}\gamma }{f}_{\pi }^{\mathrm{syn}}{f}_{\mu }^{\mathrm{syn}}{E}_{\mathrm{p}}^2{N}_{\mathrm{p}}(E_{\mathrm{p}})\frac{1+z}{4\pi d_{\mathrm{L}}^2}=\frac{1}{8}{\xi }_{\mathrm{p}}{f}_{\mathrm{p}}{f}_{\pi }^{\mathrm{syn}}{f}_{\mu }^{\mathrm{syn}}{f}_{\mathrm{p}\gamma }{S}_{\mathrm{iso}}.\end{array}$$(7)

The parameters of this equation are described below in this paragraph. The factor 1/8 accounts for the fraction of the pion energy carried by each neutrino (∼1/4) and the branching ratio of nonresonant pγ interactions into the π⁺ channel (∼1/2, as stated above). $f_{\pi }^{\mathrm{syn}}$ and $f_{\mu }^{\mathrm{syn}}$, called pion and muon synchrotron suppression factors, are included in Eq. (7) to account for a possible suppression of the high-energy neutrino component due to synchrotron emission of pions and muons that cool down before decaying. The suppression factors are defined as follows:

$$\begin{array}{c}\hfill f_i^{\mathrm{syn}}\approx 1-e^{-\frac{t_i^{\mathrm{syn}}}{t_i^{\mathrm{dec}}}},\end{array}$$(8)

where i indicates the particle type, pion or muon. The typical timescale for synchrotron cooling of a particle with energy E_i in a magnetic field B is $t_i^{\mathrm{syn}}\propto B^{-2}{E}_i^{-1}$. If this is shorter than the decay timescale $t_i^{\mathrm{dec}}\propto E_i$, the suppression factors are $f_i^{\mathrm{syn}}\simeq {(E_i^{\mathrm{syn}}/E_i)}^2$, and the neutrino flux can be suppressed at high energy. Normally, this occurs above an energy threshold $E_{\pi }^{\text{syn}}\sim {10}^{18}$ eV for pions and $E_{\mu }^{\text{syn}}\sim {10}^{17}$ eV for muons, assuming typical GRB parameters (Kimura 2022). Below these thresholds, the suppression factors are negligible, $f_i^{\mathrm{syn}}\simeq 1$. Finally, the factor f_pγ in Eq. (7) is the fraction of protons that participate in the photo-pion production. The functional dependency of f_pγ on the GRB parameters provides an important part of the physical weights of the stacking search, and will be discussed in detail for each (sub)catalog in the following of this section. It is worth noticing that numerical computations of the expected neutrino signals provide factor of ∼3 − 5 more flux in high energy neutrinos (e.g. see Hümmer et al. 2012; He et al. 2012), thus the analytical approximation we use to constrain the amount of nonthermal protons in the GRB jets is conservative.

The expected GRB neutrino flux is a double broken power law, with a low-energy and a high-energy break point, ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$ and ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$ respectively (Waxman & Bahcall 1997; Kumar & Zhang 2015):

$$\begin{array}{c}\hfill {\varphi }_{\nu }(E_{\nu })\propto \{\begin{array}{cc}{\left(\frac{E_{\nu }}{{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}}\right)}^{-s+\beta -1}\hfill & E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}\hfill \\ {\left(\frac{E_{\nu }}{{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}}\right)}^{-s+\alpha -1}\hfill & {\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}<E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}\hfill \\ {\left(\frac{{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}}{{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}}\right)}^2{\left(\frac{E_{\nu }}{{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}}\right)}^{-s+\alpha -3}\hfill & E_{\nu }>{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}\hfill \end{array}.\end{array}$$(9)

The low-energy neutrino flux break ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$ can be calculated from the break of the proton energy, assuming a typical neutrino-to-proton energy ratio of 5%. This proton break (see for instance Eq. (29) in Kimura 2022) is due to the photo-pion threshold. Hence, the neutrino low-energy break is:

$$\begin{array}{c}\hfill {\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}=60\left(\frac{\mathrm{\Gamma }}{100}\right)\left(\frac{E_{\gamma ,\mathrm{peak}}}{300\mathrm{keV}}\right){\left(\frac{2}{1+z}\right)}^2\mathrm{TeV}.\end{array}$$(10)

Given that IceCube is mostly sensitive in the energy range from few TeV to few PeV, it is very likely to observe GRB neutrinos in the region around ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$. The high-energy break ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$ is associated to the pion transition from a decay-dominated to a synchrotron-cooling-dominated regime, where the neutrino production is suppressed. Assuming a neutrino-to-pion energy ratio of 25%, the high-energy break can be estimated as ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}\simeq 0.25E_{\pi }^{\mathrm{syn}}\sim {10}^{17}$ eV. Neutrinos above the high-energy break are unlikely to be observed by IceCube, and therefore they are not considered in this analysis. The neutrino spectral index γ in each power-law region is correlated with the proton spectral index s and the photon spectral indices α and β. Below the low-energy break ($E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$), neutrinos are produced by the interactions of protons with the high-energy part of the photon spectrum (spectral index β). In this energy region, the neutrino spectrum scales as ${\varphi }_{\nu }\propto E_{\nu }^{-\gamma }=E_{\nu }^{-s+\beta -1}$. Above the low-energy break, (${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}<E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$), protons can interact with the low-energy photons with spectrum with spectral index α, producing a neutrino spectrum ${\varphi }_{\nu }\propto E_{\nu }^{-\gamma }=E_{\nu }^{-s+\alpha -1}$. Above the high-energy break ($E_{\nu }>{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$), the neutrino production is suppressed by synchrotron cooling, and the spectrum is softened by the synchrotron suppression factors, $f_i^{\text{syn}}\sim E_i^{-2}$. The high-energy part of the neutrino spectrum is thus ${\varphi }_{\nu }\propto E_{\nu }^{-\gamma }=E_{\nu }^{-s+\alpha -3}$. Assuming a proton spectral index s = 2 and typical photon parameters α ≃ 1 and β ≃ 2, the corresponding neutrino spectral indices are γ = 1 ($E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$), γ = 2 (${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}<E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$), and γ = 4 ($E_{\nu }>{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$), respectively.

3.1. Empirical correlations for the prompt emission of GRBs

The stacking weights, that are derived in Sects. 3.2 and 3.3, depend on GRB parameters, such as the isotropic-equivalent bolometric luminosity L_iso and the bulk Lorentz factor Γ. For GRBs with measured redshift, it is trivial to calculate the bolometric luminosity as $L_{\text{iso}}=4\pi d_{\text{L}}^2{F}_{\text{iso}}$ from the observed GRB flux F_iso and luminosity distance d_L. However, the majority of GRBs do not have measured redshifts. While the past IceCube analyses (e.g. Aartsen et al. 2017a; Abbasi et al. 2022) adopted the same benchmark values of the relevant parameters for all GRBs, here we propose a novel approach, that exploits empirical correlations observed in lGRBs to estimate the needed quantities for the prompt phase, as first used in Albert et al. (2020) to estimate the GRB bulk Lorentz factor. Such an approach benefits from the improvements of the last decade in the understanding of the GRB mechanisms, and it is motivated by the wide sample of GRB observations currently available.

Two solid correlations have been found between the GRB luminosity and bulk Lorentz factor (Ghirlanda et al. 2012; Lü et al. 2012), and between the GRB luminosity and the photon energy peak in the rest frame (Yonetoku et al. 2004), E_{γ,peak,300,rest} = (1+z)E_γ,peak,300, where E_γ,peak,300 = E_γ,peak/300 keV. In Fig. 3 we plot and fit the Γ₂ − E_{γ,peak,300,rest} and the L_iso,52 − E_{γ,peak,300,rest} data from the most updated samples in the literature, Ghirlanda et al. (2018) and Yonetoku et al. (2010) respectively, where²Γ₂ = Γ/10² and L_iso, 52 = L_iso/10⁵² erg s⁻¹. Although the original Yonetoku relation in Yonetoku et al. (2004) involves the peak luminosity, here we use the (time-averaged) isotropic luminosity, that is more appropriate given that the analysis described in the following of this paper is based on time-averaged GRB properties and on a search for neutrinos from the prompt phase as a whole. For both relationships, the data suggest a linear log-log dependence. In the first case, assuming homogeneous circum-burst medium around the lGRB progenitors (since Γ is mainly derived from the peak of the afterglow emission), the best-fit parameters are:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_2=(1.2\pm 0.1)E_{\gamma ,\mathrm{peak},300,\mathrm{rest}}^{0.50\pm 0.04}.\end{array}$$(11)

Fig. 3. Correlation plots of the bulk Lorentz factor Γ2 = Γ/102 (a) and luminosity Liso, 52 = Liso/1052 erg s−1 (b) with the photon energy peak in the rest frame Eγ,peak,300,rest for the prompt emission, useful for the catalog of the GRB prompts. The best exponential fit is shown as a dashed line, with the 3σ error bands represented by the blue, shaded regions. GRB data are taken from Ghirlanda et al. (2018), (a) and Yonetoku et al. (2010), (b).

For the second relationship, the best fit returns:

$$\begin{array}{c}\hfill L_{\mathrm{iso},52}=(1.0\pm 0.1)E_{\gamma ,\mathrm{peak},300,\mathrm{rest}}^{1.44\pm 0.06}.\end{array}$$(12)

These relationships are used to estimate the luminosity and Lorentz factor of prompt GRBs for which the redshift is not measured (see below). To derive the Γ, Ghirlanda et al. (2018) assumed constant efficiency of the prompt emission, that is $\eta =\frac{E_{\mathrm{iso}}}{E_{\mathrm{iso}}+E_{\mathrm{kin}}}\approx \frac{E_{\mathrm{iso}}}{E_{\mathrm{kin}}}=0.2$, where E_kin is the total kinetic energy of the jet. The baryon-loading factor in the prompt emission is ξ_p = ϵ_p/η, where ϵ_p is the fraction of total kinetic energy in nonthermal protons. The relevance of these parameters, namely the peak energy, the luminosity, and the bulk Lorentz factor, to determine the neutrino fluence from individual GRBs was first proposed for GRBs of the BATSE catalog in the AMANDA era (Guetta et al. 2004).

In principle, the bulk Lorentz factor can also be inferred by using the Γ₂ − L_iso, 52 relation (fit in Fig. 4) together with the L_iso, 52 − E_{γ, peak, 300, rest} relation in Eq. (12). While leading to the same weighting dependencies derived in Sects. 3.2 and 3.3, this approach is more convoluted, and here we prefer to fit and use the direct relation in Eq. (11).

Fig. 4. Correlation plot of the bulk Lorentz factor Γ2 = Γ/100 with the isotropic luminosity Liso, 52 = Liso/(1052 erg s−1) for the prompt emission. The best exponential fit is shown as a dashed line, with the 3σ error bands represented by the blue, shaded regions. GRB data are taken from Ghirlanda et al. (2018).

3.2. Stacking weights for γ = 2 (${\epsilon }_{\nu ,\mathsf{br}}^{\mathsf{low}}<{\mathit{E}}_{\nu }<{\epsilon }_{\nu ,\mathsf{br}}^{\mathsf{high}}$)

In this subsection, we derive the weights for the prompt and X-ray (sub)catalogs in the energy range ${\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}<E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{high}}$, where the expected neutrino spectral index is γ = 2. To do so, we estimate the factor f_pγ in terms of relevant GRB parameters mentioned above and the size of the emission region R:

$$\begin{array}{c}\hfill f_{\mathrm{p}\gamma }\approx 2\chi (\alpha ,\beta )\left(\frac{2}{1+z}\right)\frac{L_{\mathrm{iso},52}}{{\mathrm{\Gamma }}_2^2{R}_{14}}\frac{1}{E_{\gamma ,\mathrm{peak},300}},\end{array}$$(13)

where $\chi (\alpha ,\beta )=\frac{2(2-\alpha )(\beta -2)}{(\beta -\alpha )(1+\beta )}\approx 0.14$ assuming typical α = 1 and β = 2.3.

In the analysis of the prompt subcatalog, consisting of GRBs with measured redshift, the size of the emission region R can be inferred from the observed time variability δt_obs of the prompt phase (Piran 2004):

$$\begin{array}{c}\hfill R=2c{\mathrm{\Gamma }}^2\frac{\delta t_{\mathrm{obs}}}{1+z}·\end{array}$$(14)

Replacing Eq. (14) in Eq. (13), one finds for the prompt emission

$$\begin{array}{c}\hfill f_{\mathrm{p}\gamma }\approx \frac{2}{3}\chi (\alpha ,\beta )\frac{L_{\mathrm{iso},52}}{{\mathrm{\Gamma }}_2^4}\frac{1}{\delta t_{\mathrm{obs},0}}\frac{1}{E_{\gamma ,\mathrm{peak},300}},\end{array}$$(15)

where δt_obs, 0 = δt_obs/1 s. Combining Eq. (15) with Eq. (7), and neglecting the synchrotron suppression factors, the expected neutrino fluence is estimated as

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{\chi (\alpha ,\beta )}{12ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}\left(\frac{L_{\mathrm{iso},52}{S}_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak},300}}\right)\frac{1}{\delta t_{\mathrm{obs},0}}\frac{{\xi }_{\mathrm{p}}}{{\mathrm{\Gamma }}_2^4}·\end{array}$$(16)

The value of ln(E_p, max/E_p, min) is estimated from simulations. Therefore, the stacking weights for the analysis of the prompt subcatalog are:

$$\begin{array}{c}\hfill {\omega }_z^{\gamma }=\frac{L_{\mathrm{iso}}{S}_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak}}}·\end{array}$$(17)

For this subcatalog, we place constraints on ξ_p vs. Γ assuming benchmark values of the time variability δt_obs.

For the analysis of the prompt catalog, where the redshift is not measured for all GRBs and the luminosity and Lorentz factor can be unavailable, we use the correlations in Eqs. (11) and (12) to replace the dependence of F_ν from Γ₂ and L_iso, 52 in Eq. (16) with the observable E_γ, peak. We obtain the following expected neutrino fluence:

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{0.2\chi (\alpha ,\beta )}{ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}\left(\frac{S_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak},300}^{1.6}}\right)\frac{1}{{(1+z)}^{0.6}}\frac{{ϵ}_{\mathrm{p}}}{\delta t_{\mathrm{obs},0}}{\left(\frac{0.20}{\eta }\right)}^{\frac{1}{2}}·\end{array}$$(18)

Therefore, the stacking weights for the analysis of the full prompt catalog (that includes GRBs with no measured redshift) are defined as

$$\begin{array}{c}\hfill {\omega }^{\gamma }=\frac{S_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak}}^{1.6}}·\end{array}$$(19)

In this case, we constrain ϵ_p vs δt_obs for fixed values of redshift z. We have chosen the fiducial value of η of 0.20 according to the common value inferred from the X-ray and the GeV afterglow modeling of lGRBs (Nava et al. 2014; Beniamini et al. 2016).

When considering the plateau and flare subcatalogs, the relationship in Eq. (14), used to remove the dependence of f_pγ from the size R of the emission region during the prompt phase, cannot be used, as the variability timescale for the X-ray flares and plateaus is not defined. In this case, the parameter R cannot be removed from the equations, and the expected neutrino fluence for the X-ray flares and plateaus is

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{0.15\chi (\alpha ,\beta )}{ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}\left(\frac{L_{\mathrm{X},47}{S}_{\mathrm{iso}}^{\mathrm{X}}}{1+z}\right)\frac{1}{E_{\gamma ,\mathrm{peak},1}}\frac{1}{R_{14}}\frac{{\xi }_{\mathrm{p}}}{{\mathrm{\Gamma }}_1^2}·\end{array}$$(20)

Here, $L_{\mathrm{X},47}=4\pi d_{\mathrm{L}}^2{F}_{\mathrm{X}}/({10}^{47}\mathrm{erg}{\mathrm{s}}^{-1})$ is the X-ray luminosity normalized to 10⁴⁷ erg s⁻¹, F_X and $S_{\mathrm{iso}}^{\mathrm{X}}$ are the observed X-ray flux and fluence respectively, and E_{γ, peak, 1} = E_γ, peak/1 keV is the energy peak normalized at 1 keV. For the plateau and flare subcatalogs, comprising GRBs with measured redshift, the stacking weights are defined as

$$\begin{array}{c}\hfill {\omega }_z^{\mathrm{X}}=\frac{L_{\mathrm{X}}{S}_{\mathrm{iso}}^{\mathrm{X}}}{1+z}·\end{array}$$(21)

In this case, the analysis is used to place constraints on ξ_p vs Γ for benchmark values of R₁₄.

For the full plateau and flare catalogs, comprising GRBs with and without measured redshift, the basic weight is used:

$$\begin{array}{c}\hfill {\omega }^{\mathrm{X}}=S_{\mathrm{iso}}^{\mathrm{X}}.\end{array}$$(22)

In this case, no constraints are placed, as too many parameters in Eq. (20) cannot be estimated in a robust way.

3.3. Stacking weights for γ = 1 (${\mathit{E}}_{\nu }<{\epsilon }_{\nu ,\mathsf{br}}^{\mathsf{low}}$)

In this subsection, we derive the weights for the prompt and X-ray (sub)catalogs in the energy range $E_{\nu }<{\epsilon }_{\nu ,\mathrm{br}}^{\mathrm{low}}$, where the expected neutrino spectral index is γ = 1. The neutrino fluence at 1 TeV can be parametrized as

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{2.1\times {10}^{-3}\chi (\alpha ,\beta )}{ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}{L}_{\mathrm{iso},52}{S}_{\mathrm{iso}}(1+z)\frac{1}{R_{14}}\frac{{\xi }_{\mathrm{p}}}{{\mathrm{\Gamma }}_2^4}·\end{array}$$(23)

For the analysis of the prompt subcatalog, the relationship (14) is adopted to reduce the dependence of F_ν in Eq. (23) from R:

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{3.5\times {10}^{-4}\chi (\alpha ,\beta )}{ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}{L}_{\mathrm{iso},52}{S}_{\mathrm{iso}}{(1+z)}^2\frac{1}{\delta t_{\mathrm{obs},0}}\frac{{\xi }_{\mathrm{p}}}{{\mathrm{\Gamma }}_2^6},\end{array}$$(24)

and the following weights for the GRB prompt emission of GRBs with measured redshifts are defined:

$$\begin{array}{c}\hfill {\omega }_z^{\gamma }=L_{\mathrm{iso}}{S}_{\mathrm{iso}}{(1+z)}^2.\end{array}$$(25)

Constraints are then placed on ξ_p vs. Γ for benchmark values of the variability timescale δt_obs.

When considering the full prompt catalog, the luminosity L_iso, 52 and Lorentz factor Γ₂ can be replaced with E_γ, peak, using the relationships in Eqs. (11) and (12):

$$\begin{array}{c}\hfill F_{\nu }\approx \frac{2.05\times {10}^{-5}\chi (\alpha ,\beta )}{ln(E_{\mathrm{p},\max }/E_{\mathrm{p},\min })}\left(\frac{S_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak},300}^{1.6}}\right){(1+z)}^{0.4}\frac{{ϵ}_{\mathrm{p}}}{\delta t_{\mathrm{obs},0}}{\left(\frac{0.20}{\eta }\right)}^{\frac{1}{4}}.\end{array}$$(26)

Therefore, the weights for the analysis of the prompt catalog are

$$\begin{array}{c}\hfill {\omega }^{\gamma }=\frac{S_{\mathrm{iso}}}{E_{\gamma ,\mathrm{peak}}^{1.6}}·\end{array}$$(27)

As for the case with γ = 2, in this case we constrain ϵ_p vs. δt_obs considering benchmark values of the redshift z.

For the analysis of the flare and plateau subcatalogs, we go back to Eq. (23), and we replace the prompt parameters with the corresponding X-ray parameters. We define the following weights:

$$\begin{array}{c}\hfill {\omega }_z^{\mathrm{X}}=L_{\mathrm{X}}{S}_{\mathrm{iso}}^{\mathrm{X}}(1+z).\end{array}$$(28)

Constraints are set on ξ_p vs. Γ for some values of R₁₄.

For the flare and plateau catalogs, the basic weight is used:

$$\begin{array}{c}\hfill {\omega }^{\mathrm{X}}=S_{\mathrm{iso}}^{\mathrm{X}}.\end{array}$$(29)

Similarly to the case with γ = 2, no constraints are placed from these two catalogs.

So far we have been ignoring the suppresion of the high energy neutrinos by the synchrotron losses of muons and pions (see Eq. (8)). Let us introduce another parameter, ϵ_B which is the fraction of the kinetic energy of the jet in the magnetic field energy, that is

$$\begin{array}{c}\hfill \frac{{ϵ}_B{L}_{\mathrm{iso}}}{4\pi \eta R^{2}c{\mathrm{\Gamma }}^2}\approx B^2/8\pi ,\end{array}$$(30)

where B is the strength of the magnetic field in the comoving frame of the emitting region of a GRB. In the Fig. 5 we show the neutrino suppression factors by the synchrotron losses of pions and muons for the neutrino energies of 1 TeV, 10 TeV, and 100 TeV. One can notice that the magnetic field of > 10⁴ G is enough to suppress the high energy neutrinos in the IceCube energy range. To estimate ϵ_B which corresponds to a value of B ∼ 10⁴ G for a typical long GRB, we use the relation (30) together with Eqs. (11), (12), and (14):

$$\begin{array}{c}\hfill B\approx 3\times {10}^4\mathrm{G}{\left(\frac{{ϵ}_B}{0.01}\right)}^{\frac{1}{2}}{\left(\frac{\eta }{0.20}\right)}^{-\frac{1}{8}}{\left(\frac{E_{\gamma ,\mathrm{peak}}}{300\mathrm{keV}}\right)}^{-0.8}{\left(\frac{\delta t_{\mathrm{obs}}}{0.1}\right)}^{-1}{(1+z)}^{-0.2}.\end{array}$$(31)

Fig. 5. Neutrino suppression factor by the synchrotron losses of muons and pions as the function of the comoving strength of the magnetic field in the emission region of a GRB. Green, red, and blue lines correspond to the energies of neutrinos equal to 1 TeV, 10 TeV, and 100 TeV.

This means that the synchrotron losses can be safely ignored for ϵ_B < 0.01.

4. The IceCube open-access neutrino data

The IceCube Neutrino Observatory (Aartsen et al. 2017b) is a km³ sized telescope located at the South Pole and designed to detect high-energy (E > 100 GeV) astrophysical neutrinos. Since 2011, the detector has consisted of 86 strings embedded in the Antarctic ice at a depth of 1.5−2.5 km, but the data taking period already started in 2006 in partial detector configurations. The strings are equipped with a total of 5160 digital optical modules, each hosting a 10-inch photomultiplier tube. IceCube detects Cherenkov light from neutrino-induced muon tracks and cascades from neutral and charged current neutrino interactions of tau and electron neutrinos.

The analysis presented in this paper uses the data collected by IceCube in 10 years of operation, from April 6, 2008 to July 8, 2018, and publicly released as described in Abbasi et al. (2021a). The data are grouped into five independent sets that reflect the different detector configurations with 40, 59, and 79 strings (from 2008 to 2010), and the different event selections adopted in the full 86-string configuration for the years 2011 and 2012−2018. Details about the filtering and selection to final level of analyses of the five data sets can be found in Table 1 of Aartsen et al. (2020a).

These data are selected for point-source searches, and thus optimized for muon track-like events, that provide the best angular resolution. The sample spans a neutrino energy range between about a hundred GeV to above PeV energies, with median angular resolution below 0.4° above 10 TeV. The roughly 1.1 million events collected in 10 years are mostly dominated by muons from atmospheric neutrino interactions in the Northern Hemisphere (declination δ ≥ −5°), and by high-energy atmospheric muons beyond several TeV or large bundles of low-energy muons in the Southern Hemisphere (δ < −5°).

5. Analysis method

For all the searches of this work, we use an unbinned maximum likelihood method similar to the one used for previous point-source IceCube searches (e.g. Braun et al. 2008, 2010; Aartsen et al. 2016b, 2020a,b; The IceCube Collaboration 2018; Abbasi et al. 2021b,c). Such a method is implemented in a public code, known as PSLab, that has recently been released to the community by the IceCube collaboration. This method is used to search for astrophysical neutrino correlations with individual GRBs (single-source search) and with each (sub)catalog as a whole (stacking search). In the two searches, the tested background hypothesis is the absence of astrophysical neutrinos correlated with individual GRBs for the single-source search, or with the ensemble of GRBs in each (sub)catalog for the stacking search.

5.1. Single-source search

For a GRB g at direction Ω_g = (α_g, δ_g), with right ascension α_g and declination δ_g, observed at the instant t_g and lasting Δt_g, the unbinned likelihood of the single-source search reads as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{L}}_j^g(n_{\mathrm{s}})={\displaystyle \prod _{i=1}^{N_j}}& [\frac{n_{\mathrm{s}}}{N_j}{\mathcal{S}}_{i,j}^g(E_i,{\mathbf{\Omega }}_i,{\sigma }_i|{\mathbf{\Omega }}_g,t_g,\mathrm{\Delta }{t}_g,{\sigma }_{\mathrm{GRB},g},\gamma )\hfill \\ \hfill & +\frac{N_j-n_{\mathrm{s}}}{N_j}{\mathcal{B}}_{i,j}^g({\delta }_i,E_i)].\hfill \end{array}\end{array}$$(32)

As the 10-yr IceCube data consist of five independent sets, that differ for the detector configuration and event selection criteria, the index j denotes the IceCube data set (from 1 to 5) that contains the time range of the analyzed phase of GRB g. The index i runs over the events collected in the IceCube data set j, for a total of N_j events. n_s is the amount of signal-like events observed from the tested GRB g, and it is the only free likelihood parameter. ${\mathcal{S}}_{i,j}^g$ and ${\mathcal{B}}_{i,j}^g$ are the single-source signal and background probability density functions (PDFs) respectively, that depend on the GRB parameters (Ω_g, t_g, Δt_g) and on the localization error σ_GRB, g, as well as on the assumed spectral index γ of the astrophysical neutrino flux, on the energy proxy E_i, on the reconstructed direction Ω_i = (α_i, δ_i), and on the estimated angular uncertainty σ_i of the ith event.

The signal and background PDFs can be decomposed into the product of a spatial, an energy, and a time PDF. The signal spatial PDF, also called Point-Spread Function (PSF), is modeled as a symmetric bivariate Gaussian distribution, centered at the GRB location and with a variance given by the quadratic sum of the neutrino angular uncertainty and the GRB localization error: ${\sigma }^2={\sigma }_i^2+{\sigma }_{\text{GRB},g}^2$. For large standard deviations, σ > 10°, the PSF extends over a nonnegligible fraction of the sky, and following the prescription adopted in one of the previous IceCube searches (Aartsen et al. 2017a), the PSF is replaced by the first-order nonelliptical component of the Kent distribution (Kent 1982):

$$\begin{array}{c}\hfill \frac{\kappa }{4\pi sinh(\kappa )}{e}^{\kappa cos(\mathrm{\Delta }{\mathrm{\Psi }}_i^g)}.\end{array}$$(33)

Here, $\mathrm{\Delta }{\mathrm{\Psi }}_i^g$ is the angular distance between the reconstructed direction of the event i and the location of the GRB g, and κ = 1/σ² is called concentration parameter. The Kent distribution is the extension of the bivariate Gaussian distribution onto a spherical support, and the two distributions do not differ significantly for large values of the concentration parameter. The signal energy PDF is modeled as a power law with fixed spectral index, E^−γ, as also adopted in both past IceCube searches from GRBs (Aartsen et al. 2017a; Abbasi et al. 2022). The analysis is performed by independently assuming γ = 1 or γ = 2, that correspond to the expected shape of the neutrino spectrum below the low-energy and the high-energy break, respectively (Kimura 2022). The signal time PDF is constant during the period of the tested GRB, [t_g, t_g + Δt_g], and null otherwise. For the analysis of the prompt phase, Δt_g is taken as the full prompt duration T₁₀₀, while for the X-ray afterglow it is derived from the automatic analysis of the Swift/XRT products (Evans et al. 2009).

As the large majority of the IceCube events are of atmospheric origin, the background spatial PDF is determined through a data-driven method, also referred to as “scramble method”, where the time of the events is randomized within the uptime of the detector, and the right ascension is corrected accordingly, assuming fixed local coordinates (i.e., azimuth and zenith). Due to the peculiar position of the IceCube experiment at the geographical South Pole, the spatial background PDF depends only on the declination, and it is independent of the particular right ascension of the GRBs. Furthermore, since the directional reconstruction efficiency is slightly enhanced if the track events are aligned with the strings of the IceCube detector, an azimuth-dependent correction is included in the background spatial PDF, as described in Abbasi et al. (2021b). The background energy PDF is also estimated from scrambled data, and it depends on the event declination and energy proxy. The background time PDF is uniform within the time boundaries of each IceCube data set. We neglect the seasonal variations of the atmospheric neutrino and muon background, which is on the order of ∼4% for neutrino events (Heix et al. 2020) and ∼8% for downgoing muons (Tilav et al. 2019).

5.2. Stacking search

The stacking likelihood is based on the same PDFs and parameters as described for the single-source search. It is obtained by extending the single-source likelihood in Eq. (32) to include all the five IceCube data sets, and by replacing the signal PDF with the sum of single-source signal PDFs (for simplicity, the dependence on the neutrino and GRB parameters is omitted):

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{stack}}(n_{\mathrm{s}})={\displaystyle \prod _{j=1}^5\prod _{i=1}^{N_j}\left[\sum _{g\in j}(\frac{n_{\mathrm{s},j}^g{\mathcal{S}}_{i,j}^g}{N_j}+\frac{N_j-n_{\mathrm{s},j}^g}{N_j}{\mathcal{B}}_{i,j}^g)\right].}\end{array}$$(34)

The sum runs across all GRBs g occurring in the livetime of the data set j. $n_{\text{s},j}^g$ is the expected amount of signal-like events in the data set j originating from GRB g, and it is related to the “cumulative” signal-like neutrino events n_s through the following equation:

$$\begin{array}{c}\hfill n_{\mathrm{s},j}^g=\frac{w_g{R}_j^{\mathrm{IC}}({\delta }_g,\gamma )\mathrm{\Delta }{T}_j}{{\sum }_j{\sum }_{g\in j}{w}_g{R}_j^{\mathrm{IC}}({\delta }_g,\gamma )\mathrm{\Delta }{T}_j}{n}_{\mathrm{s}}.\end{array}$$(35)

Here, $R_j^{\text{IC}}({\delta }_g,\gamma )$ is the detector acceptance for a neutrino flux ∝E^−γ related to the data set j at declination δ_g, and ΔT_j = ∑_g ∈ jΔt_g is the overall time extension of all GRBs in the data set j. w_g is the theoretical weight on the relative contribution of GRB g to the overall neutrino flux, and it is computed as described in Sect. 3. By definition, $n_{\text{s}}={\sum }_j{\sum }_{g\in j}{n}_{\text{s},j}^g$, that is in the stacking search n_s is the amount of signal-like neutrinos from all the stacked GRBs and across the entire analysis period.

5.3. Test statistics and p value

In both the single-source and the stacking searches, the likelihood is maximized with respect to its parameter. The maximum-likelihood parameter is denoted with a hat, ${\widehat{n}}_{\mathrm{s}}$. By setting n_s = 0, the likelihood describes the background-only hypothesis, and it is thus referred to as background likelihood. The log-likelihood ratio is then used to define the test statistics (TS) of the analysis:

$$\begin{array}{c}\hfill \mathrm{TS}=-2ln\left(\frac{\mathcal{L}(n_{\mathrm{s}}=0)}{\mathcal{L}({\widehat{n}}_{\mathrm{s}})}\right).\end{array}$$(36)

To avoid loss of generality, in the above definition, the likelihood ℒ is written with no label and it is intended to be the single-source likelihood ${\mathcal{L}}_j^g$ (Eq. (32)) or the stacking likelihood ℒ_stack (Eq. (34)), depending on the analysis. To quantify the significance of an observed TS, the analyses are repeated for background pseudo-experiments generated with the scramble method, each returning a “background-like” TS. The pretrial p value of an analysis is then estimated as the fraction of background pseudo-experiments that produce a TS value larger than the one observed in the data. When needed, the pretrial p value is eventually corrected for the so-called “look-elsewhere effect” into a post-trial p value that takes into account the number of searches performed. The detail of such a correction are provided in the appropriate Sect. 6. In the following, we refer to cases in which TS = 0 as “under-fluctuations”, and cases in which TS > 0 as “over-fluctuations”.

The major difference of the analysis described in this paper and the previous IceCube analyses of GRB neutrinos (Aartsen et al. 2017a; Abbasi et al. 2022) lies in the stacking method. In the stacking analysis in Aartsen et al. (2017a), the TS is calculated by summing the TS of the individual GRBs from the single-source search. However, the TS distribution is generally declination-dependent, and it might differ for GRBs at different declinations. In Abbasi et al. (2022), the stacking method coincides with the one adopted for this work, but the individual GRB contributions in Eq. (35) are only calculated based on the detector acceptance, implicitly assuming an equal GRB fluence at Earth. As a matter of fact, both analyses prefer GRBs that are located in strategic declination bands for IceCube, regardless of the intrinsic properties of the GRBs. In this work, we have proposed a weighted stacking method that, based on physical motivations, estimates the GRB contribution to the expected neutrino fluence by also incorporating relevant GRB observational properties expressed by the weights.

6. Results

The results of the single-source search from the plateau catalog are summarized in Table 2. The search returns GRB 140518A as the most significant GRB, and one signal-like neutrino event identified in temporal coincidence with the X-ray plateau. With a reconstructed muon energy of 3.3 TeV and a directional angular uncertainty of 1.6°, this event lies 3.3° from the location of the source. The corresponding pretrial p value of GRB 140518A is p_loc = 1.1 × 10⁻³ (3.1σ) for a spectral index γ = 1, and p_loc = 1.6 × 10⁻³ (3.0σ) for γ = 2. The post-trial p value is calculated by repeating the analysis on many background pseudo-experiments with the corresponding spectral index, and counting the fraction of pseudo-experiments that produces a smaller p value than the one observed in the data. This translates into a post-trial p value of 69% for both spectral indices, thus not significant.

All the GRBs of the X-ray flare catalog are under-fluctuations, except for GRB 120213A when γ = 2 is considered, as reported in Table 3. In this case, the data contain one neutrino event in temporal coincidence with the X-ray flare. However, the reconstructed muon energy for this event is only 562 GeV, and the reconstructed track direction is displaced with respect to the source location by 5.8°, with an angular uncertainty of 2.9°. This results in a loose association of such neutrino with GRB 120213A, thus explaining the likelihood preference to fit only a fractional event, ${\widehat{n}}_{\mathrm{s}}=0.4$. The pretrial p value of this GRB is p_loc = 2.9 × 10⁻² (1.9σ), corresponding to a post-trial p value of 77%.

The number of under-fluctuating GRBs (∼96% for the plateau catalog and ∼99.5% for the flare catalog) is much larger than the number of over-fluctuating, and for the flare catalog only one over-fluctuating GRB is identified. This is due to the short duration of the time window in which neutrinos are searched, that reduce the expected amount of events for a fixed event rate. This is particularly true for flares, whose associated search time window ranges between few tens to few hundreds of seconds. In the flare catalog, the expected average under-fluctuation rate under the background hypothesis is ≳99%. However, the very low event rate expected for the flare catalog search makes this analysis almost background-free, suggesting that even the smaller amount of signal can be spotted in the data, potentially carrying a high (pretrial) significance. This can be easily understood by looking at the background and signal TS distributions of one example GRB from the flare catalog, notably the over-fluctuating GRB 120213A, shown in Fig. 6. These distributions are obtained by running the single-source analysis at the location of GRB 120213A for many background and signal pseudo-experiments: the former are produced with the scramble method described in Sect. 5.1, and the latter by additionally injecting signal-like neutrino events simulated with the PSLab code (The IceCube Collaboration 2022).

Fig. 6. (a) Example of TS distributions for the background (injected ns = 0, with 100 000 generated pseudo-experiments) and various signal hypotheses (injected ns > 1, with 10 000 pseudo-experiments each) at the location of GRB 120213A, for the single-source search of the flare catalog and assuming spectral index γ = 2. The signal distributions are very different from the background distribution already for ns = 1, suggesting that few signal events can lead to a high pretrial significance. (b) Background TS distribution, the same as in (a), but with different binning and axis ranges. The bin at TS = 0 corresponds to under-fluctuations: for this GRB, about 96% of the background pseudo-experiments result in under-fluctuations.

All the results of the single-source analyses are compatible with the background hypothesis. As such, 90% confidence-level (CL) upper limits on the neutrino fluence from individual non-underfluctuating sources are placed. The upper-limit neutrino fluence for the single-source searches is defined as follows:

$$\begin{array}{c}\hfill F_{\nu }(E_{\nu })=E^2\frac{\mathrm{d}{N}_{\nu }}{\mathrm{d}{E}_{\nu }\mathrm{d}A}=F_{90\%}{\left(\frac{E_{\nu }}{\mathrm{TeV}}\right)}^{2-\gamma }.\end{array}$$(37)

The factor F_90%, corresponding to the normalization of the upper-limit fluence at the benchmark value of 1 TeV, is reported in the last column of Tables 2 and 3 for the related sources.

The results of the stacking searches are summarized in Table 4. These searches produce underfluctuations in all the (sub)catalogs, except for the analysis of the prompt catalog in the Northern Hemisphere with spectral index γ = 2, with no requirements on the measured redshift. In this case, the stacking likelihood fits ${\widehat{n}}_{\mathrm{s}}\simeq 1$ event, resulting in a pretrial p value of p_loc = 3.5 × 10⁻² and a post-trial p value of 13%. The post-trial correction of the stacking search is due to performing several searches of this kind, on different hemispheres, on different catalogs, with different spectral indices (γ = 1 and γ = 2), and with different requirements on the measured GRB redshift.

As no significant astrophysical correlation is observed by the stacking searches, these analyses are used to place 90% CL upper limits on the cumulative neutrino flux from each (sub)catalog. The upper-limit neutrino flux for the stacking searches is defined as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill E_{\nu }^2{\varphi }_{\nu }^{\mathrm{Stack}}(E_{\nu })& =E^2\frac{\mathrm{d}{N}^{\mathrm{Stack}}}{\mathrm{d}E\mathrm{d}A\mathrm{d}t\mathrm{d}\mathrm{\Omega }}\hfill \\ \hfill & =F_{90\%}^{\mathrm{Stack}}{\left(\frac{E_{\nu }}{\mathrm{TeV}}\right)}^{2-\gamma }\frac{1}{\mathrm{\Delta }T}\frac{1}{\mathrm{\Delta }\mathrm{\Omega }}={\varphi }_{90\%}^{\mathrm{Stack}}{\left(\frac{E_{\nu }}{\mathrm{TeV}}\right)}^{2-\gamma },\hfill \end{array}\end{array}$$(38)

where ΔT ≃ 10 yr is the full livetime of the data used for the analysis, and ΔΩ ≃ 2π is the solid angle of each hemisphere. $F_{90\%}^{\mathrm{Stack}}$ and ${\varphi }_{90\%}^{\mathrm{Stack}}$ are the normalizations at 1 TeV on the 90% CL upper limit of the neutrino stacking fluence and the neutrino stacking flux, respectively. The 90% CL upper-limit flux normalization for each (sub)catalog is reported in Table 4. For a better visualization, such upper limits are also displayed in Fig. 7, assuming the same neutrino spectrum $\propto E_{\nu }^{-2}$ for all the GRBs. This figure additionally shows the energy-dependent stacking sensitivity of this analysis, and a comparison of the estimated IceCube-Gen2 sensitivity for the prompt (sub) catalog, assuming an $E_{\nu }^{-2}$ neutrino spectrum. The IceCube-Gen2 (Clark 2021) sensitivity is calculated by considering a 10 times larger detector volume, and by scaling the current effective area (and hence the sensitivity) by a factor 10^2/3, close to the designed factor ∼5 (Aartsen et al. 2014). The flux in Eq. (38) follows a similar definition as for the IceCube search in Abbasi et al. (2022), but it differs from the previous IceCube search in Aartsen et al. (2017a).

Fig. 7. Differential sensitivity (solid lines) and 90% CL upper limits (dashed lines) for the stacking analysis of the plateau (red), flare (yellow), and prompt (green) catalog, assuming γ = 2. Curves are shown for the Northern (left) and Southern (right) Hemisphere, and for catalogs comprising all the GRBs (first row) or the subcatalogs of GRBs with measured redshift (second row). As a comparison, the shaded green region shows the 10-yr differential sensitivity calculated for IceCube-Gen2 (Clark 2021) as explained in the text.

7. Discussion

The upper limits derived in Sect. 6 from the nonobservation of significant events in the various stacking searches are used to constrain relevant quantities of the model discussed in Sect. 3. Our assumption for the prompt emission is based on the empirical correlations between the GRB luminosities and the peak energies, and between the GRB bulk Lorentz factors and the peak energies, as detailed in Yonetoku et al. (2010) and Ghirlanda et al. (2012), respectively. Unlike previous GRB analyses performed by IceCube, that assumed all GRBs to have the same bulk Lorentz factor, fluence, or luminosity, we take into account the large range of values for these parameters (e.g. the luminosity spans five orders of magnitude). We assume the aforementioned correlations to estimate realistic values for these parameters of each GRB.

We obtain constraints on the baryon loading factor ξ_p of the prompt emission as a function of the timescale variability δt_obs and bulk Lorentz factor Γ. Such constraints are shown in Figs. 8 (assuming a neutrino spectral index γ = 2) and 9 (assuming a neutrino spectral index γ = 1), for the two hemispheres and considering measured redshifts or benchmark values (when not measured), as well as some values of the variability timescale δt_obs. Our analysis generally limits the baryon loading factor to ξ ≲ 10, unless extreme GRB values are considered (e.g. Γ > 500, δt_obs ∼ 1 s, z ∼ 10). This seems to disfavor a baryonic origin of the GRB prompt, as also indicated by the absence of significant astrophysical neutrino results. Furthermore, if the synchrotron cooling suppression factors in Eq. (8) are not neglected, we can use Eq. (7) to similarly constrain ϵ_p as a function of the magnetic equipartition parameter ϵ_B, assuming typical values of δt_obs = 0.1 s, η = 0.20 and z = 1. This is shown in Fig. 10, and indicates that for neutrino energies ≲100 TeV (where most of the IceCube events are observed), a large magnetic field B ≳ 10⁴ G (ϵ_B > 0.01, see Eq. (31)) must be considered to accommodate the observations reported in this work. In other words, considering a single-zone model and considering protons and not heavier nuclei, the neutrino limits on the prompt emission phase imply that the fraction of total kinetic energy in nonthermal protons is less than 10% or that the magnetization is so high that GRB associated neutrinos should be found in the lower energies (sub-TeV range). The latter might hint at the fact that GRB jets are mostly magnetic-dominated, as already proposed by some authors (Usov 1992; Thompson 1994; Lyutikov & Blandford 2003; Ghisellini et al. 2020).

Fig. 8. 90% CL upper limits for the baryon equipartition parameter ϵp and the baryon loading parameter ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the prompt catalog (upper plots) and subcatalogs of GRBs with available redshift (bottom plots), assuming a neutrino spectral index γ = 2. Excluded values are shown as shaded regions for different hypotheses of the redshift z (top plots) and of the timescale variability δtobs (bottom plots).

Fig. 9. 90% CL upper limits for the baryon equipartition parameter ϵp and the baryon loading parameter ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the prompt catalog (upper plots) and subcatalog of GRBs with available redshift (bottom plots), assuming a neutrino spectral index γ = 1. Excluded values are shown as shaded regions for different hypothesis of the redshift z (upper plots) and of the timescale variability δtobs (bottom plots).

Fig. 10. 90% CL upper limits for the baryon equipartition parameter ϵp in the Northern (left) and Southern (right) Sky as a function of the magnetic equipartition ϵB. These constraints are obtained with the stacking analysis of the prompt catalog, assuming a neutrino spectral index γ = 2. Excluded values are shown as shaded regions for different values of the neutrino energy Eν. Typical values of the timescale variability δtobs = 0.1 s and redshift z = 1 are used to compute these limits. To derive ϵB, we assumed a typical observed spectrum of lGRB with the peak energy of 300 keV (see Eq. (31)). It is worth notice that the asymptotic value of ϵp in the low-magnetic field regime in the Southern Sky is consistent with the limits computed in Fig. 8 in the same hemisphere, assuming the same parameters.

We additionally constrain the baryon loading factor of the plateau and flare subcatalogs (with available redshift) as a function of the Lorentz factor Γ, for some values of the normalized radius of the GRB emission site R₁₄. In this case, we cannot use the empirical correlations observed for the prompt phase and mentioned above, but we exploit all the available observations of the X-ray afterglows, namely the observed luminosity, fluence, and redshift. The constraints of our analyses are shown in Figs. 11 (assuming a neutrino spectral index γ = 2) and 12 (assuming a neutrino spectral index γ = 1). It should be noticed that in the case of X-ray afterglows, the variability range of the bulk Lorentz factor is typically lower (10 < Γ < 10²) than for prompt (10² < Γ < 10³), (Kimura 2022). The constraints on the X-ray afterglows are looser than those on the prompt phase, and the analysis does not exclude a possible baryonic origin of the plateau and flare emission.

Fig. 11. 90% CL upper limits for the baryon loading factor ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the plateau (upper plots) and flare (bottom plots) subcatalogs of GRBs with available redshift, assuming γ = 2. Excluded values are shown as shaded regions for different hypotheses on the normalized radius of the GRB emission site.

Fig. 12. 90% CL upper limits for the baryon loading factor ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the plateau (upper plots) and flare (bottom plots) subcatalogs of GRBs with available redshift, assuming γ = 1. Excluded values are shown as shaded regions for different hypotheses on the normalized radius of the GRB emission site R14 (see Sect. 3).

The nonobservation of an astrophysical neutrino signal from GRB prompt and afterglow phases is compatible with previous IceCube results, and with the estimated upper limits of GRB neutrinos to the diffuse flux of ≲1% (Aartsen et al. 2017a; Abbasi et al. 2022).

8. Conclusions

In this paper, we analyzed the 10-year IceCube data publicly available (Abbasi et al. 2021a). We used the PSLab code (The IceCube Collaboration 2022), recently released by the IceCube collaboration, to perform an unbinned maximum-likelihood search for spatial and temporal coincidence of astrophysical neutrinos from individual GRBs with X-ray flare and X-ray plateau afterglows. We also performed a stacking search for a cumulative neutrino excess from a flare, a plateau, and a prompt catalog. Unlike previous IceCube stacking searches, that assumed the same fluence at Earth for all GRBs, in our analysis we proposed a stacking scheme based on physically motivated GRB weights and observational properties. Furthermore, while past searches adopted the same benchmark values of relevant parameters for all the GRBs, we fit and used empirical relationships to estimate the GRB luminosities and Lorentz factors in the analysis of the prompt catalog when the redshift is not available. Since such parameters can vary across several orders of magnitude, our novel approach improves the physical reliability of the results.

We did not find any statistically significant neutrino correlations in our searches, consistently with previous IceCube results (Aartsen et al. 2017a; Abbasi et al. 2022). The nonobservation of neutrinos is therefore used to place upper limits on the individual and stacking flux, and to constrain general parameters of fundamental interest, such as the baryon loading, the Lorentz factor, and the magnetic field of the jet in a general single-zone fireball model. The results of the stacking analysis for the prompt phase suggest a typical baryon loading factor ξ_p ≲ 10, thus excluding the hypothesis of a baryonic origin of the jets and supporting the scenario of a magnetic-dominated ejecta. In this scenario, GeV neutrinos might be expected as a consequence of the highly efficient synchrotron cooling of pions and muons, and potentially detectable by IceCube DeepCore and KM3NeT/ORCA (Zegarelli et al. 2022). The constraints of the afterglow phase are less tight, but it is worth mentioning that this is the first analysis that specifically targets GRB X-ray plateaus and flares.

Acknowledgments

We thank the IceCube Collaboration which offered the optimal environment for the analysis code development and for its many applications and for having released the 10-year data sample. This work also made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We are thankful to Christopher Wiebusch for his detailed comments to the paper, and to the Publication Committee of IceCube for having discussed it. G.O. thanks Annalisa Celotti and Dafne Guetta for fruitful discussions. The work was financed by the Swiss National Foundation grant n. 200020_178918 and by the University of Geneva. G.O. and M.B. acknowledge financial support from the AHEAD2020 project (grant agreement n. 871158). B.B. and M.B. acknowledge financial support from MUR (PRIN 2017325 grant 20179ZF5KS).

References

All Tables

All Figures

	Fig. 1. Examples of a prompt emission light curve (a) observed by Fermi/GBM, a flare (b) and a plateau (c) of the X-ray afterglow observed by Swift/XRT.
In the text

	Fig. 2. Illustration of the photo-pion processes that lead to the production of high-energy neutrinos inside the dissipation region of GRB jets.
In the text

Fig. 3. Correlation plots of the bulk Lorentz factor Γ2 = Γ/102 (a) and luminosity Liso, 52 = Liso/1052 erg s−1 (b) with the photon energy peak in the rest frame Eγ,peak,300,rest for the prompt emission, useful for the catalog of the GRB prompts. The best exponential fit is shown as a dashed line, with the 3σ error bands represented by the blue, shaded regions. GRB data are taken from Ghirlanda et al. (2018), (a) and Yonetoku et al. (2010), (b).

In the text

	Fig. 4. Correlation plot of the bulk Lorentz factor Γ2 = Γ/100 with the isotropic luminosity Liso, 52 = Liso/(1052 erg s−1) for the prompt emission. The best exponential fit is shown as a dashed line, with the 3σ error bands represented by the blue, shaded regions. GRB data are taken from Ghirlanda et al. (2018).
In the text

	Fig. 5. Neutrino suppression factor by the synchrotron losses of muons and pions as the function of the comoving strength of the magnetic field in the emission region of a GRB. Green, red, and blue lines correspond to the energies of neutrinos equal to 1 TeV, 10 TeV, and 100 TeV.
In the text

Fig. 6. (a) Example of TS distributions for the background (injected ns = 0, with 100 000 generated pseudo-experiments) and various signal hypotheses (injected ns > 1, with 10 000 pseudo-experiments each) at the location of GRB 120213A, for the single-source search of the flare catalog and assuming spectral index γ = 2. The signal distributions are very different from the background distribution already for ns = 1, suggesting that few signal events can lead to a high pretrial significance. (b) Background TS distribution, the same as in (a), but with different binning and axis ranges. The bin at TS = 0 corresponds to under-fluctuations: for this GRB, about 96% of the background pseudo-experiments result in under-fluctuations.

In the text

Fig. 7. Differential sensitivity (solid lines) and 90% CL upper limits (dashed lines) for the stacking analysis of the plateau (red), flare (yellow), and prompt (green) catalog, assuming γ = 2. Curves are shown for the Northern (left) and Southern (right) Hemisphere, and for catalogs comprising all the GRBs (first row) or the subcatalogs of GRBs with measured redshift (second row). As a comparison, the shaded green region shows the 10-yr differential sensitivity calculated for IceCube-Gen2 (Clark 2021) as explained in the text.

In the text

Fig. 8. 90% CL upper limits for the baryon equipartition parameter ϵp and the baryon loading parameter ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the prompt catalog (upper plots) and subcatalogs of GRBs with available redshift (bottom plots), assuming a neutrino spectral index γ = 2. Excluded values are shown as shaded regions for different hypotheses of the redshift z (top plots) and of the timescale variability δtobs (bottom plots).

In the text

Fig. 9. 90% CL upper limits for the baryon equipartition parameter ϵp and the baryon loading parameter ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the prompt catalog (upper plots) and subcatalog of GRBs with available redshift (bottom plots), assuming a neutrino spectral index γ = 1. Excluded values are shown as shaded regions for different hypothesis of the redshift z (upper plots) and of the timescale variability δtobs (bottom plots).

In the text

Fig. 10. 90% CL upper limits for the baryon equipartition parameter ϵp in the Northern (left) and Southern (right) Sky as a function of the magnetic equipartition ϵB. These constraints are obtained with the stacking analysis of the prompt catalog, assuming a neutrino spectral index γ = 2. Excluded values are shown as shaded regions for different values of the neutrino energy Eν. Typical values of the timescale variability δtobs = 0.1 s and redshift z = 1 are used to compute these limits. To derive ϵB, we assumed a typical observed spectrum of lGRB with the peak energy of 300 keV (see Eq. (31)). It is worth notice that the asymptotic value of ϵp in the low-magnetic field regime in the Southern Sky is consistent with the limits computed in Fig. 8 in the same hemisphere, assuming the same parameters.

In the text

Fig. 11. 90% CL upper limits for the baryon loading factor ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the plateau (upper plots) and flare (bottom plots) subcatalogs of GRBs with available redshift, assuming γ = 2. Excluded values are shown as shaded regions for different hypotheses on the normalized radius of the GRB emission site.

In the text

Fig. 12. 90% CL upper limits for the baryon loading factor ξp in the Northern (left) and Southern (right) Sky. These constraints are obtained with the stacking analysis of the plateau (upper plots) and flare (bottom plots) subcatalogs of GRBs with available redshift, assuming γ = 1. Excluded values are shown as shaded regions for different hypotheses on the normalized radius of the GRB emission site R14 (see Sect. 3).

In the text