Hadronic Interactions and Air Showers: Where Do We Stand?

The interpretation of EAS measurements strongly depends on detailed air shower simulations. CORSIKA is one of the most commonly used air shower Monte Carlo programs. The main source of uncertainty in the prediction of shower observables for different primary particles and energies is currently dominated by differences between hadronic interaction models even after recent updates taking into account the first LHC data. As a matter of fact the model predictions converged but at the same time more precise air shower and LHC measurements introduced new constraints. Last year a new generation of hadronic interaction models was released in CORSIKA. Sibyll 2.3c and DPMJETIII.17-1 are now available with improved descriptions of particle production and in particular the production of charmed particles. The impact of these hadronic interaction models on air shower predictions are presented here and compared to the first generation of post-LHC models, EPOS LHC and QGSJETII-04. The performance of the new models on standard air shower observables is derived. Due to the various approaches in the physics treatment, there are still large differences in the model predictions but this can already be partially resolved by comparison with the latest LHC data.


Introduction
Knowing the elemental composition of cosmic ray particles arriving at Earth is of crucial importance to understanding their production and propagation. Unfortunately, cosmic rays can only be measured 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. With the operation of modern large-scale experiments, the reliability of air shower simulations has become the source of the largest systematic uncertainty in the interpretation of cosmic ray data [2]. 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 [3].
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 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. 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 experi- * e-mail: tanguy.pierog@kit.edu ments and is also the region of the highest rapidity-density of produced particles. The first data have been compared to cosmic ray models in [4]. 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 collider experiment -such as in diffractive interactions -are 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 particles and their energy spectra which have a direct impact on air shower development [5]. These latest measurements are not yet taken into account in the available hadronic interactions models, but are very important to understanding the open issues in these models and for their future developments.
At the same time, a new generation of hybrid cosmic ray detectors such as the Pierre Auger Observatory [6] (surface and fluorescence detectors), the IceCube/IceTop experiments [7,8] (low energy particles at the surface and high energy muons deep underground) or the KASCADE/ KASCADE Grande experiment [9,10] (particles 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  Figure 1. Elementary parton-parton scattering: the hard scattering in the middle is preceded by parton emissions attached to remnants. The remnants are an important source of particle production even at intermediate energies (∼100 GeV cms). models is able to give a fully consistent picture of the different observables within a given experiment [11].
In this paper, we compare the latest hadronic model predictions after LHC data and their consequences on air shower observables. In the second section, we give a general description of the models and in the third section we compare their results for the observables important for the air shower development. Using detailed Monte Carlo simulations done with CONEX [12], the new predictions for X max , X µ max and for the number of muons are finally presented.

Hadronic Interaction Models
There are several hadronic interaction models commonly used to simulate air showers. Nowadays there are four high energy models which were updated to take into account LHC data at 7 TeV: QGSJETII-03 [13] changed into QGSJETII-04 [14], EPOS 1.99 [15] replaced by EPOS LHC (v3400) [16], and more recently Sibyll 2.1 [17] updated to Sibyll 2.3c [18] in CORSIKA v7.6300 [19]. The old DPMJET2.55 [20] has been updated to a new version DPMJETIII.17-1 [21,22] whose preliminary results are presented in this paper. There is no major change in these models but, in addition to some technical improvements, some parameters were changed to reproduce TOTEM [23] cross sections.
They all are based on the simple parton model associated with the Gribov-Regge multiple scattering approach which can be seen as a multiple exchange of "parton ladders" between a projectile and a target, see Fig. 1. But they differ in their philosophy.

DPMJET model
DPMJETIII.17-1 is a minimum bias Monte Carlo hadronic generator used for both heavy ion interactions and cosmic ray air shower simulations. It is a full extension of the parton model to nuclear interactions but it does not include any final state interactions due to high density (no collective hadronization). As a consequence it is better suited to study jet production than soft particle production. For instance, charm particle production will follow perturbative quantum chromodynamic (pQCD) calculations.

EPOS model
EPOS LHC is a minimum bias hadronic generator used for both heavy ion interactions and cosmic ray air shower simulations. The goal of this model is to describe soft particle production (p t <∼ 5 GeV/c) for any system and energy in very fine details (rare particles, all possible data). To be able to describe all types of heavy ion data, nuclear effects related to Cronin transverse momentum broadening, parton saturation, and screening have been introduced into EPOS [15]. Furthermore, high density effects leading to collective behavior in heavy ion collisions are also taken into account [24].

QGSJETII model
QGSJETII-04 model is a minimum bias nuclear 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 limits the data sets to which it can be compared.

Sibyll model
Like QGSJETII, the Sibyll model is a minimum bias hadronic interaction model optimized for air shower simulations but with a different approach. It is a minimal extension of the parton model to run with light nuclei using the semi-superposition model. It produces more types of particles compared to QGSJETII but the comparison to nuclear collision data is really limited. To some extent it is a simplified version of the DPMJETIII model.
Compared to Sibyll 2.1 the new version has an improved production of baryon-antibaryon pairs, in particular from the mini-jet (hard) particle production, and a phenomenological model for the production of charm particles which is important for the production of high energy muons and neutrinos.
Even if the four 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 [25], 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 [26]. We will show that the extrapolation to high energy is not really the issue anymore.

Inelastic cross section
As shown in [25], the inelastic nuclear cross section is very important for the development of air showers and in particular for the depth of the 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 via the optical theorem in all hadronic interaction models. Therefore it is very well described by all the models up to LHC energies, where data exist. As shown in [27], thanks to the measurements at the LHC even the extrapolations up to the highest energy are now very similar.
However, plotting the prediction of these models for the proton-air and pion-air inelastic cross-sections as shown in Fig. 2 left-hand side, one notices that significant differences appear which will have direct consequences on air shower development. In all the figures DPMJETIII.17-1 is represented by a dotted (indigo) line, EPOS LHC by a full (blue) line, QGSJETII-04 by a dashed (red) line and Sibyll 2.3c by a dash-dotted (green) line. Not only do the evolutions diverge at high energy, but for Sibyll 2.3c and DPMJETIII.17-1 the relative behavior of the proton and pion-air cross-section is different from the other models (faster increase of the pion-air cross-section to reach the proton-air one).

Multiplicity
According to [25], the multiplicity plays a similar kind of role as the inelastic cross section, but with a weaker dependency (log). On the other hand the predictions from the models have larger differences for the multiplicity compared to the cross section.
As shown in [27], 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 [28] 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 a function of 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. 2 right-hand side). This effect is small compared to the pre-LHC era [26] but can change the elongation rate of the air shower maximum development. Here again Sibyll 2.3c and DPMJETIII.17-1 have 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 a 10% difference.
So, for both cross section and multiplicity, when the models are constrained by LHC data up to 7 TeV, the extrapolations 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 [25] defined as the largest energy fraction carried by a secondary particle (the leading particle).
In additions to the parton ladder, there is another source of particle production: the two off-shell remnants, see Fig. 1. This is directly related to the elasticity since  the leading particle is usually produced by the projectile remnant.
All models have some remnant from the projectile and target but the simplest scheme which does not allow more than one quark exchange with the central ladder is used except for EPOS which has more options. The simplest approach allows, for instance, the production of leading ρ 0 in π−p/A interactions while EPOS using a fully generalized scheme allows any flavor in the remnant. This is in fact needed for the consistency of the model (no difference between first scatterings and next ones) and to reproduce multi-strange baryon production at low energy [29].
The model predictions are shown in Fig. 3 for p-p, πair and p-air (as inelasticity=1-elasticity) as a function of center of mass energy. Sibyll 2.3c has the largest elasticity which 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 p-p 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 to the elasticity is the diffractive dissociation. Diffraction is a special case of interaction where no central ladder is produced and there is only some momentum exchange between excited remnants. There are only technical differences to treat diffraction in the models. 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 [30]. At the LHC various measurements related to diffraction are now available [31,32]. 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 uncertain-ties in air shower simulations at a level of 10 g/cm 2 [33]. Nevertheless the difference between models seems to be even larger as illustrated in Fig. 4 left-hand side. 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.3c to produce a rapidity gap around 2 to 4 is a direct consequence of the too narrow pseudorapidity distribution and implies a large elasticity.

