Multi-messenger gamma-ray counterpart of the IceCube neutrino signal
A. Neronov M. Kachelriess D. V. Semikoz



 





^1Astronomy Department, University of Geneva,
Ch. d'Ecogia 16, Versoix, 1290, Switzerland
^2Institutt for fysikk, NTNU, Trondheim, Norway
^3APC, Universite Paris Diderot, CNRS/IN2P3, CEA/IRFU,
Observatoire de Paris, Sorbonne Paris Cite, 119 75205 Paris, France
^4National Research Nuclear University MEPHI (Moscow Engineering Physics Institute), Kashirskoe highway 31, 115409 Moscow, Russia




A signal of high-energy extraterrestrial neutrinos from unknown source(s)
was recently discovered by the IceCube experiment. Neutrinos are always
produced together with γ-rays, but the γ-ray flux from
extragalactic sources is suppressed due to attenuation in the intergalactic
medium. We report the discovery of a γ-ray excess at high Galactic
latitudes starting at energies 300 GeV in the data of the Fermi telescope.
We show that the multi-TeV γ-ray diffuse emission has spectral
characteristics at both low and  high Galactic latitudes compatible with
those of the  IceCube high neutrino signal in the same sky regions.
This suggests that these γ-rays are the counterpart of the IceCube
neutrino signal, implying that a sizable part of the IceCube neutrino flux
originates from the Milky Way. 
We argue that the diffuse neutrino and γ-ray signal at high Galactic
latitudes originates either from  previously unknown nearby cosmic ray
"PeVatron" source(s), an extended Galactic CR halo or from decays of heavy 
dark matter particles. 



Multi-messenger gamma-ray counterpart of the IceCube neutrino signal
    D. V. Semikoz^3,4
    December 30, 2023
====================================================================







§ INTRODUCTION


The discovery of an extraterrestrial neutrino signal in the TeV–PeV energy range by the IceCube collaboration <cit.> has recently opened the era of multi-messenger astronomy. High-energy neutrinos are produced by cosmic rays (CR) interacting at their acceleration sites or during propagation through interstellar and intergalactic space <cit.>. Alternatively, neutrinos may be produced in decays of metastable heavy dark matter (DM) particles <cit.>. 
The source(s) of this  neutrino signal have remained unidentified so far because of the limited statistics of the IceCube data. Moreover, the High-Energy Starting Events (HESE) which provide the most significant contribution to the neutrino signal have a poor angular resolution. At the same time, the production of high-energy neutrinos is accompanied by γ-rays.  This implies that the neutrino sources could be identified using a "multi-messenger" approach by combining neutrino and γ-ray data <cit.>. 

The TeV–PeV γ-ray flux from distant sources is suppressed by  electron-positron pair production in interactions with low-energy photons of the extragalactic background light and the cosmic microwave background <cit.>.   Therefore, the presence or absence of a γ-ray counterpart can be used to clarify the origin of the neutrino signal: If the signal originates from extragalactic sources at cosmological distances, no γ-ray counterpart is expected in the multi-TeV to PeV band. In contrast, a Galactic origin  implies the presence of a comparable multi-TeV γ-ray flux. 

