Recent results from the GRAPES-3 experiment

The GRAPES-3 is a dense extensive air shower array operating with ∼400 scintillator detectors and a 560 m2 large tracking muon detector (E μ >1 GeV), at Ooty in India. The muon detector has been used to observe acceleration of muons during thunderstorm conditions. The muon multiplicity distribution of the EAS is used to probe the composition of primary cosmic rays below 1 PeV, with an overlap with direct measurements. More recently we have explored the possibility of using the angular distribution of >1 GeV muons to identify the best from among several lowand high-energy hadronic interaction models. We have major expansion plans to enhance the sensitivity of the GRAPES-3 experiment in all of the areas listed above.


Introduction
The GRAPES-3 (Gamma Ray Astronomy at PeV EnergieS-phase 3) is a high density extensive air shower (EAS) array designed for precision study of acceleration of charged particles in nature using a number of phenomena as detailed below.The experiment was started with 256 plastic scintillator detectors (each 1 m 2 in area) in south India in 2001 as shown in Fig. 1.The array has now been expanded to about 400 detectors that are deployed in a dense hexagonal configuration with an inter-detector separation of 8 m, at Ooty (2200 m altitude, 11.4 • N, 76.7 • E).About 100 scintillator detectors are instrumented with two fast photomultiplier tubes (PMTs) for extending the dynamic range to ∼5×10 3 particles m −2 .The scintillators, signal processing electronics and data recording systems were fabricated in-house to minimize costs and optimize performance [1].
The array also contains a large area (560 m 2 ) tracking muon detector [2] to measure the cosmic ray induced muon flux with high statistical accuracy.This high statistical accuracy has allowed us to probe the acceleration of muons in the atmospheric electric fields set up during thunderstorm activity that is also known to result in terrestrial flashes of γ-rays [3].The muon detector has also provided new information on solar flares, coronal mass ejections (CMEs) and Forbush decrease (Fd) events [4,5] besides serving as a sensitive probe of the interplanetary magnetic field and the turbulence caused in it due to the propagation of the CMEs [6].
It has been shown earlier by us that the cosmic ray energy spectrum and its composition can be extracted using the muon multiplicity distribution (MMD) and the shower size of detected EAS [7].Using the GRAPES-3 data, this a e-mail: gupta@grapes.tifr.res.intechnique was used to extract the composition of primary cosmic rays from 30 TeV to 1 PeV, allowing a sizable overlap with direct measurements using balloon and satellite based detectors.However, this interpretation of EAS data involves extrapolation of hadronic interaction properties to energies beyond those available from accelerators.This necessitates modeling of these interaction assuming certain reasonable assumptions that lead to several models of hadronic interactions depending on the nature of assumptions.Therefore, the interpretation of the EAS data in terms of composition and energy spectrum gets coupled to the choice of hadronic interaction models used.In GRAPES-3, we used the overlap in energy with the direct measurements to identify hadronic interaction models that had yielded composition and energy spectra consistent with direct measurements in the overlap energy region from 300 TeV to 1 PeV [8].
However, it will be enormously helpful if the nature of hadronic interactions could be ascertained independent of the composition.That would allow the cosmic ray composition to be later on independently extracted.Use of the angular distribution of muons detected by GRAPES-3 experiment offers one such possibility.The GRAPES-3 muon detector measures muons at a rate of 4×10 9 d −1 from a solid angle of 2.7 sr, allowing precise measurement (∼0.1%) of muon angular distribution to be made.Extensive Monte Carlo simulations using CORSIKA [9] were done to obtain muon angular distribution for GRAPES-3 experiment for all possible combinations of, three low-and three high-energy hadronic interaction models.Low energy models included GHEISHA [10], UrQMD [11] and FLUKA [12], while QGSJet-II [13], SYBILL 2.1 [14] and EPOS 1.99 [15] were high energy models used.
While simulations may demonstrate the feasibility of using muon angular distribution to probe the low-and high-energy hadronic interaction models, the final arbiter would be an actual comparison of experimentally measured muon angular distribution with predictions from simulations using different combinations of low-and highenergy hadronic interaction models.Since the expected difference in muon angular distributions are at the level of 1%, a major challenge is to measure the efficiency of muon detector to a precision better than 0.1%.We have developed a new numerical technique to model the slow variation in the gain of proportional counters that comprise the GRAPES-3 muon detector.This technique as discussed below in §2 results in measuring efficiency of muon detector modules to a precision ∼0.01%over a period of 100 days.Once the muon detector configuration has been modeled using GEANT4, we hope to be able to extract the appropriate hadronic interaction models from data in near future.

