29th International Cosmic Ray Conference Pune (2005) 10, 243-254

GRAPES-3 Experiment
Y. Hayashi on behalf of the GRAPES Collaboration
Graduate School of Science, Osaka City University, Osaka 558-8585, Japan

Presenter: Yoshio Hayashi (hayashi@sci.osaka-cu.ac.jp)
It is important to study nuclear composition and energy spectrum near knee region and search for highenergy γ-rays sources. For this purpose, we designed GRAPES-3 as high density air shower array with large
area muon detector. We study short and long term cosmic ray variation also with large area muon detector.
Here we report following results : 1) the mean mass number gradually increases through the knee region andpredictions of SIBYLL agree with JACEE results, but some discrepancy is seen between QGSJET and
JACEE, 2) we found clear trend to observe Loss-Cone precursor decrease with greater frequency before alarge Forbush decrease and 3) Our analysis shows the existence of Tail-IN sidereal anisotropy in 5 years data.
Addition to them we report new detector development also.

1. Introduction
Primary cosmic ray energy spectrum index changes around knee region, 3×1015eV, from ~-2.7 to ~-3.1. It is
important to know precise data when we think acceleration and propagation mechanism of cosmic rays. This
information will give restriction to models. Direct measurement such as satellite based experiment at knee
region has difficulty due to low flux. Extensive air shower experiment is only practical way for them. Many
air shower experiments were carried out so far. However, conclusion is not yet given.
Sun is the only one object that we can observe directly and precisely in the cosmos. Extrapolation of themodels for cosmic ray particle acceleration and propagation at sun to very high energy has meaning. For thispurpose, we can observe long and short-term variation. Mainly we can get acceleration information fromshort-term variation and propagation from long-term variation.
GRAPES-3 experiment has been designed to have high density scintillation detectors and large area muon
detectors. High density scintillation detectors allow to get precise data of electro-magnetic component toestimate primary cosmic ray energy and direction. We estimate ratio of nuclear species by comparison
observation and Monte Carlo simulation of relation between shower size and number of muons. This data isuseful to reject back ground from high energy gamma-ray data. Large area muon detectors are used formeasurement of cosmic ray variation. Rigidity cut off of primary cosmic ray is about 17GV for verticalmuon and median energy of primary proton is estimated as about 60GeV.
Monte Carlo simulation is used to reconstruct the shower in the atmosphere and estimate properties ofprimary
We reported diffuse gamma ray results at previous conference [1]. Here we report result for nuclear
composition, short and long-term variation of cosmic ray.

2. Detectors
The GRAPES-3 air shower array is located at Ooty in southern India (2200m altitude, 11.4ºN, 76.7ºE). It
consists of 288 scintillation detectors to measure the electro-magnetic component and about 3700
proportional counters to detect muon component in showers and low energy cosmic rays. In Fig. 1 aschematic layout of the GRAPES-3 array is shown. It is desirable to make reliable measurements on theenergy spectrum at lower energies, around 1013-1014eV using the EAS technique in order to have an overlap

244

Yoshio Hayashi

with direct measurements and to understand the systematic. A good agreement between direct and EAS
measurements at these energies would provide required confidence in measurements at UHE. Further,measurements over entire energy range of interest, 3×1013-3×1016eV should be made with same array tominimize possible systematic errors. Fig. 1 also shows 16 squares in the north-west region (Fig 2). Each ofthese squares represents a 4-layers muon-tracking detector with an energy threshold of 1 GeV for verticalmuons. Each layer consists of 58 proportional counters, each 6m long with 10×10 cm2 cross sectional area.
This large area, muon detector has been designed for studies on the composition of primary cosmic rays andstudies on cosmic sources of UHE gamma rays. Large area muon detector is also required for observations
on muons produced by lower energy protons that are affected by phenomena occurring on the Sun, such as
solar flares and coronal mass ejections leading to magnetic disturbances around the Earth. The arrival
direction of detected muon is categorized into 15×15 directional cells. The total area of 16 muon telescopes
is 560 m2. Statistical accuracy of hourly count in each cell is 0.1% ~ 0.3%[2][3].

Figure 1. Schematic layout of GRAPES-3 air shower detectors.
Small square shows 1m2 scintillation detectors.
Large square shows 35m2 muon detectors each.

Figure 2. Full view of GRAPES-3 Air Shower ArrayGRAPES-3 Experiment

245

Figure 3. Schematic of a muon detector module showing 4 layers of proportional countersembedded within concrete blocks.

