Reanalysis of GU miniarray data using CORSIKA

Astroparticle Physics 22 (2005) 421–429 www.elsevier.com/locate/astropart

Reanalysis of GU miniarray data using CORSIKA U.D. Goswami

a,d,*

, K. Boruah a, P.K. Boruah b, T. Bezboruah

c

a

Department of Physics, Gauhati University, Guwahati 781 014, Assam, India Department of USIC, Gauhati University, Guwahati 781 014, Assam, India c Department of Electronic Science, Gauhati University, Guwahati 781 014, Assam, India Department of Physics, Central Institute of Himalayan Culture Studies, Dahung 790 116, West Kameng, Arunachal Pradesh, India b

d

Received 19 March 2004; received in revised form 25 May 2004; accepted 9 September 2004 Available online 8 October 2004

Abstract The GU miniarray is a ultra high energy cosmic ray (UHECR) detector consisting of eight plastic scintillators of carpet area 2 m2, each viewed by fast PMTs. It is used to detect Giant EAS by the method of time spread measurement of secondary particles produced in the atmosphere. The energies of the air showers have been reestimated using CORSIKA program. As in the original analysis the Cosmic Ray energy was determined via its relation to the ground level parameter Ns, the shower size. This relation was obtained previously through a best fit relation in agreement with QGS model and Yakutsk data. In this work we use CORSIKA code with QGSJET model of high energy hadronic interactions to simulate miniarray data leading to a modified relation between primary energy and shower size. A revised energy spectrum is reported for 1017–1019 eV primary energy.
2004 Elsevier B.V. All rights reserved.

1. Introduction In the ultra high energy range, the flux of Cosmic Ray particles is so low that direct detection of the primary particles is almost impossible. The traditional method of measurement in this region has been to observe extensive air showers (EAS) pro-

* Corresponding author. Address: Department of Physics, Gauhati University, Guwahati 781 014, Assam, India. Tel.: +91 036 12674537; fax: +91 036 12570133. E-mail address: [email protected] (U.D. Goswami).

duced in the earthÕs atmosphere, by primary Cosmic Ray particles, using ground based array of particle detectors spread a over wide area (several km2). Alternatively, this can be done by a low cost method, as suggested by Linsley [1,2], requiring a few closely packed detectors capable of measuring the arrival time spread of individual shower particles. The GU miniarray detector system is based on this idea and designed specially to measure both charged particle density and arrival time at the detector level [2]. Here we present a reanalysis of our old data taken during the period from October 1996 to April 1998. Previously the primary

0927-6505/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.astropartphys.2004.09.001

422

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

energy estimation was done by using the best fit relation in agreement with QGS model and Yakutsk data [3]. Based on Gribov–Regge theory and attempt to include QCD in a consistent way, a variety of more sophisticated models for hadron production have been developed in recent years. QGSJET [4] is such a model and high energy EAS simulated using the Monte-Carlo code CORSIKA [5] with QGSJET model as generator, has been proved to be very successful in explaining a variety of experimental results from 1012 to 1020 eV. In this paper we present a reanalysis of our miniarray data using CORSIKA6019, with QGSJET01 and a revised energy spectrum.

2. The miniarray experiment The miniarray experiment is based on the Linsley effect which is the increase in the spread of arrival times in a particle sample from a given shower with the increasing distance from the shower centre. Thus the measured time spread of particles striking a localized detector system gives an estimate of the distance (r) to the shower axis.

The number of particles give the measure of the local particle density (q). From these measurements, the shower size (Ns) is estimated from the assumed lateral distribution function q(Ns, r). Finally, the primary energy (E0) corresponding to an event with the estimated shower size Ns is derived using the best fit relation obtained from simulations using CORSIKA with QGSJET. The GU miniarray was located on the rooftop of the Physics building, Gauhati University (2610 0 N, 9145 0 E and altitude 51.8 m) [6]. The array consists of eight plastic scintillators viewed by fast PMTs. Fig. 1 shows the block diagram of the miniarray data acquisition system. The signals from the eight detectors are amplified and then carried to the control room via co-axial cables. In the control room all the eight signals are discriminated to provide corresponding logic signals. The discriminated output signals are then individually shaped into narrow pulses of 20 ns width and ORÕed together to give a serial pulse train. The serial pulse train is branched into two individual channels, one fed to the digital storaged oscilloscope and the other to the trigger unit. The trigger unit senses the incoming pulse train and generates