Experimental Details
To operate the GRAPES-3 EAS array at the lowest possible threshold energy, a simple 3-line coincidence of detectors is used to generate the Level-0 trigger, that acts as the fast GATE and START for the charge to digital (ADCs) and time to digital converters (TDCs), respectively.Since this trigger selects very small local showers and also larger showers whose cores land far from the area covered by the array, we require at least 12 out of the inner 127 detectors should have triggered their discriminators within 1 µs of the Level-0 trigger.This Level-1 trigger with an observed EAS rate of 15 Hz is used to record the charge (ADC) and the arrival time (TDC) of the pulses from each detector [1].With the expansion of the array the trigger rate has now increased to 35 Hz.The pulse charge is later converted into equivalent number of minimum-ionizing particles (MIPS) using the most probable charge for a single MIP measured.The high sensitivity of scintillator detectors allows the at-mospheric radon decay products to be routinely detected during every episode of rain in Ooty.In Fig. 1 the 16 squares shown in lower left side, each represents a 4-layer tracking muon detector module with an energy threshold of 1 GeV for vertical muons.Each layer consists of a total of 58 proportional counters (PRCs), each 6 m long with 10 × 10 cm 2 cross-sectional area.The 560 m 2 muon detector consists of 16 modules as displayed in Fig. 1.The muon energy threshold of 1 GeV is achieved by placing a concrete absorber of thickness 550 g.cm −2 above the bottom layer of PRCs.Thus each detector module of a sensitive area of 35 m 2 consists of 232 PRCs arranged in four layers of 58 PRCs each, with alternate layers placed in orthogonal directions.Two successive layers of PRCs are separated by a 15 cm thick concrete layer, consisting of 60 × 60 × 15 cm 3 blocks.Fourlayer PRC configuration of the muon modules allows a 3-D reconstruction of muon track directions to an accuracy of ∼6 • .The 16 modules detect ≥1 GeV muons at a rate of 4×10 9 d −1 , allowing tiny variations in rate due to atmospheric, solar and other systematic effects to be precisely measured.Using the reconstructed muon direction the muons are segmented into nine solid angle bins, labeled vertical (V), north (N), south (S), east (E), west (W), south-east (SE), north-east (NE), south-west (SW), northwest (NW) covering a solid angle of 2.7 sr.
Daily rate of muons measured for a period of 100 days (1 Jan-10 Apr 2006) by modules 4 and 5 is shown in Fig. 2a, Fig. 2b, respectively.The rates are nearly same and display similar variation caused by atmospheric and solar effects.However, the ratio of their rates shown in Fig. 2c displays a much smaller (<0.4%), but systematic variation that represents the ratio of efficiencies of modules 4 and 5.This is caused by a gradual change in muon detection efficiency of PRCs that is modeled using a fourth order polynomial in time.Since GRAPES-3 muon detector consists of 16 modules, a total of 120 independent combinations of efficiency ratios can be formed using corresponding fourth order polynomial for each of the 100 days.The 120 series of 100 ratios represents the ratio of unique 04005-p.2efficiency profile of the corresponding two modules.Since the magnitude of variation is small (<1.0%) this problem can be solved after performing Taylor series expansion and using singular value decomposition (SVD) technique that yields the five coe ffi cients of the fourth order polynomials that describe the efficiencies of each of the 16 modules.The details of this technique will be described in a separate work under preparation.The efficiency ratio for modules 4 and 5 obtained from SVD solution yields a fit that describes the data rather well as seen from Fig. 2c.When the muon rates for modules 4 and 5 are corrected by using corresponding efficiency, resultant ratio shown in Fig. 2d almost equals unity as expected for an ideal detector.
In Fig. 3 daily rates of all 16 modules are shown for the 100 days interval after correcting for efficiency using the fourth order polynomial correction.Although, muon rates are displayed over a narrow range of only 20 (3020-3040) counts, yet the 16 rates are almost indistinguishable.The rms spread in daily muon rates is only 0.32 comparable to the statistical error of 0.26 counts.Thus we conclude that the efficiency correction was effective and an error of 0.32 counts in a rate of 3000 counts implies a precision of 0.01% in measurement of efficiency.

