Cosmic ray energy spectra derived from KASCADE-Grande data using post-LHC hadronic interaction models

TL;DR

This paper presents updated cosmic ray energy spectra from KASCADE-Grande data, analyzing the impact of various post-LHC hadronic interaction models, including a new model, on observed spectral features like the knee and ankle.

Contribution

It introduces the use of the newly released EPOS-LHC-R model for cosmic ray spectrum analysis with KASCADE-Grande data, comparing results across multiple post-LHC models.

Findings

Identification of a knee-like structure at ~100 PeV in heavy component

Detection of an ankle-like feature in the light component

Comparison of spectral results across different hadronic models

Abstract

KASCADE-Grande was dedicated to measuring the energy spectrum and mass composition of cosmic rays in the energy range of 10 PeV to 1 EeV. We observed a knee-like structure in the heavy mass component at around 100 PeV and an ankle-like structure in the light component. In this contribution, we present updated energy spectra based on shower size measurements, using the post-LHC hadronic models QGSJet-II-04, EPOS-LHC, and SIBYLL 2.3d, including accounting for shower-to-shower fluctuations. In addition, the newly released EPOS-LHC-R model is tested for the first time with KASCADE-Grande. We will compare and discuss the results obtained using the different hadronic interaction models.

Tables

Table 1. Table 1: Spectral break positions and spectral indices obtained from broken power-law fits to the unfolded spectra of electron-poor (heavy) primaries. Note that the fit results of the EPOS-LHC-R model are based on the raw (not unfolded) spectrum.

electron-poor	log${}_{10}{}^{}(E_k/\mathrm{GeV})$	${\gamma }_1$	${\gamma }_2$	$\Delta \gamma$	${\chi }^2$/ndf
QGSJET-II-04	7.77 $\pm$ 0.02	2.69 $\pm$ 0.01	3.22 $\pm$ 0.04	0.53	3.92
EPOS-LHC	7.79 $\pm$ 0.04	2.78 $\pm$ 0.01	3.38 $\pm$ 0.06	0.60	4.03
SIBYLL 2.3d	7.78 $\pm$ 0.02	2.73 $\pm$ 0.01	3.24 $\pm$ 0.03	0.51	1.68
EPOS-LHC-R	7.79 $\pm$ 0.01	2.77 $\pm$ 0.01	3.23 $\pm$ 0.13	0.46	4.46

Full text

1 Introduction

To understand the transition between galactic and extragalactic cosmic rays, experiments investigate the energy spectrum and mass composition of cosmic rays in the range from PeV to EeV. KASCADE [1] and its extension KASCADE-Grande [2] were dedicated to studying extensive air showers in this energy range, particularly their primary energy, chemical composition, and arrival direction. Both experiments were located at the Karlsruhe Institute of Technology in Germany (49.1*∘N, 8.4∘*E, 110 m a.s.l.) and measured extensive air showers produced by cosmic rays in the PeV to EeV energy range. It is necessary to note that KASCADE’s multi-detector system enabled separate measurements of the number of electrons and muons in each shower, providing for studying mass composition.

Over more than 20 years of data collection, the experiments revealed valuable results. KASCADE observed a "knee" at around 3 PeV in the all-particle energy spectrum, mainly due to a drop in the number of light particles, such as protons and helium nuclei [3]. KASCADE-Grande extended the energy range and observed a "concave" feature in the spectrum above 1016 eV as well as an "iron-like" break around 1017 eV. This second break likely results from a steepening in the spectrum of heavy elements such as iron, which fits with the expectation that heavier particles are less deflected in magnetic fields [4, 5]. Additionally, KASCADE-Grande found an "ankle-like" feature around 100 PeV in the spectrum of light cosmic rays [6]. This feature might indicate a sign of an extragalactic component, which already dominates below $10^{18}$ eV.

