The cosmic ray energy spectrum in the range 1016–1018 eV measured by KASCADE-Grande

The KASCADE-Grande experiment, located at Campus North of the Karlsruhe Institute of Technology (Ger- many) is a multi-component extensive air-shower experiment devoted to the study of cosmic rays and their interactions at primary energies 10 14 − 10 18 eV. One of the main goals of the experiment is the measurement of the all particle en- ergy spectrum in the 10 16 − 10 18 eV range, i.e. extending the range accessible by KASCADE alone. The Grande detec- tor samples the charged component (Nch) of the air shower while the original KASCADE array provides in addition a measurement of the muon component (Nµ). The combined information of Nch and Nµ is used to estimate the energy on an event-by-event basis and to derive the all-particle energy spectrum. Since the calibration of the observables in terms of the primary energy depends on Monte Carlo simulations, three different methods with partially different sources of un-


Introduction
The study of the energy spectrum and of the chemical composition of cosmic rays are fundamental tools to understand origin, acceleration and propagation of cosmic rays. The range between 10 16 eV and 10 18 eV can only be studied by means of air-shower arrays which sample specific components of the cascade of particles produced by the primary cosmic ray in the atmosphere. This range is quite important from astrophysical point of view because it is expected that in this energy range the transition between galactic and extragalactic origin of cosmic rays will occur. The results obtained at lower energies by KASCADE (Antoni et al., 2005) and EAS-TOP (Aglietta et al., 2004) as well as by other experiments suggest that the knee in the primary energy spectrum around 3−4×10 15 eV is due to the break in the spectra of light elements. Moreover, KASCADE results seem to indicate that a rigidity dependent mechanism is responsible for such knees. Therefore, a knee of the heaviest components would be expected in the range of 10 16 eV to 10 18 eV. Various theories with different assumptions try to explain the rather smooth behavior of the cosmic ray energy spectrum between 10 16 eV and 10 18 eV. The idea of the "dip model" (Berezinsky et al., 2006) is that the ankle at (4 − 10) · 10 18 marks the energy loss by electron pair production of extragalactic protons, hence, at these energies the composition is already proton dominant and the transition between galactic and extra-galactic components already happens at energies lower than 10 18 eV. On the other hand, Hillas (Hillas, 2005) proposed the idea of a "component B" of cosmic rays of galactic origin which shows a rigidity dependence, like the "component A" responsible for the well-known knee, but shifted in energy. In this case, the transition between galactic and extragalactic components will occur at the ankle, and in the region between 10 16 eV and 10 18 eV a mixed composition is expected to prevail.
In order to discriminate between the different models, a very precise measurement of the possible structures of the energy spectrum is needed. The aim of the KASCADE-Grande experiment is to provide high-quality data in exactly this energy range by using the sizes of the charged-particle and muon components.

The apparatus
The KASCADE-Grande experiment (Apel et al., 2010) is a multi-detector setup consisting of the KASCADE experiment (Antoni et al., 2003), the trigger array Piccolo and the scintillator detector array Grande. Additionally, KASCADE-Grande includes an array of digital read-out antennas, LOPES (Falcke et al., 2005;Schröder et al., 2011), to study the radio emission in air showers at E > 10 17 eV. Most important for the analysis presented here are the two scintillator arrays: KASCADE and Grande. The KASCADE array comprises 252 scintillator detector stations structured in 16 clusters. The detector stations house two separate detectors for the electromagnetic (unshielded liquid scintillators) and muonic components (shielded plastic scintillators) where muon detectors are housed in 12 clusters (or 192 stations), only. This enables to reconstruct the lateral distributions of muons and electrons separately on an event-by-event basis. The Grande array is formed by 37 stations of plastic scintillator detectors, 10 m 2 each (divided into 16 individual scintillators) spread on a 0.5 km 2 surface, with an average grid size of 137 m. All the 16 scintillators per station are viewed by a high gain photomultiplier (for timing and low particle density measurements) and the four central ones are additionally viewed by a low gain one (for high particle densities). The dynamic range of the detectors is 0.3−8000 particles / 10 m 2 . Grande is arranged in 18 hexagonal clusters formed by 6 external detectors and a central one. Grande and KASCADE arrays are both triggered when all the 7/7 stations in a hexagon have fired (rate ∼ 0.5 Hz). Full efficiency for proton and iron primaries is reached at logN ch > 5.8.
KASCADE-Grande provides the unique opportunity of evaluating the reconstruction accuracies of the Grande array by a direct comparison with an independent experiment. For a subsample of events collected by the Grande array it is possible to compare on an event-by-event basis the two independent reconstructions of KASCADE and Grande. By means of such a comparison the Grande reconstruction accuracies are found to be for the shower size: systematic uncertainty ≤ 5%, statistical inaccuracy ≤ 15%; for arrival direction: σ ≈ 0.8 • ; for the core position: σ ≈ 6 m (Apel et al., 2010). All of them are in good accordance with the resolutions obtained from simulations.