The search for the γ-ray counterpart of the  neutrino signal is challenging with both ground and space-based γ-ray telescopes. Ground-based telescopes like HESS [<https://www.mpi-hd.mpg.de/hfm/HESS/>], MAGIC [<https://magic.mpp.mpg.de>] and VERITAS [<https://veritas.sao.arizona.edu>] or air shower arrays like HAWC [<https://www.hawc-observatory.org>] and ARGO-YBJ [<http://argo.na.infn.it>] suffer from a high background of events produced by charged CRs <cit.>. The arrival directions of the CR background events are distributed over large angular scales, similar to the expected γ-ray counterpart of the neutrino signal. Space-based telescopes like the Fermi Large Area Telescope (LAT) <cit.> achieve a much better suppression of the charged CR background, but they have small collection areas which severely limit the signal statistics.   

In this Letter, we report a study of the TeV diffuse gamma-ray sky based on the data of Fermi/LAT. The small effective area of Fermi/LAT is compensated by the very long exposure time of nine years of the Fermi/LAT data we use. We show that the   flux and spectrum at low and high  Galactic latitudes are compatible with the flux of the measured neutrino signal, in the energy range where the two signals overlap. We suggest that the γ-ray in the multi-TeV band  is the counterpart of the IceCube neutrino signal. 





§ CROSS-CALIBRATION OF THE LAT DATA IN THE MULTI-TEV BAND


Our analysis uses  events from the ULRACLEANVETO class collected by Fermi/LAT  during the period
between October 28, 2008 and December 15, 2017. We calculate the spectra of large sky regions using the  "aperture photometry" approach [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/aperture_photometry.html>]. 

The energy resolution and the calibration of the telescope effective area degrade in the TeV band  <cit.>. Therefore we perform an additional cross-calibration of the Fermi/LAT flux measurements with those of ground-based γ-ray telescopes via a comparison of spectral measurements of the stacked spectra of selected calibration sources, see the Appendix for details.  We find that  a cross-calibration factor κ=1-clog(E/100), with c=0.25± 0.12  has to be applied to the LAT flux measurements above  300 GeV to achieve better consistency with the ground-based telescope measurements. We apply this factor in our analysis, following a practice common in X-ray data analysis [<https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/caldb_xcal.html>]. The uncertainty of the parameter c  is taken into account as an additional systematic error. We have verified that the cross-calibration factor also assures the consistency of the Fermi/LAT measurements of diffuse TeV emission from large regions of the sky with the measurements by the ground-based air shower arrays ARGO-YBJ <cit.> and MILAGRO <cit.>, see the Appendix for details.




§ DIFFUSE TEV  SIGNAL


Figure <ref> compares Fermi/LAT γ-ray spectra of the full sky (upper panel), of the Galactic plane |b|<10^∘ (middle panel) and at  Galactic latitudes |b|>10^∘ (lower panel) with the neutrino spectra of the same sky regions  <cit.>. In the spectra of the all-sky and the |b|>10^∘ region we remove residual CR background, while the Galactic plane spectrum is calculated by subtracting high Galactic latitude background and residual cosmic ray contributions (see the Appendix for details). The γ-ray and neutrino all-sky flux and spectral slope measurements agree   in the overlapping multi-TeV band, confirming a previous analysis based on the high-energy extrapolation of the Fermi/LAT spectrum <cit.>. Figure 1 also shows the model of diffuse γ-ray emission from pion decays derived from an all-sky analysis of the LAT data <cit.>. It is this component which is expected to have the neutrino counterpart, since pion decays produce simultaneously γ-rays (π^0 decays) and neutrinos (π^± decays). The power-law extrapolation of the pion-decay model <cit.> into the multi-TeV band agrees with the neutrino spectrum measured by IceCube. 






This agreement suggests the interpretation of the TeV γ-ray signal as the multi-messenger counterpart of the neutrino signal. However, the γ-ray flux below TeV is dominated by the emission from the Galactic plane, while only a moderate fraction of the neutrino flux in the 100 TeV range comes from the Galactic plane <cit.>. A consistent interpretation of the multi-TeV γ-ray flux should provide an explanation for this fact. If the multi-TeV γ-ray flux is the counterpart of the neutrino signal, the high Galactic latitude flux should have harder spectrum than the flux from the Galactic plane so that its relative contribution to the all-sky flux could grow with increasing energy. A hint of such a behavior can be noticed in the bottom panel of Fig. <ref> where a hardening of the flux is noticeable in the last two energy bins.

This hardening appears more pronounced in the analysis of the spectrum of the part of the sky at higher Galactic latitude, |b|>20^∘, shown in Fig. <ref>. In this figure we have removed contributions from resolved point sources, extragalactic isotropic diffuse γ-ray background (IGRB) and residual CR backgrounds thus leaving only the Galactic diffuse emission (see the Appendix).  The hardening of the  spectrum of diffuse emission  at  high Galactic latitudes starts at 300 GeV and it can not be explained by instrumental effects (see the Appendix for details).  Below 300 GeV the spectrum is well fit by a smoothly broken power-law with the slope Γ=2.906± 0.015 in the 30–300 GeV range.  The spectrum in the 0.3–3 TeV range has the slope Γ=2.09± 0.09.

The most significant excess above the extrapolation of the  power-law valid below 300 GeV is in the energy bin 1–1.7 TeV. The model prediction of the number of photon counts in this bin is 16.4.   The observed number of counts is 39. The chance coincidence probability of such an excess is 1.5× 10^-6. In the energy bin 1.7–3.16 TeV the expected number of counts is 3.8, while the observed one is 10. The chance coincidence probability of such an excess is 5.8× 10^-3. In the energy bin 0.3-1 TeV, the model prediction is <66.5 counts while the observed signal is 100 counts. The chance coincidence probability of the excess in this bin is  8× 10^-5.  The energy-binning independent combined chance probability of the excess  above 300 GeV is less than 8× 10^-10.  






§ INTERPRETATION


The multi-TeV band  flux at high and low Galactic latitude shown in Fig. <ref> originates from the Milky Way. The low Galactic latitude flux is certainly dominated by the emission from decays of pions produced in interactions of Galactic CRs with interstellar matter in the Galactic disk. Since the  and neutrino fluxes from pion decays are comparable, the  flux measurement in the multi-TeV range can be used as an estimate for the minimal possible neutrino flux from the Galactic plane. One can see from the middle panel of Fig. <ref> that this lower bound on the neutrino flux is consistent with the upper limit derived by the IceCube and ANTARES telescopes <cit.>. Combining the lower and upper limits one finds that the neutrino flux from the Galactic plane has to be  just at the level of the multi-TeV  flux from this part of the sky. 

More puzzling are the spectral characteristics of the multi-messenger signal at high Galactic latitudes.  The conventional high Galactic latitude diffuse emission components have soft spectra in the TeV range <cit.> and can not explain the observed spectral hardening above 300 GeV.  The same is true for the IGRB, which is dominated by the cumulative flux of blazars <cit.>, a special class of active galactic nuclei  which do not provide the dominant contribution to the neutrino signal <cit.>.  Thus, the observed hardening of the  spectrum has to be interpreted as due to the presence of a new Galactic γ-ray flux component above 300 GeV. It is this component which is the counterpart of the neutrino signal with comparable flux in the multi-TeV range. 

Only few source types could produce multi-TeV multi-messenger  emission on large angular scales at high Galactic latitude with a hard spectrum. One possibility is interactions of CRs forming a previously unknown component of the Galactic CR population.  If this new component would reside everywhere in the Galactic disk, an equivalent spectral hardening would be observed in the spectrum of  the Galactic plane—which is not the case. Instead, the hard spectrum CRs could either reside in our local Galactic environment, or be a part of a very large halo.








The local source of CRs with a hard spectrum reaching PeV energies (a "PeVatron") should be  a recent  and nearby  source, like e.g. the Vela supernova <cit.>. It should have injected CRs less than  10^5 year ago at a distance d not larger than several hundred parsecs. These two conditions are required for the presence of PeV CRs which produce 10–100 TeV neutrinos and the large angular extent Ω of the multi-messenger emission <cit.>. Cosmic rays with total energy  U_CR∼ 10^50 erg injected by the PeVatron and loosing their energy on the time scale t_ pp≃ 1.5× 10^8(n_ ISM/0.5^-3) yr in interactions with the interstellar medium of the density n_ ISM∼ 0.5 cm^-3 produce the γ-ray and neutrino flux F=U_CR/(4π d^2 Ω t_pp) 
with magnitude

    F∼
    2× 10^-7(Ω/2π sr)^-1n_ ISM/0.5/cm^3( d/0.3 kpc)^-2 GeV/ cm^2  s  sr .

This flux estimate matches the observed signal level, cf. with Fig. 2. Otherwise, the high Galactic latitude emission could be from a very large (hundred kiloparsec) CR "storage" around the Milky Way disk <cit.>. 

The local PeVatron model predicts strong  variability of the multi-messenger signal  across the sky. This variability is determined by the peculiarities of the energy-dependent spread of the CRs and  of the matter distribution in the local Galaxy. Low energy CRs which had no time to escape from the source region would not contribute to the large angular scale emission. This leads to a low-energy hardening of the spectrum, as shown in the top panel of Fig. 2 <cit.>.  In contrast, the signal is not expected to experience neither strong fluctuations nor a low-energy hardening in the large scale halo model <cit.>. 

An alternative possibility shown in the bottom panel of Fig. 2 is that decays of metastable DM particles X with mass m_X≃ 5 PeV generate photons and neutrinos <cit.>. The spectral shape of the decay mode X→q̅q→ hadrons is determined by Quantum Chromodynamics (QCD). Since at the end of the  QCD cascade quarks combine more easily to mesons than to baryons, mainly neutrinos and photons from pion decays are produced.  The γ-ray and neutrino flux measurements constrain the X particle lifetime  to be τ_X∼ 2× 10^27(Ω_X/Ω_DM)^-1 s, where Ω_X/Ω_DM is the fraction of the DM in the form of X particles <cit.>. Since the mass m_X  is above the unitarity limit <cit.>, the X particles were never in thermal equilibrium. They should have been  produced by gravitational interactions or other non-thermal processes  and may serve as a tool to study the earliest phases of the Universe. 

 The DM decay neutrino signal has a sizable extragalactic contribution, while its γ-ray  component in the TeV-PeV range has only the Galactic part. This leads to a systematically lower normalisation of the multi-TeV γ-ray component. The same is true for the large scale CR halo, which should be present around all galaxies, so that the neutrino flux is expected to have a significant extragalactic contribution. To the contrary, the neutrino and γ-ray components in the local PeVatron  model both originate from the Milky Way. The absence of the extragalactic component leads to similar γ-ray and neutrino fluxes   (see top panel of  Fig. 2). 
 
 The DM halo of the Galaxy is denser in the direction of the inner Galaxy. This means that in the  DM model, the flux from the inner Galaxy should be stronger than that from the outer Galaxy. However, the signal from the Galactic plane shown in Fig. <ref> contains both the direction toward the Galactic center and the anticenter, from which the strongest and the weakest DM decay signal should be observed. We have verified that the expected excess of the DM decay signal from the Galactic Plane is consistent with the IceCube upper bounds on the Galactic plane flux. The fraction of the DM decay signal from the region |b|<10^∘ is 0.22. Combining the information form Fig. <ref> and Fig. <ref>,  one can see that the neutrino flux from the high Galactic latitude region which is supposed to account for the full neutrino signal at high Galactic latitude at 100 TeV is at the level 6× 10^-7 GeV/cm^2s at this energy (cf. with the bottom panel of Fig. <ref>). Re-scaling it by a factor 0.22/0.78≃ 0.3, one could check that the expected DM decay flux from the direction of the Galactic plane is at the level of 2× 10^-7 GeV/cm^2s, i.e. marginally consistent with the IceCube upper limit on the neutrino flux from the Galactic plane (the IceCube upper limit is exactly at the level of the flux estimate, which means that the signal of DM origin should soon reveal an excess toward the inner Galaxy). There is, however, one important reservation which should be added.  The IceCube upper limit on the Galactic emission is derived assuming certain spatial template for the signal distribution. This template does not correspond to the spatial template of the DM signal. Thus, the IceCube limit on the Galactic emission is not directly comparable to the DM model prediction. 
 
 For the local PeVatron model, there is no fixed spatial template because the source morphology is not known. No  excess toward the Galactic Plane is generically expected. In this respect, the IceCube limit on the Galactic emission component does not provide constraints on the local PeVatron model.


The distinction between possible models of the multi-messenger  signal based on  spectral or spatial characteristics will be possible with next generation instruments like the IceCube-Generation II <cit.>, KM3NeT <cit.> neutrino telescopes and the space-based γ-ray telescope  HERD[<http://herd.ihep.ac.cn>] which will accumulate higher signal statistics. The detection of the γ-ray part of the signal by  ground-based telescopes like CTA[<https://www.cta-observatory.org>], LHAASO <cit.> and CARPET <cit.> will be possible provided that a sufficiently high (by a factor ∼ 10^5) rejection level of the CR background is achieved.    




§ CONCLUSIONS


We have demonstrated that the properties of the large scale diffuse Galactic  flux in multi-TeV band are compatible with the flux and spectrum of the neutrino signal in 1-100 TeV range, so that the two signals may be considered as different components of one and the same "multi-messenger" signal in the multi-TeV sky.  The  γ-ray flux at high Galactic latitude exhibits a pronounced hardening above 300 GeV, while no hardening is observed in the low Galactic latitude flux. This effect explains the lower contribution from the Galactic plane to the neutrino signal at higher energies, as observed by IceCube.  We have suggested three possible models which could explain the observed hard spectrum high Galactic latitude multi-messenger emission above 300 GeV: (i) interactions of CRs injected by a recent nearby cosmic PeVatron, 
(ii) CR interactions in a large halo around the Milky Way, or 
(iii) decays of  DM particles. 


We would like to thank T. Porter and V. Savchenko for useful discussions of the data analysis. 











§ FERMI LAT DATA SELECTION AND DATA ANALYSIS
 

Our analysis uses data of the Fermi/LAT telescope collected during the period between October 28, 2008 and December 15, 2017 [We discard the data between August 8 and October 28, 2008, as recommended by the Fermi Science Support Center for the analysis of the data above 30 GeV, see <https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html>]. We use  the ULTRACLEANVETO class <cit.> events which have the lowest residual CR contamination. 

The data are processed using the version v10r0p5 of the Fermi Science Tools[<https://fermi.gsfc.nasa.gov/ssc/data/analysis/>], via
a gtselect – gtmktime – gtbin – gtexposure chain to produce the spectra of different parts of the sky using the "aperture photometry" approach [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/aperture_photometry.html>]. The exposures for large regions of the sky are calculated averaging the exposures estimated on a grid of points with 10 degree spacing. 

We have verified that the spectra of isolated point sources extracted using this method are consistent with those extracted using an unbinned likelihood analysis [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html>]. As an example, Fig. <ref> shows a comparison of Crab spectra extracted using the two methods.  The two spectra are also extracted using two different photon selections: SOURCE for the likelihood and ULTRACLEANVETO for the aperture photometry. One can see that the  errorbars of the ULTRACLEANVETO measurements (using aperture photometry) are somewhat larger because of the lower signal statistics. 









 §.§ LAT analysis in the multi-TeV band



The spectral measurements based on the likelihood analysis do not extend into the TeV band because the public version of the gtlike tool has as analysis limit ∼ 850 GeV. The Fermi Science Support Centre provides photon data in the energy range up to 10 TeV and information on the instrument characteristics (e.g. energy resolution, effective area) derived from Monte-Carlo simulations up to 3.16 TeV [<http://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm>]. These data could still be used for the aperture photometry analysis. 

The energy resolution of the telescope degrades in this energy band because of the increasing leakage of the signals produced by particle showers in the calorimeter and because of the saturation of the calorimeter crystals <cit.>. Still, reliable estimates of the energy are achieved up to at least 3 TeV, as described in Refs.  <cit.>. The energy resolution decreases from 10% at 1 TeV to 25% at 3 TeV.   In our analysis, we bin events in wide energy bins (4 bins per energy decade) which are much wider than the energy resolution over the entire analysis energy range. The most recent analysis results extend into the multi-TeV energy <cit.> thus validating the energy calibration in the TeV band based on the real data, via a direct comparison of the spectra derived from Fermi/LAT with the measurements by the ground-based γ-ray telescopes.  

The rapid degradation of the energy resolution might also result in the effect of a "pile-up" of higher energy events which were mis-reconstructed and attributed an energy close to a characteristic energy at which the energy resolution starts to worsen. In order to explore if pile-up effects might affect our spectral measurements, we have repeated the analysis using events from the EDISP2 and EDISP3 event sub-selections. This sub-selection is characterised by increasingly better energy reconstruction quality [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_overview.html>]. Higher energy events with poor energy reconstruction which might produce a pile-up effect are naturally excluded in the high quality energy reconstruction event sub-samples. The spectra extracted using the EDISP2 and EDISP3 event sub-selections are consistent with those based on the full ULTRCLEANVETO sample. 
 
Apart from the energy resolution, the systematic uncertainty of the effective area also keeps growing from 5% at 100 GeV to ≥ 15% above 1 TeV [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html>]. Extrapolating into the multi-TeV energy range, one finds that the uncertainty exceeds 25%, i.e. the overall effective area is uncertain by a factor of 2.  Such a large uncertainty motivates a cross-calibration of the Fermi/LAT flux measurements with those of ground-based γ-ray telescopes, as  described below. 

Most of the observations of astronomical sources in the multi-TeV band are done using ground-based γ-ray telescopes, including Imaging Atmospheric Cherenkov Telescope (IACT) systems (HESS, MAGIC and VERITAS) and air shower arrays (MILAGRO, HAWC, ARGO-YBJ, Tibet-ASγ). The validation of the instrument response function of Fermi/LAT is possible via a cross-calibration of the LAT observations of selected sources with ground-based γ-ray telescope observations.   Taking into account the limited statistics of the LAT data in the multi-TeV band, we perform a comparison of spectral measurements of selected calibration sources for the stacked source signal rather than on a source-by-source basis. The sources selected for the cross-calibration purposes should have the following basic characteristics which make them suitable for the calibration analysis.

First, we require that the sources are steady in the TeV band. This excludes active galactic nuclei which are known to be strongly variable. Among the Galactic sources pulsar wind nebulae and supernova remnants are suitable.  All the pulsar wind nebulae and supernova remnants are extended sources. A further selection criterion for the calibration sources is the requirement that the source should have a well constrained spatial morphology. Uncertainties in the spatial structure of the source lead to uncertainties in the re-calculation of the flux measurements for different telescopes because of differing telescope point spread functions. The spatial morphology constraint leaves only a handful of isolated TeV γ-ray sources for the cross-calibration analysis. These sources are listed in Table <ref>.   





The stacked spectrum of selected sources is shown in Fig. <ref>. For each source, the source signal was extracted from a circle of the radius listed in the 4th column of Table <ref> (the circle radius is adjusted to cover the source extent and to include the wings of the LAT point spread function). The background is estimated from circles with radius either equal to the signal circle radius (for extended sources) or to 1^∘ (for Crab) and shifted from the source positions along constant Galactic latitude.  The reference spectral model is obtained by averaging the model spectra of individual sources derived in the References listed in Table <ref>. One can see that the Fermi/LAT measurement in the TeV band generally agrees with the spectral measurements done using ground-based γ-ray telescopes. The calibration of Fermi/LAT using ground-based measurements could be explicitly forced via a re-normalisation of the signal in the multi-TeV energy range on the model stacked source spectrum, as shown in Fig. <ref>. The comparison of the LAT and ground-based telescope data in multi-TeV band shows that the agreement of the flux measurements for the stacked source spectrum is reached if the LAT flux is renormalised by ≃ 20% (comparable to the systematic error [<https://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html>]).  We apply a cross-calibration factor κ=1-clog(E/100), with c=0.25± 0.12  at the energy E>300 GeV to achieve better consistence with the ground-based telescope measurements, following a practice common in X-ray data analysis, see e.g. <https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/caldb_xcal.html>. The uncertainty of the parameter of the cross-calibration factor  c is taken into account as an additional systematic error. 






An additional cross-check of the Fermi/LAT calibration in the multi-TeV band can be extracted from a comparison of the spectra for the part of the Galactic plane observed  by  ARGO-YBJ <cit.> and MILAGRO <cit.> air shower arrays. Figure <ref> shows the combined Fermi/LAT + ARGO-YBJ + MILAGRO spectrum of the 40^∘<l<100^∘ region of the Galactic plane. The high Galactic latitude signal discussed in the main text is used as the background estimate for this region. One can see that the Fermi/LAT measurements agree with both ARGO-YBJ and MILAGRO data for this sky region. To calculate the Fermi/LAT spectrum of an extended region of the sky, we average the exposure  calculated on a grid of points with 10^∘ spacing using the gtexposure tool.  


 






 §.§ Residual CR background estimate and systematic effects



The ULTRACLEANVETO event sample contains residual charged CR background events which are arriving from random directions on the sky and could mimic a nearly isotropic γ-ray signal. A study of the residual CR background contamination in the related sub-selection of ULTRACLEAN events with additional veto imposed to reduce CR background  was reported in Ref. <cit.> for the PASS7 event selection. This study shows that the level of residual CR background at 850 GeV is at the level of 10% of the γ-ray flux from the high Galactic latitude region in this energy range, as shown in Fig. <ref>. This study has derived the residual CR background count rate which follows a power-law dependence on equivalent γ-ray energy above 50 GeV. The recalculation of this power law into an equivalent diffuse emission flux is performed by dividing the residual CR count rate by the energy-dependent effective area. This results in the residual CR background flux shown by dashed line in Fig. <ref>. Applying the same method we extend the residual CR background flux model to the energy range above 1 TeV, assuming that the power law for the CR count rate extends with the same slope into the multi-TeV energy range. One can see that the residual CR flux could not provide the dominant contribution to the high Galactic latitude emission in the multi-TeV range.





The study of Ref. <cit.> was based on the IGRB class of the PASS7 event selection. This class is  a sub-class of the ULTRACLEAN event class with additional veto conditions applied to reduce the residual CR background level. Our analysis is based on the ULTRACLEANVETO class of the PASS8 event selection. The IGRB class of PASS7 is not publicly available and a direct comparison of the ULTRACLEANVETO/PASS8  and IGRB/PASS7 event samples is not possible. However, we have verified that the total (γ-ray + residual cosmic ray background) fluxes of the sky regions |b|>20^∘ in the two event classes are compatible within the systematic uncertainty. This is shown in Fig. <ref>.


 


There is a systematic shift between the PASS7 and PASS8 measurements. It is possible, in principle, that this shift is due to higher level of residual CR background in the PASS8 event selection rather than to a different modelling of the instrument response functions. Adopting this hypothesis, one can estimate "conservatively" the maximal possible level of residual cosmic ray background in the PASS8 ULTRACLEANVETO event sample by subtracting the  flux of the |b|>20^∘ part of the sky derived in Ref. <cit.> from the total count spectrum calculated for the PASS8 ULTRACLEANVETO event selection. This estimate of the maximal possible residual cosmic ray background is shown in Fig. <ref>. One can see that subtracting the maximal possible residual CR background does not alter the properties of the hard excess above 300 GeV.









 §.§ Point source flux subtraction








 Figure <ref> of the main text shows the spectrum of high Galactic latitude γ-ray emission after subtraction of not only the residual CR and isotropic  diffuse γ-ray backgrounds, but also of the flux from isolated catalog sources <cit.>.

For the determination  of the point-source flux we have used the method described in Ref. <cit.>. First we stack the angular distributions of photons around the brightest point sources to obtain a measurement of the point-spread function in each energy interval. Next, we take all sources from the Fermi catalog <cit.> and perform a stacking analysis of the source and background signal around them. For the calculation of the background we exclude photons in circles with radius equal to the 95% containment radius of the point-spread function and shuffle the remaining photons in Galactic longitude l according to the Fermi exposure. 

The  point-source flux for |b|>20^∘ in the 10 GeV band constitutes the same fraction of the total flux as found in Ref. <cit.>. However, we find a smaller point-source fraction in the energy range above 100 GeV, see Fig. <ref>. The signal of point sources in this energy range is dominated by the contribution of BL Lacs and unidentified blazars <cit.>, most of which are also BL Lacs. The contribution of unknown types of sources is less then 10 % at 100 GeV and reduces to zero at E>300 GeV. Taking into account that the point-spread function in this energy range has a narrow width, we have verified the point-source flux calculation by summing the photons within circles of 1 degree around the catalog source positions and estimating the remaining diffuse flux from the photon counts outside the source circles.  






§ ICECUBE UPPER LIMIT ON THE GALACTIC PLANE SIGNAL



The upper limit on the Galactic Plane neutrino signal shown in the bottom panel of Fig. 1 is plotted as a envelope curve of the 90% confidence level upper limits on the power-law type spectra as derived in Ref. <cit.>. The envelope curve is tangent to the set of straight lines representing upper bounds on the flux for different values of the power-law slope. 





