Neutrinos from Extra-Large Hadron Collider in the Milky Way

# Neutrinos from Extra-Large Hadron Collider in the Milky Way

## Abstract

Neutrino telescope IceCube has recently discovered astrophysical neutrinos with energies in the TeV – PeV range. We use the data of Fermi -ray telescope to demonstrate that the neutrino signal has significant contribution from the Milky Way galaxy. Matching the -ray and neutrino spectra we find that TeV-PeV Galactic cosmic rays form a powerlaw spectrum with the slope . This spectral slope is consistent with the average cosmic ray spectrum in the disks of the Milky Way and Large Magellanic Cloud galaxies. It is also consistent with the theoretical model of cosmic ray injection by diffusive shock acceleration followed by escape through the Galactic magnetic field with Kolmogorov turbulence. The locally observed TeV-PeV cosmic ray proton spectrum is softer than the average Galactic cosmic ray spectrum. This could be readily explained by variability of injection of cosmic rays in the local interstellar medium over the past  yr and discreetness of the cosmic ray source distribution.

## I Introduction

Our knowledge of the properties of cosmic rays produced by the star formation in the Milky Way galaxy is both precise and largely incomplete. It is precise because the details of the cosmic ray spectrum and elemental composition are measured with high precision (Olive et al., 2014). And it is incomplete because all the precision measurements refer to only single location within the Galaxy: the position of the Solar system. It is even not obvious what kind of information is contained in such single point measurements. Deflections of cosmic rays by turbulent component of Galactic magnetic field forces them into a random walk through the interstellar medium (ISM) (Berezinsky et al., 1990). As a result, cosmic rays arrive from random directions on the sky. It is clear that some information on the sources is encoded in the cosmic ray spectrum, mass composition and global anisotropy (Olive et al., 2014). However, it is not clear how the spectral slope and the break energies of the locally measured piecewise powerlaw cosmic ray spectrum are related to the properties of source population(s) and their distribution in the Galaxy. The diffusion through the ISM and escape from the Galactic Disk modify the cosmic ray spectrum, in a model-dependent way. The locally measured properties of the spectrum are not necessarily representative for those of the global distribution of Galactic cosmic rays. Instead, they could be determined by the peculiarities of recent injection of particles in the local ISM (Ptuskin et al., 2006; Thoudam & Hörandel, 2012; Sveshnikova et al., 2013; Ahn et al., 2010; Neronov & Semikoz, 2012; Eichler et al., 2013; Pohl et al., 2013; Kachelriess et al., 2015; Savchenko et al., 2015).

An illustration of these uncertainties is found in the model of modification of the cosmic ray spectrum by the propagation effects encoded in an energy-dependent diffusion coefficient (Berezinsky et al., 1990). The slope of the interstellar cosmic ray spectrum is determined by and by the slope of the injection spectrum, . The most commonly considered acceleration mechanism is diffusive shock acceleration (DSA) (Krymski, 1977; Bell, 1978; Malkov & Drury, 2001), which is expected to give a slope . Comparison with the slope of the locally measured cosmic ray spectrum, , points to a value of . This is, however, in tension with the measurements of the ratio of abundances of primary and secondary nuclei which give (Oliva, 2013) and with the measurements of anisotropy of the cosmic ray flux (Blasi & Amato, 2012). The value is also favoured by theoretical models of cosmic ray diffusion through the ISM (Berezinsky et al., 1990).

Recent the data on cosmic ray positrons and antiprotons and also the slope and anisotropy properties of the cosmic ray nuclei spectra suggest that the TeV range cosmic ray spectrum has a strong contribution a supernova which exploded approximately two million years ago within several hundred parsec distance from the Sun (Kachelriess et al., 2015; Savchenko et al., 2015). The data on deposition of isotopes in the deep ocean crust provide information on the occurrence of such nearby supernova events (Ellis et al., 1996) indicate that only one such nearby supernova event has occurred over last 14 million years (Fields et al., 2005; Fry et al., 2015). The slope of the locally observed cosmic ray spectrum in the TeV range might be largely determined by this recent supernova event. The information on the ”characteristic” slope of the Galactic cosmic ray spectrum is completely erased by this last supernova contribution.

