Open issues in hadronic interactions for air showers

. In detailed air shower simulations, the uncertainty in the prediction of shower observables for different primary particles and energies is currently dominated by differences between hadronic interaction models. With the results of the ﬁrst run of the LHC, the difference between post-LHC model predictions has been reduced to the same level as experimental uncertainties of cosmic ray experiments. At the same time new types of air shower observables, like the muon production depth, have been measured, adding new constraints on hadronic models. Currently no model is able to consistently reproduce all mass composition measurements possible within the Pierre Auger Observatory for instance. Comparing the different models, and with LHC and cosmic ray data, we will show that the remaining open issues in hadronic interactions in air shower development are now in the pion-air interactions and in nuclear effects.


Introduction
Knowing the elemental composition of cosmic ray particles arriving at Earth is of crucial importance to understand the production and propagation of cosmic rays. Unfortunately, cosmic rays can be measured only indirectly above an energy of 10 14 eV, through the cascades of secondary particles, called extensive air showers (EAS), that they produce in the atmosphere (for a recent review, see [1]). Only by simulating the generation of EAS and comparing the predictions with measurements can one draw conclusions on the primary mass composition of the arriving particles [2]. With the operation of modern largescale experiments, the reliability of air shower simulations has become the source of the largest systematic uncertainty in the interpretation of cosmic-ray data [3][4][5][6][7][8][9]. While the electroweak interaction processes are reasonably well understood, modeling of hadronic multi-particle production is subject to large theoretical uncertainties that are, moreover, difficult to estimate [10][11][12].
The Large Hadron Collider (LHC) at the CERN laboratory allows us to access, for the first time, the energy region above the cosmic ray spectral knee in the laboratory with about 10 17 eV in the laboratory frame. Therefore an analysis of inclusive particle data taken at the LHC is particularly interesting for constraining existing hadronic interaction models and for testing possible new mechanisms of hadron production [13]. The first published data from LHC experiments have mostly been taken with detectors covering the central phase space region in pseudorapidity (|η| < 2.5). This region is most easily accessible in collider experiments and is also the region of the highest rapidity-density of produced particles. The first data have been compared to cosmic ray models in [14]. On the other hand, since the number of particles in an air shower is roughly proportional to the energy of the primary particle, the most energetic outgoing particles of an interaction, emitted in the very forward region of a a e-mail: tanguy.pierog@kit.edu collider experiment -such as in diffractive interactionsare the most important ones for understanding air showers. For the first time at the LHC, collider experiments include a large variety of forward detectors to study forward particle and energy spectra which could have a direct impact on air shower development [15]. These latest measurements are not yet taken into account in the available hadronic interactions models, but are very important to understand the open issues in these models and for their future development.
At the same time, a new generation of hybrid cosmic ray detectors such as the Pierre Auger Observatory [7] (surface and fluorescence detectors), the IceCube/IceTop experiments [16,17] (low energy muons at the surface and high energy muons deep underground) or the KASCADE/KASCADE Grande collaboration [18,19] (muons of different energies and at different distances) gives access to various precise measurements of the mean logarithmic mass of cosmic rays within the same experiment. By definition the mean logarithmic mass should be independent of the measurement technique. If the physics is well described by a given hadronic model, the masses obtained from different observables should be consistent. This constraint is much stronger than the traditional test limiting the results to the range between proton and iron induced showers. This is now satisfied in most of the cases, but none of the current models is able to give a fully consistent picture of the different observable within a given experiment [20][21][22].
In this paper, we will discuss the remaining open issues in the hadronic model predictions after LHC data and their consequences on air shower observables. In the first section, we will describe the latest hadronic interaction models and then compare their results for the observables important for the air shower development in the second section. The main source of remaining uncertainties will be identified. Using detailed Monte Carlo simulations done with CONEX [23], the new predictions for X max and for the number of muons will be presented. Finally we will take the example of the muon production depth (MPD) measured by the Pierre Auger Observatory [20] to see how air shower measurements can constrain hadronic interaction physics and can be used to solve the remaining open issues.

