IOP PUBLISHING

JOURNAL OF PHYSICS G: NUCLEAR AND PARTICLE PHYSICS

doi:10.1088/0954-3899/39/2/025201

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201 (16pp)

Studies of the energy spectrum and composition of the
primary cosmic rays at 100–1000 TeV from the
GRAPES-3 experiment
H Tanaka 1 , S R Dugad 1 , S K Gupta 1 , A Jain 1 , P K Mohanty 1 , B S Rao 1 ,K C Ravindran 1 , K Sivaprasad 1 , S C Tonwar 1 , Y Hayashi 2 , N Ito 2 ,S Kawakami 2 , M Minamino 2 , T Nonaka 2 and A Oshima 2
1 Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400 005, India
2 Graduate School of Science, Osaka City University, Osaka 558-8585, Japan

E-mail: gupta@grapes.tifr.res.in

Received 25 September 2011
Published 12 January 2012
Online at stacks.iop.org/JPhysG/39/025201
Abstract
The composition and energy spectrum of primary cosmic rays (PCRs) are theonly observables at high energies to study the nature of sources accelerating
PCRs to �1000 TeV. These observables have been directly measured up to∼300 TeV with detectors aboard balloons and satellites. But measurements
at >1000 TeV have to be obtained indirectly from ground-based observationsof extensive air showers. However, their interpretation relies on an inadequate
knowledge of hadronic interactions at �1000 TeV. The GRAPES-3 experimentis designed to operate at �30 TeV providing a sizable overlap in energy
with direct measurements, enabling the selection of a suitable model ofhadronic interactions at ∼1000 TeV. We present salient features of GRAPES-3
including details of muon multiplicity distributions observed with a 560 m2
detector as a function of shower size from an analysis of data of 545 days.
These distributions were compared with expectations from Monte Carlosimulations, using some of the hadronic interaction generators in CORSIKA,
to deduce energy spectra of ﬁve nuclear groups in the 100–1000 TeV region.
A comparison of GRAPES-3 results with direct measurements indicates thatSIBYLL provides a good description of hadronic interactions for interpretingour data. These measurements extend energy spectra and composition of PCRs
that is consistent with extrapolation of direct measurements.
(Some ﬁgures may appear in colour only in the online journal)

1. Introduction
The energy spectrum of the primary cosmic rays (PCRs) is so remarkable that it can be
described by a simple power law over nearly ten decades in energy, with only two changes
0954-3899/12/025201+16$33.00 © 2012 IOP Publishing Ltd

Printed in the UK & the USA

1

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et alin the exponent at ∼3×1015 the knee and ∼3×1018 eV the ankle. Despite recent progress, anunderstanding of the nature of sources and acceleration mechanism(s) responsible for PCRs
still eludes us. In the absence of direct detection of PCR sources, owing to their deﬂection inintervening the galactic magnetic ﬁeld, precise measurements of energy spectra of different
nuclear components over a wide energy range could be helpful in deciphering the process of
their acceleration and propagation. At energies � 300 TeV, direct measurements from satellite
and balloon platforms have provided crucial information on the composition and energy spectraof PCRs [1–4]. However, results obtained from direct measurements above 100 TeV sufferfrom limited statistics, due to relatively small exposures of �103 m2 h.
The indirect measurement on some of the components of an extensive air shower (EAS)such as electrons, muons, hadrons and Cherenkov photons from ground-based detectors is the
only method available for study of PCRs, wherein exposure factors of >107 m2 h are achieved[5–12]. However, with the exception of Tibet-ASγ and ARGO-YBJ experiments [7, 8], theenergy thresholds for most EAS arrays are signiﬁcantly higher than the highest energy PCRs
detected in direct experiments. This primarily is due to the large atmospheric overburden
at observational sites, and also owing to practical constraints, such as large inter-detector
separation in EAS arrays that contribute to higher threshold energies. Equally signiﬁcant, theinterpretation of the results from EAS observations are strongly inﬂuenced by the details of
hadronic interaction generators used in simulations of EAS.
Due to the absence of data from accelerator experiments at �1000 TeV in the laboratory
frame, especially in the forward region relevant for propagation of an EAS in atmosphere, thedetails of assumptions made on the energy dependence of various interaction parameters havea major impact on the ﬁnal outcome. In the past, different characteristics of an EAS, such as
shower size (Ne ), e-μ or e-h correlations, Cherenkov photon density, etc, were used to obtainthe composition and energy spectrum at �1000 TeV. Unfortunately, energy dependence of
these components for different hadronic generators is often signiﬁcantly different. Therefore,
it is necessary to validate a combination of these observables along with a hadronic generator,by having a sizable overlap in energy with direct measurements. This objective was difﬁcult
to achieve for most of the EAS experiments, due to a variety of practical constraints, some ofwhich have been outlined above.
Here, we report GRAPES-3 observations of muon multiplicity distributions (MMDs)[13–15] for several shower size ‘Ne ’ groups that were used to obtain the energy spectra ofvarious nuclear groups in PCRs in the range ∼100–1000 TeV. Therefore, a sizable overlap
in energy ∼100–300 TeV of GRAPES-3 data with direct measurements could be achievedthat helped in the selection of an appropriate hadronic generator. The well-known correlation
between Ne and the muon content Nμ as a sensitive probe of PCR composition has been usedin this work.
Simulations have shown that the shape of the MMD measured with a large area muon
detector has excellent sensitivity to the composition of PCRs [10]. This is because lighter
nuclei (H, He etc) predominantly contribute to lower multiplicities of MMD, as opposedto heavier nuclei (Al, Fe, etc) that dominate high multiplicities. This may be explained byconsidering two PCRs, an ‘H’ and an ‘Fe’, each having the same energy E0 . The Fe nuclei
would break up in the ﬁrst interaction producing A nucleons of energy ∼E0 /A each, A is themass number of Fe. In subsequent interactions, ‘A’ nucleons would produce a larger numberof pions of systematically lower energy than those produced by H of same energy. The pions
in Fe initiated EAS would more likely decay, producing a larger number of muons than in
an H initiated EAS. These intuitive arguments are also supported by detailed simulations asshown subsequently.

