proton discrimination

Astroparticle Physics 33 (2010) 206–215

Contents lists available at ScienceDirect

Astroparticle Physics journal homepage: www.elsevier.com/locate/astropart

Spectral statistics applied to TeV cosmic gamma/proton discrimination E. Faleiro a,*, L. Muñoz b, A. Relaño b,c, J. Retamosa b a

Departamento de Física Aplicada, E.U.I.T. Industrial, Universidad Politecnica de Madrid, 28012 Madrid, Spain Grupo de Física Nuclear, Departamento de Física Atómica, Molecular y Nuclear, Universidad Complutense de Madrid, 28040 Madrid, Spain c Instituto de Estructura de la Materia, C.S.I.C. Serrano, 123, 28006 Madrid, Spain b

a r t i c l e

i n f o

Article history: Received 11 December 2009 Received in revised form 2 February 2010 Accepted 3 February 2010 Available online 10 February 2010 Keywords: Spectral statistics Cosmic ray fluctuations Extensive air showers c/hadron discrimination

a b s t r a c t Spectral statistics techniques based on random matrix theory are applied to study the correlations at ground level among secondary particle density distributions of simulated extensive air showers. Highenergy interactions of TeV cosmic c-rays and protons with Earth atmosphere and the resulting extensive air showers have been simulated by means of the CORSIKA Monte Carlo code. The statistical description of shower fronts provided by this new type of analysis allows to distinguish events generated by primary protons from those initiated by primary c-rays, and leads to a powerful discrimination method whose efficiency is evaluated using the standard quality factors. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Extensive air showers (EAS) result from the interaction of high altitude atmospheric constituents with very high energy particles, mainly atomic nuclei and photons, arriving from any direction of space and collectively called primary cosmic rays [1,2]. New particles, mostly eþ ; e ; c and l, commonly called secondaries, are created as a consequence of these interactions, and interact again with atmospheric nuclei leading to a multiplicative cascade through the atmosphere. If the energy of the incident cosmic ray is high enough, secondary particles resulting after many interactions with the atmosphere may reach the Earth surface, as a wave front originated by the first interaction. Besides electrons, positrons and low energy photons, an important component of the EAS is the front of Cerenkov photons produced when the speed of secondary charged particles is higher than the speed of light in the propagation medium. The precise identification of primary gamma rays is very important for Cosmic Ray Astronomy, since only these particles retain information on the position of the emitting source. Unfortunately, it is known that for cosmic ray energies greater than 1 TeV, 99% of the primaries are hadrons, and hence identifying this small fraction of primary gamma rays inside a hadronic bulk is a very difficult task [3,4]. From the experimental measurement of the lateral

* Corresponding author. Present address: Ronda de Valencia 3, 28012 Madrid, Spain. Tel.: +34 913367686; fax: +34 913366850. E-mail address: [email protected] (E. Faleiro). 0927-6505/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.astropartphys.2010.02.001

distributions of secondary particles on large Earth array detectors, or the Cerenkov images collected by appropriated telescopes, it is possible to determine several properties of EAS, in particular, the nature of the primary particle, i.e. whether it is a high energy gamma ray or a hadronic particle. Although several separation methods based on the information contained on the front of secondary particles as well as on the distribution of Cerenkov light at the Earth level have been proposed, the introduction of new gamma–hadron separation methods is still an important task. Several present experiments can benefit from the method that we propose here. We can point out the Pierre Auger collaboration [5], designed to detect events in the energy range of EeV; the HEGRA–MAGIC collaboration [6,7], whose experimental setup is based on a Cerenkov telescope working in an energy range from some tens of GeV to a few TeV; the KASKADE–Grande experiment, with operative energy ranges from 10 PeV to some EeV [8], and the ARGO experiment, which is briefly described below [9]. Statistical properties of EAS have been studied from several points of view. Multifractal characterization of secondary particle distribution fluctuations [10], two-dimensional power spectrum properties of fluctuations [11] and Principal Component Analysis of particle distributions [12], among others, have given rise to very efficient gamma–hadron separation methods. All these methods highlight that the existence of strong correlations among secondary particles is an essential characteristic [13,14]. From an experimental point of view, a common technique is the use of large scintillator counter arrays. These are usually widely spaced to estimate the lateral distribution of an EAS by sampling the particle density over an area as large as possible. Short and intermediate

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

