© The Authors 2026

1. Introduction

IceCube has measured a diffuse flux of astrophysical neutrinos with a spectrum roughly consistent with a power law, Φ_ν ∝ E_ν^−γ, over the range ∼10 TeV–10 PeV, with γ = 2.4 − 2.9 depending on the specific analysis considered (e.g., Abbasi et al. 2021, 2022). Since its discovery in 2013 (IceCube Collaboration 2013), improved statistics and analysis techniques have refined this picture. The latest results extend the spectrum down to ∼1 TeV and have revealed evidence of a low-energy break near 10 TeV (Basu et al. 2025), indicating that the diffuse neutrino spectrum deviates from a simple power law. The origin of this unresolved emission remains unknown.

In addition to the measurement of the diffuse flux, searches for neutrino point sources have been performed over the past decade. To date, only two astrophysical neutrino sources have been identified at a high confidence level, and both are associated with active galaxies of different types: the blazar TXS 0506+065 (IceCube Collaboration 2018; Padovani et al. 2018) and the nearby Seyfert II NGC 1068, which stands as the “hottest” spot in the neutrino sky (Abbasi et al. 2025a,b).

The detection of neutrinos from NGC 1068, in combination with the absence of very high energy γ-rays, points toward neutrino production in a compact γ-ray opaque region, the corona, which is located in the vicinity of the central black hole. Three main theoretical scenarios have been proposed to model proton acceleration in active galactic nuclei (AGN) cores, namely, diffusive shock acceleration at shock waves (e.g., Stecker et al. 1991; Inoue et al. 2020), stochastic acceleration in magnetohydrodynamic turbulence (e.g., Murase et al. 2020; Fiorillo et al. 2024; Lemoine & Rieger 2025; Saurenhaus et al. 2026), and systematic acceleration in magnetic reconnection (e.g., Fiorillo et al. 2024; Karavola et al. 2025). In a magnetically powered corona, the last two mechanisms are generally considered to accelerate protons to nonthermal energies. In turbulence-based models, neutrinos are produced by pγ and p-p interactions in a coronal region of ∼10 − 100 gravitational radii. Based on this scenario, Murase et al. (2020) and Padovani et al. (2024) determined the diffuse neutrino spectrum expected from these turbulent coronae according to the model of Murase et al. (2020). More recently, Fiorillo et al. (2025) calculated the neutrino diffuse flux according to the strongly magnetized turbulence scenario of Fiorillo et al. (2024). The general consensus is that turbulent AGN coronae may account for the ∼1–10 TeV diffuse neutrino emission but not the higher energy tail, due to strong Bethe-Heitler proton cooling.

In this paper, we calculate the diffuse neutrino flux arising from non-jetted AGNs based on the magnetic reconnection coronal model introduced by Fiorillo et al. (2024), Karavola et al. (2025). In this model, protons are accelerated in magnetospheric current sheets and on average obtain an energy of E_p ∼ σ_pm_pc² (Comisso 2024; Sironi et al. 2025; Hakobyan et al. 2025), where σ_p = B²/(4πn_pm_pc²)≫1 is the proton plasma magnetization; here, B is the magnetic field strength and n_p is the number density of non-relativistic protons in the corona. Relativistic protons interact with coronal X-ray photons, producing pions. Neutrinos with a typical energy of E_ν ∼ E_p/20 are subsequently produced through the decays of pions and muons. Our model is characterized by three key parameters: the proton plasma magnetization (σ_p), the coronal X-ray luminosity (L_X), and the coronal radius (R; or equivalently the black hole mass, M). We coupled this single-source model with an AGN mock catalog (Georgakakis et al. 2020) to compute the cumulative neutrino flux from AGN coronae up to redshift z = 4. We show that the observed diffuse neutrino flux in the 1 TeV–1 PeV range can be naturally reproduced if ∼10% of the AGN coronae have σ_p ∼ 10⁵, while the rest are distributed over a range of lower magnetizations¹.

2. Diffuse neutrino flux calculation

In this section we outline the corona model and introduce its key parameters. We then present the AGN mock catalog and describe the diffuse neutrino flux calculation.

2.1. Corona model

We assumed that neutrinos in non-jetted AGNs are produced within a compact, highly magnetized region close to the black hole coinciding with the site of nonthermal X-ray emission commonly referred to as the corona. We briefly outline the model introduced in Fiorillo et al. (2024) to explain the IceCube observations of NGC 1068 and explored more generally for AGN coronae in Karavola et al. (2025). In this model, the powerhouse of coronal emission is magnetic reconnection in current sheets that form within the magnetospheric region of an accreting black hole. The current sheet is described by a cuboid with dimensions L × L × β_recL, with $L=\overline{R}{r}_{\mathrm{g}}$, β_rec ∼ 0.1 being the reconnection rate and r_g = GM/c². For the scope of the numerical calculations, we describe the corona as a spherical environment with volume β_recL³. Under this description, one can define an effective coronal radius: R = [3β_rec/(4π)]^1/3L ∼ 0.29L. In this work, we adopted $\overline{R}=12$, motivated by global numerical simulations (Ripperda et al. 2022).

A fraction, η_X ∼ 0.5, of the dissipated magnetic energy is radiated as X-rays, with a bolometric luminosity (0.1–100 keV) given by L_X = η_Xcβ_recB²R², where B is the strength of the reconnecting magnetic field, and R is the characteristic size of the corona. A fraction, η_p ∼ 0.1, of the dissipated energy is transferred to relativistic protons with luminosity: L_p = (η_p/η_X)L_X. We discuss the specific choices of η_X and η_p in Sec. 4. The distribution function of accelerated protons is a broken power law, dN/dγ ∝ γ⁻¹ for 1 < γ < γ_br or γ^−s_p otherwise, with the post-break slope (s_p) depending on the guide-field² strength (Comisso 2024). Here we adopted s_p = 3, which corresponds to strong guide-field reconnection.

