IceCube Neutrinos from Hadronically Powered Gamma-Ray Galaxies

TL;DR

This study investigates whether gamma-ray galaxies powered by cosmic ray interactions can explain the high-energy neutrino flux observed by IceCube, suggesting they could be the dominant sources above 100 TeV.

Contribution

It demonstrates that hadronically powered gamma-ray galaxies can account for the majority of IceCube's high-energy neutrino observations within current gamma-ray background constraints.

Findings

Hadronically powered gamma-ray galaxies can explain IceCube's neutrino flux above 100 TeV.

A spectral index of ≤ 2.12 is compatible with gamma-ray background measurements.

These sources could account for about half of the neutrino flux between 10-100 TeV.

Abstract

In this work we use a multi-messenger approach to determine if the high energy diffuse neutrino flux observed by the IceCube Observatory can originate from $\gamma$-ray sources powered by Cosmic Rays interactions with gas. Typical representatives of such sources are Starburst and Ultra-Luminous Infrared Galaxies. Using the three most recent calculations of the non-blazar contribution to the extragalactic $\gamma$-ray background measured by the Fermi-LAT collaboration, we find that a hard power-law spectrum with spectral index $\alpha \leq 2.12$ is compatible with all the estimations for the allowed contribution from non-blazar sources, within 1$\sigma$. Using such a spectrum we are able to interpret the IceCube results, showing that various classes of hadronically powered $\gamma$-ray galaxies can provide the dominant contribution to the astrophysical signal above 100 TeV and about half of the contribution to the energy flux between 10-100 TeV. With the addition of neutrinos from the Galactic plane, it is possible to saturate the IceCube signal at high energy. Our result shows that these sources are still well motivated candidates.

Figures

Tables

Table 1. Table 1 : Decomposition of the 5.7 years of the 82 HESE events [ 35 ] into different source components. The background events are given by IceCube in [ 33 ] , whereas the Galactic neutrinos are computed in [ 19 ] and they are roughly 1 per year in the HESE dataset. Neutrinos from HAGS accounts for 33 of the 41 signal events, i.e. about 80% of the astrophysical signal.

	Expected	Fraction
Atmospheric muons	$25.2\pm 7.3$	$32\%\pm 9\%$
Atmospheric neutrinos	${15.6}_{-3.9}^{+11.4}$	${20\%}_{-5\%}^{+14\%}$
Extragal. neutrinos (HAGS)	$33\pm 8$	$41\%\pm 10\%$
Galactic neutrinos	$\le$ 5.7	$\le 7\%$
Total	${79.5}_{-11.7}^{+15.9}$	–

Table 2. Table 2 : Summary of our main results. In the table the values of different parameters are reported, namely: the energy cutoff E cut ⋆ superscript subscript 𝐸 cut ⋆ E_{\text{cut}}^{\star} in units of PeV, the local source density ρ 0 ​ ( α ) subscript 𝜌 0 𝛼 \rho_{0}(\alpha) in units of Mpc − 3 superscript Mpc 3 \mbox{Mpc}^{-3} , the contribution to the extragalactic γ 𝛾 \gamma -ray background as a function of the spectral index α 𝛼 \alpha . The parameter ρ 0 ​ ( α ) subscript 𝜌 0 𝛼 \rho_{0}(\alpha) denotes the local density of sources that is required to power the TGM energy flux. Here we show two different examples, assuming as a prototype source NGC 253 and ARP 220. In the last three columns on the right we write what is the tension between the contribution to the EGB expected from our model and the results of the three calculations discussed in the paper [ 41 , 42 , 43 ] .

		NGC 253	ARP 220	Contribution to		Tension with
$\alpha$	$E_{\text{cut}}^{\star }$	${\rho }_0\left(\alpha \right)$	${\rho }_0\left(\alpha \right)$	total EGB	Fermi coll. [41]	Lisanti et al. [42]	Zechlin et al. [43]
1.9	3.2	$7.7\times {10}^{-5}$	$5.1\times {10}^{-7}$	6%	no	no	no
2.0	4.1	$2.3\times {10}^{-4}$	$1.6\times {10}^{-6}$	11%	no	no	no
2.1	5.1	$7.9\times {10}^{-4}$	$5.3\times {10}^{-6}$	25%	0.8$\sigma$	no	0.3$\sigma$
2.12	5.3	$1.0\times {10}^{-3}$	$6.7\times {10}^{-6}$	28%	1$\sigma$	no	0.5$\sigma$
2.2	6.3	$2.8\times {10}^{-3}$	$1.8\times {10}^{-5}$	50%	2.6$\sigma$	2.1$\sigma$	1.7$\sigma$
2.3	7.8	$1.1\times {10}^{-2}$	$7.2\times {10}^{-5}$	100%	6.9$\sigma$	8.7$\sigma$	4.8$\sigma$

