SciPost Phys. Proc. 13, 031 (2023)

Updated results on the cosmic ray energy spectrum andcomposition from the GRAPES-3 experiment
F. Varsi1⋆ , S. Ahmad2 , M. Chakraborty3 , A. Chandra2 , S. R. Dugad3 , U. D. Goswami4 ,S. K. Gupta3 , B. Hariharan3 , Y. Hayashi5 , P. Jagadeesan3 , A. Jain3 , P. Jain1 ,S. Kawakami5 , H. Kojima6 , S. Mahapatra7 , P. K. Mohanty3 , R. Moharana8 , Y. Muraki9 ,P. K. Nayak3 , T. Nonaka10 , A. Oshima6 , B. P. Pant8 , D. Pattanaik3,7 , G. Pradhan11 ,M. Rameez3 , K. Ramesh3 , L. V. Reddy3 , R. Sahoo11 , R. Scaria11 , S. Shibata6 ,K. Tanaka12 and M. Zuberi3
1 Indian Institute of Technology Kanpur, Kanpur 208016, India
2 Aligarh Muslim University, Aligarh 202002, India
3 Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai 400005, India
4 Dibrugarh University, Dibrugarh 786004, India
5 Graduate School of Science, Osaka City University, Osaka 558-8585, Japan
6 College of Engineering, Chubu University, Kasugai, Aichi 487-8501, Japan
7 Utkal University, Bhubaneswar 751004, India
8 Indian Institute of Technology Jodhpur, Jodhpur 342037, India
9 Institute for Space-Earth Environmental Research, Nagoya University,
Nagoya 464-8601, Japan
10 Institute for Cosmic Ray Research, Tokyo University, Kashiwa, Chiba 277-8582, Japan
11 Indian Institute of Technology Indore, Indore 453552, India
12 Graduate School of Information Sciences, Hiroshima City University,
Hiroshima 731-3194, Japan⋆ fahim@iitk.ac.in
21st International Symposium on Very High Energy Cosmic Ray Interactions(ISVHECRI 2022)
Online, 23-28 May 2022
doi:10.21468/SciPostPhysProc.13

Abstract
Here, we present the updated results on the cosmic ray energy spectrum and composition
analysis from the GRAPES-3 experiment over the energy range of 50 TeV to 1000 TeV
since ICRC 2021. The simulation of cosmic ray showers was performed using the postLHC high energy hadronic interaction model QGSJetII-04 and low energy hadronic model
FLUKA. A detailed GEANT4 simulation of the GRAPES-3 muon telescope was performed.
The composition was obtained by fitting simulated muon multiplicity distributions for
proton, helium, nitrogen, aluminium, and iron primaries to the observed data. Theenergy spectrum connects to the direct measurements with a fairly good agreement influx.
Copyright F. Varsi et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

031.1

Received 30-08-2022
Accepted 22-06-2023
Check forupdatesPublished 29-09-2023
doi:10.21468/SciPostPhysProc.13.031

SciPost Phys. Proc. 13, 031 (2023)

1

The GRAPES-3 experiment

The GRAPES-3 (acronym: Gamma Ray Astronomy at PeV Energies Phase-3) experiment islocated at Ooty, South India (11.4◦ N latitude, 76.7◦ E longitude, and 2200 m altitude above
m.s.l.). It consists of a most compact array of 400 plastic scintillator detectors, each with anarea of 1 m2 and an inter-detector separation of 8 m, spreading over a total area of 25000
m2 [1–3]. Being a highly dense EAS array, the GRAPES-3 experiment can observe the primarycosmic rays (PCRs) in TeV-PeV energy range, providing a significant overlap with direct experiments [4–6]. It also contains 16 similar independent muon modules, forming the largest muon
telescope (G3MT) with an area of 560 m2 [7]. The GRAPES-3 muon telescope (G3MT) has anenergy threshold of sec θ GeV, where θ is the zenith angle of the incident muon. The muon
multiplicity distribution observed by the G3MT is sensitive to PCRs mass and an importantparameter for the precise measurements of the PRCs composition [8].

2MC simulations

