© The Authors 2026

1. Introduction

The KM3NeT Collaboration reported the observation of a neutrino with an estimated $E_{\nu }={220}_{-110}^{+570}$ PeV with the KM3NeT/ARCA detector in a partial configuration. It is the highest energy neutrino observed to date. The event represents the first neutrino of presumable astrophysical origin observed in the ultra-high-energy (UHE) regime. Subsequent studies explored how this event fits into the global neutrino landscape, taking into account the lack of positive detections reported so far by the IceCube (IC; Meier 2024) and Pierre Auger (Auger; Halim et al. 1488.) Observatories. By combining these different (non-)observations, the most likely single-flavour diffuse astrophysical neutrino flux required to produce an event such as KM3-230213A is $E^2{\mathrm{\Phi }}_{\nu +\overline{\nu }}^{1\mathrm{f}}=7.5_{-4.7}^{+13.1}\times {10}^{-10}$ GeV cm⁻² s⁻¹ sr⁻¹, assuming an E⁻² neutrino flux (Adriani et al. 2025).

Since its discovery in 2013 (Aartsen et al. 2013), the diffuse flux of astrophysical neutrinos between TeV and PeV energies has been extensively investigated (see e.g. Halzen & Kelley 2024), and a multitude of potential sources have been scrutinised. Despite the identification of a few likely sources (Aartsen et al. 2018b; Abbasi et al. 2022a; Sclafani et al. 2024), the origin of the majority of the diffuse neutrino flux remains unknown. The lack of observed neutrino multiplets (Abbasi et al. 2025) further suggests that the population of neutrino-producing astrophysical objects consists of relatively dim, abundant, and isotropically distributed sources (Murase & Waxman 2016).

Gamma-ray bursts (GRBs) are the most energetic transient events observed in the Universe in electromagnetic wavebands and are potential sources of UHE cosmic rays (E  ≳  10¹⁸ eV) (Waxman 1995; Vietri 1995). Consequent neutrino production has been predicted from interactions of cosmic-ray protons with photons within the fireball in internal shocks (Waxman & Bahcall 1997). High-energy neutrinos may also be created when the outgoing blast wave, in which particles are accelerated by external shocks, interacts with matter and radiation fields surrounding the GRB (Waxman & Bahcall 2000; Dai & Lu 2001). Previous searches for GRB neutrinos yielded no detection (Aartsen et al. 2016; Albert et al. 2017; Abbasi et al. 2022c). These searches have focused on analysing triggering GRBs – GRBs bright enough in gamma rays to initiate multi-wavelength follow-ups – and coincident neutrinos from varying time frames around the prompt emission phase of the GRBs. Although these searches have so far not resulted in direct detection, constraints have been put on the ratio of energy between protons and electrons, known as ‘baryon loading’ (Rees & Meszaros 1994). This method, albeit successful in constraining individual GRB models, is not sensitive to the potentially large population of GRBs undetected by gamma-ray satellites (see e.g. Li et al. 2025). Furthermore, the imposed constraint on the baryon loading depends on the GRB prompt emission model, which is uncertain, as well as on the model parameters. The purpose of this paper is therefore to investigate whether a larger population of GRBs can produce a significant fraction of the diffuse UHE neutrino flux, with emphasis on the undetected part of this population.

Previous constraints on the diffuse UHE neutrino flux have been set by both the IceCube and Pierre Auger observatories (Aartsen et al. 2021a; Aab et al. 2019). The lack of detection of UHE neutrinos by these observatories made the upper limits the strongest constraints available until the detection of KM3-230213A.

In this paper, we use the recent observation of a UHE neutrino event to constrain the baryon loading of GRB blast waves and the density of the medium in which the blast wave interacts to produce PeV–EeV neutrinos. We used the first-ever observation of a UHE neutrino to constrain the total contribution of GRB blast waves to the diffuse UHE neutrino flux and consequently estimated some of the relevant model parameters. This paper is structured as follows: In Section 2 we outline the model under consideration for GRB blast wave neutrino production and how a large population of these GRBs can generate a diffuse extragalactic neutrino flux capable of producing KM3-230213A. In Section 3, we present an innovative technique to calculate the diffuse UHE neutrino flux, and the statistical framework used to put constraints on the GRB model parameters. In Section 4 we report and discusses the obtained constraints and their implications on GRBs. Finally, we summarise our findings in Section 5 and highlight our conclusions therein.

2. Diffuse GRB neutrino flux in the PeV–EeV range

We investigated the possible constraints enforced by KM3-230213A on specific model parameters of long-duration GRBs (lGRBs) by considering the contribution of a large population of lGRBs, up to redshift z = 5, to the diffuse neutrino flux at UHEs. Despite a more prominent prompt emission phase at lower neutrino energies (Waxman & Bahcall 1997), a significant amount of energy in lGRBs is expected to be converted to kinetic energy in the form of protons propagating outwards from the GRB in a blast wave that subsequently interacts with surrounding matter and radiation fields (Razzaque 2013). The kinetic energy of the blast wave, E_k, is connected to the inferred gamma-ray luminosity during the prompt phase, L_γ, through the relation

$$\begin{array}{c}\hfill \frac{E_k}{L_{\gamma }}=f_b\eta t^{\ast },\end{array}$$(1)

where f_b is the baryon loading ratio, η is the efficiency of converting kinetic energy to gamma-ray energy, and t^* is the timescale for the prompt emission. In our calculations, we set η = 0.2 (Fan & Piran 2006) and adopted a normal distribution for log(t^*), fit to the distribution of lGRBs reported in von Kienlin et al. (2020), as log(t^*)∼𝒩(μ, σ²) with μ = 1.44 and σ = 0.50.

2.1. GRB blast wave models

After the prompt emission phase, relativistic ejecta from the GRB drive a blast wave interacting in the gaseous media surrounding the GRB progenitor system. The blast wave expands and cools adiabatically, and after a deceleration time (t_dec), its evolution is described by a simple similarity solution (Blandford & McKee 1976). Synchrotron radiation by electrons in the surrounding gas, accelerating and cooling in the magnetic field in the shock region (forward shock), explains the radio-to-gamma-ray afterglow emission from GRB afterglows (Mészáros & Rees 1997; Sari et al. 1998). Protons are co-accelerated with electrons to UHEs in the reverse shock, providing material to the ejecta and in the blast wave itself. These protons then interact with synchrotron photons to produce UHE neutrinos through photopion interactions (see e.g. Waxman & Bahcall 1997; Dai & Lu 2001; Murase 2007; Razzaque 2013). In the following, we consider two different scenarios: (1) an adiabatic blast wave interacting with the interstellar medium (ISM) of constant density and (2) a wind-type environment (WIND) around the GRB progenitor system with its density falling as R⁻², where R is the distance from the GRB. Protons accelerated in the forward shock interact with afterglow synchrotron photons to produce UHE neutrinos. Further details on the GRB blast wave model and neutrino production model used in this paper can be found in Razzaque (2013), Razzaque & Yang (2015).

