Advances in High Energy Physics
Volume 2016 (2016), Article ID 4676024, 7 pages
http://dx.doi.org/10.1155/2016/4676024

Research Article

On the Impact of Tsallis Statistics on Cosmic Ray Showers

M. Abrahão,¹ W. G. Dantas,¹ R. M. de Almeida,¹ D. R. Gratieri,¹ and T. J. P. Penna2,3

¹EEIMVR, Universidade Federal Fluminense, Volta Redonda, RJ, Brazil
²ICEx, Universidade Federal Fluminense, Volta Redonda, RJ, Brazil
³INCT-SC, National Institute of Science and Technology, CNPq, Rio de Janeiro, RJ, Brazil

Received 9 June 2016; Accepted 25 September 2016

Academic Editor: Sally Seidel

Copyright © 2016 M. Abrahão et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The publication of this article was funded by SCOAP³.

Abstract.

We investigate the impact of the Tsallis nonextensive statistics introduced by intrinsic temperature fluctuations in -Air ultrahigh energy interactions on observables of cosmic ray showers, such as the slant depth of the maximum and the muon number on the ground . The results show that these observables are significantly affected by temperature fluctuations and agree qualitatively with the predictions of Heitler model.

1. Introduction

The Pierre Auger Observatory [1, 2] has led to great discoveries in the field of ultrahigh energy cosmic rays (UHECRs) such as the confirmation of suppression of the cosmic ray flux at energies above  eV [3], first observed by the HiRes Collaboration [4], limits on photon [5–7], and neutrino [8–11] fluxes at ultrahigh energies and a hint of large scale anisotropies at energies above  eV [12]. Nevertheless, many questions related to these particles are still open. Particularly interesting is the behavior of the slant depth of the shower maximum with energy. Understood in terms of the LHC-tuned shower models, the HiRes [13] and Telescope Array Collaborations [14] reported a light mass composition at energies above  eV, while the Auger results suggest a gradual shift to a heavier composition, with a large fraction of protons at  eV, changing to a heavier composition at  eV [15]. However, we should interpret these results with caution, since measurements of shower properties performed by the Auger Collaboration have revealed inconsistencies between data and present shower models. For instance, the Pierre Auger Collaboration has reported the first hybrid measurement of the average muon number in inclined air showers at ultrahigh energies, suggesting a muon deficit in simulations of about to at  eV, depending on the hadronic interaction model [16]. Hence, the measured behavior of the slant depth of the shower maximum evolution could be understood as a hint of new hadronic interaction physics at energy scales beyond the reach of LHC.

In this work, we will deal with hadronic interactions in a statistical model, as first introduced by Hagedorn [17] ideas in the sixties. Recently a power-law function based on the Tsallis statistics [18] has been widely used in fitting the transverse momentum () and pseudorapidity () distributions measured in high energy collisions [19–24], while several studies have been devoted to discuss these results in the literature [25–35]. The Tsallis statistics, which is frequently present to model different branches of science, is often used to describe systems which display properties like memory effects, long range interactions, intrinsic fluctuations, (multi)fractal phase space, and so on. It consists in replacing the classical Boltzmann-Gibbs entropy by the form proposed by Tsallis:

where is the probability of a particle occupying the state and is the Tsallis index. This definition comprises the Boltzmann-Gibbs entropy as a particular case, where . On the other hand, a straight consequence from this expression is that the generalized entropy is no longer an extensive quantity, once we can verify that

with the parameter being a measure of the nonextensivity of the system. As a consequence, we must replace the usual exponential Boltzmann-Gibbs distribution, , by the Tsallis power-law distribution:

where is the state energy and is the temperature of the system.

According to [35], the behaviors presented by the transverse momentum and pseudorapidity distributions, in high energies domain, are best described using a nonexponential distribution, such as the one proposed by Tsallis. In fact, following the ideas discussed in [35], that behavior emerges from fluctuations of the thermal energy within the gas of quarks and gluons before the hadronization process. Using this approach, we can relate the parameter with those thermal fluctuations:

in which is the variance of the temperature. Obviously, when , we recover the expected result obtained in the Boltzmann-Gibbs description, where we get an equilibrium at temperature .

By assuming such scenario, in which the temperature fluctuates within each collision, the energy distribution of the particles generated in a single high energy interaction follows a power-law Tsallis distribution, given by (3). Figure 1(a) presents the Tsallis energy distribution with a fixed temperature for different values of . We can see that as values become higher, the probability for generating particles with larger energy values becomes greater. As a consequence of the total energy conservation constraint, , where is the total energy of the interaction in the center of momentum frame; it can be shown that the Tsallis statistics leads to a negative binomial multiplicity distribution given by

Such distribution has a form shown in Figure 1(b), where it is possible to see how its maximum is affected by , becoming closer to zero as grows. The value of the temperature  GeV chosen for both plots of Figure 1 is the same as that used in the simulations described in Section 3 and, as will be discussed later, corresponds to the mean energy available (at center of mass frame and inelasticity = 1) per produced particle in a collision between a proton of  eV, in the lab frame, and a nucleus of the atmosphere, according to the Sibyll hadronic interaction model. Besides, the inset plot of this figure shows that the relationship between the value for the maximum of the multiplicity distribution and is quite linear, at least in that domain of values. Therefore, one can see that the introduction of the Tsallis statistics in this context changes the energy, momenta, and multiplicity distributions of the particles generated in the hadronic interaction.

(a)

(b)

Figure 1: Tsallis energy distributions (a) and the corresponding multiplicity distributions (b) as a function of for a temperature  GeV. The mean multiplicity of all distributions presented in (b) is , according to the prediction of Sibyll model.

The transverse momentum and pseudorapidity distributions resulting from high energy collisions measured by several experiments show a large discrepancy in the values of the parameter , reflecting different physics for the transverse and the longitudinal space. The transverse distributions are thermal-like, presenting a parameter almost independent of the energy, while those from the longitudinal space have a temperature sensitive to the energy of the collision, understood as the mean energy available per produced particle [35], , where and are, respectively, the mean multiplicity and inelasticity of the interaction. Moreover, since the measured Tsallis index for the longitudinal space is much larger than that measured for the transverse space and , resulting . Also, as verified by simulations, the transverse momentum distribution has a minor contribution on the cosmic ray observables studied in this work. Therefore, from now on, we will assume a statistical equilibrium for the transverse momentum space and we will refer to the entropic index and temperature .

The goal of the present paper is to study the impact of temperature fluctuations, represented by the parameter , on the shower maximum, , and number of muons on the ground, . The simulations performed in this work are described in Section 2. In Section 3, we present the results of the simulations and discuss them in light of the Heitler model. Finally, we present the conclusions of this work in Section 4.

2. Simulations