range correlations are lost, since only characteristic lengths larger than the counter separation are available. On the other hand, there are ground-based cosmic-ray detectors that locate the individual impacts of secondary particles on a quite large detection surface. An example is the previously mentioned ARGO detector, that consists of a single layer of resistive plate counters covering an area of 6700 m2, and which provides a detailed space–time picture of the shower front [9]. In this paper we propose a new approach to study the correlations among secondary particles. The eigenvalue spectrum of a matrix is usually characterized by strong correlations. As an example, the energy spectrum of the Hamiltonian matrix of a chaotic quantum system exhibit 1=f noise, while that of an integrable system exhibit 1=f 2 noise [15,16]. In this paper we consider a formal analogy between the impact positions at ground level of EAS secondary particles and the set of complex eigenvalues of a matrix without any particular symmetry. Therefore, this allows to study the correlation properties among secondary particles using the same methods developed to analyze the spectral fluctuations of matrix eigenvalues. They will reveal clear differences for the correlations of EAS generated by primary gammas and protons. These methods are commonly known as spectral statistics and are based to a large extent on Random Matrix Theory (RMT). RMT was introduced by Wigner in the mid 1950s to deal with the statistical properties of nuclear spectra and wave functions [17]. Nowadays RMT has developed into a powerful tool in statistical mechanics that can be successfully applied to describe generic properties of different systems, like atomic nuclei, complex atoms and molecules, disordered systems, one-dimensional interacting fermion systems, quantum chromodynamics and quantum gravity. A comprehensive review of the most important concepts and developments of RMT in quantum physics can be found in [18]. RMT has also been applied to a large number of fields like multivariate statistics, combinatorics, graph theory, number theory, biology, genomics and wireless communications [19]. The paper is organized as follows. Section 2 describes the different EAS databases used in this work, and explains how spectral statistics is adapted to study the statistical properties of EAS fronts. Section 3 reports the main results of the statistical analysis and the resulting separation method. The conclusions are summarized in Section 4.

207

parameters are needed to characterize the spectral points. Using polar coordinates, we denote the spectral density gðr; hÞ, and assume that it can be decomposed as

gðr; hÞ ¼ gðr; hÞ þ g~ðr; hÞ;

ð1Þ

where gðr; hÞ stands for the smooth part and g~ðr; hÞ for the fluctuating one. The cumulative spectral density

Gðr; hÞ ¼

Z

Z

r

0

h

gðu; v Þu du dv ;

ð2Þ

0

that counts the number of impacts up to radius r and angle h can be e hÞ. When the decomposed in a similar form Gðr; hÞ ¼ Gðr; hÞ þ Gðr; system exhibits axial symmetry, as is our case, an important simplification can be made. The smooth part of the spectral density depends only on the radial coordinate, that is, gðr; hÞ ¼ gðrÞ. Therefore, the smooth cumulative density reads

GðrÞ ¼ 2p

Z

r

gðuÞ du;

ð3Þ

0

and obviously the inverse relation,

gðrÞ ¼

1 dGðrÞ 2pr dr

ð4Þ

holds. Our aim is to find a change of variables so that the smooth part of the spectral density on the new variables is equal to one. Since the smooth part of the density only depends on the radial variable, we look for an appropriate change r ! q and keep the same angular variable h. Fig. 1 sketches the spectral points ðr; hÞ before the unfolding procedure is performed and the transformed points ðq; hÞ characterized by an uniform density. If we denote the new density as

cðq; hÞ ¼ cðqÞ þ c~ðq; hÞ ¼ 1 þ c~ðq; hÞ;

ð5Þ

the cumulative density is given by

e ðq; hÞ ¼ pq2 þ C e ðq; hÞ: Cðq; hÞ ¼ CðqÞ þ C

ð6Þ

Using (6) and the fact that the total number of spectral points up to a certain radius is the same regardless of the choice of the radial variable, i.e.,

GðrÞ ¼ CðqÞ;

ð7Þ

one easily gets that 2. Spectral statistics and extensive air showers

r!q¼ We introduce an analogy between single hits, characterized by their coordinates (x, y) on the detection plane, and the complex eigenvalues of a generic matrix without any symmetry property. We shall call those single hits spectral points. Spectral statistics has been successfully adapted to deal with the complex spectra of diagonalizable non Hermitian matrices like, for instance, the Ginibre ensembles [20,21]. In order to analyze the correlation among secondary particles, we will use spectral statistics methods adapted to study ensembles with complex spectra. Very often, the fluctuations of the spectral density, i.e., the number of spectral points per unit area, are modulated by its average behavior. Therefore, it would be convenient to remove the influence of the average behavior to compare the fluctuations for different events or even for different parts of a single shower front. This task is accomplished by means of an essential procedure called unfolding [21]. The unfolding procedure is well established for one-dimensional spectra, and consists in locally mapping the actual spectrum into another one with constant mean eigenvalue density. The starting point is the separation of the eigenvalue denðxÞ and a fluctuating part g ~ðxÞ, i.e., sity into a smooth part g ~ðxÞ. For complex or bidimensional spectra, two gðxÞ ¼ gðxÞ þ g

sffiffiffiffiffiffiffiffiffi GðrÞ

p

:

ð8Þ

Thus, the unfolding procedure requires the precise knowledge of GðrÞ. There is general agreement that, on average, the secondary particle density at level ground follows the NKG distribution [22]

GNKG ðrÞ ¼ a1

 ða2 2Þ  ða2 4:5Þ r r 1þ ; a3 a3

ð9Þ

where a1 ; a2 and a3 are real parameters to be estimated. For this reason we propose the following functional form for the smooth cumulative density,

GðrÞ ¼ GNKG ðrÞð1 þ Pn ðrÞÞ;

ð10Þ

