Air showers, hadronic models, and muon production

We report on a study about some characteristics of muon production during the development of extended air showers initiated by ultra-high-energy cosmic rays. Using simulations with the recent new version of the AIRES air shower simulation system, we analyze and discuss on the observed discrepancies between experimental measurements and simulated data.


Introduction
The determination of the composition of the ultraenergetic cosmic rays is one of the most relevant current challenges of Cosmic Ray experimental physics. The hypothesis widely accepted at present is that the flux of ultra-energetic cosmic rays incident on the terrestrial atmosphere is constituted essentially by a mixture of different atomic nuclei. Under this hypothesis, the determination of the primary composition is equivalent to the determination of the masses of the incident nuclei. In the experiments that can measure the longitudinal development of particle showers initiated after the incidence of the mentioned ultra-energetic cosmic rays, the atmospheric depth of the maximum of these showers, X max , is a robust observable directly correlated with the mass of the primary particle. Such correlation can be established by comparison with numerical simulations that prove to depend on the algorithms used to model the hadronic collisions that take place during shower development [1].
Measurements of X max and its fluctuations made by the Pierre Auger Observatory [2,3] seem to indicate that the composition of ultra-energetic cosmic rays varies with their energy. The results published in [2,3] clearly show that the values of X max and its fluctuations are approximately compatible with those obtained by numerical simulations for light primaries (A ∼ 1) for primary energies around 2 × 10 18 eV, while for larger energies a progressive shift towards a heavier composition is observed. This characteristic of variable composition with energy is analyzed in several works. In particular, in [4] an adjustment of the total flux of cosmic rays as a function of primary energy is made, simultaneously with X max and its fluctuations, using a conveniently parameterized procedure, obtaining a combined setting that describes the total flux of ultraenergetic cosmic rays that arrive at the Earth as a mixture of elements (H, He, N, and Si) obtaining for each one of them the corresponding partial contribution in function of the primary energy. This work also includes a discussion * e-mail: sciutto@fisica.unlp.edu.ar on extensions of the fitted contibutions for energies below 5 × 10 18 eV that adds a fifth contribution assigned to Fe nuclei (see discussion in reference [4]).
The muonic content of showers initiated by cosmic rays is another important observable directly related to the composition of cosmic rays. However, the estimation of the primary composition of particle showers using experimental measurements has not yet been satisfactorily completed due to the discrepancies that are observed when comparing the experimental data with the results of computational simulations [5][6][7]. An example of this is the estimation of the number of muons in inclined showers published by the Pierre Auger Observatory [5], where it can be clearly observed that the simulations produce estimates that are significantly lower than the experimental results. It is important to mention that similar disagreements between experimental measurements and simulations are also obtained with alternative analyses such as those reported in references [6,7].
In this work we address the issue of the observed discrepancies between experiment and simulations of the muon production in air showers by presenting an alternative discussion of the estimates of the number of muons at ground level, under conditions similar to those of reference [5], using simulations carried out with the new version of our AIRES program, recently publicly launched [8,9].

