Berezinsky Hidden Sources: An Emergent Tension in the High-Energy Neutrino Sky?

Abstract

The IceCube Collaboration has recently reported compelling evidence of high-energy neutrino emission from NGC 1068, and also mild excesses for NGC 4151 and CGCG420-015, local Seyfert galaxies. This has increased the interest along neutrino emission from hot-corona surrounding the super massive black holes of Seyfert Galaxies. In this paper, we revisit phenomenological constraints on the neutrino emission from hot-coronae of seyfert galaxies, using an assumption of equi-ripartition between cosmic-rays and magnetic energy densities. We show that not only these sources are consistent with such an assumption but also that the data point towards low beta plasma parameters inside Seyfert Galaxies. We exploit this finding to constrain the Seyfert diffuse neutrino flux and we obtain that, in order not to overproduce neutrinos, not all the sources can be in an equi-ripartition state. We conclude (along with previous findings) that seyfert galaxies cannot explain the diffuse neutrino spectrum above $\sim 100\,\rm TeV$, allowing space for other astrophysical sources.

1 Introduction

Seyfert Galaxies are Active Galactic Nuclei (AGNs) powered by super-massive black holes (SMBHs) positioned in their core, whose accretion disk related processes produce very energetic phenomena, such as the creation of a hot corona region surrounding the black hole [1]. In the late 1970s V.S. Berezinsky proposed that the hot corona region might produce and accelerate high-energy cosmic rays (CRs) [2] (see also [3] and references therein for further details). these CRs were predicted to predominantly interact with ambient particles leading to the production of $\gamma$-rays and neutrinos, with only the latter being able to escape these very-high dense environments [2, 4, 5]. Despite this idea could have not been experimental explored for many years, current $\gamma$-ray and neutrino telescopes have strongly improved our knowledge about astrophysical messengers, making it possible to quantitative constrain this scenario. Firstly, many analyses have highlighted a tension between the diffuse $\gamma$-ray [6] and neutrino [7, 8, 9, 10] fluxes, especially because of the large neutrino flux below 100 TeV, exploring the potential role of hidden CR accelerators [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] (see also [22] for other details). Furthermore, the IceCube collaboration has recently found a 4.2 $\sigma$ excess above the background-only hypothesis of 72 high-energy neutrino events coming from the direction of NGC 1068, a nearby AGN [23], inferring a power-law differential neutrino flux $(\sim E^{-3.2})$ approximately ten times higher than expected from a $\gamma$-ray transparent source, taking into account the $\gamma$-ray emission observed by the Fermi-LAT telescope from this source [24, 25, 26]. This, from one hand, has confirmed the idea of the presence of hidden high-energy neutrino sources in the Universe, but from the other hand, it has triggered a lot of controversy for the interpretation of this observation. Indeed, a lot of models have been proposed to explain such observations [3, 27, 28, 21, 19, 29, 30, 31, 32, 33, 34, 35, 36]. At the moment, the most natural explanation resides in high-energy neutrino emission in the hot corona emission where the high x-ray photon density is able to strongly suppress the high-energy $\gamma$-ray flux accompanying the neutrino flux, leading to a natural explanation for a $\gamma$-opaque source [37]. However, there is not consensus about the acceleration mechanisms for high-energy CRs inside the coronae. For instance, Refs. [21, 30] have proposed stochastic processes as acceleration mechanism while [20, 28] have proposed diffusive shock acceleration, both normalising the spectrum imposing that the CR pressure is just a small fraction of the thermal gas pressure. On the other hand, Ref. [3] recently explored the role of magnetic re-connection, normalising the spectrum considering equi-partition between CR energy density and magnetic energy density. On the contrary, Ref. [27] ascribed a fraction of the accretion mass rate energy to the production of high-energy CRs. Nonetheless, all the models seemingly share some common properties: firstly, results suggest energy losses be very important in the hot corona regions. In fact, Beithe- Heitler pair production leads to electromagnetic cascades shifting the energy of photons in the MeV energy range [21, 19, 38, 39]. Secondly, the predicted neutrino luminosity of hot corona is expected to be proportional to the AGN X-ray luminosity. Therefore, it is fundamental to observe of other seyfert galaxies to shed light onto the processes occurring in these kind of environments [40, 41]. On this regard, IceCube has searched for other seyfert galaxies in the local Universe, both as individual emitters and performing a stacking search [42, 43, 44, 45]. The authors found a $\sim 3\sigma$ post-trial excess in the vicinity of NGC 4151 [43] with spectral index $\gamma=2.83_{-0.28}^{+0.35}$. Furthermore, an excess has been found also in the vicinity of CGCC 420-015 $(2-2.5~{}\sigma)$ [42], with a spectrum very similar to NGC 1068 spectrum. By constrast, Ref. [41], analysing 10 years of public IceCube data, found an excess even in the direction of NGC 3079, another seyfert galaxy. Seemingly, these potential observations point towards a correlation between the neutrino and hard x-ray luminosity for these sources, although the results are too premature to draw robust conclusions [33]. In fact, while Refs. [38, 39] have produced optimistic estimates for the neutrino emission of several local seyfert galaxies, Ref. [46] has, on the contrary, argued that the seyfert hot coronae should not be able to produce a neutrino flux as high as the one observed by the IceCube collaboration. Furthermore, it is not clear if the steep spectra measured are byproducts of copious p$\gamma$ interactions happening in the hot coronae or if they directly represent the cut-off of the injected CR fluxes [39]. In this paper, we derive new data-driven constraints on the seyfert neutrino emission exploiting all the state-of-the-art observations of the IceCube collaboration also exploring the potential role played by these sources to diffuse neutrino flux. To such a purpose, we employ a model which takes into account both CR escape mechanisms as well as CR energy losses and consistently solve the CR transport equation inside the hot-coronae. We normalise the spectrum assuming equi-ripartition between CR and magnetic energy densities, which represent an absolute maximum value for the CR normalization. We also ensure that the injected CR luminosity is only a fraction of the bolometric AGN luminosity which is self-consistently calculated through the background photon energy density. We show that the four sources which provide an excess into IceCube data are consistent with an equi-ripartion between CRs and magnetic density. In fact, the fit prefers small beta plasma values pointing towards high magnetic field values and a rather high energetic carried by CRs into Seyfert Galaxies (see below for further details). We then extrapolate this information to the whole seyfert population (diffuse neutrino flux), using the X-ray luminosity function (see [47]), showing that it slightly overproduce the diffuse neutrino flux inferred by the ICeCube collaboration through the starting track sample [10] in the $\sim 1-10\,\rm TeV$ range. Future observations from the upcoming neutrino telescopes, such as KM3NeT [48], IceCube gen 2 [49], P-ONE [50] and the TRIDENT [51] telescopes will be fundamental to unveil the role of seyfert galaxies to the diffuse neutrino spectrum as well as the role of other astrophysical accelerators. The paper is organised as follows: in Sec. 2, we describe the CR transport model and the neutrino emission of seyfert galaxies, in Sec. 3, we summarise the X-ray and neutrino observations for the four sources observed by IceCube. In sec. 4, we quantitative test our model using the observations and report the results. In Sec. 5, we discuss the role of seyfert galaxies into the diffuse neutrino spectrum and in Sec. 6 we draw our conclusions. Finally, in appendix A, we provide details on the dynamical timescales in the hot-coronae and the neutrino production efficiency for proton-proton collisions and photomeson interactions.