2.2. GRB population

To understand how a large population of lGRBs contributes to the diffuse astrophysical UHE neutrino flux, their distributions throughout the Universe need to be considered. We adopted the redshift-luminosity distribution obtained from the Swift and Fermi-GBM observations of lGRBs (Banerjee et al. 2021) to construct the evolution function for the lGRB population. The distributions of lGRBs in gamma-ray luminosity and redshift are independent of one another, i.e. Ψ(L_γ, z) = ψ(L_γ)ρ(z), where ψ(L_γ) and ρ(z) denote the luminosity function and redshift distribution, respectively. The luminosity function ψ(L_γ) is modelled as a broken power law and is given by

$$\begin{array}{c}\hfill \psi (L_{\gamma })=c_0[{\left(\frac{L_{\gamma }}{L_0}\right)}^{-\alpha }+ϵ{\left(\frac{L_{\gamma }}{L_0}\right)}^{-\beta }],\end{array}$$(2)

where α = 1.33, β = 1.42, and ϵ = 1, and the break luminosity is L₀ = 3 × 10⁵³ erg s⁻¹. The normalisation constant c₀ is fixed by integrating this expression over the considered luminosity range: c₀ = [∫ψ(L_γ)dL]⁻¹. The redshift distribution function is

$$\begin{array}{c}\hfill \rho (z)={\rho }_0\frac{{(1+z)}^{2.7}}{1+{\left(\frac{1+z}{2.9}\right)}^{5.6}},\end{array}$$(3)

with ρ₀ = 5.5 Gpc⁻³ yr⁻¹. The values for both expressions are the best-fit values from Banerjee et al. (2021). The luminosity (L_γ) considered is the intrinsic isotropic luminosity as inferred from observed GRBs in the 8 keV–40 MeV Fermi-GBM and Swift energy range. The lower and upper bounds we used for the luminosity function are 10⁴⁹ erg s⁻¹ and 10⁵⁴ erg s⁻¹, and we integrated the lGRB population out to redshift z_max = 5, as the number of GRBs decreases significantly at high redshift, and the contributions from far-away GRBs (z > 5) to a diffuse flux become negligible.

2.3. Diffuse neutrino flux

The final step in our model consists of combining the individual neutrino fluxes from long GRBs with their cosmological evolution and gamma-ray luminosity distribution. The total diffuse flux expected on Earth is

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{\nu }^{\mathrm{tot}}(E_{\nu },E_k,n_0;\theta )& ={\displaystyle {\int }_0^{z_{max}}}\frac{1}{1+z}\frac{\mathit{dV}}{\mathit{dz}}\hfill \\ \hfill & \times {\displaystyle {\int }_{L_1}^{L_2}}\mathrm{\Psi }(L_{\gamma },z)S_{\nu }(E_{\nu },E_k,n_0;\theta )d{L}_{\gamma }dz,\hfill \end{array}$$(4)

where the 1/(1 + z) factor corrects for the time dilation between distant lGRBs and the measured flux. The term S_ν is the neutrino fluence, i.e. integrated flux of all flavours over a time t ≥ t_dec for a GRB of luminosity (L_γ) at a redshift z, where an adiabatic blast wave is evolving in the ISM of constant density n₀. We denote all other model parameters with θ. A similar expression holds for a wind-type medium by replacing n₀ with the wind parameter A_*. A detailed calculation of S_ν can be found in Razzaque (2013) and in Appendix A. To compute the diffuse flux, we integrated over the comoving volume. The differential comoving volume element is given by

$$\begin{array}{c}\hfill \frac{\mathit{dV}}{\mathit{dz}}=\frac{4\pi c}{1+z}\left|\frac{\mathit{dt}}{\mathit{dz}}\right|d_L^2,\end{array}$$(5)

where d_L is the luminosity distance (Hogg 1999). The cosmic time, dt/dz, is given by