Complementary information on the cosmic ray source and propagation parameters is provided by secondary -rays and neutrinos. Contrary to the charged cosmic rays, neutral -rays and neutrinos go straight from their production sites to the Earth. -ray and neutrino signal from individual sources could provide information on the injection spectrum of cosmic rays, while diffuse -ray and neutrino emission from the ISM could provide the data on the propagation of cosmic rays in the Galaxy.

Analysis of the -ray signal from the Galactic Plane of the Milky Way Galaxy and from the disk of the Large Magellanic Cloud (LMC) galaxy suggest that the average slope of the cosmic ray spectra produced by the star formation in these two galaxies is , i.e. harder than that of the locally measured cosmic ray spectrum (Neronov & Malyshev, 2015; Foreman et al., 2015). The slope of the neutral pion decay -ray spectrum produced by the interactions of cosmic rays with such hard slope is about .

This slope is close to that of the slope of the spectrum of the isotropic diffuse -ray background (IGRB)(Abdo et al., 2010; Fermi Collab., 2015). This suggests that significant part of the IGRB might originate from the cosmic ray interactions in the star forming galaxies(Ackermann et al., 2012a; Murase et al., 2013, 2014) which are not individually resolved by the Fermi Large Area Telescope (LAT). The normalisation of neutrino and -ray flux in the GeV energy range is fixed by the known relation between the far infrared and -ray luminosity of galaxies, Ackermann et al. (2012a), and by the known level of extragalactic far infrared background,  erg/(cm s sr) Hauser & Dwek (2001). The most recent calculation of the -ray and neutrino flux from star forming galaxies which takes into account this relation along with the details of cosmological evolution of different types of galaxies Murase et al. (2013, 2014) shows that they provide a significant contribution to the IGRB.

An independent verification of the result on the hard average cosmic ray spectrum resulting from the star formation in the Milky Way and other galaxies could be obtained through the neutrino channel. The slope of the spectrum of neutrino emission from the Milky Way and other star forming galaxies is expected be about , close to the slope of the -ray spectrum.

The IceCube collaboration has recently reported the detection of astrophysical neutrino signal in the energy range from 10 TeV to 2 PeV (IceCube Collab., 2013; Aarsten et al., 2013, 2014, 2014a). The signal forms a powerlaw spectrum , with (Aarsten et al., 2014a).

This slope and normalisation of the neutrino spectrum are obviously inconsistent with the possibility that the observed neutrino flux is produced by cosmic rays in the Milky Way and / or other star forming galaxies, if one assumes that the typical spectrum of the cosmic rays resulting from the star formation is about the locally observed cosmic ray spectrum slope, , and it has a knee feature exactly at the same energy as the knee of the locally measured cosmic ray spectrum (Kachelriess & Ostapchenko, 2014), an assumption was explicitly adopted in most of the previous calculations of the Galactic neutrino fluxStecker (1979); Berezinskii et al. (1993); Tchernin et al. (2013); Neronov et al. (2014); Murase & Ahlers (2014). Different assumptions about the average Galactic cosmic ray slope result in the estimates of the Galactic contribution to the astrophysical neutrino flux which differ by orders of magnitude (see e.g. (Ahlers et al., 2015) for an example of calculation assuming ). However, none of the assumptions about the slope and position of the knee of the Galactic cosmic ray spectrum relying on local measurements is justified.

The measured slope of the astrophysical neutrino spectrum is consistent with the neutral pion decay -ray spectrum of the Milky Way and LMC disks (Neronov & Malyshev, 2015; Foreman et al., 2015), as expected if the typical slope of the cosmic ray spectrum resulting from the star formation is rather than . Below we show that not only the slope, but also the normalisation of the neutrino flux is consistent with this hypothesis. We argue that this suggests the validity of a simple model in which the injection of cosmic rays with the spectrum with the slope suggested by the DSA, is followed by the softening by produced by the escape through the turbulent galactic magnetic fields with Kolmogorov or Iroshnikov-Kraichnan turbulence spectrum. This simple model is valid for the spectrum of cosmic rays averaged over sufficiently large regions of galaxies. It is, however, not valid for the single point measurements of the cosmic ray spectrum, such as the measurements at the position of the Solar system in the Milky Way.

