# Hadronic interactions models beyond collider energies

## I Introduction

Extremely high energy cosmic rays (CR) are an extravagant phenomenon of nature that has baffled astrophysics for more than three decades . Ingenious installations with large effective area and long exposure times to overcome the rapidly decreasing flux, $``$ 1 event per km<sup>2</sup> per year (century) at $`10^{19}`$ ($`10^{20}`$) eV, are required to study them. Their energy spectrum beyond 1 PeV needs to be studied indirectly through the extensive air showers (EAS) they produce deep in the atmosphere. Thus, the interpretation of the observed cascades generally depends on Monte Carlo simulations which extrapolate hadronic interaction models to energies well beyond those explored at acellerators.

There is a couple of quite elaborate models (the dual parton model (DPM)  and the quark gluon string (QGS) model of the supercritical Pomeron ) that provide a complete phenomenological description of all facets of soft hadronic processes. These models, inspired on $`1/N`$ expansion of QCD are also supplemented with generally accepted theoretical principles like duality, unitarity, Regge behavior and parton structure. At higher energies, however, there is evidence of minijet production  and correlation between multiplicity per event and transverse momentum per particle , suggesting that semihard QCD processes become important in high energy hadronic interactions. It is precisely the problem of a proper accounting for semihard processes the major source of uncertainty of extensive air showers event generators.

Two codes of hadronic interactions with similar underlying physical assumptions and algorithms tailored for efficient operation to the highest cosmic ray energies are SIBYLL  and QGSJET . In these codes, the low $`p_\mathrm{T}`$ interactions are modeling by the exchange of Pomerons. Regge singularities are used to determine the momentum distribution functions of the various sets of constituents, valence and sea quarks. In the interaction the hadrons exchange very soft gluons simulated by the production of a single pair QCD strings and the subsequent fragmentation into colour neutral hadrons. In QGSJET these events also involve exchange of multiple pairs of soft strings.

As above mentioned, the production of small jets is expected to dominate interactions in the c.m. energy above $`\sqrt{s}`$ 40 TeV. The underlying idea behind SIBYLL is that the increase in the cross section is driven by the production of minijets . The probability distribution for obtaining $`N`$ jet pairs (with $`p_\mathrm{T}^{\mathrm{jet}}>p_\mathrm{T}^{\mathrm{min}}`$, being $`p_\mathrm{T}^{\mathrm{min}}`$ a sharp threshold on the transverse momentum below which hard interactions are neglected) in a collision at energy $`\sqrt{s}`$ is computed regarding elastic $`pp`$ or $`p\overline{p}`$ scattering as a difractive shadow scattering associated with inelastic processes . The algorithms are tuned to reproduce the central and fragmentation regions data up to $`p\overline{p}`$ collider energies, and with no further adjustments they are extrapolated several orders of magnitude.

In QGSJET the theory is formulated entirely in terms of Pomeron exchanges. The basic idea is to replace the soft Pomeron by a so-called “semihard Pomeron”, which is defined to be an ordinary soft Pomeron with the middle piece replaced by a QCD parton ladder. Thus, minijets will emerge as a part of the “semihard Pomeron”, which is itself the controlling mechanism for the whole interaction. After performing the energy sharing among the soft and semihard Pomerons, and also the sharing among the soft and hard pieces of the last one; the number of charged particles in the partonic cascade is easily obtained generalizing the method of multiple production of hadrons as discussed in the QGS model (soft Pomeron showers) .

The most outstanding point in connection with the shower development is certainly the incorporation of nuclear effects. Both, SIBYLL and QGSJET  describe particle production in hadron-nucleus collisions in a quite similar fashion. The high energy projectile undergoes a multiple scattering as formulated in Glauber’s approach , particle production comes again after the fragmentation of colorless parton-parton chains constructed from the quark content of the interacting hadrons. In cases with more than one wounded nucleon in the target the extra strings are connected with sea-quarks in the projectile. This ensures that the inelasticity in hadron-nucleus collisions is not much larger than that corresponding to hadron-hadron collisions. A higher inelastic nuclear stopping power yields relatively rapid shower developments which are ruled out by $`p`$-nucleus data .