For the simulations presented in this work, we have used CORSIKA 7.40 [36] with the interaction models Sibyll 2.1 [37] and GHEISHA 2002d [38] for high and low energy processes, respectively. The muon energy threshold used in the simulations is 0.3 GeV and the array detector position is at 1400 m above sea level, corresponding to the mean altitude of the Pierre Auger Observatory. The air shower simulation chain is as follows: first we simulate the secondaries generated in the collision between a cosmic ray and a nucleus of the upper atmosphere externally by assuming that the hadronization process is described by the Tsallis statistics; the resulting particle list is then inserted back into CORSIKA (using the stacking option and sampling option with thinning = ) to proceed with usual cascade development through the atmosphere. Such a procedure was performed 1000 times for each of several values of (1.01, 1.025, 1.05, 1.075, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, and 1.50) for showers with zenith angle initiated by a proton of fixed energy  eV. The reason for limiting the entropic index to in this work is that the mean value of the Tsallis distribution , given by , is well defined only for [31]. This model assumes that the Sibyll predictions are valid for lower energies, since they are tuned by accelerator data, while they fail for higher energies. This added to the parametrizations of the LHC transverse momentum and pseudorapidity distributions by the Tsallis statistics justifies its use in this work for energies above  eV. The point of the air shower first interaction is determined using the -Air cross section predicted by the Epos 1.99 model [39]. The reason for using this value instead of the one predicted by the Sibyll model is that the latter presents a large discrepancy in relation to that measured by the Pierre Auger Collaboration [40]. The mean multiplicity and inelasticity distribution of the -Air interaction used in this work were extracted from -Air interaction simulations using the Sibyll model. Since the Tsallis distribution is nonextensive, generating particle energies according to this distribution subject to the constraint is not a simple task, because the probabilities associated with each particle do not factorize [41]. Therefore, we perform the simulation process according to the following procedure: first, we select the number of particles generated in the -Air interaction using the expression given by (5) and we assign an energy to each particle according to the Tsallis distribution. After that, we normalize the energies such as ; that is, we multiply each energy by , where before the normalization. Then two particles and are randomly selected and a random fraction of the energy is given to particle in such a way that the new values of energies are and . After that, we compute the deviation of this energy distribution in relation to Tsallis distribution using estimator defined by

If the new value is smaller than the previous one, we accept the changes in energy of the particles and ; otherwise we cancel them. We keep repeating this procedure for another pair of particles. The whole process continues until value is stabilized. Generally, it takes iterations to reach such stabilization. To be conservative, the simulations presented in this work were performed using iterations for each -Air interaction. Since we verified through simulations that and are not sensitive to changes in , all simulations corresponding to the transverse space presented in this work were evaluated assuming a statistical equilibrium, that is, the Hagedorn [17] transverse momentum distribution:

with  MeV. The types of particles are randomly generated according to Sibyll predictions and once we have generated the particle masses , the longitudinal momentum is obtained as . These kinematic variables, along with the species of particles, complete all the information we need to reintroduce the secondary particles into CORSIKA and proceed with the simulation of the shower propagation through the atmosphere.

3. Results and Discussion

In the following, we describe the impact of fluctuations, represented by the parameter , on the shower maximum, , and number of muons on the ground. We will discuss them in terms of the predictions of the Heitler model [42, 43] and the results achieved in [44]. Although extremely simple, the predictions of the Heitler model are remarkable. The Heitler model assumes that the shower maximum is reached when the energies of particles become smaller than a critical energy, in which energy loss processes dominate the production of new particles in the case of electromagnetic component, or the charged pion interaction length becomes larger than the decay length of pions into muons, in the case of the hadronic one. As a consequence, the Heitler model predicts an increase of for smaller mean multiplicities, since larger multiplicities correspond to lower energy per particle. Besides, [44] describes a detailed investigation of the impact of the multiplicity, hadronic particle production cross section, elasticity, and pion charge-ratio on air shower observables with most of the predictions qualitatively understood within the simple Heitler model and its extension to hadronic component.

The Pearson coefficient was used to assess the degree of correlation between air showers observables and . It is defined by

and measures the linear correlation between two variables and , yielding a value in the interval , with meaning total positive correlation, meaning no correlation, and meaning total negative correlation. is the covariance between and and and are the standard deviations of variables and . The results for the mean depth of shower maximum, , and the fluctuations of are summarized in Figure 2.

(a)

(b)

Figure 2: (a) and (b) distributions obtained from air shower simulations initiated by the interaction between a proton of  eV and zenith angle and a nucleus of the upper atmosphere in which the hadronization process is described by the Tsallis statistics.