The break Lorentz factor is determined by the proton plasma magnetization, γ_br ∼ σ_p = B²/(4πn_pm_pc² (Comisso 2024; Hakobyan et al. 2025), where n_p is the plasma density of (non-relativistic) protons in the corona (upstream of the reconnection layer). We envisioned that the reconnection layer is formed in a region distinct from the dense accretion disk whenever steady accretion is temporarily halted due to the accumulation of magnetic flux. As a result, the proton density in this region is expected to be extremely low, allowing σ_p to reach high values³.

In our model, high-energy neutrinos are produced through interactions of relativistic protons with coronal X-rays. Karavola et al. (2025) demonstrated that the resulting neutrino emission is controlled by two parameters: the proton magnetization (σ_p) and the AGN X-ray Eddington ratio, defined as λ_X, Edd = L_X/L_Edd, where L_Edd = 4πGMm_pc/σ_T is the Eddington luminosity of the accreting black hole.

We also accounted for the synchrotron cooling of muons and pions (as described in Appendix A), which becomes important in the strong magnetic fields considered here, B ∼ 0.1 − 10³ kG. Synchrotron cooling introduces characteristic break energies to the neutrino spectrum, thus breaking the self-similarity of spectral shapes with respect to λ_X, Edd – see Appendix A – and softens the neutrino spectra, making the imprint of the proton slope less prominent, as shown in Appendix B.

We created neutrino spectral templates using the leptohadronic code ATHEνA (Dimitrakoudis et al. 2012) for a parameter grid spanning log₁₀L_X = 42 − 47, log₁₀(σ_p) = 3 − 5, and log₁₀(R) = 12.3 − 14.3 in integer steps. These ranges were motivated by the AGN mock catalogs described in the next section.

2.2. AGN mock catalog

To calculate the neutrino emission from a population of AGNs, the X-ray luminosity and the black hole mass of each source are required. For this purpose, we used the simulated AGN sample of Georgakakis et al. (2020), which provides the 2–10 keV X-ray luminosity⁴ and the stellar mass (M_*) of the host galaxy. This M_* is assigned to a black hole mass according to the redshift-independent relation M = 2 ⋅ 10⁻³M_* (Marconi & Hunt 2003), and it translates to a current sheet dimension, L = 12 r_g(M). The catalog was constructed using the X-ray luminosity function of Ueda et al. (2014) as described in Appendix C.

2.3. Diffuse flux calculation

We denote the template neutrino luminosity spectra on the discrete grid of proton magnetization (σ_p), the X-ray luminosity (L_X), and the coronal radius (R) as ${\mathcal{L}}_{(a,b,c)}\equiv E_{\mathrm{em}}L(E_{\mathrm{em}};{\sigma }_{\mathrm{p}}^{(a)},L_{\mathrm{X}}^{(b)},R^{(c)})$, where E_em is the neutrino energy in the rest frame of the galaxy, a = 1, ..., 3, b = 1, ..., 6, and c = 1, ..., 3. The neutrino flux spectrum of a simulated AGN corona with ${\sigma }_{\mathrm{p}}^{(i)}$, $L_{\mathrm{X}}^{(i)}$, and R⁽ⁱ⁾ at redshift z_i was then derived as

$$\begin{array}{c}\hfill E_{\mathrm{obs}}{F}^{(i)}(E_{\mathrm{obs}})=\frac{1}{4\pi D_{\mathrm{L}}^2(z_i)}{E}_{\mathrm{em}}L(E_{\mathrm{em}};{\sigma }_{\mathrm{p}}^{(i)},L_{\mathrm{X}}^{(i)},R^{(i)}),\end{array}$$(1)

where E_obs = E_em/(1 + z_i), D_L is the luminosity distance⁵ and ℐ[ℒ_{(a, b, c)}] is the interpolation operator. We performed a trilinear interpolation of the template spectra in logarithmic space, i.e., we interpolated log₁₀(E_emL(E_em)) as a function of log₁₀L_X, log₁₀R and log₁₀(σ_p), on a fixed logarithmic energy grid in the source rest frame (log₁₀E_em).

Fiorillo et al. (2024) and Karavola et al. (2025) showed that a value of σ_p ∼ 10⁵ was required to reproduce the IceCube spectrum of NGC 1068, which remains the brightest neutrino source observed to this day. Motivated by the aforementioned result, we initially assumed a common σ_p = 10⁵ for all AGNs, which led to an overestimation of the diffuse flux by a factor of roughly eight at E_ν = 10 TeV (see the black line in Fig. D.2). Therefore, sources with σ_p ∼ 10⁵ should constitute ∼10% of the AGN population to avoid overshooting the diffuse flux (Abbasi et al. 2025b).

Next, we assigned a random value (σ_p) to each source sampled from a power-law probability distribution, $P({\sigma }_{\mathrm{p}})\propto {\sigma }_{\mathrm{p}}^{-n}$ with n > 1 and log₁₀(σ_p)∈[3, 5]. As a result, we argue that sources with such a high proton magnetization should be rare in the population. The results shown in Sect. 3 were obtained for $P({\sigma }_{\mathrm{p}})\propto {\sigma }_{\mathrm{p}}^{-2}$. Predictions made for varying n and different probability distributions are presented in Appendix D. Finally, the neutrino energy flux from the whole sky (per steradian) using the two mock catalogs was computed as

$$\begin{array}{c}\hfill E_{\mathrm{obs}}^2{\mathrm{\Phi }}_{\nu +\overline{\nu }}(E_{\mathrm{obs}})={(4\pi )}^{-1}\sum _j\frac{4\pi }{{\mathrm{\Omega }}_{\mathrm{j}}}\sum _{i=1}^{N_{\mathrm{j}}}{E}_{\mathrm{obs}}{F}^{(i)}(E_{\mathrm{obs}}),\end{array}$$(2)

where j = {1, 2}, while Ω_j = {50  (π/180)²,  2 ⋅ 10⁴ (π/180)²}, and N_j = {1.2 ⋅ 10⁶,  7.1 ⋅ 10⁵} are the solid angles and number of sources in the 0.2 ≤ z ≤ 4 and z ≤ 0.2 sub-catalogs, respectively.

3. Results

Figure 1 shows the results for the diffuse all-flavor neutrino flux for the whole sky, marked with the solid black line, up to E_ν = 10¹⁷ eV. The solid gray line represents the maximum contribution of BL Lacs allowed by IceCube upper limits (beyond PeV energies) based on the model of Padovani et al. (2015). Moreover, with burgundy markers we show the IceCube data as given by the most recent analysis of Basu et al. (2025). Dashed lines of different colors indicate the contribution of sources with varying σ_p values. We find that the total diffuse neutrino spectrum is not particularly sensitive to the choice of P(σ_p) but to the fraction of the sources with 5 ⋅ 10⁴ ≤ σ_p ≤ 5 ⋅ 10⁵. In particular, we are able to describe the IceCube observations if the aforementioned population constitutes 10% of the total non-jetted AGNs.

Fig. 1. Diffuse all-flavor neutrino flux spectrum. Burgundy markers show the IceCube diffuse flux spectrum as obtained from the most recent analysis of Basu et al. (2025). The solid black line shows the total flux. This comprises emission from AGN coronae predicted by our model (up to ∼1015 eV) and from jetted AGNs (BL Lac sources according to Padovani et al. 2015) at higher energies. The contribution of sources with different σp values is shown with dashed colored lines – see the legend.

Our model is able to describe the observational data up to E_ν ∼ 10¹⁵ eV. Although coronae with σ_p = 10⁵ are the least common in the sample, their contribution to the total flux is the largest. This results from the dependence of the neutrino production efficiency on σ_p: As the latter increases, the neutrino production efficiency also increases until saturation is achieved, as discussed in Karavola et al. (2025). Suppression of the neutrino flux at higher energies (E_ν ≳ 10¹⁵ eV) due to muon and pion cooling leaves room for other source classes to contribute. For instance, gray solid line illustrates the maximum expected contribution from jetted AGNs – particularly BL Lac objects – based on the model of Padovani et al. (2015). The combined spectrum from both non-jetted and jetted AGNs provides a very good description of the diffuse flux across four orders of magnitude in energy.

4. Discussion and conclusions

We calculated the diffuse neutrino flux according to our reconnection-powered AGN corona model. We showed that AGN coronae could account for the majority of the diffuse neutrino background up to energies of 10¹⁵ eV, provided that approximately 10% of the population has coronae with proton magnetization in the range 5 ⋅ 10⁴ ≤ σ_p ≤ 5 ⋅ 10⁵. Moreover, we found that galaxies at 0.5 ≤ z ≤ 3.0 with 10⁴⁴ erg s⁻¹ ≤ L_X ≤ 10⁴⁶ erg s⁻¹ are the ones that mostly contribute to the diffuse neutrino flux at energies 1 TeV ≤ E_ν ≤ 1 PeV – see Appendix C. At higher energies the coronal neutrino spectrum is suppressed by pion and muon synchrotron cooling, leaving room for other source classes to contribute to multi-PeV energies.

In this study, we adopted a post-break slope of s_p = 3, which was motivated by the results of Fiorillo et al. (2024), who showed that this value is required to reproduce the neutrino spectrum of NGC 1068. Smaller values of s_p lead to a harder neutrino spectrum at energies ≳σ_pm_pc²/20. In particular, for s_p = 2, which would be expected from reconnection in the limit of a vanishing guide field, the neutrino spectral peak becomes insensitive to σ_p and is controlled by the maximum proton energy. For these reasons, zero-guide field reconnection in NGC 1068 has been disfavored by previous works (Fiorillo et al. 2024). However, the effects of pion and muon cooling were not accounted for in previous works. In Appendix B, we relax the constraint on s_p and demonstrate that once pion and muon cooling losses are included, the high-energy neutrino slope obtained for s_p = 2 closely resembles the result of s_p = 3 in the absence of heavy-particle cooling (as assumed by Fiorillo et al. 2024). As a result, reconnection in the vanishing guide-field regime appears perfectly consistent with the neutrino spectrum of NGC 1068, as long as pion and muon cooling losses are properly included.

In our coronal model, we assumed that pairs accelerated during magnetic reconnection transfer their energy to X-ray photons. Both 2D and 3D PIC simulations (Sironi et al. 2015; Werner & Uzdensky 2017) have shown that the fraction η_X of electromagnetic energy dissipated into pairs decreases with increasing guide-field strength (B_g), ranging from ∼0.5 to ∼0.1 for vanishing and strong guide fields, respectively. Because the proton energy fraction, η_p, shows a similar trend, the ratio η_p/η_X, which is crucial for our model, is rather insensitive to B_g. Instead, its value depends on the pair loading of the plasma. In particular, Petropoulou et al. (2019) found that η_p/η_X ∼ 0.1 in pair-dominated plasmas, increasing toward unity in electron-proton plasmas. Motivated by these results and given the significant pair enrichment in our model (Karavola et al. 2025), we adopted η_X = 0.5 and η_p/η_X = 0.2.

In the coronal model used in this study, reconnecting current sheets form not within the accretion disk itself but in the magnetospheric region around the black hole. Plasma from the disk can be intermittently funneled into this region if accretion is regulated by the magnetic flux that threads the inner disk (Ripperda et al. 2022). Only a small fraction, ζ, of the density of the particles in the disk is expected to enter the magnetosphere, with ζ not yet quantified in simulations. The proton magnetization scales as σ_p ∝ B²/(ζn_p, disk), and it is independent of both L_X and M since B² ∝ L_X/M² and n_p, disk ∝ L_X/M² (see Eqs. E1 and E2 in Fiorillo et al. 2024). Consequently, the distribution in σ_p needed to explain why the diffuse flux implies a corresponding dispersion in the fraction ζ across AGN coronae. Although the underlying physical driver of this dispersion is still uncertain, it may indicate that the formation and dynamics of magnetospheric current sheets differ among AGNs with varying coronal conditions.

Most studies adopt the hypothesis that the neutrino luminosity of the corona scales linearly with its X-ray luminosity (e.g., Padovani et al. 2024). A key feature of our model is that the neutrino–X-ray luminosity scaling depends on the coronal photon compactness, L_X/R. For high compactness, the corona is optically thick to pγ interactions and L_ν ∝ L_X, while L_ν ∝ L_X² otherwise (Karavola et al. 2025). This feature allowed us to loosen the commonly adopted assumption that all non-jetted AGNs follow a behavior similar to the neutrino-bright NGC 1068. As opposed to models that associate the corona with magnetized turbulence in the accretion flow (Murase et al. 2020; Fiorillo et al. 2024; Lemoine & Rieger 2025), in our scenario neutrinos are produced due to pγ interactions in the magnetospheric current sheets. In these environments, p-p collisions are not important, compared to pγ processes, as was discussed in Fiorillo et al. (2024). This is due to the fact that the plasma in the upstream region is pair-dominated, and thus the number density of non-relativistic protons is significantly lower than the density of leptons and photons.

Given that our model can reproduce nearly the entire diffuse neutrino flux at energies 1–10 TeV, it is natural to ask what implications it has for the diffuse γ-ray background. Photons with energies E_γ ≳ 10 MeV are produced due to Bethe-Heitler and pγ processes. However, these are attenuated due to the high opacity of the corona to γγ interactions with its X-ray photons. The photon energy flux at E_γ ∼ 1 − 10 MeV in the corona is comparable to the peak neutrino energy flux (Karavola et al. 2025). Therefore, the diffuse 1–10 MeV flux is expected to be comparable to the peak diffuse neutrino flux, which is a few times lower than the measured γ-ray flux in these energies (see figure 10 in Ackermann et al. 2015).

Beyond the diffuse neutrino flux, our model also has important implications for individual AGN sources. In previous work we provided source-specific predictions (Karavola et al. 2025), and a natural next step is to use the neutrino spectral templates developed here to perform a more robust statistical comparison with IceCube data. Such an analysis would allow us to infer plausible values of σ_p for AGNs with neutrino associations or to place meaningful upper limits on σ_p for non-detected sources. This is particularly relevant given that the proton magnetization of the magnetospheric environment of AGN remains one of the key unknown parameters of the model.

Acknowledgments

We thank the anonymous referee for their insightful comments that helped us clarify several points in the manuscript. D.K. thanks V.C. Karavolas for the useful conversations during the development of this paper. DFGF is supported by the Alexander von Humboldt Foundation (Germany) and, when this work was started, was supported by the Villum Fonden (Denmark) under Project No. 29388 and the European Union’s Horizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agreement No. 847523 ‘INTERACTIONS’. A.G. acknowledges funding from the Hellenic Foundation for Research and Innovation (HFRI) project “4MOVE-U” grant agreement 2688, which is part of the programme “2nd Call for HFRI Research Projects to support Faculty Members and Researchers”. L.C. acknowledges support from NSF grant PHY-2308944, NASA ATP award 80NSSC22K0667, and NASA ATP award 80NSSC24K1230. L.S. acknowledges support from DoE Early Career Award DE-SC0023015 and NASA ATP 80NSSC24K1238. This work was also supported by a grant from the Simons Foundation (MP-SCMPS-00001470) to L.S., and facilitated by Multimessenger Plasma Physics Center (MPPC, NSF PHY-2206609 to L.S.).

References

Appendix A: Synchrotron cooling of muons, pions and kaons

Neutrinos are produced by the decay of muons (μ^±), pions (π^±), and kaons (K^±, K⁰). In regions with strong magnetic fields, the synchrotron energy loss timescale can be shorter than the particle decay timescale, causing suppression (damping) of neutrino production above a characteristic energy. For a particle species i with rest mass m_i and decay timescale τ_i (defined in the particle rest frame) this critical energy is

$$\begin{array}{c}\hfill E_{\mathrm{br}}^{(i)}=\sqrt{\frac{3R}{4{\ell }_{B}c{\tau }_i}}{\left(\frac{m_i}{m_e}\right)}^{3/2}{m}_i{c}^2,\end{array}$$(A.1)

where ℓ_B = σ_TRB²/(8πm_ec²) is the magnetic compactness of the corona, a dimensionless measure of its magnetic energy density. Figure A.1 shows the corresponding break energies for muons, pions, and charged kaons in coronae with X-ray luminosities between 10⁴³ erg s⁻¹ and 10⁴⁷ erg s⁻¹, formed around black holes with masses between 10⁷ M_⊙ and 10⁹ M_⊙ and having a current sheet dimension L = 6 r_g. In our coronal model, the magnetic field strength scales as $B\propto L_{\mathrm{X}}^{1/2}/R\propto L_{\mathrm{X}}^{1/2}/M$, depending on both luminosity and black hole mass. For the parameter space considered here, synchrotron cooling of muons and pions steepens the neutrino spectrum above the corresponding break energies, which lie between 10 TeV and 10 PeV (see also Fig. A.2).

Fig. A.1. Characteristic energy where the synchrotron energy loss timescale becomes equal to the decay timescale for muons, pions, and charged kaons. The star marks the position of NGC 1068.

Fig. A.2. Neutrino spectra displaying the effect of pion and muon cooling. The blue dotted line shows the spectrum without synchrotron cooling while the red dashed and purple solid lines demonstrate the neutrino spectra when only muon cooling or both muon and pion cooling are included, respectively. The thick purple line shows the spectrum after applying a Savitzky-Golay filter to smooth out the numerical artifact appearing at ∼1014 eV. The solid and dashed vertical lines show the position of the cooling break for each species according to Eq. A.1.

To generate neutrino spectral templates, we used the numerical leptohadronic code ATHEνA (Dimitrakoudis et al. 2012). The code solves, over time, a coupled system of integro-differential equations for five particle populations–protons, neutrons, electrons/positrons, photons, and neutrinos–within a fixed volume, through a network of physical processes (see Cerruti et al. 2026, for a detailed summary). Synchrotron cooling of unstable particles (μ^±, π^±, K^±) is incorporated through an approximate treatment that avoids solving additional equations, as described in Petropoulou et al. (2014).

In this scheme, the synchrotron energy lost by each particle prior to decay is first estimated. The remaining energy is then transferred to the secondary particles produced in the decay, with their yields (production rates) precomputed using the Monte Carlo code SOPHIA (Mücke et al. 2000) ⁶. The procedure accounts sequentially for synchrotron cooling of K^± (whose decays may produce charged pions or muons), then π^± (which decay into muons), and finally muons themselves. Neutral kaons (K_S⁰, K_L⁰) are instead assumed to decay instantaneously, as are neutral pions (π⁰), with their decay products injected directly into kinetic equations of the relevant species. Notably, long-lived neutral kaons (K_L⁰) can directly produce neutrinos through semileptonic decay channels (K_L⁰ → π^± + e^∓ + ν_e, K_L⁰ → π^± + μ^∓ + ν_μ)⁷. These neutrinos (ν_e, ν_μ) are unaffected by synchrotron cooling⁸. As a result, above the characteristic pion cooling break energy, the neutrino spectrum is typically dominated by contributions from K_L⁰ decays (Kachelrieß et al. 2008).

To illustrate more clearly the effects of synchrotron cooling on the neutrino emission we show in Fig. A.2 spectra calculated for parameters similar to those of NGC 1068: L_X ≃ 3.6 × 10⁴³ erg s⁻¹, M ≃ 6.7 × 10⁶ M_⊙, B ∼ 9 × 10⁴ G and σ_p = 10⁵ (Karavola et al. 2025). The blue dotted line shows the spectrum without pion and muon cooling while the red dashed and purple solid lines demonstrate the neutrino spectra when only muon or both muon and pion cooling are included, respectively. Inclusion of charged kaon cooling does not alter the neutrino spectrum further. The solid and dashed vertical lines indicate the locations of the cooling breaks for muons and pions, respectively. The sharp feature that appears near the pion cooling energy is a numerical artifact; it is smoothed out when calculating the diffuse flux. Neutrino emission above the peak energy of the spectrum is suppressed by synchrotron cooling of pions and muons, washing out the dependence on the post-break slope of the proton spectrum (for a discussion on the fields required to affect the peak of the NGC 1068 neutrino signal see Blanco et al. 2025). The high-energy bump that is visible in the purple curve is related to the K_L⁰ decay. It is also worth noting that the suppression of neutrino flux due to synchrotron cooling of pions and muons weakens as the compactness of the corona L_X/R decreases, as illustrated in Fig. A.3.

Fig. A.3. Library of all-flavor neutrino spectral templates computed using the code ATHEνA for a discrete grid of 0.1-100 keV X-ray luminosities (in erg s−1) and coronal radii (in cm). The color of the line represents the value of σp with red, magenta and blue being 103, 104 and 105 respectively.

Figure A.3 presents the neutrino spectral templates calculated for a discrete grid of X-ray luminosity L_X and coronal radius R, accounting for the cooling of pions and muons. In each panel, we present the all-flavor neutrino spectra (in code units) computed for three different values of the proton magnetization. The suppression of neutrino emission due to synchrotron cooling of pions and muons becomes more pronounced for sources with higher X-ray luminosities and smaller radii. Notably, the neutrino spectra exhibit an approximate self-similarity with respect to λ_X, Edd, as can be seen by comparing panels along the same diagonal. This self-similarity, however, is not exact when synchrotron cooling of pions and muons becomes significant, since the corresponding characteristic break energies (defined in Eq. A.1) depend on additional parameters beyond λ_X, Edd.

Appendix B: Slopes for protons and X-rays

In our analysis, we adopted s_p = 3 for the post-break proton distribution slope, motivated by prior modeling of NGC 1068 (Fiorillo et al. 2024). This slope is controlled by the guide-field strength, B_g, with s_p = 2 in the limit of vanishing guide fields and s_p = 3 for strong fields (i.e., B_g ∼ B₀ where B₀ is the reconnecting magnetic field component Werner & Uzdensky (2017, 2024), Comisso & Jiang (2023)). In figure B.1 we show the imprint of s_p on the produced neutrino spectrum. Solid lines take into account the cooling of heavier particles, as described in Appendix A, whereas dashed lines ignore the effect. For s_p = 2, in the absence of cooling, the neutrino spectrum is drastically different than for s_p = 3, and its peak is dictated by the maximum energy of the protons rather than σ_p (compare red and blue dashed lines; see also App. C in Karavola et al. (2025)). However, when pion and muon cooling are included, the high-energy (≳30 TeV) neutrino spectrum softens, and the peak neutrino energy as well as the high-energy slope of the spectrum are essentially determined by the magnetic field.

Fig. B.1. Neutrino spectra displaying the effect of the proton distribution post-break slope sp for the NGC 1068 case– see Appendix A. Red and blue refer to sp = 2 and 3 respectively. Solid lines include π and μ cooling while dashed lines neglect it. The solid and dashed vertical lines show the position of the cooling break for each species according to Eq. A.1. The dotted line refers to the neutrino energy produced by the most energetic protons of the distribution.

To demonstrate the effect s_p has on the diffuse neutrino flux, we performed a simple calculation in which we assigned the NGC 1068 spectral template (shown in fig. B.1) to all galaxies after re-scaling it according to L_X and the pγ efficiency (see Karavola et al. (2025) for analytical estimates). We perform this calculation for both s_p = 2 and 3 – see magenta and blue lines in figure B.2 respectively. At PeV energies, we clearly see that the single-template analysis with s_p = 3 does not coincide with our calculation, because sources with lower magnetic fields than NGC 1068 dominate in that regime – see also Appendix D. As a result, if we were to adopt s_p = 2 for all sources, we would expect the diffuse flux to surpass the observations. Given that the single-template analysis is not representative, to fully account for the effects of the post-break slope, one would have to repeat our analysis for varying s_p, and we defer this to a future study.

Fig. B.2. Same as Fig. 1. Magenta and blue lines correspond to the calculation of the diffuse flux assuming that all the sources have the same spectral shape as NGC 1068 (for sp = 2 and 3 respectively – see fig. B.1) re-normalized in respect to LX and the pγ efficiency. Dashed lines refer to the coronal model itself while solid lines also include the BL Lac population.

Furthermore, Trakhtenbrot et al. (2017) showed that the X-ray photon index of non-jetted AGNs is not strongly dependent on the bolometric X-ray luminosity in the 2-10 keV band with its mean value being ∼2. Motivated by the latter result, for the scope of this work, we have assumed s_x = 2. In order to demonstrate the effect of the photon index we show, in figure B.3, the distributions of X-rays (top panel), neutrinos (middle panel) and protons (bottom panel) for a single source with s_x= 1, 2, and 3. For all three cases, we have used representative parameters close to NGC 1068– see Appendix A– for the black hole mass M and proton magnetization σ_p. However, to calculate the bolometric X-ray luminosity L_X|_sx, we considered fixed the luminosity in the 2-10 keV energy band, since it is the value provided by the mock catalog we have used (see Appendix C). We see that for s_x= 1 and 3 (red and blue lines), the neutrino spectra are more luminous by a factor of ∼5 compared to the case of s_x = 2 (black line), mostly due to the fact that L_p ∝ L_X is higher in these cases in comparison to the s_x = 2 one. As a result, if all of the sources were described by s_x = 1 or 3, then we would expect the diffuse flux to be higher by a factor of 5, as well. Moreover, it is interesting to note that high-energy protons (E ≳ σ_pm_pc²) appear to “cool down” more as s_x increases. The latter is due to the fact that for higher photon indices there are more photons toward the lower end of the distribution and thus, the efficiency of closer-to-threshold photohadronic interactions is increased.

Fig. B.3. Luminosity distributions for X-rays (top), neutrinos (middle) and protons (bottom) for 3 values of the photon index sx–see legend.

Appendix C: AGN mock catalogs

The galaxy and AGN mock catalogs used in this work are constructed using empirical relations on the incidence of accretion events in galaxies (e.g., Georgakakis et al. 2017; Aird et al. 2018). The key quantity in this calculation is the specific accretion rate defined as λ_SAR ∝ L_AGN/M_★, where L_AGN is a proxy of AGN accretion luminosity and M_★ is the stellar mass of the host galaxy. The λ_SAR measures the accretion luminosity produced per unit stellar mass of the host galaxy and under certain assumptions, namely bolometric correction and scaling relation between stellar mass and black hole mass, it can be converted to the Eddington ratio, λ_Edd. For the purpose of our analysis, a convenient feature of λ_SAR is that it can be measured observationally for large and unbiased samples of AGN. This is unlike the Eddington ratio that requires a direct estimate of the black hole mass. This latter step is observationally challenging as it relies on independent methods, such as measurement of broad optical emission lines or accretion disk SED fitting (e.g., Peterson 2014; Varisco et al. 2018; Popović 2020; Kyriopoulos et al. 2026). For certain populations, such as obscured or low-luminosity AGN with no detectable broad emission lines, black hole mass measurements turn out to be impossible.

Extragalactic survey fields that benefit from observations in many different wavebands, from X-ray to the far-infrared, have made possible the compilation of large and nearly unbiased AGN samples, for example via detection at X-rays, out to high redshift and the estimation of their host galaxies properties, such as stellar masses. Such observations allow for the measurement of the fraction of galaxies in a given stellar mass interval that host AGN with a specific accretion rate λ_SAR. These fractions can then be turned into probability distribution functions, P(λ_SAR), which describe the probability of accretion events with λ_SAR in galaxies (e.g., Georgakakis et al. 2017; Aird et al. 2018). The functional form of P(λ_SAR) can be approximated by a broken power-law with parameters that evolve with cosmic time. By construction, the P(λ_SAR) has the feature that if convolved with the stellar mass function of galaxies, it yields the AGN luminosity function

$$\begin{array}{c}\hfill \mathrm{\Phi }(L_{\mathrm{AGN}},z)=\int \psi (M_{★},z)·P[{\lambda }_{\mathit{SAR}}(L_{\mathrm{AGN}},M_{★}),z]·\mathrm{d}log{M}_{★},\end{array}$$(C.1)

where Φ(L_AGN, z) and ψ(M_★, z) are the AGN luminosity function and the galaxy stellar mass function, respectively, at redshift z and the specific accretion rate is expressed as a function of L_AGN and M_★. It is exactly this feature that we exploit in this work to populate mock galaxies randomly drawn from a stellar mass function with AGN.

In practice, parametric forms for Φ(L_AGN, z) and ψ(M_★, z) are assumed to solve the convolution equation above and estimate P(λ_SAR, z). We can then draw galaxy samples from the stellar mass function and paint AGN on them by probabilistically assigning them specific accretion rates from the empirically determined P(λ_SAR, z). The end product of this approach is a catalog of mock galaxies, each of which is associated with a redshift, stellar mass, λ_SAR and hence an AGN luminosity L_AGN ∝ λ_SAR ⋅ M_★. The last step is to assign each galaxy a black hole mass (and therefore an Eddington ratio given the AGN luminosity) using observationally inferred scaling relations between black hole and stellar masses.

We implement the above methods as described in Georgakakis et al. (2020). We adopt the X-ray AGN luminosity function of Ueda et al. (2014) and the parametric stellar mass functions of Ilbert et al. (2013). The specific accretion rate probability distribution function follows a three-segment broken power-law, motivated by the results of non-parametric studies (e.g., Georgakakis et al. 2017; Aird et al. 2018). For the black-hole versus stellar mass correlation we adopt the one of Marconi & Hunt (2003), M_BH = 2 ⋅ 10⁻³M_*. Experimentation with other scaling relations in the literature shows that our results are not sensitive (factor of < 2) to the adopted relation between M_BH and M_★. The final mock catalog includes for each galaxy the: stellar mass, redshift, specific accretion rate, black hole mass, and 2-10 keV X-ray luminosity of the AGN.

Since the generation of the mock catalog includes a sampling step for the stellar mass function, one has to assume a volume that the galaxies occupy, which translates to a fiducial sky area of the mock survey. The baseline mock catalog, which extends up to z = 4, corresponds to an area of 50 deg². However, for this solid angle the sampled volume at low redshift, z ≲ 0.2, is too small and our results are affected by poor statistics. We therefore choose to complement the baseline mock with a low redshift one with a solid angle of 20, 000 deg² and limited to z ≤ 0.2. For our analysis, we use the low redshift mock for z ≤ 0.2 and the baseline one for higher redshifts.

The modeling approach above has the advantage of producing realistic mock AGN catalogs that include host galaxy properties (e.g., stellar mass, black hole mass) and by construction reproduce fundamental statistics of the AGN population, such as space densities, host galaxy stellar mass function, large scale structure, multiwavelength characteristics (see Georgakakis et al. 2020, for details).

Figure C.1 shows the distribution of different properties of the mock galaxies: (a) bolometric (0.1 − 100 keV) X-ray luminosity, (b) X-ray Eddington ratio (defined as L_X/L_Edd), (c) black hole mass, and (d) bolometric X-ray flux. For each panel, different colors correspond to different redshift bands (see legend). The mock catalog contains AGN with black hole masses from ∼10⁷ M_⊙ up to ∼10⁹ M_⊙ with a distribution that does not depend on their distance, as seen in panel c. AGN with high bolometric X-ray luminosities (L_X ≳ 10⁴⁵ erg s⁻¹) are mostly found at large distances (z ≳ 1), while AGN with higher X-ray fluxes can be found relatively closer to us (z ≲ 1).

Fig. C.1. Histograms displaying properties of AGN from the mock catalog: (a) bolometric X-ray luminosity LX > 1042 erg s−1, (b) X-ray Eddington ratio λX, Edd, (c) black hole mass, and (d) the bolometric X-ray flux. Color shows different redshift bands, see legend.

Appendix D: Decomposition of the diffuse neutrino flux

In our model, the neutrino luminosity of a corona that is optically thin to pγ interactions scales linearly with the proton plasma magnetization (Karavola et al. 2025). This has an imprint on the diffuse flux as demonstrated in Fig. 1. Another key parameter of our model, which controls the expected neutrino luminosity, is the bolometric X-ray luminosity of the corona (see Eq. 3.6 in Karavola et al. (2025)). Moreover, the redshift of the source affects the observed neutrino flux. In Fig. D.1 we show the contribution of different X-ray luminosity bins (panel a), redshift bins (panel b) and magnetic field bins (panel c) to the total neutrino flux. Solid lines are the same as in Fig. 1, with colors representing the range of X-ray luminosities, redshifts or magnetic fields. One can see that the dominant contribution originates from galaxies at 0.5 ≤ z ≤ 3.0 with 10⁴⁴ erg s⁻¹ ≤ L_X ≤ 10⁴⁶ erg s⁻¹ and 10³ G ≤ B ≤ 10⁵ G.

Fig. D.1. Same as Fig. 1 except that lines are color coded based on the bolometric X-ray luminosity (panel a), the redshift z (panel b) and on the magnetic field (panel c) of AGN as indicated by the inset legend. Panel d shows the exact same calculation as in Fig. 1 but for σp extending up to 106.

Furthermore, panel d in Fig. D.1 shows the same calculation as in Fig. 1 but with σ_p extending up to 10⁶. We note that this solution can better describe the observations in the PeV energy range. However, a softer power-law distribution with n = 4.5 is needed to avoid overshooting the data. Nonetheless, the take-home message remains the same: only a small fraction of the non-jetted AGN population should have high proton plasma magnetizations. In particular, ∼15% of the sources should have σ_p ≥ 5 × 10⁴ in order to match the diffuse flux observations. In this case, coronae with magnetizations higher than the one describing NGC 1068 constitute a very small fraction of approximately 4% of the population. As a result, the detection probability of such high-σ_p coronae is expected to be low.

Finally, Fig. D.2 shows the impact of different probability functions P(σ_p) on the diffuse flux calculation. Black solid line assumes σ_p = 10⁵ for all sources of the catalog, while colored lines show results for $P({\sigma }_{\mathrm{p}})\propto {\sigma }_{\mathrm{p}}^{-n}$ with the color indicating the value of the slope n (see color bar). Additionally, the black dashed line is computed using a Gaussian distribution with respect to log₁₀(σ_p) with mean value μ = 3 and standard deviation σ = 1.5 so that, again, 10% of the non-jetted AGNs have a value of σ_p ∼ 10⁵. If we attribute σ_p = 10⁵ to all galaxies, the diffuse flux overshoots the data by almost one order of magnitude at E_ν ∼ 10 TeV. Also, we note that changes in the slope n of the σ_p probability distribution weakly affect the diffuse flux calculation, with 1 < n < 3 resulting in the best description of the IceCube measurements Basu et al. (2025). Interestingly, adopting a Gaussian probability distribution leads to very similar neutrino spectra as those obtained for the power-law distributions. We conclude that the shape of P(σ_p) does not have a large impact on the results as long as the fraction of galaxies with 5 × 10⁴ ≲ σ_p is about 10%.

Fig. D.2. Similar as Fig. 1 except that black solid line shows result of the assumption that all AGN have σp = 105, while colored lines occur from power-law probability distributions of different slopes n (see colorbar) and black dashed line comes from P(σp) being a gaussian distribution with a mean value μ = 3 and a standard deviation of σ = 1.5.

All Figures

Fig. 1. Diffuse all-flavor neutrino flux spectrum. Burgundy markers show the IceCube diffuse flux spectrum as obtained from the most recent analysis of Basu et al. (2025). The solid black line shows the total flux. This comprises emission from AGN coronae predicted by our model (up to ∼1015 eV) and from jetted AGNs (BL Lac sources according to Padovani et al. 2015) at higher energies. The contribution of sources with different σp values is shown with dashed colored lines – see the legend.

In the text

	Fig. A.1. Characteristic energy where the synchrotron energy loss timescale becomes equal to the decay timescale for muons, pions, and charged kaons. The star marks the position of NGC 1068.
In the text

Fig. A.2. Neutrino spectra displaying the effect of pion and muon cooling. The blue dotted line shows the spectrum without synchrotron cooling while the red dashed and purple solid lines demonstrate the neutrino spectra when only muon cooling or both muon and pion cooling are included, respectively. The thick purple line shows the spectrum after applying a Savitzky-Golay filter to smooth out the numerical artifact appearing at ∼1014 eV. The solid and dashed vertical lines show the position of the cooling break for each species according to Eq. A.1.

In the text

	Fig. A.3. Library of all-flavor neutrino spectral templates computed using the code ATHEνA for a discrete grid of 0.1-100 keV X-ray luminosities (in erg s−1) and coronal radii (in cm). The color of the line represents the value of σp with red, magenta and blue being 103, 104 and 105 respectively.
In the text

Fig. B.1. Neutrino spectra displaying the effect of the proton distribution post-break slope sp for the NGC 1068 case– see Appendix A. Red and blue refer to sp = 2 and 3 respectively. Solid lines include π and μ cooling while dashed lines neglect it. The solid and dashed vertical lines show the position of the cooling break for each species according to Eq. A.1. The dotted line refers to the neutrino energy produced by the most energetic protons of the distribution.

In the text

	Fig. B.2. Same as Fig. 1. Magenta and blue lines correspond to the calculation of the diffuse flux assuming that all the sources have the same spectral shape as NGC 1068 (for sp = 2 and 3 respectively – see fig. B.1) re-normalized in respect to LX and the pγ efficiency. Dashed lines refer to the coronal model itself while solid lines also include the BL Lac population.
In the text

	Fig. B.3. Luminosity distributions for X-rays (top), neutrinos (middle) and protons (bottom) for 3 values of the photon index sx–see legend.
In the text

	Fig. C.1. Histograms displaying properties of AGN from the mock catalog: (a) bolometric X-ray luminosity LX > 1042 erg s−1, (b) X-ray Eddington ratio λX, Edd, (c) black hole mass, and (d) the bolometric X-ray flux. Color shows different redshift bands, see legend.
In the text

	Fig. D.1. Same as Fig. 1 except that lines are color coded based on the bolometric X-ray luminosity (panel a), the redshift z (panel b) and on the magnetic field (panel c) of AGN as indicated by the inset legend. Panel d shows the exact same calculation as in Fig. 1 but for σp extending up to 106.
In the text

	Fig. D.2. Similar as Fig. 1 except that black solid line shows result of the assumption that all AGN have σp = 105, while colored lines occur from power-law probability distributions of different slopes n (see colorbar) and black dashed line comes from P(σp) being a gaussian distribution with a mean value μ = 3 and a standard deviation of σ = 1.5.
In the text