$$\begin{array}{c}\hfill \frac{\mathit{dt}}{\mathit{dz}}=\frac{-1}{H_0(1+z)\sqrt{{\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}.\end{array}$$(6)

For the calculations in this work, we used Ω_m = 0.286, Ω_Λ = 1 − Ω_m = 0.714, and H₀ = 69.32 km s⁻¹ Mpc⁻¹ (Ade et al. 2014).

3. Analysis method

To calculate the total diffuse neutrino flux from the GRB blast wave models outlined above, we used the nested sampling algorithm as implemented in the UltraNest package (Buchner 2021). In this approach, the parameter-free MLFriend algorithm is applied to sample a constrained likelihood and iteratively converge on the global maximum through a bootstrapping method (Buchner 2016, 2019). By defining the total integrand of our diffuse neutrino flux as the likelihood, the integration is performed by Bayesian inference with the aforementioned nested sampling. The computed marginal likelihood, also known as evidence, corresponds to the evaluated integral, i.e. the total diffuse neutrino flux. With this approach, we continuously varied all integration parameters associated with the diffuse flux calculation without making any a priori assumptions.

The most likely diffuse UHE neutrino flux responsible for producing KM3-230213A has been calculated and reported in Adriani et al. (2025). By considering that no events have been reported in the IceCube Extremely-High-Energy (IC-EHE) or sensitive Auger selections, the UHE neutrino flux of a single flavour was estimated to be $E^2{\mathrm{\Phi }}_{\nu +\overline{\nu }}^{1\mathrm{f}}=(7.5_{-4.7}^{+13.1})\times {10}^{-10}\mathrm{GeV}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}{\mathrm{sr}}^{-1}$ in the 90% reconstructed energy range (Aiello et al. 2025). To constrain the baryon loading and the density parameters of the GRB blast wave model, we repeated the calculation of Adriani et al. (2025) by including the total exposures from IC-EHE and Auger. The total number of expected events of a single flavour from the diffuse GRB flux for the combined three experiments is

$$\begin{array}{c}\hfill n_{\exp }(E_k,n_0;\theta )=\frac{4\pi }{3}{\displaystyle {\int }_{E_{\min }}^{E_{\max }}}{\mathcal{E}}^{\mathrm{tot}}{\mathrm{\Phi }}_{\nu }^{\mathrm{tot}}(E_{\nu },E_k,n_0;\theta )d{E}_{\nu },\end{array}$$(7)

where ℰ^tot = ∑_dT^dA_eff^d is the summed total exposure of all three experiments. The term T denotes the total lifetime associated with the selection, and A_eff is the corresponding effective area. The running index is d∈{KM3NeT, IC-EHE, Auger}, and the factor 1/3 corrects for considering all three neutrino flavours both in the diffuse flux and in the effective areas used. We set the integration limits as E ∈ [1 GeV, 100 EeV]. The KM3NeT exposure ℰ^KM3NeT is associated with the bright track selection reported in Aiello et al. (2025). The corresponding lifetime is 335 days, associated with the 19 and 21 detection unit KM3NeT/ARCA detector configuration, out of the foreseen 230 days (Aiello et al. 2024), and the effective area is all-flavour sky-averaged and averaged between neutrinos and anti-neutrinos. The IC-EHE exposure ℰ^{IC − EHE} is extracted from the 9-year analysis by Aartsen et al. (2018a) with a lifetime of 3145.5 days. The Pierre-Auger sky-averaged exposure, ℰ^Auger, is computed in Adriani et al. (2025) by considering the Earth-skimming, low-zenith downward-going, and high-zenith downward-going samples with the effective area from the data release (Aab et al. 2019). The corresponding lifetime is 6574.5 days, taken from the data release (Halim et al. 1488.). No neutrino events have been reported above tens of PeV in either the IC-EHE or Auger analyses. As we are using publicly available data, the exact treatment of systematics and the potential background of these analyses have been omitted.

We constructed the likelihood by considering the simple Poisson probability of observing one event in any of the three experiments given the predicted number of events from the considered GRB model:

$$\begin{array}{c}\hfill \mathcal{L}(E_k,n_0;\theta )=\mathrm{Poisson}(n_{\mathrm{obs}};n_{\exp }(E_k,n_0;\theta )).\end{array}$$(8)

We performed a 2D scan by varying the two parameters of interest, f_bη ∈ [10⁻², 10^1.7], which corresponds to varying the ratio of kinetic energy to gamma-ray luminosity through Eq. (1) and n₀ ∈ [1, 100] while keeping all other model parameters fixed. We then calculated the likelihood as defined in Eq. (8) at each point in this parameter space.

To constrain the two scan parameters, we adopted the Bayesian interpretation of probability and defined the posterior probability density as

$$\begin{array}{c}\hfill P(f_b\eta ,n_0;\theta )=P(E_k/L_{\gamma },n_0;\theta )=a_n\mathcal{L}(E_k/L_{\gamma },n_0;\theta ),\end{array}$$(9)

where we introduced the normalisation constant a_n such that the full posterior parameter space integrates to unity. We assumed uniform priors on both varying parameters.

The corresponding 1D marginalised probability distributions can be used to extract the best-fit and confidence intervals. Our results do not indicate well-defined constraints on the parameters. This is caused by the high degeneracy in the posterior distribution, which can be seen from Fig. 1 and the intrinsic uncertainty in the convergence of the nested sampling algorithm. Since the expected number of events increases with both parameters of interest (f_bη, n₀), we calculated the conditional posterior probability density as

$$\begin{array}{c}\hfill p(f_b\eta |n_0)=\frac{p(f_b\eta |n_0)}{{\displaystyle \int p}(f_b\eta ,n_0)d{f}_b\eta }.\end{array}$$(10)

Fig. 1. Posterior probability density of the 2D scan in fbη and n0 for the ISM model. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.

The resulting conditional probability density functions are shown in Fig. 2.

Fig. 2. Conditional posterior probability density for fbη calculated by fixing the density parameter to n0 = 1 cm−3 (A* = 0.2) for the ISM (WIND) model. The solid lines show the distribution of the conditional posterior, and the dark (light) shaded area shows the 68% (90%) confidence region. The ISM blast wave model is indicated by the teal colour and corresponding lower and left axes. The WIND blast wave model is shown in purple with the right and upper axes. For the ISM model, the best fit and 90% confidence level is $f_b\eta =5.4_{-4.4}^{+4.7}$. For the WIND model, it is $f_b\eta =29.5_{-28}^{+609}$.

4. Results and discussion

In order to calculate the diffuse flux from the population of GRB blast waves while varying both the baryon loading and the density of the surrounding medium, we fixed the remaining model parameters in θ. This procedure is described in Appendix A, and the corresponding values are listed in Table 1.

The Γ₀ parameter was fixed to a constant, as it determines the deceleration time (t_dec) of the model, which serves as the lower-limit of the time-integration over flux in Eq. (A.27). The parameter ϵ_e = 0.1 and ϵ_B = 0.1 − 0.01 are standard choices in GRB afterglow modelling (see e.g. Kumar & Zhang 2015; Miceli & Nava 2022) based on particle-in-cell simulations (Sironi & Spitkovsky 2011; Sironi et al. 2013). Additionally, the value of ϵ_B was chosen after performing a parameter space scan in f_bη versus ϵ_B similar to that described above while fixing n₀ = 1 cm⁻³. The resulting posterior is shown in Fig. B.1. This yielded a best-fit value of ϵ_B that is degenerate with f_bη, but it agrees well with broad-band modelling of very-high-energy GRBs (Barnard et al. 2025). Using the conditional posterior from Eq. (10), we estimated the best-fit value of f_bη that corresponds to one observed event in the combined exposure of KM3NeT, IceCube-EHE, and Auger. The best fit with its associated 90% confidence level is found to be $f_b\eta =5.4_{-4.4}^{+4.7}$ for fixed n₀ = 1 cm⁻³. These values were extracted from the conditional posteriors show in Fig. 2, where we observed a clear highest posterior density point defining the best-fit value and a smooth distribution covering the 68% and 90% credible interval. The diffuse neutrino flux, assuming these best-fit parameters, and the corresponding 68% containment region, is shown in Fig. 3. Fixing the efficiency parameter to η = 0.2, we constrained the baryon loading to be f_b ≤ 51 at a 90% confidence level for the ISM model.

Fig. 3. Energy-squared per-flavour diffuse astrophysical neutrino flux assuming (νe : νμ : ντ = 1 : 1 : 1) flavour equipartition. The teal (purple) dashed line corresponds to the diffuse flux from the ISM (WIND) blast wave models described in the text, with the corresponding 68% confidence interval. The reported KM3-230213A flux (Aiello et al. 2025) and corresponding 90% neutrino energy range is indicated by the grey cross. The joint fit flux considering non-observation in the IceCube-EHE and Auger samples in the same energy range (Adriani et al. 2025) is shown by the blue cross. The 68% confidence level contours from the IceCube NST (Abbasi et al. 2022b) and HESE (Abbasi et al. 2021) diffuse flux analyses are shown with the magenta and purple contours, respectively. The corresponding segmented fit analyses are shown by the magenta and purple crosses, and the IceCube Glashow resonance event (Aartsen et al. 2021b) is shown with an orange cross. The dotted lines show the upper limits from the ANTARES (Albert et al. 2024), IceCube-EHE (Aartsen et al. 2018a), and Auger (Halim et al. 1488.) analyses.

For the WIND model, we find the best-fit value and 90% confidence interval to be $f_b\eta =29.5_{-28}^{+609}$ for A_* = 0.2. This corresponds to a baryon loading of f_b ≈ 148. Conversely, when fixing f_b = 50, we found $A_{\ast }=0.{17}_{-0.10}^{+0.08}$ at a 90% confidence level. For this model, the density parameter, A_*, has a significantly larger impact on the diffuse flux than the baryon loading. We also note that increasing the two parameters simultaneously does not necessarily give a higher flux. The GRB blast wave in the WIND model loses significant amounts of energy through adiabatic expansion before the deceleration time, where UHE neutrino production proceeds photo-hadronically, and thus the flux is reduced. If the density of the surrounding material is too high or the baryon loading is significant, too much of the kinetic energy in the blast wave is lost in the initial interactions with the surrounding medium. As a consequence, the remaining shock energy will not be sufficient to accelerate the protons to UHE, reducing the subsequent UHE neutrino flux. This is not the case for the ISM mode, where the UHE neutrino fluence increases with f_b.

From the analytical expressions in Appendix A, we observed that for the ISM model, the time-dependent bulk Lorentz factor Γ is independent of Γ₀. Thus, the initial bulk Lorentz factor only enters the expression for the deceleration time; varying the value of Γ₀ has a negligible impact on the total UHE neutrino flux (see Fig. B.2). For the WIND model, the Γ₀ parameter does influence the UHE neutrino flux.

In both of the GRB blast wave models, we only considered neutrinos produced photo-hadronically by the interaction of accelerated UHE protons in the jet with photons from the afterglow. These pγ interactions dominate over the hadro-nuclear (pp, pn) interactions due to the significantly diluted particle density in the environment surrounding the GRB progenitor (Waxman & Bahcall 1997; Aartsen et al. 2016). Moreover, we only considered the primary single-pion production channel in the pγ interactions. The authors of references Murase (2007) and Razzaque (2013) have compared the full pion-production cross-section while considering higher multiplicities, and they found that the contribution is small below E_p ∼ 10¹⁸ eV.

In this analysis, we chose to fit the product of baryon loading and the efficiency parameter η rather than the baryon loading directly. They are related by Eq. (1). Thus, our results are interpretable for other values of the efficiency parameter η.

Previous analyses constrained the baryon loading from non-observations of sub-PeV neutrinos (see e.g. Aartsen et al. 2016). These analyses primarily focused on the neutrino emission from the GRB prompt phase and not from the subsequent afterglow. In the models investigated by Aartsen et al. (2016), the baryon loading factor f_b varies as a function of the Lorentz factor Γ₀. In the internal shock model considered therein, f_b is found to be ≲10 at a 90% confidence level, which is more constraining than our result. However, this is only for Γ₀ ≲ 300. Since the considered ISM model does not depend strongly on Γ₀ (see Appendix B), when comparing with the value Γ₀ = 10^2.8 used for our calculations, the 90% confidence level from Aartsen et al. (2016) is f_b ≲ 200. For the internal collision-induced magnetic reconnection and turbulence model considered therein (Zhang & Yan 2011), the constraints are consistent with our findings also for lower Γ₀ values. The photospheric model, which is severely constrained in Aartsen et al. (2016), is disfavoured by the non-observation of a dominant thermal component in the GRB spectra (Goldstein et al. 2012). Finally, we note that the baryon loading factor may be different in the initial prompt emission phase of the GRBs and the afterglow. Although the fraction of baryons in the prompt phase will contribute to the particle spectra in the afterglow emission, interactions with the surrounding environment and hadro-nuclear interactions in the prompt phase itself can alter the baryon-to-electron energy ratio during the evolution of the jet. Furthermore, if the GRB jet is dominated by the magnetic field (Lyutikov & Blandford 2003), there will be fewer baryons to produce neutrinos in the prompt phase. Constraints on f_b from the prompt (afterglow) phase may therefore not directly relate to the afterglow (prompt) emission.

Abbasi et al. (2022c) found that the prompt emission of GRBs observed in the electromagnetic spectrum contributes to less than 1% of the observed diffuse neutrino flux below 10 PeV and that emission from the afterglow on timescales up to 10⁴ s is ≲24%. Although our analysis focuses on the ≳10 PeV energy range associated with KM3-230213A, we see from Fig. 3 that the diffuse neutrino flux below 10 PeV for both models is well below that of the NST and HESE IceCube fits. Our results are consistent with Abbasi et al. (2022c) regarding the contribution of GRBs to the diffuse neutrino flux at lower energies and allow for the interactions of lGRB blast waves to explain the UHE diffuse neutrino flux.

Although our focus has been exclusively on GRBs, other classes of astrophysical objects, such as active galactic nuclei (AGNs) may also be significant contributors to the diffuse high-energy neutrino flux. Luminous AGNs, such as blazars, are expected to accelerate protons to UHEs in highly collimated jets, resulting in high-energy neutrino production when these protons interact with the varying radiation fields. As the predominant neutrino production channel is through pion decay, significant emission of high-energy gamma-rays from π⁰ decays will accompany neutrinos from such non-transient sources as blazars. The measurement of the diffuse gamma-ray sky by Fermi-LAT (Atwood et al. 2009) places strong constraints on the contribution of steady sources to the diffuse UHE neutrino flux. Since GRBs are transient sources, they are not affected by these constraints. Other types of transient sources are also candidates for UHE neutrino production, such as tidal disruption events (e.g. Lunardini & Winter 2017) and magnetar-powered super-luminous supernovae (e.g. Fang et al. 2019). However, the current sparsity of their detection in gamma rays makes the population and cosmic evolution of these transients less known and their contribution to the diffuse UHE neutrino flux uncertain (e.g. Das et al. 2025).

5. Conclusions

We have calculated the contribution of different GRB blast wave models to the diffuse neutrino flux and constrained the lGRB model parameters with respect to the diffuse UHE neutrino flux most likely associated with KM3-230213A (Aiello et al. 2025). The total neutrino fluence from individual lGRBs was calculated following the blast wave models presented by Razzaque (2013) and the lGRB luminosity function used was fitted to Swift and Fermi-GBM observations by Banerjee et al. (2021). The total diffuse flux was calculated using the nested sampling algorithm implemented in the UltraNest python package, allowing the variation of all integration model parameters, which takes into consideration a broader population of lGRBs up to z = 5 following the fitted luminosity function.

We considered two different models for UHE neutrino production from GRB blast waves: one in which the density of the surrounding matter remains constant around the GRB progenitor and another in which the density decreases radially. For the GRB blast wave model with constant ISM density n₀ = 1 cm⁻³, the baryon loading is constrained to be f_b ≤ 51 at 90% confidence. Assuming a larger value for the ISM density, f_b is significantly more constrained. The corresponding best-fit baryon loading, f_b = 27, would produce a diffuse UHE neutrino flux consistent with the detection of one UHE event (i.e. KM3-230213A) within the cumulative exposure of KM3NeT, IceCube-EHE and Auger. When assuming the GRB blast wave interacts with a wind-type medium with radially decreasing density, the baryon loading is less well defined, as the density parameter (A_*) dominates. It is constrained to be f_b ≤ 1065 at 90% confidence for A_* = 0.1. By fixing the baryon loading factor to f_b ∼ 50, i.e. typical values from prompt emission models (Aartsen et al. 2016), we found the best-fit value of the density parameter $A_{\ast }=0.{017}_{-0.10}^{+0.08}$ at 90% confidence.

Our results show that a large population of lGRBs can give rise to the diffuse UHE neutrino flux associated with KM3-230213A. Moreover, both GRB models we considered are shown to be consistent with existing limits on their contribution to the diffuse neutrino flux at lower energies. Although the true diffuse neutrino flux at UHEs may come from additional sources (see Mészáros 2017), GRB blast waves can contribute significantly to the UHE neutrino flux required for KM3-230213A whilst remaining consistent with previous limits on GRB model parameters (Aartsen et al. 2016). Future observations by upcoming large radio detectors such as GRAND (Álvarez-Muñiz et al. 2020), Askaryan detectors such as RNO-G (Aguilar et al. 2021), and combined Askaryan and Cherenkov detectors such as IceCube-Gen2 (Aartsen et al. 2021a) can contribute to better characterisation of the UHE neutrino flux. Our modelling of the UHE neutrino flux corresponding to KM3-230213A motivates further observations of the electromagnetic sky and exploration of the role of GRBs as multi-messenger sources.

Acknowledgments

We thank the anonymous reviewers for their insightful comments and constructive suggestions, which greatly improved this manuscript. The authors acknowledge the financial support of: KM3NeT-INFRADEV2 project, funded by the European Union Horizon Europe Research and Innovation Programme under grant agreement No 101079679; Funds for Scientific Research (FRS-FNRS), Francqui foundation, BAEF foundation. Czech Science Foundation (GAČR 24-12702S); Agence Nationale de la Recherche (contract ANR-15-CE31-0020), Centre National de la Recherche Scientifique (CNRS), Commission Européenne (FEDER fund and Marie Curie Program), LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001), Paris Île-de-France Region, Normandy Region (Alpha, Blue-waves and Neptune), France, The Provence-Alpes-Côte d’Azur Delegation for Research and Innovation (DRARI), the Provence-Alpes-Côte d’Azur region, the Bouches-du-Rhône Departmental Council, the Metropolis of Aix-Marseille Provence and the City of Marseille through the CPER 2021-2027 NEUMED project, The CNRS Institut National de Physique Nucléaire et de Physique des Particules (IN2P3); Shota Rustaveli National Science Foundation of Georgia (SRNSFG, FR-22-13708), Georgia; This work is part of the MuSES project which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (grant agreement No 101142396). The General Secretariat of Research and Innovation (GSRI), Greece; Istituto Nazionale di Fisica Nucleare (INFN) and Ministero dell’Università e della Ricerca (MUR), through PRIN 2022 program (Grant PANTHEON 2022E2J4RK, Next Generation EU) and PON R&I program (Avviso n. 424 del 28 febbraio 2018, Progetto NRRP), Italy; IDMAR project Po-Fesr Sicilian Region az. 1.5.1; A. De Benedittis, W. Idrissi Ibnsalih, M. Bendahman, A. Nayerhoda, G. Papalashvili, I. C. Rea, A. Simonelli have been supported by the Italian Ministero dell’Università e della Ricerca (MUR), Progetto CIR01 00021 (Avviso n. 2595 del 24 dicembre 2019); KM3NeT4RR MUR Project National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 3.1, Funded by the European Union – NextGenerationEU,CUP I57G21000040001, Concession Decree MUR No. n. Prot. 123 del 21/06/2022; Ministry of Higher Education, Scientific Research and Innovation, Morocco, and the Arab Fund for Economic and Social Development, Kuwait; Nederlandse organisatie voor Wetenschappelijk Onderzoek (NWO), the Netherlands; The grant “AstroCeNT: Particle Astrophysics Science and Technology Centre”, carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund; The program: “Excellence initiative-research university” for the AGH University in Krakow; The ARTIQ project: UMO-2021/01/2/ST6/00004 and ARTIQ/0004/2021; Ministry of Education and Scientific Research, Romania; Slovak Research and Development Agency under Contract No. APVV-22-0413; Ministry of Education, Research, Development and Youth of the Slovak Republic; MCIN for PID2021-124591NB-C41, -C42, -C43 and PDC2023-145913-I00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, for ASFAE/2022/014 and ASFAE/2022/023 with funding from the EU NextGenerationEU (PRTR-C17.I01) and Generalitat Valenciana, for Grant AST22_6.2 with funding from Consejería de Universidad, Investigación e Innovación and Gobierno de España and European Union – NextGenerationEU, for CSIC-INFRA23013 and for CNS2023-144099, Generalitat Valenciana for CIDEGENT/2020/049, CIDEGENT/2021/23, CIDEIG/2023/20, ESGENT2024/24, CIPROM/2023/51, GRISOLIAP/2021/192 and INNVA1/2024/110 (IVACE+i), Spain; Khalifa University internal grants (ESIG-2023-008, RIG-2023-070 and RIG-2024-047), United Arab Emirates; The European Union’s Horizon 2020 Research and Innovation Programme (ChETEC-INFRA – Project no. 101008324).

