High energy cosmic ray interactions and UHECR composition problem

The differences between contemporary Monte Carlo generators of high energy hadronic interactions are discussed and their impact on the interpretation of experimental data on ultra-high energy cosmic rays (UHECRs) is studied. Key directions for further model improvements are outlined. The prospect for a coherent interpretation of the data in terms of the UHECR composition is investigated.


Introduction
Among the most challenging problems of the high energy cosmic ray research is the enigma of ultra-high energy cosmic rays (UHECRs). Despite a long history of UHECR studies, one is still quite far from understanding their potential sources and acceleration mechanisms. On the experimental side, because of the very low UHECR flux, one has to rely on indirect detection techniques: measuring the properties of nuclear-electromagnetic cascades, so-called extensive air showers (EAS), initiated by interactions of the primary cosmic ray (CR) particles in the upper atmosphere [1]. In turn, a correct interpretation of the respective data requires a good quantitative understanding of the EAS development.
A substantial progress in the latter direction has been achieved with the development of contemporary Monte Carlo (MC) tools for air shower simulations, like the COR-SIKA [2], AIRES [3], SENECA [4], and CONEX [5] programs. One of the key ingredients of those numerical codes are MC models of high energy collisions of hadrons and nuclei, which describe interactions of both primary CRs and of the produced secondary hadrons with nuclei of the air. However, since such collisions generally involve nonperturbative physics, no first principle description can be employed for them and the corresponding models are largely phenomenological ones. Therefore, further progress in the understanding of ultra-high energy cosmic rays is closely linked to accelerator studies of high energy collisions of hadrons and nuclei, which are required both for calibrating parameters of such hadronic interaction models and, more importantly, for discriminating between the respective theoretical approaches.

