Air Shower Simulation and Hadronic Interactions

The aim of this report of the Working Group on Hadronic Interactions and Air Shower Simulation is to give an overview of the status of the field, emphasizing open questions and a comparison of relevant results of the different experiments. It is shown that an approximate overall understanding of extensive air showers and the corresponding hadronic interactions has been reached. The simulations provide a qualitative description of the bulk of the air shower observables. Discrepancies are however found when the correlation between measurements of the longitudinal shower profile are compared to that of the lateral particle distributions at ground. The report concludes with a list of important problems that should be addressed to make progress in understanding hadronic interactions and, hence, improve the reliability of air shower simulations.


Introduction
After the first interaction of a cosmic ray particle of very high energy in the atmosphere a multitude of subsequent interactions, leading to particle multiplication, and decay processes give rise to a cascade of secondary particles called extensive air shower (EAS) [1]. With the electromagnetic and weak interactions being well described by perturbative calculations within the Standard Model of Particle Physics, the limited understanding of strong interactions becomes the dominant source of uncertainties of shower predictions. Even though Quantum Chromodynamics (QCD) is the well-established and experimentally confirmed theory of strong interactions, only processes with very large momentum transfer can be predicted from first principles until now. It is still not possible to calculate the bulk properties of multiparticle production as needed for air shower simulation. Additional, simplifying assumptions as well as phenomenological and empirical parametrizations are needed to develop models for hadronic interactions describing various particle production processes. These additional assumptions need to be verified, parametrizations constrained, and parameters tuned by comparisons to accelerator data.
The simulation of extensive air showers forms one of the pillars on which the data analysis of modern experiments for ultra-high energy cosmic rays (UHECRs) rests. Improving the understanding and modeling of hadronic particle production is one of the important prerequisites for a reliable interpreta-tion of UHECR data. Even though calorimetric techniques based on the measurement of fluorescence and Cherenkov light have been developed for a nearly model-independent energy determination of extensive air showers, determining the type of the primary particles and, in particular, estimating the mass of the primary particle can only be done with the help of sets of simulated reference showers. 1 Much effort has been devoted to improving both the simulation of extensive air showers in the atmosphere and the understanding of the corresponding, and typically very complex detector response function. In the following we will concentrate on discussing the simulation of extensive air showers as this aspect of the overall data simulation chain is the same in all UHECR experiments. We will assume that uncertainties arising from simulating the detector response due to the shower particles are much smaller and can be disregarded here. This is certainly not guaranteed per se and it is the task of each collaboration to verify the quality of the detector simulation before attributing possible discrepancies to, for example, shortcomings in the understanding of air showers or hadronic multiparticle production.
After recalling some basic features of air showers (Sec. 2) we will give an overview of the most frequently used code packages for air shower simulation in Sec. 3. The LHC data allow us to test the hadronic interaction models employed in these code packages, for the first time, at equivalent energies beyond that of the knee in the cosmic ray spectrum. In Sec. 4 some representative model predictions are compared to LHC data and implications are discussed. A good overall bracketing of the LHC data by model predictions is found, even though each of the models will have to be improved to obtain a satisfactory description of the LHC data. Therefore it is not surprising that the distributions of most of the shower observables are well reproduced by shower simulations (Sec. 5).
Taking advantage of the hybrid detection setups of the latest generation of UHECR detectors, longitudinal and lateral particle distributions can be measured shower-by-shower. A comparison of the particle densities at ground with the ones expected according to the shower longitudinal profile reveals possible shortcomings of the shower predictions that are not yet understood. The overall data set of such comparisons is discussed in Sec. 6 and the results from the Auger, TA, and Yakutsk Collaborations are compared. There are strong indications that a significant part of the shortcomings of the shower predictions are related to the muonic shower component. Therefore we summarize recent developments aiming at better understanding muon production in air showers in Sec. 7.
Given the increased awareness of apparent deficits of current air shower predictions and related theory developments, the wealth of new data from the LHC, and the high-quality measurements of the current generation of air shower observatories, significant progress can be expected in the reliability of air shower predictions in the next few years. Open questions that need to be addressed to optimally benefit from these developments to improve the interpretation of EAS data and to solve some outstanding problems will be listed in the concluding section of this report.

