Testing of almost all the hadronic interaction models by comparing calculated muon energy spectrum with data

. Uncertainties of the model energy spectra of the most energetic secondary charged mesons are discussed. Computer simulations of the partial energy spectra of the atmospheric vertical muons induced by primary cosmic particles with various ﬁxed energies in terms of hadronic interactions models had been carried out with the help of the CORSIKA package. These partial spectra have been convolved with the contemporary spectra of the primary cosmic particles in the energy range 0.1-10 000 TeV. Results of simulations are compared with the contemporary data of the atmospheric vertical muon ﬂux. Comparison shows that all models underes-timate the production of secondary charged π ± -mesons (and K ± -mesons) by a factor of ∼ 1.4 ÷ 2 at the highest energies. This underestimation induces a more rapid development of extensive air showers in the atmosphere and results in uncertainties in estimates of energy and composition of the primary cosmic particles.


Introduction
Extensive air showers (EAS) are the only tool to understand the origin and composition of cosmic rays, their possible sources and the transport of cosmic particles in various magnetic fields on their way to the Earth at very high energies. All features of the energy spectrum, arrival directions and composition of the primary cosmic particles should be determined through an analysis of the EAS data. These data as signals in the surface and underground detectors are usually interpreted in terms of various models of hadronic interactions [1][2][3][4][5][6][7][8][9]. It happened, that such interpretation leads to some inconsistency. As an example, the energy of showers calculated in terms of the QGSJET II-03 [3] model with the help of the surface detectors signals at the Telescope Array [10] happened to be 1.27 times lager than this energy estimated with help of the fluorescence light. To ensure that results of such interpretation are as accurate as possible these models should be thoroughly tested. Usually these models are tested with the help of accelerator data at small values (∼0) of the pseudorapidity, η, where most of secondary particles (mainly mesons) are produced [11][12][13]. However, calculations have shown that the maximal energy flow carried by secondary particles occurs at much larger values (∼8-10) of the pseudorapidity η [14]. Let us also note that the longitudinal development of EAS de-pends strongly on the rate of the projectile particle energy fragmentation. The atmospheric muon flux also depends strongly on this production of the highest energy mesons. So, it is of primary importance to verify a production of the most energetic mesons simulated in terms of various models. This verification may be carried out by comparing model predictions of these muon fluxes with data. We select the classical experiments L3+Cosmic [15], MACRO [16], LVD [17] and IceCube [18] and elaborate the smooth approximation of these muon data in the energy interval of 10 2 − 10 5 GeV. The model predictions of muon flux have been estimated as follows. Showers induced by primary protons and helium nuclei with different fixed energies have been simulated with the help of the CORSIKA package [19] and the muon partial energy spectrum in each individual shower have been calculated. Results of these simulations for every type of primary particles multiplied by intensities of these particles should be integrated on the energy of particles. Thus, we also need some expressions for the energy spectra of various primary particles. Inspired by modern results and a new precise cosmic ray data base [20], we suggested new approximations of cosmic ray energy spectra for primary protons and helium nuclei. Indeed, there are results of many measurements of the fluxes of the primary cosmic nuclei (e.g.AMS-02 [21], PAMELA [22], ATIC-2 [23], CREAM [24], ARGO-YBJ [25], ARGO-YBJ & FWCTA [26], KASCADE [27], KASCADE-Grande [28], Tunka [29], IceCube [30], Telescope Array [31]). Besides there are some calculations of spectra of the primary proton and helium nuclei in SNR [32]. We will use these ap-  Energy spectra of the primary protons. Solid linesuggested approximation. All particle spectra depends on the energy per particle.
proximations for the energy spectra of the primary protons and helium nuclei to estimate convolutions with the partial muon spectra. Thus, with the help of any of the interaction models [1][2][3][4][5][6][7][8][9], the CORSIKA package and approximations of data on fluxes of the primary cosmic nuclei [21][22][23][24][25][26][27][28][29][30][31] one can predict the energy spectra of atmospheric vertical high energy muons at sea level. These predictions can be compared with data observed by the L3+Cosmic, MACRO and LVD collaborations at energies above 100 GeV. Finally, some conclusion can be drawn about the validity of various models. In fact, some low energy models with the FLUKA package [33] have been tested in such a way. We are sorry that some our results of models testing in [34][35][36] are not correct. We do apologize for our mistake in input data for the atmosphere.