Equations

$$\begin{array}{rl}\frac{d{\Phi }_{\nu }(E)}{dE}& =\int \mathrm{d}z\frac{\rho (z)}{H(z)}\frac{d{\varphi }_{\nu }}{dE}(E(1+z))\\ & =\int \mathrm{d}z\frac{\rho (z)}{H(z)}\frac{3}{2}\frac{d{\varphi }_{\gamma }}{dE}(E(1+z))\\ \approx {\Phi }_0& \left(\frac{{\rho }_0}{{\rho }_0(\alpha )}\right)\left(\frac{{\mathcal{L}}_{\text{100 GeV}}}{{\mathcal{L}}_{\text{100 GeV}}^{\text{NGC253}}}\right){\left(\frac{E}{\text{100 TeV}}\right)}^{-{\alpha }^{\star }}{e}^{-E/E_{\text{cut}}^{\star }(\alpha )},\end{array}$$

$${\rho }_0{L}_{\gamma }^{0.1-100}=2\times 10^{44}\mathrm{erg}\mathrm{Mp}{\mathrm{c}}^{-3}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\mathcal{R}}_{10-100}\equiv \frac{{\int }_{10\text{TeV}}^{100\text{ TeV}}\frac{d{\Phi }_{\nu }}{dE}dE}{{\int }_{10\text{ TeV}}^{100\text{ TeV}}\frac{d{\Phi }_{\text{CAS4}}}{dE}dE}$$

$${\int }_0^{\infty }\frac{d{\varphi }_{\gamma }}{dE}EdE\approx \frac{2}{3}{\int }_0^{\infty }\frac{d{\varphi }_{\nu }}{dE}EdE$$

$$\frac{d{\varphi }_{\nu }}{dE}(E,E_{\text{cut}},\alpha ,z)={\int }_0^z\mathrm{d}z\frac{\frac{d{\mathcal{L}}_{\nu }}{dE}[E(1+z),E_{\text{cut}},\alpha ]}{4\pi D_c(z{)}^2(1+z{)}^2{E}^2}$$

$$\frac{\mathrm{d}N}{\mathrm{d}V}(z)={\rho }_0{\left[{\left(1+z\right)}^{a\eta }+{\left(\frac{1+z}{B}\right)}^{b\eta }+{\left(\frac{1+z}{C}\right)}^{c\eta }\right]}^{1/\eta },$$

$$\frac{d{\Phi }_{\nu }}{d{E}_{\nu }}(E,E_{\text{cut}},\alpha )={\int }_0^{z_{\text{max}}}\frac{\mathrm{d}{\varphi }_{\nu }}{\mathrm{d}E}(E,E_{\text{cut}},\alpha ,z)\frac{\mathrm{d}N}{\mathrm{d}V}(z,\alpha )\frac{\mathrm{d}V}{\mathrm{d}z}dz$$

Full text

1 Introduction

The IceCube Collaboration has detected a diffuse high-energy astrophysical neutrino flux [1]. However, despite the six years since its detection, its origin remains unclear. Several candidate source populations have been proposed, such as blazars [2, 3, 4], $\gamma$-Ray Bursts (GRBs) [5, 6, 7, 8], Star Forming Galaxies and Starburst Galaxies [9, 10, 11, 12], dark matter decay [13, 14] or Galactic sources, like Galactic center [15, 16, 17], Galactic plane [18, 19, 20], Galactic halo [21]. Some models are strongly constrained by the absence of correlations between the directions of the high-energy neutrinos and known sources. This result is compatible with neutrino-bright sources being dim photon sources. For example, distant blazars that are not resolved by the Fermi satellite due to their low luminosity, could have a high fraction of interacting protons and thus be efficient neutrino sources [22]. GRBs with chocked jets are another type of $\gamma$-ray dim and neutrino-bright sources [23]. So far, only one neutrino source candidate has been tentatively identified, the blazar TXS 0506+056 [24], one of the brightest blazars located at a redshift $z=0.34$.