2 On the Neutrino Emission of Seyfert Galaxies

The hot coronae regions are generally very small and extend only for 1-100 Schwarzschild radii $(\mathcal{R}_{S})$ around the SMBHs [38]. In this work, we use spherical geometry with a radius $R=r\cdot\mathcal{R}_{S}$, fixing $r=20$ (consistent with the recent upper limits obtained by Ref. [39] with sub-GeV gamma-rays data). The CR transport equation reads

$$\frac{N_{CR}(E)}{\tau_{\rm esc}}-\frac{d}{dE}\bigg{[}\frac{E}{\tau_{\rm loss}}N_{CR}(E)\bigg{]}=Q(E)$$ (2.1)

where $\tau_{\rm esc}$ is the escape timescale, $\tau_{\rm losses}$ is the energy loss timescale (see app. A for details about the modelling on timescales) and finally $Q(E)$ is the injection rate of CRs. Eq. 2.1 is solved through the green function as (see for instance [20])

$$N_{CR}(E)=\frac{\tau_{\rm loss}}{E}\int_{E}^{+\infty}Q(E_{1})e^{-G(E,E_{1})}dE_{1}$$ (2.2)

with

$$G(E,E_{1})=\int_{E}^{E_{1}}dE_{2}\frac{\tau_{\rm loss}(E_{2})}{E_{2}\cdot\tau_{\rm esc}}$$ (2.3)

The injection CR rate is assumed to be a power-law with an exponential cut-off at $E_{\rm max}=200\,\rm TeV$, namely $Q(E)=A(\frac{E}{m_{p}c^{2}})^{-\gamma}e^{-E/E_{\rm max}}$. The normalisation parameter $A$ is fixed so that the CR energy density $U_{CR}$ is in equi-ripartion with the magnetic energy density $U_{B}$, namely $U_{CR}=U_{B}$, with

$$U_{CR}=\int_{m_{p}c^{2}}^{+\infty}EN_{CR}(E)dE$$ (2.4)