Method
To estimate the energy spectra D(E µ ) of atmospheric vertical muons in the energy range of 10 2 − 10 5 GeV we need to know the energy spectra dI p /dE and dI He /dE of the primary protons and helium nuclei within the energy interval 10 2 − 10 7 GeV and the partial energy spectra S p (E µ , E) and S He (E µ , E) of the vertical muons in EAS induced by the primary protons and helium nuclei with various fixed energies, E. Simulations of these partial muon spectra have been carried out in terms of the QGSJET01, QGSJET II-03, QGSJET II-04, DPM-JET 2.55, VENUS 4.12, EPOS LHC, SIBYLL 2.1 and SIBYLL 2.3 hadronic interaction models in the same energy range of 10 2 −10 7 GeV. The smooth approximation of the atmospheric muon data observed by the L3+Cosmic,  Functions S p (E µ , E) and S He (E µ , E) are the partial differential energy spectra of muons in showers induced by primary protons and helium nuclei with fixed values of energy E. These spectra were calculated for 24 different fixed values of energy E of the primary protons. The energy distributions of muons induced by primary helium nuclei S He (E µ , E) have also been calculated and compared with simulations based on the hypothesis of superposition [37]. As direct results coincide with simulations in terms of the hypothesis of superposition, we will use this hypothesis.
To take into account a possible change of primary spectrum above the "knee" at energies above E 1 10 6 GeV for the primary protons and helium nuclei we have used an additional exponential multiplier. The values of parameters for these new approximations are listed in Table 1. For primary protons β=5.5417 and γ=0.024. For primary helium nuclei β=4.4074 and γ=0.027.
The approximation parameters: R 1 , R 2 , α, s, β and γ are dimensionless. Parameters K 1 and K 2 are dimensional as [1/(GeV · m 2 · s · sr)]. The notation E is the kinetic energy per nucleon in GeV.
E ∈ 10 6 ÷ 10 7 GeV (1) The CORSIKA 7.4 package (CORSIKA 6.9 in the case of DPMJET 2.55 and QGSJET II-03 models and COR-SIKA 7.56 in the case of SIBYLL 2.3 model) have been used to simulate the second important ingredients -the partial energy spectra S p (E µ , E) and S He (E µ , E) of vertical muons in showers induced by the primary protons and helium nuclei with various fixed energies, E, in terms of the eight models in the energy range 10 2 − 10 5 GeV with statistics 10 6 events for the most energetic muons.
The results of these calculations in the energy range 10 2 −10 7 GeV were interpolated for 100 values of energies E with equal intervals in decimal logarithmic scale. The energy interval 10 2 − 10 5 GeV of muons was divided into 60 equal bins also in decimal logarithmic scale. So, the width of the bin was equal to h = lg(E µ,(i+1) /E µ,i ) = 0.05. Let us note that average muon energies for the 1 st , 21 st and 41 st bins we will use later are equal to 1.059 · 10 2 , 1.059 · 10 3 and 1.059 · 10 4 GeV respectively. In fact simulations for helium nuclei have been carried out only for energies 10 4 and 10 6 GeV to test the hypothesis of superposition. As results of simulations for the primary nuclei showed a good agreement with this hypothesis we have used this hypothesis to estimate the flux of the nucleons from the primary helium nuclei.
The energy spectra D p (E µ ) and D He (E µ ) of muons for primary protons and helium nuclei are calculated as integrals of products of functions S p (E µ , E) and S He (E µ , E) with corresponding intensities dI p /dE and dI He /dE of the primary protons, on energy E of primary nucleons.
The resulting energy spectrum of atmospheric muons is the sum of these energy spectra of muons produced by primary protons and helium nuclei.

