Angular correlation between IceCube high-energy starting events and starburst sources
Starburst galaxies and star-forming regions in the Milkyway, with high rate of supernova activities, are candidate sources of high-energy neutrinos. Using a gamma-ray selected sample of these sources we perform statistical analysis of their angular correlation with the four-year sample of high-energy starting events (HESE), detected by the IceCube Neutrino Observatory. We find that the two samples (starburst galaxies and local star-forming regions) are correlated with cosmic neutrinos at (pre-trial) significance level, when the full HESE sample with deposited energy TeV is considered. However when we consider the HESE sample with deposited energy TeV, which is almost free of atmospheric neutrino and muon backgrounds, the significance of correlation decreased drastically. We perform a similar study for Galactic sources in the 2nd Catalog of Hard Fermi-LAT Sources (2FHL, GeV) catalog as well, obtaining (pre-trial) correlation, however the significance of correlation increases with higher cutoff energy in the HESE sample for this case. We also fit available gamma-ray data from these sources using a interaction model and calculate expected neutrino fluxes. We find that the expected neutrino fluxes for most of the sources are at least an order of magnitude lower than the fluxes required to produce the HESE neutrinos from these sources. This puts the starburst sources being the origin of the IceCube HESE neutrinos in question.
Prepared for submission to JCAP
Angular correlation between IceCube high-energy starting events and starburst sources
Department of Physics, University of Johannesburg,
P. O. Box 524, Auckland Park 2006, South Africa
The origin of cosmic neutrinos in the sample of high-energy starting events (HESE) detected by the IceCube Neutrino Observatory is one of the outstanding puzzles in recent years . The four-year HESE sample contains 54 events, collected within 1347 days of operation from both the Northern and Southern skies, with deposited energy in the detector in the range of 20 TeV to 2 PeV [2, 3]. These events correspond to a detection of cosmic neutrinos above an estimated atmospheric neutrino background events and atmospheric muon background events [2, 3].
Galactic and extragalctic origin of the HESE events from astrophysical sources [5, 4, 6, 7, 8] and from dark matter decay [9, 10, 11] have been widely discussed in literature. Poor angular reconstruction of the cascade-type events, which dominates the HESE sample, is the main hindrance to identify the sources. Searches for significant association of sources from the same direction as the IceCube neutrino events have been done, importantly with Galactic origin, by the Galactic center [12, 13], Fermi bubbles [12, 14] and pulsar wind nebulae . In case of extragalactic sources, IceCube neutrino events have been associated with different kind of blazars [15, 16, 17, 18] and with star-forming galaxies [19, 20].
In our previous study we investigated possible correlation of the three-year HESE sample of 37 events with Ultra-High Energy Cosmic Rays (UHECRs) with energy EeV . We found a (pre-trial) significance of correlation between the HESE sample and the UHECRs with energy EeV. A similar study has also been done in Ref. , finding similar or higher significance of correlation between the UHECRs and IceCube neutrino events. These studies are independent of the astrophysical source models and useful to check a common origin of the UHECRs and cosmic neutrinos. Subsequently we found that a number of nearby () Seyfert galaxies in the Swift-BAT 70 month X-ray source catalog  are within error circles of the EeV CR directions, which in turn are correlated with the three-year HESE neutrinos .
Statistical correlations directly between the astrophysical sources and the HESE neutrinos have also been studied. Correlation with different types of blazars, such as high-energy peaked balazars (HBLs) , three-year Fermi Large Area Telescope (LAT) Active Galactic Nuclei (AGN) Catalog (3LAC) , AGNs in 70 month Fermi-LAT data  and Seyfert galaxies . No significant correlation was found in these studies. A statistically significant correlation, however, was found between the three-year HESE neutrinos and starburst galaxies and regions  in the IRAS (Infra-Red Astronomical Satellite) catalog . This is very interesting from astrophysical perspective since supernova remnants, which has a high rate in starburst sources, are widely thought to be the sources of CRs at least up to PeV . Starburst galaxies have also been proposed as high-energy neutrino sources [28, 29](see, however, [30, 31]).
In this work we study cross correlation between -ray selected starburst galaxies, as well as Galactic sources in the 2nd Fermi Hard source List (2FHL) , and the four-year HESE sample. We also calculate expected neutrino fluxes from the starburst sources which are correlated with the HESE sample, using interaction model to produce observed rays from them. We describe our data samples in Sec. 2 and results of correlation study in Sec. 3. We calculate expected neutrino fluxes and compare with required fluxes to explain neutrino data in Sec. 4. We discuss our results and conclude in Sec. 5.
2 IceCube HESE neutrino and source samples
We consider the four-year HESE sample of 54 neutrinos with deposited energy in the range TeV-2.3 PeV  as well as the highest-energy ( PeV) event which happens to be a track event . Two track events (event ID numbers 28 and 32) are coincident hits in the IceTop surface array and are almost certainly a pair of atmospheric muon background events  and we excluded them from our analysis. Fig. 1 shows sky map of these 53 events in Galactic coordinates with reported angular errors.
Figure 1 also shows the sources samples we have used to do the analysis. The first sample of sources, named as Sample-I consists of 7 Starburst galaxies (sbgs) and 3 galactic TeV detected star-forming regions. We use the same selection criteria for the 7 starburst galaxies from the IRAS catalog  as in Ref. . In particular the flux density Jy at 100 was imposed. Two of the starburst galaxies, namely M 82 and NGC 253, have been detected in TeV ray [39, 40] and appear in the -LAT 3FGL catalog as well. The other two appear only in the -LAT 3FGL catalog , but have not been detected at TeV energies. In another set of source sample, called Sample-II, we included the 4 starburst galaxies from these 7, that are detected by high energy gamma rays in addition to the 3 local starforming regions detected at TeV gamma rays.
We also used 33 Galactic sources that are supernova remnants (SNRs) and pulsar wind nebula (PWNs) from the Second Catalog of Hard -LAT Sources (2FHL) detected between 50 GeV and 2 TeV and called it Sample-III, as these are objects are related to star-forming activities. We have done the analysis for subcategories of these Galactic sources. In particular 12 PWNs called as Sample-IV, 15 SNRs as Sample-V and 6 sources which are either SNR or PWN (SPP) as Sample-VI. Table 1 lists all the source samples used for analysis.
3 Cross correlation study
We use the directional information of the neutrino events and sources to map them as unit vectors, , on the Galactic sphere and find angular separation between a pair of neutrino and source directions as . Cross correlation between the neutrino events and sources is based on the number of such pairs observed in data within a given angular error . For our study, we divided twice the reported angular error of each neutrino event () into concentric rings, each with angular width . We count the number of pairs in a given annular ring for , where This method is different from the minimum or nearest neighbor approach used in our previous paper  and allows to explore cross correlation at various angular distances.
To evaluate the significance of cross correlation between the neutrino events and sources, we generated sets of Monte Carlo source distributions by randomizing their positions (more on this later) and count the number of neutrino-source pairs in a simulated data set and in a given annular ring The p-value is calculated by counting the number of times is greater or equal to in each bin over . The average number of pairs in the Monte Carlo data sets for a given is is the average hits expected for the null distribution. We report the minimum p-value as the significance for that source sample.
We generate Monte Carlo source distributions in two different ways. In the first case, by randomizing the Galactic latitude () for the sources, keeping their longitude () fixed. We call this semi-isotropic null distribution as in Ref. . In the second case, we randomize and also within its allowed range, particularly for Galactic sources we use and for other sources we take . We call this isotropic null distribution. The agreement between these two null distributions is good.
The 54 IceCube events contain 21.6 atmospheric neutrino and muon background events. To minimize its effect on the correlation study we have taken two sets of IceCube neutrino events. One with all 53 events and another with neutrino events having energy > 60 TeV, as at higher energy the background is less. The number IceCube neutrino events in the sample having energy > 60 TeV is 33.
4 Results and Discussion
The results of cross correlation study on starburst galaxies (Sample-I) with the 53 IceCube neutrino events is shown in figure 2. The -isotropic and isotropic null distributions are shown with red and green lines, respectively in the left panels. As we have binned over the area swipe in rings increases with increasing the bin number and so does . The are shown with gray histogram for the 20 bins. For Sample-I within the given angular resolution of the IceCube events and the average expected number of hits is 7.7 and 6.1 for -isotropic and and isotropic null distributions, respectively. In case of the 33 IceCube neutrino events above 60 TeV energy, within the given angular resolution as shown with yellow histogram in the lower panel, while the is nearly 2 for the two null distributions. The lowest p-value that we got for Sample-I is nearly 0.01 for the case of all neutrino events in the bin , for isotropic null distribution. However for the correlation study for events with > 60 TeV there is no significant correlation.
A similar study on the Sample-II is shown in figure 3. For this sample of starforming region we found is 11 over an expected count of 8.1 and 6.5 from -isotropic and isotropic null distributions, respectively in case of all neutrino events within the median angular error reported for IceCube neutrino events. For the events having TeV, is 3 over an expectation of 2.77 and 1.75 from -isotropic and isotropic null distributions, respectively. The lowest p-value is 0.002 in the bin , for the correlation study with all neutrino events for the isotropic null distribution. Like Sample-I we have not found a significant correlation for these sources with the neutrino events above energy 60 TeV.
The cross correlation study of IceCube neutrino events with the 33 Galactic 2FHL sources (Sample-III) is shown in figure 4. The lowest p-value for the correlation study with all neutrino events is 0.017 in the bin, for isotropic null distribution. In this bin the is 15 while the expected hit is nearly 8 as shown in the upper panel of figure 4. However unlike the first two samples the significance of correlation of neutrino events above energy 60 TeV increased as shown in the lower panel of the figure 4. The lowest p-value is 0.0016 in the same bin (i.e., ), for the isotropic null distribution.
We have also studied cross correlation of IceCube neutrino events with different type of sources within the 2FHL Galactic sources. Correlation study of IceCube neutrino events with Sample-IV (12 Galactic PWNs listed in 2FHL) sources is shown in figure 5. The lowest p-value for all neutrino events is 0.025 for semi-isotropic null distribution in the bin shown in the upper panel. While for events above 60 TeV energy the minimum p-value is 0.009 in the same bin, shown in the lower panel.
Similarly the correlation study with 15 SNR sources, Sample-V is shown in figure 6. The significance of correlation is lower compared to Sample-IV. From upper panel, the lowest p-value is 0.11 in the bin . And for the lower panel the lowest p-value, 0.089 occurs in the bin .
For the spp sources from the 2FHL catalog (Sample-VI) we found a more significant p-value, shown in figure 7. The lowest p-value is 0.006 in the bin , from the upper panel. And from the lower panel the lowest p-value, 0.022 occurs in the bin .
We have listed the results of this study in Table 2. The post trial p-value for -isotropic and isotropic null distributions are listed in column 4 and 7 respectively. The post trial p-value is computed as, . Where p is the smaller pre-trial p-value between all events and events for energy > 60 TeV of IceCube search. And the number of trial contains, the number of bins (20) taken within the directional error reported by IceCube and the two samples of IceCube events taken with cut in energy. Although, strictly speaking, an energy cut does not produce two independent neutrino samples but we are cautious here as we first considered events TeV and then all events. This nevertheless gives more conservative results. We found post trial p-value nearly 0.07 only for the Sample-II and III. We have also computed post-trial p-value while combining different samples. Samples-IV, V and VI are not independent but sub-samples of Sample-III. The independent samples are Samples-I, II and III, therefore it may be interesting to compute a global post-trial p-value, which is , for exploring different samples. The range here depends on whether treating 60 TeV energy cut producing two independent samples or not.
A set of sources can be truly associated with the IceCube neutrinos, if the statistic of correlation increases with the increase of neutrino events. To verify whether the starfoming regions are the sources of IceCube neutrino events we have performed the cross correlation study with three year and four year IceCube neutrino sources separately. The results are shown in Fig. 8. The top panel is the result when computed with all the events while the lower panel is for neutrino events with energy more than 60 TeV. The analysis result of the Samples-I and II for 3 has also been done in , and their result is comparable with our analysis for isotropic null distribution, shown in the top-left panel of figure 8. The faded and darker colors are the minimum p-value occurred in the analysis within the 20 bins of for 3 year and 4 year data respectively. We see an increase in significance with increasing neutrino events only for the sbg sources and for all the Galactic sources combined.
5 Star-forming regions with interaction
The rate of core-collapse SN generally follows the star-formation rate , therefore SNRs in star-forming region is expected to accelerate protons and heavy ions to energy similar to the knee region ( eV) of cosmic ray spectrum. High energy gamma rays and neutrinos are suppose to be produce by the interactions between these accelerated protons and surrounding gas in the environment. The neutrino flux is nearly 2/3 of the gamma ray flux from the (proton-proton interaction) channel. Here we have fitted the high energy gamma ray flux produced by the IceCube neutrinos correlated to the star-forming sources, and then estimated the neutrino flux assuming the correlated neutrino event(s) originated from those source(s). In case of pulsar wind nebulae the neutrino production from proton photon interaction can dominate the channel. So here we have done the analysis of interaction of the sbgs, starforming regions and the spp sources that correlated with the IceCube neutrino events listed in Table 3. The parameters used for channel fitting are explained below. Total number of proton at the source per unit energy interval is,
where is the cutoff energy of protons, is the normalization constant and,
Here is the total energy in protons. The photon flux from the interaction given as in ,
Here is the number of pions produced for the given distribution of function, here we have taken it to be 1, is the density of the ambient hydrogen gas, and taken to be 1 . The inelastic scattering cross section for interaction is taken as in , and . is the luminosity distance of the source. We have fitted the high energy gamma ray data available for the individual sources that correlated with the IceCube neutrino events using equation (5.3) by varying parameters , and . The IceCube neutrino flux for the correlated sources is calculated as in  for detection time 1374 days and the HESE event effective area.
We have calculated the luminosity of CR protons, from the sources by, , where is the cooling time of proton. We used corresponding to the pion threshold energy, , for MeV and . We have also calculated the gamma-ray luminosity, from fits to the data, for the above mentioned and , and listed them in Table 3.
5.1 Sample-I and II
Sources from Sample-I and II
Sample-I and II contains in total 10 sources, 7 are starburst galaxies and 3 are local starforming regions. Even though all the 7 sbgs correlated with IceCube neutrino events, only 4 are detected with gamma-rays. So in our interaction analysis we have excluded those 3 sources (M 83, IC 342 and NGC 6946). NGC 1068 and NGC 4945 are not detected with TeV events. We have done the fit of gamma rays from interaction with the Fermi-LAT detected events , shown in the top panel of figure 9. The middle panel of figure 9 shows the Fermi-LAT and H.E.S.S detected gamma rays from NGC 253 (left)  and VERITAS detected gamma rays from M 82 (right) . It also shows the IceCube neutrino flux from their directions as well as the best fit interaction fluxes. All the three local starforming regions correlated with one and more IceCube neutrino events. W49 A is in the same location of SNR W49 B. So we have listed it in the SNR sub catalog. Cygnus Cocoon is again in the same position as Gamma cygni for SNR catalog. However we have shown the result here. Cygnus Cocoon is a starforming region in the Milky Way and detected by Fermi-LAT, ARGO-YBJ  and Milagro , shown in the left of bottom panel of figure 9. 30 Dor C is a starforming region in the Large Magellanic Cloud and detected by H.E.S.S , shown in the right of bottom panel of figure 9.
Sources from Sample-V
Sources from Sample-V cont…
Sample-V contains 15 SNRs from the 2FHL catalog, 10 of these are correlated with the direction of the IceCube neutrino events. The SNRs are surrounded by molecular cloud and high energy cosmic ray protons can interact with these clouds to produce high energy gama rays, like the Fermi-LAT detected gamma rays at the vicinity W28 are explained as in . The result of the interaction analysis for these sources are shown in figure 10-11.
Sources from Sample-VI
Sample-VI contains 6 sources from the 2FHL catalog which are not definitely identified but associated with SNR/PWN. Among these 5 are correlated with the direction of the IceCube neutrino events. W31 is detected by H.E.S.S. , but has not been detected with TeV events, so we have not done analysis for this source. The result of our interaction analysis for the 4 spp sources is shown in figure 12.
The result of the interaction analysis is given in Table 3. The gamma ray fitting parameters, and are listed in column 3 and 4. The cosmic ray and gamma ray luminosity as a result of this fitting are shown in column 5 and 6, respectively. The general cosmic ray luminosity () for sbg sources from this fitting, is within to erg/sec and the photon luminosity is within to erg/sec, which is roughly of the cosmic ray luminosity, as expected from the model. These luminosities decrease for the local starforming regions and the Galactic SNRs, however their ratio still remains at .
With the interaction fitting parameters we calculated the neutrino flux assuming the neutrino energy is 2/3 that of the gamma ray energy. We expected this flux to fit with the IceCube neutrino flux from the direction of the sources. However, for none of the sources the model is able to account for the neutrino flux from the source directions inferred from the 4 year HESE sample.
The astrophysical sources of the IceCube HESE neutrino events is one of the important puzzles since their discovery. A report of significant correlation of the 3 year IceCube neutrino events with starforming host galaxies (sbgs) and local starforming regions  motivated us to revisit the study with more events added in the 4th year of IceCube observation. We carried our study with the sample sets used in  with the cross-correlation method, and found similar result for 3 year IceCube data, shown in figure 8. However, after performing the same study with 4 year data the significance increased for sbg sources but not for the TeV detected sbgs and the local starforming regions. Moreover, we do not find any hint of correlation with the sbg sources and the local starforming regions when considering TeV HESE neutrino sample, which is almost free of atmospheric background. We therefore conclude that the starburst galaxies and the local starforming regions considered in  are most likely not the sources of the IceCube HESE neutrinos.
Furthermore we performed the cross correlation study for the highest energy gamma ray sources in the 2FHL catalog that includes 33 Galactic sources such as the SNRs, PWNs which are associated to star formation. These Galactic sources show a significant correlation with the HESE neutrinos, giving p-value (pre-trial) 0.017 for all neutrino energies and a lower p-value (pre-trial) 0.0016 for neutrinos with energy greater than 60 TeV. A global post-trial p-value for combining the 2FHL Galactic source sample with the other star-forming sources (Samples-I and II) is reduced to . We have also done correlation study with the SNRs, PWNs and spp sources listed in the 2FHL Galactic sources separately but could not find a significant p-value.
For further study of the possible gamma-ray and neutrino emissions from the star-forming sources, we have done interaction modeling for the sources that correlated with the IceCube neutrino events. We fitted the observed high energy gamma rays detected from individual sources by changing different intrinsic parameters for interaction and then calculated the neutrino flux with those parameters. While interactions can reproduce gamma-ray data in many cases, the resulting neutrino flux falls short by at least an order of magnitude from the inferred neutrino flux from the HESE sample.
Therefore, even though we found significance (post-trial) for cross correlation of the IceCube neutrino events with different starforming sources, it is still not sufficient to say that, they are the sources of the IceCube HESE sample. A further increase in number of IceCube data can verify or discard this small significance in future.
We thank Markus Ahlers, Claudio Koper and Cecilia Lunardini for helpful discussion. This work was supported by the National Research Foundation (South Africa), grant no. 87823 (CPRR).
-  L. A. Anchordoqui et al., JHEAp 1-2, 1 (2014) [arXiv:1312.6587 [astro-ph.HE]].
-  F. Halzen 08-13 July, (2015) 25th International Workshop on Weak Interactions and Neutrinos.
-  C. Kopper, ICRC (2015), Netherlands, PoS(ICRC2015)1081
-  K. Murase, M. Ahlers and B. C. Lacki, Phys. Rev. D 88 (2013) no.12, 121301 [arXiv:1306.3417 [astro-ph.HE]].
-  J. C. Joshi, W. Winter and N. Gupta, Mon. Not. Roy. Astron. Soc. 439, no. 4, 3414 (2014) [arXiv:1310.5123 [astro-ph.HE]].
-  S. Troitsky, JETP Lett. 102, no. 12, 785 (2015) [Pisma Zh. Eksp. Teor. Fiz. 102, no. 12, 899 (2015)] [arXiv:1511.01708 [astro-ph.HE]].
-  A. Palladino and F. Vissani, Astrophys. J. 826, no. 2, 185 (2016) [arXiv:1601.06678 [astro-ph.HE]].
-  A. Neronov and D. Semikoz, Phys. Rev. D 93, no. 12, 123002 (2016) [arXiv:1603.06733 [astro-ph.HE]].
-  A. Esmaili and P. D. Serpico, JCAP 1311, 054 (2013) [arXiv:1308.1105 [hep-ph]].
-  A. Bhattacharya, M. H. Reno and I. Sarcevic, JHEP 1406, 110 (2014) [arXiv:1403.1862 [hep-ph]].
-  K. Murase, R. Laha, S. Ando and M. Ahlers, Phys. Rev. Lett. 115, no. 7, 071301 (2015) [arXiv:1503.04663 [hep-ph]].
-  S. Razzaque, Phys. Rev. D 88, 081302 (2013) [arXiv:1309.2756 ].
-  K. Fang, T. Fujii, T. Linden and A. V. Olinto, Astrophys. J. 794 (2014) 2, 126 [arXiv:1404.6237 ].
-  C. Lunardini, S. Razzaque, K. T. Theodoseau and L. Yang, Phys. Rev. D 90 (2014) 2, 023016 [arXiv:1311.7188 ].
-  P. Padovani and E. Resconi, Mon. Not. Roy. Astron. Soc. 443 (2014) 1, 474 [arXiv:1406.0376 ].
-  F. Kraußet al., Astron. Astrophys. 566 (2014) L7 [arXiv:1406.0645 ].
-  S. Sahu and L. S. Miranda, Eur. Phys. J. C 75, 273 (2015) [arXiv:1408.3664 [astro-ph.HE]].
-  R. Moharana, R. Britto & S. Razzaque, ICRC (2015), Netherlands, PoS(ICRC2015)1122
-  L. A. Anchordoqui, T. C. Paul, L. H. M. da Silva, D. F. Torres and B. J. Vlcek, Phys. Rev. D 89 (2014) 12, 127304 [arXiv:1405.7648 ].
-  K. Emig, C. Lunardini and R. Windhorst, JCAP 1512, 029 (2015) [arXiv:1507.05711 [astro-ph.HE]].
-  T. Dahlen et al., Astrophys. J. 613, 189 (2004) [astro-ph/0406547].
-  R. Moharana and S. Razzaque, JCAP 1508, no. 08, 014 (2015) [arXiv:1501.05158].
-  A. Christov et. al., ICRC 2015, Netherlands, PoS(ICRC2015)1082
-  W. H. Baumgartner, J. Tueller, C. B. Markwardt, G. K. Skinner, S. Barthelmy, R. F. Mushotzky, P. Evans and N. Gehrels, Astrophys. J. Suppl. 207, 19 (2013) [arXiv:1212.3336].
-  A. M. Brown, J. Adams and P. M. Chadwick, Mon. Not. Roy. Astron. Soc. 451, no. 1, 323 (2015) [arXiv:1505.00935].
-  D. B. Sanders, J. M. Mazzarella, D. C. Kim, J. A. Surace and B. T. Soifer, Astron. J. 126, 1607 (2003) [astro-ph/0306263].
-  D. F. Torres and L. A. Anchordoqui, Rept. Prog. Phys. 67, 1663 (2004) [astro-ph/0402371].
-  A. Loeb and E. Waxman, JCAP 0605, 003 (2006) [astro-ph/0601695].
-  I. Tamborra, S. Ando and K. Murase, JCAP 1409, 043 (2014) [arXiv:1404.1189 [astro-ph.HE]].
-  K. Murase, D. Guetta and M. Ahlers, Phys. Rev. Lett. 116, no. 7, 071101 (2016) [arXiv:1509.00805 [astro-ph.HE]].
-  K. Bechtol, M. Ahlers, M. Di Mauro, M. Ajello and J. Vandenbroucke, arXiv:1511.00688 [astro-ph.HE].
-  M. Ackermann et al. [Fermi-LAT Collaboration], Astrophys. J. Suppl. 222, no. 1, 5 (2016) [arXiv:1508.04449 [astro-ph.HE]].
-  C. Kopper, ICRC (2015), Netherlands, PoS(ICRC2015)002
-  M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett. 113, 101101 (2014) [arXiv:1405.5303 [astro-ph.HE]].
-  A. Virmani, S. Bhattacharya, P. Jain, S. Razzaque, J. P. Ralston and D. W. McKay, Astropart. Phys. 17, 489 (2002) [astro-ph/0010235].
-  F. Acero et al. [Fermi-LAT Collaboration], Astrophys. J. Suppl. 218, no. 2, 23 (2015) [arXiv:1501.02003 [astro-ph.HE]].
-  M. Ajello et al. [Fermi-LAT Collaboration], Astrophys. J. 819, no. 2, 98 (2016) [arXiv:1601.06534 [astro-ph.HE]].
-  P. Tam and S. Wagner, Rome, Italy, (2011)
-  V. A. Acciari et al., Nature 462, 770 (2009) [arXiv:0911.0873 [astro-ph.CO]].
-  HESS Collaboration, Science 326 1080 (2009)
-  F. A. Aharonian and A. M. Atoyan, Astron. Astrophys. 309, 917 (1996).
-  G. Di Sciascio [ARGO-YBJ Collaboration], J. Phys. Conf. Ser. 632, no. 1, 012089 (2015)
-  A. A. Abdo et al., Astrophys. J. 753, 159 (2012) [arXiv:1202.0846 [astro-ph.GA]].
-  F. Aharonian [H.E.S.S. Collaboration], Aux. information and data points Astron. Astrophys. 503, 817 (2009) [arXiv:0906.1247 [astro-ph.GA]].
-  F. Aharonian et al. [H.E.S.S. Collaboration], Astrophys. J. 636, 777 (2006) [astro-ph/0510397].
-  F. Aharonian et al. [H.E.S.S. Collaboration], Science 307, 1938 (2005) [astro-ph/0504380].
-  A. Abramowski et al. [H.E.S.S. Collaboration], arXiv:1601.04461 [astro-ph.HE].
-  Y. Tang, C. Yang, L. Zhang and J. Wang, Astrophys. J. 812, no. 1, 32 (2015) [arXiv:1508.07712 [astro-ph.HE]].
-  Ester Aliu, for the VERITAS Collaboration, ICRC proceeding 2011.
-  S. Federici, M. Pohl, I. Telezhinsky, A. Wilhelm and V. V. Dwarkadas, Astron. Astrophys. 577, A12 (2015) [arXiv:1502.06355 [astro-ph.HE]].
-  F. Aharonian [H.E.S.S. Collaboration], Astron. Astrophys. 464, 235 (2007) [astro-ph/0611813].
-  A. A. Abdo et. al., APJ 722 2, (2010) 1303 .
-  F. Brun et. al. for the H.E.S.S Collaboration, , 1303 2011.
-  Tanaka, T., Allafort, A., Ballet, J., et al. The Astrophysical Journal Letters 740 L51 2011.
-  E. Aliu et al., Astrophys. J. 770, 93 (2013) [arXiv:1305.6508 [astro-ph.HE]].
-  A. A. Abdo et. al., [Fermi-LAT Collaboration], APJ, 718 (2010) 348
-  F. Aharonian [HESS Collaboration], Astron. Astrophys. 481, 401 (2008)
-  P. G. Tinyakov and I. I. Tkachev, JETP Lett. 74, 445 (2001) [astro-ph/0102476].
-  S. R. Kelner, F. A. Aharonian and V. V. Bugayov, Phys. Rev. D 74 (2006) 034018 [Phys. Rev. D 79 (2009) 039901] [astro-ph/0606058].
-  W. H. Baumgartner et al., 2013 APJS 207 19 , http://tevcat.uchicago.edu/
-  S. P. Wakely and D. Horan (ICRC2008) 3 1341–1344
-  A. Aab et al. [Pierre Auger Collaboration], Astrophys. J. 804, no. 1, 15 (2015) [arXiv:1411.6111 [astro-ph.HE]].
-  C. Lunardini and S. Razzaque, Phys. Rev. Lett. 108, 221102 (2012) [arXiv:1112.4799].
-  A. Abramowski et al. [HESS Collaboration], Science 347, no. 6220, 406 (2015) [arXiv:1501.06578].
-  B. Bartoli et al. [ARGO-YBJ Collaboration], Astrophys. J. 790, no. 2, 152 (2014) [arXiv:1406.6436].
-  F. Brun et al. [H.E.S.S Collaboration] (2011) [arXiv:1104.5003].
-  J. Aleksić et. al., Astronomy & Astrophysics, 541, 11, 2012.
-  A. A. Abdo et al., The Astrophysical Journal Letters, 700, Issue 2 (2009) pp. L127-L131 .
-  S. Kumar for VERITAS collaboration, 34th International Cosmic Ray Conference (ICRC2015), The Hague (The Netherlands).
-  L. Nava and S. Gabici, Mon. Not. Roy. Astron. Soc. 429, no. 2, 1643 (2013) [arXiv:1211.1668].