Except by the depth of maximum of the shower – a quantity which is well known to depend on the location and features of the first interaction – one can expect that global observables should not be affected by the above reasonable alternative physics assumptions. However, the question on the sensitivity of the free parameters of these models (which have been derived from available accelerator data) when they are extrapolated to energies essentially greater than attained with colliders, is surely an interesting one. To answer this question is the main goal of the present article.

In this work we shall present several comparative studies between SIBYLL and QGSJET . The outline of the paper is as follows. In Sec. II we proceed by first analyzing the different predictions on $`p`$-air and $`\pi `$-air cross sections, and then we compare single $`p`$-$`\overline{p}`$ and $`\overline{p}`$-nuclei hadronic interactions. In Sec. III we present results of several numerical analyses. Around 5000 air showers induced by protons with energies ranging from $`10^{14}`$ up to $`10^{20.5}`$ eV are generated with the code AIRES , a realistic air shower simulation system which includes electromagnetic interactions algorithms  and links to the mentioned SIBYLL and QGSJET models. We conclude in Sec. IV with the final remarks.

## II Hadronic collisions

Let us start our comparative analysis of SIBYLL and QGSJET by discussing briefly the $`p`$-air and $`\pi `$-air cross sections as calculated by these models and also by the simulation program AIRES.

In Fig. 1 (2) the $`p`$-air ($`\pi `$-air) cross section is plotted versus the projectile laboratory energy. In the cases of AIRES cross sections (solid lines) –which are equivalent to the so-called “Bartol cross sections”– and QGSJET (dotted lines), the mentioned laboratory energy is the input energy for the corresponding algorithms. On the other hand, the input parameter for SIBYLL is the c.m. energy in the hadron-nucleon system . Notice that similar plots can be found in Ref.  with a different behaviour of the SIBYLL cross sections. We attribute these differences to the c.m. conversion procedure. In fact, if the laboratory energy is (mistakenly) converted to the hadron-nucleus system, one may reproduce the data of Ref. .

Note that in QGSJET the growth in the cross section is fitted using a mixture of soft and hard interactions, contrariwise, SIBYLL does so just with the hard processes (minijets). Thus, it is feasible to expect a deviation in the predictions when the algorithms are extrapolated several orders of magnitude. Furthermore, SIBYLL predictions for the $`p`$-$`\overline{p}`$ cross section ought to be higher than the ones of QGSJET since hard processes overrule the soft ones with the rise of energy. As aforementioned, the extension to hadron nuclei interactions is computed in both codes in the framework of Glauber theory with minute differences, yielding no significant additional divergences. As a consequence, we attribute the different behaviours for the cross sections shown in Figs. 1 and 2, to the way in which the free parameters of both codes are fitted to reproduce $`p`$-$`\overline{p}`$ collider data. Namely, $`p_{\mathrm{min}}^2=5`$ GeV<sup>2</sup> and the multiplicative ad hoc factor $`k=1.7`$ in SIBYLL code.<sup>*</sup><sup>*</sup>*Actually one still has another parameter but hardly would have any influence at high energies. Recall that in SIBYLL the soft part of the eikonal function is taken as a constant fitted to low energy data . On the other hand, in QGSJET we have: i) parameters of the Pomeron trajectory: $`\mathrm{\Delta }=0.07`$, $`\alpha _\mathrm{P}^{}(0)=0.21`$ GeV<sup>-2</sup>, ii) the ones from the Pomeron vertices: $`R_{pp}^2=3.56`$ GeV<sup>-2</sup>, $`\gamma _p^2`$ = 3.64 GeV<sup>-2</sup>, iii) the so-called “shower enhancements coefficient $`C_{pp}=1.5`$, and iv) the parameters of semihard processes: $`p_{\mathrm{min}}^2`$ = 4 GeV<sup>2</sup> (this parameter as in the SIBYLL code represent the threshold of hard interactions), the parameter associated with the parton density $`r^2=0.6`$ GeV<sup>-2</sup>; for a survey the reader is referred again to  and references therein. It is interesting to remark that with an alternative value in the Pomeron trajectory parameter one can reproduce the cross section without hard processes, viz. with QGS model (see Table I of ).