Figure 2(a) shows as a function of . A strong correlation is observed yielding a Pearson correlation coefficient . The red line presents the best linear fit () corresponding to reduced . The comparison between the predictions from Heitler model and [44] with our results requires caution, since we did not change the mean multiplicity in our simulations of the first interaction. However, the changes in distributions of energy and momenta of the particles generated in the first interaction result in a spread of the multiplicity distribution and a shift of its peak to lower values as is shown in Figure 1. For illustration, we also show in the same picture the multiplicity distribution obtained from simulations using the Sibyll model. The reason for the strong correlation of and is that most of the showers generated with larger values are initiated with smaller multiplicities, or equivalently, with higher energy per particle.

Besides, Figure 2(b) presents the corresponding plot for as a function of . In this case, the observed correlation is not strong, with and corresponding to the best linear fit, shown by the red line. According to the Heitler model, the variance of distribution depends on the mean free path of -Air interaction, , and multiplicity, , via , where  g/cm² is the electromagnetic radiation length. Although increases with , the observed correlation between and is weak, since the increase of spread of the multiplicity distribution for larger values is dominated by the first term contribution in as a consequence of the relatively large value of . For example, using the -Air mean free path corresponding to the -Air cross section from the Epos 1.99 model,  gcm⁻²,  gcm⁻², and .

The impact of on the mean number of muons in the ground, , and the fluctuations of are summarized in Figure 3. as a function of , shown in Figure 3(a), presents a strong anticorrelation with , with and reduced corresponding to the best linear fit, marked in red. A superficial analysis of Figures 2(a) and 3(a) could indicate wrongly that and are anticorrelated. The positive correlation between and position exists but it is weak, since muons are hardly attenuated in the atmosphere. Therefore, it is not the most important factor for behavior as a function of . Indeed muons are mainly produced as a result of pion decay and their abundance in the ground, especially considering the most energetic particles, is strongly correlated with the number of pions in the shower. As a consequence of the reduction of the peak of the multiplicity distribution for larger values, most showers present lower production of pions in the first interaction, constituting the main reason for the observed anticorrelation between and .

(a)

(b)

Figure 3: (a) and (b) distributions obtained from air shower simulations initiated by the interaction between a proton of  eV and and a nucleus of the upper atmosphere in which the hadronization process is described by the Tsallis statistics.

On the other hand, Figure 3(b) shows a strong correlation of and , with , as a natural consequence of the spread of the multiplicity distribution. The red line shows the best linear fit with corresponding . The use of the Sibyll cross section instead of the one from Epos 1.99 model would reduce the mean and RMS of the first interaction depth in ~10 gcm⁻², with the same impact on . On the other hand, differently from electrons, which present a fast absorption at larger depth, the muon profile only exhibits a very moderate absorption after the maximum, and, therefore, the change in the first interaction depth is not significant in relation to the number of muons in the ground. Nevertheless, although there are slightly differences in the values of the observables, it is important to notice that the correlations coefficients obtained in this paper are not affected by such a change in the hadronic interaction cross section.

Finally, Table 1 summarizes the correlation coefficients between air shower observables and obtained in this work by using the Sibyll model. For comparison, it also shows the correlation coefficients obtained after repeating the analysis with QGSJET-II hadronic interaction model [45] instead of Sibyll. We can see that the correlations coefficients are not affected significantly by the change of the hadronic interaction model.

Table 1: Summary of the correlation coefficients between air shower observables and obtained in this work with the Sibyll and QGSJET-II models.

(Sibyll)

(QGSJET-II)

4. Conclusions