2

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

Figure 1. The GRAPES-3 experimental system of 257 scintillator detectors ( ) and 16 muon
detector modules ( ). Distances along the X- and Y-axes are in meters.

Salient features of the GRAPES-3 array are presented in section 2, details of the EAS trigger to optimize selection of lower energy EAS and muon detection efﬁciency are described in
section 3 including the data analysis to obtain the EAS direction, size Ne , core location and
MMDs for different-size group. In section 4, details of simulations are presented along with Ne
and muon size distributions for three hadronic generators. Results obtained on energy spectraof H and other nuclear groups are compared with direct measurements in section 5. This isfollowed by a discussion on the suitability of different generators used in the energy range
100–1000 TeV in section 6. Conclusions are presented in the ﬁnal section.

2. The GRAPES-3 experimental setup
The GRAPES-3 (Gamma Ray Astronomy at PeV Energies Phase-3) experiment consists of adensely packed EAS array at Ooty in India (11.4◦ N, 76.7◦ E, 2200 m). A schematic of the arrayof 257 scintillator detectors, each 1 m2 and deployed on a hexagonal grid with the separationof 8 m [13], is shown in ﬁgure 1. Each detector is shown as a ﬁlled triangle. The hexagonal
geometry in deployment of detectors ensured uniform selection of EAS, over the array.
The array also contains 16 modules of the tracking muon detector of total area 560 m2
shown as ﬁlled squares in ﬁgure 1. A cluster of four modules is housed inside a common hallforming a supermodule. The hall contains signal processing and recording system for muons.
Each 35 m2 module consists of 232 proportional counters (PCs) deployed in four layers of
58 each, placed in orthogonal directions. This permits a reconstruction of muon track to an
accuracy of 6◦ in two orthogonal vertical planes. Energy threshold of muons was sec(θ ) GeV,
where θ is the muon angle relative to the vertical direction. An arrangement of concrete layers
in the shape of an inverted pyramid provides the shielding up to 45◦ [14].
3

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

Figure 2. Muon time offset T0 in μs.

Figure 3. The number of reconstructed tracks as a function of incident muons from MC.

3. The trigger and data analysis
The 560 m2 GRAPES-3 tracking muon detector was used to identify multiple muon tracks
in each EAS. The muons due to accidental coincidences during the latency time of PCs, notassociated with EAS, were rejected by recording the time offset ‘T0 ’ of each muon relative tothe EAS trigger. The distribution of T0 for detected muons is shown in ﬁgure 2.
The peak in ﬁgure 2 is due to EAS muons, and the tail on both sides is due to thechance coincidence of unassociated muons. By limiting T0 within −0.8–2.6 μs, the numberof muons by chance was reduced to �0.2 per EAS. Since EAS muons are aligned along the
EAS direction, a further reduction was achieved by reconstructing the muon direction in two
orthogonal vertical planes (X−Z and Y −Z, where Z is vertical) and by requiring tracks to bewithin 10◦ of the projected EAS direction. A number of muons were measured for each plane
and the larger of the two was treated as the true muon number for that module. Total muons
were obtained by summing all 16 modules. The twin criteria of restricting time offset anddirection ensured a clean identiﬁcation of EAS muons. This estimate is reliable for a small
number of tracks, however, as the muon number increases, the tracks begin to overlap and their
number is underestimated. This saturation effect in a 35 m2 module was investigated throughsimulations and the results are displayed in ﬁgure 3.
For example, on average, ten incident muons actually get reconstructed as eight tracks.
For 16 modules, this effect becomes signiﬁcant for muon multiplicities of �120. But
it is automatically accounted for, since the reconstruction capability of the detector wasincorporated in EAS simulations. In analysis of real data, EAS cores landing between 60 and4

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