The multi-messenger approach can be applied in various ways, either on a “per source” basis where a single object is studied simultaneously in multiple wavelengths, in neutrinos and in gravitational waves. Alternatively, one can also study the entire populations of sources using diffuse (time and angular integrated) fluxes of neutrinos and photons, where the relative differences between observed energy spectra can constrain certain scenarios [25]. The latter approach has been recently applied to study the contribution of Star Forming Galaxies (SFG), including Starburst and Star-Forming Active Galactic Nuclei, to the diffuse neutrino flux [26]. The study concluded that these sources can only contribute at most a few percent of the observed neutrinos, essentially excluding SFGs as the dominant sources of high-energy neutrinos in IceCube. A similar conclusion is also reached in [27, 28]. Since the $\gamma$-ray radiation in SFGs predominantly originates from Cosmic Ray interactions with gas, the result may be generalized to the extent that proton-proton ($pp$) interactions can not be the main mechanism of astrophysical neutrino production.

In astrophysical environments where pp interactions dominate cosmic ray cooling, the relation between the $\gamma$-ray and the neutrino emission is fixed. With approximately equal numbers of $\pi^{+}$, $\pi^{-}$ and $\pi^{0}$ being produced in high-energy collisions, the total energy budget of $\gamma$-rays is 2/3 of the total energy budget of neutrinos. If these “hadronic sources” power the entire diffuse neutrino flux, it is sufficient that only a subset of all galaxies significantly contributes, in line with the requirement for a high-energy cutoff of $E_{\nu}\sim$ few PeV. Such sources are likely to carry characteristics of Starburst galaxies (SBG) such as enhanced stellar light and dust production, which results in infrared luminosities 10-100 times higher compared the more abundant normal galaxies [29]. In reality, this heightened level of infrared galaxies is a continuum, with even higher values observed in what are classified as Ultra-Luminous Infrared Galaxies (ULIRGs), which are more luminous (infrared luminosity 100-1000 higher than normal Galaxies, see Fig.1 of [30]) but less numerous. Fermi reveals an almost linear correlation between the infrared and the $\gamma$-ray luminosities [30, 31]. Therefore, “infrared bright” galaxies like Starburst and ULIRGs are good candidates of HAdronically powered $\gamma$-ray GalaxieS (HAGS). The neutrino and $\gamma$-ray spectrum from these objects is expected to be a power law, whose spectral index remains uncertain in the relevant energy range (above 50 GeV), since only a few members of the broad HAGS source class (NGC253, M82) have been detected up to very high energies. We will therefore consider the spectral index as a free parameter in the present work, keeping in mind that the observation of NGC253 suggests an $\sim E^{-2.15\pm 0.10}$ above 50 GeV, which is the energy region relevant for our purpose (see Fig. 2).

This work aims to re-evaluate the compatibility of the most recent neutrino observations by IceCube with an origin from HAGS. We relax the hypothesis that the entire astrophysical neutrino flux is produced by a single source class, focusing first on the high energy part (above 100 TeV). Then we apply a multi-messenger method and combine the diffuse extragalactic $\gamma$-ray background (EGB) observations by Fermi with the throughgoing muon energy flux measured by IceCube. After that we evaluate the agreement with other two IceCube dataset, where low energy events (down to $\sim$ TeV) are contained.

2 Methods

We start with a comparison of the three IceCube data analyses, namely: i) the through-going muons (TGM) [32, 33] originating from muon neutrino and anti-neutrino interactions outside the detector. This selection contains tracks from the opposite (Northern) hemisphere and at a higher energy threshold of 200 TeV. In [33] the throughgoing muon flux is fitted by using a power-law flux $\propto N\times E^{-\alpha}$, with $\alpha=2.2\pm 0.1$ and $N=1.01^{+0.26}_{-0.23}$, where $N$ is the single flavor normalization at 100 TeV in units of $10^{-18}\ \rm GeV^{-1}cm^{-2}s^{-1}sr^{-1}$;

ii) the high-energy starting events (HESE) [34, 35], characterized by an interaction vertex contained in the fiducial detector volume. HESE contain shower and track like events mostly coming from the Southern hemisphere. Concerning the HESE flux we refer to [33], in which it has been found $\alpha=2.9\pm 0.3$ and $N=2.5\pm 0.8$ for this dataset, in the same units defined above;

iii) the 4-year cascade (CAS4) [36] sample contains neutrinos of all flavors interacting with a cascade topology, i.e. mostly electron and tau neutrinos. The energy threshold is the lowest for this dataset, around $\sim$ TeV energies. The flux associated to CAS4 [33] is fitted by using $\alpha=2.48\pm 0.08$ and $N=1.57^{+0.23}_{-0.22}$.