Although the data collection completed in 2013, analysis continues. Among others, the data are now used to test and improve models of cosmic-ray interaction in the atmosphere. In this paper, we present updated energy spectra of light and heavy cosmic ray components, along with corrections for shower fluctuations based on an unfolding method. We also discuss recent studies on a new version of the hadronic interaction model EPOS.

2 All-particle and individual energy spectra

This work analyzes approximately 17 million events recorded by KASCADE-Grande over a total measuring time of 1865.62 days. After applying quality cuts including zenith angles of up to 40*∘, only events with full trigger and reconstruction efficiency were considered. Showers with more than 106* charged particles correspond to a primary energy of approximately 1016 eV.

Air shower simulations were performed using the CORSIKA program [7] with different hadronic interaction models. FLUKA was used for low-energy (E < 200 GeV) interactions, while QGSJET-II-04 [8], EPOS-LHC [9], SIBYLL 2.3d [10], and EPOS-LHC-R [11] were used for high-energy interactions. The simulations included five types of primary nuclei: proton, helium, carbon, silicon, and iron. Simulated showers cover a primary energy range from 1014 to 1018 eV and zenith angles from 0*∘* to 42*∘*. They were initially generated with a spectral index of -2 but subsequently reweighted to -3 for analysis.

The all-particle energy spectrum, along with individual spectra for heavy and light primary groups, is reconstructed based on the shower size ($N_{ch}$), which represents the number of charged particles. Measured data is divided into two mass groups, heavy and light, on an event-by-event basis. This separation is based on the correlation between the attenuation-corrected number of charged particles and the number of muons, referred to as the parameter $y_{CIC}={\rm log}_{10}(N_{\mu})^{CIC}_{CF}/{\rm log}_{10}(N_{ch})^{CIC}$. Hadronic interaction models predict different values for the shower size and muon numbers at the ground level. Consequently, the $y_{CIC}$ values are model-dependent. In general, air showers induced by heavy primaries, such as silicon and iron nuclei, contain fewer electrons at ground compared to those induced by light primaries, such as protons, helium and carbon nuclei. As a result, heavy primaries are referred to as the "electron-poor" group, while light primaries are defined as the "electron-rich" group.

Energy calibration functions are obtained by performing a linear fit to the correlation between the true primary energy and the shower size, both plotted on logarithmic scales. The slope of the fit is approximately 1. Further details on the energy calibration coefficients can be found in Ref. [13]. Using these energy calibration relations, the spectra for heavy and light primary mass group are reconstructed accordingly.

Since shower fluctuations in the energy determination are larger than the bin size of the reconstructed energy spectrum, an unfolding procedure is applied to correct for these effects. A brief summary of the method is provided below. More details can be found in Ref. [4]. Using Monte Carlo simulations, a response matrix $R_{ij}$ is constructed over the full energy interval of log${}_{10}(E/\rm GeV)$ from 6.0 to 9.0, where fluctuations can affect the measured energy spectrum. This matrix represents the conditional probability $P(E_{j}|E^{true}_{i})$, i.e., the probability that an event with true energy in bin log${}_{10}(E^{true}_{i})$ is reconstructed with energy log${}_{10}(E_{j})$. Based on this response matrix, a system of equations, $n^{exp}_{j}=\sum_{i=1}^{N}P(E_{j}|E^{true}_{i})n^{true}_{i}$, is established to relate the measured events $n^{exp}_{j}$ to the true energy distribution, $n^{true}_{i}$. This system is solved iteratively using the Bayesian unfolding algorithm [14]. The impact of this unfolding procedure on the resulting flux is estimated to be less than 10% for all energy bins and therefore does not significantly alter the overall shape of the spectra.

Figure 1 presents the unfolded energy spectra for heavy (Si+Fe) and light (H+He+CNO) primary cosmic rays, based on three different hadronic interaction models: QGSJET-II-04 (left), EPOS-LHC (middle), and SIBYLL 2.3d (right). These spectra have been corrected for shower-to-shower fluctuations using the unfolding procedure described above. The all-particle spectrum is obtained as the sum of the light and heavy mass components.

