#### 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 thatwith 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 bySuch 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)**

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 byIf 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 byand 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) 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^{2} 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^{−2}, gcm^{−2}, 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)**

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^{−2}, 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.

#### 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).