where P n is a polynomial of degree n in r. In this work n  10, and the parameters of NKG(r) and the coefficients of Pn ðrÞ are obtained by means of a least-squares fit to the numerical cumulative density calculated with the spectral points. Fig. 2 compares the spectral densities of an EAS at ground level and a Poisson process in the plane. The left panel shows N ¼ 8290 impacts of secondary photons on the detection plane, placed at a depth of 800 g. Coordinates X and Y are given in meters. These

208

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

20

200

ρ

r 100

10

θ

0

-100

-10

Y

0

-200 -200

-100

0 X

100

200

θ

-20

-20

-10

0 X

10

20

Fig. 1. Secondary photons at the detection plane level generated by a primary photon. The impact positions are fixed by the coordinates ðr; hÞ (left panel) or by the coordinates ðq; hÞ once the unfolding procedure has been performed and the spectral density is constant.

Fig. 2. Left panel: one event with N ¼ 8290 secondary photons at the detection plane level, generated by a primary photon. Center panel: the same event after unfolding the  ¼ 1.  ¼ 1. Right panel: 10000 points distributed according to the bidimensional uniform distribution with c spectral density so that c

particles have been generated by the interaction of a 20 TeV photon with atmospheric nuclei at 18 km over the detection plane. The central panel shows the same set of hits after unfolding the spectral density. Since the smooth part of the ffi unfolded density is pffiffiffiffiffiffiffiffiffi c ¼ 1, the radius of the domain is R ’ N=p ’ 51 m. The right panel displays an ensemble of N ¼ 10000 points distributed accord ¼ 1. ing the two-dimensional uniform distribution with density c The similarity between the two sets of points plotted in the central and right panels suggests that the unfolding procedure is correct. Assuming that we have a large ensemble of experimental or simulated EAS the question is how to obtain the relevant statistical properties of its members. Spectral statistics introduces a set of convenient fluctuation measures or statistics of the observed (calculated) set of eigenvalues [23]. In this paper we shall consider two different statistics: the spacing distributions and the variance of the number of spectral points. (a) The nearest-neighbor spacing (NNS) distribution, denoted P(s), gives the probability that the spacing between a spectral point and its nearest neighbor lies between s and s + ds. Since only nearest impacts are considered, the NNS

distribution is a measure of the short-range correlations in the density of secondary particles. For the so-called Poissonian random process in the plane, i.e., N points thrown onto  ¼ 1 and no correlathe plane with uniform mean density c tions between throws, P(s) is given by the Wigner law [21]

PðsÞ ¼

p 2

p 2

se4s ;

ð11Þ

to be compared with the usual Poisson distribution, PðsÞ ¼ es , corresponding to the Poissonian process in the real axis. It is interesting to note that (11) describes the distribution of spacings between consecutive unfolded energy levels of chaotic quantum systems [24]. Similarly, one can define the k th nearest-neighbor spacing (k-NNS) distribution as the probability P(k, s) that the spacing between a spectral point and its k th nearest neighbor lies between s and s + ds. As far as we know, the P(k, s) distributions have no closed expressions for the Poissonian process when k > 1. For the sake of simplicity, we shall call Wigner laws to the spacing distributions of this process, regardless of the value of k. The cumulative distributions

209

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

Iðk; sÞ ¼

Z

s

Pðk; uÞdu

ð12Þ

0

will be of practical importance since they are not affected by binning effects. (b) Let N(A) be the number of impacts of secondary particles in a surface of area A of the detection plane. If the spectral density has been unfolded the average number of spectral points per unit area is one, and thus NðAÞ ¼ A. Here, the symbol  stands for an average over the whole ensemble of observed EAS. Thus, the variance of N(A) is given by

VðAÞ ¼ N 2 ðAÞ  A2 :

ð13Þ

Since A can take small as well as large values, V(A) explores the statistical properties of the spectral points at very different scales. For the Poissonian random process in the plane we have V(A) = A.

3. Analysis We have used the version 6.20 of the CORSIKA Monte Carlo code [25] to generate large databases of simulated EAS initiated from high energy cosmic c-rays and protons. The hadronic interactions at high energies are described by the QGSJET model, based on the parton-based Gribov–Regge theory [26], that simulates ultrarelativistic heavy ion collisions. The energy of the primary particle varies from 10 GeV to 60 TeV according to the power spectral law

Table 1 Database outline of primary cosmic rays and secondary particles considered in our analysis. We show the primary energy range, the number of events for each primary type, the average number and the standard deviation of different types of secondary particles. Primary type

Primary energy

Secondary particles N ev

Nc

N eþ

N e

rc

reþ

re

c

10–40 TeV 20–60 TeV

3000 3500

8484 10800

974 1430

1232 1769

7203 8695

959 1292

1187 1572

proton

dN / Eb ; dE

ð14Þ