Impact of LHC data and further progress
Valuable information on the properties of high energy proton-proton, proton-nucleus, and nucleus-nucleus inter-⋆ e-mail: ostapchenko@fias.uni-frankfurt.de actions has been provided by experimental studies at the Large Hadron Collider (LHC). In particular, of the highest importance for the high energy cosmic ray field proved to be measurements of the total and elastic pp cross sections by the TOTEM [6][7][8][9][10] and ATLAS experiments [11,12]. First of all, the corresponding results allow one to constrain seriously model predictions for proton-air and nucleus-air inelastic cross sections and thereby, for the mean free path (m.f.p.) of primary CR particles in the atmosphere. This is related, in turn, to the position of the starting point for air shower development and, thus, impacts all EAS observables, most notably, the shower maximum depth X max . Secondly, experimental results on the energy-dependence of σ tot pp constrain the underlying theoretical approaches of hadronic interaction models and allow one to fix important model parameters.
Among the most commonly used MC generators of high energy cosmic ray interactions are EPOS [13,14], QGSJET-II [15][16][17][18], and SIBYLL [19][20][21][22]. All those models are inspired by the qualitative picture of quantum chromodynamics (QCD): hadron-proton, hadron-nucleus, and nucleus-nucleus collisions are mediated by cascades of partons [(anti-)quarks and gluons]. Since such cascades generally involve partons of small virtuality, one deals with nonperturbative phenomenological approaches. Nevertheless, the respective models have a considerable predictive power; using a particular approach, the energy dependence of the basic interaction properties, such as total and inelastic cross sections, multiplicity of secondary particles, etc., is a mere consequence of the increase of the phase space available for such cascades to develop. Moreover, moving from proton-proton to pion-proton interactions or to proton-nucleus and nucleus-nucleus collisions, one keeps the same interaction mechanism; the only thing which changes are the initial conditions for those parton cascades, i.e. the structure of the so-called constituent parton Fock states representing the interacting particles. Thus, fixing basic parameters of a particular approach, based on the measured proton-proton cross sec-tions, substantially constrains many other model predictions.
Nevertheless, despite the fact that all the abovementioned models have been updated with the data of Run I of LHC, their predictions for various characteristics of extensive air showers differ considerably. If we compare, for example, the energy-dependence of the average shower maximum depth, calculated using the newest model versions, the differences between the predicted X max values reach 40 g/cm 2 at the highest energies, as one can see in  Primary energy dependence of the average shower maximum depth for proton-and iron-initiated vertical EAS, as calculated using the QGSJET-II-04 [17,18], EPOS-LHC [14], and SIBYLL-2.3 [21] models (solid, dashed, and dashed-dotted lines respectively), compared to the data of the Pierre Auger Observatory (PAO) [23] (squares). experimental accuracy. It is thus highly desirable to reveal the reasons for those differences and to find ways to further constrain model predictions.
Starting first with the differences between the SIBYLL-2.3 results and the ones of the other two models, those appeared to be related to the very basic model assumptions concerning the structure of constituent parton Fock states in hadrons, i.e. to the above-mentioned initial conditions for parton cascades, as discussed in more detail in Ref. [24]. Like most of the hadronic event generators used in the collider field, the SIBYLL model is based on the "minijet" approach [25][26][27] which corresponds implicitly to the picture shown schematically on the left-hand side of Fig. 2. There, one starts from the same universal parton Fock state of the proton. Additional partons (sea quarks or gluons) giving rise to new branches of the parton cascade, which take part in multiple scattering processes, result from the evolution of the parton density corresponding to this initial state and their momentum fractions x are distributed as ∝ 1/x in the very high energy limit. Such a picture reflects itself in the hadron production pattern predicted by the model: multiple scattering affects mostly central particle production, while having a weak influence of forward hadron spectra. Indeed, the latter are formed by the hadronization of partons emerging from the initial part of the parton cascade, which starts from the same initial conditions and covers a short rapidity interval, being thus weakly dependent on the further development of the cascade.
A direct consequence of the above-discussed approach is a weak energy-dependence of the inelasticity K inel , i.e. the relative energy loss of leading (most energetic) hadrons produced. With increasing energy, one obtains a significant enhancement of secondary particle production in the central rapidity region only, which has a weak impact on the energy loss of leading nucleons in proton-proton and proton-nucleus collisions. As one can see in Fig. 3, the energy dependence of K inel pp is indeed almost flat for SIBYLL-2.3. In turn, a slower energy-rise of the inelasticity implies a larger EAS elongation rate, hence, a larger X max at sufficiently high energies (see, e.g. Ref. [28]), as we observed indeed in Fig. 1.
In the alternative approach, implemented in the EPOS and QGSJET-II models, any hadron is represented by an infinite sequence of Fock states containing different numbers of large x constituent partons, as shown schematically on the right-hand side of Fig. 2. Further cascading of those partons "dresses" them with low x parton clouds. As the overall multiplicity of (anti-)quarks and gluons in the central rapidity region is roughly proportional to the number of initial constituent partons, stronger multiple scattering is typically associated with larger Fock states. Thus, there is a strong long-range correlation between central and forward particle production; higher multiplicity in the central region reflects stronger multiple scattering. In turn, this implies that larger ensembles of constituent partons are involved in the process, which has a strong impact on forward hadron spectra.
This naturally leads to a substantial energy-rise of the inelasticity, which is clearly seen in Fig. 3 for EPOS-LHC and QGSJET-II-04. The reason for this rise is twofold. First, for any given Fock state, increasing multiple scattering implies that more and more constituent partons take part in the interaction, thus leaving smaller fractions of the initial proton momentum for spectator partons which finally form the leading nucleons. Additionally, Fock states containing bigger numbers of large x constituent partons come into play. Momentum sharing between these partons results in a smaller fraction of the initial proton momentum, possessed by each parton, which thus enhances the energy loss of the leading nucleons.
A crucial discrimination between the two approaches may be provided by measuring correlations between the signal strengths in central and forward-looking detectors at the LHC. For the particular case of the CMS and TOTEM experiments, this is illustrated in Fig. 4, where we plot for √ s = 8 TeV the pseudorapidity density of produced charged hadrons dn ch pp /d|η| at |η| = 6 (averaged over the interval 5.5 < |η| < 6.5) as a function of the central ηdensity (|η| < 1) of charged hadrons. Both EPOS-LHC and QGSJET-II-04 predict a strong correlation of the signal strengths in CMS and TOTEM. The respective results of the two models practically coincide with each other, apart from the tails of the multiplicity distributions. In contrast, for SIBYLL 2.3 such a correlation is twice weaker, being thus a "smoking gun" signature for the desirable discrimination.
In turn, the difference between the EPOS-LHC and QGSJET-II-04 predictions for X max is related to the treatment of forward hadron production in pion-air collisions, the main effect coming from a copious production of baryon-antibaryon pairs and from a higher rate of the inelastic diffraction in EPOS-LHC [29,30]. As discussed in detail in Ref. [29], this can be discriminated based on PAO measurements of the maximal muon production depth [31].    Generally, further progress concerning the reliability of model predictions for various properties of extensive air showers will come both from improving the calibration of model parameters (so-called model tuning), based on the presently available and forthcoming accelerator data, and from further development of the underlying theoretical approaches, aiming at a better self-consistency of the models. For cosmic ray applications, especially important is the treatment of forward hadron production in protonair and pion-air collisions. In this respect, one can benefit from recent results of the NA61 [32] and LHCf [33][34][35][36][37][38] experiments. In Ref. [18], the dominance of the pion exchange mechanism has been postulated for forward meson production in nondiffractive pion-proton and pion-nucleus collisions, in order to reach an agreement with measured spectra of π 0 and ρ 0 mesons. While this conjecture has been further confirmed by the NA61 results on ρ 0 production in pion-carbon interactions [32], which are consistent with previous measurements for π + p collisions [39], the signatures of the very same mechanism are now clearly seen in the neutron spectra for pp collisions at LHC energies: at large Feynman x, LHCf observed a pronounced "bump" in the forward neutron spectra, corresponding to the pion exchange mechanism [37,38]. Thus, a consistent model implementation of the mechanism is an urgent task. While the basic contribution to the process comes from the so-called Reggeon-Reggeon-Pomeron (RRP) diagram, the treatment of the respective absorptive corrections, which describe the suppression of this contribution with increasing energy [40][41][42], is a nontrivial problem in a Monte Carlo framework.