80 m from the center of the muon detector were selected as shown in ﬁgure 1. This was not to
reduce the saturation effect, but to minimize the probability of hadrons penetrating the muon
detector.
The GRAPES-3 trigger is optimized to select low energy EAS, by requiring a 3-foldcoincidence of detectors on three adjacent lines, formed by the inner 153 detectors locatedinside the dashed hexagon shown in ﬁgure 1. Also a minimum of �10 detectors must trigger,each with a signal of �30% of a minimum ionizing particle. The arrival time and charge in thephotomultiplier (PMT) pulse was recorded for all detectors, to obtain direction (θ , φ) and size
(Ne ) [13]. CORSIKA simulations of the EAS [16] for ﬁve nuclei, H, He, N, Al and Fe, had
shown a triggering efﬁciency of ∼90% at 50 TeV for H that increases with mass, reaching 100
TeV for Fe [13]. For each trigger, the following information was recorded for each PC in themuon detector: (i) time of the digital pulse with a resolution of 167 ns, (ii) pulse width with aresolution of 333 ns. Due to the exponential shape of the ampliﬁed pulse, its width provides ameasure of logarithm of the pulse amplitude.
Earlier simulations had shown that MMD is a sensitive probe of the relative mix ofnuclei in PCRs [10]. The energy spectrum of ﬁve nuclei is extracted through a comparisonof observed MMDs with simulation results for different generators of hadronic interactions.
Lateral distribution is used to ﬁt the EAS core, size Ne , age, etc. Errors in the core location,
Ne and age were studied as a function of Ne . The dependence of the differential size spectrumwith the zenith angle was also studied and used to derive the vertical ﬂux.
6×108 EAS collected in a lifetime of 4.7×107 s (545 days), spread over 2 years, wereanalyzed to measure EAS parameters, including direction (θ , φ), size Ne , core xc ,yc , age s.
The zenith ‘θ ’ and azimuth ‘φ’ were measured from timing data of each triggered detector.
The angular resolution of the array is known using three methods, even–odd, left–right andshadow of the moon on the PCR ﬂux. The angular resolution is 0.7◦ at Ne � 104 [18]. Here,EAS of θ < 25◦ and cores within 30 m from the array center as shown in ﬁgure 1 were used
for analysis.
The EAS direction θ, φ was calculated by a χ 2 minimization of the relative arrival
times of particles. Ne was obtained by ﬁtting the Nishimura–Kamata–Greisen (NKG) lateral
distribution function to observed densities [17]. Age (s) was calculated by the maximum
likelihood algorithm MINUIT [26] as described below. If ni particles pass detector i, locatedat (xi , yi ), of area Si , the expected density ρi was calculated as follows:�
ri = (xi − xc )2 + (yi − yc )2 − D2
(1)
D = (xi − xc ) sin θ cos φ + (yi − yc ) sin θ sin φ�� �s−2.0 �
ri s−4.5ri1+
ρi = C
rm
rm

(2)

(3)

(4.5 − s)Ne(4)
2πrm2 (s) (4.5 − 2s)
where rm is the Molière unit and rm = 103 m for atmospheric depth at Ooty. For incident
density ρi , the probability of detecting ni particles is pi , and then the probability of suchdetection is L, which may be derived from the following equations:(ρi Si cos θ )ni −ρi Si cos θ
e
pi =
(5)
ni !
C=

L=

�

pi .

(6)

i

5

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

1.8
1.6

σr (m)

1.4
1.2
1.0
0.8
0.6
0.4
4.2

4.4

4.6

4.8
LOG10(Ne )

5

5.2

5.4

Figure 4. Mean error σr in the core location as a function of size Ne .

An increase in the path length for the zenith angle θ was included when estimating ni . The
unknown parameters were determined by maximizing the value of L. Since the detector
response becomes nonlinear if the detected density ρi Si sec(θ ) is �200, such detectors were
excluded in the NKG ﬁt. The detector i was rejected, if probability pi was � 10−4 of the
expected value. MINUIT [26] was used for a successive approximation procedure to maximize
L. The initial core position was estimated from the weighted mean of locations of sevendetectors with the highest particle densities, labeled ‘High7’ where weights equaled density inthe detectors. Approximations used in equation (8) were originally obtained from simulations:�High7 ni ri(7)
Rcore = �High7 ni
Ne = 102(ni )0.97(8)

s = 1.0.

(9)

The core location accuracy was estimated by simulating proton EAS with the samereconstruction procedure as the real data. Next, the distance ‘ r’ between the original andreconstructed cores was calculated. The mean of r termed mean error σr is shown with Ne
in ﬁgure 4. The mean error σr improves with size, from ∼1.6 m at 104.3 to ∼0.6 m at 105.2 .
Similarly size accuracy was obtained by comparing Nﬁt from NKG with actual Ntrue . For size
range 104.4 � N < 105.2 , rms errors were �10%. The mean of ratio Nﬁt /Ntrue is shown inﬁgure 5 for 10 size groups, for the range 104.3 –105.3 in steps of 100.1 . Clearly, size is measured
to ∼8% at Ne = 104.4 , and it improves to ∼5% at 105.2 . However, the ﬁtted size was always
smaller than the true size when averaged over the size bin. Since the size for the measured andsimulated EAS was being underestimated in exactly the same manner, it would not inﬂuence
the ﬁnal outcome on composition. A similar analysis for age ‘s’ showed that an accuracy of∼0.1 was achieved.
Observed differential size (Ne ) spectra for different zenith angles θ , grouped in bins ofsec(θ ) of width 0.05, are shown in ﬁgure 6. Flux values shown on the Y-axis were multiplied by
Ne 2.5 for a clearer presentation. Statistical errors in data in ﬁgure 6 are very small, since each
Ne , sec(θ) bin contains >105 EAS. A comparison of the nine spectra for sec(θ ) from 1.025(median of ﬁrst bin) to 1.375 (for last bin) shows a close similarity, apart from differences in6

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

1.0(Nfit / Ntrue )

0.8

0.6

0.4

0.2

0
4.2

4.4

4.6

4.8
5.0
LOG10(Ne )

5.2

5.4

dI/dN N2.5 (m-2sec-1sr-1Ne1.5)

Figure 5. Ratio of ‘ﬁtted’ and ‘true’ size from CORSIKA for 10 size groups. The X-axis is thelogarithm of size.

400

200

100
60

30

4.5

4.7

4.9
Log10(Ne)

5.1

5.3

Figure 6. Observed Ne spectra for sec(θ ) intervals of 0.05. Top plot
1.00<sec(θ )<1.05, . . . , last plot 1.35<sec(θ )<1.40.

vertical, second plot