where the value of the spectral index b depends on the type of primary particle. Most authors use b ¼ 2:5 for primary photons and b ¼ 2:7 for primary protons. However, the effect of this power law on our analysis is relatively unimportant because the events are classified according to the number of secondary particles that arrive to the detection plane, and the large dispersion in the number of secondary particles at fixed energy together with the wide energy range erase the differences with a flat starting spectrum. The height of the first interaction is chosen randomly according to the appropriate mean free path, and the observation level is set to 2200 meters above sea level, corresponding to the altitude of the Roque de los Muchachos, where the HEGRA experiment is located. As a simplification, all the showers are generated vertically. The effect of a non vertical shower axis may be incorporated to the analysis by geometrical reconstruction of the shower front from the arrival time structure at detection level. We have generated event databases for primary gammas and primary protons, whose main features are summarized in Table 1. 3.1. Spacing distributions and c/hadron separation After computing the spectral densities on the detection area with the coordinates of the secondary particles generated with CORSIKA, their fluctuations are isolated from the secular behavior by means of the unfolding procedure according to Eqs. (8) and (10). Once the separation of the fluctuating part has been achieved, we compute the values of P(k, s) for several values of k. The distributions P(k, s) are calculated for different types of secondary particles, namely low energy photons, electrons and positrons, separately. In order to obtain these distributions we calculate the distances between a certain spectral point and all the remaining points of the EAS, and select the smallest kx distances. By repeating this procedure for all the points pertaining to this event we obtain the distributions P(k, s) with k ¼ 1; 2; . . . ; kx corresponding to this event. Finally, we take into account the whole database to obtain the averaged distributions Pðk; sÞ. This is important, since our main goal is to find appropriate characteristics that allow the separation of the EAS induced by hadrons and photons. In what follows we

γ −> γ 0.8

p −> γ 0.8

0.6

0.6

P γ (1.s)

0.3

0

0.2 0 2

0

0.2

0.4

0.4

0 2

0.2

0

0.2

0.4

0.2

1.6

p

0.8

P γ (6.s)

0.1

1.2

0

0

0.2

0.4

0.1

1.2 0.8

0

0.6

0.4 0

0

0.2

1.6

γ

0.3

p

0.4

γ

P γ (1.s)

0.6

P γ (6.s)

0.6

0

0.2

0.4 0.6

0.4 0

1

2

s

3

4

0

0

1

2

3

4

s

Fig. 3. Average distributions P c ð1; sÞ and P c ð6; sÞ for secondary gammas pertaining to EAS generated by primary c (left panel) and primary protons (right panel).

210

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

label the k-NNS distributions by a script or a double script showing the primary particle that originated the event and the selected type of secondary particle. For instance, Pab ðk; sÞ is the k th NNS distribution of b secondary particles that have been generated by an a primary ray. Similarly, P b ðk; sÞ stands for any k th NNS distribution of b secondary particles, regardless of the primary particle. Fig. 3 shows the average distributions Pc ð1; sÞ and P c ð6; sÞ. The left panels correspond to events generated by primary photons while right panels display the result for events generated by primary protons. It is seen that the distributions Pcc ð1; sÞ and Pcc ð6; sÞ are very well described by the Wigner laws, except for a small difference near the maximum. Although we only display the distributions with k ¼ 1 and 6, this result is true for all k 2 ½1; 10, considered in this work. As previously stated, for bidimensional spectra the agreement with Wigner law means that the spectral points are completely uncorrelated. On the other hand, the averaged k-NNS distributions Ppc ð1; sÞ and P pc ð6; sÞ show clear discrepancies with Wigner laws. These are noticeable through the whole range of spacing values, and are larger for k ¼ 6 than for k ¼ 1. Even though the differences between the actual distribution and the Wigner laws are larger in the central region, the relative deviation is very important already for small spacings, as it can be seen in the panel insets. The situation is very similar for the spectral density of charged secondary particles. Since the statistical properties of these particles are almost identical, we shall display only the results for electrons. Fig. 4 shows the averaged spacing distributions P e ð1; sÞ and Pe ð6; sÞ, characteristic of secondary electrons generated by primary c-rays and protons. The agreement of Pce ð1; sÞ and P ce ð6; sÞ with the Wigner laws is quite good, showing that after the unfolding these spectral points are uncorrelated. On the contrary, the behavior of P pe ð1; sÞ and Ppe ð6; sÞ show marked deviations from Wigner curves. These deviations point out the existence of correlations that are absent in the events generated by c-rays. Nevertheless, they are smaller than those observed for secondary gammas. The different behavior of the distributions Pðk; sÞ for primary photons and protons, specially from k ¼ 6 onwards, unveils that the correlations among spectral points are more noticeable for EAS generated by protons than by photons. This fact has already been reported by the authors using different methods [10–14].