The CORSIKA package version 7.6900 was used to simulate the extensive air showers (EAS)for proton (H), helium (He), nitrogen (N), aluminium (Al) and iron (Fe), using the QGSJETII-04 and FLUKA models as high and low energy hadronic interaction models, respectively. For
each simulated primary, nearly 6×107 EAS were simulated in an energy range of 1 TeV - 10
PeV divided in 20 equal logarithmic bins of a width of 0.2. The number of EAS simulated ineach energy bin is listed in Table 1. In each energy bin, EAS were generated in zenith angle
range of 0◦ - 45◦ , assuming a power law of differential spectral index −2.5. For the analysis,
each EAS was thrown ten times in a circular area of radius 150 m from the center (−13.85 m,
6.29 m) of the GRAPES-3 scintillator detectors array with a random core position to improve
the statistics. The CORSIKA simulation output was further analyzed to obtain particle density,using a detailed simulation of the GEANT4 response of γ-rays, electrons, muons, and hadrons
in the plastic scintillator detectors. A detailed GEANT4 simulation of the response of the secondaries in the G3MT had also been performed. Since the top concrete of G3MT absorbs theelectromagnetic component of the EAS, the GEANT4 response of muons and hadrons was measured in terms of PRC hits, using a discriminator threshold of 0.2 times of the energy depositedby the minimum ionizing particle (∼4 keV). This PRCs hit status was used to reconstruct themuon tracks (statistical measure of the muons).

3

Reconstruction procedure, quality cuts and data summary

For each triggered EAS, the parameters such as core location (X c , Yc ), size (Ne ), and age (s)are obtained by minimizing the observed lateral distribution of particles densities in different
scintillator detectors with a lateral density distribution function, namely Nishimura-KamataGreisen (NKG) as given below, based on a log-likelihood method using TMINUIT package,Γ (4.5 − s)ρ(ri ) =2 Γ (s)Γ (4.5 − 2s)
2πr M
Ne



ri
rM

(s−2) ri (s−4.5)
1+
,
rM

(1)

where ρ(ri ) is the particle density observed by i th detector, ri is the lateral distance of i thdetector from shower core, and r M is the Moliere radius which is 103 m for the GRAPES-3
observational site. The initial estimate of arrival direction (θ , φ) is measured by fitting therelative arrival time of secondary particles recorded by different scintillator detectors with a

031.2

SciPost Phys. Proc. 13, 031 (2023)

Table 1: Summary of the number of EAS simulated using CORSIKA for each of the
PCRs (H, He, N, Al, and Fe).

Bin

Energy range [TeV]

Number of EAS

Bin

Energy range [TeV]

Number of EAS

1
2
3
4
5
6
7
8
9
10

1.00 - 1.58
1.58 - 2.51
2.51 - 3.98
3.98 - 6.31
6.31 - 10.00
10.00 - 15.85
15.85 - 25.12
25.12 - 39.81
39.81 - 63.10
63.10 -100.00

2.5×107
1.5×107
1.0×107
5.0×106
2.5×106
1.5×106
1.0×106
5.0×105
2.5×105
1.5×105

11
12
13
14
15
16
17
18
19
20

100.00 - 158.49
158.49 - 251.19
251.19 - 398.11
398.11 - 630.96
630.96 - 1000.00
1000.00 - 1584.89
1584.89 - 2511.89
2511.89 - 3981.07
3981.07 - 6309.57
6309.57 - 10000.00

1.0×105
5.0×104
2.5×104
1.5×104
1.0×104
5.0×103
2.5×103
1.5×103
1.0×103
5.0×102

plane front, followed by correction for the shower front curvature based on shower size andage to obtain a more accurate EAS direction [9].
To reduce the systematics cause by EAS reconstruction, event selection criteria are applied. The EAS having successful parameters and arrival direction reconstruction are used in
the analysis. Shower cores are restricted to a circular area of 50 m from the center of the arrayto minimise the contamination of the improper reconstruction due to EAS with core landing
near the edge or outside the array. The shower age is restricted between 0.2 and 1.8 to avoid
improper reconstruction due to shower age converging to its limits. The zenith angle is retricted to 18◦ to minimize the systematics due to inclined EAS. To ensure the trigger efficiency> 90%, the EAS having shower size ≥ 104 are used.
Data collected during 1 January 2014 - 26 October 2015 (∼ 22 months) is used for the
analysis. The total live time of data collection is ∼463 days. The number of showers remaining
after applying all the quality cuts are 1.47×107 from a total set of 1.75×109 EAS.

4

Analysis