We notice that above 200 TeV both throughgoing muons and HESE suggest a hard power law spectrum $d\Phi_{\nu}/dE\propto E^{-2.2\pm 0.1}$. Since above this energy the contribution coming from the atmospheric background is small ($\sim 20\%$ using the signalness reported in Tab.4 of [32]), we assume that this flux is representative of the true astrophysical signal. On the other hand one cannot neglect the information coming from the spectrum measured below 200 TeV; therefore, at the end of our multi-messenger analysis, we subsequently evaluate the compatibility between our findings with low energetic HESE and CAS4 datasets.

In order to compute the diffuse neutrino and $\gamma$-ray fluxes we make use of the relation, in which the energy budget in neutrinos is 3/2 of that in $\gamma$-rays. We assume for the density of the hadronic sources an evolution $\rho(z)=\rho_{0}(1+z)^{m}$, with $m=3.4$ for $z<1$ and $m=-0.5$ up to $z=4$, like the star forming rate [37]. The diffuse all-flavor neutrino flux $\phi_{\nu}$ is, therefore, related to the product of the local source density $\rho$ and the luminosity at 100 GeV $\mathcal{L}_{\text{100 GeV}}$, i.e. to the local emissivity:

[TABLE]

where $\Phi_{0}$ is constant and equal to $3\times 10^{-18}\ \rm GeV^{-1}\ cm^{-2}\ s^{-1}\ sr^{-1}$, $\mathcal{L}_{\text{100 GeV}}^{\text{NGC 253}}=(5\pm 2)\times 10^{39}\text{erg/s}$ is the $\gamma$-ray luminosity of NGC253 at 100 GeV [38], which we use here as a HAGS prototype. We do not use any evolution of the $\gamma$-ray luminosity function, but instead assume a single prototype object for our calculation, which evolves only with an “effective” number density, such that the luminosity density evolves according to the Star Formation Rate. We determine the spectral-index dependent local ($z=0$) source density $\rho_{0}(\alpha)$ by fitting the astrophysical neutrino signal to the TGM spectrum in the range $0.1-1$ PeV. More details on the calculation are reported in the supplementary material. We assume throughout that the spectrum of neutrinos and $\gamma$-rays at the sources is described by a power-law with a spectral index $\alpha$ and an exponential cutoff energy $E_{\text{cut}}$. The spectral index remains approximately the same at Earth while the energy cutoff is slightly different and we denote it with $E^{\star}_{\text{cut}}$. The values of $\rho(\alpha)$ and $E_{\text{cut}}^{\star}(\alpha)$, given in Tab. 2, are discussed in the result section.

In the computation of the diffuse $\gamma$-ray flux, the interactions with the Extragalactic Background Light (EBL) have to be taken into account. The high-energy photons are absorbed and reprocessed through electromagnetic cascades, resulting in the re-appearance of this energy flux at lower energies. We use the same source spectrum as for the neutrinos, offset by 2/3 coming from energy budget considerations (see Eq. (A.1)). The transport of gamma rays through the intergalactic medium yields two components at the observer $\phi_{\gamma}=\phi_{\gamma}^{\text{dir}}+\phi_{\gamma}^{\text{casc}}$. The direct component here arrives from the source population, and is attenuated at high energies by the EBL. This attenuation feeds electromagnetic cascades, giving rise to the cascade component at energies below a few hundred GeV. We use [39] for the EBL model and the method given in [40] to compute the electromagnetic cascade during the propagation.

The resulting propagated $\gamma$-ray flux has to be compared with the non blazar component of the extragalactic gamma-ray background (EGB). This residual component of the EGB was estimated by the Fermi collaboration [41]. The remaining (non-blazar) fraction, is thought to be shared by all the other $\gamma$-ray emitters, such as normal and starburst galaxies and misaligned blazars. The share between the blazars and the non-blazar contribution varies between the different analyses, and are accompanied by large errors: the Fermi collaboration identifies the contribution from blazars to the EGB above 50 GeV as $86\%^{+16}_{-14}\%$ [41]; Lisanti et al. [42] as $68\%^{+9\%}_{-8\%}$ and Zechlin et al. [43] as $81\%^{+52\%}_{-19\%}$. Complementary to this measurement, the contribution from HAGS to the extragalactic $\gamma$-ray background has been discussed in [44, 45, 46]. The cumulative $\gamma$-ray flux at Earth $\phi_{\gamma}(\alpha)$ non-trivially depends on the spectral index of the source prototype since the share between the direct $\phi_{\rm dir}(\alpha)$ and the cascade $\phi_{\rm casc}(\alpha)$ components change. We evaluate the compatibility of HAGS with the non-blazar contribution by comparing the $\gamma$-ray flux integrated above 50 GeV with the integral of the EGB flux in the same energy range.