SCINTILLATION DETECTORS D1

D2

D4 D6

8 0 8 6 MICRO COMPUTER

D3

R

D5 D7

A M

8 CHANNEL

D8

C O U N T ER

MULTIPLEXER CONTROL

PRE AMPX8

I/O

PULSE WIDTH TRIGGER TRIG IN

500MHz D.S.O.

OR

8 DISCRIMINATORS

Y IN

TEKTRONIX TDS 520A GPIB

PC

COM PORT AT

COMPUTER

Fig. 1. Block diagram of the miniarray data acquisition system.

RS232C

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

423

the necessary trigger pulse. Once triggered, the oscilloscope trace is stopped and the recorded data is transferred to the PC via GPIB interface. The microprocessor monitors the status of the detectors by recording count rates at regular interval and also transfers count rate data to the PC via RS232C serial interface for storage [2]. 2.1. The detectors Each detector unit consists of a plastic scintillator block (size 50 · 50 · 5 cm3) having polyvinyltoluene base (with resolution 20% and decay time 4 ns), a fast photomultiplier tube (EMI 9807B) and a preamplifier unit enclosed into a light tight enclosure. The PMT is a round faced end window type with semitransparent bialkali photocathode having maximum sensitivity in the blue region of the spectrum. It operates at an anode potential of 1800 V and has a gain of 7.1 · 106 with rise time 2 ns. The anode pulse is coupled to the preamplifier unit which is a double stage differential amplifier with rise time 2 ns, gain 10 and overall band width 200 MHz [2]. All the eight detectors are housed in a hut constructed on the roof top of the Physics building, Gauhati University. The detector pulses after preamplification, are transmitted to the control room via co-axial cables (RG58U). In the control room pulses are processed by an eight channel discriminator with a trigger circuit, an ORing circuit and then pluses are displayed and stored in a digital storage oscilloscope with GPIB interface. A microcomputer based on 8086 microprocessor is used to monitor the operation of the detectors by recording count rate at predetermined intervals. Figs. 2 and 3 show a pulse train of a genuine shower selected by visual scan and a spurious event triggered by noise, respectively. Only genuine pulses are considered for analysis and noise pulses are effectively rejected by triggering criteria and visual inspection. The single particle rate for one detector channel is estimated to be 45 s1 from the integral Cosmic Ray flux for the area 0.25 m2. All the eight discriminator biases are adjusted at the individual single particle level. The minimum energy of the detected particle is about 100 MeV.

Fig. 2. A typical pulse train of a genuine shower recorded by the miniarray for an event with q = 1.5 m2.

Fig. 3. A spurious event triggered by noise.

For recording an event, a trigger pulse is generated under some prerequisite criteria viz, (i) a hardware trigger requiring particle in the range 2 or above within the time window, (ii) the minimum time spread between the particles must be 100 ns. Trigger rates were monitored daily over the length of the experiment and were found to be stable. The data for the period of 3.633 · 107 s (421 days) include for each event, the charged particle densities and arrival times of these

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