For each simulated primary, the trigger (ϵ T ) and reconstruction (ϵR ) efficiencies are calculated
as a function of primary energy, and total efficiency (ϵ t ot ) is determined by the product oftrigger and reconstruction efficiencies. Acceptance(Aacc ) is represented as the product of theeffective area and the effective viewing angle as,nθ

Aacc (E T ) =πA X
ϵ t ot (E T , θk ) (cos 2θk − cos 2θk+1 ) ,
2 k=1

(2)

where A is the fiducial area, nθ is the total number of angle bins and θk and θk+1 are low
and high edges of each angle bin, respectively. The trigger efficiency and acceptance of all
simulated primaries are shown in Figure 1. The trigger efficiency is >90% at 40 TeV, 45 TeV,
60 TeV, 70 TeV and 85 TeV for H, He, N, Al and Fe, respectively. The acceptance is increased
to 2300 m2 sr at 100% efficiency for θ < 18◦ .
The G3MT is dedicated to measure the number of muon tracks for each triggered shower.
With an increase in the number of incident muons in a given module, the muon tracks start
to overlap, resulting in the underestimation of reconstructed tracks. This saturation effect, for
one module, is shown in the left panel of Figure 2 for simulations. The curve is modelled with031.3

Acceptance [m2 sr]

Trigger efficiency (%)

SciPost Phys. Proc. 13, 031 (2023)

100
80
60
40

103

105

104

2000

1500

1000

protonHeliumNitrogenAluminium
Iron

20

2500

protonHeliumNitrogenAluminium
Iron

500

106
107
Primary energy [GeV]

103

104

105

106
107
Primary energy [GeV]20

15

χ2 / ndf
p0
p1
p2
p3

Normalized shower rate

Incident muons

Figure 1: Trigger efficiency (left) and acceptance of GRAPES-3 scintillator detectors
array (right) as a function of primary energy for θ < 18◦ .

2.006 / 5
− 0.2016 ± 0.4174
1.162 ± 0.2144− 0.05336 ± 0.03464
0.00759 ± 0.001746

10

5

log10 (Ne): 4.4 - 4.6

10−1

data
proton
iron

−2

10

10−3
10−4
10−5
10−6

0

2

4

6

8

10
12
Detected muons

0

10

20

30
40
Muon multiplicity

Figure 2: Left: Muon saturation for one module of the G3MT, modelled with thirdorder polynomial represented by the black curve. The blue dashed line guides thelinearity between the incident and detected muons. Right: Muon multiplicity of H
and Fe fitted with the negative binomial distribution, plotted along with observed
MMD for 4.4 ≤ log(Ne ) < 4.6.
a third-order polynomial for detected muons ≥ 4, which is used to get the correct estimate of
the muon multiplicity for both simulation and observed data. The shape of muon multiplicity
distribution is well described by the negative binomial distribution (NBD) and is given by,2 m  2m2  σ2 − m  x
Γ x + σ2m−m

σ −m
2
 2 
N BD x; m, σ =
,
(3)
2
m
σ
σ2
Γ (x + 1)Γ 2
σ −mwhere x is the muon multiplicity value and m and σ2 are the mean value and standard deviation of muon multiplicity distributions (MMDs), repectively. Thus, the normalized MMD ofeach simulated primary is fitted with NBD to model statistical fluctuations and shown in theright panel of Figure 2 for proton and iron along with observed MMD, for 4.4 ≤ log(Ne ) <4.6. The distributions of proton and iron are scaled such that the tails of distributions overlap
with the observed MMD. The low and high multiplicity of observed MMD is well described by
the proton and iron, respectively. However, primaries of intermediate mass group are required
to describe the middle range of the observed MMD. The unfolding method is used to extract
the relative composition of the simulated primaries group for each shower size bin, using theGold’s algorithm. The observed muon multiplicity vector (µ) contains the observed MMDs andthe response matrix R contains the probability values such that R i, j represents the probability

031.4

log10 (Ne): 4.4 - 4.6

3

10−1
2
10−2

Normalized shower rate

probabilityloge(µ)SciPost Phys. Proc. 13, 031 (2023)

log10 (Ne): 4.4 - 4.6

10−1

data
Resultprotonheliumnitrogenaluminium
iron

10−2
10−3
10−4

1
10−5

0

H

He

N

Al

Fe

0

1

2

3
4
Mass index (A)

−3

10

10−6
0

10

20

30
40
Muon multiplicity