Thus, the spacing distributions may be used to separate EAS generated by protons from those generated by primary c-rays. Whenever we want to perform a precise comparison of these distributions, especially for individual events, it is helpful to use the integrated or cumulative distributions I(k, s), which are independent of binning effects. They are shown in Figs. 5 and 6 for k ¼ 1 and 6. Fig. 5 depicts the cumulative distributions Ic ð1; sÞ and Ic ð6; sÞ. Solid lines represent the averaged distributions: upper lines correspond to proton databases and the lower ones to primary-c databases. It also shows integrated distributions for a few individual events, using dashed lines for events generated by protons and dash dotted lines for events generated by photons. It is evident that the distributions Ic ð6; sÞ grow more slowly than the distributions Ic ð1; sÞ, reflecting that small 6 th nearest-neighbor spacings are rather scarce. One can also observe that, though the integrated distributions for individual events are different from one another, they are all near the global distribution of the associated database. What may be really important for our purposes is that the distance between the cumulative distributions Ipc ðk; sÞ and Icc ðk; sÞ is larger for k ¼ 6 than k ¼ 1. Our calculations show that these differences increase from k ¼ 1 to k ¼ 6 (shown in the figure) and remain essentially constant for higher values of k. This result is true for averaged as well as for the individual distributions. Fig. 6 displays the cumulative distributions of spacings between secondary eþ for k ¼ 1 and 6. Although the observed pattern is very similar, the gap between the curves corresponding to events generated by protons and photons is clearly smaller. These results can be used to separate events generated by protons from those generated by c-rays. For instance, we propose a simple method that uses cut-off parameters I0 ðkÞ corresponding to the values of Iðk; s0 ðkÞÞ at adequate small spacings s0 ðkÞ, and such that

 > I0 ðkÞ ! Primary ¼ hadron; Iðk; s0 ðkÞÞ < I0 ðkÞ ! Primary ¼ photon:

After performing an analysis of the distributions Iðk; s0 ðkÞÞ for different spacings s0 ðkÞ, we selected reasonable values for I0 ðkÞ and s0 ðkÞ. They are listed in Table 2 for k ¼ 1; 2; 4; 6; 8 and 10. As it can be seen the values of s0 lie in the initial part of the

_

0.8 0.6 (1.s)

0.4 0.2

P

0

0

0.2 0 2

0.2

0.4

0

0 2

0.2

0

0.2

0.4

0.2

1.6 0.1

1.2

0.8

0

p

e

_

(6.s)

0.1

1.2

0

0.2

0.4

0.6

0.4 0

P

(6.s) e_ γ

0.2

0.4 0.2

1.6

P

0.4

p

γ

e

_

0.4

0.6

P

(1.s)

p−> e

0.6

0.6 e_

_

γ −>e

0.8

ð15Þ

0.8

0

0

0.2

0.4

0.6

0.4 0

1

2

s

3

4

0

0

1

2

s

3

4

Fig. 4. Average distributions Pð1; sÞ and Pð6; sÞ for secondary e pertaining to events generated by primary c-rays (left panel) and primary protons (right panel).

211

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

0.08

0.08

I0

S0

0.06

I γ (6,s)

I γ (1,s)

0.06

0.04

0.02

0.04

0.02

I0

S0 0 0

0.2

0

0.6

0.4

0

0.2

0.4

s

0.6

s

Fig. 5. Cumulative spacing distributions Ic ð1; sÞ (left panels) and Ic ð6; sÞ (right panels) for secondary photons. Solid lines correspond to complete databases, the upper ones for primary protons and the lower ones for primary photons. Dashed lines represent various individual events for primary protons and dash dotted lines the same for primary photons.

0.08

0.08

I0

S0

0.06

I e + (6,s)

I e + (1,s)

0.06

0.04

0.04

0.02

0.02

I0

S0 0

0 0

0.2

0.4

0.6

0

0.2

0.4

0.6

s

s

Fig. 6. Cumulative spacing distributions Ieþ ð1; sÞ (left panel) and Ieþ ð6; sÞ (right panel) for secondary positrons. Solid lines correspond to averaged distributions, the upper ones for primary protons and the lower ones for primary photons. Dashed lines represent various individual events for primary protons and dash dotted lines the same for primary photons.

Table 2 Reference values of the separation spacings s0 ðkÞ as well as those of the cutting parameters I0 ðkÞ for secondary gammas and positrons. p,c ! eþ

k

p,c ! c s0

I0

s0

I0

1 2 4 6 8 10

0.25 0.35 0.40 0.45 0.55 0.55

0.064 0.037 0.015 0.011 0.018 0.013

0.25 0.35 0.45 0.50 0.55 0.60

0.060 0.030 0.015 0.010 0.010 0.013

distributions, but we would obtain similar results studying their tails. One could envisage improved separation methods, but more complicated, using the whole I(k,s) distributions. The distributions of Iðk; so ðkÞÞ values are shown in Fig. 7 for k ¼ 1 and 6. As it can be seen, the distributions for electromagnetic and hadronic events are pretty different. The former are displaced to the left and are clearly narrower than the latter, so that the overlap between the two kinds of distributions is quite small. The distributions for primary c-rays become sharper as k increases, leading to smaller overlaps with proton distributions. Moreover, this effect is more marked for secondary gammas than for secondary charged particles. Thus, we expect the efficiency of the separation method to increase with k and secondary photons. The vertical

212

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

+

p,γ −> γ

p,γ −>e

200

200 I0

N

150 100

100

50

50

0

0

0.05

I0

150

0.1

0.15

0

0

0.1 0.05 I(6,s 0)

0.15

200 I0

150 N

0.15

I(1,s0)

200

100

50

50 0

I0

150

100

0

0.1

0.05

I(1,s0)

0.05

0.1

0.15

0

0

I(6,s 0)

Fig. 7. Distributions of Iðk; s0 ðkÞÞ values for secondary photons (left panels) and secondary eþ (right panels). Solid lines correspond to events generated by primary photons and dashed lines to events induced by primary protons.