ﬂuxes observed. Since the spectral slope is independent of θ , we combined the data for EAS
with sec(θ ) < 25◦ , after suitable normalization, for further analysis of the muon content. The
vertical ﬂux shown in ﬁgure 6 derived from the EAS of θ < 25◦ was converted into equivalentvertical size by following the relation, where vertical depth Xv = 790 g cm−2 and attenuation
length λ = 200 g cm−2 :Xv

N(0) = N(θ ) e λ (sec θ−1) .

(10)

Observed muons in 16 modules were counted by reconstruction of tracks based on hits
in the four layers of PCs as described earlier. A muon outside the time window from −0.8
to 2.6 μs, relative to the EAS or if the direction was >10◦ from the EAS in any of the twoplanes was rejected as unassociated. EAS with cores between 60 and 80 m from the center ofmuon detector and inside a 30 m radius circle at the center of the array were selected. Thiskidney-shaped area of selection is shown as shaded in ﬁgure 1. The resultant MMD is shown inﬁgure 9 for the size 105.0 � Ne � 105.2 . Similar MMDs were obtained for the remaining sizes
7

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et alas shown in ﬁgure 6. MMDs were used for measuring energy spectra of ﬁve nuclei in PCRs.
This was done by comparing observed MMDs with CORSIKA simulations after including thedetector response from Geant4 [13] as discussed in the next section.
4. Monte Carlo simulations
Most ground-based experiments record electrons, muons, hadrons and other components
deep down in the atmosphere (∼8–10λ). Thus, detected particles at the observational level
are the end product of hundreds of interactions at different depths in the atmosphere.
Thus, the extraction of information on PCRs and characteristics of the ﬁrst few ultra-highenergy interactions require a comparison of observations with results with simulations, using
generators of interactions for an assumed energy spectra of PCRs. Since the extreme forward
region of interactions in the laboratory frame relevant here cannot be studied in high energy
collider experiments, models of high energy interactions need to be extrapolated from energies
of a TeV to 1000 TeV and beyond.
Since CORSIKA parameters to describe hadronic interactions are extrapolated from lower
energy measurements, a signiﬁcant model-dependent uncertainty in simulation results is to
be expected. CORSIKA version 6.020 has several built-in hadronic generators [19]. The ones
considered here include quark–gluon string-based QGSJet01 [20], dual parton with mini-jet
based SIBYLL-2.1 [21], NEXUS-2.01 [22] generators, that hereafter would be referred to as
QGSJet, SIBYLL, NEXUS, respectively. The low energy hadrons (<80 GeV) were treatedusing GHEISHA [23]. In CORSIKA, the electromagnetic cascade is treated by two methods,
EGS4 code [24] or NKG function [17]. In EGS4, every particle is tracked but NKG calculatesdensities at the array level and runs faster than EGS4. To save computing time, we have used
NKG. But we performed simulations using EGS4 for a limited number of events, to investigate
differences between the two methods, after the folding detector response [13, 25]. Saturation
effects on a number of muons at high multiplicities were accounted for. This comparison
showed that ﬁnal results were insensitive to the speciﬁc method used.
The PCRs were assumed to be of ﬁve species, protons (H, A = 1), helium (He, A = 4),
nitrogen (N, A = 14), aluminum (Al, A = 27) and iron (Fe, A = 56). Heavier nuclei N, Al andFe were used to represent light (CNO), medium (Mg, Al, Si) and heavy (Mn, Fe, Co) masses in
PCRs. In simulations, PCR energies were randomly selected from a differential power law ofslope −2.7 for the ﬁve nuclei. The selected threshold energies for H, He, N, Si, Fe were 5, 8, 16,
16, 20 TeV, respectively. The thresholds were sufﬁciently small, so that trigger efﬁciencies were<10% [13]. The core of each EAS was thrown randomly over an area of radius 100 m centeredon the array. Each EAS was subjected to analysis identical to the one used for the observed
EAS. Once the GRAPES-3 trigger criterion was met, then EAS parameters, Ne , Nμ etc were
calculated to obtain simulated MMDs. Since MMDs for Al and Fe overlap substantially, they
were combined by assuming a ratio of Al/Fe = 0.8 based on a direct experiment [3]. Sinceabundance ratio Al/Fe was ﬁxed, effectively the number of independent species was reduced
to 4. EAS with θ � 30◦ were simulated, but only those with θ �25◦ were used in analysis.
Muon and hadrons were tracked until their energy dropped below 1 GeV. Simulated MMDs
for ﬁve nuclei were normalized per EAS and the fraction ai of each nuclei was extracted byminimizing χ 2 using MINUIT to reproduce the observed MMD for each size Ne bin:��2
�ai nik
n0k −
�
i
2
χ =
(11)
�
2
(ai ik )2
k
0k +
i

8

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

1000
(a)

(c)
H (1000 TeV)

Frequency

500

H (100 TeV)

0
(b)

1000

(d)
Fe (100 TeV)

Fe (1000 TeV)

500
0
3.5

4

4.5

5

5.5

6

LOG10(Ne )Figure 7. Size distribution at Ooty (790 g cm−2 , θ � 30◦ ), (i) QGSJet continuous line, SIBYLL
dashed line, NEXUS dotted line. Primary particle and energy: (a) H 100 TeV, (b) Fe 100 TeV,(c) H 1000 TeV, (d) Fe 1000 TeV.