Results of simulations
The partial energy spectra S p (E µ , E) of the atmospheric vertical muons simulated for various fixed energies E of the primary protons in terms of the EPOS LHC model are shown in figure 3. It is seen that statistics of ∼ 10 6 at the higher energy end (tail) of the spectra is not enough. Table 2 displays the total number of muons with energies above 10 2 and 10 3 GeV in showers induced by primary protons with energies 10 5 and 10 6 GeV estimated in   terms of the eight models in our simulations and in [38]. A very reasonable agreement is seen. The next figure 4 demonstrates a comparison of the partial muon energy spectra S p (E µ , E) calculated in terms of the eight models for a fixed energy E = 10 4 GeV of the primary protons. The results for the SIBYLL 2.1 model are the highest, and the QGSJET01 model values are the lowest at muon energy E µ 10 3 GeV. The result of the rest of models are in between these two models.
It is important to estimate energy intervals of the primary protons which contribute into various bins of the muon energy spectrum. But first some dependence of muon numbers inside the partial bins on energy E should be illustrated. Figure 5 demonstrates the distributions of the primary proton energy E for the bins of muon energy spectra for the SIBYLL 2.3 model. Let us remember that average energies of muons for these bins are equal to 1.059 · 10 2 , 1.059·10 3 and 1.059·10 4 GeV accordingly. It is possible to note that nearly two orders of energy E are of importance   for any fixed muon energy E µ . The maximal contributions occur at energies ∼ 5 · 10 2 GeV (to the 1 st bin), ∼ 5 · 10 3 GeV (to the 21 st bin) and ∼ 5 · 10 4 GeV (to the 41 st bin). The final results of the muon energy spectra D(E µ ) calculated in terms of the eight hadronic interaction models and data [15][16][17][18] are shown in figure 6. The ratios of MC simulation to data are shown in figure 7.
It should be noted that at energies above ∼ 100 GeV both the simulated spectra and data are steepened. It is because the decay constant B for the charged mesons is equal to B ∼ 100 GeV and the probability of decay for charged mesons is decreasing at higher energies.
It should be noted that contributions to the generation of muons into the 21 st bin with an average energy of E µ 10 3 GeV are most considerable from primary protons with energy E = 10 4 GeV as shown in figure 5. So, various models may be effectively compared at energies at which their contributions into the muon spectrum is most considerable. Such comparison of relative contributions of the eight models to the generation of muons into the 21 st bin are shown in figure 4 and corresponds to the final result in the same energy region shown in figure 7.
It is seen that calculated spectra are ∼ 2 times below data in the case of the QGSJET01 model and ∼ 1.4   times below data for the SIBYLL 2.1 model. The results of DPMJET 2.55 model is ∼ 1.86 below data at energies E µ 10 2 GeV and ∼ 1.36 below data at energies E µ 10 4 GeV. The results of QGSJET01 model is ∼ 2 below data at energies E µ 10 2 GeV and ∼ 1.7 below data at energies E µ 10 4 GeV. The results of SIBYLL 2.1 model is ∼ 1.5 below data at energies E µ 10 2 GeV and ∼ 1.35 below data at energies E µ 10 4 GeV. The result of the rest of models are in between these limits and slightly different from each other by ∼10-15% . The main conclusion is quite clear. All considered models demonstrate the valuable deficit of muons [39][40][41][42][43].
In this paper the muon fluxes intensities for the DPM-JET 2.55, QGSJET II-03, SIBYLL 2.3 and VENUS 4.12 models differ by ∼ 2 ÷ 15% with an increase in energy from ∼ 10 2 GeV to ∼ 10 4 GeV compared with the work of Ref. [42,43] due to slightly different modified approximations of the primary cosmic ray spectrum.

Conclusion
Muons which contribute much to the muon energy spectra are produced in decays of the most energetic π ± -mesons and K ± -mesons generated in first interactions of the primary particles with nuclei in the atmosphere. As calculated vertical muon energy spectra in the case of the  3) hadronic interaction models are ∼ 1.4 ÷ 2 times below data we can conclude that production of the most energetic π ± -mesons and K ± -mesons in these models is considerably suppressed. This suppression may induce smaller values of signals in the surface scintillation detectors and will result in larger values of the calculated energy estimates. So, the coefficient 1.27 used by the TA collaboration [10] to decrease the energy estimates of showers calculated on the basis of signals in the surface detectors may be understood as a result of this suppression. The increased intensity of the primary particle flux observed at the Yakutsk array at super high energies [44] may also be a result of smaller values of calculated signals in surface scintillation detectors.