Composition of UHECRs: coherent interpretation of experimental data?
At present, we have a very puzzling situation with the studies of the composition of ultra-high energy cosmic rays. On the one side, there exists a wealth of experimental results of high statistical significance, concerning measurements of various properties of UHECR-initiated EAS, which are sensitive to the primary composition. The list is quite extensive and includes, in particular, the primary energy dependence of the average shower maximum position X max and its fluctuations RMS(X max ), of the maximal muon production depth X µ max , and of the muon density at ground ρ µ . On the other hand, a coherent interpretation of all those data sets, in terms of a unique UHECR composition, proves to be a very difficult challenge. This is related, in particular, to the fact that UHECR composition studies necessarily involve models of high energy hadronic interactions and that the relevant predictions of those models are yet rather different from each other. Moreover, none of the presently available models appeared to be able to provide an overall self-consistent interpretation of the available experimental results on UHECRs. Most gravely, they can not even describe specific data sets, like ρ µ for the highest energy events [43], for any reasonable assumptions concerning the primary composition.
One possible explanation of the above-mentioned problems is that the present models of high energy interactions do not provide an adequate enough description of the corresponding physics. In the following, we are going to test this assumption by checking whether a coherent interpretation of the various experimental results can be achieved applying some physically meaningful changes of the relevant model predictions. We shall concentrate on the data of the Pierre Auger Observatory, which have high statistics of UHECR events and contain measurements of the different EAS observables of interest.

