Cosmic Ray Interaction Models: an Overview

I review the state-of-the-art concerning the treatment of high energy cosmic ray interactions in the atmosphere, discussing in some detail the underlying physical concepts and the possibilities to constrain the latter by current and future measurements at the Large Hadron Collider. The relation of basic characteristics of hadronic interactions to the properties of nuclear-electromagnetic cascades induced by primary cosmic rays in the atmosphere is addressed.


Introduction
Experimental studies of high energy cosmic rays (CR) are traditionally performed using the so-called extensive air shower (EAS) techniques: the properties of the primary CR particles are reconstructed from measured characteristics of nuclear-electro-magnetic cascades induced by their interactions in the atmosphere. This naturally implies the importance of detailed Monte Carlo simulations of EAS development, particularly, of its backbone -the cascade of nuclear interactions of both the primary particles and of the secondary hadrons produced. Thus, the very success of experimental studies depends crucially on the validity of hadronic interaction models used in the analysis.
Typically, one chooses between two main experimental procedures [1]. In the first case, one deals with the information obtained with scintillation detectors positioned at ground. The energy of the primary particle is reconstructed from the measured lateral density of charged particles (mostly, electrons and positrons) while the particle type is inferred from the relative fraction of muons, compared to all charged particles at ground. Alternatively, one may study the longitudinal EAS development by measuring fluorescence light produced by excited air molecules at different heights in the atmosphere, for which purpose dedicated fluorescence telescopes are employed. In the latter case, the primary energy is related to the total amount of fluorescence light emitted. In turn, the particle type may be determined from the measured position of the so-called shower maximum X max -the depth in the atmosphere (in g/cm 2 ) where the number of ionizing particles reaches its maximal value (which is thus the brightest spot on the shower image).
Not surprisingly, the observables used to determine the primary particle type -the lateral muon density ρ μ and the EAS maximum position X max -appear to be very sensitive to details of high energy hadronic interactions. More precisely, X max depends strongly on the properties of the primary particle interaction with air nuclei: the inelastic cross section and the forward spectra of secondary hadrons produced. In turn, the EAS muon content is formed in a multi-step cascade process, driven mostly a e-mail: sergei@tf.phys.ntnu.no by interactions of secondary pions and kaons with air. Hence, ρ μ depends strongly on the properties of pion-air collisions, most importantly, on the multiplicity and spectral shape for charged hadrons produced.
The discussion above helps one to understand the basic requirements to cosmic ray interaction models: a treatment of various hadron-nucleus and nucleus-nucleus interactions over a wide energy range, including the respective inelastic cross sections and particle production, with a special importance of forward particle spectra. By far, the most crucial requirement is the one for a significant predictive power, which stems from the fact that there is a little possibility for model tuning based on air shower data. Rather, cosmic ray data allow one to check the overall consistency of a model approach in question.
Among the models most frequently used in cosmic ray studies have been QGSJET [2] and SIBYLL [3], as well as more recent ones: QGSJET-II [4,5] and EPOS [6,7]. Overall, CR interaction models behaved rather well when confronted with the first set of data from Run 1 of LHC [8]. QGSJET-II and EPOS had been further retuned with LHC data, and an update of the SIBYLL model is in progress [9].