Atmospheric Acceleration
During thunderstorm activity, the muon rate displays a sudden reduction and recovery over time-scales of about 20-40 minutes.However, this reduction does not occur uniformly but along a narrow solid angle.In Fig. 4 the percent variation in rate of muons in nine directions is shown during a thunderstorm on 18 April 2011 during 12 successive 10 minutes exposures.A significant deficit appears in panels 5 and 6 for about 20 minutes.The episodes of reduced muon flux are accompanied by presence of large electric fields on ground and its polarity indicates that the atmosphere is at a large negative potential relative to the ground.Variation in the rate of muons from 14 to 18 h IST on 18 April 2011 in SE out of 9 directions is shown in Fig. 5a.A sharp reduction in muon the rate is clearly visible.The electric field at a location 6 km east of GRAPES-3 too shows a positive electric field on ground concurrent with the deficit in muon rate also lasting about 20 minutes is seen from Fig. 5b.It is well known that the ratio of positive to negative muons is ∼1.2 at a few GeV and gradually increases to ∼1.4 at 1 TeV [16][17][18].Thus a negatively charged cloud in atmosphere would cause the ground to appear at a relatively positive potential resulting in deceleration of positive muons, causing a drop in their detected rate.On the other hand acceleration of negative muons would result in an increase in their detected rate.Since the GRAPES-3 muon detector operates without magnetic field, it records the sum of both types of muons.However, as the positive muons outnumber negative ones, a net decrease would be observed in the total rate of detected muons by GRAPES-3 exactly as seen during the thunderstorm episode on 18 April 2011.

Cosmic Ray Composition Below 1 PeV
The composition of primary cosmic rays below 1 PeV has been measured using the MMD as explained in §1.These measurements have provided evidence of a steeper proton spectrum, consistent with direct measurements from balloon and satellite experiments.GRAPES-3 spectra for He, N, Al, Fe agree well with direct results, extending their energy and confirming their harder nature relative to protons as described elsewhere [8].The consistency of GRAPES-3 results with direct measurements require SYBILL to be the suitable model of hadronic interactions in overlap energy of 100-300 TeV.The extension of energy spectra by a factor of ∼3 up to 1000 TeV for five nuclei was consistent with extrapolation of directly measured spectra [8].

Hadronic Interaction Models
As explained in §3.2 the interpretation of EAS data to obtain cosmic ray composition is closely coupled with the choice of hadronic interaction models for low-and highenergy interactions.Ideally, the hadronic interaction models should be determined independent of composition to have greater confidence in results.The precision measurements of angular distribution of >1 GeV muons in GRAPES-3 offers one such possibility as outlined in §1.
A comparison of the observed muon angular distribution from GRAPES-3 experiment with the predictions from simulations using CORSIKA with specific combinations of low-and high-energy hadronic interaction models should in principle determine the most appropriate models of hadronic interactions.However, as a pilot study we compared the predictions of different combinations of low-and high-energy hadronic interaction models to estimate the sensitivity of this technique.We used three low energy models including GHEISHA [10], UrQMD [11] and FLUKA [12].QGSJet-II [13], SYBILL 2.1 [14] and EPOS 1.99 [15] were three high energy models used in simulations.
Since primary cosmic rays responsible for production of low energy muons detected by GRAPES-3 are influenced by geomagnetic field, the simulations also have to take account of it.Depending on location the geomagnetic field causes a sizable bending of particles.Since primary cosmic ray are mostly positively charged, the minimum A Lego plot of difference of normalized angular distributions using low energy models FLUKA and UrQMD with QGSJet-II as high energy model.An excess of +1.5% in E and a deficit of -1.1% in W directions, respectively.energy per unit charge, termed cutoff rigidity has a lower magnitude from west as compared to east.This effect is very significant for us due to proximity of geomagnetic equator to Ooty.Since charged particles are bent in presence of magnetic field, an accurate model of geomagnetic field "International Geomagnetic Reference Field (IGRF) 11th generation model" [19] was used.
The properties of detected muons including their energy and direction are governed by the interaction properties of parent hadrons.The angular distribution of detected muons is a result of superposition of several effects including, atmospheric attenuation, various cross sections, production heights, multiplicities and pseudorapidity distributions.Most of these properties are effectively parametrized in a given interaction model following certain prescriptions that may include extrapolation of properties from lower energy accelerator measurements.
Simulations have been carried out using COR-SIKA 6990 and all nine possible combinations of the three low-and three high-energy hadronic interaction models listed above were used for proton primaries with a differential spectral slope of -2.7 at the top of atmosphere.Primary protons were subjected to geomagnetic field (IGRF-2011) by employing the back-tracking method.Protons passing the geomagnetic cutoff and subsequent secondary charged particles were subjected to bending in geomagnetic field using multi-platform cosmic ray trajectory program [20].A total of 10 9 primary protons above 10 GeV were simulated.Muons above 1 GeV were tracked down to GRAPES-3 location.These muons then were subjected to detection threshold of sec(θ) GeV and resulting hits in PRCs were histogrammed into nine solid angle bins (V, N, E, W, S, NE, SE, NW, SW) exactly as is done for real data.Typically, 2×10 7 simulated muons were recorded for 10 9 primary protons.Actual number varied for each of the nine different combinations of models within ±10%.
To compare the predictions of different models, total number of muons in each of nine sets corresponding to the three low-and three high-energy models were normalized to 2×10 7 .Next, to compare two angular distributions, dif- ference in number of muons in a given bin, say 'N' was divided by mean number of muons in N and expressed in percent.A comparison of the predictions of low energy models FLUKA and UrQMD using SYBILL 2.1 as common high energy model is shown in Fig. 6.Largest excess of +1.9% in SE and biggest deficit of -1.1% in V bins, respectively are observed.A similar exercise where SYBILL 2.1 was replaced by QGSJet-II is shown in Fig. 7.
Here, largest excess of +1.5% in E and biggest deficit of -1.1% in W directions, respectively are seen.Clearly, the difference in the angular distributions for FLUKA and UrQMD with SYBILL 2.1 is unique and is significantly different when SYBILL 2.1 is replaced by QGSJet-II as shown in Figs. 6 and 7.
As an additional example, we examined the difference in angular distributions using GHEISHA and UrQMD with SYBILL 2.1 as shown in Fig. 8. Largest excess of +3.1% in NE and biggest deficit of -2.1% in V directions, respectively are seen.Interestingly, although the structure of the plots shown in Figs. 6 and 8 are distinctly different but they do display some similarity possibly indicating the common influence of two hadronic interaction models (SYBILL 2.1, UrQMD) used there.The unique shapes of the angular distributions shown in Figs.6,7,8 indicate that they represent the characteristics of underlying hadronic interaction models in a subtle but definite manner.Typical differences ∼1% in angular distributions should be possible to measure given that the efficiency of muon detector modules is known to a precision of 99.99% as discussed in §2 and that the exact configuration of the PRCs can be modeled with a precision of 1%.It is also important to note that despite a low muon threshold of 1 GeV, the angular distributions shown above are sensitive to both low-and high-energy hadronic interaction models.This is possibly because of proximity of GRAPES-3 location near geomagnetic equator resulting in a large cutoff rigidity of 17 GV in vertical direction.Consequently the median energy of primary protons is also much higher ∼80 GeV and falling at the boundary of the domains of the low-and high-energy models.

