Testing of the VENUS 4.12, DPMJET 2.55, QGSJET II-03 and SIBYLL 2.3 hadronic interaction models via help of the atmospheric vertical muon spectra

Uncertainties of the model energy spectra of the most energetic secondary π± -mesons (and K±-mesons) are discussed. Computer simulations of the partial energy spectra of the atmospheric vertical muons induced by primary cosmic particles with various fixed energies in interval 102−107 GeV in terms of the VENUS 4.12, DPMJET 2.55, QGSJET II-03 and SIBYLL 2.3 models had been carried out with help of CORSIKA package. These partial spectra should be convoluted with the contemporary spectra of the primary cosmic particles. Results of simulations are compared with the contemporary atmospheric vertical muon flux data. Comparison shows that all models underestimate production of secondary π±-mesons (and K±-mesons) by factor of ∼ 2 at the most high energies. This underestimation induces 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
The extensive air showers (EAS) are the only a tool to understand the origin and composition of cosmic rays, their possible sources and a 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 some 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, energy of showers calculated in terms of the QGSJET II-03 [3] model with help of the surface detectors signals at Telescope Array [10] happened to be 1.27 times lager than this energy estimated with help of the fluorescence light. To be sure 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 the 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 depends strongly on the rate of the projectile particle energy fragmentation. The atmospheric muon flux also depends strongly on this production of the most high energy mesons. So, it is of the 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 this muon fluxes with data. We select the classical experiments L3+Cosmic [15], MACRO [16] and LVD [17] 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 estimate as follows. Showers induced by the primary protons and helium nuclei with different fixed energies have been simulated with help of the COR-SIKA package [18] and the muon partial energy spectrum in each individual shower have been calculated. Results of these simulations for every type of the primary particles multiplied by intensities of these particles should be integrated on energy of particles. Thus, we need also some expressions for the energy spectra of various primary particles. Inspired by modern results and new precision cosmic rays data base [19], we suggested new approximations of cosmic ray energy spectra for primary protons and helium nuclei. Indeed, there are results of many measurements (e.g.AMS-02 [20], PAMELA [21], ATIC-2 [22], CREAM [23], ARGO-YBJ [24], ARGO-YBJ & FWCTA [25], KAS-CADE [26], KASCADE-Grande [27], Tunka [28], IceCube [29], Telescope Array [30]) of the fluxes of the primary cosmic nuclei. Besides there are some calculations of spectra of the primary proton and helium nuclei in SNR [31]. We will use these approximations for energy spectra of the primary protons and helium nuclei to estimate convolutions with the partial muon spectra. Thus, with the help of any interaction models [1][2][3][4][5][6][7][8][9], the package CORSIKA and approximations of data on fluxes of the primary cosmic nuclei [20][21][22][23][24][25][26][27][28][29][30] one can predict the energy spectra of atmospheric vertical high energy muons at sea level. These predictions can be compared with data observed by collaborations L3+Cosmic, MACRO and LVD at energies above 100 GeV. Finally, some conclusion can be drawn about validity of various models.
In fact, some low energy models with the package FLUKA [32] have been tested in such a way. We are sorry that some our results of models testing in [33][34][35] are not correct. We do apologize for our mistake in input data for the atmosphere.
In this paper models QGSJET II-03, DPMJET 2.55 [8], VENUS 4.12 [9] and SIBYLL 2.3 [5] have been tested. A comparison of muon data observed in [15][16][17] with results of simulations allows to draw a conclusion about the most energetic mesons production described by these models.

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 the various fixed energies E. Simulations of these partial muon spectra have been carried out in terms of the VENUS 4.12, DPMJET 2.55, QGSJET II-03 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 collaborations L3+Cosmic, MACRO and LVD had been used for comparison with results of simulations.  -IceCube (All particles); -KASCADE (All particles); -KASCADE-Grande (All particles). All particles spectra depends on energy per particle.
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 been also calculated and compared with simulations based on the hypothesis of superposition [36]. As direct results coincides with simulations in terms of the hypothesis of superposition, we will use this hypothesis.
The energy spectra of the primary particles are important ingredients of simulations. As the energy per nucleon is of importance only the energy spectra of the primary protons and helium nuclei should be taken into account. We had used approximations (1) for (dI p /dE) and (dI He /dE) suggested in [37]. Figures 1 and 2 demonstrate how these approximations fit data [20][21][22][23][24][25][26][27][28][29][30].
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 had used additional exponential multiplier. The values of parameters for these new approximations are listed in Table 1.
The results of these calculations in the energy range of 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 accordingly. 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 had 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 the 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.
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 VENUS 4.12, DPMJET 2.55, QGSJET II-03 and SIBYLL 2.3 models are shown in figure 3. It is seen that statistics of ∼ 10 6 at the higher energy end of the spectra is not enough.   Table 3 displays the total number of muons with energies above 10 2 and 10 3 GeV in showers induced by the primary protons with energies 10 5 and 10 6 GeV estimated in terms of the VENUS 4.12, DPMJET 2.55, QGSJET II-03 and SIBYLL 2.3 models in our simulations and in [38]. The 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 SIBYLL 2.  It is of importance 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.
These dependence on energy E for the 1-st, 21-st and 41-st bins of muon energy spectra S p (E µ , E) are displayed in figure 5 for the QGSJET II-03 and SIBYLL 2.3 model. The good agreement is also seen. Figure 6 demonstrates distributions of the primary proton energy E for the bins of muon energy spectra for the SIBYLL 2.3 model. Let us remind 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 energy E µ of muons. 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 VENUS 4.12, DPM-JET 2.55, QGSJET II-03 and SIBYLL 2.3 models and data [15][16][17] are shown in figure 7 (left panel) and ratios of MC simulation to data -in the same figure (right panel).
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 a probability of decay for charged mesons is decreasing at higher energies.
It is seen that calculated spectra are ∼ 2 times below data in case of the QGSJET II-03 model and ∼ 1.7 times below data for the SIBYLL 2.3 model. The result of the rest of models are in between of these limits. The main conclusion is quite clear. All considered models demonstrate the valuable deficit of muons.

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 case of the VENUS 4.  JET 2.55, QGSJET II-03 and SIBYLL 2.3 models are ∼ 1.7 ÷ 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 base 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 [39] may be also a result of smaller values of calculated signals in surface scintillation detectors.