particles with reference to the leading particle (first arriving particle hitting the miniarray) [2]. 2.2. Calibration of detectors The integral Cosmic Ray flux of the secondary charged particles is 1.8 · 102 m2 s1, and thus the number of charged particles crossing the scintillator block of area 0.25 m2 is 45 s1. Therefore, the single particle rate for one channel of the detector array can be considered as 45 Hz. The calibration of the detector for single particle pulse height is performed using a standard single channel analyzer (ECIL, SC604B) and a counter. All the eight discriminator biases are adjusted at the individual single particle peak position. The single particle peak of a scintillation detector corresponds to a mean count rate of about 50 s1. This corresponds to a flux of 200 particles m2 s1 [2]. 2.3. Theoretical background 2.3.1. Simulation with CORSIKA A Monte-Carlo code to simulate the EAS development in the atmosphere, taking into account all knowledge of high energy hadronic and electromagnetic interactions involved is a good procedure to study the cascade development in the atmosphere. In the present work the extensive air showers are generated by CORSIKA [5] version 6019 with QGSJET01 [4] for high energy hadronic interactions. CORSIKA is an open program package for performing a complete 4-dimensional simulation of air shower with primary energies from 1012 to 1020 eV. It was originally developed to perform simulations for the KASCADE experiment at Karlsruhe [7] and it has been refined over the past years. An important task in EAS Monte-Carlo simulation is to take into account all knowledge of high energy hadronic and electromagnetic interactions, although the most serious problem is the extrapolation of hadronic interaction to higher energies, which is not covered by experimental data in accelerators. For these reasons CORSIKA provides different hadronic interaction models at high energies. In this work we invoked the QGSJET01 [4], to simulate high energy hadronic interactions.

This model is very successful in explaining experimental results in the high energy range. For our Monte-Carlo data library we have simulated 500 vertical showers in the energy range from 1017 to 1020 eV with an energy slope of 2.65, proton and iron as primary particles. Statistical thinning [8–10] of shower particles was applied at the level of 104E0 with a maximum particle weight limit of 1030. The threshold energies are 0.1 GeV for hadrons, muons and electromagnetic component. 2.3.2. Numerical equations The relation between shower disk thickness r (ns) and the core distance r (m) has been derived by Capdevielle et al. [11] from their simulation with CORSIKA for near vertical shower as given by
r b rðrÞ ¼ B 1 þ ; ð1Þ c where B = 2.6, c = 25 and b = 1.4. Fig. 4 shows the comparison of shower disk thickness at different distances from shower axis as calculated from Eq. (1) with results from LinsleyÕs original equation [2,12] for the same. This figure shows a reasonable agreement between these two equations with slight discrepancy in the higher core distances. The particle density distribution for large shower and medium core distances (100 m < r < 1000 m) has been parameterized by using our

900 Linsley Corsika

800 700 600 σ (ns)

424

500 400 300 200 100 200

400

600

800 1000 1200 1400 R (m)

Fig. 4. Shower disk thickness versus distance from the shower axis. LinsleyÕs relation is shown in comparison with CORSIKA simulation.

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

Density (m-2)

old q1 = 1.5 m2 and its inner radius determined by the minimum time spread r1 as selected. From Eqs. (1) and (2),

Proton Iron Old Relation

1000

100

10

rmin ¼ c½ðr1 =BÞ1=b  1;

ð7Þ

 1=n N s ¼ : q1

ð8Þ

rmax

1 100

425

The minimum detectable shower size is given by the condition

1000 Distance (m)

AðN min Þ ¼ 0;

Fig. 5. Lateral particle density distribution obtained from simulation with CORSIKA and from old relation.

which gives N min ¼ ðq1 =Þ½cfðr1 =BÞ1=b  1gn :

n

qðN s ; rÞ ¼ N s r ;

ð2Þ

where Ns is the size of the shower. For proton primary  = 0.00053, n = 1.5 and for iron primary  = 0.0086, n = 1.9. Fig. 5 shows the lateral density distribution for proton and iron initiated showers simulated for the miniarray using CORSIKA as compared with earlier relation [2,13] used for miniarray analysis. There is a noticeable disagreement between CORSIKA and the earlier relation. The integral and differential shower size spectra [2,3] are J ðN s Þ ¼ DN c s ;

jðN s Þ ¼ cDN c1 ; s

ð3Þ

where the constants have values D = 318 and c = 1.65. With the help of the above relations we derived the frequencies of LinsleyÕs events as a function of the minimum time spread (r1), the threshold density (q1) and the shower size (Ns) as [2] Z 1 F L ð> r1 ; q1 Þ ¼ AðN s ÞjðN s Þ dN s ; ð4Þ N min

F L ð> N s > r1 ; q1 Þ ¼

Z

1

AðN s ÞjðN s Þ dN s ;

ð5Þ