and $U_{B}=B^{2}/8\pi$. The magnetic field is set as [21]

$$B=\sqrt{\frac{8\pi n_{p}K_{B}T_{p}}{\beta}}$$ (2.5)

where $K_{B}$ is the Boltzmann constant, $T_{p}$ is the thermal proton virial temperature and finally $\beta$ is the plasma parameter which represents the ratio between the thermal and magnetic pressure [21]. The proton virial temperature is $m_{p}c^{2}/(K_{B}\cdot r)\simeq 6.0\cdot 10^{10}(r/30)^{-1}\,[\rm K]$ [21]. It is important to take into account that the luminosity injected in CR is just a fraction of the bolometric luminosity of the AGN disk $(L_{\rm ph})$. Therefore, we also constrain the normalization parameter $A$ with

$$\int_{m_{p}c^{2}}^{+\infty}EQ(E)dE\leq\eta L_{\rm ph}$$ (2.6)

The bolometric luminosity gets contribution both from the accretion disk in the ultraviolet (UV) regime (usually also called the blue bump of AGNs) and from the coronae in the x-ray luminosity [20, 21, 3]. In this work, we utilise the model put forward by Ref. [20], where the number density of background photons can be modelled in terms of the x-ray luminosity of the AGN in the $[2-10]\,\rm KeV$ band (here forth denominated simply $L_{X}$) (see [20] for more details). We fix $\eta=25\%$ for the maximal efficiency conferred into CRs. In order to calculate the pp neutrino production rate, we employ the analytical prescription of ref. [52]

$$Q_{\nu+\bar{\nu}}^{pp}(E)=\frac{1}{3}c\cdot n_{p}\int_{10^{-3}}^{1}\sigma_{pp}\big{(}\frac{E}{x}\big{)}N_{CR}\big{(}\frac{E}{x}\big{)}\tilde{F}_{\nu}\big{(}x,\frac{E}{x}\big{)}\frac{dx}{x}$$ (2.7)

where $\tilde{F}_{\nu}(x,\frac{E}{x})$ is defined in Ref. [52] and it takes into account all the neutrinos produced in the interaction, while the factor 1/3 takes into account neutrino oscillation. Since for $E_{CR}\lesssim 200\,\rm TeV$, the neutrino production is dominated by pp collisions (see appendix A for details), for photomeson neutrino production, we estimate the neutrino production rate with multi-messenger relations [19, 53]

$$E_{\nu}^{2}Q_{\nu+\bar{\nu}}^{p\gamma}\simeq\frac{1}{8}\bigg{[}\frac{E_{\rm CR}^{2}N_{\rm CR}(E)}{\tau_{p\gamma}}\bigg{]}|_{E_{\rm CR}\simeq 20E_{\nu}}$$ (2.8)

we evaluate the neutrino luminosity of the sources between $1-10\,\rm TeV$ as

$$L_{\nu+\bar{\nu}}^{1-10\,\rm TeV}=\int_{1\,\rm TeV}^{+\infty}E_{\nu}Q_{\nu+\bar{\nu}}^{\rm tot}(E_{\nu})dE_{\nu}$$ (2.9)

with $Q_{\nu+\bar{\nu}}^{\rm tot}(E)=Q_{\nu+\bar{\nu}}^{pp}(E)+Q_{\nu+\bar{\nu}}^{p\gamma}(E)$. The differential neutrino flux at Earth is

$$\phi_{\nu+\bar{\nu}}(E)=\frac{(1+z)^{2}}{4\pi D_{l}(z)^{2}}Q_{\nu+\bar{\nu}}^{\rm tot}(E(1+z))$$ (2.10)

where $z$ and $D_{l}(z)$ are respectively the redshift and the luminosity density.

3 Neutrino vs X-ray Luminosity

In this section, we summarise the observations regarding the four neutrino sources which give an excess into IceCube data. For consistency and in order to avoid any bias, we use the values reported [54] in the BASS AGN catalogue for all the sources for redshifts and the x-ray luminosities. In particular, regarding the x-ray flux, we use the intrinsic x-ray (absorption-corrected) in the 2-10 KeV band. Tab. 1 summarises all the values considered.

Source	redshift (z)	$D_{l}(z)(\rm Mpc)$	$F_{X}(10^{-12}\,\rm erg\,\rm cm^{-2}\,\rm s^{-1})$
NGC 1068	$3.03\cdot 10^{-3}$	13.4	$268.30$
NGC 4151	$3.14\cdot 10^{-3}$	13.9	84.80
CGCG 420-015	$2.959\cdot 10^{-2}$	133.9	50.50
NGC 3079	$3.40\cdot 10^{-3}$	15.1	6.60

Table 1: The table summarises the value we use for redshift, luminosity distance and the intrinsic x-ray fluxes of the sources.

