Constraints on the flux of $\sim (10^{16} - 10^{17.5})$ eV cosmic photons from the EAS-MSU muon data
Yu. A. Fomin N. N. Kalmykov I. S. Karpikov G. V. Kulikov M. Yu. Kuznetsov G. I. Rubtsov V. P. Sulakov S. V. Troitsky




D.V. Skobeltsyn Institute of Nuclear Physics,
M.V. Lomonosov Moscow State University, Moscow 119991, Russia

D.V. Skobeltsyn Institute of Nuclear Physics,
M.V. Lomonosov Moscow State University, Moscow 119991, Russia

Institute for Nuclear
Research of the Russian Academy of Sciences, 60th October Anniversary
Prospect 7a, Moscow 117312, Russia

D.V. Skobeltsyn Institute of Nuclear Physics,
M.V. Lomonosov Moscow State University, Moscow 119991, Russia

Institute for Nuclear
Research of the Russian Academy of Sciences, 60th October Anniversary
Prospect 7a, Moscow 117312, Russia

Institute for Nuclear
Research of the Russian Academy of Sciences, 60th October Anniversary
Prospect 7a, Moscow 117312, Russia

D.V. Skobeltsyn Institute of Nuclear Physics,
M.V. Lomonosov Moscow State University, Moscow 119991, Russia

st@ms2.inr.ac.ru
Institute for Nuclear
Research of the Russian Academy of Sciences, 60th October Anniversary
Prospect 7a, Moscow 117312, Russia









Results of the search for ∼ (10^16 - 10^17.5) eV primary
cosmic-ray photons with the data of the Moscow State University (MSU)
Extensive Air Shower (EAS) array are reported. The full-scale reanalysis
of the data with modern simulations of the installation does not confirm
previous indications of the excess of gamma-ray candidate events. Upper
limits on the corresponding gamma-ray flux are presented. The limits are
the most stringent published ones at energies ∼ 10^17 eV.



Constraints on the flux of ∼ (10^16 - 10^17.5) eV cosmic
photons from the EAS–MSU muon data
    S. V. Troitsky
    December 30, 2023
===========================================================================================





§ INTRODUCTION


Searches for astrophysical gamma rays with energies 10^15 eV≲
E_γ≲ 10^20 eV attracted considerable attention for
years <cit.>. At these energies, photons interact with the
atmosphere and produce extensive air showers (EAS) which may be detected
by installations studying cosmic rays. The primary motivation for the
studies includes multimessenger astronomy. Photons work as a diagnostic
tool to distinguish between various models of the origin of energetic
cosmic rays and neutrinos. In particular, because of the pair production on
cosmic background radiation, the photon flux from extragalactic sources is
strongly suppressed at sub-PeV to sub-EeV energies, and non-observation of
gamma rays in this energy band would strongly favor the extragalactic
origin of IceCube-detected astrophysical neutrinos over the Galactic
models <cit.>. At higher energies,
non-observation of cosmogenic photons <cit.> would strongly
constrain models with proton composition of cosmic rays with E ≳
10^19.5 eV, see e.g. Refs. <cit.>. Contrary,
observation
of high-energy photons may be a smoking-gun signal of new physics,
including superheavy dark matter <cit.> (see also
Ref. <cit.> and, for a recent reanalysis,
Ref. <cit.>), axion-like particles <cit.>, or
an ultimate test of Lorentz-invariance <cit.>.

For air-shower experiments, the main problem in the photon search is to
separate gamma-ray induced events from the usual, hadron-induced air
showers. One of the best discriminating variables is the muon number of a
EAS, which is much lower in photon-induced events, compared to the bulk of
showers. Indeed, a photon-induced shower develops by means of
electromagnetic interactions mostly, and the only source of muons is
provided by photonuclear reactions, which have a relatively small cross
section. Contrary, in a hadron-induced shower, lots of muons are produced
in decays of π mesons, born in hadronic interactions. Unfortunately,
muon detectors are small or missing in many modern cosmic-ray experiments.
Here, we take advantage of having a large muon detector in the EAS-MSU
experiment and use its data to search for air showers poor in muons.