Baryon and resonance production
Another important observable for EAS is the number of muons reaching the ground. It has been shown in [35] 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 [36].
Recent measurements by NA61 [37] show that the ρ 0 production in π-C interactions seems to be underestimated by a relatively large amount (from 20% to 100% for all models but Sibyll 2.3c which was tuned to the data) potentially leading to a large increase of muon production [38]. Furthermore in [39] 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 section 4.3). And indeed from [34,40] it can be concluded that the excess of protons seen in [35,41] is not due to newly produce baryons but is due to some baryons 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 implication on muon production in air showers as shown in the next section.

Depth of shower maximum
As shown in Fig. 5, the mean depth of shower maximum, X max , for proton and iron induced showers simulated with CONEX is different for DPMJETIII.17-1, EPOS LHC, QGSJETII-04 and Sibyll 2.3c as a direct consequence of the differences shown in section 3. However, the elongation rate (the slope of the 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 [26]. 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 section 3 it is likely that on the one hand Sibyll 2.3c predicts too large values of the X max since the multiplicity is already too low and the elasticity too high at the LHC. On the other hand QGSJETII-04 is at the lower edge of possible predictions compatible with LHC data since the multiplicity is at the higher limit, the cross-section is low and the rapidity gap (diffraction) is also low. 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 [42].
In fact, a further study using the fluctuations of X max around the mean can be used to test model consistency. Indeed both 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 [25]) than the mean value, it can be checked that the composition corresponding to a given X max is consistent with the observed fluctuations. In [42] the Pierre Auger Collaboration shows that while it is possible to describe the observed data with EPOS LHC, QGSJETII-04  Refs. to the data can be found in [1] and [42]. is in tension with data at a 1 sigma level ( X max too shallow by ∼15g/cm 2 ) confirming that this model is the lower edge of the allowed X max region.