References

Appendix A: Individual GRB neutrino fluence

We provide a more detailed analytical summary of the UHE neutrino flux produced in a single GRB blast wave and the total energy released in UHE neutrinos. The deceleration timescale is

$$\begin{array}{c}\hfill t_{\mathrm{dec}}(z,E_k)={\left[\frac{3E_k{(1+z)}^3}{64\pi n_0{m}_p{c}^5{\mathrm{\Gamma }}_0^8}\right]}^{1/3},\end{array}$$(A.1)

where E_k is the kinetic energy in the blast wave, n₀ is the number density of the ISM, m_p is the proton mass, c is the speed of light, and Γ₀ is the initial bulk Lorentz factor of the outflow. The bulk Lorentz factor after the deceleration timescale evolves in the constant density ISM as

$$\begin{array}{c}\hfill \mathrm{\Gamma }(t)={\mathrm{\Gamma }}_0{\left(\frac{t_{\mathrm{dec}}}{4t}\right)}^{3/8},\end{array}$$(A.2)

and the radius and magnetic field strength of the blast wave increases correspondingly as

$$\begin{array}{c}\hfill R(t)=\frac{2{\mathrm{\Gamma }}^2(t)act}{1+z},\end{array}$$(A.3)

and

$$\begin{array}{c}\hfill B^{\prime }(t)={[32\pi {ϵ}_B{n}_0{m}_p{c}^2]}^{1/2}\mathrm{\Gamma }(t),\end{array}$$(A.4)