3. Nuclear composition
A total of 6×108 showers observed during two years of run in 2000 and 2001 (about 560 live days) havebeen analyzed. Each shower parameter, i.e. shower size, age parameter and core location, are estimated fromlateral distribution of particles with a Maximum Likelihood Method fit to NKG function. Arrival direction of
shower is estimated from the timing data. The muon tracks are reconstructed using tracking muon detectors
to agree with shower direction and arrival time. They help to reduce the background of muon. Reconstructed
showers are selected to be within a fiducial area and 80 m away from center of muon detectors, which isshown as yellow (or gray) area in Fig. 4, to reject showers that have unreliable information due to situatedoutside or edge of the array and to prevent saturation of muon detectors. In Fig. 5, average number ofdetected muon is shown for a given shower size. CORSIKA MC(v6.02)[4] results are also shown for proton
and iron nuclei. The Ne-Nμ relation features the primary nuclear composition [5]. Results bases on hadronicinteraction models, QGSJET and SIBYLL, are presented here. The energy spectrum of each component canbe extracted from the shapes of muon multiplicity distributions. Fig. 6 shows the muon multiplicity
distributions from observation and MC’s after fitting them to the observations to estimate relative abundancefor each component for Ne range 6.0<log(Ne)<6.2 as a sample. Here, the ratio of abundance of aluminum toiron is fixed to 0.8, because it is difficult to separate them. This value was decided from extrapolation of low
energy observations by RUNJOB. The relations between shower size and primary energy are calculated by
MC for each mass groups and the primary energies are converted from shower size as explained below.

Figure 4. Accepted area of
shower core location shownin yellow (or gray) area islocated inside the array but>80 m away from center ofmuon detectors

Figure 5. Relation between shower
size and average number of muons.
Red and blue lines show MC resultsfor proton and iron primaries.
Dependenceonthehadronic
interaction models is represent-ed by
QGSJET (fat lines) and SIBYLL (finelines)

Figure 6. Muon multiplicity distributionsclassi-fied by shower size. Observations
are shown black. MC (QGSJET) results for
proton, helium, nitrogen, aluminum andiron are shown in red, green, light blue,
violet and blue respectively. MC distributions are fitted to observations toestimate relative abundance for each.

246

Yoshio Hayashi

The MC is carried out using E-2.7 spectrum and then the MC showers are subjected to same cuts asobservational data. They are binned according to shower size and then the mean energy < log(E) > for eachbin is calculated.
The energy spectra and mean mass values derived from GRAPES-3 data are shown and compared with otherobservations in Fig. 7, 8, 9 and 10.

Figure 7. A comparison proton and helium energy
spectra with other observations, RUNJOB, Tibet,
EAS-TOP etc. from Kobayakawa et al.[6].

Figure 8. A comparison of energy spectra of CNO,
medium heavy and iron groups with other observations.
The flux of CNO and Fe-groups is multiplied by 10 and
1/100, respectively [6].

Figure 9. A comparison of all-particle energy
spectra.
GRAPES-3
comparedwith
RUNJOB, Tibet, Akeno, CASA-DICE,CASA-MIA, CASA-BLANCA etc. [6].

Figure 10. A comparison of mean mass (<lnA>) for
GRAPES-3 with JACEE, RUNJOB, CASA-DICE,CASAMIA, CASA-BLANCA, EAS-TOP, KASCADE,
Chacaltaya etc. [6].

In Fig. 9, the effect of interaction models is relatively small and flux from QGSJET is only 20% higher thanSIBYLL result. No disagreement with other experiments is found for both models with all-particle spectra.
However, there are significant differences in the flux for each component. In Fig. 7 and 8, QGSJET flux ishigher for proton, and SIBYLL is higher in heavier component. Since direct observations are good checks on
the interaction model, a comparison of direct results from JACEE with QGSJET and SIBYLL shows good
agreement with SIBYLL and some disagreement with QGSJET results.

GRAPES-3 Experiment

247

A comparison of the mean masses is shown in Fig. 10. SIBYLL has 0.7 ~ 0.8 larger <lnA> than QGSJET.
Both have tendency to increase <lnA> above 106 GeV. Here, SIBYLL has good agreement with direct
measurements, too. Some models are suggested to reproduce knee in the energy spectrum. In Fig. 11, plotsof mean mass number that are expected from various knee models are shown [7] and results from SIBYLL,
which is one of the most reliable model of hadronic interactions are also shown. Results from the following
models of knee, (A) Acceleration in supernova remnants, (B) Acceleration by supernova shocks, (C) The
single-source model, (D) The minimum-path-length model, (E) Drift in the global regular magnetic field of
the Galaxy and (F) Diffusion and drift are shown in Fig. 11. (A) ~ (C) expect knee appears due toacceleration mechanisms, (D) ~ (F) it is due to propagation and leakage. In these, similar results are obtained
from “(F) Diffusion and drift” model by comparison of energy spectrum for each mass component. In Fig.
12, CNO flux is 1/3 smaller than SIBYLL, though it has good agreement in proton and helium spectra. In theabove, there is no model that agrees with our results completely, although “Diffusion and drift” model has
almost same results as our data.