The most direct way to analyze the differences between the models is to study the characteristic of the secondaries generated under similar conditions. For each hadronic code we generate sets of $`10^5`$ collisions in order to analyze the secondaries produced by SIBYLL and QGSJET in $`\overline{p}p`$ and $`\overline{p}A`$ ($`A`$ represents a nucleus target of mass number $`A=10`$) at different projectile energies.

In all the considered cases we found that the number of secondaries coming from QGSJET collisions is larger than the ones corresponding to the SIBYLL case. This shows up clearly in the energy versus total number of secondaries two-dimension distributions. Although differences are not quite obvious at 100 TeV (Figs. 3 and 4) they grow up dramatically when the projectile energy reaches $`10^{20}`$ eV (Figs. 5 and 6). With the exception of protons and neutrons, the normalized energy distributions possess similar shapes for both models.

Again we remark that the case of $`p`$-nuclei collisions do not show major differences with respect to $`p`$-$`\overline{p}`$ case.

## III Air shower simulations

Proton induced air showers are generated using AIRES+SIBYLL and AIRES+QGSJET. Primary energies range from $`10^{14}`$ eV up to $`10^{20.5}`$ eV. To put into evidence as much as possible the differences between the intrinsic mechanism of SIBYLL and QGSJET we have always used the same cross sections for hadronic collisions, namely, the AIRES cross section plotted in Figs. 1 and 2.

All hadronic collisions with projectile energies below 200 GeV are processed with the Hillas Splitting algorithm , and the external collision package is invoked for all those collisions with energies above the mentioned threshold. It is worthwhile mentioning that for ultra-high energy primaries, the low energy collisions represent a little fraction (no more than 10% at $`10^{20.5}`$ eV) of the total number of inelastic hadronic processes that take place during the shower development. It is also important to stress that the dependence of the shower observables on the hadronic model is primarily related to the first interactions which in all the cases are ultra high energy processes involving only the external hadronic models. All shower particles with energies above the following thresholds were tracked: 500 keV for gammas, 700 keV for electrons and positrons, 1 MeV for muons, 1.5 MeV for mesons and 80 MeV for nucleons and nuclei. The particles were injected at the top of the atmosphere (100 km.a.s.l) and the ground level was located at sea level.

We have analyzed in detail the longitudinal development of the showers. The number and energy of different kind of particles have been recorded as a function of the vertical depth for a number of different observing levels (more than 100).

The charged multiplicity, essentially electrons and positrons, is used to determine the number of particles and the location of the shower maximum by means of four-parameter fits to the Gaisser-Hillas function .

In Fig. 7 $`X_{\mathrm{max}}`$ is plotted versus the logarithm of the primary energy for both the SIBYLL and QGSJET cases. It shows up clearly that SIBYLL showers present higher values for the depth of the maximum, and that the differences between the SIBYLL and QGSJET cases increase with the primary energy. This is consistent with the fact that SIBYLL produces less secondaries than QGSJET –as discussed in the previous section– and as a result, there is a delay in the electromagnetic shower development which is strongly correlated with $`\pi ^0`$ decays. The fluctuations, represented by the error bars, decrease monotonously as long as the energy increases, passing roughly from 95 g/cm<sup>2</sup> at $`E=10^{14}`$ eV to 70 g/cm<sup>2</sup> at $`E=10^{20.5}`$ eV.  At this stage it must be stressed that the mean values we have obtained for depths of shower maximum are slightly different to those recently presented by Pryke and Voyvodick . The main difference arise from the fact that in our treatment the mean free paths for hadron-nucleus collisions are the same in both models. Besides, one cannot say that interactions lengths are quite similar in these two models at $`10^{19}`$ eV (see Figs. 1 and 2), unless both are seen in the framework of QGSJET (see discussion in Sec. II).