Hadronic interaction models
There are several hadronic interaction models commonly used to simulate air showers. Here we will focus on the three high energy models which were updated to take into account LHC data at 7 TeV: QGSJETII-03 [24,25] changed into QGSJETII-04 [26], EPOS 1.99 [27,28] replaced by EPOS LHC (V3400) [29], and Sibyll 2.1 [30][31][32] updated to Sibyll 2.3 [33] all available since CORSIKA V7.5600 [34]. There is no major change in these models but in addition to some technical improvements, some parameters were changed to reproduce TOTEM [35] cross sections. These are based on Gribov-Regge multiple scattering, perturbative QCD and string fragmentation.

EPOS model
EPOS LHC is a minimum bias monte-carlo hadronic generator used for both heavy ion interactions and cosmic ray air shower simulations. It is based on EPOS 1.99 retuned to reproduce LHC data on a higher precision level. As with most of high energy hadron-hadron interaction models, it is based on the simple parton model which can be seen as an exchange of a "parton ladder" between the two hadrons.
In EPOS, the term "parton ladder" is actually meant to contain two parts [36]: the hard one, as discussed above, and a soft one, which is a purely phenomenological object, parameterized in Regge pole fashion. This is the so called Pomeron used as elementary parton-parton interaction in EPOS.
In addition to the parton ladder, there is another source of particle production: the two off-shell remnants, see Fig. 1. We showed in Ref. [37] that this "three object picture" can solve the "multi-strange baryon problem" of conventional high energy models, see Ref. [38].
Hence EPOS is a consistent quantum mechanical multiple scattering approach based on partons and . c) The most simple and frequent collision configuration has two remnants and only one cut Pomeron represented by two q − q strings. d) One of the q string ends can be replaced by a qq string end. e) With the same probability, one of the q string ends can be replaced by a qq string end.
strings [36], where cross sections and the particle production are calculated consistently, taking into account energy conservation in both cases (unlike other models where energy conservation is not considered for cross section calculations [39]). The main consequence of this energy sharing process is that the number of Pomerons generated event-by-event do not follow a simple Poissonian distribution. As a consequence it is much less likely to produce events with a very large number of Pomerons (large multiplicity) compared to the standard Gribov-Regge approach like in QGSJETII. Nuclear effects related to Cronin transverse momentum broadening, parton saturation, and screening have been introduced into EPOS [27]. Furthermore, high density effects leading to collective behavior in heavy ion collisions are also taken into account [40].
Within a Monte Carlo, first the collision configuration is determined: i.e., the number of each type of Pomerons exchanged between the projectile and target is fixed and the initial energy is shared between the Pomerons and the two remnants. Then particle production is accounted from two kinds of sources, remnant decay and cut Pomeron. A Pomeron may be regarded as a two-layer (soft) parton ladder attached to projectile and target remnants through its two legs. Each leg is a color singlet, of type qq , qqq or qqq from the sea, and then each cut Pomeron is regarded as two strings, Cf. Fig. 2a)

and b).
It is a natural idea to take quarks and anti-quarks from the sea as string ends for soft Pomerons in EPOS, because an arbitrary number of Pomerons may be involved.
Thus, besides the three valence quarks, each remnant has additional quarks and anti-quarks to compensate the flavors of the string ends, as shown in Fig. 2-c. According to its number of quarks and anti-quarks, to the phase space, and to an excitation probability, a remnant decays into mesons and/or (anti)baryons [37]. Furthermore, this process leads to a baryon stopping phenomenon in which the baryon number can be transferred from the remnant to the string ends (for instance in 2-d, depending on the process, the 3q + 3q can be seen as 3 mesons or a baryonantibaryon pair).  In the case of a meson projectile, this kind of diquark pair production at the string ends leads to an increase of the (anti)baryon production in the forward production in agreement with low energy pion-nucleus data [41] as shown Fig. 3. Comparing to the QGSJETII model which does not have diquark as string ends or using only qq as string end in EPOS, we can clearly see that this process is needed to reproduce experimental data. As a consequence it contributes to the larger number of muons in air shower simulations with EPOS.
Energy momentum sharing and remnant treatment are the key points of the model concerning air shower simulations because they directly influence the multiplicity and the inelasticity of the model.

