Revision of the energy calibration of the Yakutsk EAS array
aa r X i v : . [ a s t r o - ph . H E ] A ug Revision of the energy calibration of the Yakutsk EAS array
A. V. Glushkov, M. I. Pravdin, and A. Sabourov ∗ Yu. G. Shafer Institute of cosmophysical research and aeronomy and677980, Lenin Ave. 31, Yakutsk, Russia
Responses of surface and underground scintillation detectors of Yakutsk arrayare calculated for showers initiated by primary particles with energy E ≥ eVwithin the frameworks of QGSJet01d, QGSJet-II-04, SIBYLL-2.1 and EPOS-LHChadron interaction models. A new estimation of E is obtained with the use ofvarious methods. The resulting energy is lower compared to the obtained withearlier method by factor ∼ . PACS numbers: 96.50.sb, 96.50.sdKeywords: cosmic rays, extensive air showers, energy spectrum
I. INTRODUCTION
Energy spectrum of ultra-high energy cosmic rays (UHECR), cosmic rays (CR) with E ≥ eV energy, is one of the main links in the chain of complex problems associatedwith understanding the nature of primary particles with such energy. Mechanisms of theirproduction and acceleration in astrophysical sources and various effects they experience dur-ing the propagation in the Universe have direct effect on the observed primary CR spectrum.Recently, a significant progress has been achieved in interpretation of its structural featuresin the ultra-high energy domain. The black-body cutoff at E ∼ × eV predictedby Greisen [1] and Zatsepin and Kuz’min [2] (the GZK cutoff) has been confirmed [3, 4],thus pointing at extragalactic origin of the most energetic CR particles. The second knee(at ∼ eV) and the ankle (at ∼ eV) are commonly associated with a transitionbetween galactic and extragalactic CR components; and though its exact location on theenergy scale is not known precisely, there are plenty of theoretical scenarios compatible withexisting experimental observations (e.g. [5–7]). ∗ tema@ikfia.sbras.ru However, there is certain discrepancy in the world array of experimental results. CRspectra measured by various UHECR experiments [8–12] confirm such spectral features asthe ankle and the second knee, but they differ from each other in absolute intensity by factorof almost 2 [13, 14]. In particular, the spectrum measured by Yakutsk experiment lies aboveall the world data. In this context, the data published by the Yakutsk group signify theupper limit of the spectrum intensity, and data from the Pierre Auger Observatory (PAO) —its lower limit.Such situation to a large extent stems from the fact that the only available method ofUHECR observation is indirect, conducted by registering cascades of secondary particlesproduced by primary UHECRs in the Earth’s atmosphere: extensive air showers (EAS).Most of largest experiments employ differing observational techniques and, consequently,rely on different methods to reconstruct the energy of primary particles. Hence, one cannotdo without theoretical notion of EAS development.The Yakutsk EAS array stands out from other large arrays for its complexity: since itis equipped with detectors of three types, it simultaneously registers several shower compo-nents. Charged particles (electrons, positrons and muons) are recorded with 2 m surfacescintillation detectors (SSD). Muon component arising from nuclear interactions is regis-tered with detectors of the same type placed below the ground level, in order to preventelectromagnetic contamination by creating a shield with 1 × sec θ GeV threshold. Cherenkovlight emitted by EAS charged particles is recorded with integral Cherenkov detectors basedon the FEU-49 photomultiplier tube.Cherenkov component carries information about ∼
80 % of primary energy dissipated inthe atmosphere and, thus, enables to determine E with calorimetric method [15–19]. Thismethod defines the E as a sum of energies of all EAS components and connects it withexperimentally measured value ρ s, (it will be discussed in greater detail in Section III).Originally, it was introduced in [20] for energies ∼ eV. In Yakutsk experiment it wasapplied to showers with E ≃ (1 . − × eV at zenith angles θ ≤ ◦ [15, 16] andresulted in the following approximation for primary energy reconstruction: E =(4 . ± . × · ( ρ s, (0 ◦ )) . ± . (eV), (1) ρ s, (0 ◦ ) = ρ s, ( θ ) × exp (cid:16) (sec θ − · x λ ρ (cid:17) , (2) λ ρ = 400 ±
45 (g/cm ), (3)where x = 1020 g/cm , ρ s, ( θ ) is the density of charged particles (m − ) measured bySSDs at the distance R = 600 m from shower axis and λ ρ is attenuation length. Later, therelations (1) and (3) were changed slightly (see [17–19]): E = (4 . ± . × · ( ρ s, (0 ◦ )) . ± . , (4) λ ρ = (450 ±
44) + (32 ± · log ρ s, (0 ◦ ). (5)The intensity of the CR energy spectrum estimated with the use of (4) turned out to besignificantly higher than the world data (see e.g. [21]). In [22] estimation of E for Yakutskdata was presented obtained for primary protons within the framework of QGSJet01 model,which was 1.6 times lower than (4). Here we consider energy calibration of registered showersbased on modern CORSIKA code (version 6.7370) [23]. II. LATERAL DISTRIBUTION OF THE DETECTOR SIGNAL
Basic EAS parameters measured in Yakutsk experiment (arrival direction, coordinatesof the shower axis, primary energy) are reconstructed with lateral distribution function(LDF) of all particles (electrons, muons and high-energy photons) which are registered bySSDs. These particles pass through a multilayer shield consisting of snow, iron, wood andaluminium (total thickness amounts to ∼ . ) and then — through a 5 cm thickscintillator (with the density ∼ .
06 g/cm ). The energy deposit in a scintillator ∆ E s ( R ) isproportional to the number of particles passed though a detector and is measured in relativeunits: ρ s ( R ) = ∆ E s ( R ) E (m − ), (6)where E = 11 .
75 MeV, which is the energy released in a scintillator during the passage ofa vertical relativistic muon (the response unit).Scintillation detectors are calibrated and controlled with the use of amplitude densityspectra from background cosmic particles [24]. Herewith, the integral spectra of two typesare used. First one — a spectrum from a single detector, which is controlled by a nearbydetector mounted in the same station (the so-called “spectrum of double coincidence” withthe frequency ≃ (2 −
3) s − ). Second one — uncontrolled spectrum with the frequency ∼
200 s − , which is used to calibrate muon detectors. Both spectra are described by apower law: F ( > ρ ) ∼ ρ − η ∼ U − η , (7)where values of η for both spectra were obtained experimentally. For spectrum of the firsttype η = 1 . η = 3 . ρ = U/U — particle densitymeasured in units of signal amplitude U from a reference detector during the passage ofa vertical relativistic cosmic muon. The procedure of calibration and control consists ofcontinual monitoring of the U value in all detectors by periodical measurements of theirdensity spectra. The procedure is performed once a two days. Spectra of double coincidenceare collected for 2 hours, uncontrolled spectra — for 30 minutes.Within the framework of models QGSJet01d [25], QGSJet-II-04 [26], SIBYLL-2.1 [27]and EPOS-LHC [28] we calculated LDF of the SDD response in showers, initiated by pri-mary protons and iron nuclei with energies 10 . − . eV arriving at different zenithangles. FLUKA package [29, 30] was chosen for treatment of lower energy interactions.At first, the response u m ( ǫ, θ ) from a single particle of a type m (where m is electron,muon or gamma-photon) with energy ǫ was calculated. During the calculation, all the mainprecesses occurring in the detector during energy release/consumption with correspondingcross-sections were put into consideration: ionization and bremsstrahlung – for charged par-ticles; pair production and delta-electrons from Compton effect – for gamma-photons. Thenthe development of air shower was simulated with CORSIKA code. For each set of primaryparameters (mass of primary particle, its energy and incident zenith angle) 500 showers weresimulated. In order to speed-up the simulation, the “thin-sampling” mechanism, introducedin [31], was activated in the CORSIKA code [23, 32]. The thinning level ǫ th. = ǫ min /E , con-trolling the minimal energy of secondary particles ǫ min treated by CORSIKA, was defined inthe interval 3 . · − , − and weight limit of secondary particles w max — in the interval10 , . · , depending on the primary energy. This was done in order to limit the growthof artificial fluctuations induced by thin-sampling in showers with lower energies.During conversion to density, the number of particles was calculated in the detector of agiven area. Resulting showers were averaged together and mean energy spectra d m ( ǫ, R, θ )were calculated for all particle types in intervals (log R j , log R j + 0 . R was defined as a sum of responses: ρ s ( R ) = X m I m X i =1 u m ( ǫ i , θ i ) d m ( ǫ i , R, θ i ), (8) FIG. 1. Energy dependence of log ( ρ s, ( θ ) /E ) for primary protons and iron nuclei accordingto predictions of QGSJet01d for vertical showers. where I m — the number of particles of a type m hitting a detector at a distance R .On Fig.1 the dependence of the value log ( ρ s, (0 ◦ ) /E ) from E is shown for primaryprotons (open circles) and iron nuclei (black circles) as predicted by QGSJet01d model.They satisfy the relation: E = (3 . ± . × · ( ρ s, (0 ◦ )) . , (9)Other models — QGSJet-II-04, SIBYLL-2.1 and EPOS-LHC — give the following esti-mations correspondingly: E = (3 . ± . × · ( ρ s, (0 ◦ )) . , (10) E = (3 . ± . × · ( ρ s, (0 ◦ )) . , (11) E = (3 . ± . × · ( ρ s, (0 ◦ )) . , (12)Averaging over all models gives the dependence: E = (3 . ± . · ( ρ s, (0 ◦ )) . , (13) FIG. 2. Zenith-angular dependence of log ( ρ s, ( θ ) /E ) according to QGSJet01d predictions forprotons and iron nuclei with energy E = 10 , and 10 eV. which resulted in a lower estimated value of E by factor 1 .
20 when compared to (1) and by1 .
41 when compared to (4).Zenith-angular dependences of log ( ρ s, ( θ ) /E ) according to QGSJet01d model areshown on Fig.2. They satisfy a linear dependency with λ ρ = 415 ±
15 g/cm at any compo-sition of primary CR when sec θ is lesser than:sec θ lim. = a + b log ρ s, ( θ ), (14)where a = 1 .
26 and b = 0 . a = 1 . E = 10 eV, with attenuation length λ ρ = 415 g/cm and for θ ≤ ◦ . In other cases the dependency is more complex. III. CALORIMETRIC METHOD
We considered the energy balance starting from the example of experimental datafrom [15, 16]. Earlier, these data had provided a basis for the calorimetric method of
TABLE I. Observables of EAS with E = 10 eV and cos θ = 0 .
95 from primary nuclei ( A )according to CORSIKA [23] simulation and experiment [16]. k γ ( θ ) k ion. ( θ ) F ( θ ) h N s ( θ ) i ρ s, ( θ ) h N µ ( θ ) i model A ( × ) ( × ) ( × ) ( × ) ( × )eV eV eV − m − QGSJet01d p 0 .
341 2 .
846 2 .
104 2 .
178 2 .
312 5 . .
224 2 .
910 2 .
148 1 .
250 2 .
432 7 . .
364 2 .
816 2 .
070 2 .
296 2 .
438 5 . .
246 2 .
894 2 .
148 1 .
358 2 .
636 7 . .
345 2 .
822 2 .
100 2 .
512 2 .
193 4 . .
224 2 .
910 2 .
228 1 .
384 2 .
249 4 . .
377 2 .
815 2 .
023 2 .
355 2 .
655 5 . .
230 2 .
894 2 .
133 1 .
419 2 .
917 8 . .
357 2 .
825 2 .
074 2 .
335 2 .
400 5 . .
231 2 .
902 2 .
164 1 .
353 2 .
558 7 . .
294 2 .
864 2 .
119 1 .
844 2 .
479 6 . .
700 2 .
510 1 .
793 2 .
656 6 . E estimation adopted for the Yakutsk array. The observables and main components consti-tuting the primary energy are given in Tables I and II for E = 10 eV and cos θ = 0 .
95. The F column in the Table I is the flux of Cherenkov photons measured with integral Cherenkovlight detectors. The values for k γ and k ion. in the same table were obtained in simulationwith CORSIKA. Mean values of N s and N µ were obtained from the LDFs averaged overenergy interval. The row entitled “average p-Fe”corresponds to values averaged over allmodels and compositions. The energy dissipated in the atmosphere by electromagneticcomponent equals to E i = E γ + E ion. , (15)where E γ is energy of gamma-photons on observation level, E ion. — summary ionizationlosses of all electrons and positrons. It is proportional to the total flux F of Cherenkov TABLE II. Energy balance of EAS with E = 10 eV and cos θ = 0 .
95 from primary ( A ) accordingto CORSIKA [23] simulation and experiment [16]. E γ E ion. E el E µ ∆ E E model A ( × ) ( × ) ( × ) ( × ) ( × ) ( × )eV eV eV eV eV eVQGSJet01d p 0 .
806 6 .
620 1 .
469 0 .
517 0 .
565 9 . .
529 6 .
660 1 .
306 0 .
785 0 .
798 9 . .
859 6 .
476 1 .
474 0 .
547 0 .
624 9 . .
582 6 .
430 1 .
302 0 .
844 0 .
866 9 . .
909 6 .
625 1 .
523 0 .
428 0 .
491 9 . .
528 6 .
679 1 .
340 0 .
702 0 .
716 9 . .
891 6 .
412 1 .
482 0 .
524 0 .
657 9 . .
543 6 .
415 1 .
305 0 .
794 0 .
898 9 . .
866 6 .
533 1 .
487 0 .
504 0 .
584 9 . .
546 6 .
531 1 .
313 0 .
781 0 .
820 9 . .
706 6 .
532 1 .
400 0 .
643 0 .
702 9 . .
287 0 .
947 0 .
636 0 .
860 11 . .
926 0 .
947 0 .
618 0 .
702 10 . radiation in the atmosphere: E i = k · F , (16)where k (eV/photon eV − ) is the scaling factor: k = k γ + k ion. = E γ + E ion. F . (17)On Fig.3, the dependence of the scaling factor (17) from the path from x max to observationlevel ( x obs. = x × sec θ g/cm ) is shown. The flux F is determined with respect to itsattenuation by factor 1 .
15 due to Rayleigh scattering in clean atmosphere and degradationof the relative transparency in sampling events [15, 16] by factor 1 .
1. It is given for radiationinterval 1 eV: F = 1 . F obs. ∆ ǫ , (18) FIG. 3. Dependence of the scaling ratio (17) from the path between x max and x obs. for two CRcompositions. Lines represent approximations. where F obs. is the flux measured in experiment with integral Cherenkov light detectors and∆ ǫ = 12400 · (cid:18) λ − λ (cid:19) ≃ .
58 (eV). (19)Here λ = 3000 ˚A, λ = 8000 ˚A. The energy E el. is the amount of primary energy carriedby electrons and positrons to the observation level. It was estimated by integrating thedifferential energy deposit over the cascade curve N e ( x ) below the observation level x obs. : E el. = Z ∞ x obs. (cid:18) d E d x (cid:19) ion. · N e ( x )d x ≃≃ . × · N e ( x obs. ) × Z ∞ x obs. exp (cid:18) x obs. − x h λ N i (cid:19) d x , (20)where h λ N i ≃
240 g/cm . N e ( x obs. ) is the number of electrons at observational level, whichwas determined from the relation: N e ( x obs. ) ≃ h N s ( x obs. ) i − . · h N µ ( x obs. ) i , (21)where h N s ( x obs. ) i and h N µ ( x obs. ) i are mean values of the total number of responses fromall particles and muons with 1 GeV threshold, obtained by integrating of experimentally0measured corresponding LDFs [15, 16]. The ratio 1 . E µ = h E µ i · h N µ ( x obs. ) i , (22)where h E µ i = 10 . E i + E el. + E µ amounts to ≃
93 % from primary energy. The rest of it (∆ E ) is not controlled by the array. Itincludes energy of neutrinos, energy transferred to nuclei in various reactions and ionizationlosses of muons and hadrons in the atmosphere. In [15, 16] this value was obtained fromearlier calculations and is roughly consistent with predictions obtained with CORSIKA. IV. DISCUSSION
Summary values of all constituents are given in the rightmost column of the Table II.The value E = 1 . × eV in the “experiment” column exceeds the mean value h E i = 0 . × eV obtained in simulation by factor ≃ . k , occurred in [15, 16], where it was deter-mined as k = 3 . × eV/photon eV − , while simulation with CORSIKA gave h k i =3 . × eV/photon eV − .The new estimation of primary energy obtained with the use of calorimetric methoddescribed above is given in the bottom row of the Table II. The value E = 1 . × eVwas determined with corrected values E i = h k i · F , h E µ i = 10 . E . It is shownon Fig.4 together with other data from [16] with black circles. White circles represent thedata from [17, 18] reprocessed with the revised values of F and E ion. with the account ofthe adjusted atmosphere transparency and with introduction of a new scaling factor k (seeFig.3). Solid line represents the dependency: E = (3 . ± . × · ( ρ s, (0 ◦ )) . ± . , (23)which describes all the experimental data when ρ s, (18 . ◦ ) is converted to vertical with theuse of (2) with λ ρ = 415 g/cm . Dotted and dashed lines reflect the relations (11) and (12)1 FIG. 4. Energy E reconstructed from the parameter ρ s, ( θ ) in showers with h cos θ i = 0 .
95. Sym-bols represent the data given in [16] and [17, 18] reprocessed with the new calorimetric estimation.Solid line represents the best fit to all data. which signify limits of the interval containing predictions of all the abovementioned models.The closest to experiment are QGSJet-II-04 and EPOS-LHC though one cannot exclude thecredibility of two others.On Fig.5 the CR energy spectra are shown measured by modern giant EAS arrays. Circlesand squares (showers selected by master triangles with 500 m and 1000 m sides correspond-ingly) represent the data of the Yakutsk experiment [33]. Energy E , estimated with the useof expressions (4) and (5), are shown with red symbols; open symbols represent the samedata with energy estimated according to (23); green symbols — according to (13). Blacktriangles represent TA SD data [34], blue triangles — PAO [35]. V. CONCLUSION
Application of the CORSIKA code to the Yakutsk EAS array data provided an opportu-nity to critically examine the experiment’s energy calibration which for a long time has beena subject of debates and controversy among our colleagues from other world EAS arrays.2
FIG. 5. Differential energy spectrum of CR according to the data from different experiments.
This became possible thanks to the availability of modern EAS development models to awide range of researchers. With these models we have managed to calculate the responsesof scintillation detectors and obtain a set of probable estimations for primary energy (9-12). Calculations have revealed, that in relations (1) and (4) the energy dissipated in theatmosphere in the form of electromagnetic component, was overestimated by (12 − x max (see Fig.3). This was made worse in (4) due to un-derestimation of the atmosphere transparency by ≃ E in comparison with (4) by factor ≃ .
33 and in decreasedintensity of the CR energy spectrum measured on the Yakutsk EAS array (see Fig.5).Independent techniques of E estimation from SSD LDFs (13) and with the use of calori-metric method (23) gave close results, which agree with simulations within (10 − E ≥ × eV they do not contradict the TA data [34] and consistently point at thesteepening of the primary CR spectrum in the region of extreme energies ( E ≥ × eV).This steepening does not contradict to GZK cutoff but, probably, has a different astrophys-ical reason. As for the difference in spectral intensities at E ≤ × eV, it could haveother reasons. Probably, it is the effect of systematical errors in primary energy reconstruc-3tion techniques adopted by different experiments. But one cannot exclude that the saiddifference is caused by geographical locations of arrays observing different regions of the sky.In [13] such a correlation had been noticed. Our current plan is to continue the elaborationof E estimation at the Yakutsk array with a more detailed analysis of the Cherenkov lightdata. ACKNOWLEDGMENTS
This work is financially supported by Russian Academy of Science within the program“Fundamental properties of matter and astrophysics”, by RFBR grant 13-02-12036 ofi-m-2013 and by grant of President of the Republic Sakha (Yakutia) to young scientists, special-ists and students to support research. [1] K. Greisen, Phys. Rev. Lett. , 748 (1966).[2] G. T. Zatsepin and V. A. Kuz’min, JETP Lett. , 78 (1966).[3] R. U. Abbasi, T. Abu-Zayyad, M. Allen, J. F. Amman, et al. (High Resolution Fly’s EyeCollaboration), Phys. Rev. Lett. , 101101 (2008).[4] J. Abraham, P. Abreu, M. Aglietta, C. Aguirre, et al. (The Pierre Auger Collaboration), Phys.Rev. Lett. , 061101 (2008).[5] R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, et al. , Astroparticle Physics , 76 (2007).[6] E. G. Berezhko, S. P. Knurenko, and L. T. Ksenofontov, Astropart. Phys. , 31 (2012).[7] V. N. Zirakashvili and V. S. Ptuskin, Bull. of the Russ. Acad. Sci.:Phys. , 555 (2009).[8] D. M. Edge, A. C. Evans, H. J. Garmston, R. J. O. Reid, A. A. Watson, et al. , J. Phys. A:Math. Nucl. Gen. , 1612 (1973).[9] A. V. Glushkov, V. M. Grigoriev, M. N. Dyakonov, T. A. Egorov, V. P. Egorova, et al. , in Proc. of 20th ICRC, Moscow , Vol. 5 (1987) p. 494.[10] N. Sakaki, M. Chikawa, M. Fukushima, N. Hayashida, K. Honda, et al. , in
Proc. of 27th ICRC,Hamburg , Vol. 1 (2001) p. 329.[11] R. U. Abbasi, T. Abu-Zayyad, J. F. Amman, G. C. Archbold, J. A. Bellido, et al. , Astropart.Phys. , 157 (2005), arXiv:0208301 [astro-ph]. [12] J. Abraham, P. Abreu, M. Aglietta, E. J. Ahn, D. Allard, et al. (Pierre Auger Collaboration),Physics Letters B , 239 (2010), arXiv:1002.1975 [astro-ph].[13] A. V. Glushkov and M. I. Pravdin, JETP Lett. , 345 (2008).[14] B. R. Dawson, I. C. Maris, M. Roth, F. Salamida, et al. , EPJ Web of Conferences , 01005(2013).[15] A. V. Glushkov, V. M. Grigoriev, O. S. Diminshtein, N. N. Efimov, L. I. Kaganov, et al. ,“Phenomenology of extensive air showers and primary cosmic rays,” Preprint by the Instituteof cosmophysical research and aeronomy, Yakutsk (1978), A contribution for VI ECRS, Kiel.[16] A. V. Glushkov, Lateral distribution and total flux of Cherenkov radiation from EAS withprimary energy E ≥ eV , Ph.D. thesis, SINP MSU (1982), (in Russian).[17] A. V. Glushkov, M. N. Dyakonov, T. A. Egorov, N. N. Efimov, N. N. Efremov, et al. , Bull.Acad. Sci. USSR: Phys. Ser. , 713 (1991), (in Russian).[18] Egorov, T. A. (Yakutsk Group), in Proceedings of the Tokyo Workshop on Techniques for theStudy of Extremely High Energy Cosmic Rays (ICRR, U. of Tokyo, 1993) p. 35.[19] A. V. Glushkov, V. P. Egorova, A. A. Ivanov, S. P. Knurenko, V. A. Kolosov, et al. , in
Proc.of the 28th ICRC, Tsukuba , Vol. 1 (2003) p. 393.[20] S. I. Nikolsky, in
Proc. of the 5th Intern. Seminar on Cosmic Rays, La Paz , Vol. 2 (1962) pp.48–52.[21] N. Budnev, D. Chernov, O. Gress, E. Korosteleva, L. Kuzmichev, et al. , Astropart. Phys. , 18 (2013).[22] L. G. Dedenko, D. A. Podgrudkov, T. M. Roganova, G. F. Fedorova, E. Y. Fedunin, and G. P.Shozieev, Phys. of At. Nucl. , 1806 (2007).[23] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, “CORSIKA: A Monte CarloCode to Simulate Extensive Air Showers,” FZKA 6019 (Forschungszentrum Karlsruhe, 1998).[24] A. V. Glushkov, O. S. Diminshtein, T. A. Egorov, N. N. Efimov, L. I. Kaganov, et al. , in Collection of thesises (Yakutsk: Ya.Fil. SB Acad. Sci. USSR, 1974) p. 43, (in Russian).[25] N. N. Kalmykov, S. S. Ostapchenko, and A. I. Pavlov, Nucl. Phys. B - Proc. Suppl. , 17(1997).[26] S. Ostapchenko, Phys. Rev. D , 014018 (2011).[27] E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari, and T. Stanev, Phys. Rev. D , 094003 (2009). [28] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner(2013),arXiv:1306.0121v3 [hep-ph].[29] G. Battistoni, S. Muraro, P. R. Sala, F. Cerutti, A. Ferrari, et al. , in Proceedings of theHadronic Shower Simulation Workshop 2006 , Vol. 896, edited by M. Albrow and R. Raja,Fermilab (AIP Conference Proceeding, 2007) pp. 31–49.[30] A. Fass´o, A. Ferrari, J. Ranft, and P. R. Sala,
FLUKA: a multi-particle transport code (CERN-2005-10 (2005), INFN/TC 05/11, SLAC-R-773).[31] W. R. Nelson, H. Hirayama, and D. W. O. Rogers, Report SLAC 265 (Stanford Linear Accel-erator Center, 1985).[32] D. Heck and T. Pierog,
Extensive Air Shower Simulation with CORSIKA: A User’s Guide,Version 7.37xx (Karlsruhe Institut f¨ur Technologie, 2013).[33] V. P. Egorova, A. V. Glushkov, A. A. Ivanov, S. P. Knurenko, V. A. Kolosov, et al. , Nucl.Phys. B Proc. Suppl. , 3 (2004).[34] Y. Tsunesada (Telescope Array Collaboration)(2011), arXiv:1111.2507v1 [astro-ph.HE].[35] P. Abreu et al.et al.