New Constraints on the Origin of Medium-Energy Neutrinos Observed by IceCube

Any conceivable scenario for the origin of high-energy cosmic neutrinos, observed by the IceCube neutrino detector at the South Pole, predicts the generation of accompanied high-energy $\gamma$ rays. Propagation of high-energy photons over cosmological distances initiates an electromagnetic cascade that degrades their energy to $\lesssim1$ TeV energies, contributing to the extragalactic $\gamma$-ray background (EGB) observed by the Fermi-LAT experiment. By taking into account various established components of the EGB, as well as the latest IceCube shower data in the 10-100 TeV range, we derive new multimessenger bounds on the properties of their sources. We find that conventional scenarios suffer from a tension at $\gtrsim3\sigma$ level for a neutrino flux extending down to $10$ TeV, which is interpreted as evidence for a new class of sources that prohibit $\gamma$ rays from escaping them or have redshift evolution significantly different from the cosmic star formation rate.

It is critically important to identify and understand the sources of the medium-energy neutrinos in the 10-100 TeV range. First, the large neutrino flux compared to the extragalactic γ-ray background (EGB) flux [11] is naturally explained by γ-ray hidden sources [23]. If es-tablished, the IceCube data will enable us to utilize neutrinos as a unique probe of particle acceleration in dense environments. Second, the large neutrino flux implies that the energy generation rate density of high-energy neutrinos is significant as one of the nonthermal energy budgets in the Universe [24]. Not many source candidates can satisfy the energy budget requirement, and possible candidate sources include the cores of active galactic nuclei (AGNs) (see a review [25]) and chokedjet supernovae (SNe) [26][27][28][29][30]. Revealing the sources is also important for us to understand the multi-messenger connection among neutrinos, γ rays and UHECRs. For example, the medium-energy neutrino flux cannot be explained by conventional γ-ray transparent sources such as galaxy clusters and starburst galaxies [31][32][33], so that a multi-component model may be required for the IceCube data from TeV to PeV energies [23,[34][35][36][37][38][39][40].
With the latest IceCube data in the 10-100 TeV range [19,20], and the EGB data from Fermi-Lat [11], this work provides the first quantitative constraints on the parameter space allowed by intent neutrino sources, from which GeV-TeV γ rays escape. We show that the conventional γ-ray transparent scenario suffers from the > ∼ 3σ tension with EGB data, which is regarded as evidence for hidden cosmic-ray accelerators or unknown Galactic sources.

MODELING OF ν AND γ-RAY SPECTRA
Production of high-energy neutrinos in astrophysical sources requires hadronic processes creating π ± and K ± which subsequently decay to neutrinos; e.g, π + → µ + ν µ → e + ν µνµ ν e . The pions can be produced via interactions of accelerated protons with ambient protons (pp scenarios) or photons (pγ scenarios). The resulting neutrino flux has different characteristics in each sce-nario: while in the pp scenario the neutrino spectrum extends to lower energies and increases with the decrease in energy, in the pγ scenario a large fraction of produced neutrinos have energies larger than the threshold energy ∼ 4 × 10 −2 m π m p /ε t ∼ 6 × 10 6 GeV (eV/ε t ) (where ε t is the energy of target photons), below which the neutrino spectrum drops rapidly. For this work, we parametrize the neutrino spectrum in the pγ scenario by introducing a break energy, ε br , where for ε ν < ε br the spectrum hardens due to pion-decay kinematics, and we take it ∝ ε −s l ν with s l = 0 [23]. Theoretical calculations of ε br require detailed knowledge on the source characteristics. We take the following energy spectrum: where ε ν Q εν = n s ε ν dLν dεν is the differential energy generation rate density of neutrinos with energy ε ν for neutrino luminosity L ν and the number density of the sources n s . The neutrino flux is conservatively set to zero for ε ν > 10 PeV since, so far, there is no observed neutrino flux at this energy range [41,42].
We emphasize that the neutrino spectrum in Eq. (1) is the minimal assumption about the neutrino production in the source(s) that can accommodate the diffuse neutrino flux observed by IceCube. Extending the energy range either to lower energies, as in pp scenario, or to higher energies, by increasing the assumed 10 PeV cutoff, increases the accompanying γ-ray flux.
The energy flux observed at the Earth from a source at redshift z is ε ν (dL ν /dε ν )| εν =(1+z)Eν /(4πd 2 L ), where d L is the luminosity distance, H(z) is the z-dependent Hubble parameter and E ν is the observed neutrino energy at Earth. Knowing the differential energy generation rate density of neutrinos at redshift z = 0, which is ε ν Q εν in Eq. (1), and the dimensionless redshift evolution of the sources, F(z), which we take to be the cosmic start formation rate discussed in the Appendix, the all-flavor diffuse flux of neutrinos at the Earth from the distribution of sources is given by where dV c /dz = 4π[c/H(z)]d 2 L /(1 + z) 2 and V c is the comoving volume.
The γ-ray flux accompanied by the neutrino flux is calculated using the following argument: from isospin symmetry, not only π ± but also π 0 have to be produced at the sources. The subsequent decay of π 0 to photons (π 0 → 2γ) generates a γ-ray spectrum given by where K ≈ 1 for pγ sources. The γ-ray flux should generally be larger than Eq. (3) because the charged pions can lose part of their energies before the decay, adiabatically or radiatively, and also ambient electrons and positrons can enhance the γ-ray production via cascades inside the sources. The calculation of the γ-ray flux at Earth is more complicated than the neutrino flux. Even for the most conservative setup, in which the sources are optically thin to γ rays, the Universe is opaque to γ rays with energy > ∼ 1 TeV that propagate distances z > ∼ 10 −2 , due to the absorption by pair production on the Cosmic Microwave Background (CMB) and Extragalactic Background Light (EBL) photons [43]. At > ∼ 1 PeV energies, this absorption is significant even at the Galactic scale [16,44]. However, the pairs produced in the pair production process inverse-Compton scatter off the CMB and EBL, creating new γ rays at slightly lower energies than the original γ rays. These successive processes initiate an electromagnetic cascade which ceases at the pair-production threshold ∼ m 2 e /ε t . For the CMB, this cutoff appears at ∼ 100 TeV, while for the EBL it is ∼ 100 GeV. Thus, although the Universe is opaque to high-energy photons, the initial high-energy γ-ray flux will be redistributed in the GeV-TeV range due to the electromagnetic cascade. Although the approximate spectrum can be calculated analytically [45,46], the exact energy dependence of the flux needs numerical calculations taking into account the z-dependence of the EBL and CMB. In this work, we use the public code γ-Cascade for this purpose [47], which agrees well with results of the previous literature [15].

MULTIMESSENGER ANALYSES
The diffuse and isotropic γ-ray flux arising from cascades induced by high-energy photons in the intergalactic space contributes to the EGB. A conservative limit on the neutrino sources has been derived in Ref. [15] by requiring that the resulting flux should not overshoot the isotropic diffuse γ-ray background (IGRB) part of the then measured EGB [48] (which was extending to ∼ 100 GeV) at any energy. In our analysis we consider the latest measured whole EGB data [11] which is extending to 820 GeV. However, there are other contributions to the EGB originating from populations of unresolved sources at low energies, i.e., < ∼ 1 TeV, which should be taken into account, such as the guaranteed contributions from jetted AGNs (including blazars and radio galaxies), star-forming galaxies and cosmogenic γ rays [49]. One approach is to subtract such point-source contributions from the EGB and use the remaining flux to set limits on any additional diffuse γ-ray contribution, including the cascaded flux that we are interested in. A more appropriate approach is to perform a χ 2 analysis by taking into account the various contributions.
In the following we describe two different analyses performed in this work: A. χ 2 analysis: The contribution of blazars to the EGB, including BL Lac objects and Flat Spectral Radio Quasars (FSRQs), has been calculated in Ref. [50]. We use their luminosity-dependent density evolution (LDDE) model for the luminosity function of the blazars. The emissions from starforming galaxies [51] and radio galaxies [52] have also been taken into account. Using these contributions and the EGB data we set a limit on any extra contribution to the EGB by defining the following χ 2 function: where F i,EGB , F i,a and F i,cas are, respectively, the observed EGB flux, the astrophysical contribution (blazars, star-forming galaxies and radio galaxies) and the cascaded flux contribution to the i-th energy bin. Also σ i is the uncertainty on EGB flux and the last term is the pull-term, taking into account the normalization uncertainty of the astrophysical contribution, given by σ A ≈ 35% [50].
B. Integrated flux above 50 GeV: It has been shown in Ref. [13] that 86 +16 −14 % of the total EGB above 50 GeV can be accounted for by the contribution from the sources in the 2FHL catalog, mainly consisting of blazars. The total EGB integrated flux above 50 GeV is J EGB >50 GeV = 2.4 × 10 −9 ph/cm 2 /s/sr. So, by requiring 820 GeV we can derive limits on the cascaded γ-ray flux Φ cas γ . In the above relation, q is the percentage of total EGB intensity (above 50 GeV) which can be explained by the blazars, with central value q = 86%.
Method B is an independent analysis which is not sensitive to the spectral shape of the cascaded flux as in method A. Also, the majority of sources in the 2FHL catalog are blazars (and among those, 74% are BL Lac objects). This means that the constraints derived from method B are very conservative and based on the contribution of a single source population to the EGB. Although the principal result of this Letter comes from the method A, we perform the analysis of method B as a sanity check.
The cascaded γ-ray flux from the distribution of sources responsible for the neutrino flux observed in Ice-Cube depends on s h and ε br via Eq. (1) (through Eq. (3)). Using the EGB data, we can derive constraints on the s h and ε br parameters, or equivalently, on s h and E br , as well as on the normalization of the corresponding neutrino flux. Here the E br is the redshifted observed energy break at Earth (see Appendix).

RESULTS
In analysis method A, constraints in the (s h , ε br ) plane, shown in Figure 1, are derived by defining ∆χ 2 = χ 2 − χ 2 min , where χ 2 min is the minimum value of χ 2 in Eq. (4) without the pull-term (free A) and a free normalization for the cascaded γ-ray flux. The s h range for each IceCube data set, see Appendix, is depicted and the gray shaded regions show the excluded ε br from IceCube data (by translating the E br to ε br for each data set). The color-shaded regions show the excluded ε br values from the EGB data at 90% C.L. limits derived from the condition ∆χ 2 < 4.61 (for 2 d.o.f.). For each color (corresponding to a different IceCube analysis), the upper and lower curves correspond, respectively, to the highest and lowest IceCube allowed normalizations, Φ astro defined in Eq. (7) in the Appendix, at 1σ (shown in Figure 2). Clearly, from Figure 1, the HESE and through-going ν µtrack data sets of IceCube are compatible with the EGB data, while the measured neutrino flux in the cascade data set leads to a diffuse γ-ray flux that is incompatible with EGB data. For each color, the upper and lower curves respectively correspond to the maximum and minimum allowed flux normalizations, Φ astro , at 1σ reported by IceCube.
To quantify the tension in Figure 1, using method A, we derive constraints in the (s h , Φ astro ) plane for fixed values of E br . The color-shaded regions in Figure 2 show the allowed regions in each IceCube data set in the (s h , Φ astro ) plane. The solid curves show the limits, at 2σ C.L., from method A of analyzing the EGB data for the depicted E br values, where the arrows point toward the allowed regions. We can see that having astrophysical neutrinos down to ∼ 10 TeV, as the 6-year cascade data set indicates [20], leads to a tension with the EGB data. As in Figure 1, the HESE and through-going ν µ -track analyses rely on the data above ∼ 60 TeV and ∼ 120 TeV, respectively, so they are compatible with the EGB data. Both the 4-year [53] and 6-year [20] cascade data sets are essential for the tension. From Figure 2, we can also conclude that extending the astrophysical neutrino flux to energies < ∼ 20 TeV results in tensions with the EGB data for all the three sets of IceCube data. The present shower data with E br ≈ 10 TeV is in tension with the EGB data at > ∼ 3σ C.L., whereas for E br ≈ 20 TeV, it drops to ≈ 2σ. The statistical significance of this tension increases in a more realistic setup. FIG. 2: Constraints in the (s h , Φ astro ) plane from method A of analyzing the EGB data. The solid black curves depict the allowed regions, for fixed E br , from EGB data. The green shaded regions show the allowed regions for the 4-year cascade events [53] which are similar to the 6-year cascade [20] allowed regions.
As an independent analysis, Figure 3 shows the results based on method B. The solid (dashed) curves correspond to the highest (lowest) allowed normalization of astrophysical neutrinos at 1σ level. The label on each curve shows the q value in Eq. (5). Consistent with method A, Figure 3 shows the tension between the IceCube cascade data set and the EGB data for q > ∼ 80%. Obviously, method B is less constraining since the analysis is based on just the integrated flux of EGB above 50 GeV and is independent of the spectral shape of the cascaded flux, which in fact is important at ∼ 100 GeV.
The redshift evolution slightly affects the tension quantitatively but not qualitatively and the conclusions remain the same for redshift evolution of the most of the source classes including galaxy clusters, star-forming galaxies, and AGNs [15].  5)).