Previous preliminary studies indicated some excess of gamma-ray candidates
in the EAS-MSU data with respect to the expected background <cit.>. These results invited for a
detailed study of the data with state of art methods. This study has been
performed: a full Monte-Carlo model of the data and the detector was
constructed <cit.> which gives an excellent description of
the surface-detector results <cit.> as well as of the
bulk of muon-detector data <cit.>. Here, we present the
ultimate results of this study and consider muonless photon-candidate
events. We will see that, within the present approach, no excess of these
photon-like events is seen in the data. This allows us to put upper limits
on the primary gamma-ray flux which are world-best at energies E_γ∼ 10^17 eV.

The rest of the paper is organized as follows. In Sec. <ref>, we
briefly describe the installation, the data and the Monte-Carlo
simulations used in this work; references to more detailed descriptions
are given. We derive the main result of the paper, the limits on the
gamma-ray flux, in Sec. <ref>. In Sec. <ref>, we
estimate systematic uncertainties of the result and compare it with other
studies, including previous photon searches in the EAS-MSU experiment. We
briefly conclude in Sec. <ref>.



§ DATA AND SIMULATIONS


The EAS-MSU experiment and the Monte-Carlo model are described in
detail in the previous papers <cit.>.
The installation consists of 76 charged-particle detector stations, in
which multiple Geiger-Mueller counters are located which determine the
number of charged particles (N_e) in a shower. The total number of
Geiger-Mueller counters in the surface part of the installation is
about 10,000. The total area of all Geiger-Mueller counters is ∼
250 m^2; the total area of the installation is ∼
0.5 km^2. In this work, we use the data of the large muon
detector located at the center of the array at a depth of 40 meters of
water equivalent underground corresponding to the threshold energy of
10 GeV for vertical muons. The muon detector had the total area of
36.4 m^2 and consisted of Geiger-Mueller counters with area of
0.033 m^2 each.

For the present study, we select events with the following cuts:

 
  * the event passes the reconstruction procedure and the reconstruction
quality criteria are satisfied <cit.>;
 
  * the age parameter of the EAS is in the range 0.3<S<1.8;
 
  * the reconstructed zenith angle θ is below 30^∘;

  * the distance between the reconstructed shower axis and the array center
where the muon detector is located is R<240 m;

  * the reconstructed shower size is N_e>10^7;

  * out of 32 sections of the muon detector, at least 28 were operational.

The only difference with the cuts used in Ref. <cit.> is
the lower N_e threshold which is chosen to cover a wider range of
energies of potential photons. The extension of the N_e range
  requires a new comparison of data to Monte-Carlo to validate the
  simulations. The corresponding distributions of N_e, shower age
  S, core distance of the array center R, zenith angle θ
  and muon density at 100 meters ρ_μ(100) are shown in
  Appendix.  With the account of the muon detector operation cut, we
are left with 1204 days of data taking within 1984–1990. The data set
contains 3148 events.

Following the previous studies, we consider the following criteria to
select photon-candidate events: the muon detector
is not triggered by the shower.
However, low-energy protons do not always produce a sufficient number of
muons in EAS to activate the muon detector, in particular, at large
distances between the detector and the shower axis.
To evaluate the background of muonless events from hadronic primaries, we
make use of the full Monte-Carlo (MC) model of the EAS-MSU array described
in Ref. <cit.>.
It includes simulations of artificial air showers by the CORSIKA
 7.4001 <cit.> package
with the  QGSJET-II-04 <cit.> high-energy hadronic
 interaction model,
FLUKA 2011.2c <cit.> low-energy hadronic interaction model
and EGS4 <cit.> electromagnetic model.
Artificial EAS are recorded and processed identically to the experimental
data.
The Monte-Carlo simulations have been performed for proton and iron
primaries, see Ref. <cit.> for details; a realistic
composition was assumed to be a mixture of the two. The muon component of
EAS is highly dependent on the primary composition, therefore the number
of muonless background events depends on the assumed proton fraction.
For this study, the primary composition has been
determined  in Ref. <cit.> by fitting the observed distribution
of the muon densities.
The assumed fraction of protons is 46 ± 6 %. We note that, as discussed in
Ref. <cit.>, it agrees well with the composition obtained from
the EAS-MSU surface-detector data <cit.>. More details on
the MC set with hadronic primaries may be found in
Ref. <cit.>.

We also need a MC simulation with gamma-ray primaries which is used to
determine the efficiency of the installation for the gamma-ray detection,
to relate the reconstructed N_e to the primary-photon energy
E_γ and to account for the amount of photon showers which do not
produce photon-candidate (that is, muonless) events. The simulation and
reconstruction of the artificial photon showers was performed in a way
similar to that for hadron-induced events <cit.>. The
total number of simulated independent gamma-induced showers is 300, their
thrown energies follow the E_γ^-1 spectrum with
  lower bound of MC 10^16<E_γ≤