2.3.3. Detector response Simulation of detector response for all the charged particles in the vertical direction of miniarray detector has been performed. Particles are collected for 500 simulated showers within the annular area ranging from Rmin to Rmax and for primary energies from 1017 to 1020 eV with minimum detectable energy of the particles as 100 MeV. From this, mean particle density at different core distances have been calculated. For the minimum detectable shower size and the threshold density, events are counted for all those simulated showers for different primary energy bins. Fig. 6 shows the variation of effective area of miniarray with primary energy. The effective

8

10 Effective annular area (m2)

simulated miniarray data for proton and iron initiated showers separately is as follows:

ð9Þ

Proton Iron Old Relation

7

10

6

10

5

10

Ns

1017

where AðN s Þ ¼ pðr2max  r2min Þ:

ð6Þ

Here, the effective area A(Ns) is an annular ring with outer radius determined by density thresh-

1018 Energy (eV)

Fig. 6. Effective area of miniarray versus primary energy for proton and iron shower simulation using CORSIKA compared with the results from old relation.

Number of events

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

Shower size

426

108

7

10

1017

Proton Iron Old Relation 1018 Energy (eV)

100

10

100

1000 σ (ns)

Fig. 7. Shower size (Ns) versus primary energy for proton and iron shower simulation using CORSIKA. Solid line represents energy calibration curve for proton and the dotted line for iron.

Fig. 8. Event number distribution as a function of shower disk thickness (r).

area increases with the increasing energy. The simulated data is not in close agreement with the data from old relation [2].

ground soft radiations and small air showers. In order to eliminate the large number of small air showers a minimum time spread has to be assigned. For our miniarray of 2 m2 area, minimum acceptable shower sizes are 1.57 · 107 for proton and 9.67 · 106 for iron with a minimum time spread r = 100 ns. In view of the small particle density encountered, each scintillator is not expected to receive more than one particle at a time from a shower. Fig. 8 shows the event distribution as a function of shower disk thickness (r).

2.3.4. Energy calibration The primary energy for each shower is reconstructed using the estimated value of Ns. The calibration curve for proton and iron primaries are plotted in Fig. 7 and are compared with the simulation and the best fitted relation with Yakutsk data [2,3], that was used for previous miniarray data analysis. The energy calibration is parameterized as E0 ðeVÞ ¼ aN bs ;

ð10Þ

where a = 2.217 · 1011, b = 0.798 with an approximate percentage error of 12% for proton and a = 6.194 · 1011, b = 0.898 with an approximate percentage error of 8% for iron. There is a noticeable difference between the old and the new calibration, which leads to the reconstruction of energies. 2.4. Data selection criteria Numerical calculations show that for a given threshold density (q1 = 1.5 m2), the minimum detectable shower size increases and the shower rate decreases with increasing time spread. The miniarray should be able to pick out at least few large air shower events from a swarm of irrelevant events including the counter noises, the back-

3. Results and discussion We present here the reanalysis of data that were collected from October 1996 to April 1998. Most of the data collected do not belong to true large air shower events. The data are reduced by the selection process and by visual inspection. True large shower events with a time spread of shower front r P 100 ns and with the local particle densities q P 1.5 m2 are collected and analyzed. In the present reanalysis using CORSIKA, the following new results are derived. 3.1. Shower rate and shower size spectrum The integral shower rate spectrum of the selected events with q1 = 1.5 m2 as function of time spread is shown in Fig. 9. The errors are estimated by considering the ±10 ns instrumental error. Here

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

427

4 3

FL(>N) (m sr eV day )

-1

2

10

101 100 10-1 10-2 10-3

10

2

1

10

100

10-4 100 200 300 400 500 600 700 σ (ns)

Fig. 9. Integral shower rate spectrum with q1 = 1.5 m2 and r1 = 100 ns. Data for proton and iron are obtained from reanalysis of experimental data using CORSIKA assuming proton and iron as primary Cosmic Ray particles.

unfolding of data is done by first parameterizing the integral flux as a function of time spread (r) from the simulation with the CORSIKA. Then the parameterized equation is used to unfold the values of the integral flux for different experimental values the r. The corresponding parameterized equation is

107

108 Shower size

