Tests of hadronic interaction models with the KASCADE-Grande muon data

The KASCADE-Grande experiment is an air-shower ground-based observatory designed to study in detail the energy spectrum and composition of primary cosmic rays in the region of 10 − 10 eV. These analyses are based on precise measurements of the charged, electron and muon numbers of the cosmic ray air-showers performed through different detector systems which come into play simultaneously in KASCADE-Grande during the data acquisition. Due to the quality of the data and the number of air-shower observables at disposal through the experiment the collected data proves to be also useful to test hadronic interaction models used for air-shower simulations. In this contribution, predictions of the QGSJET II-2, SIBYLL 2.1 and EPOS 1.99 hadronic models are confronted with the KASCADEGrande muon data. Besides, the influence of these models on the all-particle energy spectrum derived from the muon size is also investigated.


Introduction
Cosmic-ray simulations are a useful tool to reconstruct the properties of primary cosmic rays at high energies, when studying the extensive air showers that the primary radiation induce in the Earth's atmosphere. Hadronic interaction models are an important part of these simulations. They are also the main source of uncertainties in cosmic ray studies, which arises due to the lack of an accurate description of the physical phenomena ocurring in the kinematical region of small transverse momenta, the most important one for EAS development. The difficulty comes from the fact that in the very forward region (small p t ), QCD can not be applied perturbatively and, even worse, almost no data is available. This situation demands the a e-mail: arteaga@ifm.umich.mx employement of phenomenological models tuned up with accelerator data at low energies to describe hadronic interactions. At the high-energy regime models are extrapolated to be used in EAS simulations [1]. The distinct phenomenological approaches and parameterizations used in the models result in very important differences in relevant EAS quantities at high-energies, such as the inelastic cross section, the inelasticity of hadron-hadron interactions and the total number of charged particles, which can be measured with dedicated air-shower observatories. Combining precise measurements of several EAS parameters, airshower data can be used to test and improve hadronic interaction models. In this way, EAS facilities can serve also as particle physics laboratories to explore kinematic and energy regions not available for present collider experiments.  In general, EAS studies show that predictions of hadronic models employed in EAS simulations present discrepancies from the observed data. For example, analyses performed with the multi-detector set-up KASCADE have shown that models such as QGSJET 98 and 01, EPOS 1.66, SIBYLL 1.6 and 2.1, VENUS, NEXUS and DPM-JET are not able to describe simultaneously all the data on the hadron, electron and muon contents of air showers (see [2] and references therein) in the primary cosmic ray energy interval 10 14 − 10 16 eV. These studies in KASCADE have been extended up to 10 18 eV with its extension, the KASCADE-Grande observatory [3]. Some of these analyses are focused on the muon content. Muons, as a penetrating component of air-showers, are an important tool to get insight into the hadronic interactions happening during the development of the EAS. With KASCADE-Grande, in [4] it has been shown that predictions with QGSJET II about the muon production height (H µ ) distributions for EAS with zenith angles below θ < 18 • show a discrepancy compared to the measured data (specifically, a disagreement at large muon production heights, H µ > 3.5 km). Investigations of muon pseudorapidities with the Muon Tracking Detector of KASCADE-Grande [5] have also shown the deficiencies of QGSJET II. In addition, in [6] the lateral distribution function (LDF) of muons was studied with KASCADE-Grande. In general, it was observed that the slope of the muon LDF is not in agreement with results of the Monte Carlo simulations (QGSJET II and EPOS 1.99). In this work, predictions on the muon content of EAS from several current hadronic interaction models are confronted with measurements from KASCADE-Grande in the energy interval 10 16 − 10 18 eV. In particular, the muon attenuation length is extracted and studied as a tool to investigate the dependence of the muon content on the zenith angle in the atmosphere. The study is performed in an energy-independent way by using the so-called constant intensity cut method.