10^17.5 eV, and they are selected at the reconstruction stage to
reproduce the primary spectrum ∼ E_γ^-2, as is customary in
high-energy photon searches. The dependence of the resulting limits on the
assumed spectrum is through the efficiency only and is weak. The
total number of MC realizations of gamma-induced events is 27310 (see
Ref. <cit.> for description of the sampling). Of them,
3898 passed all cuts.




§ THE GAMMA-RAY FLUX LIMIT


The total number of muonless events in the set is 86, while the
expected number of background muonless events from primary hadrons is
80.1. The muon detector core distance distribution of the observed and
expected muonless events is shown in Fig. <ref>.


To study various energy ranges, we consider certain subsamples of the
data. The quantity reconstructed for each shower is N_e, not the
primary energy E; the relation between the two quantities may be
obtained from simulations. The (E-N_e) relations are different for
photons and for hadrons, and it is important to keep track of this
difference for gamma-ray searches, see e.g. Refs. <cit.>.
The mean gamma-ray energy E_γ is related to N_e as

    N_e(E_γ)= 4.1 × 10^-10E_γ/eV,

see Fig. <ref>

for the MC simulated points and the fit.
The condition used to select the data for the search of photons with
E_γ> E_γ^ min is defined as

    N_e>max{ 10^7,   a  N_e(E_γ^ min)   },

where the coefficient a=0.56 was chosen in such a way that at least 90%
of MC photon-induced showers with E_γ>E_γ^ min are
reconstructed with these N_e.

For each energy cut E_γ^ min, we determine the number of
observed muonless events, n_ obs, in the sample, as well as the
expected number of muonless events from the background of
hadron-induced showers, n_ b. No excess of muonless events is
seen, and we estimate the maximum number of photon-induced
events in the sample, n_γ^ FC, by means of the standard
Feldman–Cousins method for the Poisson distribution <cit.>, for
the 90% and 95% confidence levels (CL). The upper limits on the
  photon flux are then estimated as a ratio of n_γ^ FC
  and the effective exposure of the experiment to photons of the given
  energy range. The effective exposure accounts for the fraction of
photon events lost in the reconstruction and, also, of reconstructed
photon events which are not muonless. We estimate the effective exposure as
follows.

The geometrical exposure for the conditions we use (R ≤ R_
  max=240 m, θ≤θ_ max=30^∘) is given by
A_ geom=Ω× S × T, where
Ω=2π(1-cosθ_ max) is the solid angle, S=π
R_ max^2 is the area and T=14060.7 hours is
the on-time of the installation corresponding to the data set
used. Note that the R cut is defined in the plane orthogonal to the
shower axis and therefore Ω is calculated differently from the
conventional case when exposure is determined by the area of the array.

MC photon-induced showers are thrown
in a square with area S_ MC=(280  m)^2=0.3136 km^2 and
with zenith angles up to θ_ MC=35^∘. The corresponding MC
geometrical exposure is then A_ MC=Ω_ MC× S_ MC× T, where Ω_ MC=πsin^2θ_ MC and
θ_ MC=35^∘. We calculate the number n_ pass, 0μ
of events from the MC set which passed all cuts (that is, in particular,
were reconstructed with geometrical properties corresponding to A_
geom), satisfy the criterion (<ref>) and are muonless and divide it
by the number n_ MC of thrown MC events (corresponding to A_
MC). The effective exposure is then given by

    A_ eff=n_ pass, 0μ/n_ MC A_ MC .

 The flux limit is then obtained as

    I_γ=n_γ^ FC/A_ eff,

where n_γ^ FC corresponds to the required confidence level.
Next, we define the exposure correction as a ratio of the effective
exposure to geometrical one:

    ξ = A_ eff/A_ geom .


Note that the exposure correction factor may
exceed unity because Monte-Carlo events are thrown to the area larger
than the installation.
Figure <ref>

presents the exposure correction factor ξ as a function of energy. For
comparison, the reconstruction efficiency for primary photons and protons
is also shown (it is determined in a similar way as ξ but without the
muonlessness condition and with the criterion N_e>10^7 instead of
Eq. (<ref>)). Decline of ξ at higher energies reflects the
fact that the probability for a primary photon to produce a EAS which is
not muonless grows with energy.