"Marrying" X max and X µ max in terms of UHECR composition
Among the most composition-sensitive EAS observables is the average shower maximum depth X max . It is also a good anchor point for our analysis since the respective measurements of the Pierre Auger Observatory and the Telescope Array experiment are consistent with each other [44,45]. However, comparing the PAO data with model predictions for X max in Fig. 1, it is quite evident that the inferred UHECR composition will be highly modeldependent. To discriminate between the different models, we can use another EAS observable measured by PAOthe maximal muon production depth X µ max . The respective data clearly disfavor EPOS-LHC and SIBYLL-2.3 as the results of the two models require the primary cosmic rays to be much heavier than iron nuclei 1 [30,31]. Neverthe-1 A dedicated study performed in Ref. [29] showed that the deep X µ max of EPOS-LHC is caused by a too abundant (anti-)baryon production in pion-air interactions at high energies, predicted by that model [46]. In addition, the model is characterized by a rather high rate of inelastic diffraction for pion-air collisions, which also contributes to the elongation of the profile for muon production [30]. less, even using QGSJET-II-04, a consistent interpretation of the PAO data on X max and X µ max is not possible: while a comparison with the measured shower maximum position indicates a rather light primary composition, an agreement of the model prediction for X µ max with the measurements can only be reached if UHECRs are iron nuclei [31].
Let us now try to "marry" the PAO data on X max and X µ max composition-wise by changing the model predictions for high energy proton-air and pion-air collisions. It is not obvious a priori if this is possible since any modification of the model treatment, which affects X max , will also impact X µ max predictions, and vice versa, and both quantities will change in the same direction.
Let us first try to increase the predicted X max , in order to change the interpretation of the respective PAO data towards a heavier cosmic ray composition. It is worth reminding that the shower maximum position is most sensitive to the properties of the primary particle interaction in air, notably, to the respective inelastic cross section, to the rate of the inelastic diffraction, and to the inelasticity of the collisions [28]. It is important to stress, however, that potential changes of these characteristics of the first interaction will only impact the initial stage of EAS development, hence, give rise to an approximately parallel shift of the electromagnetic and hadronic longitudinal shower profiles, without any serious modification of the profile shapes. As a result, one will obtain a similar effect on X max and X µ max : ∆X max ≃ ∆X µ max . As discussed in Section 2, the inelastic proton-air cross section σ inel p−air is seriously constrained by measurements of the total and elastic proton-proton cross sections at LHC. In turn, changing the treatment of the inelastic diffraction, within present uncertainties on the diffractive pp cross section, allows one to change X max by less than 10 g/cm 2 [47]. Thus, the only significant freedom left is related to the energy-dependence of K inel p−air . In order to check the impact of a slower energy-rise of K inel p−air on model predictions for X max and X µ max , we shall use the "cocktail" method [29]: treat the primary particle interaction with the air, using the SIBYLL-2.3 model, but keep QGSJET-II-04 for the description of interactions of all secondary hadrons in the cascade. The so-obtained energy-dependences of X max and X µ max are compared in Fig. 5 with the respective results of QGSJET-II-04 and with the PAO data. 2 As one can see in the left panel of Fig. 5, a slower energy-rise of K inel p−air allows one to shift the position of the shower maximum somewhat deeper in the atmosphere and, thus, to change the interpretation of the PAO data on X max towards a slightly heavier UHECR 2 While the PAO results for X max have been corrected for the detector acceptance and the experimental reconstruction effects [23], this is not the case for X µ max . Therefore, the QGSJET-II-04 predictions for X µ max (solid lines in the right panel of Fig. 5) are taken from the original PAO publication [31], while the respective results of the above-described "cocktail" model (dashed lines in the right panel of Fig. 5) are obtained assuming that relative differences between X µ max results of any two models remain unaffected by the acceptance and reconstruction effects. Such a conjecture is supported by the fact that relative differences between X µ max values of EPOS-LHC and QGSJET-II-04 for p-and Fe-induced EAS in Ref. [31], the experimental acceptance and reconstruction effects included, coincide, within few g/cm 2 , with the ones obtained from pure CONEX simulations. QGSJET-II-04 SIBYLL for p-air, QGSJET-II -rest p Fe Figure 5. Primary energy dependence of X max (left) and X µ max (right) for proton-and iron-induced vertical EAS, as calculated using the QGSJET-II-04 model (solid lines) or applying SIBYLL-2.3 for the treatment of proton-air collisions and QGSJET-II-04 for the rest of the hadron cascade (dashed lines). PAO data [23,31] are shown by squares.
composition. But at the same time, the obtained maximal muon production depth, plotted by the dashed lines in the right panel of Fig. 5, also moves deeper and the respective interpretation requires primaries heavier than iron. Therefore, proceeding this way, we can not reach the desirable consistency between the two sets of experimental data.
Let us now try to shift the predicted X µ max higher in the atmosphere and thus change the interpretation of the PAO data on the maximal muon production depth towards a lighter UHECR composition. As discussed in Ref. [29], the most efficient way to do so is to change the treatment of pion-air and kaon-air interactions, e.g., by increasing the respective inelastic cross sections, by decreasing the rate of the inelastic diffraction, or by increasing the inelasticity K inel π−air and the multiplicity of secondary hadrons produced. Very importantly, whatever option we choose, it will impact every step in the hadronic cascade development, hence, its effect on X µ max will be cumulative over all those steps. For example, higher σ inel π−air and σ inel K−air would shorten the mean free path of every pion and kaon in the cascade, thus considerably speeding up the cascade development and giving rise to a significant shrinkage of the longitudinal hadron and muon production profiles, hence, to sizable upward shifts of the maxima of those profiles. A similar effect can be obtained by using a smaller diffractive cross section for pion-air collisions. Indeed, since in diffractive interactions the incident hadron typically looses a small portion of its energy, those may be regarded, in the first approximation, as quasielastic processes. Therefore, using a smaller σ diffr π−air is essentially equivalent to having a shorter pion m.f.p. Finally, an increase of the inelasticity and multiplicity of pion-air collisions would speed up the decrease of the average pion energy in the cascade, hence, most of the pions will reach the pion critical energy, for which the probabilities of pion interaction and decay become equal, at a smaller depth.
Very remarkably, the above-discussed changes of the treatment of pion-air collisions would have a much weaker impact on the average shower maximum depth since X max is only sensitive to pion-air and kaon-air interactions in the beginning of the hadron cascade, before most of the energy of the primary particle is channelled into secondary electromagnetic cascades [29].
Let us now check these conjections by applying again the "cocktail" method: keeping QGSJET-II-04 for the description of (anti-)nucleon interactions but treating pionair and kaon-air collisions by means of the old QGSJET model [48][49][50]. As σ inel π−air , K inel π−air , and, especially, the multiplicity of pion-air interactions have a steeper energy-rise in QGSJET, compared to QGSJET-II-04, we obtain indeed a much smaller X µ max , as one can see in Fig. 6 (right). At the same time, the shower maximum position goes just slightly higher in the atmosphere, as shown in Fig. 6 (left), such that the two PAO measurements can be consistently interpreted in terms of a light UHECR composition.
Thus, in this subsection we obtained two important results. First, changing the treatment of high energy pion-air collisions, one is able to reach a consistent interpretation of the PAO data on the average shower maximum depth and the maximal muon production depth. Secondly, proceeding in such a way, we inevitably arrive to a very light, predominantly proton-dominated, UHECR composition.