where k is the multiplicity bin number, i the mass group, i.e. i = 1−5 for H, He, N, Al, Fe,respectively, and i = 0 is the observed data. nik and ik are the mean frequency and statisticalerror for the kth bin. The initial value of parameter ai was set to 0.2 and the effect of the initialvalue on the ﬁnal result was found to be negligible.
It was pointed out by the CORSIKA team [16] that the fragmentation of the incident
nucleus, from full fragmentation in the ﬁrst interaction (superposition model) to partial
fragmentation into other nuclei, led to essentially the same results for most observables deep
down in the atmosphere, where EAS are observed. Any distinctive features of fragmentation
are smeared out by ﬂuctuations in the huge number of interactions, during development andpropagation of EAS in the upper atmosphere. Therefore, we use a simpler superposition modeland treat interactions of an incident nucleus of mass A, energy E as the superposition of A
nucleons, each of energy ‘E/A’ from the ﬁrst interaction.
Before a detailed comparison of observed MMDs for different Ne with simulations, it
is useful to examine the broad features of predictions of different generators. Distributionsof simulated Ne for H and Fe, for an energy of 100 TeV, are shown in ﬁgures 7(a) and (b)by continuous, dashed and dotted lines for QGSJet, SIBYLL and NEXUS, respectively. On
average, QGSJet predicts a smaller Ne at Ooty compared to the other two generators. The
ratio of Ne for H and Fe at 100 TeV is ∼3. A similar result is obtained at 1000 TeV as seen inﬁgures 7(c) and (d). Compared to QGSJet, Ne is ∼10% larger for SIBYLL and is the largest
for NEXUS. Comparison of Ne distributions at 100 and 1000 TeV shows a narrower one for
Fe than H for all three generators. This behavior may be understood as follows: ﬂuctuations
in depth of the ﬁrst interaction are smaller for Fe- than H-induced EAS due to superpositionof a large number of cascades for Fe that reduces the ﬂuctuations because of averaging.
A similar comparison of muons of energy > 1 GeV, the so-called muon size Nμ , in H and
Fe EAS for the same generators, QGSJet, SIBYLL and NEXUS, is shown in ﬁgures 8(a) and(b) for 100 TeV and in ﬁgures 8(c) and (d) for 1000 TeV. Nμ on average is signiﬁcantly larger
9J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

1000
(a)

(c)
H (100 TeV)

H (1000 TeV)

Frequency

500

0
(b)

1000

(d)
Fe (100 TeV)

Fe (1000 TeV)

500
0
0

1000

2000

10000

20000

Number of Muons
Figure 8. Distributions of �1 GeV muons at Ooty (790 g cm−2 , θ � 30◦ ) for (i) QGSJet, continuousline, SIBYLL, dashed line, NEXUS, dotted line. Primary particle and energy: (a) H 100 TeV, (b)Fe 100 TeV (c) H 1000 TeV, (d) Fe 1000 TeV.

for QGSJet than SIBYLL or NEXUS. This reversal of trend further enhances the differencesin MMDs when EAS of same Ne are compared, thus increasing the sensitivity of the MMD
technique. At 100 TeV, muon distribution is signiﬁcantly narrower for Fe- than H-induced
EAS as shown in ﬁgures 8(a) and (b). The same trend is repeated at 1000 TeV as shown inﬁgures 8(c) and (d), although this size-based selection does reduce the number of Fe EAS in
the sample. In the next section, we compare experimental data with simulations using only
QGSJet and SIBYLL since they were reported to provide better agreement with observational
results [11, 27].
5. Results
GRAPES-3’s objective is to compare indirect spectra with direct results. Thus, the numberof simulated EAS were normalized to the number expected during the observed time of
4.71×107 s. The detected EAS of a size Ne are more heavily dominated by H than Fe, as
discussed in section 4 and shown in ﬁgure 7 and their ratio is ∼3 for H and Fe at 100–1000
TeV. The observed MMD for Ne 105.0 –105.2 is shown as ﬁlled circles in ﬁgure 9(a). The error
bars represent only statistical errors in data. Also shown are results of the simulated EAS
initiated by H, He, N, Al and Fe using QGSJet. Relative abundances of four primary groups(Al/Fe ratio being ﬁxed at 0.8) were obtained from a χ 2 minimization by ﬁtting 4+1 simulated
MMDs to the observed MMD displayed with error bars in ﬁgure 9(a). The summed distribution
shown as the top histogram in ﬁgure 9(a) gives an excellent ﬁt to the observed MMD for thefull range of multiplicities from 0 to 60 muons. A comparison of the same observed MMD
with the summed MMD for SIBYLL is shown in ﬁgure 9(b). This summed MMD also yields
an equally good ﬁt, although with different relative abundances of ﬁve nuclei.
A comparison of ﬁgure 9(a) and ﬁgure 9(b) shows that far fewer heavier nuclei are
required to ﬁt large (nμ > 40) multiplicities for QGSJet than SIBYLL. Smaller (nμ < 10)10

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

1.0
10-1
Normalized Shower Rate

(a)

H

10-2

He
N

10-3
10-4

Fe
Al

10-5
10-1

(b)

10-2
N

10-3

Al

He
Fe

10-4

H

-5

10

0

10

20
30
40
Detected Muons

50

60

Figure 9. A comparison of the observed MMD for Ne 105.0 –105.2 with simulations for H, He, N,
Al and Fe for (a) QGSJet, (b) SIBYLL. Note enhanced Al, Fe in SIBYLL.