Figure 11. A comparison of GRAPES-3 results withSIBYLL using six models of the knee from J.R.
Hoerandel[7].

Figure 12. Diffusion and drift (F) knee model by
E. Roulet[8] and GRAPES-3 results using
SIBYLL.

4. The Loss-Cone Precursor Decrease (LCPD)
From 1 March 2001 to 31 December 2002, a total 659 days of data were analyzed. Forbush decrease (FD)events were selected from this dataset. For the analysis, FOV was divided into 9 directional cells as shown in
Fig. 13. From hourly count rate data for these 9 cells, the FD events were searched for. FD candidate wereselected under the condition that any one of the 9 cells showed a deficit greater than 0.4% between minimum
intensity of current 24 hours and the up-coming 24 hours period. The start-time of FD is defined by the time
when intensity decreased lower than lowest intensity in current 24 hours. The amplitude of FD is defined by
maximum deficit among the data of 9 cells. The empty histogram in black, in Fig 15 corresponds to an FD
candidate as described above. Some of those events with amplitude smaller than 0.6% were rejected during
this analysis due to a gradual slope.

248

Yoshio Hayashi

Figure 13. 15×15 cells in FOV. The FOV wasreconfigured into 9 (NW,N,NE,W,V,E,SW,S,SE)
also for this analysis.

Figure 15. Distribution of the amplitude ofobserved Forbush decreases.

Figure 14. Asymptotic Viewing direction at themedian rigidity of muon telescope in geographicalcoordinate.

Figure 16. An example of subtraction of backgroundanisotropy. (2001/3/25 – 26)

It is already known from earlier works [9] that some of the LCPD anisotropy seems to have a relatively
broad 2-dimension structure. In such a case, the limited FOV of the telescope does not permit an accurate
measurement of amplitude of the anisotropy. To evaluate amplitude correctly, it is necessary to estimate thesolar diurnal anisotropy. For this purpose the data are separated into units of 24 hours duration (start from 5h
in UTC) and each unit was categorized into one of the three types “A”, ”B” and “N” using the following
criteria.
(A): Data unit include start of FD and following two data units were defined as category “A”(B): Data unit includes period before the start of the FD. Two data units prior to the units that include start of
FD were defined as category “B”. If a unit is included in “A” also, the unit is categorized into “B”.
(N): Data units those were not included neither “A” nor “B”. Number of data sets in each category is 136
units in “A”, 84 units in “B”, 27 units in “B” and 412 units in “N” data.
We have used the central 169 out of 225 directional cells as the FOV. Loss-Cone (LC) search were done forall category of data. In this analysis, relative intensity in the data unit is used hereafter. The observing
direction of each cell is defined as asymptotic arrival direction of median rigidity particles. Shape of large
scale anisotropy such as solar diurnal anisotropy or Swinson flow can be estimated using data points in the
direction of GSE longitude 0º-270º. The estimations were done by fitting a spherical harmonic function tothe observed data. Up to second order harmonics and linear function of time were used for fitting. The
second harmonics is to reduce the residuals due to approximate knowledge of the observing direction.

GRAPES-3 Experiment

249

Common variation was assumed as a linear function of time. In the case of data in category “A”, a total of
17 sets could not be fitted due to complicated intensity variation caused by FD. After obtaining the shape of
diurnal anisotropy, residual between best fit function and the data are calculated for each data point. At each
5º, mean residuals are calculated using data points within the opening angle of 5º.
The top two panel of Fig. 16 are an example of observed anisotropy and best of fit large scale anisotropy
estimated by the method described above. Bottom two panels on Fig. 16 show a residual after subtracting
best of fit anisotropy from observed data. Bottom right panel shows the histogram of residual values in the
direction of GSE 270º ~ 360º. From eye-scan analysis, it is known that for LCPD events, the histograms
have significant distribution below -0.35%. Thus in this analysis we defined the integrated value of thishistogram below -0.35% as Σd 0.35 . Using this Σd 0.35 parameter, it becomes possible to evaluate local
anisotropy near IMF more objectively even with data from a single station. Following table lists frequencyof observed Σd 0.35 value in each category.
Table 1: Frequency of observing Σd 0.35 under several condition.