Our limits on the integral gamma-ray flux, which represent the main result
of this work, are presented in Table <ref>.




§ DISCUSSION




 §.§ Systematic uncertainties


The main part of systematic uncertainties in the study
comes from the simulation of the background of muonless events
from hadronic showers. Indeed, it is known that hadronic-interaction
models used in the air-shower simulations are not perfect, in
particular in the part related to the description of the muon content
of EAS. While our previous study indicates <cit.> that the
bulk of E>10 GeV muon data of the EAS-MSU experiment is well
described by the QGSJET-II-04 simulations, assuming the primary
chemical composition implied by the surface-detector studies, this
is not directly tested for muon-poor showers. The uncertainties of
the hadronic model, in principle, may reveal itself in
the incorrect estimation of the background. Indeed, we note
(cf. Table <ref>) that for certain energy ranges, the number
of observed muonless events in our sample is smaller than expected
under the background-only hypothesis. Therefore, in addition to the
standard statistical estimate of the upper limit on the gamma-ray
flux, we calculated also the “expected” flux upper limits which
would be obtained if the number of observed muonless events followed
the simulations under the background-only hypothesis. Alternatively,
one may estimate the flux limits which might be obtained under
assumption that the MC model doesn't provide a reliable prediction of the
background. To this end, we use the “data-driven background”, that
relies on the assumption that the correct background is equal to the number of the
observed muonless events. These “expected” and
“data-driven-background” limits on the gamma-ray flux are presented
in Table <ref> and compared in Fig.<ref>.

We also estimate the systematic errors associated with the uncertainty
of the chemical composition. The change of the proton fraction within
it's error ± 6% results in the energy dependent correction of flux
limits. The variation of limits is 21% for the minimum energy and 4%
for maximum energy.




A careful look at the relation between n_ b and n_ obs
reveals one peculiarity which is most probably related to the modelling of
hadronic interactions but, in principle, might also be explained in terms
of the presence of a certain amount of primary photons. The ratio n_
obs/n_ b remains constant, ∼ 0.4, for energies E_γ≳ 8 × 10^16 eV, but quickly raises to ∼ 1 below this
energy. If, due to some systematics in the modelling of the background,
the real n_ b is indeed ∼ 0.4 of the MC one, then we have a
certain excess of muonless events at 2× 10^16 eV≲ E
≲ 8 × 10^16 eV, which may correspond to an excess
of primary photons expected, for instance, in the heavy dark-matter
decay scenario.



 §.§ Comparison with previous EAS-MSU results


Previous analyses of the EAS-MSU data in the same energy range suggested
an excess of muonless events which might be explained by the presence of a
certain amount of photons in the primary cosmic-ray flux at E∼
10^17 eV <cit.>. The
present study does not confirm that claim and puts strong upper limits on
the gamma-ray flux, see Table <ref>. It is therefore important to
understand the differences between this analysis and the previous ones.

With respect to the previous preliminary analysis, this final study
has several important advantages. First, it is based on the new
full Monte-Carlo description of the air showers and of the
installation <cit.>. This uses modern simulation tools
and an analysis technique with the real and simulated events processed
by one and the same reconstruction program. In such a way, we take
into account all possible biases introduced at the reconstruction
step, as well as keep track of rare fluctuations in the EAS
development and registration. Second, the reconstruction program
  has been slightly revised for this study. The main overall effect
of the reconstruction update is that, while the muonless events
remain in the data sample, their reconstructed energies become, for
most of them, lower than before. For these lower energies, the
background of muonless hadron-induced showers is higher, and the same
amount of muonless events does no longer represent an
excess. Third, an additional check of the quality of muon data
  was performed: at least 28 out of 32 sections of the muon detector
  are required to be operational.

 To be
specific, let us consider 48 muonless events with N_e≥ 2× 10^7
studied in Ref. <cit.>. Of them, 28 events have
N_e< 2× 10^7 in the new analysis; 4 events arrived before
1984 (not included in the present data set); 10 events arrived at the
days


excluded from the present analysis because of stricter criteria on the
quality of muon data;


6 muonless events remained in the
data set. In addition, 3 new muonless events joined the data set in the
new reconstruction (they did not pass the cuts in the old one), so the
total number of muonless events with N_e≥ 2× 10^7 is now 9.
This reduced number of photon candidates is below the MC background
expectation of 18.9 events, so no excess is present in the updated data
set. Our new results are compared to previous studies in
Fig.<ref>.