Summary
Measurements using high efficiency tracking muon detector in GRAPES-3 EAS array have been used to observe acceleration of muons during thunderstorms.The muon multiplicity distribution of EAS is used to probe the composition of primary cosmic rays below 1 PeV, with an overlap with direct measurements.Using a new approach we have shown that the angular distribution of >1 GeV muons can be used to determine most appropriate models from among several known low-and high-energy hadronic interaction models.

Figure 1 .
Figure 1.A schematic layout the GRAPES-3 air shower array shown here with 256 detectors (triangles).Each of the 16 squares represent a 35 m 2 area tracking muon detector with E µ ≥ 1 GeV.

Figure 2 .
Figure 2. Muon rates for modules, (a) 4, (b) 5 for 100 days from 1 Jan-10 Apr 2006.Ratio of muon rates, (c) before efficiency correction and its fourth order polynomial model, (d) after efficiency correction.

Figure 3 .
Figure 3.The muon rates of all 16 modules for 100 days from 1 Jan-10 Apr 2006 after efficiency correction using a fourth order polynomial model.

Figure 4 .
Figure 4. 12 panels from top left to right show muon rate in nine directions on 18 April 2011 from 15 to 17 h (IST) during a thunderstorm during successive 10 minutes.Significant deficit is seen in SE direction for ∼20 minutes.

Figure 5 .
Figure 5. Top panel shows muon rate during a thunderstorm on 18 April 2011 from 14 to 18 hours.Bottom panel shows the electric field at a location 6km east of GRAPES-3.The arrows indicate duration of activity to be about 20 minutes.

Figure 6 .
Figure 6.A Lego plot of difference of normalized angular distributions using low energy models FLUKA and UrQMD with SYBILL 2.1 as high energy model.An excess of +1.9% in SE and a deficit of -1.1% in V directions, respectively.

Figure 8 .
Figure 8.A Lego plot of difference of normalized angular distributions using low energy models FLUKA and UrQMD with QGSJet-II as high energy model.An excess of +3.1% in NE and a deficit of -2.1% in V directions, respectively.