Figure 17. Frequency of Fds which follows LCPD evaluated by Σd 0.35
parameter as displayed in the same histogram.

Figure 18. Relation between FD
amplitude
andamplitudeof
LCPD( Σd 0.35 <-10.0)

During 13 days which show Σd 0.35 < −10.0 were found from 84 days categorized “B”. With same conditions
13 days were found from category “N”. Thus it is clear that days showing large Σd 0.35 value are concentrated
before an FD rather than “N”(Normal) days. The chance probability to get this result is less than 10-5. By
combining 4 events observed in “B” data, the correlation between FD amplitude and LCPD were obtained as
Fig. 17. In Fig. 17, one can see a clear trend to observe LCPD with greater frequency before a large FD.
These results which were expected from LCPD model proposed by Nagashima et.al [10][11]. Fig. 18 shows
correlation between maximum of LC depth and following FD amplitude. This trend is also expected for themodel in [10][11].

250

Yoshio Hayashi

5. Sidereal anisotropy of Cosmic rays
Using data of five years (2000-2004) from GRAPES-3 muon telescopes, we obtained the sidereal time
variation of primary cosmic rays. Our muon telescopes are able to measure the incidence direction of thecosmic ray with rather good accuracy (7 degrees). The data were classified according to the direction(polarity) (Toward and Away) of IMF. The polarity of each day was determined from the satellite data of
ACE level 2. Calculating the 72 hours running average we adopted the value of 22h UT for the polarity of
the corresponding day. The muon telescope consists of a total of 16 modules (each module has an effective
area of 35 m2, total area of 560 m2) resulting in smallest statistical error among all instruments operating in
the world. The data for each direction is summed for one hour and the resulting counting rate is normalizedby average value for the day and also corrected for barometric pressure [3]. In addition, in order to remove
various systematic errors, we adopted the so call E-W method (refer to [12] for details) to obtain the siderealtime variation [13]. For a median energy of 60 GeV the sidereal time variation, there would becontamination due to the daily variation of pseudo sidereal time anisotropy called Swinson flow [14]. Thisanisotropy arises due to the gradient of cosmic rays in the vicinity of the Earth and the polarity (Toward,
Away) of IMF. The amplitude (about 0.07%) of anisotropy due to Swinson flow is larger than the expected
genuine sidereal time anisotropy (about 0.03%). The average anisotropy by Swinson flow should show thesame variation with a phase difference of 12 hours with respect to the polarity of IMF. So, if we take averageof the intensity for both polarities (T+A)/2, then pseudo anisotropy due to Swinson flow can be cancelled
and only a genuine anisotropy would be left in the data. In case of (T-A)/2 this effect of the original Swinson
flow would be enhanced by a factor of 2. These two sets of data namely, (T-A)/2 and (T+A)/2 observed by
the GRAPES-3 muon telescopes, are shown in Fig. 19 and Fig. 20 for 5 years from 2000 to 2004. Each panel
shows the sidereal time variation in an individual year. As one can see in these two figures that for each year
the value fluctuates considerably. In the case of (T-A)/2 “Swinson flow” intensity variation should follow asinusoidal wave with a maximum at a phase at 18 hours and a minimum at 6 hours for a long-term time
average. The daily variation for each year is distorted in due to effect of “Swinson flow” in the case of single
year data. It implies that single year data cannot be used to judge the genuine sidereal anisotropy.

Figure 19. (T-A)/2 data for each year

Figure 20. (T+A)/2 data for each year

Therefore, the average of five years of data is taken and is shown for (T-A)/2 in Fig. 21 and (T+A)/2 in Fig.
22. As can be seen in Fig. 19, (T-A)/2, the phase of the minimum appears around 6 hours and the maximum
near 18 hours. Thus averaging over five years has resulted in effect of Swinson flow being smoothed into asine wave and the effect of pseudo sidereal anisotropy could be largely cancelled. Same thing is true for the
genuine sidereal anisotropy as well. As seen in Fig. 20 for (T+A)/2, there exist a clear excess at 6 hours, but
no other clear excess or deficit can be seen in present analysis.

GRAPES-3 Experiment

Figure 21. Averages of five years data for(T-A)/2 for north(N), vertical(V) andsouth(S) directions

251

Figure 22. Averages of five years data for(T-A)/2 for north(N), vertical(V) andsouth(S) directions