3 Results

We check the compatibility of our model with the TGM energy flux by scanning over the source spectral index $\alpha$ and by performing fits of the local density $\rho_{0}(\alpha)$ of HAGS that maximize the throughgoing muon energy flux. We remark that the throughgoing muon sample is less contaminated by atmospheric muons and the 200 TeV energy threshold reduces the contribution of atmospheric neutrinos. Moreover no contamination from the Galactic plane is expected.

In Fig.1 we report our result for different spectral indices and different luminosity densities. The bands represent forbidden regions at 1$\sigma$. The gray region indicates the exclusion of the model by its contribution to the most pessimistic non-blazar estimation of the EGB. The orange band limits the allowed neutrino energy flux integrated between 200 TeV and 5 PeV with respect to what is maximally allowed by the TGM measurements. The dashed green lines represent partial contributions to the TGM energy flux. The result is valid for generalized $pp$ sources, since the luminosity density is given by the product between the source density at redshift $z=0$ and the luminosity of the sources at 100 GeV. More details are reported in Tab. 2 of the supplementary material. For harder spectrum and the two HAGS prototypes, the starburst galaxy NGC253 and the ULIRG Arp220 (with a luminosity 100-150 times higher than NGC [30]), we find the local densities (given in Tab.2) to be compatible with the expectation for such sources. Concerning Arp220, we remark that although this object is expected to be close to the calorimetric regime and ideal for neutrino production [47, 48], present measurements indicate a spectrum softer than $E^{-2}$ [49]. In addition to this, there are also works discussing the possibility that this object is contamined by AGN activity [12, 50, 51], rather than a pure Starburst Galaxy.

As demonstrated in the Supplementary Material, normal galaxies contribute very little to the diffuse $\gamma$-ray and the neutrino fluxes due to their soft spectral indices and low luminosities. Our result places an upper limit on the source spectral index, with values of $\alpha>2.3$ being excluded at 5$\sigma$, thus confirming the findings by [26] who based their arguments on the soft single power-law fit to the HESE data available at that time. On the contrary, we find that a source spectral index $\alpha\leq 2.12$ is found to be compatible with the three estimations of the non blazar contribution within $1\sigma$. Moreover an index of $\alpha=2.12$ allows the saturation of the throughgoing muon energy flux. The first constraint on the spectral index of the neutrino flux is reported in [25], where the authors find $\alpha<2.1-2.2$. Similar results were also reached in [52, 53, 54]. Our result slightly improves the previous constraints, although the uncertainty in the redshift evolution of the sources is not included and this may have some impact on the spectral index constraint, since the source evolution affects the diffuse $\gamma$-ray flux associated to neutrinos and, as a consequence, the comparison with the EGB.

We note here that this discussion on the contribution to the non-blazar of the EGB neglects the inevitable contribution from misaligned AGN. An estimation of the level of this component is in the range between 4% and 40% [55], with a best fit value of 12%. Taking this into account the remaining non AGN contribution would sit at level of 16%. Even for such a reduced level we still find compatibility with the throughgoing muon energy flux using a hard spectrum ($\alpha=2$) (see Appendix, Tab.2). For the sake of completeness, we point out that in [56] it is shown that a hard spectrum would be the best scenario to reconcile the IceCube neutrinos above 100 TeV in the context of star-forming galaxies modeling, although it is not possible to saturate the flux under this hypothesis. This work is based on new theoretical model to predict the luminosity and spectrum of gamma-ray and neutrino emission from star-forming galaxies. The authors find that an injected $E^{-2}$ spectrum is able to power $\sim 22\%$ of the HESE events and $\sim 50\%$ of the throughgoing muon energy flux. An additional contribution from AGN-Starburst, rather than pure Starburst Galaxies, may help to saturate the observed astrophysical neutrino flux [26]. In addition to, it has been discussed in the literature the possibility that high energy neutrinos can be produced by Galaxy Clusters [25, 57] and by Radio Galaxies [58, 59].