Although the simulations presented in this work are a very simple description of ultrahigh energy interactions, the results presented here show that intrinsic fluctuations of the system with respect to , given by the parameter , change the energy, momenta, and multiplicity distributions of the particles generated in the interaction between a cosmic ray and a nucleus of the atmosphere, impacting air shower observables such as the slant depth of the maximum and the muon number on the ground . The results show that the higher the temperature fluctuations, the greater the values of the mean slant depth of maximum and variance of the number of muons on the ground , with Pearson correlation coefficients of and , respectively. This results from the spread and shift of the maximum of the multiplicity distribution to lower values for larger temperature fluctuations. Besides, as muons are mainly produced by the decay of charged pions and the shift in the peak of multiplicity distribution reduces the number of such particles generated in the first interaction, the mean number of muons on the ground presents a negative correlation with , producing . On the other hand, the variance of the slant depth distribution presents a weak correlation with the temperature fluctuations, with , because the contribution of the hadronic -Air interaction fluctuations dominates the one originated from the multiplicity distribution. These results agree qualitatively with the Heitler model and [44] predictions. The correlation coefficients are not affected by changing the Sibyll by the GQSJET-II hadronic interaction model. Although the results presented in this paper have been obtained for a specific nonextensive hadronic interaction model, we believe that they capture the essential features related to the presence of Tsallis statistics in UHECR showers and can shed light on the understanding of their properties in addition to particle interactions at these energies. Studies regarding different nonextensive particle interaction models are out of the scope of this work and will be addressed in a future work.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Acknowledgments

The authors thank Professor J. C. Anjos for useful discussions in the beginning of this work. This work was partially supported by the “Conselho Nacional de Desenvolvimento Científico e Tecnológico” (CNPq) and by the “Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro” (FAPERJ).

References

J. Abrahama, M. Agliettab, and I. C. Aguirre, “Properties and performance of the prototype instrument for the Pierre Auger Observatory,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 523, no. 1-2, pp. 50–95, 2004.

Pierre Auger Collaboration, A. Aab et al., “The pierre auger cosmic ray observatory,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 798, pp. 172–213, 2015.

J. Abraham, P. Abreu, M. Aglietta et al., “Observation of the suppression of the flux of cosmic rays above  eV,” Physical Review Letters, vol. 101, no. 6, Article ID 061101, 2018.

R. U. Abbasi, T. Abu-Zayyad, M. Allen et al., “First observation of the Greisen-Zatsepin-Kuzmin suppression,” Physical Review Letters, vol. 100, no. 10, 5 pages, 2008.

J. Abraham, M. Aglietta, C. Aguirre et al., “An upper limit to the photon fraction in cosmic rays above 10¹⁹ eV from the Pierre Auger Observatory,” Astroparticle Physics, vol. 27, no. 2-3, pp. 155–165, 2007.

J. Abraham, P. Abreu, M. Aglietta et al., “Upper limit on the cosmic-ray photon fraction at EeV energies from the Pierre Auger Observatory,” Astroparticle Physics, vol. 31, no. 6, pp. 399–406, 2009.

J. Abraham, P. Abreu, M. Aglietta et al., “Upper limit on the cosmic-ray photon flux above 10¹⁹ eV using the surface detector of the Pierre Auger Observatory,” Astroparticle Physics, vol. 29, no. 4, pp. 243–256, 2008.

J. Abraham, P. Abreu, M. Aglietta et al., “Upper limit on the diffuse flux of ultrahigh energy tau neutrinos from the Pierre Auger Observatory,” Physical Review Letters, vol. 100, no. 21, Article ID 211101, 7 pages, 2008.

P. Abreu, M. Aglietta, M. Ahlers et al., “A search for ultra-high energy neutrinos in highly inclined events at the Pierre Auger Observatory,” Physical Review D, vol. 84, Article ID 122005, 2011.

P. Abreu, M. Aglietta, M. Ahlers et al., “Search for point-like sources of ultrahigh energy neutrinos at the Pierre Auger Observatory and improved limit on the diffuse flux of tau neutrinos,” The Astrophysical Journal, vol. 755, p. L4, 2012.

P. Abreu, M. Aglietta, M. Ahlers et al., “Ultrahigh energy neutrinos at the Pierre Auger Observatory,” Advances in High Energy Physics, vol. 2013, Article ID 708680, 18 pages, 2013.