QGSJETII model
QGSJETII-04 [26,42] model is a minimum bias hadronic interaction model optimized for air shower simulations. It has a minimum set of parameters to reduce the uncertainty due to the extrapolation to high energy and as a consequence has a less detailed description of the final stage of hadronic interactions (no final state effect, no rare particle production, etc.) which limit the data sets to which it can be compared. It is the current last evolution of the series of models based on Quark-Gluon and Strings model with Jet [43].
The elementary interaction used in QGSJETII is a semi-hard Pomeron very similar to the one used in EPOS. But to take into account non-linear effects at high energy, a very complex resummation scheme has been developed to take into account any type of Pomeron-Pomeron interactions resulting in a net fan diagram [42] which can be used for both cross-section and particle production. But then only the standard Gribov-Regge type of calculation can be used, meaning that energy sharing is not taken into account in cross-section calculation and when the total number of multiple scattering is calculated.
Like in EPOS, QGSJETII has some remnant from the projectile and target but in a simplified scheme which does not allow more than one quark exchange with the string ends. In particular diquarks are not allowed as string ends which gives quite different forward baryon production compared to EPOS as shown Fig. 3.

Sibyll model
Like QGSJETII, the Sibyll model is a minimum bias hadronic interaction model optimized for air shower simulations but with a different approach. The crosssection calculation is based on the same Gribov-Regge calculation without energy sharing and diffraction is now done with a two state approach like in QGSJETII. However the eikonal is based on the minijet model with an energy dependent p t cutoff [44][45][46][47] and the string hadronization is using the Lund model [48] producing all type of particles.
The new Sibyll includes the same simplified remnant treatment as in QGSJETII. This allows for instance the increase of leading ρ 0 in π −p/A interactions. Compared to Sibyll 2.1 the new version has an improved production of baryon-anti-baryon pairs, in particular from the minijet (hard) particle production, and a phenomenological model for the production of charm particles which is important for the production of high energetic muons and neutrinos.
Unlike the two other models, Sibyll uses the semisuperposition model to treat the nuclear interactions meaning that for a nuclear projectile the result is the sum of a given number of nucleon-nucleus interactions without additional nuclear effects. Another main difference from the other models is a true scaling of forward particle production as discussed in [49].
So, even if the three models are based on very similar approaches (parton ladder and Gribov-Regge based multiple interactions), the detailed treatments of energy sharing, non-linear effects, nuclear effects and remnant production lead to different extrapolations in both proton and pion interactions and thus for air shower observables as shown in the next section.

Model comparison
A toy-model, as described in [50], only gives a very much over-simplified account of air shower physics. However, the model allows us to qualitatively understand the dependence of many air shower observables on the characteristics of hadronic particle production. Accordingly the parameters of hadron production which are most important for air shower development are the cross section (or mean free path), the multiplicity of secondary particles of high energy, the elasticity and the production ratio of neutral to charged particles. Until the start of the LHC, these parameters were not well constrained by particle production measurements at accelerators. As a consequence, depending on the assumptions of how to extrapolate existing accelerator data, the predictions of hadronic interaction models were very different [51]. We will show that the extrapolation to high energy is not really the issue anymore.

Cross section
As shown in [50], the cross section is very important for the development of air showers and in particular for the depth of shower maximum. As a consequence, the number  of electromagnetic particles at ground level is strongly correlated to this observable (if the shower maximum is closer to ground, the number of particles is higher).
The inelastic cross section of proton-proton scattering is usually used as an input to fix basic parameters in all hadronic interaction models. Therefore it is very well described by all the models up to the LHC, where data exist [52]. As shown on Fig. 4, thanks to the measurements at the LHC, even the extrapolations up to the highest energy are now very similar. In all the figures EPOS LHC is represented by a full (blue) line, QGSJETII-04 by a dashed (red) line and Sibyll 2.3 by a dashed-dotted (green) line.
However, using these models to predict proton-air and pion-air inelastic cross-sections as shown on Fig. 5, one can notice that significant differences appear which will have direct consequences on air shower development. Not only do the evolutions diverge at high energy, but for Sibyll 2.3 the relative behavior of the proton and pionair cross-section is different from the other models (faster increase of the pion-air cross-section).