## Ii γ-ray and neutrino data analysis

Our analysis uses all publicly available Fermi / Large Area Telescope (LAT) data (Atwood et al., 2009) collected over the period from August 2008 till June 2014. We have processed the data using the Fermi Science Tools v9r32p51 - a standard software package, provided by the Fermi collaboration to reduce the data, obtained by the Fermi/LAT. We have used the Pass 7 ”reprocessed” event selection.

We have filtered the event lists using the gtselect tool with parameter evclass=3, which leaves the -ray events and rejects most of the cosmic ray background events. To produce the all-sky spectrum shown in Fig. 2, we have used the aperture photometry method2, applied to the full sky. An estimate of the exposure in each energy bin was done using the gtexposure tool with the option apcorr=no (there is no need to correct the exposure for the point spread function for the full sky). To separate the diffuse emission from the point source contributions we subtract the flux in the circles of the radius around point sources from the four-year Fermi catalogue Fermi Collab. (2015). This is sufficient for the analysis of the diffuse emission in the energy band above 10 GeV. In this energy band the point spread function of the LAT is smaller than Atwood et al. (2009).

We have verified that the aperture photometry approach gives the result which is consistent with the results obtained using the likelihood analysis. In particular, the spectrum of the part of the sky, calculated using the aperture photometry method is identical to that reported by Fermi Collab. (2015) over the entire energy range from 100 MeV up to 1 TeV.

To produce the Galactic latitude profiles for the neutrino data, we have binned the neutrino events reported by Aarsten et al. (2014) in Galactic latitude and have estimated the IceCube exposure in each Galactic latitude bin. This was done via a Monte-Carlo simulation of the Galactic latitude distribution of neutrino events driven from an isotropic sky distribution, with account of the declination dependence of the IceCube effective area, derived from the information reported by Aarsten et al. (2014). Having estimated the exposure in each Galactic latitude bin, we have multiplied the model fluxes of the Galactic and isotropic components by the Galactic latitude dependent exposure to estimate the expected number of neutrino counts in each latitude bin.

## Iii γ-ray and neutrino all-sky spectrum

Fig. 1 shows the combined -ray and neutrino all-sky spectrum in a broad GeV-PeV energy range. The statistical errors of the -ray signal are small in all energy bins up to  GeV. The uncertainty of the -ray flux measurement is dominated by the systematic errors (Ackerman et al., 2012b). The IceCube neutrino spectrum is derived from the analysis of (Aarsten et al., 2014a). The dark grey shaded band shows the 68% uncertainty of the flux and slope of the neutrino signal.

From Fig. 1 one could see that the neutrino spectrum lies at the high-energy extrapolation of the -ray spectrum of the entire sky. Thus, not only the slopes, but also the normalisation of the two spectra agree with each other. The uncertainty of the neutrino flux at 100 TeV is just by a factor of . Assuming a negligible uncertainty of the -ray flux at  GeV, one could find that a powerlaw fit to the combined -ray plus neutrino spectrum gives a very precise measurement of the slope of the powerlaw, , with an error . This is due to a very large dynamic range of the energy on which the powerlaw is observed and to the moderate uncertainty of the neutrino flux measurement.

A more precise characterisation of consistency of the -ray and neutrino spectra is given by the model calculation shown in Fig. 1. Red and blue thin solid curves show a model of the neutrino and -ray emission from interactions of cosmic rays with the ISM. The cosmic rays have a powerlaw spectrum with the slope . The -ray and neutrino spectra are calculated using the parametrizations of the pion production spectra by Kelner et al. (2006); Kamae et al. (2006). The model does not fit the -ray data in the energy band below  GeV. This is expected, because in this energy band significant contribution from electron Bremsstrahlung is expected (Ackermann et al., 2012; Neronov & Malyshev, 2015).