A. Aab, P. Abreu, and M. Aglietta, “Large scale distribution of ultra high energy cosmic rays detected at the Pierre Auger Observatory with zenith angles up to 80°,” The Astrophysical Journal, vol. 802, no. 2, p. 111, 2015.

R. U. Abbasi, T. Abu-Zayyad, G. Archbold et al., “A study of the composition of ultra-high-energy cosmic rays using the high-resolution Fly's eye,” The Astrophysical Journal, vol. 622, no. 2, p. 910, 2005.

R. U. Abbasi, M. Abe, T. Abu-Zayyad et al., “Study of Ultra-High Energy Cosmic Ray composition using Telescope Array’s Middle Drum detector and surface array in hybrid mode,” Astroparticle Physics, vol. 64, pp. 49–62, 2015.

A. Aab, P. Abreu, M. Aglietta et al., “Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above 1017.8 eV,” Physical Review D, vol. 90, no. 12, Article ID 122005, 25 pages, 2014.

A. Aab, P. Abreu, M. Aglietta et al., “Muons in air showers at the Pierre Auger Observatory: mean number in highly inclined events,” Physical Review D, vol. 91, Article ID 032003, 2015, Erratum in: Physical Review D, vol. 91, Article ID 059901, 2015

R. Hagedorn, “Statistical thermodynamics of strong interactions at high energies,” Supplemento Al Nuovo Cimento, vol. 3, pp. 147–186, 1965.

C. Tsallis, “Possible generalization of Boltzmann-Gibbs statistics,” Journal of Statistical Physics, vol. 52, no. 1-2, pp. 479–487, 1988.

B. I. Abelev, J. Adams, M. M. Aggarwal et al., “Strange particle production in p+p collisions at  GeV,” Physical Review C, vol. 75, Article ID 064901, 2007.

A. Adare, S. Afanasiev, C. Aidala et al., “Measurement of neutral mesons in collisions at  GeV and scaling properties of hadron production,” Physical Review C, vol. 83, no. 5, Article ID 052004, 26 pages, 2011.

A. Adare, S. Afanasiev, C. Aidala et al., “Identified charged hadron production in collisions at and 62.4 GeV,” Physical Review C, vol. 83, no. 6, Article ID 064903, 29 pages, 2011.

K. Aamodt, N. Abel, U. Abeysekara et al., “Production of pions, kaons and protons in pp collisions at  GeV with ALICE at the LHC,” The European Physical Journal C, vol. 71, aricle 1655, 2011.

G. Aad, B. Abbott, J. Abdallah et al., “Charged-particle multiplicities in pp interactions measured with the ATLAS detector at the LHC,” New Journal of Physics, vol. 13, Article ID 053033, 2011.

V. Khachatryan, A. M. Sirunyan, A. Tumasyan et al., “Strange particle production in pp collisions at and 7 TeV,” Journal of High Energy Physics, vol. 2011, no. 5, pp. 1–40, 2011.

J. Cleymans, G. I. Lykasov, A. S. Parvan, A. S. Sorin, O. V. Teryaev, and D. Worku, “Systematic properties of the Tsallis distribution: energy dependence of parameters in high energy p-p collisions,” Physics Letters, Section B: Nuclear, Elementary Particle and High-Energy Physics, vol. 723, no. 4-5, pp. 351–354, 2013.

J. Cleymans, “The tsallis distribution for p-p collisions at the LHC,” Journal of Physics: Conference Series, vol. 455, no. 1, Article ID 012049, 2013.

M. D. Azmi and J. Cleymans, “Transverse momentum distributions in p-Pb collisions and Tsallis thermodynamics,” https://arxiv.org/abs/1311.2909

M. D. Azmi and J. Cleymans, “Transverse momentum distributions in proton-proton collisions at LHC energies and Tsallis thermodynamics,” Journal of Physics G: Nuclear and Particle Physics, vol. 41, no. 6, Article ID 065001, 2014.