What about the small RMS(X max )?
Let us now check if small X max fluctuations observed by PAO for highest energy UHECRs [23] can be made compatible with a light primary composition. Obviously this is not the case for the present high energy interaction models, see Fig. 7. Thus, we have to address the question whether there exists a physically sensible way to substantially reduce the predicted RMS(X max ). QGSJET-II-04 QGSJET for π-air p Fe p Fe Figure 6. Primary energy dependence of X max (left) and X µ max (right) for proton-and iron-induced vertical EAS, as calculated using the QGSJET-II-04 model (solid lines) or applying QGSJET-II-04 for the treatment of (anti-)nucleon-air collisions and QGSJET for the description of pion-air and kaon-air interactions (dashed lines). PAO data [23,31] are shown by squares.  In case of proton-induced EAS, the largest contribution to RMS(X max ) comes from fluctuations of the free path of the primary proton [51], which is related to the inelastic proton-air cross section. As already discussed above, there is very little freedom for changing the predicted σ inel p−air since it is strongly constrained by measurements of the total and elastic pp cross sections. Additionally, RMS(X max ) depends somewhat on the rate of the inelastic diffraction in proton-air collisions. However, modifying the diffraction treatment within present uncertainties on σ diffr pp , one obtains only few g/cm 2 effect for RMS(X max ) [47]. Finally, the quantity depends on the Glauber(-Gribov) geometry for nondiffractive proton-air collisions [52,53], mainly, on fluctuations of the impact parameter and the related variations of the number of target nucleons taking part in the interaction. Here, one also has little freedom, given the known nuclear transverse profiles and the effective "size" of the proton, constrained by measurements of the differential elastic pp cross section at LHC. The overall uncertainty for calculations of RMS(X max ) for p-induced EAS is well represented by the differences between the current model predictions, which are quite small, as one can see in Fig. 7.
Coming now to the case of primary nuclei, fluctuations of the free path of primary particles are of secondary importance, since the corresponding m.f.p. is short due to the large nucleus-air cross sections, especially, for relatively heavy nuclei. The same is true for the inelastic diffraction since the diffraction rate in nucleus-nucleus interactions is much smaller than in the pA case [48]. Thus, the dominant contribution to RMS(X max ) comes here from the Glauber(-Gribov) geometry of nucleus-nucleus collisions, notably, from fluctuations of the impact parameter and, correspondingly, of the number of projectile nucleons taking part in the interaction [48,54]. While this is well-constrained, given the known nuclear profiles, a potential source of uncertainty for X max fluctuations is related to the fragmentation of the "spectator part" of the incident nucleus, i.e. to the treatment of the projectile nucleons not taking part in the inelastic interaction. For the two extreme cases, when the whole spectator part is kept as a single secondary nucleus or when it is broken into separate nucleons, the respective RMS(X max ) differs by a factor of two [48].
Since experimental data on the production of nuclear fragments are available at fixed target energies only, this may pose a problem for predicting nuclear fragmentation at very high energies. Fortunately, relative fragment yields scale above few GeV per nucleon (see, e.g., Ref. [55] for a review), which allows one to calibrate fragmentation pro-cedures with low energy data and then to extrapolate them reliably to ultra-high energies. Thus, also in the case of nuclear primaries, there is very little freedom to change the model predictions for RMS(X max ). This is well illustrated by the excellent agreement between the respective results of QGSJET-II-04 and SIBYLL-2.3 for different primary nuclei, see In view of the discussion above, it is quite surprising that the EPOS-LHC model predicts substantially smaller RMS(X max ) for nuclei-initiated EAS, compared to QGSJET-II-04 and SIBYLL-2.3, especially, for intermediate mass nuclei. To clarify this problem, let us compare in Fig. 9 the default EPOS-LHC results for Fe-induced EAS to calculations of QGSJET-II-04 and SIBYLL-2.3 for the two above-mentioned extreme cases: when the nuclear spectator part is kept as a single nucleus or when it is broken into separate nucleons. It is easy to see that the EPOS-LHC predictions are very close to the full break-up option, irrespective the primary energy, which clearly indicates that the small RMS(X max ) predicted by that model originates from an underestimation of the production of relatively heavy nuclear fragments. 3 The main conclusion of this subsection is that the actual uncertainties for calculations of RMS(X max ) are very small, as, in fact, has been stated earlier in Ref. [57]. Thus, "marrying" experimental data on X max fluctuations with the ones on the average shower maximum depth and the maximal muon production depth is a purely experimental challenge.