Simulations and results
We recall that AIRES [8,9] is a system for the realistic simulation of particle showers in the Earth's atmosphere, also containing a series of tools for the analysis of the corresponding results. In its current version, 18.09.00, recently launched [9], a large series of new features have been incorporated. Among these new features, we find important to mention: (1) Links with the hadronic interaction models EPOS [10], QGSJET [11], and SIBYLL [12], in both their pre-and post-LHC versions. 100 Auger data AIRES+Epos LHC AIRES+Sibyll 2.3c simulation of decays of unstable particles, including up to many tens of relevant branches in the case of particles with complex decay patterns. (4) Extended sets of data that are recorded in the listings of particles reaching the ground or of longitudinal survey, including, for example, identity and energy of the mother particle, identity and energy of the projectile of the last hadronic event in the history of the particle, etc. A complete description of the AIRES system, including access to the corresponding user manual, is placed at [9].
As a first step of our analysis with AIRES, we have carried out numerous simulations in order to reproduce the most relevant results that have already been published. In the case of X max and its fluctuations, we have verified that there is an excellent agreement between the simulations with the current version of AIRES and those published in references [2,3], which were made with CONEX [13], with links to the already mentioned EPOS, QGSJET, and SIBYLL hadronic models in their pre and post-LHC versions. For reasons of brevity, we omit presenting in this work the data in graphic form.
A similar check can be made in the case of the number of muons at ground level. Following reference [5], we run simulations of air showers with an inclination from the vertical of 67 degrees, and for each set of showers we evaluate the average number of muons at ground level that have a total energy larger than 300 MeV, N MC µ . Then we calculate R MC µ can be directly compared with the experimental measurements, as explained in reference [5].
In figure 1 we present our results for R MC µ /E prim corresponding to a first series of simulations performed with AIRES linked to EPOS-LHC or SIBYLL 2.3c, and proton or iron nuclei primaries. Experimental data of the Pierre Auger Observatory published in [5] is also displayed for comparison purposes. We verified that the simulations performed with AIRES + EPOS-LHC (blue lines) are virtually coincident with those corresponding to simulations with CORSIKA [14] + EPOS-LHC plotted in the similar figure 4 of reference [5]. The green lines in figure 1, correspond to simulations with AIRES + SIBYLL 2.3c, and correspond to larger values of R MC µ than the CORSIKA + SIBYLL 2.1 values represented in figure 4 of the aforementioned reference [5], as expected for this pre-LHC version of SIBYLL which produces a smaller number of muons during the development of the showers (for a detailed discussion on this issue see reference [1]).
In both simulations shown in figure 1 it is evident that the number of muons predicted by them is significantly lower than the corresponding experimental data. As already mentioned, this lack of agreement between experimental data and simulations is presently one of the most relevant enigmas of ultra-energetic cosmic rays, and is receiving most attention in current research. It should be noted that this "excess" of muons in the experimental data is also reported in other analysis, either from the Pierre Auger Observatory [6], or the Telescope Array Collaboration [7].
The curves of figure 1 corresponding to the simulations for proton and iron clearly show the well-known characteristic that for showers of a given primary energy, the pro- duction of muons is larger the larger the mass of the primary. Connecting this with the already mentioned result that measurements of the depth of the shower maximum, X max , and its fluctuations provide favorable evidence for a mixed primary composition of cosmic rays that reach the Earth's atmosphere, with proportions that depend on energy, we consider it worthwhile to investigate what the impact of this hypothesis would be in the case of muon production.
To this end, we take into account the adjustment published in reference [4], where the flux of cosmic rays with primary energies greater than 5×10 18 eV is the sum of four components, namely, protons (H 1 ), He 4 , N 14 , and Si 28 , whose relative abundance is estimated with an elaborate procedure, described in detail in reference [4], which takes into account the whole process of creation, acceleration and propagation of the primaries through the intergalactic space, and allows the simultaneous adjustment of the energy spectrum, X max , and its fluctuations σX max . This combined adjustment can be extended to energies slightly below 5 × 10 18 eV by adding an additional contribution to the cosmic ray flux, of galactic origin and including a new heavy fraction that is considered to be composed of iron (Fe 56 ) nuclei (see reference [4], figure 14 and the related discussion). From these results published by the Pierre Auger Observatory, we calculated the normalized abundances of the five components, protons (H), He, N, Si and Fe, in the range of energies relevant to our study of muon production in inclined showers, which are represented in figure 2.
With these primary fractions plotted in figure 2 we have carried out numerous simulations with AIRES. First, we have verified the correct adjustment of X max and its fluctuations with the results published by the Pierre Auger Observatory [2,3], which are also used for the combined adjustment of the reference [4]. Our results are shown in figure 3 (curves in red line), where the very good degree of agreement between the simulations and the experimental values of both X max and its fluctuations can be appreciated.
Finally we have carried out simulations to determine R MC µ under the hypothesis of a combined composition, with the abundances as represented in figure 2. The results are shown in figure 4. The red squares with solid lines represent the results with AIRES + EPOS-LHC, this being the same hadronic model as used in the combined fit of reference [4]. The most outstanding aspect of this curve is that it is approximately horizontal, that is, with a similar dependence on energy as the experimental data points. Thus, the difference between simulations and experimental data virtually reduces to a constant that accounts for the muon deficit that arises from comparing the experimental data with the simulations. In effect, if we add a constant, properly adjusted, to the values of R MC µ /E prim of this red curve, we can raise it and make it coincide with all the experimental points within the statistical errors range (red dashed curve figure 4). Figure 4 also includes the results of the simulations with AIRES + SIBYLL 2.3c (green circles with solid line) which show a behavior similar to the case of EPOS-LHC, with slightly higher muon production. As in the case of simulations with EPOS-LHC, the curve corresponding to SIBYLL 2.3c can be moved towards higher values of R MC µ /E prim with the addition of a global constant, and get a one sigma adjustment to all experimental data points.

Final remarks
A new study of the production of muons in showers initiated by ultra-energetic cosmic rays has been presented, using simulations with the recently released version of AIRES [9].
To carry out this study it has been necessary to verify the proper simulation of usual observables, such as X max and its fluctuations, for example, in a series of well-known cases, allowing to ensure the correctness of the new interfaces with hadronic models that are currently included in AIRES.
In the case of the number of muons in inclined showers, it has been possible to corroborate that the simulations with AIRES linked to the EPOS-LHC, QGSJET-II-04, and SIBYLL 2.3c models, predict a significantly lower number of muons compared to the experimental data, in the same direction as the results reported in the references [5][6][7].
In our analysis of the number of muons, we have also made the comparison of the experimental data with simulations taking into account a combined primary composition [4], with abundances dependent on energy. In this case, it is noteworthy that although a muon deficit persists in the simulated data with respect to the experimental ones, this difference can be virtually reduced to a constant, independent of the primary energy, shift in R µ /E which added to the results of the simulations allows to locate the displaced curve in agreement with the experimental data within error bars. We understand this as a new indication   Figure 3. X max and σX max versus primary energy. The experimental points reproduce the Auger Observatory data published in reference [3]. The error bars (horizontal square brackets) correspond to statistical (systematic) uncertainties. The red curves (labelled CFitC) correspond to simulations with AIRES + EPOS-LHC using an energy dependent combined composition (see text).  supporting the hypothesis of the combined composition of ultra-energetic cosmic rays. The results of this analysis call for a more complete study in order to understand the origin of this displacement, and to improve the hadronic models to obtain predictions compatible with the experimental results. We are currently investigating in this direction.