Finally, the efficiency of the muonless detection of photons is
estimated in the present work with the full photon Monte-Carlo. The
numerical values of limits became apparently weaker due to account of
the fact that only about 40% of the photon-induced showers are
registered as muonless.



 §.§ Comparison with other results and possible applications


Many experiments searched for primary photons with the EAS technique and
none has yet found any. Our flux limits are compared to others in
Fig. <ref>.

We see that our limits are similar to those of
the KASCADE-Grande experiment, being world-best for certain energies. At
the same time, one can see that the energy range (10^16-10^18) eV
discussed here is one of the least studied bands. While future studies to
improve the limits are important, already the present data may be used to
constrain astrophysical models with Galactic sources of PeV neutrinos <cit.>
 or decaying dark matter with appropriate mass and sufficiently hard spectrum (see <cit.> for a review).




§ CONCLUSIONS


This work presents the results of the search for primary photons in the
EAS-MSU data. The photon candidate events are defined as ones
  giving
 no signal in the muon detector
of the installation. We made use of the full Monte-Carlo simulation of the
installation and of the updated reconstruction of EAS parameters. Contrary
to the previous analysis of the same data, no evidence was found
for an excess of photon-candidate events, and this fact allowed us to put
upper limits on the diffuse flux of primary gamma rays at energies
above ∼ (10^16-10^17) eV. For certain energies, the limits are
world-best ones. The difference with the
previous study is, mainly, due to the change of reconstruction: the
energies of muonless events moved downwards, while the background of
muonless showers from primary hadrons is higher at lower energies. The
limits obtained in this work may be used in multimessenger astrophysics as
well as for constraining exotic particle-physics models.



§ ACKNOWLEDGEMENTS

ST thanks J.-M. Frère for a discussion on the data-driven background and
the Service de
Physique Théorique at Université Libre de Bruxelles for
hospitality at the final stages of the work. Monte-Carlo simulations have
been performed at the computer cluster of the Theoretical Physics
Department, Institute for Nuclear Research of the Russian Academy of
Sciences. The experimental work of the EAS-MSU group is supported by the
Government of the Russian Federation (agreement 14.B25.31.0010) and the
Russian Foundation for Basic Research (project 14-02-00372). Development
of the analysis methods and application of them to the EAS-MSU data is
supported by the Russian Science Foundation (grant 14-12-01340).




38

1975
G.B. Khristiansen, G.V. Kulikov, Yu.A. Fomin,
 “Ultra-high-energy cosmic radiation”, Moscow, Atomizdat, 1975.

Risse
M. Risse and P. Homola,
  “Search for ultrahigh energy photons using air showers,”
  Mod. Phys. Lett. A 22 (2007) 749
  [astro-ph/0702632].
  
  

Ahlers
 M. Ahlers and K. Murase,
 “Probing the Galactic Origin of the IceCube Excess with Gamma-Rays,”
  Phys. Rev. D 90 (2014)  023010

  [arXiv:1309.4077 [astro-ph.HE]].
  
  

Kalashev-ST-gamma
O. E. Kalashev and S. V. Troitsky,
“IceCube astrophysical neutrinos without a spectral cutoff and 10^15-10^17 eV cosmic gamma radiation,”
  JETP Lett.  100 (2015)   761
   [Pisma Zh. Eksp. Teor. Fiz.  100 (2014)  865]
  
  [arXiv:1410.2600 [astro-ph.HE]].
  
  

BerZats
 J. Wdowczyk, W. Traczyk, C. Adcock and A. W. Wolfendale,
“The possibility of detectable fluxes of cosmic ray γ rays with
energy above 10^19 eV”,
J. Phys. A: Gen. Phys. 4 (1971) L37.

Gelmini
G. Gelmini, O. E. Kalashev and D. V. Semikoz,
  “GZK photons as ultra high energy cosmic rays,”
  J. Exp. Theor. Phys.  106 (2008) 1061
[Zh. Exp. Teor. Fiz. 133 (2008) 1214]
  [astro-ph/0506128].
  
  

Sarkar
D. Hooper, A. M. Taylor and S. Sarkar,
  “Cosmogenic photons as a test of ultra-high energy cosmic ray composition,”
  Astropart. Phys.  34 (2011) 340
  [arXiv:1007.1306 [astro-ph.HE]].
  
  