"Muon excess" problem
The most outstanding problem with the interpretation of UHECR data is a significant excess of the EAS muon 3 A likely origin of the problem is an incorrect matching between the interaction and nuclear fragmentation procedures in EPOS-LHC [56]. content for very high energy primaries, compared to current EAS simulation results: the observed muon density at ground level is substantially higher than model predictions for iron-induced air showers [43,58,59]. Let us first address the question if this contradiction can be eliminated by a change of the treatment of the primary particle interaction. To obtain a substantial enhancement of the muon number N µ , one has to consider a very radical change of the interaction mechanism, producing an order of magnitude increase of the number of secondary hadrons [28,60]. Such a change necessarily involves physics beyond the Standard Model. However, to have a significant impact on the average EAS muon content, such new physics should impact the bulk of UHECR-induced air showers, hence, to emerge with barn-level cross sections [61]. It is noteworthy that such a speculative scenario would give rise to a substantial increase of fluctuations of muon density at ground [61], which can be easily discriminated by current UHECR experiments. Without invoking an exotic scenario, to obtain a substantial muon excess at very high energies, one has to modify the treatment of pion-air collisions over a wide energy range [60]. In such a case, the final effect on N µ is a cumulative one over many generations of hadrons in the atmospheric nuclear cascade. Thus, the process is described qualitatively by the Geitler's picture [62,63] and we have the usual power-like energy-dependence of N µ : with E ref being some reference energy. The desirable modification amounts to an increase of the corresponding effective exponent: α µ →α µ . Since no muon excess has been observed at energies around the knee of the CR spectrum [64], we can safely set E ref = 1 PeV. Now, if at 10 EeV one has an excess R (10 EeV) µ = 10000α µ −α µ , there should be R (10 EeV) µ excess at 10 17 eV, which is at variance with the data of the KASCADE-Grande [65] and Ice-Top [66,67] experiments.
While we thus see no sensible way to explain the observed muon excess by changing the treatment of high energy hadronic interactions, one may investigate other possibilities. In particular, since most of the relevant muon density measurements have been performed at large distances from the shower core, a somewhat flatter shape of the muon lateral distribution function (LDF) could improve the agreement with the data. It is worth reminding that the muon LDF shape depends mostly on muon multiple scattering in the atmosphere and on the transverse momentum p t spectra of pions produced in pion-air collisions at relatively low energies (E 100 GeV) [68]. In particular, at large core distances (R core > 1 km), one is quite sensitive to the tails of the above-mentioned p tdistributions, which gives rise to a substantial dependence of predicted muon LDF on low energy hadronic interaction models [69]. Apart from that, one may think about improving the treatment of the Landau (large angle) muon scattering in current EAS simulation programs, since this mechanism may also influence significantly the shape of muon LDF at large distances.