The analysis
A sample of Monte Carlo data was simulated including the full air shower development in the atmosphere, the response of the detector and its electronics as well as their uncertainties. In this way, the reconstructed parameters from simulation are obtained in the same way as for real data. The EAS events were generated with an isotropic distribution with spectral index γ = −3 and were simulated with CORSIKA (Heck et al., 1998) and the hadronic Monte Carlo generators FLUKA (Fassò et al., 2005) and QGSJet II (Ostapchenko, 2006a,b). Sets of simulated events were produced in the energy range from 10 15 eV to 10 18 eV with high statistics and for five different representative mass groups: H, He, C, Si and Fe (≈ 353.000 events per primary). Few events up to 3 · 10 18 eV were also generated in order to cross-check the reconstruction behavior at the highest energies.
Grande stations are used to provide core position and angle-of-incidence as well as the total number of charged particles in the shower, by means of a maximum likelihood procedure comparing the measured number of particles with the one expected from a modified NKG lateral distribution function (Apel et al., 2006) of charged particles in the EAS.
The total number of muons is calculated using the core position determined by the Grande array and the muon densities measured by the KASCADE muon array detectors. The total number of muons N µ in the shower disk (above the energy threshold of 230 MeV) is derived from a maximum likelihood estimation assuming the locally detected muons to fluctuate according to a Poisson distribution. Here, we have chosen a lateral distribution function based on the one proposed by Lagutin and Raikin (Lagutin and Raikin, 2001). Since the muon densities are very low, except for the highest energy showers, stable fits on shower-by-shower basis are only obtained if the shape of the lateral distribution function is kept constant and only the muon number N µ is taken as a free parameter. By means of Monte Carlo simulations a resolution of the total muon number in the showers of ≈ 25% for energies above 10 16 eV is expected. The reconstruction procedures and obtained accuracies of KASCADE-Grande observables are described in detail in Apel et al. (2010).
For the reconstructed events, we restricted ourselves to events with zenith angles lower than 40 • due to increasing muon number and shower size uncertainties coming from the reconstruction procedure. Additionally, only air showers with cores located in a central area on KASCADE-Grande, represented in Fig. 1, were selected. With this cut on the fiducial area, border effects are discarded and the under-and overestimations of the muon number for events close to and far away from the center of the KASCADE array are reduced. This is motivated by the fact that the muon detector is asymmetrically located with respect to the array. Therefore, in case of small and large distances, the total number of muons obtained by fitting the local muon densities with the lateral distribution function in a limited range of distances, would introduce some bias for very small and large core distances. All of these cuts were applied also to the Monte Carlo simulations to study the effects and to optimize the cuts. Full efficiency for triggering and reconstruction of air-showers is reached at primary energy of ≈ 10 16 eV, slightly depending on the cuts needed for the reconstruction of the different observables (Apel et al., 2010).
The analysis presented here is finally based on 1173 days of data and the cuts on the sensitive central area and zenith angle correspond to a total acceptance of A = 1.976 · 10 9 cm 2 sr, and an exposure of N = 2.003 · 10 17 cm 2 s sr.
In KASCADE-Grande, the all-particle energy spectrum is derived by means of different techniques: using only one component, such as the size of the charged (Kang et al., 2009) or muon (Arteaga et al., 2009) components, a combination of the two (Bertaina et al., 2009), or the particle density at 500 m, S(500) (Toma et al., 2009). In this way, by combining and comparing the results of the different methods, which are sensitive to different parameters of the showers, and have intrinsically partially different systematics or sources of uncertainty it is possible to: (i) cross-check the measurements by different sub-detectors; (ii) cross-check the reconstruction procedures; (iii) cross-check the influence of systematic uncertainties; (iv) test the sensitivity of the observables to the elemental composition; (v) test the validity of hadronic interaction models underlying the simulations.
The conversion between the observed quantities (charged particle and muon sizes as well as the combination of both) of the EAS to the energy of the primary particle requires the assumption of a specific hadronic interaction model, whose general suitability has to be verified beforehand. In this work, the energy estimations are based on the CORSIKA-QGSjetII model, motivated by the fact that this model reproduces fairly well the distributions of the ratio of the muon and electron sizes measured by KASCADE-Grande as a function of both, the electron size and the atmospheric depth (up to 40 • zenith angle) (Cantoni et al., 2009). In the following, results will be presented on the basis of three among these methods: N ch or N µ alone and their combination. In the N ch -method and N µmethod, the reconstructed charged particle or muon shower size per individual event is corrected for attenuation in the atmosphere by the constant intensity cut method and calibrated by Monte Carlo simulations under the assumption of a dependence E 0 ∝ N α ch ch (E 0 ∝ N α µ µ ) and a particular primary composition (Kang et al., 2009;Arteaga et al., 2009). The energy conversion relation (calibration function) is obtained by selecting events of zenith angle ranges of 17 • ≤ θ ≤ 24 • , i.e. around the reference angle (20 • for N ch and 22 • for N µ distributions), to reduce systematic effects. The relation of the primary energy as a function of the number of charged particles (muons) is obtained assuming a linear dependence in logarithmic scale logE 0 = a N ch,µ + b N ch,µ · logN ch,µ . The same procedure is performed for H, He, C, Si and Fe primaries as well as for a uniform mixed composition to examine the dependence of the calibration on the assumed primary particle type. The energy resolution is estimated from the difference between simulated and derived energy. In case of N ch (N µ ) the energy resolution for proton and iron are about 31% (25%) and 15% (12%), respectively, over the whole energy ranges, where the reconstruction uncertainty for N ch gives the largest contribution. In addition, the systematic uncertainties on the reconstructed energy spectrum are investigated considering various possible contributions: a) the associated errors to the fit of the attenuation curve and of the calibration function, b) uncertainties due to a possible misdescription of the attenuation in the Monte Carlo simulation are estimated by using two different reference zenith angles (10 and 30 degrees), c) shower fluctuations, d) uncertainty in the spectral index. All these individual systematic contributions were considered to be uncorrelated, and combined thus in quadrature to obtain the total systematic uncertainty, where the composition dependence was not taken into account. The systematic uncertainty (i.e. sum in quadrature of all terms discussed above except the energy resolution) on the flux for proton and iron are for N ch (N µ ) 21% (7%) and 17% (13%), respectively, at energies of 10 17 eV. Finally, to reduce the effects of the fluctuations on the determined energy spectrum, an unfolding procedure is applied which uses a response matrix built using Monte Carlo simulations. This matrix represents the probability that an event being reconstructed with an energy inside the bin log 10 (E i ) has a true energy in the channel log 10 (E T rue j ). By unfolding such matrix it is possible to deconvolve the effect of shower fluctuations.
In the N ch − N µ method, the information provided by the two observables is combined. By help of Monte Carlo simulations a formula is obtained to calculate the primary energy per individual shower on the basis of N ch and N µ . The formula takes into account the mass sensitivity in order to minimize the composition dependence. Advantage of this technique is that the correlation between the two particle components is taken into account on an event-by-event basis. The attenuation is corrected for by deriving the formula for different zenith angle intervals independently and combining the energy spectrum afterwards (Bertaina et al., 2009). The energy assignment is defined as E = f (N ch ,k) (see Eq. (1)), where N ch is the size of the charged particle component and the parameter k is defined through the ratio of the sizes of the N ch and muon (N µ ) components: k = g(N ch ,N µ ) (see Eq. (2)). The main aim of the k variable is to take into account the average differences in the N ch /N µ ratio among different primaries with same N ch and the shower to shower fluctuations for events of the same primary mass: log 10 (E[GeV ]) = [a p + (a F e − a p ) · k]· log 10 (N ch ) (1) k = log 10 (N ch /N µ ) − log 10 (N ch /N µ ) p log 10 (N ch /N µ ) F e − log 10 (N ch /N µ ) p (2) log 10 (N ch /N µ ) p,F e = c p,F e · log 10 (N ch ) + d p,F e .
The k parameter is, by definition of Eq. (2), a number centered around 0 for a typical proton shower and 1 for a typical iron shower. Also in this case a response matrix is applied afterwards to take into account the migration of low energy events to higher energies. After applying the procedure to the 5 angular bins independently, a combined spectrum is obtained by using all the events of the 5 bins calibrated with their own calibration function. Also in this case systematic uncertainties due to the following reasons were considered: a) different intensity in the five angular bins, b) capability of reproducing an, a priori assumed, single primary spectrum with slope γ = −3, c) uncertainty on the index of the energy spectrum used to derive the energy relation E(N ch ) and the response matrix, d) uncertainty on the core location. The total uncertainty is around 15% and it is quite stable in the entire energy range. Regarding the energy resolution, it varies between 26% at the threshold and 20% at the highest energies due to the decrease of the intrinsic shower fluctuations.