SUMMARY AND DISCUSSIONS
The neutrino flux observed in IceCube should be accompanied by the γ-ray flux, which provides a powerful diagnostic in the search for their possible sources. Assuming a minimal model for high-energy cosmic neutrinos, for the first time, we showed that the new IceCube data extended down to < ∼ 10 TeV leads to > ∼ 3σ tension with the EGB data from Fermi-LAT. We stress that the derived limits and reported tension are based on very conservative assumptions. The tension is ≈ 3σ for a break energy of E br ≈ 10 TeV, and larger for more realistic setups. First, the neutrino spectrum is modified by the cooling of mesons and muons, which yields a larger ratio of γ rays to neutrinos. Second, additional γ rays must be produced by the Bethe-Heitler process; for example, these Bethe-Heitler-induced γ rays are dominant in the AGN core scenario [54]. Third, γ rays should also be produced by leptonic processes which do not produce any neutrinos. GeV-TeV γ rays of blazars are conventionally explained by the leptonic components.
The reported tension suggests an additional population of the sources, which are different from conventional cosmic-ray reservoirs. Hidden (γ-ray opaque) cosmicray accelerators are among the promising sources of the medium-energy IceCube neutrinos. Candidate classes include choked GRB jets [30,55,56], AGN cores [54,[57][58][59][60], and MeV blazars [23]. Alternatively, high-redshift source population that do not exist in the local universe can alleviate the tension. For example, with the redshift evolution of POP-III stars, the EBL cutoff can be down to 10 GeV energies [61]. Finally, in principle, Galactic sources that lead to quasi-isotropic emission, such as the Galactic halo [16,62,63], may give a significant contribution. Although the 10-100 TeV neutrinos come from both hemispheres and there is a tension with some of the upper limits from air-shower experiments [16,23], further multimessenger studies are necessary [64]. In fact, the claimed upturn in the IGRB [65] can support such Galactic halo scenarios. Our results also impact nonastrophysical scenarios that explain the shower data with physics beyond the Standard Model (BSM) (see reviews [66,67]). For example, decaying dark matter has been invoked as an interpretation of the IceCube data [68,69] (see Refs. [70,71] for recent analyses). Final states involving quarks, charged leptons and gauge bosons are accompanied by a comparable γ-ray flux [44,72], which gives strong constraints especially for models explaining the medium-energy neutrino data [73][74][75][76]. Other BSM explanations, such as neutrino decay [77], increase the ratio of γ rays to neutrinos, which further strengthens the results of this work [78].
Further observations of the medium-energy range (by more efficient rejection of background events) to lower energies is of crucial importance. IceCube-Gen2 will give us more statistics, but the threshold energy should not be far from ∼ 10 TeV. KM3NeT [79] will be able to give us information on the northern sky, which is complementary, and the detection of showers with a better angular resolution will be particularly useful. In addition, stacking searches with source catalogues at different wavelengths are strongly encouraged. Intriguingly, a hidden cosmic-ray accelerator with a steep neutrino spectrum is independently indicated from the recent ∼ 3σ observation of NGC 1068 [8,80]. Searching for lowerenergy γ-ray counterparts in the MeV energy range will also be important.
We thank Ali Kheirandish and Pasquale Serpico for useful discussions and comments. This work has been supported by the Alfred P. Sloan Foundation and NSF

Redshift evolution
For F(z) in Eq. (2) and the corresponding computations of the γ-ray flux, we use the cosmic star formation rate (SFR) [81,82] given by where a = 3.4, b = −0.3 and c = −3.5. The constants B 5000 and C 9 correspond to breaks at z 1 and z 4, respectively, and η = −10 smooths the transition between the breaks. The normalization in Eq. (1) is fixed by the observed IceCube neutrino flux. Notice that the abrupt break in the injection spectrum at 10 PeV is smoothed at the Earth due to cosmological redshift. The same effect causes the position of the break in the energy spectrum at the Earth, E br , to be shifted with respect to the energy break at the source, ε br . In our case, this shift depends only on the spectral index s h and on our choice of adopting the SFR evolution. Figure 4 shows E br in terms of ε br for various s h values, where the diagonal gray line depicts E br = ε br , that is, no redshift. We can see that an increase in s h results in a decrease in the ratio E br /ε br , and for s h 3, it reaches approximately 50%. The astrophysical neutrino flux has been measured by IceCube in several channels. The channels can be char-acterized by the event topology, either cascade or ν µtrack events, and the location of neutrino-nucleus vertex, which can be either inside or outside the fiducial volume of IceCube, leading to starting or through-going ν µ -track events, respectively. The measured differential flux from the data in each channel can be parametrized by (in units of [GeV −1 cm −2 s −1 sr −1 ]) where Φ astro and s ob are the (observed) normalization and energy index of the flux. Since the background events for each channel are different, the minimum observed energy, or threshold energy E thr , which depends on the efficiency of background rejection at low energies, varies among the data sets.
We consider the following three data sets: i) 7.5-years of High Energy Starting Events (HESE) over the fullsky [83], consisting of both cascade and ν µ -track events with the interaction vertex inside the fiducial volume of IceCube and with the threshold energy E thr = 60 TeV. The reported all-flavor normalization and energy index are (1σ error) Φ astro = 6.45 +1. 46 −0.46 and s ob = 2.89 +0.2 −0. 19 . ii ) 6-years cascade events [20] over the entire sky with E thr = 16 TeV, one-flavor normalization Φ astro = 1.66 +0. 25 −0.27 and energy index s ob = 2.53 ± 0.07. The precedent 4-year cascade data set [53] has almost the same normalization and index. iii ) 9.5-years of through-going ν µ -track events over the northern hemisphere [84] with E thr = 119 TeV, one-flavor normalization Φ astro = 1.44 +0. 25 −0.24 and energy index s ob = 2.28 +0.08 −0.09 . All the reported normalization and energy index values in the three data sets come from single power-law fits to data. For all data sets, a broken power-law fit also has been performed showing no preference over the single power-law fit.
The γ-ray data set consists of the Extragalactic γ-ray Background (EGB) measured by the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) [11]. The EGB is the sum of contributions from all the extragalactic γ-ray sources, including individual sources (faint and unresolved sources) and diffuse ones such as the Galactic foreground and (possible) contributions from electromagnetic cascades and dark matter annihilation/decay. The latest EGB data set covers the energy range 100 MeV to 820 GeV.