C.-Y. Wong and G. Wilk, “Tsallis fits to pT spectra and multiple hard scattering in pp collisions at the LHC,” Physical Review D-Particles, Fields, Gravitation and Cosmology, vol. 87, no. 11, Article ID 114007, 2013.

C. Y. Wong and G. Wilk, “Tsallis fits to spectra for pp collisions at the LHC,” Acta Physica Polonica B, vol. 43, pp. 2047–2054, 2012.

G. Wilk and Z. Włodarczyk, “Interpretation of the nonextensivity parameter in some applications of Tsallis statistics and Lévy distributions,” Physical Review Letters, vol. 84, no. 13, pp. 2770–2773, 2000.

L. Marques, J. Cleymans, and A. Deppman, “Description of high-energy pp collisions using Tsallis thermodynamics: transverse momentum and rapidity distributions,” Physical Review D—Particles, Fields, Gravitation and Cosmology, vol. 91, no. 5, Article ID 054025, 2015.

B.-C. Li, Y.-Z. Wang, F.-H. Liu, X.-J. Wen, and Y.-E. Dong, “Particle production in relativistic and AA collisions at RHIC and LHC energies with Tsallis statistics using the two-cylindrical multisource thermal model,” Physical Review D, vol. 89, no. 5, Article ID 054014, 2014.

B.-C. Li, Z. Zhang, J.-H. Kang, G.-X. Zhang, and F.-H. Liu, “Tsallis statistical interpretation of transverse momentum spectra in high-energy P A collisions,” Advances in High Energy Physics, vol. 2015, Article ID 741816, 10 pages, 2015.

G. Wilk and Z. Włodarczyk, “Power laws in elementary and heavy-ion collisions: a story of fluctuations and nonextensivity?” The European Physical Journal A, vol. 40, no. 3, pp. 299–312, 2009.

D. Heck, G. Schatz, T. Thouw, J. Knapp, and J. N. Capdevielle, “CORSIKA: a Monte Carlo code to simulate extensive air showers,” Report FZKA 6019, Forschungszentrum Karlsruhe, Karlsruhe, Germany, 1998, http://www-ik.fzk.de/corsika/physics_description/corsika_phys.html.

E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari, and T. Stanev, “Cosmic ray interaction event generator SIBYLL 2.1,” Physical Review D—Particles, Fields, Gravitation and Cosmology, vol. 80, no. 9, Article ID 094003, 2009.

H. Fesefeldt, “The simulation of hadronic showers: physics and applications,” Report PITHA-85/02, RWTH Aachen, 1985.

K. Werner, F.-M. Liu, and T. Pierog, “Parton ladder splitting and the rapidity dependence of transverse momentum spectra in deuteron-gold collisions at the BNL relativistic heavy ion collider,” Physical Review C, vol. 74, no. 4, Article ID 044902, 2006.

P. Abreu, M. Aglietta, E. J. Ahn et al., “Measurement of the proton-air cross section at  TeV with the Pierre Auger Observatory,” Physical Review Letters, vol. 109, no. 6, Article ID 062002, 9 pages, 2012.

C. Beck, “Non-extensive statistical mechanics and particle spectra in elementary interactions,” Physica A: Statistical Mechanics and its Applications, vol. 286, no. 1, pp. 164–180, 2000.

W. Heitler, “Theory of meson production,” Reviews of Modern Physics, vol. 21, no. 1, pp. 113–121, 1949.

J. Matthews, “A Heitler model of extensive air showers,” Astroparticle Physics, vol. 22, no. 5-6, pp. 387–397, 2005.

R. Ulrich, R. Engel, and M. Unger, “Hadronic multiparticle production at ultrahigh energies and extensive air showers,” Physical Review D, vol. 83, no. 5, Article ID 054026, 2011.

S. Ostapchenko, “QGSJET-II: towards reliable description of very high energy hadronic interactions,” Nuclear Physics B, vol. 151, no. 1, pp. 143–146, 2006.