Physics of Air Showers
Reviews of the physics of extensive air showers and hadronic interactions can be found in, for example, [3][4][5][6] and a very instructive extension of Heitler's cascade model to hadronic showers is given in [7]. Here we give only a qualitative introduction to some aspects of shower physics that will be needed in the discussions later on.
With pions being the most abundant secondary particles of hadronic interactions at high energy, the difference in lifetime and decay products of neutral and charged pions lead to fundamentally different ranges of interaction energies that give rise to the electromagnetic and muonic shower components. In air showers, neutral pions decay almost always before interacting again (except E π 0 > 10 19 eV) and provide high-energy photons that feed the electromagnetic shower component. In contrast, the energy of charged pions has to be degraded down to some 30 − 100 GeV before the probability to decay will exceed that of another interaction. In first approximation, the depth profile of charged shower particles, being dominated by e + e − of the electromagnetic component, is linked to the secondary particles of the first few interactions in the hadronic core of the shower. In contrast, the muons at ground stem from Shown are also the sub-showers produced by the secondary particles of the first 100 highest-energy interactions and the cumulative profile of these sub-showers [8].
The simulation was done with a modified version of CONEX [9] for a proton shower of 10 19 eV. a chain of 8 − 10 successive hadronic interactions with the energy of the last interaction producing pions or kaons that lead to observed muons being in the range of 20 − 200 GeV [10,11]. This is illustrated in Fig. 1 showing the sub-showers produced by the secondary particles of the first 100 highest-energy interactions in an air shower [8]. While more than 50% of the electromagnetic shower particles are produced by decaying neutral pions of the first 100 interactions in an air shower, only a negligible fraction of the muonic component stems from decaying charged pions and kaons of the same interactions.
The longitudinal profile of showers as well as the energy and angular distributions of charged particles in a shower exhibit universality features [12][13][14][15][16][17] if considered as a function of shower age, here approximated by with X being the slant depth along the shower axis and X max the depth of shower maximum. Universality features are routinely used to estimate the Cherenkov light signal of showers, see, for example, [15,18]. Not only the electromagnetic particles but also the production depth as well as the energy and transverse momentum distributions of muons can be described by functions [19,20] that depend only weakly on the assumed hadronic interaction model. Breaking down the overall shower signal at ground into different muonic and electromagnetic components, the description of showers can be improved considerably [21][22][23].
While the bulk of charged particles is considered for longitudinal profiles of the electromagnetic and muonic shower components, lateral particle distributions at ground are only measured at large distances from the shower axis. Particle densities at, for example, 600 or 1000 m from the shower core are only indirectly related to the overall particle multiplicity of a shower at ground. This is reflected in the fact that, for example, the last interaction for producing a muon observable at ground has a median energy of ∼ 100 GeV for the KASCADE array (typical primary energy 3 × 10 15 eV; lateral distance 40 − 200 m, see [11]) and only 20 − 30 GeV for the Auger array (typical primary energy 10 19 eV; lateral distance 1000 m, see [24]). Still universality arguments can also be applied to particle densities far away from the core [21] but the model-related differences between the predictions are larger. In particular, low-energy interactions can lead to a model-dependent violation of universality profiles that is reflected in both the muonic and electromagnetic shower content at large lateral distance [16,23].