we compute the luminosity distance through redshift and using the standard cosmological parameters $\Omega_{M}=0.31$, $\Omega_{\Lambda}$ and $H_{0}=67.74\,\rm Km\,\rm s^{-1}\,\rm Mpc^{-1}$ [55]. We compute $L_{X}$ as [54]

$$L_{X}=F_{x}4\pi D_{l}^{2}(z)(1+z)^{2-\Gamma}\simeq F_{x}4\pi D_{l}^{2}(z)$$ (3.1)

where we can neglect the contribution coming from redshift because these sources are very near and their x-ray spectrum is near to a $E^{-2}$ dependence of the x-ray flux.

It is paramount to correctly assess the uncertainty affecting the values reported in Tab. 1. Firstly, there is uncertainty in the distance value; for instance, Refs. [32, 3] report $10.1\,\rm Mpc$ as the best quoted value for the distance of NGC 1068. We conservatively consider a 0.13 dex uncertainty in the distance in order to be consistent with this estimate. Furthermore, there is an uncertainty on the intrinsic x-ray flux, since the geometry of the absorption region in AGNs is unknown [56]. In fact, estimates on $L_{X}$ for NGC 1068 span from $\sim 5\cdot 10^{42}\,\rm erg\,\rm s^{-1}$ [54] to $\sim 5\cdot 10^{43}\,\rm erg\,\rm s^{-1}$ [56]. Therefore, we consider an uncertainty of 0.7 dex for $F_{X}$ in order to obtain an overall range for $L_{X}$ consistent with both estimates (we sum the uncertainty quadratically) leading to a total 0.75 dex uncertainty on $\log L_{X}$.

For the neutrinos, we consider the energy range $1-10\,\rm TeV$ which is interval where IceCube is most sensitive to very steep spectra observed for these sources [23]. Therefore, we compute the luminosity as

$$L_{1-10\,\rm TeV}=4\pi D_{l}(z)^{2}(1+z)^{2-\gamma}\int_{1\,\rm TeV}^{10\,\rm TeV}E\phi_{\nu}(E,\gamma)dE\simeq 4\pi D_{l}(z)^{2}\int_{1\,\rm TeV}^{10\,\rm TeV}E\phi_{\nu}(E,\gamma)dE$$ (3.2)

where we neglect the dependence over z, because all the sources reside in the very local Universe. For NGC 1068, NGC 4151 and CGCG420-015, we consider the best-fit and the 68% CL contours for the power-law fit recently published by the IceCube collaboration [43, 42] (see also [23, 44, 45]); while for For NGC 3079, we utilise the spectral energy distribution (SED) shown by [41]. For the uncertainty, we consider the statistical $68.3\%$ CL uncertainty as well as the distance uncertainty. Fig. 1 shows $L_{\nu}^{1-10\,\rm TeV}$ in terms of $L_{X}$ for all the sources. We report NGC 1068, NGC 4151 and CGCG420-015 as black datapoints, while NGC 3079 is reported as grey data point.

Figure 1: $L_{\nu}^{(1-10\,\rm TeV)}$ as a function of the $L_{X}$. The uncertainty on $L_{\nu}$ comprises distance and statistical uncertainty, while the uncertainty on $L_{X}$ is set at 0.75 dex in order to be consistent with all the $L_{X}$ values inferred in the literature. NGC 1068, NGC 4151, CGCG 420-015 are shown as black datapoints while NGC 3079 as a grey datapoint. The blue band corresponds to the $1\sigma$ uncertainty band from our statistical analysis, while the golden dashed line corresponds to the best-fit scenario.

4 Statistical Analysis and Results

In this section, we test if our theoretical model (2) can follow the experimental results outlined in the previous section 3. For this purpose, we define the following chi-square

$$\chi^{2}(\gamma,\beta)=\sum_{i}\frac{(\ logL_{\nu}^{\rm obs,i}-\rm logL_{\nu}^{\rm mod}(\rm logL_{X}^{\rm obs,i}\gamma,\beta))^{2}}{\sigma_{\ logL_{\nu}}^{2}+\bigg{(}\frac{\partial\rm logL_{\nu}^{\rm mod}(\rm log(L_{X}^{\rm obs,i}),\gamma,\beta))}{\partial\ logL_{X}}\bigg{)}^{2}\sigma_{\ log_{L_{X}}}^{2}}$$ (4.1)

The $\partial\ logL_{\nu}/\partial L_{X}$ term allows us to take into account also the uncertainty on $L_{X}$. We use $\Delta\chi^{2}=\chi^{2}(\gamma,\beta)-\rm min\,\chi^{2}$ as a mean to probe the parameter space for $\gamma\in[0,3]$ and $\beta\in[10^{-2},10]$. Consistently with the Wilks’s theorem [57], the $\Delta\chi^{2}$ distribution is a chi-squared with a degree of freedom equal to the number of free parameters in the fit (2).