The KASCADE-Grande observatory
In KASCADE-Grande, measurements of the total muon number in EAS (N µ , number of muons with energy greater than 230 MeV) are performed with an array of 192×3.2 m 2 shielded scintillator detectors belonging to the KASCADE part of the experiment. The procedure implies the calculation of the number of penetrating particles traversing each KASCADE muon station from the energy deposits obtained by sampling the shower front at distances larger than 40 m from the core [3]. The estimation is done by using a lateral energy correction function (LECF), derived from simulations based on CORSIKA. The number of muons in the EAS is estimated from a maximum likelihood estimation, assuming that locally the detected muons fluctuate according to a Poisson distribution [3]: where n i is the number of muons measured at a core distance r i in the i−th KASCADE muon station with sensitive area A i . Here, θ is the zenith angle of incidence of the EAS and f (r) is a lateral distribution function for muons, which has a Lagutin-Raikin form [7], The coefficients p 1 = −0.69, p 2 = −2.39, p 3 = −1.0, and r 0 = 320 m were obtained from fits to COR-SIKA/QGSJET01 simulations (protons and iron nuclei with energies of 10 16 and 10 17 eV). As it can be seen, the shape of the lateral distribution function is fixed and is not fitted event by event [3]. That comes from the fact that the number of muon data points is too low to produce a stable fit. On the other hand, arrival times and charged particle densities, employed for estimations of the EAS arrival direction, charged particle content and core position, are measured with an extension called the Grande array [3], which is composed by 37 × 10 m 2 plastic scintillator detectors distributed on a surface of 0.5 km 2 .

MC simulations
Three high-energy hadronic interaction models were tested in this study: QGSJET II-02 [8], EPOS 1.99 [9] and SIBYLL 2.11 [10]. EAS development in the atmosphere was simulated with CORSIKA [11] and the response of the KASCADE-Grande detectors, with a GEANT 3.21 based code. The low-energy hadronic interactions were treated with FLUKA [12]. MC simulations for single primaries: H, He, C, Si and Fe were performed. Additionally, a set with mixed composition was created (five primaries of equal abundance). Data was generated following a E −2 distribution. Proper weights were added to produce data sets with spectral indexes: γ = −2.8, −3.0 and −3.2.
Selection cuts were applied to both experimental and MC data. They were chosen according to MC studies to avoid as much as possible the influence of systematic uncertainties in the measurements of the EAS parameters. Data sets were composed of events with more than 11 triggered stations in Grande, shower cores inside a central area of 1.52 × 10 5 m 2 and arrival directions confined to the zenith angle interval of ∆θ = 0 • − 40 • . These events were registered during stable periods of data acquisition and passed successfully the standard reconstruction procedure of KASCADE-Grande. Additionally, only showers with log 10 N µ > 5.1 were considered for this work. Both the experimental and simulated data were analyzed and reconstructed with the same algorithms. With the above quality cuts, the effective time of observation with KASCADE-Grande was equivalent to 1424 days. The threshold for full efficiency was found at log 10 N µ ≈ 5.4.
The models tested predict different muon contents for EAS at a fixed energy. In general, EPOS 1.99 gives more muons than QGSJET II-2 and SIBYLL 2.1, while the latter produces less muons for the same energy than QGSJET II-2. For example, comparing simulated data for air showers in the zenith angle interval θ = 20 • − 26 • , at different energies (see Fig 1) it is found that the mean number of muons from EPOS 1.99 is higher by 14% than QGSJET II-2, and 21% higher than SIBYLL 2.1. The evolution of  Figure 4. Integral muon number spectra for two different zenith angular ranges for measured and simulated data in KASCADE-Grande. In these plots, N µ has been corrected for systematic effects according to QGSJET II. Simulated fluxes are normalized according to experimental data using the first angular bin. N µ in the atmosphere is also different in each case. This point will be analysed later in section 5.