SHDM
V. Berezinsky, M. Kachelriess and A. Vilenkin,
 “Ultrahigh-energy cosmic rays without GZK cutoff,”
  Phys. Rev. Lett.  79 (1997) 4302
  [astro-ph/9708217].
  
  

SHDM1
 V. A. Kuzmin and V. A. Rubakov,
“Ultrahigh-energy cosmic rays: A Window to postinflationary reheating epoch of the universe?,”
  Phys. Atom. Nucl.  61 (1998) 1028
  [astro-ph/9709187].
  
  

KRT-SHDM
O. E. Kalashev, G. I. Rubtsov and S. V. Troitsky,
   “Sensitivity of cosmic-ray experiments to ultra-high-energy photons:
reconstruction of the spectrum and limits on the superheavy dark matter,”
  Phys. Rev. D 80 (2009) 103006 [arXiv:0812.1020 [astro-ph]].

  

Kuznetsov
 O. K. Kalashev and M. Y. Kuznetsov,
“Constraining heavy decaying dark matter with the high energy gamma-ray limits,”
  Phys. Rev. D 94 (2016)   063535
  
  [arXiv:1606.07354 [astro-ph.HE]].
  

axion
M. Fairbairn, T. Rashba and S. V. Troitsky,
“Photon-axion mixing and ultra-high-energy cosmic rays from BL Lac type
objects – Shining light through the Universe,”
Phys. Rev. D 84 (2011) 125019 [arXiv:0901.4085 [astro-ph.HE]].

  

Lorentz-violation
 M. Galaverni and G. Sigl,
  “Lorentz Violation in the Photon Sector and Ultra-High Energy Cosmic Rays,”
  Phys. Rev. Lett.  100 (2008) 021102
  [arXiv:0708.1737 [astro-ph]].
  
  

GrishaLIV
 G. Rubtsov, P. Satunin and S. Sibiryakov,
“Prospective constraints on Lorentz violation from ultrahigh-energy photon detection,”
  Phys. Rev. D 89 (2014) 123011
  
  [arXiv:1312.4368 [astro-ph.HE]].
  
  

Khorkhe
 N. Kalmykov, J. Cotzomi, V. Sulakov and Yu. Fomin,
“Primary cosmic ray composition at energies from 10^17 to 10^18 eV
according to the EAS MSU array data”,
Bull. Russ. Acad. Sci. Phys. 73 (2009) 547
[Izv. Ross. Akad. Nauk Ser. Fiz.  73 (2009) 584]

EAS-MSU_gamma_2013
Yu. A. Fomin et al.,


“Estimate of the fraction of primary photons in the
  cosmic-ray flux at energies ∼ 10^17 eV from the EAS-MSU
experiment data,” J. Exp.   Theor. Phys.  117 (2013) 1011

[Zh. Exp. Teor. Fiz. 144 (2013) 1153]
[arXiv:1307.4988 [astro-ph.HE]].
  
  

EAS-MSU_gamma_2014
Yu. A. Fomin et al.,


“Estimates of the cosmic gamma-ray flux at PeV to EeV
  energies from the EAS-MSU experiment data,” JETP Lett.  100
  (2015)   699 [Pisma Zh. Eksp. Teor. Fiz.  100 (2014)  797]

  [arXiv:1410.2599 [astro-ph.HE]].
  
  

EAS-MSU_Full_MC
Yu. A. Fomin et al.,


 “Full Monte-Carlo description of the Moscow State University Extensive Air Shower experiment,”
  JINST 11 (2016)  T08005

  [arXiv:1607.00309 [astro-ph.HE]].
  
  

EAS-MSU_mu
Yu. A. Fomin et al.,


“No
  muon excess in extensive air showers at ∼ 10^17 eV primary energy:
  EAS-MSU muon versus surface detector data,” arXiv:1609.05764
  [astro-ph.HE].


EAS-MSU
S. N. Vernov



et al.,
“New Installation Of Moscow State University For Studying
Extensive Air Showers With Energies To 10^18 eV,”
  Bull. Russ. Acad. Sci. Phys.  44 (1980) 80
   [Izv. Ross. Akad. Nauk Ser. Fiz.  44 (1980) 537].
  
  
book
G.B. Khristiansen, G.V. Kulikov, Yu.A. Fomin. Cosmic Rays of
Superhigh Energies, Karl Tiemig Verlag, Munchen, 1980, 246 p.  