Fig. 10. Integral shower size spectrum with q1 = 1.5 m2. Data for proton and iron are obtained from reanalysis of experimental data using CORSIKA assuming proton and iron as primary Cosmic Ray particles.

data obtained using CORSIKA is compared with the same from earlier analysis. It is observed that disagreement is more for small shower size. From the new analysis it is observed that shower flux decreases almost linearly with increasing shower size for both proton and iron primaries. 3.2. Energy spectrum

ð11Þ

where F0 = 0.00015, b = 974.284, c = 0.91 with an approximate percentage error of 19% for proton flux and F0 = 0.0029, b = 72 866.7, c = 1.78 with an approximate percentage error of 16% for iron flux. This figure also shows the comparison of the result of the old analysis with the present reanalysis under proton and iron assumptions. A considerable disagreement is observed for higher values of time spreads. This indicates a more rapid decrease of shower rate with increasing time spread than as observed in the earlier analysis [2]. Integral shower size spectrum considering events with threshold density q1 = 1.5 m2 is shown in Fig. 10. The same procedure as explained above is also adopted to plot Fig. 10. Here the corresponding parameterized equation is ð12Þ

where F0 = 2.557 · 108, b = 0.94 with an approximate percentage error of 22% for proton flux and F0 = 3.773 · 107, b = 0.83 with an approximate percentage error of 24% for iron flux. Reanalyzed

Fig. 11 shows differential energy spectrum derived from the reanalysis of the miniarray data assuming proton and iron primaries as compared

1018 j(E) E 2.5 (m-2sec-1sr-1eV1.5)

F L ¼ F 0 expðb=rc Þ;

F L ¼ F 0 N bs ;

Proton Iron Old analysis

-1

Proton Iron Old analysis

10

-2 -1

-2 -1

-1

-1

FL(>σ1) (m sr eV day )

10

17

10

1016

Proton Iron Old analysis Yakutsk AGASA HiRes

15

10

1014 1017

1018 1019 Energy (eV)

Fig. 11. Miniarray differential energy spectrum. The differential flux is multiplied by E2:5 0 . Data for proton and iron are obtained from reanalysis of experimental data using CORSIKA assuming proton and iron as primary Cosmic Ray particles. Region between solid lines give the flux as compiled from Akeno and Haverah Park data by Nagano and Watson (2000) [15].

428

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

with old analysis and results of Yakutsk, AGASA and HiRes [14]. A compilation of spectrum derived from Akeno and Haverah Park data by Ave et al. [15] is also shown in Fig. 11 for comparison. It is observed that (1) new analysis gives estimates of primary energy significantly higher than the previous analysis. Energy spectrum after reanalysis is found to span from 1017 to 1019 eV for proton primary and from 1017 to 1019.4 eV for iron primary and further confirms the irregular behavior of energy spectrum at ultra high energy region with a prominent dip; (2) the differential energy spectrum shows structure similar to that observed by other world groups. Spectral breaks are found to occur at higher energies compared with old method of analysis. The spectrum becomes steeper around 1017.5 and 1017.7 eV and flattens around 1018.7 and 1019.1 eV for proton and iron primaries, respectively forming a dip. Earlier analysis showed a dip around 1018.2 eV. A comparison of different features with other world data are shown in Tables 1 and 2; (3) there is a significant difference between spectra predicted by pure proton and pure iron assumptions. However, proton assumption results agree more with Nagano and Watson compilation [15] within spectral range from 1017.5 to 1018.5 eV. Beyond 1019 eV results of different giant arrays are contradictory. AGASA provides strong evidence for the existence of Cosmic Rays with energies beyond Greisen–Zatzepin–Kuzmin (GZK) Table 1 A comparison of slope before the dip of the differential energy spectrum Experiment

Slope before the dip

Energy range (eV)

Yakutsk AGASA HiRes-I Miniarray (old analysis) Miniarray (new analysis, p) Miniarray (new analysis, Fe)

3.195 ± 0.009 2.981 ± 0.058 3.185 ± 0.059 3.468 ± 0.131 2.733 ± 0.094 2.538 ± 0.081