dashed line gives the selected value of the cut-off I0 ðkÞ, and thus events lying to the left or to the right of the separation line are considered of electromagnetic or hadronic origin, respectively. The gain obtained by a certain c/proton separation method is usually measured by the so-called quality factor Q, defined as [27,28]

Long-range correlations among spectral points are studied by means of the variance V(A). For A ¼ pq2 we obtain the values of this statistic as

VðqÞ ¼

jc jp

Q ¼ pffiffiffiffiffiffi ;

ð16Þ

where jc and jp are the fractions (efficiencies) of photons and procut tons kept by the algorithm, namely jc ¼ Ncut c =N c and jp ¼ N p =N p , where Ncut is the number of remaining particles after the cutting procedure. Table 3 shows the evolution of the quality factors Q I and their errors in terms of k. It can be observed a growing trend for k < 6 and then a saturation from k ¼ 6 onwards. This is due to the fact that the NNSD become essentially the same as k increases, once they are normalized and sðkÞ ¼ 1. Such quality factors make it possible to detect and to reject the majority of primary protons. For instance, the averaged quality factor for secondary gammas and electrons, at k ¼ 6, is Q Ik¼6  4:4, meaning that about 95.0% of background protons are rejected meanwhile the 98.4% of photons are kept. This gives rise to a significant enrichment in primary photon events.

Table 3 Quality factors Q I associated to the cumulative spacing distributions for secondary gammas and positrons in terms of the neighbor order k. QI

k

1 2 4 6 8 10

3.2. The variance of the number of spectral points

p,c ! c

p,c ! eþ

3.2 ± 0.2 3.7 ± 0.2 4.3 ± 0.3 4.9 ± 0.5 4.4 ± 0.4 4.4 ± 0.3

1.9 ± 0.1 2.2 ± 0.1 2.7 ± 0.2 2.6 ± 0.3 2.7 ± 0.2 2.6 ± 0.2

N 1X ðCi ðqÞ  pq2 Þ2 ; N i¼1

ð17Þ

where N is the total number of spectral points in a certain event, Ci ðqÞ is the number of points contained in a disk of radius q centered at the spectral point ðxi ; yi Þ, and pq2 is the expected number of spectral points in the disk. We recall, that after the unfolding pro ¼ 1, and thus the average number of cedure the average density is c spectral points in a disk of radius q is just pq2 . Then, repeating the calculation for the whole database we can obtain the averaged variance VðqÞ. Since this average as well as its fluctuations grow very quickly with q, we found more convenient to use the widths

rðqÞ ¼

pffiffiffiffiffiffiffiffiffiffiffi VðqÞ;

 ðqÞ ¼ and r

qffiffiffiffiffiffiffiffiffiffiffi VðqÞ:

ð18Þ

Fig. 8 shows the behavior of rðqÞ for two different types of secondaries; c particles in the left panel and e in the right one. Solid lines represent averaged widths: upper lines correspond to proton databases and lower ones to photon databases. The widths for a few individual events are also shown using different symbols: squares for events generated by protons and triangles for events generated by photons. Similarly to the behavior of the cumulative spacing distributions, the widths for individual events fluctuate, but they tend to be near the average curves. It should be noted that rp ðqÞ increases quicker than rc ðqÞ; moreover, for secondary gammas, the curves corresponding to proton events and primary-c events are usually apart and may be used to distinguish between both types. This behavior is less substantial for secondary electrons, and it is likely to find crossings between curves corresponding to the two types of events. This result confirms again that the correlations among spectral points are more marked for primary protons than for photons.

213

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

100

100 _

p,γ −> γ

p,γ −> e

80

80

ρ0

ρ0

60

σ

60

40

40

σ0

20

σ0

20

0

0 0

2

4

ρ

6

8

0

10

2

4

ρ

6

8

10

Fig. 8. Widths rp ðqÞ and rc ðqÞ for secondary photons (left panel) and electrons (right panel). Solid lines correspond to complete databases, the upper ones for primary protons and the lower ones for primary photons. Squares represent individual events generated by primary protons and triangles events initiated by primary photons. The dot dashed line correspond to widths of the Poissonian process in the plane.

0.5

0.5 _

p, γ −> γ

0.4

p, γ −> e

0.4 σ0

0.3

σ0

N

0.3

0.2

0.2

0.1

0.1

0

0

20

40 60 σ(ρ0 )

80 100

0 0

20

40 60 σ(ρ0 )

80 100

Fig. 9. Distributions of the rðq0 Þ values for secondary gammas (left panels) and secondary e (right panels). Solid lines correspond to events generated by primary photons and dashed lines to events induced by primary protons.

A convenient separation method based on this result is to use a cut-off parameter r0 corresponding to the value of the width rðqÞ at a certain radius q0 , and such that



rðq0 Þ

> r0 ! Primary ¼ hadron; < r0 ! Primary ¼ photon:

ð19Þ