respectively. The constant a = 4; ϵ_B denotes the fraction of the forward-shock energy from the blast wave converted into magnetic energy; and the timescale is t > t_dec. In the WIND model, we instead consider the deceleration timescale

$$\begin{array}{c}\hfill t_{\mathrm{dec}}(z,E_k)=\frac{E_k(1+z)}{16\pi A{m}_p{c}^3{\mathrm{\Gamma }}_0^4},\end{array}$$(A.5)

and the bulk Lorentz factor as

$$\begin{array}{c}\hfill \mathrm{\Gamma }(t)={\mathrm{\Gamma }}_0{\left(\frac{t_{\mathrm{dec}}}{4t}\right)}^{1/4},\end{array}$$(A.6)

with A = 3.02 × 10³⁵A_* cm⁻¹. Consequently, the density of the surrounding medium in the WIND model is n(R) = AR⁻², whereas it is constant for the ISM model.

Within this blast wave model, electrons accelerated in the shock exhibit distinct behaviours in different regimes characterised by different electron Lorentz factors

$$\begin{array}{cc}& {\gamma \prime }_m={ϵ}_e\frac{m_p}{m_e}\mathrm{\Gamma }(t),\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill & {\gamma \prime }_c=\frac{6\pi m_{e}c(1+z)}{{\sigma }_{T}tB\prime {}^2(t)\mathrm{\Gamma }(t)},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill & {\gamma \prime }_s={\left[\frac{6\pi e}{{\sigma }_{T}B\prime (t)\varphi }\right]}^{1/2},\hfill \end{array}$$(A.9)