Multiplicity
According to [50], the multiplicity plays a similar kind of role as the cross section, but with a weaker dependence (log). On the other hand the predictions from the models have larger differences for the multiplicity compared to the cross section. As shown in Fig. 6 (left-hand side) and in a more general way in [59], the average multiplicity is well reproduced by all the models up to 1 TeV and even up to 13 TeV for EPOS LHC and QGSJETII-04 [60] and a difference appears between these two models only at the highest energy (beyond 100 TeV). However, in the case of a nuclear target the slope of the rise of the multiplicity as function of the energy is different for all three models leading to a difference of about 20-30% at the highest energies in p or π -air interactions (Fig. 6 right-hand side). This effect is small compared to the pre-LHC era [51] but can change the elongation rate of the air shower maximum development. Here again Sibyll 2.3 has a different behavior than the other models with a smaller slope and the same multiplicity for p or π -air interactions while other models have about 10% difference.
On Fig. 7 left-hand side it can be seen on the pseudorapidity distribution of charged particles at 7 TeV that even if the central density of particles is well reproduced by all models, the width of the distribution is too narrow in the case of Sibyll 2.3 which leads to a reduced total multiplicity as seen on Fig. 6. Furthermore on the right-hand side of Fig. 7, we can observe that the fluctuations are very similar for QGSJETII-04 and EPOS LHC but again Sibyll 2.3 seems to have problems to reproduce the shape of the distribution. This can be important for the fluctuations of the air shower maximum.
Finally by comparing the multiplicity of p-air, heliumair and iron-air interactions we can check that the difference between models increases due to nuclear effects. On Fig. 8 it is clear that from a maximum difference of about 30% for p-air interactions, one can reach more than a factor of 3 in the case of Fe-air. This is of course important for the simulation of air showers with a primary mass heavier than protons.
So, for both cross section and multiplicity, when the models are constrained by LHC data up to 7 TeV, the extrapolation to the highest energy in p-p are very similar but differences remain in nuclear and pion interactions because of the lack of data at high energy and with light ions (only heavy ion data available from RHIC and LHC at high energy).