The distributions of rðq0 Þ values corresponding to q0 ¼ 5 are shown in Fig. 9. The observed behavior is very close to that shown by cumulative spacing distributions. The curves associated to electromagnetic events are displaced towards small r values and are very narrow. On the contrary, those corresponding to protons are spread along a wider domain. The overlap between the two curves is small, specially for the case of secondary photons. The vertical

dashed line gives the selected value of the cut-off parameter, that in the present work is r0 ¼ 13:5. Events lying to the left or to the right of the separation line are considered of electromagnetic or hadronic origin, respectively. Using the aforementioned values of q0 and r0 we obtain Q Vc ¼ 4:9 and Q Ve ¼ 1:96, for secondary gammas and electrons, respectively. Fig. 10, compares different Q values, that have been obtained using the cumulative spacing distributions I(k, s) for different neighbor orders, and the variance of the number of spectral points VðAÞ. A first important feature is that Q Ic > Q Ie for all the k values used in this work. For secondary charged particles, it can be observed a growing trend for k < 4 and then a saturation from k ¼ 4 onwards. Something similar occurs for secondary photons, but the saturation seems to start at k ¼ 6. The horizontal lines in

214

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215 Table 4 Comparison of averaged quality factors evaluated by different methods. Note that methods that use the Cerenkov light have a range of energy different than those using EAS secondary particles as incoming data.

6 V γ

Q

5

4.9 4.4

4.3

4

4.4

3.7

Q

3.2

3

2.7

2.6

2.7

2.6

2.2 1.9

2

V _ e

Q

I Qγ I Q e_

1 0

0

2

4

6

8

10

k Fig. 10. Comparison of the quality factors Q V associated to the variance VðAÞ (dot dashed horizontal lines) to the quality factors Q I extracted from the cumulative spacing distributions (circles and diamonds).

both panels indicate the values of the quality factors Q V . It can be seen that Q Ve K Q Ie , for any value of k, and thus V(A) does not contribute significantly to the separation method in this case. On the contrary, Q Vc J Q Ic and it contributes a great deal to the discrimination of c-ray/proton events. This result suggest that there may exist a correlation among the values of the two statistics. To explore further this point Fig. 11 plots the re ðq0 Þ values versus those of the cumulative distribution Ie ð4; s0 Þ. Similar plots are obtained for other cumulative distributions proving that there exists a non simple correlation between these statistics. To end this section it may be interesting to compare these quality factors with those reported in the literature using different methods. Table 4 shows the values of several quality factors together with a brief description of the evaluation method, the energy rank involved in, and the corresponding bibliographic reference. Although there is a large dispersion of energy ranges, the simulations were performed under very similar conditions for both, the methods that use as input the Cerenkov light and also for those using the front of secondary particles. With the exception of the method proposed in[30], where the effect of the detector has been also partially simulated, the results are quite similar in all cases. Nevertheless, the values obtained in this work are slightly better probably because our procedure deals with individual impacts instead of particle densities, and therefore is able to put in evidence correlations that were missed by applying other procedures.

Fig. 11. Bidimensional plot representing the re ðq0 Þ values versus those of the cumulative distribution Ie ð4; s0 Þ for exploring correlations between these statistics.

Method

Reference

Energy rank

Q

Light-electron-slope (LES) method Fluctuations of Cerenkov densities Imaging Cerenkov Telescopes parameters Secondary densities distribution asymmetries Multifractal features of secondary densities 2D power spectrum of secondary densities PCA separation of secondary densities PCA separation of Cerenkov densities Cumulative spacing distributions ðk ¼ 6Þ Variance of the spectral point number

[29]

1 TeV–1 PeV

4.0

[30]

500 GeV–10 TeV

2.0

[31]

30 GeV–10 TeV

3.9

[10]

10 TeV–50 TeV

3.3

[11]

10 TeV–70 TeV

3.7

[12]

20 TeV–70 TeV

4.0

[14]

20 TeV–70 TeV

3.3

[13]

10 GeV–500 GeV

3.6

This work

10 Gev–60 TeV

4.4

This work

10 GeV–60 TeV

4.3

4. Conclusions The main purpose of this work is to point out the existence of measurable correlations in the EAS front that can be studied by considering individual impacts of secondary particles at Earth level. In order to achieve this goal we have applied RMT spectral statistics methods to simulated extensive air showers initiated by protons and c-rays with energies ranging from 10 to 60 TeV. As far as we know this method has not been applied so far to cosmic rays studies. For each EAS, the set of impacts is considered analogous to the complex eigenvalue spectrum of a generic matrix, without any type of symmetry, and we have called these hits spectral points. We have used spectral statistics measures like the cumulative spacing distributions I(k, s) or the variance of the number of spectral points V(A) to discriminate between events generated by the two types of primary rays. More precisely, we have defined discriminative parameters Iðk; s0 Þ and rðq0 Þ associated to these statistics, and whose distributions are clearly different for primary protons and c-rays. This fact allows to establish cutting values I0 and r0 (see Table 2) so that whenever Iðk; s0 Þ < I0 or rðq0 Þ < r0 we consider that the event has been generated by a primary c-ray. This method leads to an enrichment of primary photons in the initial data sample. The c/proton discrimination is important for high energy astronomy, since only electromagnetic EAS provide information on the direction of the source of the cosmic rays. The enrichment in primary cosmic c-rays is measured by quality factors Q associated to the cutting values of the separation parameters. Table 3 shows the Q values associated to the cumulative spacing distributions and Fig. 10 compares these values with the quality factors obtained using the variance V(A). The present values of the quality factors are quite similar or slightly better that those obtained using very different techniques, but it is worth mentioning that, with some exception, they represent upper bounds of the values that will be obtained in real experimental situations. In order to estimate the degradation of the separation power of a given method it would be necessary to know detailed characteristics of the experimental setup, like the type of detectors, their efficiency, the effective detection area, etc. Although further work is needed to study the combination of different methods, to explore the role played by other primary types, and the influence of factors such as not centered showers or experimental limitations, we conclude that the application of