multiplicities are entirely accounted for by lighter nuclei in both cases. The dominance ofeach mass group at different multiplicities underlies the sensitivity of this technique for thecomposition. Figures 9(a) and (b) also indicate evidence of heavier nuclei that generate the EAS
in the 105.0 –105.2 region. This shows that interpretation of the same data requires signiﬁcantly
different compositions for the two generators. Similar behavior is seen for other sizes. Thus,
the derived composition is intimately linked to the hadronic generator used. Since it was notpossible to identify a unique hadronic generator from the GRAPES-3 data, we resorted to otherconstraints such as its consistency with direct results to identify the correct generator. Using
the above procedure, relative abundances of four nuclei were obtained for QGSJet simulations
for ﬁve sizes, 104.2 –104.4 , 104.4 –104.6 , 104.6 –104.8 , 104.8 –105.0 and 105.0 –105.2 . The variationin abundance of four nuclei for Ne 104.2 –105.2 is shown by ﬁlled squares in ﬁgure 10. Resultsfor SIBYLL are shown by ﬁlled circles. Signiﬁcantly different abundances for two generators
are seen in ﬁgure 10.
The data show different trends for composition dependence of Ne . H is more abundant for
QGSJet than SIBYLL, due to ∼15% more muons in QGSJet than SIBYLL. This is reﬂectedin the results shown in ﬁgure 10, where the H fraction for QGSJet is >70% at Ne = 104.3 .
Thereafter, the H fraction decreases with increasing Ne for both generators. The abundanceof H sharply decreases by ∼30% from Ne = 104.3 –105.1 for both generators. This suggests
a steeper spectra for H, than other nuclei. A fraction of He increases with Ne for QGSJet,
but stays almost constant for SIBYLL. A different dependence is observed for the N fraction,that increases with Ne for SIBYLL, but is almost constant for QGSJet. The heavy fraction(Al+Fe) is much larger for SIBYLL than QGSJet. But a small fraction of ∼5% (Al+Fe) at
an Ne corresponds to a larger fraction (�15%) for the corresponding primary energy. Thisis because Fe-EAS develop faster, reach maximum Ne higher and then attenuate faster inlower atmosphere when detected. This leads to ∼30% size of the H EAS for the same energy.
However, the heavy fraction is much smaller at larger Ne for QGSJet. These results may havesigniﬁcant implications for the nature of the sources of PCRs and the acceleration processes
dominating over the region 100–1000 TeV [28]. Next, observed Ne spectra were converted into11

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

90
H

70

Q

Relative Abundance (%)

50
30

S

40
30
20
10

S
Q
He

N

30

S

20
10

Q
Al+Fe

10

S
Q

0
4.2

4.4

4.6

4.8

5.0

5.2

Log10(Ne)
Figure 10. Relative abundance of four groups with Ne by
SIBYLL(S).

for QGSJet(Q) and by

for

3.0
Fe
Al
N

2.7
LOG10 (ETeV)

He
H
2.4

2.1

1.8

1.5
4.2

5

5.2

Figure 11. Size and energy relationship for ﬁve nuclei, (a) H, (b) He(e) Fe

4.4

4.6
4.8
LOG10(Ne ), (c) N, (d) Al,for SIBYLL. The dotted lines are there only to guide the eye.

energy spectra for ﬁve nuclei. Taking account of the exposure factor for data used in MMDs
in ﬁgure 9 for a size and for the two generators, it is easy to convert relative abundances
in ﬁgure 10 to ﬂuxes for ﬁve nuclei at corresponding primary energy. The energy of eachlogarithmic Ne bin is computed from the logarithm of energy of the EAS in the bin, after
imposing selection cuts discussed above.
An energy is assigned to each Ne , using size-to-energy relation from simulations for ﬁve
nuclei and for two generators. The dependence of energy E0 on Ne for ﬁve nuclei, namely
H, He, N, Al, Fe for SIBYLL is shown in ﬁgure 11. The dotted lines joining data for individual
nuclei are only to guide the eye. The ratio of energies of He and H at Ne = 104.3 is 1.5. Similarly
12

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

104
3

10

H

102

10
CNO

1

E

2.5

dI/dE [m-2sr -1s-1 (GeV/n)

1.5

]

He

Al

10-1

1
10-1

Fe
2

10

103

104

105

106

Energy E (GeV/n)
Figure 12. H, He, CNO, Al, Fe spectra from direct experiments: Ryan [29],
JACEE [2],RUNJOB [3], CREAM [4],TIBET [7], KASCADE [11]QGSJet and GRAPES-3 SIBYLL, QGSJET.

SOKOL [30],SIBYLL,ratios for N, Al, Fe relative to H are 2.1, 2.5, 3.0, respectively, at Ne = 104.3 . The magnitudeof these ratios decreases slightly with increasing Ne , reducing by ∼10% at 105.1 as seen inﬁgure 11. Next, Ne was converted into median energy using this relation. Energy spectra soderived for GRAPES-3 for H and He are shown in the top and second panel of ﬁgure 12 for
QGSJet by ﬁlled squares and SIBYLL by ﬁlled circles. Spectra for N, Al, Fe are shown in thethird, fourth and ﬁfth panels in ﬁgure 12, respectively. These are compared with the results of
the other groups: H and He data are from [2, 3, 29, 30] and N, Al and Fe are from [1, 3, 4, 30].
Also shown are results of indirect experiments, Tibet-ASγ on H, He [7] and KASCADE on H
[11]. It may be recalled that GRAPES-3 spectra for Al and Fe are not independent since theirrelative abundance was ﬁxed at 0.8, but are plotted separately for comparing them with direct
results that are displayed in the ﬁve panels of ﬁgure 12.
A good agreement on composition and energy spectra of PCRs is seen between
GRAPES-3 and direct results using SIBYLL, especially with JACEE, SOKOL and CREAM
results [2, 30, 4]. This comparison shows that SIBYLL provides a better description of hadronicinteractions in the overlap region of 100–300 TeV. Ours are absolute ﬂuxes that were derived
from Ne and MMD directly, without recourse to any normalization. The extension of energy
spectra by a factor of ∼3 for ﬁve nuclei from GRAPES-3 is consistent with extrapolationfrom direct methods. Although Ne spectra shown in ﬁgure 6 can be ﬁtted to a single powerlaw, when combined with MMDs, different spectral slopes emerge for the 5 nuclei as seen13

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