Shower Simulation Packages and Hadronic Interaction Models
The most frequently used code packages for the simulation of air showers for the Auger, TA and Yakutsk experiments are AIRES [25], CORSIKA [26], and COSMOS [27]. In these packages, the Monte Carlo method is applied to simulate the evolution of air showers with all secondary particles above a user-specified energy threshold. The large number of secondary particles that would have to be followed in the simulation of ultra-high energy showers makes the simulations so time consuming that thinning techniques are applied. Only a representative set of particles is included in the simulation below a pre-defined energy threshold. The discarded particles are accounted for by increasing the statistical weight of the particles remaining in the simulation [28,29]. The main drawbacks of this method are artificial fluctuations related to particles with large statistical weight [30,31] and the need for treating weighted particles in the detector simulations. The Auger and TA Collaborations have developed different methods for de-thinning simulated showers, see [32] and [33], respectively. Also a limited number of fully simulated showers of the highest energies are available in each collaboration (see, for example, [34]) but the statistics is very limited and there is a lack of methods to work with such large amounts of data in an efficient way.
The classic Monte Carlo codes are complemented by hybrid simulation programs that combine the Monte Carlo technique with either numerical solutions of cascade equations (CONEX [9] and SENECA [35]) or libraries of pre-simulated air showers (COSMOS [27]). The advantage of the hybrid codes is the much reduced CPU time needed for simulating showers at very high energy and the possibility to simulate showers either without thinning or with very low statistical weights. But there is always a remaining risk that rare shower fluctuations or local correlations might not be correctly described.
Different codes are applied for the simulation of the electromagnetic shower component. The code used in AIRES is based on the one developed by Hillas for MOCCA [28,36]. CONEX, CORSIKA, and SENECA are interfaced to modified versions of EGS4 [37]. Similar to AIRES, COSMOS comes with a custom-developed code for electromagnetic interactions called EPICS [27]. The suppression of high-energy interactions of photons and electrons due to the Landau-Pomeranchuk-Migdal (LPM) effect (for example, see [38]) is taken into account in all code packages.
Several comparisons of the predictions of shower simulation packages for the same interaction models are available in the literature, see [4] and [61,62,59,60,63]. In the ideal case, the predictions for air showers of a given energy and primary particle depend only on the user-selected energy thresholds and on the choice of hadronic interaction models selected for the simulation. While it is possible to chose a high-energy interaction model that is supported in all shower simulation codes, less flexibility is available for the low-energy interaction models, which are often different in each of the simulation packages. Therefore the results of these comparisons typically show a good agreement of the longitudinal shower profiles, which are mainly related to hadronic interactions at high energy. Differences are found for observables that are sensitive to low-energy interactions, which are, for example, the total number of muons and particle densities at large lateral distances. The differences between the predictions of different simulation packages are typically of the order of 5%, increasing in some phase space regions up to ∼ 10%. For example, longitudinal shower profiles from a recent comparison between CORSIKA and COSMOS [59,60] are shown in Fig.2.

Performance of Hadronic Interaction Models
Most of the interaction models used in air shower simulations are not commonly applied in high energy physics (HEP) simulations and, conversely, HEP models are not used for air shower simulations. This is related to the fact that HEP models are typically limited to a set of primary particles that are available in accelerator experiments and optimized only for collider energies. Similarly, cosmic ray interaction models are not routinely used by HEP collaborations for acceptance correction calculations including  Fig. 3. Pseudorapidity distribution of charged particles measured at LHC. CMS data [64] are shown together with model predictions. Similar data sets are available from the ATLAS and ALICE Collaborations [65][66][67] (not shown). Below 2 TeV center-of-mass energy the models were tuned to previously available collider data from Tevatron and SPS. the comparison of raw, uncorrected data as this can be done with the HEP event generators as well. Furthermore dedicated HEP event generators describe high-p ⊥ and electroweak physics processes in much more detail and contain many more parameters for improving the description of particular distributions than any of the cosmic ray or general purpose models. In contrast, interaction models for cosmic ray physics are tuned to describe data of accelerator measurements over a wide range of energies, often sacrificing a perfect description of some distributions in favour of a good overall reproduction of the whole data set. Hence comparisons of model predictions to data from fixed-target and collider experiments are typically made by the authors of the models and using only fully acceptance corrected data. However, over the last years much progress has been made regarding the interaction between the cosmic ray and high energy physics communities. A number of workshops and meetings (including this one) and the direct involvement of cosmic ray physicists in HEP collaborations have lead to a much higher level of awareness and intensified communication of both communities.
Comparisons of high-energy interaction models used for air shower simulations with LHC data can be found in, for example, [69][70][71]. Here only some representative examples can be included for illustration. The pseudorapidity distribution of charged particles is shown together with model predictions in Fig. 3. Not only most central particle distributions are reasonably well bracketed by the model predictions, also the predicted energy flow at larger pseudorapidity is in good agreement with the LHC data. This can be seen in Fig. 4 where the energy flow measured by CMS [68] is compared to model predictions. Note that all these model results are true predictions as the models were tuned years before LHC data became available. Nevertheless, in many cases, the minimum bias data are better described by interaction models developed primarily for cosmic ray physics than the various tunes of standard HEP models [68,69].
On the other hand, there are also some important distributions measured at LHC that are not well described or bracketed by these pre-LHC interaction models. Examples of large deviations from measurements are multiplicity distributions of charged particles and particle production spectra at large pseudorapidity and Feynman-x. The possible impact of these deviations on air shower predictions is subject to ongoing research and not yet fully understood. In Fig. 5   photons produced in forward direction is shown [72]. The data have been obtained within the LHCf experiment that is specifically built for measuring particles in the very forward direction.
In general, the predictions of the different cosmic ray interaction models bracket the LHC data reasonably well [69]. This observation supports the hope that simulating air showers with different interaction models should also bracket the correct, but unknown predictions for air showers. On the other hand, this success does not mean that there is no need to improve the interaction models. There is not a single interaction model that reproduces most of new LHC data well and, hence, would be much better than the others. Each of the models needs to be re-tuned or underlying ideas be extended to obtain a better description of the data.
Many new data sets measured in LHC experiments and also the fixed-target experiments NA49 [73] and NA61 [74,75] became available within the last 2−3 years and the process of tuning and modifying the interaction models to obtain an improved description is ongoing. At the time of writing these proceedings the new model versions QGSJET II.04 [70] and EPOS LHC [71] have been released as the first versions of interaction models tuned to LHC data. In addition, work is in progress to develop improved versions of SIBYLL and DPMJET that are also tuned to the new LHC data.
Finally it should be mentioned that, in addition to accelerator data, also air shower measurements provide important information on hadronic interactions. For example, cross section measurements made with air shower observables [76] can be used to estimate not only the rise of the cross section but also to verify the model approximations used to calculate the proton-air cross section from the data on proton-proton cross sections. Even though the systematic uncertainties of air shower measurements of this kind are much higher than that of accelerator data, such measurements are the only way to study particle production at interaction energies well beyond the range of colliders [77][78][79][80].