Outlook
We discussed here the state of the art concerning the description of cosmic ray interactions at very high energies. The main conclusion may be formulated as follows: while we need hadronic interaction models for interpreting UHECR data, those models will remain phenomenological ones and, thus, improvements of the model treatment will always be an issue. The present Monte Carlo generators of high energy hadronic interactions considerably benefited from experimental studies of proton-proton, proton-nucleus, and nucleus-nucleus interactions at the Large Hadron Collider, especially, from LHC data on the total and elastic proton-proton cross sections. Nevertheless, there remain substantial differences between model predictions for various EAS characteristics. While some of those differences can be explained by evident drawbacks of certain models, additional dedicated accelerator studies are needed to discriminate between alternative theoretical approaches and thereby reduce other sources of uncertainty.
We also investigated the possibility to obtain a consistent interpretation of the data on various air shower properties, in terms of the UHECR composition, by considering some plausible changes of the predictions of high energy hadronic interaction models. Changing the treatment of pion-air collisions, one is able indeed to "marry" the PAO data on the average shower maximum depth and the maximal muon production depth, which leads one to a very light, predominantly proton-dominated, cosmic ray composition at very high energies. However, no sensible changes of the models of high energy hadronic interactions would allow one to bring the PAO data on fluctuations of the shower maximum depth and on the EAS muon density in agreement with the ones on X max and X µ max . Thus, reaching an overall self-consistent composition interpretation of the UHECR data requires to address systematic issues in the current experimental measurements and, possibly, to improve other parts of air shower simulation procedures, e.g., concerning the treatment of hadronic interactions at relatively low energies.