Concerning the $\gamma$-ray spectrum from NGC253, it was expected to be close to $\sim E^{-2.3}$ from theoretical arguments [60]. However the update measurements provided by Fermi and HESS [38] suggest an $E^{-2.15\pm 0.10}$ spectrum above 50 GeV, i.e. our energy region of interest (see Fig.2). Therefore we rely on the observations, instead of considering theoretical predictions for the spectrum. Let us notice that such a hard spectrum can be explained from the theoretical point of view, assuming a small diffusion coefficient [61]. We re-fitted the data above 50 GeV, since we are particularly interested in the high-energy spectrum. The spectral indices obtained by Fermi and HESS [38] are $\alpha=2.09\pm 0.07^{\text{stat}}\pm 0.05^{\text{sys}}$ above 60 MeV and $\alpha=2.22\pm 0.06^{\text{stat}}$ above 3 GeV, which both are compatible with our fit above 50 GeV. Therefore as a baseline spectrum we choose $\alpha=2.12$, i) being the softest spectrum that is able to saturate the TGM energy flux, ii) without producing any tension with the non blazar contribution, iii) being compatible with the recent update on the observation result of NGC 253. We also remark that a spectral index $\alpha=2.12$ is perfectly consistent with the theoretical prediction presented almost 15 years ago in [10], in which a spectral index $\alpha=2.15\pm 0.10$ was expected for this kind of objects.

Moreover we also checked the compatibility of the local density resulting from our paper with the measured one. The infrared luminosity of NGC253 corresponds to $2\times 10^{10}L_{\odot}$ (see Fig. 1 of [30]), where $L_{\odot}$ is the luminosity of the Sun. For Starburst Galaxies with this infrared luminosity the local density is $1\times 10^{-4}\rm\ Mpc^{-3}$ according to Fig. 10 of [62], i.e. 10 times smaller than that adopted in our benchmark model (see Tab. 2). However, one has to consider various elements: i) this value depends on the spectral index; for example, the local density required for an $E^{-2}$ spectrum is $2\times 10^{-4}\rm\ Mpc^{-3}$ (see Tab. 2); ii) Starburst Galaxies are not the only sources contributing to the neutrino flux, since HAGS denote a larger class of objects, in which also Star Forming - AGN with Starburst characteristics are included. These objects are likely to be close to calorimetric, providing a hard neutrino spectrum. According to Tab. 1 and Tab. 2 in [12], the local density of Star Forming - AGN (SB) is $1.5\times 10^{-4}\rm\ Mpc^{-3}$; iii) the relevant quantity is not the local density of the considered benchmark object, but the product between the local density and the $\gamma$-ray luminosity. In our benchmark model, with $\rho_{0}=10^{-3}\times\rm\ Mpc^{-3}$ and $L_{\gamma}^{0.1-100}=6\times 10^{39}\rm erg/s$ in the energy range 0.1-100 GeV (see Fig. 1 of [30]), the required local luminosity density is

[TABLE]

to power the throughgoing muon flux. This value can be compared with Ref. [63], in which the luminosity function of Starburst galaxies has been explicitly included. From Fig. 3 of [63] relates a source density $\rho_{0}=10^{-5}\rm\ Mpc^{-3}$ to the luminosity $L_{\nu_{\mu}}^{10^{4}-10^{7}}=3\times 10^{40}~{}{\rm erg~{}s}^{-1}$ in the energy range $10^{4}$ GeV and $10^{7}$ GeV. For what is discussed above, $L_{\gamma}\simeq 2/3L_{\nu}$ in $pp$ interaction and considering the neutrino flavor equipartition coming from neutrino oscillations $L_{\gamma}=2L_{\nu_{\mu}}$. To convert the luminosities computed in the different energy ranges, one obtains for an $E^{-2.2}$ spectrum $L_{\gamma}^{0.1-100}=1/10\ L_{\gamma}^{10^{4}-10^{7}}$. Therefore our result in Eq. 3.1 is not in contradiction with [63], demonstrating consistency of our simplified approach.

Figure 3 confirms that a neutrino spectrum of the form $dN/dE_{\nu}\propto E_{\nu}^{-2.12}$ with a source energy cutoff at 10 PeV is well in agreement with the TGM measurement and with the high-energy part of HESE. Smaller cutoff values are generally found to yield a better compatibility with the EGB, gradually falling short in explaining the highest energy neutrinos. An energy cutoff of 10 PeV for neutrinos at the source corresponds to an average cutoff of 200 PeV for protons. According to [9], SBGs may be capable of accelerating protons to such high energies, supporting the compatibility of HAGS with the multi-messenger observations. Moreover the possibility to accelerate protons up to the 100 PeV scale has been discussed in various works, invoking cosmic-ray accelerations due to hypernovae [65, 66], strong magnetic fields in Arp220-like galaxies [25], Galaxy mergers [67, 68] and AGN winds [12, 69, 70, 71]. Complementary to these expectations, recent measurements by the two Ultra-High Energy Cosmic Ray (UHECR) observatories see first indications for a directional correlation between the arrival directions of cosmic rays above 39 EeV and nearby SBG [72, 73]. Should HAGS be the sources of these UHECR, the local abundance of sources capable of reaching 200 PeV should be rather high.