Therefore, the exclusion limits at 1,2 and $3\sigma$ are respectively defined as $\Delta\chi^{2}=2.28,6,12$. Fig. 2 shows the best-fit scenario as well as the 1,2 and $3\sigma$ contours. In Fig. 1, we also report the $1\sigma$ band for $L_{\nu}$ as well as the best-fit scenario obtained. We also compare our predicted SEDs with the ones obtained by IceCube (see Figs. 3 and 4). In particular, we show the best-fit scenario (golden dashed line) and the $1\sigma$ band (orange) allowing for $L_{X}$ to vary within the uncertainty reported in Fig. 1 fixing $D_{L}$ to the values shown in Tab. 1. For NGC 1068, we also report the expected KM3NeT differential sensitivity after 10 years of full operation [48].

Figure 2: $\Delta\chi^{2}$ contours for our fit. the color scheme represents, respectively, the contours at $1,2$ and $3$ $\sigma$. The red star represents the best fit scenario.

Despite all the uncertainties, the results allow us to derive some general conclusions: the four sources are consistent with the assumption of equi-ripartition between CRs and the magnetic energy density. In fact, low values for $\beta~{}(\lesssim 1)$ ([39]) are preferred by the data in order to increase the magnetic field and so the CR flux normalisation. However, Ref. [46] has argued that in seyfert galaxies might be nonphysical to get low beta parameters since they are not usually characterised by strong jets. Along with previous findings [3, 19, 38, 39], very hard spectral indexes are needed in order to explain the measurements. Indeed, $\gamma\simeq 0.7$ as a best-fit value is very near to the value predicted by stochastic CR acceleration mechanism [30, 19] and magnetic reconnection [3]. In general, $\gamma\gtrsim 1.9$ are in tension at $\sim 2\sigma$ with current observations, because softer injected spectra would need a higher normalisation to fit the observations. On this regard, we stress that with lower $\eta$ values in Eq. 2.6, the model fails to completely fit all the data especially for NGC 1068 and CGCG420-015, which means that if the neutrino production indeed takes place in the hot corona, an enormous energetics is carried by CRs. Recently, Ref. [39] has proposed that CRs are only accelerated in specific energy ranges $(\sim 10-200\,\rm TeV)$ in order to reduce the tension with energetics. We leave this scenario for future exploration. Finally, we stress that for CGCG420-015 the model can barely be consistent with the observations.

Figure 3: Left: $1\sigma$ SED band measured by IceCube [23, 43] (blue band) for NGC 1068 compared with predictions of our statistical analysis. In particular, we show the best-fit scenario (golden dashed line) and the $1\sigma$ band allowing $L_{X}$ to vary within the uncertainty shown in Fig. 1 and fixing $D_{L}$ to the values shown in Tab. 1 (orange band). We also report the expected KM3NeT differential sensitivity after 10 years of full operation [48]. Right: the same as on the left but for NGC 4151. The SED is directly taken from [43] (see also [44, 45]).

Figure 4: Left: $1\sigma$ SED band for CGCG420-015 (blue band) [42, 44, 45])) compared with predictions of our statistical analysis. In particular, we show the best-fit scenario (golden dashed line) and the $1\sigma$ band allowing $L_{X}$ to vary within the uncertainty shown in Fig. 1 and fixing $D_{L}$ to the values shown in Tab. 1 (orange band). Right: the same as on left but for NGC 3079. The SED is taken from [41].

5 Diffuse Spectrum

Equipped with the previous results, in this section, we extrapolate the information to the whole population constraining the diffuse neutrino fluxes of seyfert galaxies.

For the distribution, we follow the distribution of Ref. [47] (see [20] for further details) and define the comoving density of sources as

$$\rho(L_{X},z)=\frac{d\Phi_{X}(L_{X},z)}{d\ logL_{X}}=\frac{d\Phi_{X}(L_{X},0)}{d\ logL_{X}}e(z,L_{X})$$ (5.1)

with

$$\frac{d\Phi_{X}(L_{X},0)}{d\ logL_{X}}=A\bigg{[}\bigg{(}\frac{L_{X}}{L_{*}}\bigg{)}^{\gamma_{1}}+\bigg{(}\frac{L_{X}}{L_{*}}\bigg{)}^{\gamma_{2}}\bigg{]}^{-1}$$ (5.2)

and