Figure 23. Sidereal anisotropy. Result is matching with Tail-in excess for phase and width.

Present analysis shows the existence of Tail-IN excess, but the Loss-cone deficit is not clearly seen. Since
three years out of total of five years observation period happen to fall in the near maximum phase of solar
activity, a strong disturbance due to this solar activity might have been hiding the Loss-cone type anisotropy.
The energy of the cosmic rays connected with Tail-IN is considered to be lower than Loss-cone. The factthat the anisotropy is seen only at low energy in present analysis is somewhat difficult to understand. One ofthe possible explanation for this kind of Tail-IN anisotropy could be the difference between the “Nose” and“Tail” structure of the heliosphere which exists considerably deeper in the heliosphere than usually expected.

252

Yoshio Hayashi

6. Further development
1.0E+08
1.0E+07

H igh G ain P M T
Low G ain P M T

1.0E+06
1.0E+05
1.0E+04
1.0E+03
1.0E+02
1.0E+01
0.1

1

10

100

1000

10000

N um ber of P articles

Figure 24. Density Spectrum of Dual PMT detector. Dynamic range extended to more than 1000 particles

We developed new type of scintillation detectors and circuit modules that have better performance for future
extension of air shower array. Here we show 1) Dual PMT scintillation detectors with wave length shifting
fiber[15] and 2) TDC module with CERN TDC32 chip[16] some of them.
1) Dual PMT scintillation detectors with wave length shifting fiber
We made dual PMT detector to avoid PMT output saturation. New detector configuration is followed. a) 4
blocks of plastic scintillators, each 2 cm thick and 50×50 cm2 in area, b) 12 straight grooves for each blocks,c) 18 wave length shifting fibers for high gain PMT, d) 6 wave length shifting fibers for low gain PMT. Eachfiber length is 1.3m. We could improve compare with old cone type scintillation detectors, photon yield
increased by about factor 2, position response uniformity becomes within 10% RMS include PMT responseand dynamic range. Fig.24 shows relation between high and low gain PMT density distribution. Linear range
achieves more than 1000 particles.

Figure 25. Commercial TDC module (Phillips Scientifics TDC 7186) linearity check result.

GRAPES-3 Experiment

253

Figure 26. Home made TDC module with CERN TDC32 chip linearity check result.

2) TDC module with CERN TDC32 chip
New TDC module can handle 31 channels and wide dynamic range. It has some advantages compare withcommercial one, i.e. wider dynamic rage, better linearity, more freedom for timing between start and stop
signal. Fig. 25 and 26 shows sample of commercial TDC module and home made module linearity checkdata. Commercial one has little deviation from linear.

4. Acknowledgements
We thank D.B. Arjunan, A.A. Basha, G.P. Francis, I.M. Haroon, V. Jeyakumar, K. Manjunath, B. Rajesh, K.
Ramadass, C. Ravindran and V. Viswanathan for help in construction and operation of GRAPES-3
experiment. Japanese members acknowledge partial financial support from Ministry of Education,
Government of Japan. We thank N.K. Mondal and colleagues for loan of proportional counters.

References[1] Y. Hayashi et al. Proc. 28th ICRC, 4, 2273, Tsukuba (2003).
[2] S.K. Gupta et al. Nucl. Instr. & Meth. A, 540, 311 (2005).
[3] Y. Hayashi et al. Nucl. Instr. & Meth. A 545, 643 (2005).
[4] D. Heck et al. Forschungszentrum Karlsruhe, FZKA 6019 (1998).
[5] H. Tanaka et al. Proc. 28th ICRC, 1, 155, Tsukuba (2003).
[6] K. Kobayakawa et al. Phys. Rev. D, 66(8), 083004 (2002).
[7] J.R. Hoerandel. Astropart. Phys., 21, 241 (2004).
[8] E. Roulet. Int. J. Mod. Phys. A, 19(7), 1133 (2004).
[9] T. Nonaka et.al.,Proc.28th ICRC, 6, 3569, Tsukuba (2003)[10] A. Belov and J.W Bieber et al., Proc. 27th ICRC, Hamburg (2001)[11] K. Munakata et.al., J. Geophys. Res A 105, 27457 (2000)[12] K. Nagashima et al., Nuovo Cimento, C(6), 550 (1983)[13] H. Kojima et al., Proc 28th ICRC, 7, 3957, Tsukuba (2003)[14] JD.B. Swinson, J. Geophys. Res., 745591 (1969)[15] PSF catalog, Kuraray[16] J. Christiansen, CERN/ECP – MIC