Figure 3: Left: Response matrix R for the MMD for the EAS induced by all simulated primaries. Right: Comparision of the resultant MMD with the observed
MMD after scaling the simulated MMD with the relative compositional weights;
for 4.4 ≤ log(Ne ) < 4.6.

Relative abundance (%)

of an EAS initiated by the j th simulated primary having the muon multiplicity value µi . The
response matrix for 4.4 ≤ log(Ne ) < 4.6 is shown in the left panel of Figure 3, where the color
gradient represents the probability values. The Gold algorithm is an iterative method and therelative composition vector for (k + 1) th iteration (Ak+1 ) is estimated from the estimate of thek th iteration (Ak ) as,
Ak (C R) T (C µ)
Ak+1 =,(4)(C R) T (C R)Ak
p
where C is the error matrix for the observed data such that Ci, j = δi, j / µi . The relative
composition assumed by the H4a models for the all simulated primaries used as the priorfor the unfolding. The algorithm tries to minimize the χk2 for the muon multiplicity vector((C R) T C µ) and the forward unfolded muon multiplicity vector ((C R) T (C R)Ak ). The number
2of iterations are stoped where the chi-square improvement (∆χk2 = χk2 − χk−1
) becomes less−3
than 10 . The result of the unfolding is shown in the right panel of the Figure 3 for 4.4 ≤log(Ne ) < 4.6, where each NBD curve corresonding to the simulated primaries is scaled with
its relative composition and the resulting MMD exhibiting a good agreement with the observed
MMD.
70

Preliminary

proton

60

50

40
4.0

4.5

5.0

5.5

6.0
log10 (Ne)

Figure 4: Relative composition of proton as a function of shower size.

031.5

SciPost Phys. Proc. 13, 031 (2023)

5

Results

Flux x E2.7 [m2 sr s GeV-1.7]-1

Flux x E2.7 [m2 sr s GeV-1.7]-1

The preliminary results for the measurements of relative composition and energy spectrum for
proton component are presented. The relative abundance of the proton component decreases
from ∼65% at shower size 104.1 to ∼42% at shower size 106.1 and is shown in the left panelof Figure 4. In order to get the size distribution for the proton component, the fraction of thenumber of showers is selected from the observed size distribution based on the relative abundance of the proton component for each shower size bin in the shower size range from 104.0
to 106.2 . Gold’s algorithm is used to unfold the proton energy distribution from corresponding
size distribution, assuming the H4a proton spectrum as the prior. The differential cosmic-ray
spectrum (dI/dE) can be expressed as follows,
Preliminary

104

103
103

ATIC-2 (2009)

KASCADE (2005, QGSJET-I-01)

CREAM-I+III (2017)

KASCADE (2005, SIBYLL-2.1)

CALET (2019)

KASCADE (2013, QGSJET-II-02)

NUCLEAON KLEM (2019)

ICETOP (2019, SIBYLL-2.1)

DAMPE (2019)

GRAPES-3 (This work, QGSJET-II-04)

104

105

Preliminary
Eb = 208.5 +/- 1.6 TeV
γ = -3.14 +/- 0.002
1

γ = -2.54 +/- 0.006
2

104

5

106
107
Energy [GeV]

10

6

10
Energy [GeV]

Figure 5: Left: Elemental spectrum of proton; compared with other experiments[4–6, 10–14]. Right: Elemental spectrum of proton fitted with a smoothly broken
power law (solid curve) showing an energy break at ∼210 TeV.
dI
1
=
dE
Tobs



N
∆E · Aacc



,

(5)

i

where N , Aacc and ∆E are the number of EAS, acceptance and width of i th energy bin, respectively, and To bs is the live time of the data. The measured preliminary energy spectrum of
proton plotted with direct [4–6, 10, 11] and indirect observations [12–14] is shown in the left
panel of Figure 5. The statistical error bars are smaller than the marker size. The flux of themeasured proton spectrum in this work is consistent with CREAM I+III at lower energy andis consistent with KASCADE (QGSJET-I-01) at higher energy. The measured proton spectrum
is fitted with a broken power law in an energy range from 50 TeV to 800 TeV as shown in theright panel of Figure 5. An energy break (E b ) is observed at nearly 210 TeV with the spectral
slope of −3.14 and −2.54 before (γ1 ) and after (γ2 ) the energy break.