1017.5–1019.2 1017.6–1018.9 1017.2–1018.5 1017.4–1018.2 1017.5–1018.7 1017.7–1019.1

Table 2 A comparison of overall slope of the differential energy spectrum Experiment

Slope

Energy range (eV)

Yakutsk AGASA HiRes-I Miniarray (old analysis) Miniarray (new analysis, p) Miniarray (new analysis, Fe)

3.031 ± 0.047 2.884 ± 0.059 3.151 ± 0.036 2.938 ± 0.108 2.360 ± 0.075 2.207 ± 0.067

1017.5–1019.8 1017.6–1019.2 1017.2–1019.2 1017.0–1018.8 1017.0–1019.0 1017.0–1019.4

cut-off. By contrast, data from the HiRes detector in the US are compatible with the existence of the GZK cut-off. More statistics is therefore necessary to resolve this important feature related to Astrophysics; (4) the differential energy spectrum corresponding to best least square fit in the energy region 1017.0–1019.0 eV is derived as jðE0 Þ ¼ j0 Ep 0 :

ð13Þ 12

For proton primary j0 = 7.107 · 10 , p = 2.360 ± 0.075 and for iron primary j0 = 2.772 · 1010, p = 2.207 ± 0.067. (5) A relatively higher flux is seen in the reanalysis of both for proton and iron, beyond the dip, as compared with other world data. This overestimation of event rate at higher energy range may be due to inclusion of some delayed particles, small detector area and associated triggering criterion. Our efforts are on to minimize the errors of miniarray detectors by increasing the area of the miniarray such that more particles can be detected per event and unwanted events may be distinguished by careful investigations.

Acknowledgements One of the authors (UDG) is immensely thankful to Dr. Dieter Heck of Institut fu¨r Kernphysik, Forschungszentrum Karlsruhe, Germany, Dr. Varsha R. Chitnis, Prof. P.N. Bhat and Prof. B.S. Acharya of Tata Institute of Fundamental Research, Mumbai, India for their painstaking

U.D. Goswami et al. / Astroparticle Physics 22 (2005) 421–429

help in various aspect of using CORSIKA code and analyzing data with the same.

References [1] J. Linsley, G. Phys. G 12 (1986) 51. [2] T. Bezboruah, K. Boruah, P.K. Boruah, Astropart. Phys. 11 (1999) 395. [3] A.M. Hillas, Phys. Report C 20 (1975) 79. [4] N.N. Kalmykov, S.S. Ostapchenko, A.I. Pavlov, Nucl. Phys. B 52 (Proc. Suppl.) (1997) 17. [5] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, T. Thouw, Report FZKA 6019, Forschungszentrum Karlsruhe, 1998. [6] T. Bezboruah, K. Boruah, P.K. Boruah, Nucl. Inst. Meth. Phys. Res. A 410 (1998) 206.

429

[7] K.H. Kampert et al., KASCADE Collaboration, in: Proc. 26th Int. Cosmic Ray Conf., Salt Lake City, USA, vol. 3, 1999, p. 159. [8] W.R. Nelson, H. Hirayama, D.W.O. Rogers, Report SLAC 265, Stanford Linear Accelerator Centre, 1986. [9] D. Heck, J. Knapp, Report FZKA 6019, Forschungszentrum Karlsruhe, 1998. [10] M. Hillas, Nucl. Phys. B 52 (Proc. Suppl.) (1997) 29. [11] J.N. Capdevielle et al., in: Proc. 28th Int. Cosmic Ray Conf. Tsukuba, Japan, vol. 2, 2003, p. 217. [12] J. Linsley, L. Scarsi, Phys. Rev. Lett. 9 (1962) 123. [13] T. Hara et al., in: Proc. 25th Int. Cosmic Ray Conf., Durban, vol. 6, 1997, p. 229. [14] A.V. Glushkov et al., in: Proc. 28th Int. Cosmic Ray Conf., Tsukuba, Japan, vol. 1, 2003, p. 389. [15] M. Ave et al., in: Proc. 27th Int. Cosmic Ray Conf., Hamberg, Germany, vol. 1, 2001, p. 381.