E

2.5

-2

dI/dE [m sr -1s -1 GeV

1.5

]

104

103
102

3

10

5

104

10

6

10

Energy E (GeV)

Figure 13. All-particle spectra from direct measurements:
JACEE [2],RUNJOB [3], and GRAPES-3

SIBYLL,

Grigorov [31],

SOKOL [30],QGSJET, TIBET �.

in ﬁgure 12. The GRAPES-3 data are consistent with Tibet-ASγ in the overlap region of
200–700 TeV/n for H [7]. GRAPES-3 overlaps with KASCADE at 900 TeV and the H ﬂux
from KASCADE for SIBYLL and QGSJet [11] agrees closely with corresponding ﬂuxes for
GRAPES-3 as shown in ﬁgure 12. The consistency of our results using SIBYLL with direct
methods underscores the utility of the energy overlap in the selection of a suitable generator.
A comparison of the GRAPES-3 all-particle spectra with direct results [2, 3, 30, 31] shownin ﬁgure 13 also display good agreement, except for RUNJOB [3]. However, this comparisonof all-particle spectra allows no distinction to be made between QGSJet and SIBYLL by the
GRAPES-3 data.
6. Discussion
The analysis of 6×108 EAS and their muon content measured by GRAPES-3 is used to obtain
MMDs as a function of Ne (104.2 –105.2 ), corresponding to the energy range 100–1000 TeV.
Simulated data in ﬁgure 7 show that for the same energy (100, 1000 TeV) H produces larger
Ne EAS than Fe. Thus, if EAS are selected on Ne , then the sample contains H of ∼3 times
lower energy than Fe. For ﬁxed energy (100, 1000 TeV), Nμ is a factor of ∼2 larger for Fe
than H as shown in ﬁgure 8. Thus, if EAS are binned in narrow Ne intervals and their MMDs
examined, a muon multiplicity of approximately six times larger would be observed for Fe
than H. However, this Ne based cut results in the selection of Fe of approximately three times
higher energy and thus their fraction substantially reduces in the Ne bin, a limitation of thistechnique. However, even in the presence of a relatively small fraction of Fe, the MMD gets
substantially enhanced at higher multiplicities. Such high multiplicities cannot be producedby lighter elements as shown in ﬁgure 9 for 105.0 < Ne < 105.2 . This fact summarizes thesensitivity of GRAPES-3 to the mass composition of PCRs.
The above argument especially applies to the GRAPES-3 muon detector because(i) the identiﬁcation of each muon is unambiguous, (ii) the muon statistics in an EAS are
good due to the large area (560 m2 ), (iii) the muon detector is compact, wherein the 560 m2
detector is conﬁned in a geographical area of ∼1600 m2 . This allows the simulation of muons
to be largely independent of details of the assumed shape of lateral distribution. As seen from14

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

ﬁgure 12, the abundances of Al and Fe are signiﬁcantly lower for QGSJet than for SIBYLL.
This difference is caused by the production of ∼15% more muons in QGSJet. Thus, QGSJet ﬁts
observed MMD even with a relatively lighter composition rather well as explained in section 5.
It is to be reiterated that the derived composition is strongly dependent on the nature of assumed
hadronic interactions. However, this is not a limitation unique to the GRAPES-3, but holds
true for all EAS experiments. As our knowledge of hadronic interactions improves, especiallywith the start of the LHC, measurements of composition should also become progressively
more precise.
A glance at ﬁgure 12 shows that predictions of both QGSJet and SIBYLL seem toagree with direct results for lighter nuclei within errors. However, a closer examination of
GRAPES-3 data reveals that the ﬂux of H with QGSJet is systematically higher by a factorof ∼2 than with SIBYLL. Despite the errors, the ﬂux obtained using SIBYLL is much closer
to the direct results than QGSJet and that is also true for the KASCADE ﬂux at 900 TeV.
This observation is especially relevant for comparison of H spectra for the following two
important reasons: (i) ﬂuxes from direct methods are intrinsically more reliable, especiallyat the lower end of the overlap energy region, and (ii) no major assumptions are needed fornuclear break-up in the ﬁrst H air–nucleus collision.
QGSJet ﬂux of CNO is systematically lower by a factor of ∼3 compared to SIBYLL. Hereagain SIBYLL is closer to more reliable direct results from JACEE and CREAM. Only in thecase of He do both QGSJet and SIBYLL yield nearly the same ﬂuxes consistent with the direct
results. However, for heavier Al and Fe, the agreement with the direct results for SIBYLL issigniﬁcantly better than QGSJet as seen in ﬁgure 12. In short, comparison of spectra for ﬁve
nuclei obtained from MMDs with direct results shows the inadequacy of QGSJet to provide
good agreement with data. But there is no signiﬁcant difference in the spectra of He for bothgenerators. However, differences between them widen for the three heavier nuclei.
Among the four direct experiments, JACEE, SOKOL, CREAM, RUNJOB, the
GRAPES-3 results agree well with the ﬁrst three experiments. But there is a discrepancy inﬂuxes of He reported by RUNJOB and JACEE/SOKOL/CREAM. GRAPES-3 results provide
continuity to measurements from JACEE/SOKOL/CREAM. Another compelling inference
from comparisons shown in ﬁgure 12 relates to the difference in steepness of the H spectrum
relative to the other four nuclei that are decidedly ﬂatter up to ∼1000 TeV. This difference inspectral slopes may have important implications for the nature of cosmic ray sources and theenvironment surrounding them in the range 100–1000 TeV within the galaxy. The GRAPES-3,
H spectrum agrees well with the results of other indirect experiments, KASCADE [11] andTibet AS/γ [7]. However, the He ﬂux from Tibet AS/γ is in disagreement, with GRAPES-3
as well as JACEE, SOKOL and CREAM, probably showing the inadequacy of the Ne spectrum
determined from an EAS array alone, to provide reliable spectral information on He and otherheavier nuclei in PCRs.
7. Conclusions
The observations of MMDs with the large GRAPES-3 tracking muon detector have provided
evidence of a steeper proton spectrum, consistent with direct measurements from balloon andsatellite borne experiments. GRAPES-3 spectra for He, N, Al and Fe agree well with the direct
results, extending the energy reach and conﬁrming their ﬂatter nature relative to protons. Our
results indicate SIBYLL to be a suitable generator for describing hadronic interactions in theoverlap energy region of 100–300 TeV. The extension of energy spectra by a factor of ∼3
up to ∼1000 TeV for ﬁve nuclei is consistent with the extrapolation of the spectra measured
by direct methods. With greater conﬁdence in the suitability of SIBYLL at these energies,
15