For all three models, the energy spectrum of the heavy component exhibits the knee-like feature around 1017 eV, characterized by a change in the spectral slope. The all-particle spectra shows smoother transition and less pronounced slope changes compared to the heavy components. In the light-component spectrum, a hardening, i.e., flattening of the spectrum, is observed above about 1017 eV, with the spectral slope changing gradually rather than sharply.

The analysis method and procedure introduce systematic uncertainties in the reconstructed energy spectrum and its extracted fit parameters. Some uncertainties may also be introduced by the reconstruction of the air showers. All these individual systematic contributions were considered and combined in quadrature to obtain the total systematic uncertainty in the flux, which is is of the order of about 13% for light and 10% for heavy primaries at the primary energy of $10^{17}$ eV. Further details on systematic uncertainties are discussed in Ref. [13].

Figure 2 compares the resulting all-particle energy spectra derived from different post-LHC interaction models with results from other experiments: ESA-TOP, IceCube/IceTop, Pierre Auger Observatory (PAO) and Telescope Array (TA). The all-particle energy spectra from KASCADE-Grande are obtained by summing the individual heavy and light component spectra. As shown in Fig. 2, KASCADE-Grande results on the all-particle spectrum show relatively weak dependence among post-LHC hadronic interaction models in use. However, the absolute flux measured by KASCADE-Grande is up to 10% lower than that measured by other experiments. This is possibly due to the observation altitude near sea level, as a different part of the shower evolution is taken into account in the interpretation by the hadronic interaction models. Furthermore, these differences are related to the absolute normalization of the energy scale by the various models. Despite differences in observation levels, analysis techniques and hadronic interaction models among the various experiments, there is good agreement in the energy range from PeV to EeV. At higher energies around $10^{18}$ eV, KASCADE-Grande results are statistically consistent with the result of the Pierre Auger Observatory.

3 New hadronic interaction model EPOS-LHC-R

A new version of the EPOS model, EPOS-LHC-R [11, 12], has recently been released. Improvements to the proton-proton cross-sections, multiplicity of proton-proton and proton-Air interactions, and nuclear fragmentation processes are applied. These changes impact both the shower maximum $X_{max}$ and the muon production in extensive air showers. An improved hadronization according to LHC data enhances the accuracy of muon production predictions, while the updated cross-section and nuclear fragmentation processes affect $X_{max}$. Compared to the previous version, EPOS-LHC, the EPOS-LHC-R model predicts an increased number of muons and a deeper $X_{max}$ by approximately 25 g/cm2, and thus even deeper than the predictions from the SIBYLL 2.3d model. This is expected to lead to an interpretation of a heavier composition for a given measured $X_{max}$. More details on EPOS-LHC-R are provided in Ref. [11, 12].

Full Monte-Carlo simulations, including the detector response of KASCADE-Grande, have been performed, using EPOS-LHC-R. Figure 3 (left) shows the correlation between the number of charged particles and the number of muons, $y_{CIC}$. The selection criteria separating heavy and light primaries is based on a fit to the average of silicon and CNO, yielding a separation value of 0.8496, indicated by the dashed line in Fig. 3. The right panel of Fig. 3 displays the energy calibration functions, described by: ${\rm log}_{10}(E_{true}/{\rm GeV})=a\cdot{\rm log}_{10}(N_{ch})+b$, with coefficients $a=0.897$, $b=1.705$ for heavy primaries, and $a=0.953$, $b=1.109$ for light primaries. The $y_{CIC}$ value for the new EPOS-LHC-R model is slightly lower than that of the previous EPOS model and closer to the values from the QGSJET-II-04 and SIBYLL 2.3d models. This indicates that EPOS-LHC-R predicts a slightly smaller relative muon content in air showers, aligning more closely with the predictions of QGSJET-II-04 and SIBYLL 2.3d models. The fact that the relative muon content is smaller means that for the observation with KASCADE-Grande the effect of the deeper $X_{max}$ (the electromagnetic component) outweighs the increased generated muon number in the model.

