Estimate of the fraction of primary photons in the cosmic-ray flux at energies ~10^17 eV from the EAS-MSU experiment data
Yu.A. Fomin, N.N. Kalmykov, G.V. Kulikov, V.P. Sulakov, S.V. Troitsky
EEstimate of the fraction of primary photons in the cosmic-ray fluxat energies ∼ eV from the EAS-MSU experiment data Yu.A. Fomin a , N.N. Kalmykov a , G.V. Kulikov a , V.P. Sulakov a and S.V. Troitsky b a D.V. Skobeltsyn Institute of Nuclear Physics,M.V. Lomonosov Moscow State University, Moscow 119991, Russia b Institute for Nuclear Research of the Russian Academy of Sciences, Moscow 117312, Russia
Abstract
We reanalyze archival EAS-MSU data in order to search for events with anomalouslylow content of muons with energies E µ >
10 GeV in extensive air showers with numberof particles N e (cid:38) × . We confirm the first evidence for nonzero flux of primarycosmic gamma rays at energies E ∼ eV. The estimated fraction of primary gammarays in the flux of cosmic particles with energies E (cid:38) . × eV is (cid:15) γ = (cid:0) . +0 . − . (cid:1) %,which corresponds to the intensity of I γ = (cid:0) . +0 . − . (cid:1) × − cm − s − sr − . The studyof arrival directions does not favour any particular mechanism of the origin of thephoton-like events. The study of the primary mass composition of ultra-high-energy (UHE) cosmic rays (CR)is one of the topical problems of astroparticle physics because these experimental resultsare of crucial importance for understanding the theory of both cosmic-ray generation intheir sources and their subsequent propagation to the Earth. Low UHECR intensity makestheir study by direct methods impossible, so that the only available method is the study ofextensive air showers (EAS).The dominant part of EAS is caused by primary nuclei (from protons to iron), however,there is a considerable interest to possible presence of very different particles, e.g. UHEgamma rays, among them. First works on the subject appeared already a half-centuryago (see e.g. Ref. [1]) but definitive quantitative results are still missing (cf. a review [2]and references therein). Indeed, the highest-energy cosmic photons firmly detected had theenergy of ∼
50 TeV [3]. The searches for gamma rays in the energy ranges 3 × eV (cid:46) E (cid:46) × eV (the EAS-TOP [4], CASA-MIA [5] and KASCADE [6] experiments) aswell as at E (cid:38) eV (the Haverah Park [7], AGASA [8, 9, 10], Yakutsk [11, 12], PierreAuger [13, 14] and Telescope Array [15] experiments) did not find any signal and resultedin upper limits on the photon flux only. A few claims of the experimental detection of1 a r X i v : . [ a s t r o - ph . H E ] N ov eV (cid:46) E (cid:46) eV photons (Mt. Chacaltaya [16], Tien Shan [17], Yakutsk [18] andLodz [19]) had low statistical significance. At the same time, a certain flux of UHE photonsis predicted in many models of both conventional and “new” physics. In particular, the fluxof secondary photons from interactions of extreme-energy particles with cosmic backgroundradiation, the so-called Greizen-Zatsepin-Kuzmin (GZK) photons, may serve as a tool todistinguish various models of cosmic rays at energies (cid:38) × eV because the photonflux is very sensitive to the primary composition at these energies: predominantly lightcomposition at GZK energies results in a much higher flux of secondary photons. Giventhe present contradictory situation with the mass composition at UHE (see e.g. Ref. [20]for a detailed review and Ref. [21] for a brief update), searches for GZK photons are nowconsidered very important. Also, a significant contribution to the UHE gamma-ray fluxis predicted in particular top-down mechanisms of CR origin ([22] and references therein),in particle-physics models with Lorentz-invariance violation [23] and in models with axion-photon mixing [24].One of the most promising approaches to the search of primary gamma rays is the studyof the EAS muon component. The number of muons in a gamma-ray induced EAS is anorder of magnitude smaller than in a usual hadronic shower. Therefore, one may hope tofind photon showers by selecting those which have unusually low muon content.In the present work, we study the muon content of showers with the estimated numberof particles N e > × and zenith angles θ < ◦ detected by the EAS-MSU array [25]in 1982 – 1990. We demonstrate that the number of muonless events exceeds significantlythe background expected from random fluctuations in the development of showers causedby primary hadrons. This result may be interpreted as an indication to the presence ofgamma rays in the primary cosmic radiation with energies of order 10 eV which confirmsand strengthens the first evidence for UHE cosmic photons [26].The rest of the paper is organized as follows. In Sec. 2, we briefly review the experimentalsetup (Sec. 2.1), then discuss the data set we study, and muonless events in particular(Sec. 2.2). Sec. 3 is devoted to the estimate of the number of background muonless eventsfor hadronic showers (Sec. 3.1) and to the derivation of the estimated photon flux underthe assumption that all muonless events not accounted for by the hadronic background arecaused by primary gamma rays (Sec. 3.2). Possible systematic errors in the determinationof the flux are discussed in Sec. 3.3. In Sec. 4, we present a detailed study of the distributionof the arrival directions of muonless events on the celestial sphere and test various modelsof the origin of primary photons. We put our results in the context of the present-day stateof the art and briefly conclude in Sec. 5. The description of the EAS-MSU array is given in [25]. The array had the area of 0.5 km andcontained 77 charged particle density detectors (consisted of the Geiger-Mueller counters) fordetermination of the EAS size N e employing the empirical lateral distribution function [27]2igure 1. The EAS-MSU array setup. Muon detectors (circles) are denoted by µ i , i =1 , . . . ,
4; surface detector stations are represented by squares.and 30 scintillator detectors which measured particle arrival times necessary for determi-nation of the EAS arrival direction. In addition to the surface detectors which recordedmostly electron-photon component of an EAS, the array included also four undergroundmuon detectors, also consisted of Geiger-Mueller counters, located at the depth of 40 metersof water equivalent. These detectors recorded muons with energies above 10 GeV. A muondetector with the area of 36.4 m was located at the center of the array while other threestations had the area of 18.2 m and were located at the distances between 150 m and 300 mfrom the center (see Fig. 1). To select the sample of showers with the number of particles N e > × which we use in this work, 22 scintillator detectors, each of the area of 0.5 m ,were used. The scintillator detector threshold was set at the level of 1 / ∼ ∼ µ s, of at least one ofthe 4-fold coincidence systems.To reduce the number of sub-threshold events which still satisfy the master conditions,the express-analysis of the number of fired Geiger-Mueller counters was invoked. In eachcase, it was required that at least 4 of 22 counters in the selection system recorded thedensity exceeding 1 particle per square meter. With these selection criteria implemented,the probability of detection of a shower with N e > × falling to any place of the arraywas not less than 95%. The position of the shower axis was determined with the precisionof ∼
10 m. The precision of determination of the arrival direction was ∼ ◦ . The number ofparticles in the shower was determined with the accuracy ∼ (15 − R between the shower axisand the muon detector. Line: data; shadow: expectation for hadronic primaries. The presence of muon detectors in the EAS-MSU array allows to search for primary gammarays. The method is based on the fact that, for N e (cid:38) and for an hadronic primary, it ishighly unprobable to have zero muons in the central, 36.4 m , detector if the shower axis iswithin ∼
240 m from it. At the same time, these muonless events are fully consistent withthe conjecture of primary gamma rays. The total number of events with N e ≥ × in thedata set is 1679; of them 48 are muonless.Fig. 2 presents the distribution of muonless events in the distance R between the showeraxis and the muon detector. Most of the muonless events correspond naturally to large R ;however, there are a certain number of events close to the axis which are very difficult toexplain by random fluctuations of the hadronic background. One should note that the realnumber of muonless events is larger than observed because of the non-EAS background whichresults in firing of each counter in the central muon detector with average frequency of 4.6 Hz.In three other muon detectors, the frequency of random firing was 2 to 3 times higher, andin this work, we use only the data of the central detector. It consisted of 1104 counters. Forthe time of EAS detection ∼ µ s, one expects 0.076 random firings. Therefore we assumethat the probability of absence of the random firing was 0.93.To obtain a very rough estimate of the probability to have a muonless hadronic event,one may start with the (experimentally known) mean muon lateral distribution function [27]and estimate the expected muon density ρ µ ( N e , R ) for a given core distance R . Then, bymaking use of the Poisson distribution, one may calculate the probability P ( m = 0) to haveno muons in the detector at this distance. In Fig. 3, the distribution of m = 0 events in P ( m = 0) is shown. The tail at low P ( m = 0) indicates that there might be a problem inexplaining the observed number of muonless events within the standard model of the showerdevelopment. 4igure 3. Cumulative distribution of muonless events in the LDF-based Poisson probability P ( m = 0). To quantify the observed discrepancy more precisely, we performed Monte-Carlo simulationsof proton-induced showers and compared the number of muonless events in data and insimulations.
For the shower simulations, we used the AIRES v. 2.6.0 [28] simulation code, whose choicewas determined primarily by its speed. We used the high-energy hadronic interaction modelQGSJET-01 [29]. The primary protons were thrown with zenith angles 0 ◦ ≤ θ ≤ ◦ andwith energies between 3 × eV and 2 × eV, assuming the integral spectral index of2.0. Without the account of fluctuations, the energy of a N e = 2 × proton shower wouldbe equal to E ∼ eV; however, the fluctuations reduce this value. For the study, showerswith N e ≥ × have been selected; Fig. 4 gives the distribution of the primary energiesof the selected artificial showers. In this way, the total amount of 15000 artificial showerswere simulated. The general assumption behind our estimate of the gamma-ray flux is that all muonlessevents, not accounted for by fluctuations of hadronic showers, are caused by primary gammarays. Therefore, the central moment of the estimate is the calculation of the expected numberof background muonless events from the simulated proton-induced showers.The probability of a zero muon detector reading, m = 0, was estimated under assumption(see Ref. [30] for its motivation) that the fluctuations of muon density in EAS may berepresented as a superposition of (a) fluctuations of muon density at a given distance from5igure 4. Contribution of various primary energies to the proton showers with N e ≥ × .Points: results of the simulation; shadow: naive estimate without the account of fluctuations. E ≈ eV.the shower axis, determined purely by the EAS development, and (b) the Poisson fluctuationsof the number of particles which hit the detector station. In this approach, the probability P (∆ R, S, m ) to have m muons in the detector of area S located in the ring ∆ R (distancebetween R and R + ∆ R from the shower axis), is given by P (∆ R, S, m ) = (cid:90) P EAS (∆ R, M ) · P P ( M, S, m ) dM, where P EAS (∆ R, M ) is the function of the muon number density distribution in the ring asdetermined by the shower development and P P ( M, S, m ) is the Poisson probability to recordexactly m muons in a detector of the area S for the total number M of muons in the ring.Suppose that a shower axis came within the ring ∆ R k = R k +1 − R k from the muondetector. Then the muon density in the ring is determined as ρ µ (∆ R k , i ) = N µ (∆ R k , i ) (cid:0) π (cid:0) R k +1 − R k (cid:1)(cid:1) , where N µ (∆ R k , i ) is the number of muons in this ring and i = 1 , . . . n tot is the number ofthe selected artificial showers. Then the probability of a zero detector reading in the ∆ R k ring is P ( m = 0 , ∆ R k ) = 1 n tot n tot (cid:88) i .
93 exp ( − S cos θ ρ µ (∆ R k , i )) , where θ is the zenith angle of the shower.The total probability of a m = 0 event is P tot ( m = 0) = k max (cid:88) k =1 P ( m = 0 , ∆ R k ) R k +1 − R k R k max , R , m Observed number of P ( m = 0 , ∆ R ) Expected number ofmuonless events muonless events60–90 3 1 . × − × − . × − . × − . × − . × − . × − . × − k max gives the total number of rings considered and the last factor accounts for theprobability for the shower axis to hit the ∆ R k ring. The results of the calculation of prob-ability to observe a muonless event are given, for various distances from the shower axis, inTable 1 together with the number of observed and predicted muonless events in our sampleof 1679 showers.The total probability to have a muonless proton-induced event within 240 m between thedetector and the shower axis is 1 . × − which corresponds to ≈
23 expected muonlessevents in the sample, to compare with 48 observed. As expected, the dominant part of thebackground muonless events should appear in two outer rings we considered, the same beingtrue also for the observed events. However, the total number of the observed events is almosttwice the expected one. This allows one to estimate, based on the Poisson distribution, thenumber S of signal photon-like events in the sample as S = 25 . +7 . − . , which transforms intothe fraction (cid:15) = (cid:0) . +0 . − . (cid:1) % of anomalous muonless events in the sample with N e ≥ × and θ ≤ ◦ .We want to identify the anomalous muonless showers with showers initiated by primaryphotons. To determine the fraction of these events in the energy spectrum of cosmic rays, oneneeds to take into account the difference in the development of showers caused by photonsand protons of the same energy. The gamma-ray showers develop slower in the atmosphereand arrive younger to the surface level (the vertical atmospheric depth for EAS-MSU is1025 g · cm − ). On average, for the primary energies ∼ eV, the number of particles ina gamma-ray shower detected by the EAS-MSU experiment should be ≈ .
86 times largercompared to the proton shower. The cut in N e we use thus corresponds, on average, tothe gamma-ray energy of 5 . × eV. Knowing the total cosmic-ray flux measured by theEAS-MSU array [31], we determine the main result of the present work: the photon fraction, (cid:15) γ = (cid:0) . +0 . − . (cid:1) % for E (cid:38) . × eV , and the photon flux intensity, I γ = (cid:0) . +0 . − . (cid:1) × − cm − s − sr − for E (cid:38) . × eV . (1)7odel N µ N µ (QGSJET-01) Expected number of Excess number ofmuonless events muonless eventsSIBYLL 2.1 [34] 0.70 38.0 10.7QGSJET II-03 [35] 0.89 27.6 21.0QGSJET-01 [29] 1.00 23.5 25.2EPOS 1.99 [36] 1.03 21.9 26.8experiment [37] 1.33 10.9 37.8Table 2. Effect of the choice of hadronic interaction models on the result. The last linecorresponds to experimental results on the muon content of EAS. The systematic uncertainty of our result, within the method we use, is related to the estimateof the number of background muonless events from hadronic showers.
Hadronic interaction models.
The largest uncertainty comes from the variety ofmodels of shower development which predict different values of muon number in EAS. Fur-thermore, this difference is sensitive to the muon threshold energy, which is 10 GeV in ourcase. The change of the mean expected muon density in EAS by ±
10% would result in thechange of the number of background muonless showers in the sample by ±
4. The results wequote are based on the QGSJET-01 model [29] which gives a good description of the LHCand Pierre Auger Observatory measurements of the high-energy hadronic cross section (cf.Fig. 5 of Ref. [21]) and of the LHC multiplicity distributions (see e.g. Ref. [32]); the choiceof the model was also motivated by its computational efficiency. The amount of model-to-model variations of the number of >
10 GeV muons in EAS may be estimated from Ref. [33]and from our own simulations. The effect of the change of the interaction model on our re-sults is summarized in Table 2. We shall note that, according to experimental data on EASdevelopment, all hadronic-interaction models currently in use underestimate the number ofmuons in a shower significantly. In particular, several independent indirect analyses of thePierre Auger Observatory data indicate [37] that the real number of muons is approximately1.5 times larger than predicted by the QGSJET II-03 model. This number is used in Table 2and for the estimate of the systematic error; a similar result was obtained with the help ofmuon detectors of the Yakutsk EAS array [38, 39]. The systematic error in the resultinggamma-ray flux due to the uncertainty of hadronic models is ± Primary composition.
Assumption of the purely proton composition gives a conserva-tive (i.e. large) estimate of the expected background of the muonless events because primaryheavier nuclei produce more muons in EAS. For primary iron, the corresponding number ofmuons is larger by a factor of ∼ . Large fluctuations.
Since no model gives a perfect description of hadron-induced airshowers, and in particular there are large uncertainties in predictions of muon number,one cannot exclude that the fluctuations of the EAS muon content might be much larger8han suggested by simulations. Among theoretical approaches, the probability of occasionalvery low muon density in a proton shower is the highest in the model of Ref. [40] wherethe energy equipartition between positive, negative and neutral components of the cascadewas postulated. As it has been shown in Ref. [41], in the frameworks of this model, it ispossible to obtain the probability of ∼
1% of imitation of a gamma-ray shower by a primaryproton. However, this model is much less physically motivated as compared to those whichare currently used in simulation codes.To summarize the discussion of systematic uncertainties, current experimental and the-oretical understanding of the EAS properties suggests that the flux values we obtain areconservative, though they could become lower if physically less motivated models were usedfor hadronic showers.
In this section, in order to find some hints about the origin of the events we observed, weperform various searches for deviations from isotropy in the distribution of arrival directionsof the photon-like events. All the tests are performed by comparison, by means of a certainstatistical procedure, of the real distribution of arrival directions with the simulated onewhich assumes isotropy. In all cases, the result of a test is given by the probability P that the actual distribution of events is a fluctuation of the isotropic distribution, that isfor small P , the isotropic distribution is excluded at the confidence level of (1 − P ). Fortests of the global (large-scale) isotropy, we use the Kolmogorov-Smirnov method (see e.g.Ref. [42]) which compares one-dimensional distributions of real and simulated events in someobservable (e.g., a celestial coordinate). For searches of the local (small-scale) anisotropy,we rely on the correlation-function method which estimates how often the number of paircoincidences of directions from two catalogs (e.g., one of the arrival directions of cosmic raysand another of particular astronomical objects) in simulated samples exceeds the similarnumber obtained from the real data. The notion of the “pair coincidence” depends on theangular distance ∆ between the directions, so the probability P (∆) is often quoted for acertain range of ∆. The clustering properties of the sample of the directions are estimatedby the same method with both catalogs being identical cosmic-ray lists. More details on themethod may be found e.g. in Ref. [43].In both approaches, we need to simulate sets of arrival directions in the assumption ofthe isotropic flux. These sets should take into account the experimental selection effects.For continuously operating surface detector arrays with the efficiency close to 100%, theexposure is uniform in the azimuth angle and depends on the zenith angle θ by a purelygeometric factor sin θ cos θ , assuming that the incoming flux is isotropic (this is the casewhen one studies the energy-limited sample of cosmic rays). However, our sample is limitedby N e instead of energy and, due to different age of showers coming at different zenith angles,the exposure becomes non-geometric. Based on the observed distribution of θ , we determinethe acceptance factor as ∼ sin θ cos θ . The distribution of events in the azimuth angle isperfectly consistent with uniform as expected. The distribution of the arrival directions on9igure 5. The distribution of arrival directions of muonless events in the sky (equatorialcoordinates). Grades of shade represent the distribution expected for the isotropic flux.the sky, together with the one expected from exposure for the isotropic flux, is shown inFig. 5. In the study of the arrival directions, we do not include 3 of 48 events observed in1982 for which the determination of geometry is uncertain. Among possible mechanisms of the origin of UHE gamma rays, we consider separately thosewhich do not require deviations from the standard particle-physics and astrophysical concepts(we will call these scenarios conventional) and those which require the presence of particlesand/or interactions beyond the Standard Model of particle physics (these will be called“new-physics” scenarios). Note that high-energy photons interact with cosmic backgroundradiation efficiently. Assuming standard physics, the energy attenuation length for a ∼ eV photon is as low as ∼
35 kpc due to efficient e + e − pair production on the cosmicmicrowave background (CMB). This means that the observed photons were created in theGalaxy unless some new physics is assumed. Conventional scenarios.
Scenario 1. Cosmogenic photons.
UHE cosmic particles experience intense interactionswith cosmic background radiation. For protons with energies above ∼ × ‘eV, theseare dominated by the GZK [44, 45] process of pion production through the ∆ resonance; forlower energies the dominant mechanism is the e + e − pair production. For heavier primariesat E ∼ eV, photodisintegration effectively reduces the propagation effects to those ofprotons of lower energy. The secondary particles from all these interactions (pions, electronsand positrons) are the source of the so-called cosmogenic photons which appear either fromsubsequent pion decays or from inverse Compton scattering of e ± . There are a lot of workson the GZK photons (e.g. Refs. [46, 47]); the key point of interest here is the possibilityto use these ∼ (10 − ) eV gamma rays as a tool to determine the composition ofthe bulk of E ∼ eV cosmic rays; due to the GZK process, the flux of the secondaryphotons would be much higher for super-GZK protons than for heavy nuclei. Given thepresent-day uncertainty in the primary composition at the very end of the CR spectrum, seee.g. Refs. [20, 21], this approach attracts considerable attention, though no sign of the GZK10hotons have been observed yet. The expected flux of the GZK photons at E (cid:46) eVis far too low to explain our result; we are not aware of a calculation of the flux at lowerenergies, nor of the distribution of their arrival directions which however should be close tothe isotropic one. Scenario 2. Direct photons from point sources.
UHE astrophysical accelerators are ex-pected to emit energetic photons born in interactions of charged particles with ambientmatter and radiation. The energy of accelerated particles should therefore exceed the energyof the photons, roughly by an order of magnitude. It is presently unknown whether theacceleration of particles up to ∼ (10 − ) eV may happen in any single object in theGalaxy (that is, within the propagation length of ∼ eV photons). In any case, theseobjects are not expected to be numerous; we therefore expect a certain degree of clusteringof the arrival directions of photons in this scenario. Galactic TeV gamma-ray sources mayrepresent plausible candidates for the UHECR accelerators; in this case, the arrival directionswould concentrante around them. “New-physics” scenarios. Scenario 3. Superheavy dark matter.
While the Large Hadron Collider failed to discovereasily any dark-matter candidate, models of dark matter which are beyond the reach ofthis machine are becoming more and more popular. In particular, the superheavy (mass M (cid:38) eV) dark matter (SHDM) scenario, originally put forward [48] to explain theapparent excess of E (cid:38) eV cosmic rays (presently disfavoured), has its own cosmologicalmotivation. Its important prediction is a significant fraction of secondary photons amongthe decay products of these superheavy particles; these energetic photons contribute to theUHECR flux. For M (cid:38) eV, the scenario is constrained, but not killed [49], by the UHEphoton limits; constraints for lower M have not been studied. A characteristic manifestationof this mechanism is the Galactic anisotropy [50] of the arrival directions of photons relatedto a non-central position of the Sun in the Galaxy. Scenario 4. Axion-like particles and BL Lac correlations.
The UHECR data set with thebest angular resolution ever achieved (0 . ◦ ), that of High Resolution Fly’s Eye (HiRes) in thestereo mode, demonstrated hard-to-explain correlations of arrival directions of E (cid:38) eVevents with distant astrophysical sources, BL Lac type objects [51, 52], which suggest that ∼
2% of the CR flux at these energies are neutral particles arriving from these objects. Theonly self-consistent explanation of this phenomenon [24] which does not require violationof the Lorentz invariance suggests that the observed events are caused by the gamma rayswhich mix with hypothetical new light particles (axion-like particles, ALPs) in the cosmicmagnetic fields. This would allow them to propagate freely through the cosmic photonbackground in the form of the inert ALP and then to convert back to real photons in aregion of the magnetic field close to the observer. This approach may also explain someother astrophysical puzzles. A test of this scenario may be performed by cross-correlation ofthe arrival directions with the same BL Lac catalog as in Refs. [51, 52].
Scenario 5. Lorentz-invariance violation.
There is no lack of theoretical models with tinyviolation of the relativistic invariance on the market. In some of them, this effect results inefficient increase of the mean free path of an energetic photon through CMB [23]. Thoughthese models have many free parameters, with no particularly motivated choice, one may11igure 6. The autocorrelation test: the probability P (∆) to have the observed or highernumber of pairs of events within the angular bin ∆ as a fluctuation of the isotropic distri-bution.expect that a possible effect of this change of the attenuation length would be to increasethe cosmogenic-photon flux at E (cid:46) eV by orders of magnitude. There is no evidentsignature of this scenario in arrival directions.
1. Point sources or diffuse?
The test of the presence of a relatively small number ofpoint sources is provided by the autocorrelation function. We present the results in Fig. 6where the probability that the observed excess of pairs of events in the angular bin ∆ isplotted as a function of ∆. No significant clustering of events is found.
2. Test of scenario 2: Galactic TeV sources.
Fig. 7 represents the P (∆) function forcross-correlations of the arrival directions of the photon-like events with positions of GalacticTeV sources from the TeVCat catalog [53], as of May 2013. No sign of correlation is seen.
3. Test of scenario 2: Galactic-plane correlation.
Galactic UHECR acceleratorsof a yet unknown type are still expected to concentrate along the Galactic plane, and thedistribution of events in the Galactic latitude b is a model-independent test of this scenario.Fig. 8 illustrates that the distribution of the photon-like events in b is consistent with thatexpected for an isotropic flux (the Kolmogorov-Smirnov probability P KS ≈ .
4. Test of scenario 3: Galactic anisotropy.
The SHDM-related Galactic anisotropyshould reveal itself in the dipole excess seen in distribution of events in the distance to theGalactic Center. Fig. 9 demonstrates that no such excess is seen ( P KS ≈ .
5. Test of scenario 4: BL Lac correlations.
The HiRes BL Lac correlations [51]appeared as an excess of events close to positions of 156 bright BL Lac type objects selectedfrom the catalog [54] by the cut on the optical magnitude
V < m . A subsequent study [52]suggested also a correlation with TeV-selected BL Lacs. In Fig. 10, we present the results12igure 7. The test of correlation with Galactic TeV sources: the probability P (∆) to havethe observed or higher number of events within the angular distance ∆ from TeVCat [53]Galactic TeV sources as a fluctuation of the isotropic distribution.Figure 8. The distribution of the observed photon-like events (line) and Monte-Carloisotropic events (shadow) in the Galactic latitude b .13igure 9. The distribution of the observed photon-like events (line) and Monte-Carloisotropic events (shadow) in the angular distance to the Galactic Center.Figure 10. The test of correlation with BL Lac type objects: the probability P (∆) to havethe observed or higher number of events within the angular distance ∆ from bright V`eronBL Lacs (sample of Ref. [51], full line) and TeVCat [53] TeV BL Lacs (dashed line) as afluctuation of the isotropic distribution. 14igure 11. The diffuse cosmic photon integral flux versus the photon minimal energy. Theresult of the present work is shown as a cross whose vertical line represents the error bars.Tentative detections and upper limits from other experiments are indicated by symbols: star(Tien Shan [17], detection), open star (Lodz [19], detection), gray triangle (EAS-TOP [4]),gray squares (CASA-MIA [5]), gray diamonds (KASCADE [6, 55]), triangles (Yakutsk [12]),open diamonds (Pierre Auger [13, 14]), boxes (AGASA [8]), large squares (Telescope Ar-ray [15]).of a similar analysis for our photon-like sample, with the same catalog of 156 BL Lacs andwith an updated list of TeV BL Lacs from TeVCat [53]. No significant correlation is seen. The place of our result among others is rather specific. All previous studies put upper limitson the photon flux or fraction for the primary energy intervals ∼ (10 − × ) eV and (cid:38) eV. The EAS-MSU result, first reported in Ref. [26], therefore represents the firstever statistically significant detection of cosmic photons with energies above ∼
100 TeV. Inthe present work, we performed the first estimate of the gamma-ray flux in the previouslyunstudied energy window (5 × − ) eV and estimated statistical and systematic errorsfor its value. The result is compared with limits obtained by other experiments in Fig. 11(flux) and Fig. 12 (fraction). The fraction estimates should be interpreted with great carebecause they are sensitive to the energy determination of the bulk of hadronic primarieswhich is known to suffer from large systematic uncertainties due to lack of understandingof high-energy hadronic interactions. Contrary, the photon flux estimates are more robustbecause they use the primary gamma-ray energy determination and the exposure of the arrayonly, both quantities being well understood. Therefore, the main result of this paper is theflux estimate, Eq. (1).The interpretation of the result is problematic. Leaving aside the discrepancy with theCASA-MIA result in terms of the (uncertain) gamma-ray fraction, the more robust flux15igure 12. The fraction of gamma-ray primaries in the diffuse cosmic-ray integral flux versusthe photon minimal energy. Notations are the same as in Fig. 11; in addition, more Yakutskresults [11], results from Haverah Park ([7], open squares), from reanalysis of the AGASAdata ([9], the highest-energy square like AGASA) and from a combination of AGASA andYakutsk data ([10], large open square) are shown.estimate, Eq. (1), does not formally contradict to any existing experimental constraint butis clearly in conflict with the general trend observed both at lower and higher energies , seeFig. 11. Within conventional scenarios, these photons cannot travel for longer than a fewdosen Mpc and should therefore be born in the Galaxy. However, we do not see any significantGalactic (nor any other) anisotropy in the distribution of the arrival directions. To add tothe troubles, in some scenarios it would be difficult to avoid a conflict with measurements ofthe ∼ N e (cid:38) shower We note in passing that Eq. (1) agrees well with the early Yakutsk estimate [18] of the ∼ eV photonflux, but that claim was based on observation of one event only. N e limited one), which wouldincrease the statistics for inclined events, together with more precise determination of thearrival directions and reduction of the background by means of two-parametric (e.g. bothmuon number and shower-front curvature) selection of photon-like showers. We leave thesequestions for a future study.We are indebted to O. Kalashev and G. Rubtsov for helpful discussions. The work ofN.K., G.K., V.S. and Yu.F. was supported in part by the RFBR grant 11-02-00544 and bythe Ministry of Science and Education of the Russian Federation, agreement 14.518.11.7046.The work of S.T. was supported in part by the RFBR grants 11-02-01528, 12-02-01203, 13-02-01311, 13-02-01293 and by the Ministry of Science and Education of the Russian Federation,agreements 8525 and 14.B37.21.0457. References [1] G.B. Khristiansen, G.V. Kulikov, Yu.A. Fomin,
Ultra-high-energy cosmic radiation ,Moscow, Atomizdat, 1975.[2] M. Risse and P. Homola, Mod. Phys. Lett. A (2007) 749 [astro-ph/0702632 [ASTRO-PH]].[3] T. Tanimori, K. Sakurazawa, S.A. Dazeley et al. [CANGAROO Collaboration], Astro-phys. J. (1998) L33 [astro-ph/9710272].[4] M. Aglietta, B. Alessandro, P. Antoni et al. [EAS-TOP Collaboration], Astropart. Phys. (1996) 71.[5] M. C. Chantell, C.E. Covault, J.W. Cronin et al. [CASA-MIA Collaboration], Phys.Rev. Lett. (1997) 1805 [astro-ph/9705246].[6] G. Schatz, F. Fessler, T. Antoni et al. [KASCADE collaboration], Proc. 28th ICRC,Tsukuba (2003) 2293[7] M. Ave, J. A. Hinton, R. A. Vazquez et al. Phys. Rev. Lett. (2000) 2244 [astro-ph/0007386].[8] K. Shinozaki, M. Chikawa, M. Fukushima et al. [AGASA Collaboration], Astrophys. J. (2002) L117. 179] M. Risse, P. Homola, R. Engel et al. , Phys. Rev. Lett. (2005) 171102 [astro-ph/0502418].[10] G. I. Rubtsov, L. G. Dedenko, G. F. Fedorova et al. , Phys. Rev. D (2006) 063009[astro-ph/0601449].[11] A. V. Glushkov, D. S. Gorbunov, I. T. Makarov et al. , JETP Lett. (2007) 131[astro-ph/0701245].[12] A. V. Glushkov, I. T. Makarov, M. I. Pravdin et al. , Phys. Rev. D (2010) 041101[arXiv:0907.0374 [astro-ph.HE]].[13] J. Abraham, P. Abreu, M. Aglietta et al. [Pierre Auger Collaboration], Astropart. Phys. (2008) 243 [arXiv:0712.1147 [astro-ph]].[14] P. Abreu, M. Aglietta, E.J. Ahn et al. [Pierre Auger Collaboration], arXiv:1107.4805[astro-ph.HE].[15] T. Abu-Zayyad, R. Aida, M. Allen et al. [Telescope Array Collaboration],arXiv:1304.5614 [astro-ph.HE].[16] K. Suga, Y. Toyoda, K. Kamata et al. , Astrophys. J. (1988) 1036.[17] S. I. Nikolsky, I. N. Stamenov and S. Z. Ushev, J. Phys. G (1987) 883.[18] A.V. Glushkov, N.N. Efimov, N.N. Efremov et al. , Proc. 19th ICRC, La Jolla, (1985)186.[19] J. Gawin, R. Maze, J. Wdowczyk and A. Zawadzki, Canad. J. Phys. (1968) 75.[20] K. -H. Kampert and M. Unger, Astropart. Phys. (2012) 660 [arXiv:1201.0018 [astro-ph.HE]].[21] S. Troitsky, Phys. Usp. (3) (2013), Uspekhi Fiz. Nauk (2013) 323[arXiv:1301.2118 [astro-ph.HE]].[22] P. Bhattacharjee and G. Sigl, Phys. Rept. (2000) 109 [astro-ph/9811011].[23] M. Galaverni and G. Sigl, Phys. Rev. Lett. (2008) 021102 [arXiv:0708.1737 [astro-ph]].[24] M. Fairbairn, T. Rashba and S. V. Troitsky, Phys. Rev. D (2011) 125019[arXiv:0901.4085 [astro-ph.HE]].[25] S. N. Vernov, G. B. Khristiansen, V. B. Atrashkevich et al. , Bull. Russ. Acad. Sci. Phys. (1980) 80 [Izv. Ross. Akad. Nauk Ser. Fiz. (1980) 537].[26] N. Kalmykov, J. Cotzomi, V. Sulakov et al. , Izv. Ross. Akad. Nauk Ser. Fiz. (2009)584 1827] Yu.A. Fomin, N.N. Kalmykov, V.M. Kalmykov et al. , Proc. 28th ICRC, Tsukuba, (2003) 119[28] S.J. Sciutto, AIRES: A system for air shower simulations. Version 2.6.0 , 2002.[29] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, Nucl. Phys. Proc. Suppl. (1997) 17.[30] A. A. Lagutin, V. V. Uchaikin and G. V. Chernyaev, Yad. Fiz. (1987) 757.[31] N.N. Kalmykov, L.A. Kuzmichev, G.V. Kulikov et al. , Moscow Univ. Phys. Bull. (2010) 275.[32] D. d’Enterria, R. Engel, T. Pierog et al. , Astropart. Phys. (2011) 98 [arXiv:1101.5596[astro-ph.HE]].[33] R. Engel, Talk at the International Symposium on Future Directions in UHECR Physics,CERN, 13-16 February 2012.[34] E. -J. Ahn, R. Engel, T. K. Gaisser et al. , Phys. Rev. D (2009) 094003[arXiv:0906.4113 [hep-ph]].[35] S. Ostapchenko, Nucl. Phys. Proc. Suppl. (2006) 143 [hep-ph/0412332].[36] T. Pierog and K. Werner, Nucl. Phys. Proc. Suppl. (2009) 102 [arXiv:0905.1198[hep-ph]].[37] R. Engel [Pierre Auger Collaboration], arXiv:0706.1921 [astro-ph].[38] A. V. Glushkov, I. T. Makarov, M. I. Pravdin et al. , JETP Lett. (2008) 190[arXiv:0710.5508 [astro-ph]].[39] L. G. Dedenko, G. F. Fedorova, T. M. Roganova et al. , J. Phys. G (2012) 095202.[40] V.V. Uchaikin, V.V. Ryzhov, The Stochastic Theory of High Energy Particle Transport ,Novosibirsk, Nauka, Siberian Branch (1988) (in Russian).[41] J. Cotzomi Paleta, PhD Thesis, SINP MSU, Moscow, 2010.[42] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
Numerical Recipes:The Art of Scientific Computing , Cambridge University Press, 2007.[43] P. Tinyakov and I. Tkachev, Phys. Rev. D (2004) 128301 [astro-ph/0301336].[44] K. Greisen, Phys. Rev. Lett. (1966) 748.[45] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. (1966) 78 [Pisma Zh. Eksp. Teor. Fiz. (1966) 114].[46] G. Gelmini, O. E. Kalashev and D. V. Semikoz, J. Exp. Theor. Phys. (2008) 1061[astro-ph/0506128]. 1947] D. Hooper, A. M. Taylor and S. Sarkar, Astropart. Phys. (2011) 340 [arXiv:1007.1306[astro-ph.HE]].[48] V. Berezinsky, M. Kachelriess and A. Vilenkin, Phys. Rev. Lett. (1997) 4302 [astro-ph/9708217].[49] O. E. Kalashev, G. I. Rubtsov and S. V. Troitsky, Phys. Rev. D (2009) 103006[arXiv:0812.1020 [astro-ph]].[50] S. L. Dubovsky, P. G. Tinyakov and I. I. Tkachev, Phys. Rev. Lett. (2000) 1154[astro-ph/0001317].[51] D. S. Gorbunov, P. G. Tinyakov, I. I. Tkachev and S. V. Troitsky, JETP Lett. (2004)145 [Pisma Zh. Eksp. Teor. Fiz. (2004) 167] [astro-ph/0406654].[52] R. U. Abbasi et al. [HiRes Collaboration], Astrophys. J. (2006) 680 [astro-ph/0507120].[53] TeVCat, http://tevcat.uchicago.edu/ .[54] M. P. V´eron-Cetty and P. V´eron, Astron. Astrophys. (2001) 92.[55] M. Risse and G. Rubtsov, Talk at the International Symposium on Future Directionsin UHECR Physics, CERN, 13-16 February 2012.[56] S. Ogio et al. [Telescope Array Collaboration], Talk at the International Symposium onFuture Directions in UHECR Physics, CERN, 13-16 February 2012.[57] A. Etchegoyen [Pierre Auger Collaboration], arXiv:0710.1646 [astro-ph].[58] M. Tluczykont, D. Hampf, U. Einhaus et al. , AIP Conf. Proc.1505