Overall Description of Shower Characteristics
Simulating air showers with a realistic energy spectrum, primary mass composition, and arrival direction distribution, the predicted and observed distributions of the various observables of reconstructed showers can be compared. Such comparisons are very important as they are end-to-end tests of the simulation chain for the experiments. Only if good agreement is found one can cross-check the impact of quality cuts applied in the process of data analysis and avoid unexpected biases. The results of such end-to-end simulations depend, of course, on the assumed mass composition of the primary particles and the employed hadronic interaction models.
Examples of such comparisons for fluorescence detectors can be found in [81][82][83][84]. A number of unpublished comparisons of typical observables measured with the Auger and TA surface detector arrays (zenith and azimuth angles, number of stations per event, signal size distribution) have been presented at this meeting. They all show very good agreement between the measured distributions and the corresponding Monte Carlo predictions. While the TA simulations were done for proton primaries, the Auger simulations were based on a 50/50 proton-iron mixture. Both collaborations used QGSJET II.03 as high-energy interaction model. The TA data are better described by a primary composition of only protons than using 100% iron. In case of the Auger simulations, which were also made assuming only proton or iron primaries, it was found that most of the surface detector observables exhibit only a very limited sensitivity to the primary mass composition.
It should be considered as an important success of hadronic interaction models and modern air shower simulation packages that such a good overall description of the general features of the observed showers is reached. However, different combinations of primary mass composition and shower energy can lead to very similar surface detector signals (for example, the signal at 600 or 1000 m in units of that of vertical muons). Therefore it is very important to use additional information to determine, for example, the primary energy in a composition-independent way. This is done with fluorescence telescopes in the case of the Auger Observatory [85] and Telescope Array [86], and with non-imaging Cherenkov light measurements in the Yakutsk setup [87].