Results
In Fig. 2 the resulting spectra are compiled, where the flux is multiplied by a factor of E 3.0 . Owing to the different procedures, the results for the first two methods are shown under proton and iron assumption, respectively, only, whereas for the other method the resulting all-particle energy spectrum is displayed with a band for systematic uncertainties. Taking into account the systematic uncertainties, there is a fair agreement between the all-particle energy spectra of the different applications.
Of particular interest is the fact that by using N ch , the iron assumption predicts a higher flux than the proton assumption, whereas using N µ the opposite is the case. That means that the "true" spectrum has to be a solution inside the range spanned by the two methods. If there is only the possibility of applying one method, then there is a large variance in possible solutions (everything within the range spanned by proton and iron line, not even parallel to these lines). However, more detailed investigations have shown, that a structure in the spectrum or a sudden change in composition would be retained in the resulting spectrum, even if the calibration is performed with an individual primary, only. All variations (within reasonable physics scenarios, i.e. everything between pure proton and pure iron, including other mass groups) of the composition would lead to a spectrum within the given error band. Interestingly, over the whole energy range there is only little room for a solution satisfying both ranges, spanned by N ch and N µ , and this solution has to be of relatively heavy composition -in the framework of the QGSJet-II hadronic interaction model. The narrower range for a solution provided by the N µ method compared to N ch confirms the finding of KASCADE that at sea-level the number of mostly low-energy muons N µ is a very good and composition insensitive energy estimator. The results of the composition independent N ch -N µ method lie inside the area spanned by the composition dependent methods, again confirming that the reconstructions of the shower size are under control. In addition, this finding expresses a cross-check of the intrinsic consistency of the hadronic interaction model. Future investigations will perform a similar analysis to test the validity of the interaction model EPOS (Werner et al., 2006) to describe the air shower development.
Despite the overall very smooth power law behavior of the resulting all-particle spectrum, there are some smaller structures observed, which does not allow to describe the spectrum with a single slope index. Figure 3 shows the resulting all-particle energy spectrum multiplied with a factor that the middle part of the spectrum is flat. There is a clear evidence that just above 10 16 eV the spectrum shows a "concave" behavior, which is significant with respect to the systematic and statistical uncertainties. Such a hardening of the spectrum is expected when a pure rigidity dependence of the galactic cosmic rays is assumed. In this case the gap between light primaries (proton, helium) and the CNO (Z = 6 − 12) group in their knee positions requires a hardening of the spectrum (De Donato and Medina, 2009). But there are also other astrophysical scenarios possible for a concave behavior of the cosmic ray spectrum. In general, a transition from one source population to another one should result in a hardening of the spectrum. In this aspect, the KASCADE-Grande result could be a first experimental hint to the "component B" of galactic cosmic rays as proposed by Hillas (2005). Another feature visible in the spectrum is a small break at around 10 17 eV. A fit with a power law spectrum above 10 17 eV gives a spectral index of γ = −3.24 ± 0.08. A preliminary statistical analysis with the F-test shows that the significance associated with this change of the spectral index is at the level of 99.8%. It is noteworthy that a somewhat different spectral feature at 50 − 80 PeV has been reported by the GAMMA experiment (Garyaka et al., 2008), which rather observed a "bump" instead of the spectral break seen by KASCADE-Grande. Figure 4 shows the all-particle energy spectrum together with the results of other experiments. At the threshold, there is a good overlap with the results obtained by KASCADE and EAS-TOP. This is somehow expected as the analysis performed on the subsample of events reconstructed by both, KASCADE and KASCADE-Grande, showed no significant systematical effects. The spectrum doesn't follow a pure power-law on the entire energy range. The features observed around few 10 16 eV and close to 10 17 could give interest- ing hints to the astrophysical processes in the transition region from galactic to extragalactic origin of cosmic rays. As the change of spectral index is small compared to the main knee, a significant conclusion is not possible without investigating the composition in detail in this energy range. These composition studies are presently underway at KASCADE-Grande adopting different approaches based on N ch and N µ observables, among them: a) the reconstruction of the local muon densities at a certain distance to the shower core as it gives a sensitivity to changes in the elemental composition (De Souza et al., 2009); b) the electron/muon number ratio (Cantoni et al., 2009); c) the evolution of the variable k as defined in Eq.

Conclusions
(2); d) an unfolding procedure, similar to what was already performed with the KASCADE data (Antoni et al., 2005). Finally, at the highest energies, the spectrum quite naturally overlaps, inside the statistical and systematic uncertainties, with Hi-Res (Abbasi et al., 2008) and Auger (Schüssler, 2009) results.

Fig. 4.
Comparison of the allparticle energy spectrum obtained with KASCADE-Grande data to results of other experiments. The band around the KASCADE-Grande spectrum denotes the systematic uncertainties in the flux estimation.