Underlying physical approaches
By construction, CR interaction models are designed to cope with both nonperturbative "soft" and hard processes. QGSJET(-II) and EPOS employ the so-called "semihard Pomeron" approach [10,11], treating both contributions within the Reggeon Field Theory framework [12]. Basically, one introduces a "soft-hard" separation scale Q 2 0 , describing parton evolution in the perturbative range of high virtuality |q 2 | > Q 2 0 by means of DGLAP equations. In turn, for "soft" (|q 2 | < Q 2 0 ) parton cascades (and for a "soft pre-evolution" for hard cascades) a phenomenological soft Pomeron amplitude is employed. The total amplitude of a "general Pomeron" which contains both kinds of processes is used as a basic building block to develop a Pomeron calculus and to derive all partial cross sections and probabilities for particular final states [4,6,11].
On the other hand, SIBYLL is similar to most of the Monte Carlo generators used in the collider physics, being based on the "minijet" approach [13]. The latter employs an eikonalization of the inclusive cross section σ jet pp (s, p cut t ) for the production of parton jets of transverse momenta p t > p cut t . Partial probabilities for having a given number of jet pairs for some impact parameter b are expressed via the eikonal function defined by a product of σ jet pp (s, p cut t ) and a parton overlap function A(b). This is further supplemented by a simplified treatment of soft processes.
There are substantial similarities between the results of the two approaches concerning the description of inclusive hadron spectra in the central rapidity region. This is, first of all, due to a calibration of the models to similar sets of accelerator data. Moreover, in the very high energy limit, the bulk of central particle production comes from a hadronization of (mini-)jets, thus being driven by the input parton momentum distribution functions and DGLAP evolution.
However, the predictions of the two approaches differ significantly in the fragmentation region. In the "semihard Pomeron" scheme, multiple scattering affects both central and forward production of hadrons, creating thus long-range rapidity correlations. This gives rise to a noticeable violation of Feynman scaling in the fragmentation region and of the limiting fragmentation. Also the plateaulike behavior of the pseudorapidity η spectra of secondary hadrons extends to large values of η. In contrast, in minijet-like models, multiple scattering affects mostly central hadron production, being almost decoupled from the fragmentation region. As a consequence, the above-discussed effects are practically absent in that case.
served η-dependence of dn/dη in inelastic and nondiffractive event samples has been correctly reproduced by EPOS and QGSJET-II while other Monte Carlo generators in the study, e.g. various tunes of PYTHIA [15], fell short of describing the high yield of secondary hadrons in the TOTEM detector. Another important piece of evidence comes from studies of forward production of neutral pions by the LHCf experiment [16]: due to their characteristic scaling-like behavior, π 0 spectra of the SIBYLL and PYTHIA models appear to be too hard compared to the experimental observations. However, a more reliable discrimination between the two basic approaches requires to study correlations between central and forward particle production, which may be realized in combined measurements by e.g. ATLAS and LHCf detectors. Indeed, triggering different levels of activity in ATLAS, for example, measuring at least 1, 5, or 10 charged hadrons of p t > 0.5 GeV at |η| < 2.5, one could be able to discriminate between predictions of models of the two discussed types, as illustrated in Fig. 1 for a particular case of QGSJET-II-04 [5] and SIBYLL 2.1 [3] models. Similar studies may be performed with the CMS and TOTEM detectors.