Detailed Comparison with Shower Data
Detecting air showers simultaneously with fluorescence telescopes and an array of surface detectors at ground is often referred to as hybrid measurements. One of the first setups of this type was the Conference Title, to be filled prototype instrument of the High Resolution Fly's Eye (HiRes) [88] that was operated in coincidence with the CASA-MIA array [89] in Utah. Already in this prototype experiment strong indications for a discrepancy between the measured density of muons at 600 m from the shower core and the one expected from simulations were found, if the composition inferred from the measurements of the longitudinal shower profiles was used [90]. While the data on the depth of shower maximum suggested the conclusion that a transition from a heavy to a light, almost proton dominated composition is seen, the density of muons stayed above or at the level of the predictions for iron primaries. Covering energies higher than those of the early HiRes-MIA measurements, the data of the Auger Collaboration also indicate a discrepancy between the observed and predicted number of muons in air showers [91]. Different methods have been applied to obtain either the expected surface detector signal S (1000) or the muon density at 1000 m from the shower core for a given primary energy or even an observed longitudinal shower profile (see [92] for a review given at this meeting). This discrepancy is most directly displayed in a shower-by-shower comparison of simulated profiles with those obtained from hybrid measurements, see Fig. 6. Compared to proton showers as reference, the measured values of S (1000) are about 1.5 to 2 times higher that the simulated ones. Studying inclined showers of about 10 19 eV with zenith angles larger than 60 • , whose particle content at ground is completely dominated by muons, gives a ratio of [91] N µ,data N µ,MC QGSJET,p = 2.13 ± 0.04(stat) ± 0.11(sys), where the systematic uncertainty due to the 22% energy scale uncertainty is not included. This ratio is about 1.2 for iron showers simulated with EPOS 1.99. The data on inclined showers imply that the observed discrepancy is closely related to, if not dominated by a deficiency in simulating the muon component of air showers. The detectors of the Auger, TA and Yakutsk air shower arrays are differing in their sensitivities to muons and electromagnetic particles, and the Yakutsk array even contains a number of shielded muon detectors. For example, the contribution of muons to the scintillator signal of the TA surface detectors will never exceed 30% for showers with zenith angles less than 45 • . In comparison, a signal fraction of 30−80% is expected due to muons for the Auger water Cherenkov tanks in the zenith angle range from 0 − 60 • . Therefore similar studies of the TA and Yakutsk Collaborations can provide a very important, independent information on a possible muon deficit of air shower simulations.
The Yakutsk Collaboration compared the predicted and measured charged particle and muon lateral distributions for a number of high energy showers. Within the uncertainties implied by the unknown primary composition, no significant discrepancies between the predictions obtained with both QGSJET II/FLUKA and EPOS/UrQMD and the data are found [93].
The TA Collaboration studied the difference between the energy one would assign to a shower based on the calorimetric fluorescence light measurement and that derived from the comparison with simulations of the surface detector signal. It was found that the fluorescence-assigned energy is about 27% lower than the energy estimated from the comparison with simulated showers [94]. Given the limited sensitivity of the TA surface detector stations, no direct study of the muon contribution to the overall signal has been done so far.
The observations of all three collaborations can be brought into qualitative agreement if the different energy scales of the experiments are taken into account. The factor with which the energy scales need to be re-scaled are listed in the report of the Spectrum Working Group [95]. For example, for the same shower, the Yakutsk simulations are done at an energy about two times higher than that of the Auger Observatory. Hence, also about a factor two more muons are predicted in the Yakutsk simulations in comparison to Auger simulations. Given that the Yakutsk simulations agree with their measurements, this implies two times more muons in the Yakutsk data than predicted by simulations if the Auger energy scale is used. A similar estimate can be done for the TA-Auger comparison, where one has to keep in mind that only ∼ 20% of the scintillator signal stems from muons. If the discrepancy observed by the Auger Collaboration would be attributed to the muon component only, the expected energy rescaling would be somewhat smaller than the ∼ 27% reported by the TA Collaboration [94].
It remains to be shown how much of the apparent discrepancy between the surface detector signals and the predictions based on calorimetric shower energy determination can be cured by adopting different energy scales in the experiments. A reduction of the systematic uncertainties of the energy assignment to air showers will be needed to quantify more reliably possible discrepancies between shower data and predictions.

Hadronic Interactions and Muon Production in Air Showers
The observation of a possible muon deficit in simulated showers in comparison to Auger measurements, first reported already in 2007 [96], has triggered a number of theoretical studies searching for modifications of hadronic interactions that could result in an enhancement of the muonic shower component. Changes of the inelastic cross section, inelasticity of interactions, secondary particle multiplicity at high energy, and many other parameters lead only to moderate changes of the predicted number of muons (see, for example, [97,98]). Also scenarios with new physics processes such as string percolation [99] or the drastic change of interaction properties due to, for example, chiral symmetry restoration [100] have been discussed.
So far all proposed changes that lead to a significant increase of the number of muons in air showers are directly or indirectly based on either one or both of the following effects: (i) An increase of the production rate of particles that do not decay, for example baryon-antibaryon pairs, leads to higher muon multiplicities of showers since these particles will stay being part of the hadronic shower component and loose their energy only by producing further hadronic particles. This effect has been discussed first in [101] and is one of the differences of the EPOS model with respect to the other models [102].
(ii) A change of the type of the leading particle produced in inelastic interactions can also be very efficient in reducing the energy that is transferred to the electromagnetic shower component [103]. The majority of sub-showers in an extensive air shower are initiated by charged pions. The chance probability of producing a charged pion or a neutral pion as leading particle in charged pion interactions is about 2 : 1. Replacing all leading π 0 by ρ 0 mesons, which have the same valence quarks but are spin 1 particles, leads to a drastic enhancement of muon production since neutral ρ mesons decay immediately into two charged pions [104]. Indeed, fixed-target data [105][106][107] indicate that, in contrast to conventional model predictions, the production of ρ 0 dominates that of π 0 for Feynman-x larger than 0.5.
Both effects are the reason for the higher muon multiplicities predicted with EPOS. Leading rho mesons are explicitly generated for pion interactions in the new version II.04 of QGSJET, boosting the muon multiplicity in air showers significantly. Depending on the relative importance of effect (i) and energy distribution of the baryon pairs, the energy spectrum of the additionally produced muons is very soft and might be in contradiction to shower attenuation data.
Finally it should be mentioned that also an increase of the production of kaons could lead to a higher muon multiplicity and also a harder muon energy spectrum [103,99].