3.1 Comparison with low energy events

We check the compatibility of our result with the other IceCube dataset. Concerning the low energy measurements, the expected number of HESE events are reported in Table 1. Using the HESE effective areas [1] and our baseline spectrum, we find that HAGS can account for 33 out of about 41 astrophysical signal events, obtained by subtracting the $\sim$ 41 expected background events from the 82 detected events in HESE [33]. While the event counts are well described by our model, the hard spectrum undershoots the second HESE bin, leaving space for other small contributions from other sources or background in this energy range. A possible additional contribution may come from the Milky Way’s galactic disc can be present among the events below 100 TeV from the Southern hemisphere. Galactic neutrinos are expected to give a contribution below $\sim 150$ TeV, reflecting the possible $\sim 3$ PeV knee of the primary proton spectrum [18, 74, 20]. The current estimates predict $\sim$1 neutrino/year in the HESE dataset (above 30 TeV), in line with the latest experimental limits [75]. Together with the galactic component, our model saturates the HESE signal event count to 94% at the best fit.

We also compare our result with the flux measured using the CAS4 sample. As mentioned at the beginning of Sec. 2, the single flavor flux suggested by the low energy data is $N=1.57^{+0.23}_{-0.22}$ (the parameter $N$ is defined at the beginning of Sec. 2) and the spectral index is $\alpha=2.48\pm 0.08$. On the other hand from our benchmark model we obtain a harder spectrum, having $N=1$ and $\alpha=2.12$. We compare the two energy fluxes, performing the ratio between the following integrals:

[TABLE]

obtaining $\mathcal{R}=0.40^{+0.07}_{-0.09}$ in the energy range between 10 TeV and 100 TeV, i.e. the sensitive energy range for the CAS4 dataset. This suggests that even below 100 TeV the contribution from HAGS can be relevant, accounting for about half of the energy flux. For the remaining 50% of the energy flux, various scenarios are plausible, such as a contamination from atmospheric prompt neutrinos [76], from neutrinos produced in the Galactic plane [18, 19] and/or from a second astrophysical component. However no conclusive interpretation can be provided with the present measurements.

4 Conclusion

We applied a multi-messenger approach to study the contribution of HAGS to IceCube’s diffuse neutrino flux combining the constraints from current $\gamma$-ray and neutrino observations. HAGS are (typically) infrared bright sources with a hard $\gamma$-ray spectrum from proton-proton interactions, for which the relation to the expected neutrino flux is well defined. The strongest constraint comes from the non-blazar contribution to the extragalactic $\gamma$-ray background.

We find that a hard power-law spectrum with an index $\alpha\leq 2.12$ and an energy cutoff at 10 PeV at the source to be compatible with the currently available estimations of the non-blazar contribution. In particular, due to the large uncertainties of these estimates it is impossible to derive tighter constraints solely from $\gamma$-ray observations and exclude HAGS, which include Starburst Galaxies and ULIRGs, as the dominant sources of diffuse neutrinos. Moreover, this conclusion remains valid even when an estimation of the misaligned AGN are further removed from the extragalactic $\gamma$-ray background. Such a spectrum is also found to be consistent with the spectrum of the NGC253 (a prototype of HAGS) above 50 GeV.

Concerning neutrinos, we find that our model can fully power the throughgoing muon energy flux detected by IceCube and the high energy part of the HESE flux, contributing to 80% of the total signal contained in the HESE dataset and to more than 50% of the signal contained in the CAS4 dataset. In conclusion, we do not find any contradiction between the hypothesis that the dominant contribution of the IceCube neutrino energy flux comes from HAGS, since they can provide the dominant contribution in the energy range between 10 TeV and 10 PeV, where the neutrino flux is measured nowadays. The remaining fraction can be attributed to neutrinos from the Milky Way’s galactic plane plus other possible source contributing below 100 TeV. Due to the soft spectral index and low maximal energies, the contribution from an entire population of “normal galaxies” (like ours) is at the few % level.

This result shows that relaxing the hypothesis that one single source class produces the entire astrophysical neutrino flux (as in [26]), it is still allowed for sources characterized by $pp$ interaction to provide the dominant contribution to the IceCube neutrinos above 100 TeV and to produce half of the observed energy flux between 10-100 TeV. Therefore this result shows that HAGS are still worthy to be investigated in the context of astrophysical neutrinos, in agreement with the findings of previous works [79, 80, 81, 82, 63, 83].