Consistency of the combined -ray and neutrino signal with a straightforward model of -ray and neutrino emission from a powerlaw distribution of the parent protons /nuclei suggests s simple model in which both the -ray and neutrino signal are generated by the interactions of cosmic rays produced by the star formation process.

From Fig. 1 one could immediately conclude that in such a scenario contribution of the Milky Way galaxy in the neutrino flux could not be negligible. Indeed, the neutrino flux from other star forming galaxies is constrained by the measurement of an upper limit on the -ray flux of star forming galaxies, which is given by the measured IGRB flux. Assuming that the -ray emission from cosmic ray interactions in star forming galaxies saturates the IGRB measurement, one could estimate the neutrino flux from star forming galaxies via extrapolation of the IGRB spectrum toward higher energies (the red hatched range in Fig. 1. This gives an upper limit of % of the contribution from star forming galaxies other than Milky Way.

## Iv γ-ray and neutrino all-sky anisotropy

If a large part of the astrophysical neutrino signal is of Galactic origin, one expects to find higher signal level at low Galactic latitudes. The neutrino signal indeed shows a hint of anisotropy in the direction of the Galactic Plane (Aarsten et al., 2014), which is consistent with the -ray – neutrino signal correlation. The distribution of neutrinos along the Galactic Plane is consistent with the distribution of the -rays, with higher event statistics observed around the region of Galactic Ridge (Aarsten et al., 2014; Neronov et al., 2014).

To find the anisotropy properties of the Galactic component of the neutrino flux, we use the observed Galactic latitude profile of the -ray emission in the energy band above 300 GeV as a template. This is possible because Fig. 1 suggests that the -ray and neutrino signals are both of hadronic origin and the neutrino signal above 10 TeV could be directly calculated from the -ray signal via s simple powerlaw extrapolation. In the energy band above 300 GeV the -ray signal is free from the extragalactic contribution which is suppressed by the effect of gamma-gamma pair production on the Extragalactic Background Light (Fermi Collab., 2015). Thus, the neutrino signal above 10 TeV is a sum of the anisotropic Galactic component which follows the -ray anisotropy template plus an isotropic extragalactic component.

Fig. 2 shows the Galactic latitude profile of the -ray signal compared to the ”template” profiles of pion decay and inverse Compton (IC) emission calculated via powerlaw extrapolation of the model of Ackermann et al. (2012) to the energy band above 100 GeV. Assuming that the emission from the Galactic Plane (the first bin ) is dominated by the pion decay emission (Ackermann et al., 2012; Neronov & Malyshev, 2015), one could find that the the pion decay component provides a good fit for the profile in the whole latitude range up to , without an additional contribution from the IC or isotropic flux. The sub-dominance of the IC contribution is also supported by the modelling of the diffuse emission in the 100 GeV energy band in the Galactic latitude range reported by Fermi Collab. (2015). This analysis finds that the IC component is sub-dominant even compared to the pion decay component calculated assuming the slope of the cosmic ray spectrum. The dashed histogram in the two panels of Fig. 2 shows the total flux in the latitude bins. One could see that the point source contribution to the all-sky flux is minor. Thus, the observed Galactic latitude profile of the -ray flux above 300 GeV could serve as a proxy for the Galactic latitude profile of the Galactic component of the neutrino flux.

Fig. 3 shows the distribution of the neutrino events in Galactic latitude. It is compared to the model calculation taking into account the dependence of the IceCube exposure on the declination taken from the calculation of Aarsten et al. (2014). One could see that a range of models with large Galactic contribution to the neutrino flux (assumed to follow the Galactic latitude profile of the -ray flux above 0.3 TeV, shown in Fig. 2) is consistent with the data. The statistics of the signal is currently not sufficient for the discrimination of the Galactic and isotropic contributions.

From Fig. 3 one could see that the reduced of the fits of both Galactic latitude distribution models the with 50% and 90% Galactic contribution is . The same would be true also for the isotropic models with 0% Galactic contribution. Discrimination between the dominant Galactic or dominant extragalactic flux models will be possible already with IceCube in its present configuration. Fig. 3 shows that a decrease of the size of the errorbars by a factor of is needed for this. This could be achieved with the increase of signal statistics by a factor of 4, i.e. with a long exposure of IceCube,  yr (compared to the available 3 yr). Precision measurement of each component would, however, require a more powerful detector like IceCube-Gen2 (Aarsten et al., 2014a) or km3net (Bagley et al., 2009).

At the first sight, this result seems to be in tension with the analysis of Ahlers et al. (2015) who finds that the Galactic diffuse emission could not contribute more than 50% of the observed neutrino signal. However, the discrepancy stems from model dependence of the results. One could notice that the result of Ahlers et al. (2015) depends on the choice of the spatial template for the hadronic component of the all-sky -ray flux. It relies on a spatial template derived by the Fermi collaboration (Ackermann et al., 2012) assuming the reference cosmic ray spectrum with the slope 2.7 (rather than 2.58 assumed by Ahlers et al. (2015)). The uncertainty of this template and its dependence on the assumed spectral shape of the template are difficult to estimate. In our analysis we assume that the Galactic longitude profile of the Galactic component of neutrino flux is identical to the profile of the -ray emission in the energy band above 300 GeV. This is motivated by the fact that the isotropic and the inverse Compton components of the -ray flux are suppressed at these energies. In such a way in our analysis we avoid using a model-dependent spatial template of the Galactic component of neutrino flux.

Detection of neutrino and -ray emission from the Galactic cosmic rays with powerlaw index provides an important clue for understanding of formation of Galactic cosmic ray spectrum (Kachelriess et al., 2015; Neronov & Malyshev, 2015). It is consistent with the possibility that the cosmic rays are injected with the powerlaw spectrum with the slope predicted by the shock acceleration models Krymski (1977); Bell (1978); Malkov & Drury (2001), with the energy dependence of the diffusion coefficient as expected for the Kolmogorov turbulence spectrum of the ISM (Armstrong et al., 1995).

Although the slope and normalisation of the Galactic neutrino spectrum are fixed by the -ray data at lower energies, there is no certainty about its high-energy end. The Galactic cosmic ray spectrum should have a knee like suppression at high energies produced by the escape from the Galactic Disk and / or high-energy cut-off due to the absence of sources capable of particle acceleration beyond certain energy range (Berezinsky et al., 1990; Giacinti et al., 2014, 2015). Similarly to the slope of the spectrum, the energy of the knee measured in the local Galaxy is not necessary the same as the average knee energy for the Galactic cosmic rays. The knee feature could be explained within the escape model of formation of cosmic ray spectrum as the energy at which the scattering length of cosmic rays becomes comparable to the correlation length of the Galactic magnetic field Ginzburg & Syrovatskii (1964). This model fits the shapes of the spectra of all elemental groups around the knee energy (Giacinti et al., 2014, 2015). Within the escape model, the knee energy varies across the Galaxy and from one galaxy to another in response to the variations of the properties of the galactic magnetic field. Presence of the knee of the cosmic ray spectrum leads to a suppression of the neutrino flux at the energies of the knee energy. Uncertainty of the average position of the knee in the Galactic cosmic ray spectrum introduces an uncertainty in the shape of the spectra of Galactic -ray and neutrino emission in the PeV energy band (the dashed model line ranges in Fig. 1).

The locally observed TeV – PeV cosmic ray proton spectrum appears ”peculiar” in the sense that its slope is different from the slope of the average Galactic cosmic ray spectrum. This peculiarity could not be related to the specific of the propagation process in the local ISM, because in this case also the spectra of atomic nuclei would be affected. It could also hardly be related to the process of injection from the Galactic cosmic ray sources in general, because in this case the proton spectrum would be systematically softer than the nuclei spectrum everywhere in the Galaxy and this would be visible in the softer -ray / neutrino spectrum of the Galaxy.

The most probable reason for the peculiarity of the locally observed proton spectrum is in the change of the cosmic ray injection rate over the last  yr (Neronov & Semikoz, 2012; Eichler et al., 2013; Kachelriess et al., 2015; Savchenko et al., 2015). Massive injection of cosmic rays in the local ISM in a star formation episode  yr ago could have resulted in an increase of the cosmic ray energy density at all energies. Faster diffusion of higher energy particles should have led to faster ”wash out” of the excess density at high-energies and, as a consequence, to a temporary softening of the local cosmic ray spectrum. Several candidate past events in the local Galaxy, like e.g. the event which produced an expanding ring of molecular clouds, the Gould Belt (Gehlers et al., 2000; Neronov et al., 2012), could be considered. In a similar way, a nearby supernova explosion which occurred several million years ago could have produced an additional distortion in the cosmic ray spectrum. Evidence for such an event could be found through the analysis of the secondary positrons and antiprotons in the cosmic ray spectrum (Kachelriess et al., 2015) and anisotropy (Savchenko et al., 2015). Discreteness of the distribution of the cosmic ray sources and time variability of the star formation rate generically leads to a situation in which the locally measured cosmic ray spectrum is not identical to the average Galactic cosmic ray spectrum (Kachelriess et al., 2015; Savchenko et al., 2015; Neronov & Malyshev, 2015).

## V Conclusions

To summarise, we have shown that the the all-sky neutrino spectrum in the TeV-PeV energy range coincides with the high-energy extrapolation of the all-sky -ray spectrum in the energy band below 1 TeV. Such a consistency of the -ray and neutrino spectra is not unexpected in a simple model where both -rays and neutrinos originate in the decays of neutral and charged pions produced by interactions of cosmic rays in the Milky Way and other galaxies.

Noticing that the -ray all-sky spectrum in the energy band just below 1 TeV is dominated by the Galactic contribution, we were able to estimate the flux and anisotropy properties of the Galactic component of the neutrino spectrum.

The isotropic neutrino flux from the cosmic ray interactions in star forming galaxies other than Milky Way is constrained by the level of the IGRB and could not constitute more than % of the observed neutrino flux. Presence of significant Galactic component of the neutrino flux (at least 50% of the flux) is also consistent with the anisotropy properties of the neutrino signal.

Combination of the -ray and neutrino data suggests that the Galactic cosmic ray spectrum is harder than previously thought, with the slope in the TeV – PeV energy range. The discrepancy between the slopes of the average Galactic cosmic ray spectrum and the locally measured cosmic ray spectrum is readily explained by the discreteness of the cosmic ray source distribution and time variability of the local star formation rate.

### Footnotes

1. http://fermi.gsfc.nasa.gov/ssc/data/analysis/
2. http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/
aperture_photometry.html

### References

1. Olive K.A. et al. (Particle Data Group), Chin. Phys. C, 38, 090001 (2014).
2. Berezinskii V.S., Bulanov S.V., Dogiel V.S., GInzburg V.L., Ptuskin V. S., The Astrophysics of Cosmic Rays, Nothr-Holland, Amsterdam (1990).
3. Ptuskin V.S., Jones F.C., Seo E.S., Sina R., Ad.Sp.Res. 37, 1909 (2006).
4. Thoudam S., Hörandel J.R., MNRAS, 421, 1209 (2012).
5. Sveshnikova L.G., Strelnikova O.N., Ptuskin V.S., Astropart.Phys. 50, 33 (2013).
6. Ahn H.S., et al., Ap.J., 714, L89 (2010).
7. Neronov A., Semikoz D.V., Phys. Rev., D85, 083008 (2012).
8. Eicher D., Kumar R., Pohl M., Ap.J., 769, 138 (2013).
9. Pohl M., Eichler D., Ap.J., 766, 4 (2013).
10. Kachelriess M., Neronov A., Semikoz D., arXiv:1504.06472 (2015).
11. Savchenko V., Kachelriess M., Semikoz D.V., arXiv:1505.02720.
13. Bell A.R., Royal Astronomical Society, 182, 147 (1978)
14. Malkov M., Drury L.O.C. Rep. Progr. Phys. 64, 429 (2001).
15. Oliva A., Proc. 33rd International Cosmic Ray Conference, Rio de Janeiro, id. 1266 (2013).
16. P.Blasi, E.Amato, JCAP 01 11 (2012).
17. Ellis J., Fields B.D., Schramm D.N., Ap.J., 470, 1227 (1996).
18. Fields B.D., Hochmuth K.A., Ellis J., Ap.J., 6210, 902 (2005).
19. Fry H.J, Fields B.D., Ellis J., Ap.J., 800, 71 (2015).
20. Neronov A., Malyshev D., arXiv:1505.07601 (2015).
21. Foreman G., Chu Y.-H., Gruendl R., Hughes A., Fields, B., Ricker P., arXiv:1502.04337 (2015).
22. Abdo A.A., et al., Phys. Rev. Lett., 104, 101101 (2010).
23. The Fermi LAT collab., Ap.J., 799 86 (2015).
24. Ackermann M. et al., Ap.J., 755, 164 (2012).
25. Murase, K., Ahlers, M., Lacki, B.C., Phys.Rev. D88, 121301 (2013).
26. Tambora I., Ando S., Murase K., JCAP, 09, 043 (2014).
27. Hauser M.G., Dwek E., ARA&A, 39, 249 (2001).
28. IceCube Collaboration, Science, 342, 947 (2013).
29. Aartsen M.G. et al., Phys. Rev. Lett., 111, 021103 (2013).
30. Aartsen M.G. et al. [IceCube Collaboration], Phys.Rev.Lett., 113, 101101 (2014).
31. Aarsten M.G. et al., arXiv:1410:1749 (2014).
32. M. Kachelriess and S. Ostapchenko, Phys. Rev. D 90, 083002 (2014)
33. Stecker F.W., Ap.J., 228, 919 (1979).
34. Berezinskii V.S., Gaisser T.K., Halzen F., Stanev T., Astropart.Phys., 1, 281 (1993).
35. Tchernin C., Aguilar J.A., Neronov A., Montaruli T., A&A, 560, A67 (2013).
36. Neronov A., Semikoz D., Tchernin C., Phys.Rev., D89, 103002 (2014).
37. Ahlers, M., Murase K, Phys.Rev. D90, 023010 (2014).
38. Ahlers M., et al., arXiv:1505.03156 (2015).
39. Atwood W.B., et al., Ap.J., 697, 1071 (2009).
40. Fermi Collab., arXiv:1501.02003 (2015).
41. Ackermann M. et al., Ap.J.Supp. 203, 4 (2012a).
42. Kelner S.R., Aharonian F.A., Bugayov V.V., Phys.Rev., D74, 034018 (2006).
43. Kamae T., Karlsson N., Mizuno T., Abe, T., Koi T., Ap.J., 647, 692 (2006).
44. Ackermann M., et al., Ap.J., 750, 3 (2012).
45. Aarsten M.G., et al., arXiv:1412.5106 (2014).
46. Bagley P. et al., km3net Technical Design Report http://www.km3net.org/TDR/TDRKM3NeT.pdf (2009).
47. Armstrong J.W., Rickett B.J., Spangler S.R., Ap.J., 443, 209 (1995).
48. Giacinti G., Kachelriess M., Semikoz D.V., Phys.Rev. D90, 041302 (2014).
49. Giacinti G., Kachelriess M., Semikoz D.V., Phys.Rev. D91, 083009 (2015).
50. Ginzburg V.L., Syrovatskii S.I., The Origin of Cosmic Rays (Pergamon Press, Oxford, 1964).
51. Gehlers N., et al., Nature, 404, 363 (2000).
52. Neronov A., Semikoz D.V., Taylor A.M., Phys. Rev.Lett., 108, 051105 (2012).
72294