corresponding to minimum, cooling, and saturation electron states, respectively. The primed notation indicates the blast wave comoving frame. The ϵ_e parameter describes the kinetic energy in the blast wave converted into random electron energy, m_e is the electron mass, σ_T is the Thomson cross section, and ϕ is the number of gyroradii for the accelerated electrons. These electrons produce a flux of synchrotron photons as

$$\begin{array}{c}\hfill F_{\nu }=\frac{N_e}{4\pi d_L^2}\frac{P(\gamma \prime )}{h\nu }\frac{{\mathrm{\Gamma }}^2(t)}{{(1+z)}^2},\end{array}$$(A.10)

where N_e = (4/3)πR³(t)n(t) is the total number of electrons in the blast wave and

$$\begin{array}{c}\hfill d_L(z)=\frac{c(1+z)}{H_0}{\displaystyle {\int }_0^z}\frac{d{z}^{\ast }}{\sqrt{{\mathrm{\Omega }}_m{(1+z^{\ast })}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}}},\end{array}$$(A.11)

is the luminosity distance. The synchrotron power of electrons with Lorentz factor γ′ is given by

$$\begin{array}{c}\hfill P(\gamma \prime )=\frac{c{\sigma }_T}{6\pi }B\prime {}^2\gamma \prime {}^2,\end{array}$$(A.12)

and the characteristic synchrotron frequency is

$$\begin{array}{c}\hfill h\nu =\frac{3}{2}\frac{B\prime (t)}{B_Q}\gamma \prime {}^2(t)m_e{c}^2\frac{\mathrm{\Gamma }(t)}{1+z},\end{array}$$(A.13)

with normalising magnetic field strength B_Q = 4.41 × 10¹³ G.

A.1. pγ interactions

The density of observed synchrotron photons is parametrised in the different regimes and calculated as