Compared to the brightest $\gamma$-ray sources (like blazars) HAGS are comparatively dim steady emitters, with their detection as neutrino sources requiring neutrino detectors with at least an order of magnitude more sensitivity compared to IceCube’s current sensitivity. In the $\gamma$-ray domain, the upcoming Cherenkov Telescope Array (CTA) will be able to discover more nearby HAGS, as well as extend the energy spectral energy range up to higher energies of those already detected. Being hadronic sources, HAGS are natural candidates for this neutrino emission. We find good compatibility of HAGS being the dominant neutrino source class.

Acknowledgments. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 646623).

Appendix A Supplementary material

The next sections are dedicated to more detailed calculations of i) the diffuse neutrino flux, ii) the expected number of events in the HESE and CAS4 data sets and iii) the contribution from Normal Galaxies.

A.1 Computation of diffuse neutrino flux

For hadronic, i.e. proton-proton, interactions the $\gamma$-ray and the all-flavor neutrino spectrum obey the following energy budget relation:

[TABLE]

where $\frac{dN}{dE}$ is the differential flux at the source. The diffuse neutrino flux expected from HAGS, we start from the $\gamma$-ray luminosity spectrum from a single source. We assume that it is a power law with a spectral index equal to $\alpha+2$, an energy cutoff $E_{\text{cut}}$ and we denote it as $\frac{d\mathcal{L}_{\gamma}}{dE}(E,E_{\text{cut}},\alpha)$. For each value of $\alpha$, the related neutrino spectrum is normalized in order to reproduce $L_{\nu}=3/2\ L_{\gamma}^{\text{NGC 253}}$ in the energy range between 0.1 GeV and 3 TeV, since we take the Starburst Galaxy NGC 253 as a benchmark object in our calculation. This object has been recently observed by both Fermi-LAT and HESS [38]. Under the previous assumptions the neutrino flux from a single source is given by:

[TABLE]

where $d\mathcal{L}/dE$ is in units of energy/sec. We denote the neutrino flux from a single source with $\phi(E_{\nu})$. where $z_{\text{max}}=6$, $D_{c}(z)=D_{H}\times d(z)$ is the comoving distance, $D_{H}$ is the Hubble distance and $d(z)=\int_{0}^{z}d\tilde{z}\ h(\tilde{z})^{-1}$. In the previous equation the terms $h(z)=\sqrt{\Omega_{\lambda}+\Omega_{m}(1+z)^{3}}$ with $\Omega_{\lambda}=0.73$ and $\Omega_{m}=0.27$.

In order to obtain the cumulative neutrino flux we parameterize the evolution of the sources in redshift as in [37]:

[TABLE]

where $a=3.4$, $b=-0.3$, $c=-3.5$, $B=5000$, $C=9$, $\eta=-10$. The parameter $\rho_{0}$ denotes the local density of sources at redshift ($z=0$). In our case $\rho_{0}(\alpha)$ is obtained by fitting data and depends on the spectral index (values are given in Tab.2). The diffuse neutrino flux is indicated using $\Phi(E_{\nu})$ and it is obtained from:

[TABLE]

where $dV/dz$ gives the relation between the comoving volume and the redshift as follows, $\frac{dV}{dz}=4\pi D_{H}^{3}\frac{d^{2}(z)}{h(z)}$. We have chosen $z_{\text{max}}=6$ in our calculations but we have checked that larger choices of $z_{\text{max}}$ do not produce any significant impact on the calculations.

A.2 Gamma-rays and neutrinos from Normal Galaxies

Normal galaxies, such as ours, are more abundant but less luminous than HAGS. Therefore it might be important to estimate their contribution to the diffuse $\gamma$-ray and neutrino fluxes, using the same procedure as for HAGS. As a typical luminosity of normal galaxy we use that of the Milky Way $\sim 10^{39}\rm\ erg/s$ in the 0.1 GeV - 3 TeV energy range [30]. The local density is estimated as $\rho_{0}=10^{-2.7}\ \rm Mpc^{-3}$ from HERSCHEL data (Figure 4 of [62]). The spectral index is varied in the range $\alpha\in[2.4-2.7]$, which includes the rather extreme KRA-$\gamma$ model [84] and that of galactic cosmic-rays. The resulting diffuse fluxes from normal galaxies are presented in Figure 4 for the range of spectral indices. Even for the hardest spectral indices the contribution to the EGB and the neutrino flux is at level of few %.