$$e(z,L_{x})=\begin{cases*}(1+z)^{p_{1}}&if $z\leq z_{c_{1}}(L_{X})$\\ (1+z_{c_{1}})^{p_{1}}\big{(}\frac{1+z}{1+z_{c_{1}}}\big{)}^{p_{2}}&if $z_{c_{1}}(L_{X})\leq z\leq z_{c_{2}}(L_{X})$\\ (1+z_{c_{1}})^{p_{1}}\big{(}\frac{1+z_{c_{2}}}{1+z_{c_{1}}}\big{)}^{p_{2}}\big{(}\frac{1+z}{1+z_{c_{2}}}\big{)}^{p_{3}}&if $z\geq z_{c_{2}}$\end{cases*}$$ (5.3)

where

$$p_{1}(L_{X})=p^{*}_{1}+\beta_{1}(\ logL_{X}-44)$$ (5.4)

while

$$z_{c_{1}}(L_{X})=\begin{cases*}z_{c_{1}}^{*}\big{(}\frac{L_{X}}{L_{a_{1}}}\big{)}^{\alpha_{1}}&if $L_{X}\leq L_{a_{1}}$\\ z_{c_{1}}^{*}&otherwise\end{cases*}$$ (5.5)

and

$$z_{c_{2}}(L_{X})=\begin{cases*}z_{c_{2}}^{*}\big{(}\frac{L_{X}}{L_{a_{2}}}\big{)}^{\alpha_{2}}&if $L_{X}\leq L_{a_{2}}$\\ z_{c_{2}}^{*}&otherwise\end{cases*}$$ (5.6)

Tabs. 2 and 3 summarise all the parameters used in the distribution (see [47] for further details).

$A\small{(10^{-6}\,h_{70}^{3}\,\rm Mpc^{-3})}$	$\log L_{*}$	$\gamma_{1}$	$\gamma_{2}$	$p_{1}^{*}$	$\beta_{1}$
$2.91\pm 0.07$	$43.97\pm 0.06$	$0.96\pm 0.04$	$2.71\pm 0.09$	$4.78\pm 0.16$	$0.84\pm 0.18$

Table 2: Summary of the distribution parameters taken from [47]. All the luminosity are expressed in units of $\rm erg\,\rm s^{-1}$.

$z_{c_{1}}^{*}$	$\log L_{a_{1}}$	$\alpha_{1}$	$p_{2}$	$p_{3}$	$z_{c_{2}}^{*}$	$\log L_{a_{2}}$	$\alpha_{2}$
$1.86\pm 0.07$	$44.61\pm 0.07$	$0.29\ \pm 0.02$	-1.5	-6.2	3.0	44	-0.1

Table 3: Continuation of Tab. 2 for the distribution parameters. All the luminosity are expressed in units of $\rm erg\,\rm s^{-1}$.

Along with Ref. [20], we multiply the A parameter for 1.5 in order to take into account the fraction of compton-thick AGNs into the x-ray luminosity function. The final diffuse neutrino spectrum per solid angle reads

$$\Phi_{\nu}(E_{\nu},\gamma,\eta,E_{\max})=\frac{c}{4\pi H_{0}}\int_{0}^{z_{\rm max}}\frac{dz}{E(z)}\int_{10^{41}}^{10^{47}}\rho(L_{X},z)Q_{\nu+\bar{\nu}}(E_{\nu}(1+z),L_{X},\gamma,\eta,E_{\rm max})$$ (5.7)

with $E(z)=\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}$ and the $L_{X}$ integration limits are set to be $10^{41}\,\rm erg\,\rm s^{-1}$ and $10^{47}\,\rm erg\,\rm s^{-1}$. Fig. 5 shows our final best-fit scenario (golden dashed line), the $1\sigma$ band according to our stastical analysis fixing the source distribution density to the best-fit scenario (dark red band) and finally the maximal $1\sigma$ band taking also into account the uncertainty on the density distribution (orange band). We compare the results with the latest diffuse IceCube data (6 year cascade [8] and 10 year of starting tracks [10]) and the expected KM3NeT differential sensitivity after 10 years of full operation considering all-sky shower events [48].

Figure 5: best-fit diffuse neutrino flux from seyfert galaxies (golden dashed line), the $1\sigma$ band fixing the source density distribution to the best-fit scenario (dark red scenario) and the total $1\sigma$ band considering an uncertainty over the density distribution (orange band) compared compared with the diffuse IceCube fluxes (10 years of starting tracks [10] and 6 years of cascades flux [8]) and with the expected KM3NeT differential sensitivity for shower events [48].