Conclusions and Outlook
Over the last two decades the quality and predictive power of air shower simulations has improved significantly and a very good overall description of most of the shower features observed in experiments has been reached.
The predictions of independently developed shower simulation packages agree reasonably well with each other if the same hadronic interaction models are used for the simulations. While there is very good agreement for the electromagnetic shower component, the differences between the predictions for muon multiplicities and lateral distributions can be as large as ∼ 5 − 10%. These differences are most likely related to the use of different low-energy interaction models. A comprehensive comparison of the different shower simulation packages, similar to the recent study of CORSIKA and COSMOS predictions [60], should be made to quantify the systematic uncertainties of the predictions and possibly also to identify and address shortcomings of the simulation packages.
The limited theoretical understanding of modeling hadronic multiparticle production together with the limitations of accelerator measurements in energy, covered phase space, and projectile-target combinations are the dominating source of systematic uncertainties of air shower predictions. Because of this, the systematic uncertainties of the model predictions cannot be estimated reliably. More work is needed to improve QCD calculations in the low-p ⊥ domain, in particular to understand screening and saturation effects, and to develop alternative models to describe particle production.
The new LHC data provide extremely useful input for tuning hadronic interaction models. Even though the first LHC data on multiplicities and cross sections were well bracketed by the predictions of interaction models used for air shower simulations, the comparison to data revealed the need for further model developments and tuning. Improved and re-tuned versions of EPOS and QGSJET are already available and similar versions of DPMJET and SIBYLL are in preparation.
The interaction between the cosmic ray and high energy physics communities has intensified and the direct engagement of cosmic ray physicists in LHC and fixed-target experiments has lead to a much better understanding of the needs of the cosmic ray community. There is also very large interest from the side of LHC communities to have cosmic ray physicists being involved in the analysis of accelerator data. One example of the achieved progress is the use of cosmic ray interaction models by LHC collaborations to compare data with predictions in publications.
Both the Auger and TA Collaborations have found indications for a discrepancy between the expected and observed surface detector signals for showers with fluorescence energy measurement. Accounting for the different energy scales of the Auger, TA, and Yakutsk experiments the observed discrepancies are consistent with each other. Most likely, the dominating sources of the discrepancies are shortcomings of the simulation of muon production in air showers. The systematic uncertainties of the overall energy scales of the experiments of 20 − 30% will have to be reduced significantly to be able to combine data from different experiments for a stronger test of the shower predictions.
While the production of electromagnetic particles is clearly dominated by the electromagnetic cascade induced by photons of very high energy due to neutral pion decay early on in the shower evolution, the muonic component receives contributions from all high-energy interactions above ∼ 20 GeV lab. energy. Different modifications of the simulation of particle production in low-and high-energy interactions have been considered to increase the number of muons in air showers. Progress has been made by understanding that both production of stable or long-lived hadrons such as baryon pairs and modifications to the particle types generated as leading particles are efficient mechanisms to increase the muon multiplicity without changing significantly the longitudinal shower profile. It remains to be shown whether the discrepancies between the observed and predicted surface detector signals for a given shower energy can be resolved by improving low-and high-energy interaction models using conventional physics assumptions or whether these observations indeed indicate a change of the characteristics of hadronic interactions at energies beyond that of LHC.