Diffraction and elasticity
Another important observable determining air shower development is the elasticity [50] defined as the largest energy fraction carried by a secondary particle (the leading particle). The model predictions are shown on Fig. 9 for p-p, π -air and p-air (as inelasticity=1-elasticity) as a function of center of mass energy. Sibyll 2.3 has the largest elasticity and it is probably related to the fact that the multiplicity is lower (less energy taken from the leading particle). In the cases of EPOS LHC and QGSJETII-04 the difference is smaller for an air target compared to pp interactions. This opposite behavior compared to the other observables can be explained by the fact that this quantity is very difficult to measure in collider experiments since the latter cannot cover 100% of the phase space. As a consequence there are only indirect constraints on the different contributions to the elasticity leading to a larger uncertainty in the models.    One contribution is the diffractive dissociation. Indeed diffractive events are producing the largest elasticities and are important for air shower development, not only for the position of the shower maximum but also for the muon production [63]. At the LHC various measurements related to diffraction are now available [55,[64][65][66][67][68]. Due to the difficulties of measuring very forward particles, the compatibility between the results is not as good as it is for the mid-rapidity measurements. This leads to some uncertainties in air shower simulations at a level of 10 g/cm 2 [69]. Nevertheless the difference between models seems to be even larger as illustrated on Fig. 10. The rapidity gap (range in pseudorapidity without particle detection in triggered events) cross-section measurement is poorly described by the model while it is directly related to the elasticity in general and diffraction in particular (the large rapidity gaps come from single diffractive events). For instance, the large probability for Sibyll 2.3 to produce rapidity gap around 2 to 4 is a direct consequence of the too narrow pseudorapidity distribution and implies a large elasticity.
It is interesting to notice that hadronic interaction models used for air shower simulations underestimate the diffraction dissociation as observed at LHC (as on Fig. 10 for large rapidity gaps). In fact in EPOS LHC this was improved, leading to a contradiction with air shower data as shown in Sect. 4.3.

Baryon and resonance production
Another important observable for EAS is the number of muons reaching the ground. It has been shown in [70] that the production of particles which are not π 0 (for  instance baryon-antibaryon pairs or ρ 0 resonance) plays an important role in the muon production rate especially if we take into account the leading particle effect [71]. Recent measurements by NA61 [72] show that the ρ 0 production in π -C interactions seems to be underestimated by a relatively large amount (from 20% to 100%) potentially leading to a large increase of muon production [73]. On the other hand in [74] it is demonstrated that increasing the muon production by increasing the forward baryon pair production like in EPOS leads to a very deep muon production which seems to be in contradiction with data (see Sect. 4.3). Indeed from [75] it can be concluded that the excess of protons seen in [41] (Fig. 3) is not due to newly produce baryons but is due to some baryon stopping (protons from the nuclear target). As a consequence this effect does not lead to an increase of muon production by energy transfer as in EPOS LHC. Both results imply a change in the hadronic interactions models with strong implications on muon production in air showers as shown in the next section.

Depth of maximum shower development
As shown in Fig. 11, the mean depth of shower maximum, X max , for proton and iron induced showers simulated with CONEX is different for EPOS LHC, QGSJETII-04 and Sibyll 2.3 as a direct consequence of the differences shown in Sect. 3. However the elongation rate (the slope of the mean X max as a function of the primary energy) is almost the same for all models since the difference between models is now much lower than it was in the past [51]. The difference between the models is a constant shift of about + / − 20 g/cm 2 around the value given by EPOS LHC. From the results shown in Sect. 3 it is likely that Sibyll 2.3 predicts too large values of the mean X max since the multiplicity is already too low and the elasticity too high at the LHC.
Nevertheless the very similar elongation rate is very important for the study of the primary cosmic ray composition. If the models converge to a similar elongation rate, it will allow us to have a more precise idea on possible changes in composition at the "ankle", for instance, where the Pierre Auger Observatory measures a break in the elongation rate of the data [76].
In fact, further studies using the fluctuations of X max around the mean can be used to test model consistency. Indeed both mean X max and X max fluctuations depend on the mass composition and since fluctuations are less dependent on the details of hadronic interactions (superposition model [50]) than the mean value, it can be checked that the composition corresponding to a given mean X max is consistent with the observed fluctuations. In [76] the Pierre Auger collaboration shows that while it is possible to describe the observed data with EPOS LHC, QGSJETII-04 is in tension with data at a 1 sigma level (mean X max too shallow by ∼ 15 g/cm 2 ).

Muons at ground level
Concerning the number of muons at ground level (for 40 • inclined showers at a height of 1500 m), the difference between QGSJETII-04, EPOS LHC and Sibyll 2.3 is relatively small. We can see on Fig. 12 that model predictions differ only by about 10%. The studies by the Pierre Auger Observatory show that the absolute number of muons observed in vertical shower differ from the model predictions by 1.33 ± 0.13 ± 0.09 [22] in the best case. This is a 2 sigma effect and in case of inclined showers the effect is less than 2 sigmas too [21]. Taking into account the ρ 0 measurement as explained in Sect. 3.4, it is not unlikely that the next generation of hadronic interaction models can reproduce the absolute number of muons, at least for vertical showers.
Even if the number of muons is much more similar now for all recent hadronic interactions models, and not so different compared to the data, there is still a large uncertainty related to the energy spectrum of the produced muons. This is an important factor for the attenuation length of the muons in the atmosphere [77] and for the muons at ground level in general [78]. As a consequence, one of the most sensitive measurement of how muons are produced in an air shower is the muon production depth and this is, in fact, not well reproduced by the current models.

Muon production depth (MPD)
We have seen in the previous section how LHC data could improve the description of EAS using updated hadronic interaction models. In fact, in one particular case, the update of EPOS leads to inconsistent results: the muon production depth measured by the Pierre Auger Observatory [20]. As shown on Fig. 13 the mean logarithmic mass ln A calculated from X µ max is incompatible with the one extracted from X max and even out of the range defined by the proton and iron primary mass when EPOS LHC is used for the simulation. With QGSJETII-04 the resulting ln A from X µ max is below the iron line but not consistent with the one from X max . In a previous analysis [79], EPOS 1.99 was giving results lighter than iron, so the important shift observed in the MPD simulated with EPOS LHC can partially be explained by the change in elasticity due to the corrections in diffractive interactions needed to reproduce the rapidity gap distributions measured by the ATLAS collaboration [68]. We can see on Fig. 10 that EPOS LHC gives reasonable results while QGSJETII-04 is too low.
The change of the parameters needed to describe the rapidity gap correctly (the diffractive cross-section and the diffractive mass distribution) affected both proton and pion interactions because the same parameters were used for both type of projectile. While the change of diffraction and thus of elasticity in proton interaction has very little impact on X µ max , it appears that MPD are extremely sensitive to the elasticity of pion interactions. This can be understood by the fact that muons are produced at the end of the hadronic cascade after many generations of mainly pion-air interactions. As a consequence of this cumulative effect, even a small increase of only about 10% of the elasticity of pion-air interactions can lead to large shift in X µ max .
To check this hypothesis, the diffractive cross-section for pion interactions has been reduced in EPOS LHC to get a reduction of about 10% of the elasticity of the pionair interactions. As a result X µ max is reduced by about 20 g/cm 2 as shown on Fig. 14 with the pink line with open circles. The diffraction has not been changed for proton interactions to keep full compatibility with LHC data and then the change in X max is limited to less than 5-10 g/cm 2 . Another consequence is the increase of the number of muons at ground level by few percents.
Such a small change is compatible with all pionnucleus data that are available at low energy and thus these two versions of EPOS cannot be discriminated from accelerator data. But the effect on MPD is so strong that data from the Pierre Auger Observatory can be used to constrain diffraction in pion interactions to get consistent results between the mean logarithmic mass which can be extracted from X µ max and the one deduced from X max which has very little dependence on pion hadronic interaction [74]. From EAS development we can thus say that elasticity of pion-air interactions should be lower than the elasticity of proton-air interactions.
The second factor explaining the large shift in MPD was identified in [74] as the too large production of forward baryons in pion interactions (which was indeed extended from low energy only in EPOS 1.99 to all energies in EPOS LHC to improve model consistency). Simply suppressing the production of diquark in string ends and thus the forward baryon pair production as in Fig. 3, the resulting X µ max is shown on Fig. 14 as a thin black line (on top of Sibyll 2.3 predictions) which is again about 20 g/cm 2 lower than the original EPOS LHC predictions.
The electromagnetic X max is increased by less than 5 g/cm 2 by the change of forward baryon production (more energy in the π 0 ). The muons production is reduced at the level of QGSJETII-04.
Since these two effects are cumulative, changing both diffraction and forward baryon production leads to a value of X µ max very similar to the one from QGSJETII-04 and thus compatible with Auger X µ max [20] with a X max reduced by only 5 g/cm 2 still compatible with the Auger X max [76]. The decrease of muon production could be compensated by a larger ρ 0 production which is in fact related to the increase of diffractive dissociation in pionnucleus interactions. This is the way to follow for future model development.

Summary
In [74] the uncertainty in the first proton(nucleus)-air interaction has been identified as the source of 70% of the uncertainty in the simulated X max . The remaining 30% being linked to the pion-air interactions. Concerning the muon production, 90% is coming from the pion interactions and only 10% from the first interaction. In Sect. 3 we have shown that for the first interaction the uncertainty is not in the basic p-p interaction any more, very well constrained by LHC data, but by the nuclear effects which can not be tested properly with current model and data combination (data with heavy ions only at high energy and only EPOS LHC can treat heavy ion collisions properly). These nuclear effects, being important

ISVHECRI 2016
for both the air target and in the case of heavier primaries, are the main source of the systematic shift in X max of about 20 g/cm 2 around EPOS LHC predictions. This uncertainty is comparable to the experimental uncertainty in the measurement of X max and the elongation rate is now the same for all models for a constant composition. As a consequence the interpretation of the data using a post-LHC model will be more reliable, especially concerning the possible change in mass composition with energy as summarized in [80].
To further reduce these uncertainties and improve the description of air showers by hadronic interaction models, in particular the observables based on muons, it is crucial to improve the description of pion-nucleus interactions in general and the diffractive dissociation in particular which is likely to be different than in proton interactions. Forthcoming studies of diffraction at the LHC including with a nuclear target [15,81] will reduce the model uncertainty for the first interaction at its minimum. To further improve the models it is important to take into account that air shower measurements, such as the muon production depth, can also give very strong constraints on hadronic interactions in particular for pion interactions [74] for which cumulative effects due to the hadronic cascade can be observed. This should give qualitative input to improve the models which can then be quantitatively tested against past and future NA61 measurements for instance [82].
To conclude, we can say that LHC data contribute a lot to reducing the uncertainties in air shower simulations, providing better tools to analyze cosmic ray data. The differences between the hadronic models have been reduced but one should keep in mind that there are still uncertainties in the models themselves which have to be better quantified and transferred to the calculation of the systematic errors in EAS analysis. Consistency of different EAS observables can and should be used to test the hadronic interaction models. The open issues are now mainly in the treatment of pion interactions which have a direct influence on the geometry and energy of the muons in air showers. The next generation of models taking into account more detailed LHC data and what has been learned from the MPD study should significantly improve their description of air showers.
The author would like to thank the Pierre Auger Collaboration and Sergey Ostapchenko for useful discussions.