We notice that if we assume that all the seyfert galaxies have the same relation between $L_{\nu}$ and $L_{X}$ as imposed by the four sources used in the statistical analysis, we overproduce the diffuse neutrino flux at $\sim 1-10\,\rm TeV$. On this regard, the latest diffuse flux from IceCube starting tracks strongly constrain the astrophysical diffuse neutrino spectrum at $\sim 1\,\rm TeV$ [10]. We stress that this result is independent on the possible energetics tension given by the low $\beta$ plasma parameter because $L_{\nu}$ and $L_{X}$ are data-driven and any model should be able to reproduce these, leading to almost the same result. This has already been emphasised by Ref. [32] where the authors has evaluated the cosmic neutrino background from non-jetted AGNs, assuming that all AGNs behave like NGC 1068. Our best-fit result is slightly higher than the one reported in Ref. [32] because our model self-consistently exploits the luminosity information of all the four sources providing an excess into IceCube data rather than assume all the sources to be NGC 1068-like. However, within all uncertainties the results are consistent. The CR cut-off at $200\,\rm TeV$ naturally suppresses the contribution of seyfert galaxies at $E_{\nu}\geq 100\,\rm TeV$, leaving room for other astrophysical sources such blazars (see [32, 58] for instance) or starburst galaxies [22, 26, 59]. However, the chosen cut-off is only driven by the point-like IceCube observations which sets an upper limit of neutrino emission at $\sim 10-50\,\rm TeV$ rather than first-principle calculations. Therefore, future theoretical studies directly probing the high-energy CR cut-offs in these sources are fundamental. In fact, a higher $E_{\rm max}$ value would lead to a higher diffuse neutrino flux at $E_{\nu}\sim 100\,\rm TeV$, providing further constraints to seyfert galaxies emissions. We also highlight that the upcoming KM3NeT telescope [48] will be fundamental in order to probe the role of each of these components into the high-energy neutrino sky especially in the $1-10\,\rm TeV$ energy range. Indeed, shower-like events, having a reduced background rate will be a perfect sample to investigate the diffuse emission of these galaxies with.

6 Conclusions

In this paper, we have revisited the constraints on the neutrino emission of seyfert galaxies, initially proposed by V.S. Berezinsky [2], exploiting the latest observations by the IceCube collaboration [43, 42, 44, 45, 41]. We develop a theoretical model accounting both for CR escape and energy losses mechanisms, assuming equi-ripartition between CR and magnetic energy densities. The four point-like sources providing an excess into IceCube data are consistent with an equi-ripartion hypothesis but the beta plasma parameter required to fit the data is rather low $(\beta\lesssim 1)$, pointing to a high energetics carried by CRs inside these sources. Therefore, future dedicated analyses will be crucial to investigate if such condition can be met in environment of seyfert galaxies. Furthermore, extrapolating the neutrino emission to the whole source population might overestimate the diffuse neutrino flux at $\sim 1-10\,\rm TeV$ energies leading to a potential tension in the high-energy neutrino sky. We might argue that all the sources might not be in an equi-ripartition state and if we take a back-of-the-envelope estimate of $\sim 1/10$ of the sources (considering that IceCube has observed an excess for 3 sources out of a catalogue of the 30 most luminous sources in the northern hemisphere), this would sensibly reduce the tension (see also Ref. [32] for further remarks). All in all, all these results point to the fact that regions sourrounding the SBMHs in AGNs might sensibly produce high-energy neutrinos.

Appendix A Details on the Dynamical Timescales

In this section, we report details on the dynamical timescale in the coronae. The escape timescale is given by the in-fall timescale onto the black hole $\tau_{\rm esc}=R/V_{\rm fall}$, with $V_{\rm fall}=\alpha\sqrt{2GM_{\rm BH}/R}$. Therefore $\tau_{\rm esc}$ reads

$$\tau_{\rm esc}=3.54\cdot 10^{6}\bigg{(}\frac{\alpha}{0.1}\bigg{)}\bigg{(}\frac{r}{20}\bigg{)}^{1/2}\bigg{(}\frac{M_{\rm BH}}{10^{8}\,\rm M_{\odot}}\bigg{)}\,[s]$$ (A.1)

we fix the friction coefficient $\alpha=0.1$ [46]. For the mass of the SMBHs, we make use of [20, 60]

$$M_{\rm BH}=2\cdot 10^{7}\,M_{\odot}\bigg{(}\frac{L_{X}}{1.155\cdot 10^{43}\,\rm erg\,\rm s^{-1}}\bigg{)}^{0.746}$$ (A.2)

in this way, we can express $M_{\rm BH}$ in terms of $L_{X}$. $\tau_{\rm loss}$ is given by the competition between inverse compton (IC), Synchrotron, Bethe-Heitler pair production, photomeson production and pp collisions. For IC and Synchrotron the timescales read [20]

$$\tau_{\rm IC,\rm syn}(E)=\frac{3}{4}\bigg{(}\frac{m_{p}}{m_{e}}\bigg{)}^{3}\frac{m_{e}c^{2}}{c\sigma_{T}U_{\rm ph,B}}\bigg{(}\frac{E}{m_{p}c^{2}}\bigg{)}^{-1}$$ (A.3)