$$\begin{array}{c}\hfill n{\prime }_{\gamma }(ϵ\prime )=\frac{2d_L^2(1+z)F_{\nu ,m}}{r^{2}c\mathrm{\Gamma }ϵ\prime }\times \{\begin{array}{c}{(ϵ{\prime }_c/ϵ{\prime }_m)}^{-3/2}{(ϵ\prime /ϵ{\prime }_c)}^{-2/3}\mathrm{for}ϵ{\prime }_a<ϵ\prime <ϵ{\prime }_c\hfill \\ {(ϵ\prime /ϵ{\prime }_m)}^{-3/2}\mathrm{for}ϵ{\prime }_c<ϵ\prime <ϵ{\prime }_m\hfill \\ {(ϵ\prime /ϵ{\prime }_m)}^{-k/2-1}\mathrm{for}{ϵ}_m\prime <ϵ\prime <ϵ{\prime }_s,\hfill \end{array}\end{array}$$(A.14)

where k is the spectral index of the electron distribution, m, c, s correspond to the different regimes, and ϵ′=(hν)′ = hν(1 + z)/Γ. The rate of pγ interactions for protons with Lorentz factor γ′_p in the blast wave with this synchrotron photon field is

$$\begin{array}{c}\hfill K_{p\gamma }(\gamma {\prime }_p)=\frac{c}{2\gamma \prime {{}^2}_p}{\displaystyle {\int }_{ϵ{\prime }_{\mathrm{th}}}^{\infty }}dϵ{\prime }_rϵ{\prime }_r{\sigma }_{p\gamma }(ϵ{\prime }_r){\displaystyle {\int }_{ϵ{\prime }_r/(2\gamma {\prime }_p)}^{\infty }}dϵ\prime \frac{n{\prime }_{\gamma }(ϵ\prime )}{ϵ\prime {}^2},\end{array}$$(A.15)

where ϵ′_r is the photon energy in the rest frame of the proton and ϵ′_th = m_πc² + m_π²c²/2m_p is the threshold for pion production. The cross-section for pγ → nπ⁺ interactions is

$$\begin{array}{c}\hfill {\sigma }_{p\gamma }(ϵ{\prime }_r)={\sigma }_0{\mathrm{\Gamma }}_{\mathrm{\Delta }}^2{s}^2ϵ\prime {{}^{-2}}_r{[{\mathrm{\Gamma }}_{\mathrm{\Delta }}^{2}s+{(s-m_{\mathrm{\Delta }}^2)}^2]}^{-1},\end{array}$$(A.16)

with s = m_p²c⁴ + 2ϵ′_rm_pc², σ₀ = 3.11 × 10⁻²⁹ cm², and Γ_Δ = 0.11 GeV is the width of the Delta resonance. Finally, we calculate the optical depth for pγ interactions as

$$\begin{array}{c}\hfill {\tau }_{p\gamma }({\gamma }_p^{\prime })=K_{p\gamma }({\gamma }_p^{\prime })\frac{R(t)}{2ac\mathrm{\Gamma }}=K_{p\gamma }({\gamma }_p^{\prime })\frac{t\mathrm{\Gamma }}{1+z},\end{array}$$(A.17)

where a is the same constant as in Eq. (A.3).

A.2. Proton population

The majority of the multi-messenger emissivity from an individual GRB blast wave is driven by the interactions of protons accelerated by the shocks in the blast wave. In this model, the differential number density of protons as a function of proton energy is given by

$$\begin{array}{c}\hfill n(E_p)=\frac{{\mathcal{E}}_{\mathrm{CR}}}{V{E}_p^{2}ln({\gamma }_{p,s}^{\prime }/\mathrm{\Gamma })},\end{array}$$(A.18)

where γ′_p, s is the saturation proton Lorentz factor given by

$$\begin{array}{c}\hfill {\gamma }_{p,s}^{\prime }=\frac{e{B}^{\prime }(t)}{\varphi m_{p}c}\frac{t\mathrm{\Gamma }}{1+z}.\end{array}$$(A.19)

Here, ϕ is the number of gyroradii for the protons. The energy of the cosmic rays in the blast wave after the deceleration timescale t_dec is

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{CR}}=\frac{4}{3}\pi {ϵ}_p{n}_0{r}^3(t)m_p{c}^2[{\mathrm{\Gamma }}^2-1],\end{array}$$(A.20)

with ϵ_p being the fraction of the blast wave energy that goes into proton acceleration. If this population of blast wave protons could escape their source environment as cosmic rays, their energy would be

$$\begin{array}{c}\hfill E_{p,s}=\frac{m_p{c}^2{\gamma }_{p,s}^{\prime }\mathrm{\Gamma }}{1+z}.\end{array}$$(A.21)

To avoid having a population of runaway high-energy protons, we introduce an exponential cut-off term ∝exp(−E_p/E_p, s) to suppress the proton spectrum after saturation. Including this cut-off, if the protons could freely escape the blast wave environment and make it to Earth, their flux would be

$$\begin{array}{c}\hfill J_p(E_p)=\frac{c}{4\pi }{\left(\frac{R}{d_L}\right)}^{2}n(E_p)e^{-E_p/E_{p,s}}.\end{array}$$(A.22)

A.3. Neutrino flux

The protons accelerated within the GRB blastwave shocks will interact with the synchrotron photons and produce neutrinos through pion production and consequent decays

$$\begin{array}{cc}& {\pi }^{+}\to {\mu }^{+}+{\nu }_{\mu },\hfill \\ \hfill & {\mu }^{+}\to e^{+}+{\nu }_e+{\overline{\nu }}_{\mu },\hfill \end{array}$$(A.23)

resulting in the 1 : 2 : 0 flavour ratio associated with photopion interactions. Although neutrinos are also produced through neutron beta decay, this process is subdominant and consequently ignored in our calculations. We assume an equal probability of producing π⁺ and π⁰ with mean inelasticity ⟨x⟩≈0.2, resulting in a pion flux

$$\begin{array}{c}\hfill J_{\pi }(E_{\pi })\approx \frac{1}{⟨x⟩}{J}_p\left(\frac{E_{\pi }}{⟨x⟩}\right){\tau }_{p\gamma }\left(\frac{E_{\pi }(1+z)}{m_p{c}^2⟨x⟩\mathrm{\Gamma }}\right).\end{array}$$(A.24)

The consequent flux of muons is found by integrating over the kinematics of the process governed by the mass ratio as

$$\begin{array}{c}\hfill J_{\mu }(E_{\mu })={\displaystyle {\int }_0^1}\frac{\mathit{dx}}{x}{f}_{{\pi }^{+}\to {\mu }^{+}}(x)J_{\pi }\left(\frac{E_{\mu }}{x}\right),\end{array}$$(A.25)

where x = E_ν/E_π. The resulting neutrino fluxes from all relevant decay channels are calculated by considering similar mass ratios and integrating over the kinematics of the interactions

$$\begin{array}{cc}& J_{{\nu }_{\mu }}(E_{\nu })={\displaystyle {\int }_0^1}\frac{\mathit{dx}}{x}{f}_{{\pi }^{+}\to {\nu }_{\mu }}(x)J_{\pi }\left(\frac{E_{\nu }}{x}\right);x=\frac{E_{\nu }}{E_{\pi }};\hfill \end{array}$$(A.26)