The reconstructed energy spectra using EPOS-LHC-R are shown in Figure 4 (left), where shower-to-shower fluctuations are not yet taken into account. A comparison with the EPOS-LHC model (unfolded) is presented in Fig. 4 (right). For heavy primaries, EPOS-LHC-R shows a slightly higher flux than EPOS-LHC, though spectral features remain similar. For light primaries, the EPOS-LHC-R model is less dominant compared to EPOS-LHC, resulting in a slightly lower all-particle spectrum for EPOS-LHC-R. Interestingly, both models show negligible difference above $10^{17}$ eV, while more notable differences appear in the lower energy range.

Figure 5 compares the energy spectra of heavy (Si+Fe) and light (H+He+CNO) mass groups for four different interaction models: QGSJET-II-04, EPOS-LHC, SIBYLL 2.3d, and EPOS-LHC-R. The first three models show unfolded spectra, whereas the EPOS-LHC-R spectrum is still raw. The heavy components spectra show similar knee-like structure across all models. QGSJET-II-04 and SIBYLL 2.3d yield comparable heavy spectra. In contrast, both EPOS-LHC and EPOS-LHC-R present slightly lower heavy fluxes, mainly due to differences in the $N_{ch}/N_{\mu}$ ratio. Light (electron-rich) spectra are consistently more abundant due to the separation boundary near the CNO mass group. All models exhibit similar hardening of light primaries with a smooth change of the slope around 1017 eV. The EPOS-LHC-R model shows the lowest flux of light primaries of all models, but this probably will slightly increase after applying shower-to-shower fluctuation corrections, potentially aligning more closely with QGSJET-II-04 and SIBYLL 2.3d.

To quantify the spectral features, a broken power-law fit is applied to the spectra of heavy particles. The break position and the spectral indices before and after the break are summarized in Table 1. Spectral breaks appear just below 1017 eV, with EPOS-LHC showing the most pronounced slope change.

4 Conclusion

Using the shower size measured by KASCADE-Grande and the $y_{CIC}$ technique, the energy spectra of heavy and light mass groups have been reconstructed based on the various post-LHC hadronic interaction models: QGSJET-II-04, EPOS-LHC, SIBYLL 2.3d, and the newly released EPOS-LHC-R. The EPOS-LHC-R model, tested with KASCADE-Grande data for the first time, show spectral behavior consistent with other models and confirms previous findings. To correct bin-to-bin migrations caused by shower-to-shower fluctuations, response matrices to each interaction model were constructed and applied to unfold the reconstructed spectra. The impact of this unfolding procedure on the resulting flux is estimated to be less than 10% for all energy bins and thus does not significantly change the spectral shapes. All prominent spectral features observed in previous measurements are confirmed: observation of a heavy knee at around $10^{17}$ eV accompanied by a flattening of the light component at the same energy. This flattening may be an indication of an extragalactic component, which appears to be dominant already below the energy of $10^{17}$ eV for all post-LHC models.

Additionally, the muon content in cosmic-ray air showers was investigated in the primary energy from $10^{16}$ and $10^{18}$ eV for three zenith angle intervals. While a slight zenith angle dependence is observed, the measured muon contents by KASCADE-Grande shows a reasonable agreement with the MC expectations. Further details can be found in Ref. [15].

Finally, the KASCADE Cosmic-ray Data Centre (KCDC) [16] provides a public accessible, web-based platform (https://kcdc.iap.kit.edu), where scientific data of the completed KASCADE and KASCADE-Grande (and more) experiments are made available for the astroparticle community and interested public (for latest news on KCDC see also Ref. [17]).