Uncertainties of model predictions for X max
By far, the most suitable EAS parameter for studying primary CR composition is the shower maximum depth X max . Apart from the possibility to measure it reliably by modern air fluorescence detectors, the uncertainties of the respective model predictions have been greatly reduced with the start of the Large Hadron Collider -primarily, thanks to the precise measurements of the total and elastic proton-proton cross sections by the TOTEM and ATLAS experiments [17,18].
Yet another potential source of uncertainty for X max is related to its sensitivity to the rate of inelastic diffraction in proton-proton and proton-nucleus collisions. Meanwhile, there exists presently a serious tension between the results of diffraction measurements at LHC by the TOTEM and CMS experiments, as discussed in [19,20]. The potential impact of the latter on the calculated shower  maximum position has been studied in [19] in the framework of the QGSJET-II-04 model. In more detail, one considered two alternative model parameter tunes which allowed one to reproduce either the TOTEM or CMS results on single and double diffraction cross sections, while leaving essentially unchanged the model results for σ tot/el pp and for central particle production. Applying these alternative model tunes to EAS simulations, one observed that the difference between the respective X max predictions does not exceed 10 g/cm 2 , being thus comparable to the characteristic uncertainty of shower maximum measurements by modern EAS experiments.
However, present differences between various calculations of X max are substantially bigger, e.g. reaching 20 g/cm 2 in the ultra high energy limit in the particular case of the QGSJET-II-04 [5] and EPOS-LHC [7] models, the corresponding X max predictions being shown by red and blue solid lines respectively in Fig. 2 (left). This is rather surprising, taking the fact that both models have been recently updated using LHC data. Thus, the question arises if the analysis of Ref. [19] was not general enough or the position of the shower maximum depends on some other characteristics of hadronic interactions, not well constrained by the present LHC data.
To reveal the interaction features which are responsible for the above-discussed differences in X max predictions, we are going to use a "cocktail" model approach, using QGSJET-II-04 results for some selected interactions of hadrons in the atmospheric cascades or for some particular features of the primary interaction, while treating the rest with the EPOS-LHC model. As the first step, we apply QGSJET-II-04 to determine the position of the primary particle interaction in the atmosphere and to describe the production of secondary nucleons in this interaction; all other characteristics of the first proton-air collision and all the subsequent interactions of secondary hadrons in the cascade are treated using EPOS-LHC. This way we check the sensitivity of the calculated shower maximum depth to the model differences concerning the proton-air cross section and the predicted spectra of leading nucleons (the latter thus comprise the effects of the inelastic diffraction). The obtained X max shown by the green dot-dashed line in Fig. 2 (left) differs from the original EPOS-LHC results by not more than 5 g/cm 2 , which is well within the uncertainty range obtained in [19]. Further, we apply QGSJET-II-04 to describe all the characteristics of the primary interaction while treating the rest of the hadron cascade using EPOS-LHC. The obtained X max shown by the yellow dashed line in Fig. 2 (left) is shifted further towards the QGSJET-II-04 results by up to 5 g/cm 2 . This additional shift is explained by somewhat softer spectra of secondary mesons, most importantly, of secondary pions in QGSJET-II-04, compared to the EPOS-LHC model. Here we actually observe an important change in the physics of the hadronic cascade in the atmosphere. At lower energies, there is a very pronounced "leading nucleon" effect, i.e. the most energetic secondary particles in proton-air collisions are typically protons and neutrons (produced either directly or via decays of hyperons and resonances). On the other hand, in the very high energy limit the energy loss of leading nucleons is noticeably higher and the most energetic secondary hadron may well be a pion or a kaon.
What is the reason for the remaining difference with the QGSJET-II-04 results? One potential explanation could have been a difference between the two models concerning the predicted pion-air cross section or/and concerning the multiplicity and spectral shape for secondary mesons. However, the real reason is a different and a more interesting one. If we apply QGSJET-II-04 to describe both the primary interaction and the production of nucleons and antinucleons in all the secondary hadronair collisions, while treating the rest with EPOS-LHC, the obtained X max shown by the black dotted line coincides with the QGSJET-II-04 results. Thus, it is the more copious production of baryonantibaryon pairs and harder (anti-)baryon spectra in EPOS-LHC 1 [21] that gives rise to the largest differences in X max predictions of the two models in the very high energy limit.