Muons at ground level
Concerning the number of muons at the ground (for 40 o inclined showers at a height of 1500 m), the difference between DPMJETIII.17-1, EPOS LHC, QGSJETII-04 and Sibyll 2.3c is relatively small. We can see in Fig. 6 lefthand side 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 showers differs from the model predictions by 1.33±0.13±0.09 [11] in the best case. This is a 2 sigma effect and in the case of  Figure 7. Mean logarithmic mass from X µ max (red) and from X max (black) as a function of energy from [46] using the QGSJETII-04 (left-hand side) or the EPOS LHC (right-hand side) hadronic interaction model. inclined showers the effect is less than 2 sigmas too [43]. Taking into account the ρ 0 measurement as explained in section 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 muons in the atmosphere [44] and for the muons at the ground in general [45]. As a consequence, one of the most sensitive measurements 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 [46]. As shown in Fig. 7 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 [47], EPOS 1.99 was giving a mean composition 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 [32]. We can see in Fig. 4 left-hand side that EPOS LHC gives reasonable results while QGSJETII-04 is too low. Sibyll 2.3c, which overestimate the fraction of large rapidity gaps (high elasticity), predicts deep MPD as well (and probably incompatible with the data since it is very close to EPOS).
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 types of projectile. While the change of diffraction and thus of elasticity in proton interactions has very little impact on X µ max , it appears that the MPD is 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 a large shift in X µ max . This is confirmed by the results of DPMJETIII.17-1. The rapidity gap distribution in p-p shown in Fig. 4 lefthand side is larger than the measured one at the LHC, but on the other hand the diffractive cross-section for pion interactions is very low in comparison to other models (and EPJ Web of Conferences 208, 02002 (2019) https://doi.org/10.1051/epjconf/201920802002 ISVHECRI 2018 data). In Fig. 4 right-hand side (the elasticity of the pionair interactions is about 10% lower than other models at low energy). As a result X µ max is lower than QGSJETII-04 by about 30 g/cm 2 as shown in Fig. 5 right-hand side while X max is in the same range as the other models.
Hence we can say that the dependence of the MPD on the pion elasticity is so strong that the 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 [39]. From the EAS development we can thus say that the 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 [39] 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). As explained in section 3.4, new accelerator data confirm that the forward baryon production should be reduced in EPOS, leading to shallower X µ max .

Summary
In [39] 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% is 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 section 3 we have shown that for the first interaction the uncertainty is not in the basic p-p interaction anymore, very well constrained by LHC data, but by the nuclear effects which cannot be tested properly with current model and data combinations (data with heavy ion only at high energy and only EPOS LHC can treat heavy ion collisions properly). These nuclear effects being important for both the air target [48] and in the case of a heavier primary, are the main source of the systematic shift in X max but which is limited to about ±20 g/cm 2 around EPOS LHC predictions (or the band defined by the predictions of Sibyll 2.3c and QGSJETII-04). It is very unlikely that a model compatible with accelerator data up to LHC energies could predict a X max outside this range. 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 [49].
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. Upcoming studies of diffraction at the LHC, including those with a nuclear target [5,50], will reduce the model uncertainty for the first interaction to 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 [39] for which cumulative effects due to the hadronic cascade are observed. This should give qualitative input to improve the models which can then be quantitatively tested against past and future NA61 measurements for instance.
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 concern now mainly 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.