As a representative case we are going to consider in more detail the behaviour of $`10^{20.5}`$ eV proton showers. In Fig. 8 the total number of pions, muons, gammas, and charged particles are plotted versus the vertical atmospheric depth. The observed behaviour of the electromagnetic shower is consistent with the discussion of the previous paragraphs: A shift in the maximum of the shower. It is worth to mention that even if the longitudinal development shows important differences at first stages, they decrease monotonously as far as the shower evolves. The smallest difference between models corresponds to the case of pions whose number at the ground level is very similar in both SIBYLL and QGSJET cases. On the other hand, the number of ground muons does present significant (even if not critical) differences.

With the particle data recorded we have evaluated lateral and energy distributions not only at ground altitude but also at predetermined observing levels. To the best of our knowledge this is the first time that the evolution of lateral and energy distributions along the longitudinal shower path is studied in such detail. In this paper we present the distributions corresponding to a subset of all the levels considered, taking into account particles whose distances to the shower axis are larger than 50 m.

The high-altitude lateral distributions (Figs. 9, 10 and 11) show important differences between SIBYLL and QGSJET; such differences diminish as long as the shower front gets closer to the ground level. The behaviour can be explained taking into account the differences between the number of SIBYLL and QGSJET secondaries reported in the previous section. Due to the fact that SIBYLL produces less number of secondaries, they have –in average– more energy and therefore the number of generations of particles undergoing hadronic collisions is increased with respect to the QGSJET case. As a result, during the shower development SIBYLL is called more times than QGSJET, and this generates a compensation that tends to reduce the difference in the final number of hadronic secondaries produced during the entire shower, and consequently in the final decay products, that is, electrons, gammas and muons.

The lateral distributions of electromagnetic particles are remarkably similar at both $`X_{\mathrm{max}}`$ and ground level.We want to stress that $`X_{\mathrm{max}}`$ is different in each model. However, it comes out from a more detailed analysis of the ground distributions that they are not strictly coincident and that the ratio between SIBYLL and QGSJET predictions does depend on $`r`$, the distance from the core. In fact, for electrons, this ratio runs from 1.25 for small $`r`$ to 0.73 for $`r1000`$ m, being equal to 1 at $`r350`$ m. A similar behaviour is observed for gammas where the lateral distributions intersect at $`r1000`$ m.

In the case of lateral muon distributions, QGSJET predicts a higher density for all distances, but the SIBYLL/QGSJET ratio is not constant, ranging from 0.74 near the core to 0.56 at 1000 m.

The energy distributions at varying altitudes (Figs. 12, 13 and 14) permit following the dynamics of the shower in great detail. The plots corresponding to 200 g/cm<sup>2</sup> clearly indicate that at such altitude there is an important fraction of high energy particles (more than 1 TeV). In particular, this is more evident for electrons and gammas (Figs. 12 and 14, compare the 200 g/cm<sup>2</sup> with the respective sea level distributions).

Finally, analyzing in detail the energy distributions of muons at ground level, we observe that the ratio of $`dN/d\mathrm{log}E`$ between SIBYLL and QGSJET is not constant: At the low (high) end of the spectrum it takes the value 1.0 (0.8) reaching a minimum of 0.54 around 250 GeV.

## IV Conclusions

Addressing the theoretical issues surrounding high energy hadronic collisions is intrinsically complicated since many variables are involved. However, this is crucial in understanding the data being recorded by present extremely high energy CR experiments (like the Akeno Giant Air Shower Array ) as well as CR next generation experiments (the future Pierre Auger Observatory  –fluorescence detector plus ground array– and the “eyes” of the OWL  that will deeply watch into the CR-sky).