$$\begin{array}{cc}& J_{{\nu }_e}(E_{\nu })={\displaystyle {\int }_0^1}\frac{\mathit{dy}}{y}{\displaystyle {\int }_0^1}\frac{\mathit{dx}}{x}{f}_{{\mu }^{+}\to {\nu }_e}(x,y)f_{\pi \to \mu }(x)J_{\pi }\left(\frac{E_{\nu }}{\mathit{xy}}\right);\hfill \\ \hfill & J_{{\overline{\nu }}_{\mu }}(E_{\nu })={\displaystyle {\int }_0^1}\frac{\mathit{dy}}{y}{\displaystyle {\int }_0^1}\frac{\mathit{dx}}{x}{f}_{{\mu }^{+}\to {\overline{\nu }}_{\mu }}(x,y)f_{\pi \to \mu }(x)J_{\pi }\left(\frac{E_{\nu }}{\mathit{xy}}\right),\hfill \end{array}$$

where x = E_μ/E_π and y = E_ν/E_μ in the last two expressions. The scaling functions f are parametrisations of the decays as derived in Lipari (1993). The final step is to integrate the total neutrino flux from a single GRB blastwave over its duration. The emission happens only after the deceleration time t_dec and decreases as time increases,

$$\begin{array}{c}\hfill S_{\nu }(E_{\nu })={\displaystyle {\int }_{t_{\mathrm{dec}}}^{\infty }}{J}_{\nu }(E_{\nu })dt.\end{array}$$(A.27)

This is the fluence of an individual GRB blastwave and gives the total energy released in neutrinos from photopion interactions.

Appendix B: Dependence on other model parameters

In the main text, we focused on what happens to the predicted neutrino flux as we vary the baryon loading and the density of the surrounding medium, while fixing the remaining model parameters. It is also worth investigating how changing the other parameters influences the UHE neutrino flux. We see from Eq. (1) that modifying the efficiency of energy conversion from kinetic energy to gamma-ray luminosity η is directly proportional to changing the ratio between kinetic energy and gamma-ray luminosity E_k/L_γ itself. Thus, we chose to constrain the product of the baryon loading and this efficiency parameter in the main text.

Since previous searches and literature focus on fitting the bulk Lorentz factor Γ, we also consider how this parameter influences our model and flux prediction. From the analytical expressions in Appendix A, we see that the model depends on the initial bulk Lorentz factor Γ₀, which determines the deceleration timescale t_dec and consequently the bulk Lorentz factor Γ. For the ISM model, the time-dependent bulk Lorentz factor Γ(t) is independent of the initial bulk Lorentz factor Γ₀, and, since the model GRB flux only depends on Γ(t) and not on t_dec, the effects of varying Γ₀ are minimal. The deceleration timescale t_dec only enters the calculation via the lower integration bound of the single-GRB fluence of Eq. (A.27), and the impact on the diffuse UHE neutrino flux is small. This can be seen in Fig. B.2, where we fixed the value of the ISM density to n₀ = 1 cm⁻³, and fitted instead for varying Γ₀. The same relation holds for the WIND model. In this case, the deceleration time t_dec depends more strongly on the initial bulk Lorentz factor, t_dec ∝ Γ₀⁻⁴, and the time elapsed before UHE neutrino production is shorter. However, as this dependence only enters in Eq. (A.27), it does not significantly affect the diffuse UHE neutrino flux, especially considering the competing effect of reduction in the UHE neutrino flux with increasing f_b as described in the main text.

A similar fit was performed for the magnetic energy parameter ϵ_B. In this scenario, we once again fixed n₀ = 1 cm⁻³ and performed a scan for f_bη and ϵ_B. The resulting posterior probability density is shown in Fig. B.1. We see from this parameter scan that the variables are quite degenerate, so choosing a different value for ϵ_B than the one used in the main analysis will significantly change the constraints on f_b.

Fig. B.1. Posterior probability density of the 2D scan in fbη and ϵB for the ISM model with n0 = 1 cm−3. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.

Fig. B.2. Posterior probability density of the 2D scan in fbη and Γ0 for the ISM model with n0 = 1 cm−3. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.

All Tables

All Figures

	Fig. 1. Posterior probability density of the 2D scan in fbη and n0 for the ISM model. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.
In the text

Fig. 2. Conditional posterior probability density for fbη calculated by fixing the density parameter to n0 = 1 cm−3 (A* = 0.2) for the ISM (WIND) model. The solid lines show the distribution of the conditional posterior, and the dark (light) shaded area shows the 68% (90%) confidence region. The ISM blast wave model is indicated by the teal colour and corresponding lower and left axes. The WIND blast wave model is shown in purple with the right and upper axes. For the ISM model, the best fit and 90% confidence level is $f_b\eta =5.4_{-4.4}^{+4.7}$. For the WIND model, it is $f_b\eta =29.5_{-28}^{+609}$.

In the text

Fig. 3. Energy-squared per-flavour diffuse astrophysical neutrino flux assuming (νe : νμ : ντ = 1 : 1 : 1) flavour equipartition. The teal (purple) dashed line corresponds to the diffuse flux from the ISM (WIND) blast wave models described in the text, with the corresponding 68% confidence interval. The reported KM3-230213A flux (Aiello et al. 2025) and corresponding 90% neutrino energy range is indicated by the grey cross. The joint fit flux considering non-observation in the IceCube-EHE and Auger samples in the same energy range (Adriani et al. 2025) is shown by the blue cross. The 68% confidence level contours from the IceCube NST (Abbasi et al. 2022b) and HESE (Abbasi et al. 2021) diffuse flux analyses are shown with the magenta and purple contours, respectively. The corresponding segmented fit analyses are shown by the magenta and purple crosses, and the IceCube Glashow resonance event (Aartsen et al. 2021b) is shown with an orange cross. The dotted lines show the upper limits from the ANTARES (Albert et al. 2024), IceCube-EHE (Aartsen et al. 2018a), and Auger (Halim et al. 1488.) analyses.

In the text

	Fig. B.1. Posterior probability density of the 2D scan in fbη and ϵB for the ISM model with n0 = 1 cm−3. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.
In the text

	Fig. B.2. Posterior probability density of the 2D scan in fbη and Γ0 for the ISM model with n0 = 1 cm−3. The solid (dashed, dotted) lines show the 1σ (2σ, 3σ) contours.
In the text