J. Phys. G: Nucl. Part. Phys. 39 (2012) 025201

H Tanaka et al

analysis of larger size (Ne > 105.4 ) showers will be done for extending the measurements to
the knee region and beyond, which may help in better understanding the physical process(es)responsible for this critically important feature in the energy spectrum of PCRs.
Acknowledgements
We thank D B Arjunan, G P Francis, I M Haroon, A Iyer, P Jagadeesan, V Jeyakumar,
K Manjunath, S D Morris, S Murugapandian, B Rajesh, K Ramadass, C Ravindran, V SanthoshKumar, R Suresh Kumar and V Viswanathan for their help in the installation and operation of
GRAPES-3. This work is dedicated to the memory of our brilliant colleague S Karthikeyan
who passed away in a tragic road accident.
References[1] Kobayakawa K et al 2002 Phys. Rev. D 66 083004[2] Asakimori K et al 1998 Astrophys. J. 502 278
Takahashi Y 1998 Nucl. Phys. B 60 83[3] Shibata T et al 2001 Astropart. Phys. 16 13
Derbina V A et al 2005 Astrophys. J. Lett. 628 L41[4] Ahn H S et al 2009 Astrophys. J. 707 593
Ahn H S et al 2010 Astrophys. J. Lett. 714 L89[5] Glasmacher M A K et al 1999 Astropart. Phys. 10 291[6] Aglietta M et al 1999 Astropart. Phys. 10 1
[7] Amenomori M et al 2000 Phys. Rev. D 62 112002[8] Bacci C et al 2000 Nucl. Instrum. Methods Phys. Res. A 443 342[9] Fowler J W et al 2001 Astropart. Phys. 15 49[10] Gupta S K et al 2003 Phys. Rev. D 68 052005[11] Antoni T et al 2005 Astropart. Phys. 24 1
[12] Chilingarian A et al 2007 Astropart. Phys. 28 58[13] Gupta S K et al 2005 Nucl. Instrum. Methods Phys. Res. A 540 311[14] Hayashi Y et al 2005 Nucl. Instrum. Methods Phys. Res. A 545 643[15] Hayashi Y 2005 Proc. 29th Int. Cosmic Ray Conf. (Pune) vol 10 p 243[16] Heck D et al 1998 Report FZKA 6019 Forschungszentrum Karlsruhe[17] Greisen K 1960 Annu. Rev. Nucl. Sci. 10 63
Kamata K et al 1958 Suppl. Prog. Theor. Phys. 6 93
[18] Oshima A et al 2010 Astropart. Phys. 33 97[19] Knapp J et al 1996 Report FZKA 5828 Forschungszentrum Karlsruhe
Heck D et al 2001 Proc. 27th Int. Cosmic Ray Conf. (Hamburg) vol 1 p 233[20] Kalmykov N N et al 1997 Nucl. Phys. B 52 17[21] Fletcher R S et al 1994 Phys. Rev. D 50 5710
Engel R et al 1999 Proc. 26th Int. Cosmic Ray Conf. (Salt Lake City) vol 1, p 415[22] Drescher H J et al 2001 Phys. Rep. 350 93[23] Fesefeldt H 1985 Report PITHA-85/02 RWTH Aachen[24] Bielajew A F et al 1994 National Research Council Canada Report NRC-PIRS-0436[25] Yoshikoshi T et al 2001 Proc. 27th Int. Cosmic Ray Conf. (Hamburg) vol 2 p 612[26] James F 1998 MINUIT CERN Prog. Library D506
[27] Milke J et al 2004 Acta Phys. Polon. B 35 341[28] Meli A and Biermann P L 2006 Astron. Astrophys. 454 687[29] Ryan M J et al 1972 Phys. Rev. Lett. 28 985[30] Ivanenko I P et al 1993 Proc. 23rd Int. Cosmic Ray Conf. (Calgary) vol 2 p 17
[31] Grigorov N L et al 1971 Proc. 12th Int. Cosmic Ray Conf. (Hobart) vol 5 p 1746

16