In this work we have studied the sensitivity of parameters of hadronic interactions models (fitted to low energy data) when the algorithms are extrapolated several order of magnitudes. Perhaps the most oustanding difference between SIBYLL and QGSJET, as we had expected from our theoretical analysis, is reflected in the predicted number of secondaries after single $`p`$-$`\overline{p}`$ and $`\overline{p}`$-nuclei collisions. Such a difference increase steeply with rising energy. Our investigation on air showers throws up various other points of interest. In particular, we have reported that the different number of secondaries predicted remains noticeable during the first stages of the shower development. It is, of course, immediately evident that this follows again as a direct consequence of the lower inelasticity implemented in the SIBYLL generator when compared with the one in QGSJET. The study of the evolution of lateral and energy distributions along the longitudinal shower path allows us to clearly observe how the differences in the distributions become monotonously damped, yielding rather similar shapes when reaching the ground. Further, we have shown that the differences observed at ground level do depend on the distance to the shower core. Consequently, we are convinced that it will be possible to obtain relevant information about the hadronic interactions in air showers from the measurement of particle densities at distances far from as well as close to the shower core. This can be achieved if CR experiments are designed with appropiate dynamic ranges.

On the other hand, in our opinion most of the model discrepancies discussed in Sec. II will be naturally reduced with the help of data obtained from future accelerator experiments like the well known Large Hadron Collider (LHC).

###### Acknowledgements.

It is a pleasure to thank Tom Gaisser for valuable comments on technical aspects of SIBYLL. We are also indebted to Alberto Etchegoyen for granting us computing facilities. This work has been partially supported by CONICET and the FOMEC program.     

## V Legends

Fig. 1: In the figure we have plotted the corresponding $`p`$-air cross sections of SIBYLL (dashed line), QGSJET (dots), and AIRES (solid line). We have also superimposed experimental data obtained from collider experiments $``$, $`\mathrm{}`$, together with the ones obtained in cosmic ray experiments $``$, $``$, $`\mathrm{}`$.

Fig. 2: The corresponding $`\pi `$-air cross sections of SIBYLL, QGSJET and AIRES. The conventions adopted are the ones of Fig. 1.

Fig. 3: The figure displays the two–dimension distributions (energy vs number of secondaries) obtained from $`p`$-$`\overline{p}`$ scatterings (incident energy 100 TeV). In the left hand side we present the results of QGSJET while the right hand side corresponds to the ones of SIBYLL.

Fig. 4: $`p`$-nuclei ($`A`$=10) scatterings with incident energies of 100 TeV. We have used the conventions of Fig. 3.

Fig. 5: $`p`$-$`\overline{p}`$ scatterings with incident energies of 100 EeV. The conventions are the same as Fig. 3.

Fig. 6: $`p`$-nuclei ($`A`$=10) scatterings with incident energies of 100 EeV. We have used the conventions of Fig. 3.

Fig. 7: Simulation results for the average slant depth of maximum, $`X_{\mathrm{max}}`$, for proton induced showers, plotted versus the logarithm of the primary energy. The error bars indicate the standard fluctuations (the RMS fluctuations of the means are always smaller than the symbols). The squares (circles) correspond to SIBYLL (QGSJET).

Fig. 8: Longitudinal development of $`10^{20.5}`$ vertical proton showers. The average numbers of particles are plotted versus the atmospheric depth. The solid (dashed) line stands for the QGSJET (SIBYLL) case.

Fig. 9: Comparison between the recorded electron lateral distributions displayed by SIBYLL (grey) and QGSJET (black) at different atmospheric altitudes including the depth where the shower developes its maximum and the predictions at the ground level.

Fig. 10: Same as Fig. 9 for the case of muons.

Fig. 11: Same as Fig. 9 for the case of gammas.

Fig. 12: Electron energy distributions obtained with SIBYLL (grey) and QGSJET (black) at different atmospheric altitudes, see level and depth of shower maximum.

Fig. 13: Same as Fig. 12 for the case of muons.

Fig. 14: Same as Fig. 12 for the case of gammas.