where $U_{\rm ph}=L_{\rm ph}/(4\pi R^{2}c)$ is the total background photon background with $L_{\rm ph}$ is the bolometric luminosity. $U_{B}=B^{2}/8\pi$ is the magnetic energy density. The timescale for Bethe-Heitler pair production reads [20]

$$\tau_{\rm BH}^{-1}=\frac{7(m_{e}c^{2})^{3}\alpha_{f}\sigma_{T}\cdot c}{9\sqrt{2}m_{p}c^{2}}\bigg{(}\frac{E}{m_{p}c^{2}}\bigg{)}^{-2}\int_{m_{e}m_{p}c^{4}/E}^{+\infty}d\epsilon\frac{n_{\rm ph}(\epsilon)}{\epsilon^{3}}\cdot\bigg{[}\bigg{(}\frac{2E\cdot\epsilon}{m_{e}m_{p}c^{4}}\bigg{)}^{3/2}\bigg{(}\rm log\bigg{(}\frac{2E\cdot\epsilon}{m_{e}m_{p}c^{4}}\bigg{)}-\frac{2}{3}\bigg{)}+\frac{2}{3}\bigg{]}$$ (A.4)

where $n_{\rm ph}(\epsilon)$ is the background photon density from Ref. [20]. For photomeson interaction, we use [61]

$$\tau_{p\gamma}^{-1}=\frac{c^{5}m_{p}^{2}}{2E^{2}}\int_{0}^{+\infty}dE^{{}^{\prime}}\frac{n_{\rm ph}(E^{{}^{\prime}})}{E^{{}^{\prime}3}}\int_{E_{\rm th}}^{2\frac{E\cdot E^{{}^{\prime}}}{m_{p}c^{2}}}d\epsilon\epsilon\sigma_{p\gamma}(\epsilon)K_{p\gamma}(\epsilon)$$ (A.5)

where $\sigma_{p\gamma}(\epsilon)$ is total inelastic $p\gamma$ cross section in the proton rest frame [62]. $K_{p\gamma}$ is approximated with a step function [61]

$$K_{p\gamma}(\epsilon)=\begin{cases*}0.2&if $\epsilon<1\,\rm GeV$\\ 0.6&otherwise\end{cases*}$$ (A.6)

finally, $E_{\rm th}=145\,\rm MeV$ is the kinetic threshold for the process. The pp timescale reads [20, 52]

$$\tau_{pp}=\frac{1}{k_{p}n_{p}\sigma_{pp}c}$$ (A.7)

where $k_{p}=0.5$ is the mean inelasticity of the process. $n_{p}$ is the gas density and it is set considering charge neutrality of the corona as [20]

$$n_{p}=\frac{\tau_{\tau}}{R\sigma_{T}}\simeq 2.8\cdot 10^{9}\frac{\tau_{\tau}}{1.1}\bigg{(}\frac{r}{20}\bigg{)}^{-1}\bigg{(}\frac{M_{\rm BH}}{10^{8}\,\rm M_{\odot}}\bigg{)}^{-1}\,[\rm cm^{-3}]$$ (A.8)

where $\tau_{\tau}=1.1$ is the Thompson scattering optical depth [20]. We can analytically estimate the efficiency to produce neutrinos as $F_{\rm cal}\simeq(\tau_{\rm loss}^{-1}+\tau_{esc}^{-1})^{-1}/\tau_{pp}$, $F_{\rm cal}\simeq(\tau_{\rm loss}^{-1}+\tau_{esc}^{-1})^{-1}/\tau_{p\gamma}$ respectively for $pp$ and $p\gamma$ interactions [20].

Fig. 6 shows the efficiency for three different values of $L_{X}$, $10^{42}\,\rm erg\,\rm s^{-1}$, $10^{44}\,\rm erg\,\rm s^{-1}$, $10^{46}\,\rm erg\,\rm s^{-1}$ as a function of the CR energy, fixing $\beta=0.32$ as the best-fit scenario obtained in the main text.

Figure 6: efficiency for neutrino production as a function of the CR energy. The blue line corresponds to the $pp$ collision channel, the orange line corresponds to $p\gamma$ interactions and the green line corresponds to the calorimetric limit. We also report a dashed red vertical line corresponding to the CR cut-off assumed in the analysis. The top left, top right and the bottom panels respectively corresponds to $L_{X}=10^{42}\,\rm erg\,\rm s^{-1}$, $L_{X}=10^{44}\,\rm erg\,\rm s^{-1}$ and $L_{X}=10^{46}\,\rm erg\,\rm s^{-1}$.

The results show that for CR energy below the cut-off, the neutrino production is dominated by pp collisions.