Implications of gamma-ray observations on proton models of UHECR
The origin of ultra high energy cosmic rays (UHECR) is still unknown. However, great progress has been achieved in past years due to the good quality and large statistics in experimental data collected by the current observatories. The data of the Pierre Auger Observatory show that the composition of the UHECRs becomes progressively lighter starting from eV up to eV and then, beyond that energy, it becomes increasingly heavier. These analyses are subject to important systematic uncertainties due to the use of hadronic interaction models that extrapolate lower energy accelerator data to the highest energies. Although proton models of UHECRs are disfavored by these results, they cannot be completely ruled out. It is well known that the energy spectra of gamma rays and neutrinos, produced during propagation of these very energetic particles through the intergalactic medium, are a useful tool to constrain the spectrum models. In particular, it has recently been shown that the neutrino upper limits obtained by IceCube challenge the proton models at 95% CL. In this work we study the constraints imposed by the extragalactic gamma-ray background, measured by Fermi-LAT, on proton models of UHECRs. In particular, we make use of the extragalactic gamma-ray background flux, integrated from 50 GeV to 2 TeV, that originates in point sources, which has recently been obtained by the Fermi-LAT collaboration, in combination with the neutrino upper limits, to constrain the emission of UHECRs at high redshits (), in the context of the proton models.
Despite a big experimental effort done in the past years, the origin of the ultra high energy cosmic rays (UHECRs), i.e. with energies larger than eV, is still unknown. The Pierre Auger Auger:15 () and Telescope Array TA:12 () observatories, are currently taking data in this energy range. The Southern hemisphere Pierre Auger Observatory, in the Province of Mendoza, Argentina, is the largest cosmic ray observatory in the world. The Telescope Array observatory, placed in Utah, USA, is the largest in the northern hemisphere. The UHECR energy spectrum presents two main features, the ankle, found at an energy of eV, which consists in a hardening of the flux, and a suppression at eV Valino:15 (); Jui:15 (). Note that, even though the spectra observed by The Pierre Auger and Telescope Array observatories present some differences, these two features are observed in both measurements.
UHECRs generate atmospheric air showers when they interact with the molecules of the atmosphere. These air showers can be observed by using ground detectors and/or fluorescence telescopes. The ground detectors detect the secondary particles, produced during the shower development, that reach the ground. Whereas, the fluorescence telescopes detect the fluorescence light emitted by the interaction of the secondary charged particles with the molecules of the atmosphere. The Pierre Auger and Telescope Array observatories have both types of detectors, and a subsample of the events are observed by both detector types. The arrival direction, the energy, and the composition of the primary UHECR have to be inferred from the air shower observations. In particular, the composition is obtained by comparing observables sensitive to the primary mass with simulations of the atmospheric air showers. This is subject to large systematic uncertainties because the hadronic interactions at the highest energies are unknown. There are models that extrapolate the lower energy accelerator data to the highest energies. Some of these models have been recently updated with the Large Hadron Collider data. As a consequence, the differences among the predictions of the different models decreased but did not disappear.
One of the parameters most sensitive to primary mass is the atmospheric depth of the shower maximum, . This parameter can be reconstructed with the data taken by fluorescence telescopes. The composition analyses based on observed by Auger, done by using the updated hadronic interaction models, show a decrease of the primary mass from eV up to eV and then, beyond that energy, the mass becomes increasingly heavy Porcell:15 (); AugerXmax:14 (). On the other hand, the data obtained by Telescope Array is consistent with protons TAXmax:15 (); Fujii:15 (). However, it has been recently shown that the composition data observed by the two experiments are consistent Abbasi:15 (). It is worth mentioning that, the statistics of Telescope Array is smaller than the one corresponding to Auger.
The interpretation of the UHECR spectrum depends strongly on the mass composition. If protons are the dominant component, the ankle is originated by pair production in the interaction of the protons with the low energy photons of the extragalactic background light (EBL) and cosmic microwave background (CMB), during propagation through the intergalactic medium. The suppression observed at the highest energies can originate by the photopion-production of protons interacting with the CMB photons, by the intrinsic cutoff of in the sources, or by a combination of both effects. The proton model (known as the dip model) of the UHECR spectrum have been extensively studied in the literature (see Refs. Hill:85 (); DeMarco:03 (); Bere:04 (); Bere:05 (); Bere:06 (); Aloisio:07 (); Aloisio:08 ()).
On the other hand, if the UHECR spectrum is composed by heavy nuclei besides protons, the ankle can be interpreted as the transition between the galactic and extragalactic cosmic rays (see Ref. AloisioRev:12 () for a review). However, this possibility is disfavoured by the Auger data Auger:12 (), large scale anisotropy studies show that the extragalactic component should continue below the ankle. Recently, a new light extragalactic component, that dominates the flux below the ankle, originated in different sources than the ones responsible of the flux at the highest energies, has been proposed Aloisio:14 (). This low energy component could also originate in the photodisintegration of nuclei with the photon fields in the vicinity of the acceleration region or the source Unger:15 (); Allard:15 (). In these scenarios the suppression is due to the photodisintegration process undergone by the nuclei when they interact with photons of the EBL and CMB, to an intrinsic spectral cutoff in the injected spectrum, or to a combination of both effects.
Besides composition information, the observation of secondary gamma rays and neutrinos, generated in the interactions undergone by the cosmic rays during propagation through the universe, can constrain the different models of the UHECR spectrum. In particular, a smaller number of gamma rays and neutrinos are predicted in scenarios with larger fractions of heavy nuclei at the highest energies. In Ref. Heinze:15 () the dip model has been rejected at 95% CL by using the upper limit on the neutrino spectrum obtained by the IceCube experiment Ishihara:15 (). In that work the source evolution is assumed to be the one corresponding to the star formation rate times where is the redshift and is a parameter fixed by fitting Telescope Array data. However, it is also shown that models with no emission at cannot be rejected with present neutrino data.
Extragalactic gamma-ray background (EGB) observations impose quite restrictive constrains on models of the UHECR spectrum Bere:75 (); Fodor:03 (); Semikoz:04 (); Ahlers:10 (); Decerprit:11 (); Bere:11 (); Hooper:11 (); Gelmini:12 (). The EGB has recently been measured, from 100 MeV to 820 GeV, by Fermi-LAT Ackermann:15 (). Part of the EGB originates in gamma-ray point sources. The gamma rays generated by the propagation of UHECRs can contribute to a diffuse component. By using the 2FHL catalog Ackermann:16 () of 360 point sources (mostly blazars) detected by Fermi-LAT, it has been shown that of the integrated EGB spectrum from 50 GeV to 2 TeV originates in point sources AckermannPS:16 (). By using this result, in Ref. Liu:16 () it is found that only a group of nearby sources, contributing in the energy range below the ankle, can be responsible for the light component observed in that energy range. Also, in Ref. Kalashev:16 (), a more restrictive upper limit to the energy density of the electromagnetic cascades, that are developed in the intergalactic medium, is found by using this new result.
In this work we study the impact of this new estimation on proton models of UHECRs. Instead of using the upper limit to the energy density found in Ref. Kalashev:16 (), we calculate an upper limit to the integrated gamma-ray flux, in the energy range from 50 GeV to 2 TeV, that does not originate in point sources. This upper limit is used to study the constraints that this new result imposes on proton models of the UHECRs spectrum, assuming a more relaxed source evolution than the one used in Ref. Heinze:15 (). We also compare the results obtained by using the gamma-ray information with the ones obtained by using the neutrino upper limit.
Ii Fit of the UHECR spectrum
The Telescope Array measurement is considered to fit the UHECR flux assuming a pure proton component. This is because, as mentioned above, the Telescope Array data are compatible with a proton composition at the highest energies. Following Refs. Kido:15 (); Heinze:15 () the fits are performed for energies above eV.
It is assumed that the sources are uniformly distributed in the Universe and that the injected spectrum follows a power law with an exponential cutoff,
where is a normalization constant, is the spectral index, is the cutoff energy, and parametrize the source density and luminosity evolution of the injection with redshift . The evolution function is not known, it depends on the types of sources responsible of the acceleration of the cosmic rays, which are still unknown. The most common options found in the literature correspond to the evolution of the star formation rate (SFR), the gamma-ray bursts (GRB), and the active galactic nuclei (AGN) (see for instance Refs. Heinze:15 (); Gelmini:12 (); Kalashev:13 ()). In this work a broken power law of is assumed,
where and are free parameters. The break point at is motivated by the SFR, GRB, and some special cases of AGN evolution functions Kalashev:13 (). Note that just the parameter is important for the fit of the spectrum because the contribution of the sources at is negligible for eV. The evolution corresponding to is important for the generation of gamma rays and neutrinos.
The energy spectrum of UHECR at Earth is calculated by using the program TransportCR Kalashev:14 (), which has been developed to solve numerically the transport equations that governs the propagation of the cosmic rays in the intergalactic medium. One of the advantages of this approach is the decrease of the computational time compared with the one corresponding to the methods based on the Monte Carlo technique. This makes it very useful for fitting the cosmic ray energy spectrum. Besides nuclei, the TransportCR code propagates gamma rays and neutrinos produced in the interactions of them with the photon backgrounds. The photopion production is based on the SOPHIA code sophia () and there are several models of the EBL available.
Models with fixed cutoff energy are considered. It is possible to fit also as it is done in Ref. Heinze:15 (), but the computational time increases considerably. Note that the processing time is larger in the case considered in this work because of the calculation of the propagation of gamma rays and neutrinos. For a given cutoff energy the proton flux at Earth is calculated for a discrete grid of the parameters , , and . As mentioned before, the fit of the spectrum depends on the parameters and in the energy range under consideration, then the parameter just modifies the flux of secondary particles without affecting the fit. That is why only some particular values of are considered. The spectrum for any value of energy, , and is interpolated by using a trilinear interpolation algorithm, for fixed values of .
The fit of the Telescope Array spectrum is performed by minimizing the chi-square, which is given by,
where is the interpolated spectrum, with the energy of the -th data point and is a systematic shift in the energy scale (see Ref. Heinze:15 (); Kido:15 ()), is the flux measured corresponding to the -th bin, and is the uncertainty in the determination of the flux also for the -th bin. The last term in the takes into account the systematic uncertainty on the determination of the energy, where (see Ref. Heinze:15 ()). Following Ref. Heinze:15 () and are considered as nuisance parameters and then the profile likelihood technique Agashe:14 () is used to get rid of of them. In this technique the nuisance parameters are replaced by their maximum likelihood estimators. Therefore, only and are the parameters of interest.
The EBL is still quite uncertain, principally because a direct measurement is subject to several foregrounds Cooray:16 (). However, there are different models, available in the literature, for the spectrum as a function of the photon energy and redshift. In Ref. Kido:15 () it was shown that the fit of the cosmic ray energy spectrum, assuming a pure proton composition, depends on the EBL model considered. Moreover, the gamma-ray and neutrino spectra also depend on the EBL model assumed (see Refs. Kalashev:16 () and Stanev:05 () for gamma rays and neutrinos, respectively). In this work we consider the EBL model of Ref. Kneiske:04 (), which is commonly used to model the UHECR flux and also, following Ref. Kalashev:16 (), the one of Ref. Inoue:12 ().
The top panel of Fig. 1 shows the fit of the Telescope Array energy spectrum considering the EBL model of Ref. Kneiske:04 (), no cosmic ray emission above , i.e. , and eV. The Telescope Array spectrum is taken from Kido:15 () and correspond to the spectrum obtained by the surface detectors111Instead of the combined spectrum, reported in Ref. Jui:15 (), the one obtained by the surface detectors, reported and fitted in Ref. Kido:15 (), is considered for the fit because the error bars at low energies can be extracted from the plot in that paper, which is not the case for the plot in Ref. Jui:15 ().. The energy scale of the proton spectrum has been shifted according to the value obtained. The best fit is obtained for , , and . Note that this results are compatible with the ones obtained in Ref. Kido:15 () even though they do not use the penalty term in the . The bottom panel of Fig. 1 shows the proton, gamma-ray, and neutrino spectra corresponding to the best fit. Also in this case, the energy scale of the proton spectrum has been shifted according to the value obtained. The plot also shows the upper limit on the neutrino flux obtained by IceCube at 90% CL Ishihara:15 (). The data points at low energies correspond to the spectrum of the EGB observed by Fermi-LAT Ackermann:15 (). Note that the best fit is compatible with the neutrino upper limit an also with the EGB. It is worth mentioning that the Telescope Array energy spectrum is extended to energies below eV with the combined spectrum taken from Ref. Jui:15 (). It can also be seen that the fit overshoots the data points below eV. This low energy contribution can be suppressed due to the diffusion of the particles in the intergalactic magnetic field Berezinsky:07 () or by the presence of a low energy cutoff in the injection spectrum (see Ref. Heinze:15 () for details).
As mentioned in the introduction, the EGB measurement has important implications on proton models of the UHECR spectrum. By using the results obtained in Ref. AckermannPS:16 (), an upper limit on the integrated flux corresponding to the component of the EGB spectrum that do not originate in point sources is obtained,
at 90% CL. Here the energy range considered for the integral goes from 50 GeV up to 2 TeV. The method developed by Feldman & Cousins Feldman:98 () is used to obtain the upper limit (see Appendix A for the details of the calculation).
For each pair and considered to fit the Telescope Array data, after applying the profile likelihood method to the parameters and as described above, the integral of the gamma-ray energy spectrum is calculated and the region in the plane, defined by,
The top panel of Fig. 2 shows the best fit and the regions of 68.27%, 95.45%, and 99.73% CL for the fit in Fig. 1. Also shown are the allowed regions inferred from the upper limits obtained form Eq. (5) for the case of gamma rays and from the condition , where is the upper limit on the flux obtained by IceCube (see bottom panel of Fig. 1), for the case of neutrinos. It can be seen from the plot that, while the best fit is compatible with the neutrino data, it is in tension with the upper limit obtained from the non-point sources component of the EGB. Moreover, the best fit is also in tension with the upper limit obtained at 99% CL (see Appendix B). In this case the upper limit obtained from gamma-ray data is more restrictive than the one coming from the neutrino data. The bottom panel of the figure shows the result obtained in the same conditions as the ones corresponding to the top panel but for . It can be seen that the best fit and the allowed regions are, as expected, unaltered, but in this case even the region corresponding to 99.73% CL is in tension with the gamma-ray upper limit. Any value of the parameter larger than has associated a larger production of secondary gamma rays and neutrinos. Therefore, the region corresponding to 99.73% CL of the fit is also in tension with at least the gamma-ray upper limit for any value of larger than 1.5.
From Fig. 2 it can be seen that for the case in which the UHECRs are generated in the redshift range from to , the gamma-ray upper limit imposes more restrictive conditions than the ones corresponding to the neutrino upper limit. However, when the sources produce UHECRs beyond the restrictions obtained from gamma-ray and neutrino observations are complementary, i.e. by using both, the neutrino and gamma-ray upper limits, it is possible to enlarge the rejection region. The loss of restrictive power of the gamma-ray upper limit is due to the fact that the attenuation of gamma rays, owing to pair production with the EBL photons, that originate in is larger than the one corresponding to the gamma rays originated in . This can be seen from Fig. 3 where the contributions from and to the total spectra are shown. The case considered corresponds to , , and , which is a point inside the 99.73% CL region of the fit for . This point is also in the rejection region obtained from the gamma-ray data but in the allowed region obtained from the neutrino data. In fact while the neutrino flux increases by a factor of at GeV, when the component corresponding to is added, the integral of the gamma-ray spectrum increases in a factor of .
Figure 4 shows the best fit and the associated gamma-ray and neutrino spectra for the EBL model of Ref. Inoue:12 (), for no UHECR emission above , and for eV. The best fit is obtained for , , and . Although, the best fit parameters are different from the ones obtained for the other EBL model considered before, the differences between Figs. 4 and 1 are small.
The top panel in Fig. 5 shows the best fit parameters and the regions of 68.27%, 95.45%, and 99.73% CL for the fit of Fig. 4. Also shown are the allowed regions inferred from the gamma-ray and neutrino upper limits. Also in this case, the best fit is compatible with the neutrino upper limit but it is in tension with the gamma-ray upper limit at both 90% and 99% CL (see Appendix B). Comparing with the top panel of Fig. 2, it can be seen that, in this case, the upper limits coming from gamma-ray and neutrino data are less restrictive. This is because the spectra of these secondary particles, predicted by using the second EBL model considered, are smaller in almost the entire energy range relevant for each type of particle.
The bottom panel of Fig. 5 shows the confidence regions of the fit and the allowed regions inferred from the gamma-ray and neutrino upper limits for . In this case the 99.73% CL region of the fit is in tension with the neutrino upper limit and marginally with the gamma-ray upper limit. Also in this case the gamma-ray upper limit loss restrictive power due to the attenuation in the EBL, in such a way that, the neutrino upper limit is more efficient to reject all models whose points fall in the 99.73% CL region of the fit.
Although eV is compatible with the Telescope array data Heinze:15 (), the best fit value obtained in Ref. Heinze:15 () is eV with . Therefore, the same analysis done for eV is performed for this new value of the cutoff energy. Figure 6 shows the results obtained by using the EBL model of Ref. Kneiske:04 () and for no UHECR emission beyond . The best fit parameters obtained in this case are , , and . The values obtained in Ref. Heinze:15 () are , , and which are close to the ones obtained in this work. In any case, the differences can be explained by the different EBL model and propagation code used for the two calculations. Comparing the bottom panel of Fig. 6 with the ones corresponding to Figs. 1 and 4 it can be seen that for eV the production of gamma rays and neutrinos is smaller than for the case of eV. This is due to the fact that eV is very close to the threshold energy of photopion production undergone by protons interacting with the CMB photons. Therefore, this channel for the generation of these secondary particles is quite suppressed.
The top panel of Fig. 7 shows the best fit parameters and the regions of 68.27%, 95.45%, and 99.73% CL for the fit of Fig. 6. Also shown are the allowed regions inferred from the gamma-ray and neutrino upper limits. In this case, the best fit is compatible with the neutrino upper limit and it is marginally compatible with the gamma-ray upper limits. As expected, in this case the gamma-ray and neutrino upper limits are less restrictive than for eV. It is worth mentioning that the constraints imposed by the neutrino upper limit are less affected than the ones corresponding to the gamma rays, by the decrease of the cutoff energy. This is due to the fact that, as can be seen from the bottom panel of Fig. 6, the neutrino energy distribution becomes narrower than for the eV case and then the neutrino flux at the peak does not decrease considerably. The bottom panel of Fig. 7 shows the results corresponding to in which it can be seen that, for this value, the region of 99.73% CL is in tension with the neutrino upper limit. In this case the neutrino upper limit is more restrictive than the one corresponding to the gamma-ray observations.
Also for eV a similar behavior is obtained when the EBL model of Ref. Inoue:12 () is considered. In this case, the best fit parameters are , , and . For no UHECR emission above the best fit parameters fall in the allowed regions inferred from the gamma-ray and neutrino upper limits and the 99.73% CL region of the fit is in tension with the neutrino upper limit for .
In this work we have studied the constraints, imposed by the results recently obtained by the Fermi-LAT collaboration about the EGB origin, on proton models of UHECRs. For that purpose, we have first calculated an upper limit to the component of the integrated EGB flux, from 50 GeV to 2 TeV, that does not originate in point sources. The obtained value at 90% CL is .
We have assumed that the UHECR emission as a function of redshift follows a broken power law in , with a breaking point at and an end point at . We have found, by fitting the TA data, a very fast increase of the UHECR emission for , such that the index of the power law takes the values , depending on the cutoff energy and EBL model considered. This result is consistent with previous work. We have also found that for eV the best fit, corresponding to the case in which there is no UHECR emission beyond , is rejected at 99% CL by using the inferred upper limit from gamma-ray observations. However, part of the 68.27% CL region of the fit is in the allowed region inferred from that upper limit. For eV the best fit is in the allowed regions obtained from the upper limits inferred from gamma-ray and neutrino observations.
When the UHECR emission beyond is included the gamma-ray information becomes less restrictive. This is due to the larger attenuation of the gamma-ray flux originated at in the energy interval considered for the integration. The process responsible for this attenuation is pair production on the EBL. In this case, the gamma-ray and neutrino observations become complementary. By using both types of observations we have found that the index of the evolution function for should be smaller than in order not to be in tension with the 99.73% CL region of the fit. Therefore, only proton models of the UHECRs with a much slower redshift evolution in the interval , compared with the the one corresponding to , are still compatible with the neutrino and gamma-ray upper limits. A more precise determination of the different EGB components or a more restrictive neutrino upper limits are required to reject, at a given significance level, all proton models of UHECRs by using these types of analyses, which are independent of composition measurements.
Note added: as our paper was completed another manuscript appeared on the arXiv, considering the impact of the EGB measurements on proton models of UHECR, which contains complementary results to ours BereKala:16 ().
Acknowledgements.A. D. S. is a member of the Carrera del Investigador Científico of CONICET, Argentina. This work is supported by CONICET PIP 2011/360, Argentina. The author thanks to O. Kalashev and E. Kido for their help with the TransportCR program and the members of the Pierre Auger Collaboration for useful discussions.
Appendix A Upper limit on the integrated EGB flux that do not originate in point sources
In this appendix the calculation of the upper limit on the EGB flux component that does not originate in point sources, integrated from 50 GeV to 2 TeV, is presented. In Ref. AckermannPS:16 () it is shown that the integrated EGB flux, corresponding to the energy range from 50 GeV to 2 TeV and originated in point sources is cmssr. The total integrated EGB flux corresponding to the same energy range is calculated in Ref. Bechtol:15 () and it is given by cmssr. Because has asymmetric errors the method developed in Ref. Barlow:04 () is used to find an approximated likelihood for the estimators and , which is given by,
where cmssr, cmssr, and cmssr. Here,
where cmssr and cmssr. The integrated EGB flux that do not originate in point sources is obtained by subtracting the contribution of the point sources to the total one, i.e. . Therefore, introducing in Eq. (6) and applying the profile likelihood method Agashe:14 () to get rid of the nuisance parameter an approximated likelihood function for is obtained. As expected, the minimum of is attained at cmssr. Figure 8 shows as a function of , where is the value of the likelihood at the maximum and cmssr. As can be seen from that figure the likelihood obtained is quite symmetric, this is due to the fact that the difference between and is not too big and that the likelihood corresponding to the point sources is combined with the one corresponding to the total EGB integrated flux, which is symmetric. In order to simplify the calculations a Gaussian approximation is developed where the mean value is given by cmssr and cmssr. The Gaussian approximation of the profile likelihood is also shown in the figure, it differs in less than in the two sigma region. In any case, it is chosen such that the Gaussian approximation is always smaller than the profile likelihood, which means that the upper limits obtained below are conservative.
The Feldman & Cousins Feldman:98 () method is applied in order to decide whether an upper limit or an interval containing the true value of , corresponding to a given CL, is obtained. By considering the following variable,
the problem is reduced to the one corresponding to a Gaussian distribution with unknown mean, , and standard deviation equal to 1. Figure 9 shows the confidence belts for 90% and 99% CL.
In this case,
From Fig. 9 it can be seen that, for this value, an upper limit is obtained for both, 90% and 99% CL. The results are,
Appendix B Rejection regions inferred from the gamma-ray upper limit at 99% CL
Figure 10 shows the rejection regions obtained from the upper limit on the integrated EGB flux, corresponding to the energy range from 50 GeV to 2 TeV, that do not originate in point sources at 90 and 99% CL (see Appendix A). The cases considered correspond to eV and no UHECR emission beyond . From the plot it can be seen that for both EBL models considered the best fit of the UHECR flux is rejected at 99% CL. Also shown are the rejection regions inferred from the neutrino upper limit.
- (1) A. Aab et al., Nucl. Instrum. Meth. A 798, 172 (2015).
- (2) T. Abu-Zayyad, et al., Nucl. Instrum. Meth. A 689, 87 (2012).
- (3) I. Valiño for the Pierre Auger Collaboration, International Cosmic Ray Conference, The Hague, Netherlands, (2015) (PoS (ICRC2015) 271).
- (4) C. Jui for the Telescope Array Collaboration, International Cosmic Ray Conference, The Hague, Netherlands, (2015) (PoS (ICRC2015) 035).
- (5) A. Porcelli for the Pierre Auger Collaboration, International Cosmic Ray Conference, The Hague, Netherlands, (2015) (PoS (ICRC2015) 420).
- (6) A. Aab et al., Phys. Rev. D 90, 122005 (2014).
- (7) R. Abbasi et al., Astropart. Phys. 64, 49 (2015).
- (8) T. Fujii for the Telescope Array Collaboration, International Cosmic Ray Conference, The Hague, Netherlands, (2015) (PoS (ICRC2015) 320).
- (9) R. Abbasi et al. for the Pierre Auger Collaboration and the Telescope Array Collaboration, JPS Conf. Proc. 9, 010016 (2016).
- (10) C. Hill and D. Schramm, Phys. Rev. D 31, 564 (1985).
- (11) D. De Marco, P. Blasi, and A. Olinto, Astropart. Phys. 20, 53 (2003).
- (12) V. Berezinsky, S. Grigorieva, B. Hnatyk, Astropart. Phys. 21, 617 (2004).
- (13) V. Berezinsky, A. Gazizov, S. Grigorieva, Phys. Lett. B 612, 147 (2005).
- (14) V. Berezinsky, A. Gazizov, and S. Grigorieva, Phys. Rev. D 74, 043005 (2006).
- (15) R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S.I. Grigorieva, and B. Hnatyk, Astropart. Phys. 27, 76 (2007).
- (16) R. Aloisio, V.S. Berezinsky, P. Blasi, and S. Ostapchenko, Phys. Rev. D 77, 25007 (2008).
- (17) R. Aloisio, V. Berezinsky, and A. Gazizov, Astropart. Phys. 39, 129 (2012).
- (18) P. Abreu et al., Astrophys. J. Lett. 762, L13 (2012).
- (19) R. Aloisio, V. Berezinsky, and P. Blasi, JCAP 10, 020 (2014).
- (20) M. Unger, G. Farrar, and L. Anchordoqui, Phys. Rev. D 92, 123001 (2015).
- (21) N. Globus, D. Allard, and E. Parizot, Phys. Rev. D 92, 021302 (2015).
- (22) J. Heinze, D. Boncioli, M. Bustamante, and W. Winter, arXiv:1512.05988.
- (23) A. Ishihara [IceCube Collaboration], “High Energy Neutrino Astronomy”, Talk at the TeVPA Conference (2015).
- (24) V. Berezinsky and A. Smirnov, Astrophys. Sp. Sci. 32, 461 (1975).
- (25) Z. Fodor, S. Katz, A. Ringwald, and H. Tu, JCAP 0311, 015 (2003).
- (26) D. Semikoz and G. Sigl, JCAP 0404, 003 (2004).
- (27) M. Ahlers, L. Anchordoqui, M. Gonzalez-Garcia, F. Halzen, and S. Sarkar, Astropart. Phys. 34, 106 (2010).
- (28) G. Decerprit and D. Allard, Astron. Astrophys. 535, A66 (2011).
- (29) V. Berezinsky, A. Gazizov, M. Kachelrieß, and S. Ostapchenko, Phys. Lett. B 695, 13 (2011).
- (30) D. Hooper, A. Taylor, S. Sarkar, Atropart. Phys. 34, 340 (2011).
- (31) G. Gelmini, O. Kalashev, D. Semikoz, and V. Dmitri, JCAP 01, 044 (2012).
- (32) M. Ackermann et al., Astrophys. J. 799, 86 (2015).
- (33) M. Ackermann et al., Astrophys. J. Suppl. Ser. 222, 5 (2016).
- (34) M. Ackermann et al., Phys. Rev. Lett. 116, 151105 (2016).
- (35) Ruo-Yu Liu, Andrew M. Taylor, Xiang-Yu Wang, Felix A. Aharonian, arXiv:1603.03223.
- (36) V. Berezinsky and O. Kalashev, arXiv:1603.03989.
- (37) E. Kido and O. Kalashev for the Telescope Array Collaboration, International Cosmic Ray Conference, The Hague, Netherlands, (2015) (PoS(ICRC2015)258).
- (38) O. Kalashev, A. Kusenko, W. Essey, Phys. Rev. Lett. 111, 041103 (2013).
- (39) O. Kalashev and E. Kido, Journal of Experimental and Theoretical Physics 120, 790 (2015).
- (40) A. Mücke, R. Engel, J. Rachen, R. Protheroe, and T. Stanev, Comp. Phys. Commun. 124, 290 (2001).
- (41) K. Olive et al. [Particle Data Group Collaboration], Review of particle physics, Chin. Phys. C 38, 090001 (2014).
- (42) A. Cooray, Royal Society Open Science 3, 150555 (2016).
- (43) T. Stanev, D. de Marco, M. Malkan, F. Stecker, Phys. Rev. D 73, 043003 (2005).
- (44) T. M. Kneiske, T. Bretz, K. Mannheim, and D. H. Hartmann, Astron. Astrophys. 413, 807 (2004).
- (45) Y. Inoue et al., Astrophys. J. 768, 197 (2013).
- (46) V. Berezinsky and A. Gazizov, Astrophys. J. 669, 684 (2007).
- (47) G. Feldman and R. Cousins, Phys. Rev. D 57, 3873 (1998).
- (48) K. Bechtol et al., arXiv:1511.00688.
- (49) R. Barlow, arXiv:physics/0406120.
- (50) V. Berezinsky, A. Gazizov, and O. Kalashev, arXiv:1606.09293.