E. Faleiro et al. / Astroparticle Physics 33 (2010) 206–215

spectral statistics to study the correlations among secondary particles of extensive air showers provides a new and promising c/proton separation method. Acknowledgments This work is supported in part by Spanish Government grants for the research projects FIS2006-12783-C03-02, CSPD-200700042-Ingenio 2010, and by the Universidad Complutense de Madrid grant UCM-910059. A. Relaño is supported by the Spanish program CPAN Consolider-imagenio 2010. References [1] P. Sokolsky, Introduction to Ultra-high Energy Cosmic Ray Physics, Addison Wesley, 1989. [2] T.K. Gaisser, Cosmic Rays and Particle Physics, Cambridge University Press, 1990. [3] D.J. Fegan, J. Phys. G 23 (1997) 1013. [4] R. Atkins et al., ApJ 595 (2003) 803. [5] D.V. Camin et al., The Auger Collaboration, Nucl. Instrum. Meth. A518 (2004) 172. [6] A. Karle et al., The HEGRA Collaboration, Astrop. Phys. 4 (1995) 1. [7] F.A. Aharonian et al., The HEGRA Collaboration, Astron. and Astrophys. 390 (2002) 39. [8] T. Antoni et al., Astrop. Phys. 14 (2001) 245. [9] G. Aielli et al., Nucl. Instrum. Meth. A 562 (2006) 511. [10] E. Faleiro, J.L. Contreras, J. Phys. G 24 (1998) 1795.

215

[11] E. Faleiro, J.M.G. Gomez, A. Relaño, Astrop. Phys. 19 (2003) 617. [12] E. Faleiro, J.M.G. Gomez, A. Relaño, J. Retamosa, Astrop. Phys. 21 (2004) 369. [13] J.L. Contreras, E. Faleiro, J.M.G. Gomez, R. Molina, L. Muñoz, A. Relaño, J. Retamosa, Astrop. Phys. 26 (2006) 50. [14] E. Faleiro, J.M.G. Gomez, L. Muñoz, A. Relaño, J. Retamosa, ApJS 155 (2004) 167. [15] A. Relaño, J.M.G. Gómez, R.A. Molina, J. Retamosa, E. Faleiro, Phys. Rev. Lett. 89 (2002) 244102. [16] E. Faleiro, J.M.G. Gómez, L. Muñoz, R.A. Molina, A. Relaño, J. Retamosa, Phys. Rev. Lett. 93 (2004) 244101. [17] E.P. Wigner, Interational Conference on the Neutron Interactions with the Nucleus, Columbia University, Columbia University Report CU-175, 1957. [18] T. Guhr, A. Müller-Groeling, H.A. Weidenmüller, Phys. Rep. 299 (1998) 189. [19] I. Dimitriu, Eigenvalue statistics for Beta-Ensembles, PhD thesis, Massachussetts Institute of Technology, 2003. [20] J. Ginibre, J. Math. Phys. 6 (1965) 440. [21] F. Haake, Quantum Signatures of Chaos, Springer, Berlin, 2001. [22] K. Greisen, Phys Prog. Cosmic Ray. 3 (3) (1965); Kamata, K. Nishimura, J. Theoretical Phys Prog. (6) (1958) 93. [23] M.L. Mehta, Random Matrices, Academic, New York, 1991. [24] O. Bohigas, M.J. Giannoni, C. Schmit, Phys. Rev. Lett. 52 (1984) 1. [25] D. Heck, J. Knapp, Extensive Air Shower Simulation with CORSIKA, Forschungszentrum Karlsruhe GmbH, Karlsruhe, 2002. [26] H.J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog, K. Werner, Physics Reports 350 (2) (2001) 93. [27] B.M. Schafer et al., Nucl. Instrum. Meth. A 465 (2001) 394. [28] R.K. Bock et al., Nucl. Instrum. Meth. A 516 (2001) 511. [29] C.E. Portocarrero, F. Arqueros, F.C. Andres, D.M. Borque, S. Martinez, J. Phys. G 25 (1999) 1249. [30] V.R. Chitnis, P.N. Bhat, Exper. Astron. 13 (2002) 77. [31] H. Krawczynski, D.A. Carter-Lewis, C. Duke, J. Holder, G. Maier, S. Le Bohec, G. Sembroski, V.R. Chitnis, Astrop. Phys. 25 (2006) 380.