Acknowledgements
We are grateful to D.B. Arjunan, A.S. Bosco, V. Jeyakumar, S. Kingston, N.K. Lokre, K. Manjunath, S. Murugapandian, S. Pandurangan, B. Rajesh, R. Ravi, V. Santoshkumar, S. Sathyaraj,
M.S. Shareef, C. Shobana, R. Sureshkumar for their role in efficient running of the experiment
Funding information We acknowledge support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI4002. This work was partially supported
by grants from Chubu University, ISEE of Nagoya University and ICRR of Tokyo University,
031.6

SciPost Phys. Proc. 13, 031 (2023)

Japan. Numerical computations were in part carried out on PC cluster at Center for Computational Astrophysics, NAOJ, Japan. The first author is thankful to the University Grant
Commission, New Delhi, for the support of his fellowship.

References[1] S. K. Gupta et al., GRAPES-3 – A high-density air shower array for studies on the structure in
the cosmic-ray energy spectrum near the knee, Nucl. Instrum. Methods Phys. Res. A: Accel.
Spectrom. Detect. Assoc. Equip. 540, 311 (2005), doi:10.1016/j.nima.2004.11.025.
[2] P. K. Mohanty et al., Measurement of some EAS properties using new scintillator
detectors developed for the GRAPES-3 experiment, Astropart. Phys. 31, 24 (2009),
doi:10.1016/j.astropartphys.2008.11.004.
[3] P. K. Mohanty et al., Monte Carlo code G3sim for simulation of plastic scintillator detectors with wavelength shifter fiber readout, Rev. Sci. Instrum. 83, 043301 (2012),
doi:10.1063/1.3698089.
[4] Q. An et al., Measurement of the cosmic ray proton spectrum from 40 GeV to 100 TeV with
the DAMPE satellite, Sci. Adv. 5, eaax3793 (2019), doi:10.1126/sciadv.aax3793.
[5] E. V. Atkin et al., Energy spectra of cosmic-ray protons and nuclei measured in
the NUCLEON experiment using a new method, Astronomy Rep. 63, 66 (2019),
doi:10.1134/S1063772919010013.
[6] Y. S. Yoon et al., Proton and helium spectra from the CREAM-III flight, Astrophys. J. 839,
5 (2017), doi:10.3847/1538-4357/aa68e4.
[7] Y. Hayashi et al., A large area muon tracking detector for ultra-high energy cosmic ray
astrophysics – The GRAPES-3 experiment, Nucl. Instrum. Methods Phys. Res. A: Accel.
Spectrom. Detect. Assoc. Equip. 545, 643 (2005), doi:10.1016/j.nima.2005.02.020.
[8] H. Tanaka et al., Studies of the energy spectrum and composition of the primary cosmic rays
at 100-TeV - 1000-TeV from the GRAPES-3 experiment, J. Phys. G: Nucl. Part. Phys. 39,
025201 (2012), doi:10.1088/0954-3899/39/2/025201.
[9] V. B. Jhansi et al., The angular resolution of GRAPES-3 EAS array after improved timing
and shower front curvature correction based on age and size, J. Cosmol. Astropart. Phys.
24 (2020), doi:10.1088/1475-7516/2020/07/024.
[10] O. Adriani et al., Direct measurement of the cosmic-ray proton spectrum from 50 GeV to 10
TeV with the calorimetric electron telescope on the international space station, Phys. Rev.
Lett. 122 (2019) 181102, doi:10.1103/PhysRevLett.122.181102.
[11] A. D. Panov et al., Energy spectra of abundant nuclei of primary cosmic rays from thedata of ATIC-2 experiment: Final results, Bull. Russ. Acad. Sci. Phys. 73, 564 (2009),
doi:10.3103/S1062873809050098.
[12] M. G. Aartsen et al., Cosmic ray spectrum and composition from PeV to EeV using 3 years of data from IceTop and IceCube, Phys. Rev. D 100, 082002 (2019),
doi:10.1103/PhysRevD.100.082002.
[13] W. D. Apel et al., KASCADE-Grande measurements of energy spectra for elemental groups ofcosmic rays, Astropart. Phys. 47, 54 (2013), doi:10.1016/j.astropartphys.2013.06.004.
031.7

SciPost Phys. Proc. 13, 031 (2023)[14] T. Antoni et al., KASCADE measurements of energy spectra for elemental groupsof cosmic rays: Results and open problems, Astropart. Phys. 24, 1 (2005),
doi:10.1016/j.astropartphys.2005.04.001.

031.8