Constraining model predictions for X max by cosmic ray data
As demonstrated in Section 3, a large part of the model uncertainty for X max is related to the treatment of proton-air interactions at very high energies, in particular, to the predicted diffraction rate and the forward spectra of secondary mesons. There are serious hopes that this will be strongly constrained by further studies of proton-proton and proton-nucleus collisions at LHC. Of particular importance are the measurements of very forward production of neutral pions by the LHCf experiment [16], which already demonstrated that the spectra predicted by EPOS-LHC are somewhat harder than the observed ones. However, another serious uncertainty is related to the production of baryon-antibaryon pairs in pion-air interactions. While more precise measurements of (anti-)baryon spectra at fixed target energies may be performed by the NA61 experiment [23], the prospects for extending such studies to much higher energies are rather dim.
Nevertheless, the discussed mechanism can be constrained indirectly by characteristics of EAS induced by very high energy cosmic rays. Recently, Pierre Auger collaboration reported measurements of the so-called maximum muon production depth in air showers, X μ max -the depth in the atmosphere (in g/cm 2 ) where the rate of muon production via decays of pions and kaons reaches its maximal value [24]. One observed a strong contradiction between the results of EAS simulations with EPOS-LHC and the experimental data: the measured maximum of the muon production profile appeared to be reached much higher in the atmosphere than predicted by the model, even when the heaviest possible primary cosmic ray nuclei were considered.
There are both similarities and differences concerning the relation of X max and X μ max to the properties of hadron-air interactions. Obviously, both characteristics are sensitive to the position X 0 of the primary particle interaction in the atmosphere, which depends on the respective inelastic cross section: fluctuations of X 0 shift the whole cascade upwards and downwards in the atmosphere and thus do so for X max and X μ max for a particular shower. However, in contrast to X max , maximum muon production depth is much less sensitive to the production of secondaries in the primary interaction. EAS muon content rather depends on the multi-step hadronic cascade in which the number of pions and kaons increases in an avalanche way until the probability of their decays becomes comparable to the one for interactions. This happens when energies for most of the hadrons approach the so-called critical energy: around 100 GeV for charged pions and an order of magnitude higher for charged kaons [25]. The position of the maximum of the muon production profile is actually close to this turning point.
From the above-discussion, it is clear that X μ max is very sensitive to the forward spectral shape of secondary mesons in pion-air collisions: producing in each cascade step a meson of a slightly higher energy would mean that a larger number of cascade branchings is required for reaching the critical energy, with the result that the maximum of the muon production profile will be observed deeper in the atmosphere. However, there exists another potential mechanism which may influence model predictions for X μ max . If the production of baryon-antibaryon pairs in pion-air interactions is copious enough, it will also influence the obtained X μ max . Indeed, (anti-)nucleons do not decay, hence, they continue to interact even when their energies fall below 100 GeV, producing additional generations of secondary hadrons in the cascade. Muons emerging from decays of secondary pions and kaons created in interactions of such low energy (anti-)nucleons contribute to the elongation of the muon production profile and give rise to larger values of X μ max . It is noteworthy that the respective effect is noticeable if (and only if) the yield of baryon-antibaryon pairs in pion-air interactions is comparable to the one of secondary pions.
To reveal the reasons for the large difference of the values of X μ max , predicted by QGSJET-II-04 and EPOS-LHC, shown respectively by red and blue solid lines in Fig. 2 (right), we are going to use the same "cocktail" model approach as in the previous Section. First, we apply QGSJET-II-04 to describe all the characteristics of the primary interaction while treating the rest of the hadron cascade using EPOS-LHC [yellow dashed line in Fig. 2 (right)]. As expected, the obtained X μ max deviates only slightly from the pure EPOS-LHC calculation (by less than 7 g/cm 2 ): the bulk of the difference between QGSJET-II-04 and EPOS-LHC predictions for the maximum of the muon production profile is due to secondary (mainly pion-air) interactions in the cascade. Further, we apply QGSJET-II-04 to describe both the primary interaction and the production of nucleons and antinucleons in all the secondary hadron-air collisions, while treating the rest with EPOS-LHC. The obtained X μ max shown in Fig. 2 (right) by the green dot-dashed line appears to be right in the middle between the QGSJET-II-04 and EPOS-LHC results. Finally, we have to check if the remaining half of the difference between the two models is indeed due to harder spectra of secondary mesons in EPOS-LHC for pion-air and kaon-air collisions. To this end, we use QGSJET-II-04 to describe the production of mesons for x F > 0.1 in all the secondary hadron-air interactions, apart from treating the primary interaction and the production of nucleons and antinucleons in secondary collisions. It is easy to see that the result shown by the black dotted line is indeed very close to the pure QGSJET-II-04 calculation.

Conclusions
We have demonstrated that both the maximum of the muon production profile and the average shower maximum depth are very sensitive to properties of very high energy pion-air and kaon-air interactions. In particular, a large part of the difference between the predictions of the QGSJET-II-04 and EPOS-LHC models for the two EAS characteristics is due to the copious production of baryon-antibaryon pairs in pion-air collisions in the latter model. As a matter of the fact, the mechanism is thus disfavored by the measurements of X μ max by the Pierre Auger experiment, which puts serious constraints on potential variations of model predictions for the shower maximum depth. Yet a large part of the model dependence for calculated X max is related to the properties of primary proton-air interactions. These will be constrained by further studies of proton-proton and protonnucleus collisions at the Large Hadron Collider, in particular, by measurements of forward particle production. Especially promising for testing the underlying physical mechanisms of the present interaction models seem to be combined studies by central and forward-looking particle detectors.