Muon systematics
Systematic errors on the muon number were studied in detail with MC simulations. From these analyses muon correction functions were built as functions of the arrival direction, core position and muon content of the EAS for each hadronic interaction model, assuming a mixed composition and γ = −3. Experimental data was corrected with the aforementioned functions to study also the effect of the hadronic interaction model in the interpretation of the muon data. Each MC muon data set was treated with the correction function of the respective hadronic interaction model.
In Fig. 2 the mean value of the muon correction function for different hadronic interaction models is plotted as a function of the uncorrected N µ . In general, after correction the systematic error on the muon number above threshold is found to be almost independent of the corrected muon size, N ′ µ , and smaller than 6%.

Analysis and results
To test the hadronic interaction models with the KASCADE-Grande muon data, predictions on the evolution of the muon content with the arrival zenith angle of the EAS were confronted with observations. As a first step, the muon fluxes were reconstructed for five different zenith angle intervals, each with the same acceptance (c.f. Fig 3). Then, the integral muon fluxes, J(N µ ), are calculated for each ∆θ range. If MC fluxes are normalized in such a way that vertical showers agree with the experimental values around log N µ = 5.1 − 5.6, one observes that MC values for more inclined showers deviate from the measured fluxes. The differences increase for higher zenith angles, as can be seen from Fig. 4. In general, it is observed that the values of the experimental number of 07002-p.3 muons in EAS are larger than the predictions by MC simulations and that the difference between the measured and predicted N µ increases for more inclined showers. A more detailed comparison between the expected and observed muon measurements can be done by calculating the muon attenuation length, Λ µ . This quantity is extracted by applying the Constant Intensity Cut (CIC) method to the J(N µ ) data as described in reference [13] but using a global fit to the attenuation curves, log 10 N µ (θ), with the known formula where X 0 = 1023 g/cm 2 is the average atmospheric depth for vertical showers and N µ is a normalization parameter to be determined for each attenuation curve. The results for Λ µ are presented in Table 1 for simulated data and Table 2 for experimental measurements (corrected with appropriate correction functions). Comparison of the data and simulations is shown in Fig. 5. Discrepancies between the experimental values and the simulation results can be observed for the studied models. The differences do not disappear when modifying the primary composition or spectral index. As a consequence, the predicted evolution of the muon component with the zenith angle, N µ (θ) (see equation 2) shows also a disagreement with the observations as shown in Fig. 6, where the evolution of the mean muon number from MC data (γ = −3) is compared with that from experimental data (corrected for systematic effects) in the framework of the QGSJET II hadronic interaction model. Data has been normalized at θ = 22 • , where the maximum of the experimental zenith angular distribution is found. Similar results are found when the EPOS and SIBYLL hadronic interaction models are employed. The percentage of deviation for N µ (θ) between experiment and simulations depends on the zenith angle of normalization. The maximum deviation for inclined showers is found when normalizing at θ = 0 • . Such a differences are plotted in Fig. 7 assuming a mixed composition scenario with γ = −3 for MC data, for the models QGSJET II-2, EPOS 1.99 and SIBYLL 2.1. Formula 2 was employed to calculate the curves of Fig. 7 using the values for Λ µ from Tables 1 and 2. Several factors, which are model dependent, may come into play in the observed differences. Beginning from predicted muon correction function (see graphs in Fig. 2), up to the description of the production, evolution and fluctuations of the shower, and systematic errors from the assumed shape of the muon lateral distribution function and the lateral energy correction function for muons in KAS-  CADE. Therefore, one should be cautious when extracting conclusions from these differences. Systematic studies are underway to investigate these individual possibilities.

Conclusions
The predictions on the muon content of EAS from three hadronic interaction models: QGSJET II-2, EPOS 1.99 and SIBYLL 2.1 were tested in this work. In particular, the muon attenuation lengths of the penenetrating component was calculated from MC simulations and compared with the values extracted from KASCADE-Grande observations. It was found that the above hadronic interaction models do not describe this aspect of the muon component in air-showers. More tests are under way. These tests include a more detailed analyses of the effect of the muon systematics on Λ µ . 07002-p.4   Figure 7. Relative deviation of measured data from model predictions for the dependence of the muon content with the zenith angle in the framework of three hadronic interaction models. The normalization angle is chosen at θ = 0 • .