Heck:1998vt
  D. Heck, G. Schatz, T. Thouw, J. Knapp and J. N. Capdevielle,
  “CORSIKA: A Monte Carlo code to simulate extensive air showers,”
  FZKA-6019.
  
  

Ostapchenko:2010vb
  S. Ostapchenko,
  “Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: I. QGSJET-II model,”
  Phys. Rev. D 83 (2011) 014018
  
  [arXiv:1010.1869 [hep-ph]].
  
  

Fluka
 G. Battistoni et al.,
“The FLUKA code:  Description and benchmarking”,
AIP Conf.Proc.896(2007) 31.

EGS4
W. R. Nelson, H. Hirayama, D. W. O. Rogers,
“The EGS4 code system”, Tech. Rep. 0265, SLAC (1985).

AGYak
G. I. Rubtsov, L. G. Dedenko, G. F. Fedorova  et al.,



  Phys. Rev. D 73 (2006) 063009
  [astro-ph/0601449].
  
  

FC
 G. J. Feldman and R. D. Cousins,
“A Unified approach to the classical statistical analysis of small signals,”
  Phys. Rev. D 57 (1998) 3873

  [physics/9711021 [physics.data-an]].
  
  

EAS-TOP
M. Aglietta, B. Alessandro, P. Antoni et al.  [EAS-TOP Collaboration],
 “A Limit to the rate of ultrahigh-energy gamma-rays in the primary cosmic radiation,”
  Astropart. Phys.  6 (1996) 71.
  
  

CASA-MIA
M. C. Chantell, C.E. Covault, J.W. Cronin et al.  [CASA-MIA
Collaboration],
“Limits on the isotropic diffuse flux of ultrahigh-energy gamma radiation,”
  Phys. Rev. Lett.  79 (1997) 1805
  [astro-ph/9705246].
  
  

KASCADE
 Z. Feng et al. [KASCADE-Grande Collaboration],
“Limits on the isotropic diffuse γ-rays at ultra high energies measured with KASCADE,”
  PoS ICRC 2015 (2016) 823.
  
  

KASCADE-Grande
 D. Kang et al. [KASCADE-Grande collaboration],
  “A limit on the diffuse gamma-rays measured with KASCADE-Grande,”
  J. Phys. Conf. Ser.  632 (2015)   012013.
  
  
  

Yak2
A. V. Glushkov, I. T. Makarov, M. I. Pravdin  et al.,

  “Constraints on the flux of primary cosmic-ray photons at energies E >
10^18 eV from Yakutsk muon data,” Phys. Rev. D 82 (2010) 041101
  [arXiv:0907.0374 [astro-ph.HE]].

  

PAO1
A. Aab et al. [Pierre Auger Collaboration],
“The Pierre Auger Observatory: Contributions to the 34th International Cosmic Ray Conference (ICRC 2015),”
  arXiv:1509.03732 [astro-ph.HE].
  
  

PAO2
P. Abreu, M. Aglietta, E.J. Ahn  et al.
[Pierre Auger Collaboration],
“The Pierre Auger Observatory III: Other Astrophysical Observations,”
  arXiv:1107.4805 [astro-ph.HE].
  
  

AGASA
K. Shinozaki, M. Chikawa, M. Fukushima et al. [AGASA Collaboration],

“Upper limit on gamma-ray flux above 10^19 eV estimated by the Akeno
  Giant Air Shower Array experiment,” Astrophys. J.  571 (2002)
  L117.
  
  

TA
 G. I. Rubtsov et al. [Telescope Array Collaboration],
“Telescope Array search for photons and neutrinos with the surface detector data,”
  PoS ICRC 2015 (2016) 331.
  
  






§ APPENDIX. DISTRIBUTION OF EAS PARAMETERS FOR DATA AND MONTE-CARLO


  In this appendix, we provide the data to Monte-Carlo comparison
  of the EAS parameter distribution for the extended cut
  N_e>10^7 used is the present work. The primary composition is
  the same as determined in Ref. <cit.> by fitting the
  observed distribution of the muon densities: 46% protons and
  54% iron. The distributions of S, N_e, R, θ and
  muon density at 100 meters ρ_μ(100) are shown in
  Figures <ref>-<ref>.  The validity of the N_e
  cut extension is verified by reasonable agreement of the data and
  Monte-Carlo.











