On the Impact of Tsallis Statistics on Cosmic Ray Showers

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


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 4 × 10 19 eV [3], first observed by the HiRes Collaboration [4], limits on photon [5][6][7], and neutrino [8][9][10][11] fluxes at ultrahigh energies and a hint of large scale anisotropies at energies above 8 × 10 18 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 10 18 eV, while the Auger results suggest a gradual shift to a heavier composition, with a large fraction of protons at 10 18 eV, changing to a heavier composition at 10 19.5 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 30% to 80 +17 −20 (sys)% at 10 19 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][20][21][22][23][24], while several studies have been devoted to discuss these results in the literature [25][26][27][28][29][30][31][32][33][34][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 ( BG ) 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  = 1.On the other hand, a straight consequence from this 2 Advances in High Energy Physics 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, exp(−/), 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  2  is the variance of the temperature.Obviously, when  = 1, 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, ∑ = =1   =  CM , where  CM 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  = 462 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 10 18 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.
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],   =  CM /⟨⟩, 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,  max , 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.

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 = 10 −6 ) 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  = 38 ∘ initiated by a proton of fixed energy  = 10 18 eV.The reason for limiting the entropic index to  = 1.5 in this work is that the mean value of the Tsallis distribution (), given by ⟨⟩ = /(3 − 2), is well defined only for 1 ≤  ≤ 1.5 [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  ∼ 10 18 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 ∑ = =1   =  CM 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 ∑ = =1   =  CM ; that is, we multiply each energy   by  CM / tot , where  tot = ∑ = =1   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  ,new =   − Δ  and  ,new =   + Δ  .After that, we compute the deviation of this energy distribution in relation to Tsallis distribution () using  2 estimator defined by If the new  2 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  2 value is stabilized.Generally, it takes 10 4 iterations to reach such stabilization.To be conservative, the simulations presented in this work were performed using 10 5 iterations for each -Air interaction.Since we verified through simulations that  max 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 ⟨  ⟩ = 133 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   = √ 2 −  2 −  2  .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.

Results and Discussion
In the following, we describe the impact of   fluctuations, represented by the parameter  =   , on the shower Figure 2:  max (a) and ( max ) (b) distributions obtained from air shower simulations initiated by the interaction between a proton of  = 10 18 eV and zenith angle  = 38 ∘ and a nucleus of the upper atmosphere in which the hadronization process is described by the Tsallis statistics.
maximum,  max , 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 ⟨ max ⟩ 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 chargeratio 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 [−1, +1], with 1 meaning total positive correlation, 0 meaning no correlation, and −1 meaning total negative correlation.cov() is the covariance between  and  and () and () are the standard deviations of variables  and .The results for the mean depth of shower maximum, ⟨ max ⟩, and the fluctuations of  max are summarized in Figure 2. Figure 2(a) shows ⟨ max ⟩ as a function of .A strong correlation is observed yielding a Pearson correlation coefficient  = 0.90.The red line presents the best linear fit (() = 0 + 1 ⋅ ) corresponding to reduced  2 = 0.64.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 ⟨ max ⟩ 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 ( max ) as a function of .In this case, the observed correlation is not strong, with  = 0.40 and  2 = 2.04 corresponding to the best linear fit, shown by the red line.According to the Heitler model, the variance of  max distribution depends on the mean free path of -Air interaction,   , and multiplicity, , via ( max ) ∝  2   + ( 0 ln 2) 2 (ln ), where  0 ∼ 37 g/cm 2 is the electromagnetic radiation length.Although (ln ) increases with , the observed correlation between ( max ) and  is weak, since the increase of spread of the multiplicity distribution for larger  values is dominated by the first term contribution in ( max ) 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,   ∼ 48 gcm −2 ,  0 ln 2 ∼ 26 gcm −2 , and (ln ) ∼ 1.
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  = −0.67 and reduced  2 = 2.02 corresponding to the best linear fit, marked in red.A superficial analysis of Figures 2(a    and 3(a) could indicate wrongly that ⟨  ⟩ and ⟨ max ⟩ are anticorrelated.The positive correlation between ⟨  ⟩ and ⟨ max ⟩ 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 .
On the other hand, Figure 3(b) shows a strong correlation of (  ) and , with  = 0.81, as a natural consequence of the spread of the multiplicity distribution.The red line shows the best linear fit with corresponding  2 ] = 0.96.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 ⟨ max ⟩.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.

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  max 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 ⟨ max ⟩ and variance of the number of muons on the ground (  ), with Pearson correlation coefficients of  = 0.90 and  = 0.81, 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  = −0.67.On the other hand, the variance of the slant depth distribution ( max ) presents a weak correlation with the temperature fluctuations, with  = 0.40, 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.

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

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

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