Table Of ContentUpper limit on the ultra-high-energy photon flux from AGASA and Yakutsk data
G.I. Rubtsov1, L.G. Dedenko2,3, G.F. Fedorova3, E.Yu. Fedunin3, A.V. Glushkov4,
D.S. Gorbunov1, I.T. Makarov4, M.I. Pravdin4, T.M. Roganova3, I.E. Sleptsov4 and S.V. Troitsky1
1Institute for Nuclear Research of the Russian Academy of Sciences, Moscow 117312, Russia
2Faculty of Physics, M.V. Lomonosov Moscow State University, Moscow 119992, Russia
3D.V. Skobeltsin Institute of Nuclear Physics, M.V. Lomonosov Moscow State University, Moscow 119992, Russia and
4Yu.G. Shafer Institute of Cosmophysical Research and Aeronomy, Yakutsk 677980, Russia
(Dated: January 13, 2006)
Wepresenttheinterpretationofthemuonandscintillationsignalsofultra-high-energyairshowers
observedbyAGASAandYakutskextensiveairshowerarrayexperiments. Weconsidercase-by-case
ten highest energy eventswith known muon content and conclude that at the95% confidencelevel
(C.L.) noneof them was inducedby aprimary photon. Takinginto account statistical fluctuations
and differences in the energy estimation of proton and photon primaries, we derive an upper limit
6
of 36% at 95% C.L. on the fraction of primary photons in the cosmic-ray flux above 1020 eV. This
0
result disfavors the Z-burst and superheavy dark-mattersolutions to theGZK-cutoff problem.
0
2
PACSnumbers: 98.70.Sa,96.40.De,96.40.Pq
n
a
J I. INTRODUCTION lows) [4] with reconstructed energies above 8×1019 eV
9 and measured muon content. We reject the hypothe-
1 sis that any of showers considered was initiated by a
One of the most intriguing puzzles in astroparticle
1 physics is the observation of air showers initiated by photon primary at the 95% confidence level (C.L.). We
then derive as our main result an upper limit of 36%
v particles with energies beyond the cutoff predicted by
(at 95% C.L.) on the fraction ǫ of primary photons
9 Greisen and by Zatsepin and Kuzmin [1]. Compared γ
4 to lower energies, the energy losses of protons increase with original energies above 1020 eV (the difference be-
4 sharplyat≈5×1019eVsincepionproductiononcosmic tween original and reconstructed energies is discussed in
1 Sec. II).
microwavebackgroundphotonsreducestheprotonmean
0 Therestofthepaperisorganizedasfollows. InSec.II
free path by more than two orders of magnitude. This
6
we discuss the experimental data set which we use for
0 effect is even stronger for heavier nuclei, while photons
our study. In Sec. III, the details of the simulation of
/ are absorbed due to pair production on the radio back-
h the artificial showerlibraries and comparisonof the sim-
ground with the mean free path of a few Mpc. Thus,
p ulatedandrealdata aregiven. This sectioncontains the
the cosmic-ray (CR) energy spectrum should dramati-
-
o cally steepen at ≈7×1019 eV for any homogeneous dis- description of our method and the main results. We dis-
r cusshowrobusttheseresultsarewithrespecttochanges
tribution of CR sources. Despite the contradictions in
t
s the shape of the spectrum, the existence of air showers in assumptions, to analysis procedure, and to variations
:a with energies in excess of 1020 eV is firmly established intheexperimentaldata,inSec.IV. InSec.V,wediscuss
v the differences between our approachand previous stud-
by severalindependent experiments using different tech-
i ies, which allowed us to put a significantly more strin-
X niques (Volcano Ranch [2], Fly’s Eye [3], Yakutsk [4],
gent limit on the gamma-ray fraction. Our conclusions
AGASA[5],HiRes[6]andPierreAuger[7]experiments).
r
are briefly summarized in Sec. VI.
a Some explanations for these showers,like the Z-burst or
top-down models, predict a significant fraction of pho-
tons above typically 8×1019 eV (for reviews see, e.g.,
Refs. [8]). Indications for the presence of neutral par- II. EXPERIMENTAL DATA
ticles at lower energies were found in Refs. [9]. Thus,
the determination of the photon fraction in the CR flux AGASAwasoperatingfrom1990to2003andconsisted
is of crucial importance, and the aim of this work is to of 111 surface scintillation detectors (coveringan area of
derive a stringent limit on this fraction in the integral about100km2)and27muondetectors. Theareasofthe
CR flux above 1020 eV. To this end, we compare the re- AGASA muon detectors varied between 2.8 and 20 m2.
ported information on signals measured by scintillation The detectors consisted of 14–20 proportional counters
and by muon detectors for observed showers with those aligned under a shield of either 30 cm of iron or 1 m of
expectedbyairshowersimulations. Wefocusonthesur- concrete and were placed below or close to scintillation
facedetectorsignaldensityat600metersS(600)(known detectors. The threshold energy was 0.5 GeV/cosθ for
µ
aschargedparticledensity)andthemuondensityat1000 muons with zenith angle θ [11]. During 14 years of op-
µ
m, ρµ(1000), which are used in experiments as primary eration, AGASA had observed 11 events with reported
energy and primary mass estimators, respectively. energiesabove1020 eV andzenith anglesθ <45◦ [5, 12].
We study individual events of AGASA [10] and of the Among them, six events had ρµ(1000) determined [11].
Yakutskextensiveairshowerarray(Yakutskinwhatfol- Yakutsk is observing CRs of highest energies since
2
1973, with detectors in various configurations. With 5×1019 eVand5×1020 eVtotakeintoaccountpossible
θ < 60◦, it has observed three events above 1020 eV, deviations of E from E . The arrival directions of the
est 0
all with measured muon content. Before 1978, only one simulated showers were the same as those of the corre-
muon detector with the area of 8 m2 and threshold en- sponding real events. The simulations were performed
ergy0.7GeV/cosθ wasinoperation. Later,ithasbeen with CORSIKA v6.204 [18], choosing QGSJET 01c [19]
µ
replacedby sixdetectors withareasup to 36m2 andthe as high-energy and FLUKA 2003.1b [20] as low-energy
threshold energy of 1.0 GeV/cosθ [13]. hadronic interaction model. Electromagnetic shower-
µ
In our study, we combine the AGASA and Yakutsk ing was implemented with EGS4 [21] incorporated into
datasets,motivatedbythefollowing. First,bothdatasets CORSIKA. Possible interactions of the primary pho-
are obtained from surface array experiments operated tons with the geomagnetic field were simulated with the
with similar plastic scintillation detectors. Second, the PRESHOWER option of CORSIKA [22]. As discussed
energy estimation procedures of the two experiments inSec. IVB,this choiceof the interactionmodels results
are compatible, within the reported systematic errors at in a conservativelimit on gamma-rayprimaries. As sug-
∼ 1020 eV, if differences in the observational conditions gested in Ref. [23], all simulations were performed with
are taken into account [14]. Finally, the values of the thinninglevel10−5,maximalweight106forelectronsand
CR flux at 1020 eV reported by the two experiments are photons, and 104 for hadrons.
consistent within their 1σ errors. For eachsimulated shower,we determined S(600)and
The shower energy estimated by an experiment (here- ρ (1000). For the calculation of S(600), we used the de-
µ
after denoted as Eest) is in general different from the tectorresponsefunctions fromRefs.[24, 25]. Foragiven
true primary energy (denoted as E0) because of natu- arrival direction, there is one-to-one correspondence be-
ral shower fluctuations, etc. Moreover, the energy es- tween S(600) and the quantity called estimated energy,
timation algorithms used by surface-array experiments E . Therelationisdeterminedbythestandardanalysis
est
normally assume that the primary is a proton. While procedureofthetwoexperiments[10,27]. Thisallowsus
the estimated energy for nuclei depends only weakly on toselectsimulatedshowerscompatiblewiththeobserved
their mass number, the difference between photons and onesbythe signaldensity. ThequantityS(600)isrecon-
hadrons is significant. For photons, the effects of geo- structed not precisely. In terms of estimated energy, for
magneticfield[15]resultindirectionaldependenceofthe AGASA events, the reconstructed energies are are dis-
energy reconstruction. Thus, the event energy reported tributedwithaGaussianinlog E /E¯ ;thestandard
(cid:0) est rec(cid:1)
by the experiment should be treated with care when we deviation of E is σ ≈ 25% [14]. For Yakutsk events,
est
allow the primary to be a photon. In this study we in- the correspondingσ has beendeterminedevent-by-event
clude events with Eest ≥ 8 × 1019 eV because of pos- and is typically 30–45% [28]. To each simulated shower,
sible energy underestimation for photon-induced show- we assigned a weight w proportional to this Gaussian
1
ers; these events contribute to the final limit, derivedfor probability distribution in logE centered at the ob-
est
E0 >1020 eV, with different weights. served energy E¯rec =Eobs. Additionally, each simulated
For AGASA, we use the events given in Ref. [5] that shower was weighted with w to reproduce the thrown
2
pass the “cut B” defined in Ref. [11], that is having energy spectrum ∝ E−2 which is typically predicted by
0
at least one [46] muon detector hit between 800 m and non-accelerationscenarios(see Sec. IVC for a discussion
1600 m from the shower axis. The ρµ(1000) of the in- of the variations of the spectral index). For each of the
dividual events can be read off from Fig. 2 of Ref. [11]. ten observed events, we obtained a distribution of muon
Yakutsk muon detectors have larger area and are more densities ρ (1000) representing photon-induced showers
µ
sensitive both to weak signals far from the core and to compatiblewiththe observedonesby S(600)andarrival
strongsignalsforwhichAGASAdetectorsmightbecome directions. To this end, we calculated ρ (1000) for each
µ
saturated. This allowedthe Yakutsk collaborationto re- simulated shower by making use of the same muon lat-
lax the cuts, as compared to AGASA, and to obtain re- eral distribution function as used in the analysis of real
liable values of ρµ(1000) using detectors between 400 m data [11, 13]. To take into account possible experimen-
and 2000 m from the shower axis [16, 17]. Providing tal errors in the determination of the muon density, we
these cuts, six AGASA and four Yakutsk events entered replacedeachsimulatedρ (1000)byadistributionrepre-
µ
the dataset in our study (see Table I for the event de- sentingpossiblestatisticalerrors(50%and25%Gaussian
tails). for AGASA cut B [29] and Yakutsk [17], respectively).
The distribution of the simulated muon densities is the
sum of these Gaussians weighted by w w .
1 2
III. SIMULATIONS AND RESULTS Atypicaldistributionofsimulatedρ (1000)isgivenin
µ
Fig.1,forgamma-andproton-inducedsimulatedshowers
In order to interpret the data, for each of the ten compatiblewiththeevent3byS(600)andthearrivaldi-
events, we generated a shower library containing 1000 rection. We will see below that this particular event has
showers induced by primary photons [47]. Thrown ener- the largest probability of gamma interpretation among
gies E of the simulated showerswere randomly selected all ten events in the data set; still the proton interpreta-
0
(see below the discussion of the initial spectra) between tion looks perfect for it. This is the case for all events
3
TABLEI: Descriptionoftheindividualeventsusedinthiswork. Columns: (1),eventnumber;(2),experiment;(3),dateofthe
event detection (in the format dd.mm.yyyy);(4), the reported energy assuming a hadronic primary (in units of 1020 eV); (5),
thezenithangle(indegrees);(6)theazimuthangle(indegrees,φ=0correspondstoaparticlecomingfromtheSouth,φ=90◦
–fromtheWest);(7)numberofmuondetectorsusedtoreconstructmuondensity;(8)muondensityat1000mfromtheshower
axis (in unitsof m−2); (9), probability that this eventwas initiated by a photon with E >1020 eV; (10), probability that this
event was initiated by a non-photon with E >1020 eV, assuming correct energy determination. The sum p(i)+p(i) gives the
1 2
weight of thisevent in thefinal limit on ǫγ. The probability that theprimary had the energy E <1020 eV is 1−p(1i)−p(2i).
i Experiment Date Eobs θ φ ndet ρ(µi)(1000) p(1i) p(2i)
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
1 AGASA 10.05.2001 2.46 36.5 79.2 3 8.9 0.000 1.000
2 AGASA 03.12.1993 2.13 22.9 55.5 1 10.7 0.001 0.998
3 AGASA 11.01.1996 1.44 14.2 27.5 >1 8.7 0.013 0.921
4 AGASA 06.07.1994 1.34 35.1 234.9 1 5.9 0.003 0.887
5 AGASA 22.10.1996 1.05 33.7 291.6 >1 12.6 0.000 0.581
6 AGASA 22.09.1999 1.04 35.6 100.0 >1 9.3 0.000 0.565
7 Yakutsk 18.02.2004 1.60 47.7 180.8 5 19.6 0.000 0.876
8 Yakutsk 07.05.1989 1.50 58.7 230.6 5 11.8 0.000 0.868
9 Yakutsk 21.12.1977 1.10 46.1 346.8 1 8.0 0.000 0.645
10 Yakutsk 02.05.1992 0.85 55.7 163.0 5 4.7 0.000 0.303
0.2 Let us consider the ith observed event. Denote by
M the weighted number of showers contributed to the
ρ (1000) distribution for the simulated photon-induced
µ
showerscompatiblewiththeitheventbyarrivaldirection
and S(600) (throughout this paragraph, the weighted
0.1 number is the sum of corresponding weights, that is
M is the sum of weights of all 1000 showers simu-
lated for the ith event). Some of the simulated show-
ers contributed to the part of the distribution for which
ρ (1000) > ρ(i)(1000), where ρ(i)(1000) is the observed
µ µ µ
valueforthisevent. Theweightednumberoftheseshow-
5 10 15 20
′ ′
ersisM . Somepartl ofthis M correspondstoshowers
Log ΡΜH1000L with E0 > 1020 eV, the rest (M′−l) to E0 < 1020 eV.
The probability p(i) of case (i) is p(i) = l/M, while the
1 1
FIG.1: Weighteddistributionsofmuondensityρµ(1000)for probability that the event is consistent with a photon
the simulated events compatible with the event 3 by S(600) of E < 1020 eV is p′(i) = (M′ − l)/M. Moreover,
and the arrival direction. Units in the vertical axis are ar- 0 1
cboitrrraersyp,onρdµs(1t0o00p)riimsamryeapsuhroetdonisn; imt−is2.thTehdeisttrhiibnutdioanrkulsinede tphreimparroybaisbil1it−y tph1(ai)t−thpe1′(eiv)e=nt i1s−deMscr′i/bMed. byWaenyasosuthmeer
for our analysis. The thick grey line is the distribution ob- that the experimental energy estimation works well for
tained in the same way but for 500 proton-induced showers. non-photon primaries and determine the fraction ξ of
The arrow indicates the observed value of ρµ(1000) for the events with E > 1020 eV simply from the Gaussian
event 3. The disributions include 50% Gaussian error of the log(E ) distribution, so the probability of the case (ii)
est
detector. is p(i) = ξ(1−M′/M). The values of p(i) are presented
2 1,2
in Table I. Note that p(i) +p(i) < 1 because of a non-
1 2
exceptevent7,whichhastoohighρ (1000)foraproton; zero probability that a simulated shower is initiated by
µ
possible nature of its primary particle will be discussed a primary with E0 < 1020 eV. This happens especially
elsewhere. for events with reported energies close to 1020 eV and
reduces considerably the effective number of events con-
Toestimatetheallowedfractionǫ ofprimaryphotons
γ
amongCRswithE >1020 eV,wecompare,foreachob- tributing to the limit on ǫγ: since we are interested in
0 the limit for E > 1020 eV only, each event contributes
servedevent,twopossibilities: (i)thatitwasinitiatedby 0
a photonprimarywith E >1020 eVand(ii) thatitwas to the result with the weight (p(i) +p(i)). Inspection of
0 1 2
initiated by any other primary with E > 1020 eV for Table I demonstrates that the total effective number of
which the experimental energy estimat0ion works prop- events with E0 > 1020 eV (the sum of p1(i) and p(2i) over
erly. all ten events) is 7.67.
4
If the ith primary particle was a photon with E >
0
Z-burst
1020 eV with the probability p(i) and a non-photon with
1 00..88
E > 1020 eV with the probability p(i), one can eas- A RH
0 2
ily calculate the probability P(n1,n2) to have n1 pho- 00..66 SHDM
tons and n non-photons in the set of N = 10 observed HP HP
2 ΓΓ
events (0 ≤ n1 +n2 ≤ N, the rest N −n1 −n2 events ΕΕ 00..44 A AY TD
have E < 1020 eV). From the set of N events, one A
0
should take all possible non-overlapping subsets of n1 00..22 PA
and n events and sum up probabilities of these real-
2
isations (since p(i) 6= p(j), these probabilities are dif- 00
1,2 1,2
ferent for different realisations with the same n1 and 1199 1199..55 2200 2200..55
n ). Now, suppose that the fraction of the primary pho-
2 LLoogg EE
tons at E > 1020 eV is ǫ . Then, the probability to
0 γ
have n photons and n non-photons at E > 1020 eV
1 2 0
is ǫn1(1−ǫ )n2, and the probability that the observed FIG. 2: Limits (95% C.L.) on the fraction ǫγ of photons in
γ γ the integral CR flux versus energy. The result of the present
muon densities were obtained with a given ǫ is
γ work (AY) is shown together with limits previously given in
Refs. [32](HP),[11](A),[30](RH)and[33](PA).Alsoshown
N
are predictions for the superheavy dark-matter (thick line)
P(ǫγ)= X ǫγn1(1−ǫγ)n2P(n1,n2) andtopological-defect(necklaces,betweendottedlines)mod-
n1,n2=0 els [34] and for theZ-burst model (shaded area) [35].
(cf.Ref.[30]foraparticularcasen +n =N; notethat
1 2
the combinatorial factor is included in the definition of
IV. ROBUSTNESS OF THE RESULTS
P(n ,n )). The cases n +n <N reflectthe possibility
1 2 1 2
that some of the N events correspond to primaries with
E < 1020 eV. In our case, the probability P(ǫ ) is a In this section, we discuss systematic uncertainties of
0 γ
monotonically decreasing function of ǫ . Thus the up- our limit that are related to the air shower simulations,
γ
per limit on ǫ at the confidence level α′ is obtained by to the data interpretation and to selection cuts.
γ
solvingtheequationP(ǫ )=1−α′. Forourdataset,the
γ
95%C.L.upperlimitonthephotonfractionisǫ <0.33.
γ
The limit on ǫγ is rather weak compared to the individ- A. Systematic uncertainty in the S(600) and energy
ualvaluesofp(i) becauseofthesmallnumberofobserved determination
1
events.
However,someofthe photon-inducedshowersmayes- The systematic uncertainty in the absolute energy
capefromourstudybecausetheymaynotpassthemuon determination is 18% and 30% for AGASA [10] and
measurementqualitycutsortheirestimatedenergyisbe- Yakutsk [4], respectively. These systematic errors origi-
low 8×1019 eV. Possiblereasonsfor an underestimation nate from two quite different sources: (a) the measure-
of the energy may be either the LPM effect [31] or sub- ment of S(600) and (b) the relationbetween S(600)and
stantial attenuation of gamma-induced showers at large (i)
primary energy. The probabilities p that a particu-
1
zenith angles. To estimate the fraction of these “lost”
lar event may allow for a gamma-ray interpretation are
events, we have simulated 1000 gamma-induced showers
not at all sensitive to the S(600)-to-energy conversion
for each experiment with arrival directions distributed
because we select simulated events by S(600) and not
according to the experimental acceptance. We find that
by energy. These probabilities may be affected by rela-
the fraction of the “lost” events is ∼ 3.5% for AGASA
tivesystematicsindeterminationofρ (1000)andS(600).
µ
and ∼ 15% for Yakutsk. The account of these fractions, Onthe other hand, in the calculationof p(i) we assumed
weightedwiththerelativeexposuresofbothexperiments, 2
thattheexperimentalenergydeterminationiscorrectfor
results in the final upper limit,
non-photon primaries; the values of p(i) and the effec-
2
tive number of events contributing to the limit on ǫ at
ǫ <36% (95% C.L.). γ
γ
E >1020 eVwouldchangeiftheenergiesaresystemati-
0
In Fig. 2, we present our limit on ǫ (AY) together callyshifted. Inourcase(allp(i) ≈0),thereportedvalue
γ 1
with previously published limits on the same quantity. of ǫγ would be applicable to the shifted energy range in
Also, typical theoretical predictions are shown for the that case.
superheavy dark-matter, topological-defect and Z-burst Thus, the 95% C.L. conclusion that none of the ten
models. Our limit on ǫ is currently the strongest one events considered here was initiated by a photon is ro-
γ
at E > 1020 eV. It disfavours some of the theoretical bustwithrespecttoanychangesintheS(600)-to-energy
0
models such as the Z-burst and superheavy dark-matter conversion. As for the limit on ǫ we report, instead of
γ
scenarios. E > 1020 eV, it would be applicable to a different en-
0
5
ergyrangeifallexperimentalenergiesaresystematically D. Width of the ρ distribution
µ
shifted. Oneshouldnotethattheoreticalpredictions,e.g.
the curves shown in Fig. 2, would also change because
Clearly, the rare probabilities of high values of
they are normalised to the observed AGASA spectrum.
ρ (1000) in the tail of the distribution for primary pho-
µ
tons depend on the width of this distribution. The fol-
lowing sources contribute to this width:
B. Interaction models and simulation codes
• variations of the primary energy compatible with
the observed S(600) (larger energy correspond to
Our simulations were performed entirely in the COR-
larger muon number and ρ (1000));
SIKAframework,andanychangeintheinteractionmod- µ
els or simulation codes, which affects either S(600) or
ρ (1000), may affect our limit. We have studied the • physical shower-to-shower fluctuations in muon
µ
model dependence of our results by comparing differ- density for a given energy (dominated by fluc-
ent low- and high-energy hadronic interaction models tuations in the first few interactions, including
(GHEISHA [36] versus FLUKA, SIBYLL 2.1 [37] ver- preshowering in the geomagnetic field);
sus QGSJET). Our result is quite stable with respect to
these changes. In all cases, individual values of p(i) are • artificial fluctuations in S(600) and ρµ(1000) due
1 to thinning;
always close to zero, thus the limit on ǫ is not affected.
γ
The change of the low energy model does not at all af-
• experimental errors in ρ (1000) determination.
fect the reported values. In use of SIBYLL compared µ
with QGSJET, ρ (1000) is ∼ 20% smaller for photon-
µ While the first two sources are physical and are fully
induced showers. While S(600) is almost unchanged,
controlled by the simulation code, the variations of the
events in our dataset are better explained by showers
last two may affect the results.
initiatedbyheaviernucleiandtheprobabilityofphoton-
induced showers is even smaller. A similar effect is ex-
pectedforthecominginteractionmodelQGSJETII[38].
We also performed simulations with the help of the 1. Artificial fluctuations due to thinning
hybrid code [39] which reproduced the CORSIKA re-
sults with high accuracy. Another popular simulation It has been noted in Ref. [43] that the fluctuations in
code, AIRES [40], differs from CORSIKA mainly in the ρ (1000) due to thinning may affect strongly the preci-
µ
low-energy hadronic interaction model (which is fixed sionofthe compositionstudies. Forthe thinning param-
in AIRES to be the Hillas splitting algorithm), hence eters we use, the relative size of these fluctuations is [44]
we hope that simulations with AIRES would not signif- ∼ 10% for ρ (1000) and ∼ 5% for S(600). Thus with
µ
icantly affect our results. Comparison with AIRES will more precise simulations, the distributions of muon den-
be presented elsewhere. sities should become more narrow, which would reduce
The values presented here were obtained for the stan- the probability of the gamma-ray interpretation of each
dard parameterizationof the photo-nuclearcross section of the studied events even further.
given by the Particle Data Group [41] (implemented as
default in CORSIKA). The muon content of gamma-
induced showers is in principle sensitive to the extrap-
2. Experimental errors in ρµ(1000) determination
olation of the photonuclear cross section to high ener-
gies. The hybrid code [39] allows for easy variations of
Thedistributionsofρ (1000)weuseaccountedforthe
the cross section; we checked that the results are stable µ
errorinexperimentaldeterminationofthisquantity. The
for various reasonable extrapolations, in agreement with
sizeoftheerrorswastakenfromtheoriginalexperimental
Ref. [42].
publications [17, 29]. In principle, this error depends on
the event quality and on the muon number itself, which
is lower for simulated gamma-induced showers than for
C. Primary energy spectrum
the observed ones. However, e.g. Ref. [11] states that
for the AGASA cut A (two or more muon detectors),
For our limit, we used the primary photon spectrum the error is 40%, lower than 50% we use [29]. Note that
E−αforα=2. Whiletheindividualprobabilitiesp(i) are Ref.[11]discussesmuondensitiesaslowas0.04m−2 and
0 1,2
notaffectedbythechangeofthespectralindexαbecause even 0 m−2, much lower than ∼ 1 m−2 typical for our
the simulated events are selected by S(600) anyway, the simulated gamma-induced events. Still, we tested the
value of α changes the fraction of “lost” photons and, stability of our limit by taking artificially high values of
correspondingly, the final limit on ǫ . Variations of 1 ≤ experimental errors in muon density: 100% for AGASA
γ
α ≤ 3 result in the photon fraction limits between 36% and 50% for Yakutsk. The limit on ǫ changes to 37%
γ
and 37% (95% C.L.). (95% C.L.) in that case.
6
E. Data selection cuts 21
Since all events in the data set are unlikely to be ini- 20.75
tiated by primary photons (all p(i) ≈ 0), the limit on
1 20.5
ǫ is determined by statistics only and is affected if the
γ
number of events is changed. Here, we discuss possible c20.25
e
variationsofthedatasetcorrespondingtomorestringent Er
quality cuts which reduce the event number and weaken g 20
o
the limit. L
19.75
19.5
1. Zenith angle
19.25
All Yakutsk events in the data set have zenith angles
◦ ◦ ◦ 20.1 20.2 20.3 20.4 20.5 20.6
45 < θ < 60 , so the cut θ < 45 imposed by AGASA
reducesthesampletosixAGASAeventswhichresultsin LogE0
thelimitǫ <50%(95%C.L.). Oneshouldnotehowever
γ
thatAGASAmuondetectorsarenotsensitivetoinclined FIG. 3: Direction dependence of the reconstructed energy
showers, which is not the case for Yakutsk. forgamma-rayprimaries. Plottedisthereconstructedenergy
(determinedbytheAGASAmethod from S(600)) versusthe
primary energy. Dark boxes: arrival direction of theevent 1;
2. Core inside array crosses: arrival direction of the event 3; grey circles: arrival
directionsrandomlydistributedaccordingtotheAGASAac-
ceptance (0 < θ < 45◦). Straight line represents E = E .
rec 0
AnothercutimposedontheAGASApublisheddataset Both E and E are measured in eV.
0 rec
isthelocationofthecoreinsidearray. Theeventnumber
7 does not satisfy this criterion; its exclusion from the
data set results in ǫ <40% (95% C.L.).
γ In Refs. [11, 16], no conclusion was derived about ǫ
γ
at E > 1020 eV, and the data points corresponding to
highest-energyeventswerefoundtobe quite closeto the
3. More than one muon detector gamma-ray domain. To our opinion, the main source of
thiseffectisaveragingoverarrivaldirectionswhichintro-
Reconstruction of the muon density at 1000 m from duced additional fluctuations for gamma-ray primaries
a single muon detector reading requires extrapolation of due to direction-dependent preshowering (see Fig. 3 for
the lateral distribution function with an averaged slope. an illustration). In Ref. [30] which discussed the same
Though it is well-studied, the data points corresponding six AGASA events, all simulated showers for an event
toeventswithasinglemuondetectorhitmightbeconsid- with the observed energy Eobs had energies 1.2Eobs (up
eredlessreliablethanthosewithtwoormorehits. With to the energy reconstruction uncertainty of 25%). This
◦
the account of the events with two or more hits only, conversionhadbeenobtainedastheaverageoverθ <36
we are left with seven events (four AGASA and three in Ref. [11] using AIRES simulation code [40]. That is,
Yakutsk)whichweakensthe95%C.L.limittoǫ <48%. not only the average results were applied to individual
γ
showers, but effectively muon densities were simulated
with CORSIKA while energies – with AIRES, though
V. COMPARISON WITH OTHER STUDIES the two codes result in a systematically different rela-
tions between energy and S(600). Artificially high ener-
giesresultedinhigher,closertoobserved,muondensities
Some of the previous studies used the AGASA [11,
for simulated photonic showers. In our event-by event
30] and Yakutsk [16] muon data to limit the gamma-ray
simulations with CORSIKA, the energies of gamma-ray
primaries at high energies. Our results differ from the
primarieswhoseS(600)werecompatibletoobservedval-
previous ones not only because we join the data sets of
ues,werenothigherbyafactor1.2,butinfactevenlower
the two experiments. Two major distinctive features of
than E for some of the events: besides the difference
our approach allowed us to put the stringent limit: obs
in simulation codes, this is partially due to non-uniform
• both ρ (1000) and S(600) were tracked for simu- distribution of the highest-energy AGASA events on the
µ
lated showers within framework of a single simu- celestial sphere [12, 45] which makes the usage of aver-
lation code (CORSIKA in our case); aged energies poorly motivated.
The impact of two other sources of difference between
• each event was studied individually, without aver- our approach and that of Ref. [30] is less important for
aging over arrival directions. the final result: (i) Ref. [30] does not account for the
7
2 in quadrature. We see that the main source of the dis-
agreementis inthe values ofE whichpush, forthe case
0
1.5 of Ref. [30], the simulated muon densities closer to the
observed observed one.
1
L
0 0.5
0
0
1
HΜ 0 VI. CONCLUSIONS
Ρ
og -0.5 simulated RH
L To summarize, we have studied the possibility that
-1 the highest-energy events observed by the AGASA and
Yakutsk experiments were initiated by primary photons.
-1.5
Comparing the observed and simulated muon content of
these showers, we reject this possibility for each of the
20.1 20.2 20.3 20.4 20.5 20.6 teneventsatE >8·1019 eVatleastatthe95%C.L.An
LogE important ingredient in our study is the careful tracking
0
of differences between the originaland reconstructeden-
ergies. This allows us to put an upper bound of 36% at
FIG. 4: Illustration of the difference between our study and
Ref. [30]. Plotted is the muon density at 1000 m versus the 95%C.L.onthefractionǫγ ofprimaryphotonswithorig-
primaryenergy. Smallgreyboxes: simulatedgamma-induced inalenergiesE0 >1020eV,assuminganisotropicphoton
eventswitharrivaldirectionoftheevent1. Filledbox,marked fluxandE−2spectrum. Thislimitisthestrongestoneup
0
“simulated”: simulated events compatible with the event 1 to date. It strongly disfavorsthe Z-burst and constrains
byS(600). Openbox,marked“RH”:simulatedshowerswith severely superheavy dark-matter models. The method
average E0 = 1.2Eobs from Ref. [30]. The observed value of thatwehaveusedisquitegeneralandmaybeappliedat
ρµ(1000) for the event 1 is represented by a horizontal line, other energies and to other observables.
marked “observed”. E0 is measured in eV, ρµ(1000) in m−2.
We are indebted to M. Kachelrieß, K. Shinozaki and
See thetext for more details.
M. Teshima for numerous helpful discussions and col-
laboration at initial stages of this work. We thank
“lost” photons and (ii) the detector error is applied in L. Bezroukov, V. Bugaev, R. Engel, D. Heck, A. Ring-
our study to the simulated events while in Ref. [30] – to wald,M.Risse,V.Rubakov,D.SemikozandP.Tinyakov
the observed ones. forhelpful discussionsandcomments onthe manuscript.
The difference with Ref. [30] is illustrated in Fig. 4, This study was performed within the INTAS project
where ρ (1000) is plotted versus E for simulated 03-51-5112. We acknowledge also support by fellow-
µ 0
gamma-inducedshowerswith the arrivaldirectionof the ships of the Russian Science Support Foundation and of
event#1. For simulatedevents compatible with the real the Dynasty foundation (D.G. and S.T.), by the grants
eventbyS(600),theaveragepointisshowntogetherwith NS-2184.2003.2 (D.G., G.R. and S.T.), NS-1782.2003.2,
onesigmaerrorbars. Horizontalerrorbarscorrespondto RFFI 03-02-16290 (L.D., G.F., E.F. and T.R.), NS
variations in E compatible with S(600). Vertical error 748.2003.2, RFFI 03-02-17160, RFFI 05-02-17857, FASI
0
bars include variations in simulated ρ (1000) and 50% 02.452.12.7045(A.G.,I.M.,M.P.andI.S.). G.R.andS.T.
µ
detector error. The point corresponding to simulated thank the Max-Plank-Institut fu¨r Physik (Munchen),
showers with E = 1.2E from Ref. [30] has a larger where a significant part of this work was done, for
0 obs
ρ (1000). Horizontal error bars correspond to the en- warm hospitality. Computing facilities of the Depart-
µ
ergy reconstructionaccuracy. Vertical error bars include ment of Theoretical Physics, Institute for Nuclear Re-
variationsinsimulatedρ (1000)reportedinRef.[30]and search (Moscow), were used to perform the simulations
µ
40% detector error applied to the observed value, added of air showers.
[1] K. Greisen, Phys. Rev. Lett. 16 (1966) 748; G. T. Zat- [8] P. Bhattacharjee and G. Sigl, Phys. Rept. 327 (2000)
sepin and V.A. Kuzmin,JETP Lett. 4, 78 (1966). 109;M.Kachelrieß,ComptesRendusPhys.5,441(2004).
[2] J. Linsley, Phys.Rev. Lett. 10, 146 (1963); [9] D. S. Gorbunov et. al., JETP Lett. 80, 145 (2004);
[3] D.J. Bird et al.,Astrophys.J. 441, 144 (1995); R. U.Abbasi et al.,astro-ph/0507120.
[4] V.Egorovaetal.,Nucl.Phys.Proc.Suppl.136,3(2004). [10] M. Takedaet al.,Astropart. Phys. 19, 447 (2003).
[5] N. Hayashida et al., astro-ph/0008102, see also [11] K. Shinozaki et al.,Astrophys.J. 571, L117 (2002).
http://www-akeno.icrr.u-tokyo.ac.jp/AGASA/ [12] M. Takedaet al. Phys.Rev.Lett. 81, 1163 (1998).
[6] R.U. Abbasiet al., Phys.Rev.Lett. 92, 151101 (2004); [13] A. V.Glushkov et al.,JETP Lett. 71, 97 (2000).
[7] J. Matthews et al.,FERMILAB-CONF-05-276-E-TD. [14] N. Sakaki,PhD thesis, Univ.of Tokyo,2002.
8
[15] B. McBreen and C.J. Lambert, Phys. Rev. D24, 2536 [33] M.Risse[PierreAugerCollaboration],astro-ph/0507402.
(1981). [34] R. Aloisio, V. Berezinsky, M. Kachelrieß, Phys. Rev. D
[16] S.P. Knurenkoet al. astro-ph/0411683. 69, 094023 (2004).
[17] Yakutskcollaboration, to be published. [35] D.Semikoz,G.Sigl,JCAP0404,003(2004);G.Gelmini,
[18] D. Heck et al., Report FZKA-6019 (1998), O. Kalashev and D. V. Semikoz, astro-ph/0506128;
Forschungszentrum Karlsruhe Z. Fodor, S. D. Katz, A. Ringwald, JHEP 0206, 046
[19] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, (2002). A. Ringwald, T. J. Weiler and Y. Y. Y. Wong,
Nucl.Phys. Proc. Suppl.52B, 17 (1997). Phys. Rev.D 72 (2005) 043008.
[20] A. Fass´o, A. Ferrari, P. R. Sala, in Proc. of the Mon- [36] H. Fesefeldt, PITHA 85/02.
teCarlo 2000 Conference, A. Kling el. (eds.), p. 159-164 [37] R. Engel et al., Phys. Rev. D 46, 5013 (1992);
(SpringerBerlin 2001); A.Fass´o et al., ibid. p.955. R. S.Fletcher et al.,Phys.Rev. D 50, 5710 (1994).
[21] W.R.Nelson,H.Hirayama,D.W.O.Rogers,SLAC-0265 [38] S. Ostapchenko and D. Heck, Proc. 29th ICRC (Pune),
[22] P.Homola et al., Comp. Phys. Comm. 173 (2005) 71. 2005.
[23] M. Kobal et al., Astropart. Phys., 15, 259 (2001). [39] L. G. Dedenkoet al., Proc. 29th ICRC(Pune), 2005.
[24] N.Sakakiet al.,Proc. ICRC2001, 1, 333. [40] S.J. Sciutto, astro-ph/9911331.
[25] E. Fedunin,PhD thesis, Moscow StateUniv., 2004; [41] S. Eidelman et al.,Phys.Lett. B592, 1 (2004).
[26] L. G. Dedenko et al., Nucl. Phys. Proc. Suppl. 136, 12 [42] M. Risse et al.,astro-ph/0512434.
(2004). [43] D.BadagnaniandS.J.Sciutto,Proc.29thICRC(Pune),
[27] A.V.Glushkovetal.,Phys.Atom.Nucl.63,1477(2000). 2005.
[28] M.I. Pravdin, Proc. 29th ICRC(Pune), 2005. [44] G.I. Rubtsovet al., toappear.
[29] K.Shinozakietal.,Proc.27thICRC(Hamburg)1(2001) [45] D. S. Gorbunov and S. V. Troitsky, JCAP 0312, 010
346. (2003).
[30] P.Homolaetal.,Nucl.Phys.B(Proc.Suppl.)151(2006) [46] WethankK.ShinozakiforbringingamisprintinRef.[11]
116; M. Risse et al.,Phys.Rev.Lett. 95 (2005) 171102. to our attention: “more than one” was written there.
[31] L.Landau,I.Pomeranchuk,Dokl.Acad.NaukSSSR,92, [47] FortheillustrationinFig.1,500proton-inducedshowers
535,735(1953);A.Migdal,Phys.Rev.103,1811(1956). were simulated and processed in a similar way.
[32] M. Aveet al.,Phys. Rev.D 65, 063007 (2002).
